8

I am trying to write a function that calculates all odd prime numbers from 1..n using the "Sieve of Sundaram" algorithm.

Here is my try:

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

But it gives some wrong numbers like 9,15,21,25, etc.

*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]

What am I doing wrong?

4

1 回答 1

12

这个怎么运作

Sundaram 的 seive 专注于奇数 2n+1，排除那些是数字的乘积。

如果两个数相乘得到一个奇数，它们一定都是奇数，所以我们的数 2n+1 = (2i+1)(2j+1)。如果我们将其相乘，我们得到 2n+1 = 4ij + 2i +2j + 1，我们可以将其简化为 2n=4ij+2i+2j，它再次简化为 n=2ij+i+j。所以如果我们可以把它写成 2ij+i+j，我们就不想要 n。这对于任何数字 i 和 j 都是正确的，但是只需去掉 i<=j 的数字就可以了，因为否则你肯定会两次排除同一个数字。

修复你的代码

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

首先生成排除列表要简单得多：

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

在这里，我决定只允许iandj介于1and之间n，因为否则 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]]

使用一些数学来加快速度

这些的问题是它们产生了太多的数字然后把它们扔掉。只生成我们需要的数字会更快。

实际上，我认为最好计算一下要走多远。j我们将使用的最小值是，i所以 2ij+i+j 可以是 2i 2 +2i 的最小值。如果我们不希望它超过 n，我们需要 2i 2 +2i<=n，我们可以将其重写为 2i(i+1)<=n。正确性比效率更重要，所以可以超过 na 一点，但重要的是不要错过 n 以下的数字，所以我们可以说 2i 2 <=n。这可以表示为i <= floor (sqrt (fromIntegral n / 2))(floor截断小数，所以floor 35.7是 35，fromIntegral在这里用于转换n为浮点数（允许非整数），因此我们可以进行除法和平方根。

这是很多工作，但现在我们可以计算一次应该有多大i：

sSundDelete n = filter (<= n) [i+j+2*i*j|i<-[1..floor (sqrt (fromIntegral n / 2))], j<-[i..n]]

我们可以在 j 上做类似的工作。我们想要 2ij+i+j<=n，我们可以将其重新排列为 (2i+1)j<=ni，这可以用j<=floor( (n'-i')/(2*i'+1))wherei'=fromIntegral i 和来表示n'=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计算第二个素数！

于 2013-04-26T23:54:15.600 回答