logo标签列表

推迟筛法与起始逻辑

javascriptalgorithmprimessieve-of-eratosthenessieve
3

基于这篇Python答案，作者是Will Ness。我使用了该答案的JavaScript版本来实现延迟筛法算法：

function * primes() {
    yield 2;
    yield 3;
    yield 5;
    yield 7;
    const sieve = new Map();
    const ps = primes();
    ps.next() && ps.next();
    for (let p = 3, i = 9; true; i += 2) {
        let s = sieve.get(i);
        if (s !== undefined) {
            sieve.delete(i);
        } else if (i < p * p) {
            yield i;
            continue;
        } else {
            s = 2 * p;
            p = ps.next().value;
        }
        let k = i + s;
        while (sieve.has(k)) k += s;
        sieve.set(k, s);
    }
}

但现在我需要在其中添加start点，我很难理解它，因为这里的逻辑并不简单。

当start是质数时，我需要它成为第一个值。而当start不是质数时，我需要序列以start之后的第一个质数开始。

Will Ness在其中一条评论中建议：

您需要为起始点设计一个有效的筛法字典。获取所需的素数很容易——只需运行此算法，直到sqrt（start），然后您需要找到每个核心素数的下一个（或是上一个？）刚好大于start的倍数，同时注意重复。

将其实现为现实可能并不那么简单（至少对我来说不简单）:|

有人能帮忙更新此算法以实现*primes（start）（最好是像上面的JavaScript代码）吗？

function * primes(start) {

  // how to amend it to support 'start' logic???

}

结论

在Will Ness的出色回答之后，我决定通过一个公共库 - primes-generator，分享我使用的最终代码。所有主要算法都可以在src/sieve.ts中找到。

为什么代码中没有注释？难道每一步都是如此微不足道和清晰吗？i是什么？s是什么？p是什么？为什么代码中没有注释？ - Will Ness
@WillNess 因为我不是这个适配程序的作者，我从过去某个在线资源中复制了 JavaScript 版本。至少，它与您编写的内容非常相似，我提供了链接，并且有很好的文档;) - vitaly-t
有趣，感谢澄清。 :) 我正在撰写一个答案... - Will Ness
@WillNess 我已经无法追溯我从哪里复制的了，但我记得它是完全按原样复制的，没有注释。谢谢你的关注。我本来想尽可能地为它提供500赏金的；我认为这个问题非常值得。 - vitaly-t
1
感谢您的慷慨奖励，以及让我写下这个答案的机会。 :) - Will Ness
显示剩余2条评论
2个回答
3
(更新：在这个答案的底部添加了可行的JavaScript代码)

这是埃拉托色尼筛法：

primes = {2, 3, ...} \ ⋃(p ← primes) {p², p²+p, ...} = {2} ∪ oddPrimes, oddPrimes = {3, 5, ...} \ ⋃(p ← oddPrimes) {p², p²+2p, ...}

其中 \ 表示集合减法，∪ 表示集合并，⋃ 表示集合的大并。

下面通过一个实例来说明：

{2,3,4,5,6,7,8,9,10,11,12,13,14,15,...}
 \             |
   { 4,  6,  8,| 10,   12,   14,   16,   18, ...}
    \          .
             { 9,      12,      15,      18,    21,  ...}
       \ 
                                               { 25, 30, 35, ...}
          \                                     { 49, 56, 63, ...}
            \                                     { 121, 132, 143, ...}
              \  ........

或者对于奇素数，

{3,5,7,9,11,13,15,17,19,21,23,25,27,29,31, ...}
 \                    |
     { 9,      15,    | 21,      27,      33, ...}
   \                          .
                            { 25,            35, ...}
     \                                        { 49, 63, ...}
      \                                         { 121, 143, ...}
       \  ........

问题中提到的代码实现了这种方法。在考虑某个候选者时，sieve 以及循环中的其他变量都处于特定状态。因此，我们只需要直接重新创建这个状态。

假设我们正在考虑 i = 19 作为当前候选者。此时，sieve = { (21, 6) }，p = 5。

这意味着对于一个候选者 i，sieve 包含所有满足 q^2 < i 的质数 q 的倍数，并且 p 是 q 后面的下一个质数。

每个倍数是不小于 i 的最小倍数，并且筛子中没有重复项。那么它就处于一致的状态，并且可以从该点继续执行。

因此，在伪代码中：

primes() = ..... // as before

primesFrom(start) =
  let
  { primes.next()
  ; ps = takeWhile( (p => p^2 < start) , primes() )
  ; p = primes.next_value()
  ; mults = map( (p => let 
                       { s = 2*p 
                       ; n = (start-p)/s  // integer division NB!
                       ; r = (start-p)%s
                       ; m = p + (r>0 ? n+1 : n)*s
                       }
                       ( m, s) )
                , ps)
  }
  for each (m,s) in mults
    if m in sieve
       increment m by s until m is not in sieve
    add (m,s) to sieve

然后执行与之前相同的循环。

根据要求，以下是JS代码：

function *primesFrom(start) {
    if (start <= 2) { yield 2; }
    if (start <= 3) { yield 3; }
    if (start <= 5) { yield 5; }
    if (start <= 7) { yield 7; }
    const sieve = new Map();
    const ps = primesFrom(2);
    ps.next();                 // skip the 2
    let p = ps.next().value;   // p==3
    let psqr = p * p;          // p^2, 9
    let c = psqr;              // first candidate, 9
    let s = 6;                 // step value
    let m = 9;                 // multiple

    while( psqr < start )      // must adjust initial state
    {
         s = 2 * p;            
         m = p + s * Math. ceil( (start-p)/s );  // multiple of p
         while (sieve.has(m)) m += s;
         sieve.set(m, s);
         p = ps.next().value;
         psqr = p * p;
    }
    if ( start > c) { c = start; }
    if ( c%2 === 0 ) { c += 1; }
    
    for ( ;  true ;  c += 2 )     // main loop
    {
        s = sieve.get(c);
        if (s !== undefined) {
            sieve.delete(c);      // c composite
        } else if (c < psqr) {
            yield c;              // c prime
            continue;
        } else {                  // c == p^2
            s = 2 * p;
            p = ps.next().value;
            psqr = p * p;
        }
        m = c + s;
        while (sieve.has(m)) m += s;
        sieve.set(m, s);
    }
}

正确地在ideone上瞬间找出首10个大于500,000,000的质数：链接。

Success #stdin #stdout 0.03s 17484KB
500000003,500000009,500000041,500000057,500000069,
500000071,500000077,500000089,500000093,500000099

显然，它以惊人的递归深度 5（五次）调用来实现。

重复平方的威力！或者它的反操作，log log运算：

log2( log2( 500000000 )) == 4.85

1
我需要一个简单的东西，只需点击运行并查看它的运行情况。实际上，我可以在ideone.com上自己做到这一点，我看到它有JS。 :) 如果我有任何问题，我会在这里问你。 :) --- 至于准备好的东西，这里有一个。 - Will Ness
期待您的完整解决方案。我已经尝试过跟随那个基于Haskell的实现的链接，甚至尝试将其编译成JavaScript，但目前没有太大的成功。我认为很多人现在都希望看到它在JavaScript中的实现 ;) - vitaly-t
我开始认为“筛法”算法并不像它被赞扬的那样高效。我能够创建一个高度优化的isPrime版本，它比筛法更快地给出质数。 - vitaly-t
我已经修改了我的代码，使其返回5000亿后的前500个质数，在ideone上只需要0.3秒。对于2500亿，只需要0.15秒。因此，它的增长是线性的。我怀疑你的代码不会表现得那么好。 - Will Ness
看看我添加的测试，就像你要求的那样 - 它轻松地处理那些高数字，比筛法快10倍。让我知道你的想法 ;) - vitaly-t
显示剩余17条评论
0

我修改了Will Ness的惊人答案

  1. 允许任何起始值进行打印。如下所述，如果您想要更高的质数的无限序列，则我认为无法避免必须计算所有早期质数。

  2. 调整变量名称以帮助理解算法。

  3. 将存储在映射中的内容从2倍质因子改为仅质因子，以使读者更容易跟随算法。

  4. 将代码的一部分移动到子函数中，以便读者理解。

  5. 更改算法中间三选一的控制流，并添加注释以简化理解。它可能稍微慢一些，因为它不再首先测试最常见的真条件，但对于读者来说更容易。

function* primeGenerator() {
  yield 2;
  yield 3;
  yield 5;
  yield 7;
  const multiplesWithAFactor = new Map();

  // We only need to consider potential factors  that are themselves prime
  // Start with 3 and get further prime factors on demand from a child instance of ourself.
  let biggestFactorConsidered = 3;
  const childStreamOfPrimes = primeGenerator();
  childStreamOfPrimes.next(); // Discard the 2
  childStreamOfPrimes.next(); // Discard the 3

  let candidate = 7; // We have already yielded 7 as a prime above, so on the first iteration of the while loop we will be testing 9.
  while (true) {
    candidate += 2;

    // Assess candidate, into one of three outcomes: square, non-square multiple, prime. This order is arranged for ease of understanding, but for maximum speed efficiency you should test in the order in Will Ness's version.
    const factorOfCandidate = multiplesWithAFactor.get(candidate);
    
    if (candidate === biggestFactorConsidered * biggestFactorConsidered) {
    // A square, so from now on we will need to consider the next factor, too.
      markNextUnmarkedMultiple(candidate, biggestFactorConsidered);
      biggestFactorConsidered = childStreamOfPrimes.next().value;
    } else if (factorOfCandidate) {
    // A non-square multiple. Can forget that fact now, and instead mark the next unmarked multiple of the same prime factor
      multiplesWithAFactor.delete(candidate);
      markNextUnmarkedMultiple(candidate, factorOfCandidate);
    } else {
    // Prime
      yield candidate;
    }
  }

  // This is a subfunction for ease of understanding, but for maximum efficiency you should put this in the while loop above, and have the "yield candidate" avoid calling it, via a "continue" statement.
  function markNextUnmarkedMultiple(candidate, factor) {
    let nextMultiple = candidate + 2 * factor;
    while (multiplesWithAFactor.has(nextMultiple)) {
      nextMultiple += 2 * factor;
    }
    multiplesWithAFactor.set(nextMultiple, factor);
  }
}

const maxItems = 20;
const minimumPrintablePrime = 5e8;
console.log("Running");
const primeObject = primeGenerator();
let count = 0;
do {
  const prime = primeObject.next().value;
  if (prime > minimumPrintablePrime) { // If you don't like seeing this process of reading non-printed primes out here in the parent code, you can do it within the primeGenerator itself. Either way, _someone_ has to call the prime generator for all primes up to the square root of the starting value.
    console.log(prime);
    count++;
  }
} while (count < maxItems);

这个算法在存储方面非常经济

它使用了几个素数生成器的实例。您调用“主”生成器，它沿着检查候选项运行，同时存储一组未来的候选项，这些候选项是较小素数因子的倍数。

在考虑候选项N时，它存储的映射包括每个素数因子的条目，直到sqrt(N)，但存储在下一个素数因子的倍数的映射中，该算法将从当前候选项向前探索。

每当候选项出现在此多重映射中时，它都可以拒绝该候选项。该映射告诉我们，对于该倍数，是哪个素数因子告诉我们该数字是倍数。例如，15被标记为3的倍数。当算法到达15时，它意识到它是3的倍数，因此拒绝了15。

从现在开始，它不再需要记住15是一个倍数，因为在此算法实例中，15将永远不会再成为候选项。它可以从其多重映射中删除15条目，并用3的下一个倍数替换它。然而，下一个倍数总是偶数，因为所有的候选项（例如15）和所有的因子（例如3）都是奇数。因此，它只需要告诉映射15+2x3是一个倍数，即它可以通过2x因子向前跨越。此外，如果结果数字21已经标记为倍数，它不需要标记该数字：它可以继续前进2x3到15+2x3+2x3等，直到找到一个未标记为倍数的数字。
巧妙的是，这个过程确保每个因子（例如3）都永远在映射中标记一个倍数。为了评估候选N，映射将仅包含每个小于sqrt(N)的质数的一个条目。

当候选数超过目前考虑的最大因子的平方时，算法需要获取下一个因子。它只需使用自身的第二个实例来获取序列中的下一个质数。起初我担心这会创建巨大的内存需求，但事实并非如此。评估一个约为2^64的数字，需要所有第二个实例的~2^32，这反过来又会调用第三个实例的~2^16，以此类推。即使对于巨大的质数，也只需要创建少量生成器的实例，如果一个实例的映射具有高达F的因子，则下一个实例仅需要具有高达sqrt(F)的因子，这很快就变得很小。

即使是Ness算法仍然需要构建包含所有因子的映射，直到sqrt(N)

它需要映射包含所有这些因子，以便可以正确地拒绝N附近的候选数。

但它最终也需要sqrt(N)和N之间的因子

只有这样，它才能确信从N到N^2的数字是质数。

我的结论

我认为需要迭代才能使其工作。我无法想象在任意候选数（例如10^8）开始而不预先填充地图与（所有）素数直到10^4。

最有效的预填充方法是运行算法本身，如果您只是要求它产生（但不打印）所有素数直到10^8，则实际上就是这样做的。

我的怀疑是该算法非常高效，以至于将其置于任意候选数（例如1e8）的位置的最佳方法是通过运行算法本身，即没有捷径。这是一个非常特殊的情况！

这似乎在威尔更新的答案中得到了证实，在其中算法调用自身来预填筛子。在预填充结束时，筛子包含每个小于sqrt（start）的素数的一个条目。我的版本要慢得多，但这不是因为它收集比威尔更多的素数；而是因为它（a）按不同顺序测试3个条件（使其更易于人类跟随），并且（b）将一些代码提取到子函数中（再次使其更易于理解）。显然，威尔的方法对于生产代码更好！

@will Ness - 谢谢！我的修改后的答案正确吗？您似乎非常擅长质数生成。如果您想要一个可以无限生成的生成器，我是否正确地认为它需要明确或隐含地知道所有起始值之前的所有质数？ - ProfDFrancis
隐含地，宇宙中的每一个真实的事物都是真实的，仅仅因为它是真实的。但我们从事的是显式知识的业务。埃拉托色尼筛法是 primes = [2,3,...] \ U [[p*p,p*p+p,...] for p in primes]。或者更详细地说，primes = nats2 \ U composites, composites = map multiples primes, multiples(p) = [p*p,p*p+p,...] -- 每个质数 p 的倍数从 p*p 开始枚举。所以要找到一个质数 **q**，只需要知道显式的小于等于 sqrt(q) 的质数即可。所有小于等于 q 的质数的知识都隐含在小于等于 sqrt(q) 的质数中。 - Will Ness
我们需要知道的是，第一个质数是2，其他的都会随之而来。@Eureka - Will Ness
我也在这里发布了一个答案。特别是它描述了一种快捷方式，并且表明筛法中没有考虑任何因素。那将是“试除法”筛法的方式，而不是“埃拉托斯特尼筛法”。 :) - Will Ness
FYI，我已经在我的答案中添加了一个可工作的JS代码。毕竟有捷径。 :) - Will Ness
显示剩余5条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接