Sundaram 筛 - 列表理解

Question

我正在尝试编写一个函数，使用 "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]

我做错了什么？

Answer 1

工作原理

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

修复您的代码

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

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

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

这里我决定只允许i 和j 在1 和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²+2i。如果我们不希望它超过 n，我们需要 2i²+2i2i <= 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+jj<=floor( (n'-i')/(2*i'+1)) 其中i'=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 计算第二个素数！