我正在尝试编写一个函数,使用"Sieve of Sundaram" algorithm从1..n计算所有奇数质数。

这是我的尝试:

sSund :: Integer -> [Integer]
sSund n = [ i * 2 + 1 | i <- [1..n], j <- [f i], (i + j + 2 * i * j) > n ]
  where f 1 = 1
        f y = y + 1 --use function f because i don't know how insert 1 into j's list

但是它给出了一些错误的数字,例如9,15,21,25等。
*Main> sSund 30
[7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59,61]

我究竟做错了什么?

最佳答案

怎么运行的

Sundaram的seive的工作方式是着眼于2n + 1的奇数,并排除那些由数字乘积而来的数。

如果两个数字相乘得到一个奇数,则它们都必须都是奇数,因此我们的数字2n + 1 =(2i + 1)(2j + 1)。如果将其相乘,我们得到2n + 1 = 4ij + 2i + 2j + 1,可以将其简化为2n = 4ij + 2i + 2j,再次将其简化为n = 2ij + i + j。因此我们不希望n可以写为2ij + i + j。对于任何数字i和j都是如此,但是可以摆脱i
修正程式码

在您的代码中,您生成了一些要排除的数字i + j + 2 * i * j,但实际上您只是排除了i而不是i + j + 2 * i * jj<-[f i]仅在列表中为您提供一个j值,而不是从in的所有数字,您应该将其写为[i..n]

首先生成排除列表会更简单:

sSundDelete :: Integer -> [Integer]
sSundDelete n = [i+j+2*i*j|i<-[1..n], j<-[i..n]]

在这里,我决定只允许ij1n之间,因为否则2ij + i + j肯定大于n。

现在我们可以列出不包含这些数字的x数字,然后使用2*n+1公式使它们成为奇数:
sSund :: Integer -> [Integer]
sSund n = let del = sSundDelete n in
     2:[2*x+1 | x <- [1..n], not (x `elem` del)]

哪个正确地给你
> sSund 30
[2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61]

加快速度

但是,它的速度不如预期的快,因为如果您查看
> sSundDelete 10
[4,7,10,13,16,19,22,25,28,31,12,17,22,27,32,37,42,47,52,24,31,38,45,52,59,66,73,40,49,58,67,76,85,94,60,71,82,93,104,115,84,97,110,123,136,112,127,142,157,144,161,178,180,199,220]

它的数字比我们需要的大得多-sSund 10只能达到2 * 10 + 1 = 21。这意味着我们要一次又一次地根据我们从未考虑过的数字检查我们的数字!

要做的最简单的事情是重写sSundDelete
sSundDelete n = [i+j+2*i*j|i<-[1..n], j<-[i..n],i+j+2*i*j<=n]

非常像您一样,或者
sSundDelete n = filter (<= n) [i+j+2*i*j|i<-[1..n], j<-[i..n]]

使用一点数学来加快速度

这些问题是它们产生太多的数字,然后丢掉它们。仅生成我们需要的数字会更快。

实际上,我认为最好算出要走多远。我们将使用的最小ji,因此2ij + i + j可以最小的是2i2 + 2i。如果我们不希望其超过n,则希望2i2 + 2i i <= floor (sqrt (fromIntegral n / 2))(floor截断小数,因此floor 35.7为35,此处使用fromIntegraln转换为浮点数(允许使用非整数),以便我们进行除法和平方根运算。

可以解决很多问题,但是现在我们只需计算一次i应该多大:
sSundDelete n = filter (<= n) [i+j+2*i*j|i<-[1..floor (sqrt (fromIntegral n / 2))], j<-[i..n]]

我们可以对j做类似的工作。我们想要2ij + i + j j<=floor( (n'-i')/(2*i'+1)),其中i'=fromIntegral in'=fromIntegral n。这给了我们
sSundDelete n = [i+j+2*i*j|let n'=fromIntegral n,
                           i<-[1..floor (sqrt (n' / 2))],
                           let i' = fromIntegral i,
                           j<-[i..floor( (n'-i')/(2*i'+1))]]

这样一来,我就不必放弃等待sSund 5000来计算第二个质数的速度了!

07-28 06:51