埃拉托色尼的分段筛？
Posted
技术标签:
【中文标题】埃拉托色尼的分段筛？【英文标题】：Segmented Sieve of Eratosthenes? 【发布时间】：2012-05-02 05:23:56 【问题描述】：制作一个简单的筛子很容易：
for (int i=2; i<=N; i++)
if (sieve[i]==0)
cout << i << " is prime" << endl;
for (int j = i; j<=N; j+=i)
sieve[j]=1;
cout << i << " has " << sieve[i] << " distinct prime factors\n";
但是当 N 非常大而我无法在内存中保存这种数组时怎么办？我查找了分段筛方法，它们似乎涉及在 sqrt(N) 之前找到素数，但我不明白它是如何工作的。如果 N 非常大（比如 10^18）怎么办？
【问题讨论】：
您在对 larsmans 的回答中提到，您真的有兴趣找到大 N 的质因子数。在这种情况下，假设 N 对于每个 k 到 N，而不仅仅是 N 【参考方案1】：分段筛的基本思想是选择小于 n 平方根的筛分素数，选择一个相当大但适合内存的分段大小，然后筛选每个分段反过来，从最小的开始。在第一段，计算段内每个筛素的最小倍数，然后以正常方式将筛素的倍数标记为合数；当所有的筛选素数都用完后，该段中剩余的未标记数是素数。然后，对于下一个片段，对于每个筛分素数，您已经知道当前片段中的第一个倍数（它是结束前片段中那个素数的筛分的倍数），所以您在每个筛分素数上进行筛分，依此类推直到你完成。
n 的大小无关紧要，只是较大的 n 比较小的 n 需要更长的筛分时间；重要的是段的大小，它应该尽可能大（例如，机器上主内存缓存的大小）。
您可以看到分段筛here 的简单实现。请注意，分段筛将比另一个答案中提到的 O'Neill 的优先队列筛快得多；如果你有兴趣，有一个实现here。
编辑：我写这个是为了不同的目的，但我会在这里展示它，因为它可能有用：
虽然埃拉托色尼筛法非常快，但它需要 O(n) 空间。通过在连续的段中执行筛选，这可以减少到筛选素数的 O(sqrt(n)) 加上位数组的 O(1)。在第一段，计算段内每个筛素的最小倍数，然后以正常方式将筛素的倍数标记为合数；当所有的筛选素数都用完后，该段中剩余的未标记数是素数。然后，对于下一段，每个筛分素数的最小倍数是前一段筛分结束的倍数，所以筛分一直持续到结束。
考虑在 20 段中从 100 到 200 筛分的示例。五个筛分素数是 3、5、7、11 和 13。在从 100 到 120 的第一段中，位数组有十个槽，槽 0对应101，slot k对应100+2k+1，slot 9对应119。segment中3的最小倍数是105，对应slot 2；槽 2+3=5 和 5+3=8 也是 3 的倍数。槽 2 处 5 的最小倍数是 105，槽 2+5=7 也是 5 的倍数。7 的最小倍数是 105在槽 2 处，槽 2+7=9 也是 7 的倍数。以此类推。
函数 primesRange 接受参数 lo、hi 和 delta； lo 和 hi 必须是偶数，lo
function primesRange(lo, hi, delta)
function qInit(p)
return (-1/2 * (lo + p + 1)) % p
function qReset(p, q)
return (q - delta) % p
sieve := makeArray(0..delta-1)
ps := tail(primes(sqrt(hi)))
qs := map(qInit, ps)
while lo < hi
for i from 0 to delta-1
sieve[i] := True
for p,q in ps,qs
for i from q to delta step p
sieve[i] := False
qs := map(qReset, ps, qs)
for i,t from 0,lo+1 to delta-1,hi step 1,2
if sieve[i]
output t
lo := lo + 2 * delta
当调用 primesRange(100, 200, 10) 时，筛选素数 ps 为 [3, 5, 7, 11, 13]； qs 最初是 [2, 2, 2, 10, 8] 对应于最小的倍数 105, 105, 105, 121 和 117，并且对于第二段重置为 [1, 2, 6, 0, 11] 对应于最小123、125、133、121 和 143 的倍数。
您可以在http://ideone.com/iHYr1f 看到这个程序的运行情况。除了上面显示的链接之外，如果您对使用素数进行编程感兴趣，我在我的博客上谦虚地推荐 essay。
【讨论】：
你看了吗？我指出的实现包括一个很好的解释。 您要求提供示例。引用的站点精确地显示了如何在五个段中筛选 100 到 200 的范围，包括如何选择筛分素数以及如何重置每个段的筛分素数。您是否亲自为自己制定了示例？你还有什么不明白的？ 看例子。小于 200 的平方根的筛选素数是 3、5、7、11 和 13。让我们考虑第一段，它有十个值 101 103 105 107 109 111 113 115 117 119。段中 3 的最小倍数是 105，所以打 105 和后面的每三个数字：111、117。段中 5 的最小倍数是 105，所以打 105 和后面的第五个数字：115。段中的 7 是 105，所以打 105 和后面的第七个数字：119。段中没有 11 的倍数，所以无事可做。 13的最小倍数 在segment中是117，所以敲它。剩下的数 101 103 107 109 113 是素数。第二段121 123 125 127 129 131 133 135 137 139每个素数的最小倍数为123、125、133（段外）、121和143（段外），都可以通过计数来计算第一段结束后的下一个筛分素数的倍数。这有帮助吗？ +1 以获得对分段筛和编程实践链接的出色描述。【参考方案2】：有一个基于优先级队列的 Sieve 版本，可以根据您的要求产生尽可能多的素数，而不是所有素数都达到上限。它在经典论文"The Genuine Sieve of Eratosthenes" 中进行了讨论，并且在谷歌上搜索“eratosthenes 优先队列的筛子”会在各种编程语言中找到相当多的实现。
【讨论】：
我遇到过这些实现，但问题是我不理解它们。那些纸总是很密。我主要是在寻找示例，因为我认为这些示例最容易使用和理解。从技术上讲，我正在使用筛子为大 N 获取每个值 k 的唯一素因子 # 个。 Melissa O'Neill 在链接论文中使用的增量筛与基于数组的筛相比非常慢，而且更糟糕的是，其渐近计算复杂度的增长速度远快于随范围线性增长的速度，所以可能不适合这个问题。至于“不需要上限”限定条件，如果基素数小于当前范围的平方根，则页面分段筛也不必具有指定的上限）被实现为“可扩展数组”或作为一种列表形式。 @gordonbgood 对我来说，基于迭代器和优先级队列的筛子“与基于数组的筛子相比非常慢”，这显然不是正确的。运行时间是，据我所知：O（从 i=2 到 n 的总和 log(π(i)-1) ω(i)) （其中 π 是小于或等于给定数的素数数, ω 是给定数的唯一素因数的数量）。基于数组的筛子的一个相对简单的实现是 O（从 i=2 到 n 的总和（如果 i 是素数，则为 n/i，如果 i 不是素数，则为 1））。 @gordonbgood 基本上，我没有数学能力来快速思考，但目前我认为你是对的，前者比后者慢，而前者有比后者更糟糕的无症状生长。但这不是一个非常微不足道的分析，我最初的直觉是怀疑你的评论。我必须像这样明确地说明每次迭代的复杂性，以便我看到您似乎总体上是正确的（除了主观模糊强化词，例如“相当”）。 @gordonbgood 但经过进一步分析，它变得更加不清楚。让我们看看总和中的这些术语：n/i in array-based vs log(π(i)-1) ω(i)。前者趋势从 n 除以一个小常数，趋向于 1。后者从 1 趋向于 log(π(n)-1) ω(n)。这将直觉引诱为“前者缩小，后者增长，所以显然前者更快，后者更慢”。但我发现如果我们对给定的 n 取所有这些和的所有项，并从最小到最大对它们进行排序，那么我发现它是有用的，它们都从 1 开始并爬升到 n/2 和 log(π(n)-1) ω(n)，分别。【参考方案3】：只是我们正在用我们现有的筛子进行分段。 基本思想是假设我们必须找出 85 到 100 之间的素数。 我们必须使用传统的筛子，但方式如下所述：
所以我们取第一个素数 2，将起始数除以 2(85/2)，然后四舍五入到较小的数，我们得到 p=42，现在再次乘以 2，得到 p=84，从这里开始开始添加 2 直到最后一个数字。所以我们所做的是我们已经删除了范围内 2(86,88,90,92,94,96,98,100) 的所有因子。
我们取下一个素数 3，将起始数除以 3(85/3)，然后四舍五入到较小的数，得到 p=28，现在再次乘以 3，得到 p=84，从这里开始将 3 添加到最后一个数字。所以我们所做的就是删除了该范围内 3(87,90,93,96,99) 的所有因子。
取下一个素数=5，以此类推...... 继续做上面的步骤。你可以通过使用传统的筛子直到 sqrt(n) 得到素数 (2,3,5,7,...)。然后将它用于分段筛子。
【讨论】：
【参考方案4】：如果有人想看 C++ 实现，这里是我的：
void sito_delta( int delta, std::vector<int> &res)
std::unique_ptr<int[]> results(new int[delta+1]);
for(int i = 0; i <= delta; ++i)
results[i] = 1;
int pierw = sqrt(delta);
for (int j = 2; j <= pierw; ++j)
if(results[j])
for (int k = 2*j; k <= delta; k+=j)
results[k]=0;
for (int m = 2; m <= delta; ++m)
if (results[m])
res.push_back(m);
std::cout<<","<<m;
;
void sito_segment(int n,std::vector<int> &fiPri)
int delta = sqrt(n);
if (delta>10)
sito_segment(delta,fiPri);
// COmpute using fiPri as primes
// n=n,prime = fiPri;
std::vector<int> prime=fiPri;
int offset = delta;
int low = offset;
int high = offset * 2;
while (low < n)
if (high >=n ) high = n;
int mark[offset+1];
for (int s=0;s<=offset;++s)
mark[s]=1;
for(int j = 0; j< prime.size(); ++j)
int lowMinimum = (low/prime[j]) * prime[j];
if(lowMinimum < low)
lowMinimum += prime[j];
for(int k = lowMinimum; k<=high;k+=prime[j])
mark[k-low]=0;
for(int i = low; i <= high; i++)
if(mark[i-low])
fiPri.push_back(i);
std::cout<<","<<i;
low=low+offset;
high=high+offset;
else
std::vector<int> prime;
sito_delta(delta, prime);
//
fiPri = prime;
//
int offset = delta;
int low = offset;
int high = offset * 2;
// Process segments one by one
while (low < n)
if (high >= n) high = n;
int mark[offset+1];
for (int s = 0; s <= offset; ++s)
mark[s] = 1;
for (int j = 0; j < prime.size(); ++j)
// find the minimum number in [low..high] that is
// multiple of prime[i] (divisible by prime[j])
int lowMinimum = (low/prime[j]) * prime[j];
if(lowMinimum < low)
lowMinimum += prime[j];
//Mark multiples of prime[i] in [low..high]
for (int k = lowMinimum; k <= high; k+=prime[j])
mark[k-low] = 0;
for (int i = low; i <= high; i++)
if(mark[i-low])
fiPri.push_back(i);
std::cout<<","<<i;
low = low + offset;
high = high + offset;
;
int main()
std::vector<int> fiPri;
sito_segment(1013,fiPri);
【讨论】：
【参考方案5】：基于 Swapnil Kumar 的回答，我用 C 语言做了以下算法。它是用 mingw32-make.exe 构建的。
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
int main()
const int MAX_PRIME_NUMBERS = 5000000;//The number of prime numbers we are looking for
long long *prime_numbers = malloc(sizeof(long long) * MAX_PRIME_NUMBERS);
prime_numbers[0] = 2;
prime_numbers[1] = 3;
prime_numbers[2] = 5;
prime_numbers[3] = 7;
prime_numbers[4] = 11;
prime_numbers[5] = 13;
prime_numbers[6] = 17;
prime_numbers[7] = 19;
prime_numbers[8] = 23;
prime_numbers[9] = 29;
const int BUFFER_POSSIBLE_PRIMES = 29 * 29;//Because the greatest prime number we have is 29 in the 10th position so I started with a block of 841 numbers
int qt_calculated_primes = 10;//10 because we initialized the array with the ten first primes
int possible_primes[BUFFER_POSSIBLE_PRIMES];//Will store the booleans to check valid primes
long long iteration = 0;//Used as multiplier to the range of the buffer possible_primes
int i;//Simple counter for loops
while(qt_calculated_primes < MAX_PRIME_NUMBERS)
for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
possible_primes[i] = 1;//set the number as prime
int biggest_possible_prime = sqrt((iteration + 1) * BUFFER_POSSIBLE_PRIMES);
int k = 0;
long long prime = prime_numbers[k];//First prime to be used in the check
while (prime <= biggest_possible_prime)//We don't need to check primes bigger than the square root
for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
if ((iteration * BUFFER_POSSIBLE_PRIMES + i) % prime == 0)
possible_primes[i] = 0;
if (++k == qt_calculated_primes)
break;
prime = prime_numbers[k];
for (i = 0; i < BUFFER_POSSIBLE_PRIMES; i++)
if (possible_primes[i])
if ((qt_calculated_primes < MAX_PRIME_NUMBERS) && ((iteration * BUFFER_POSSIBLE_PRIMES + i) != 1))
prime_numbers[qt_calculated_primes] = iteration * BUFFER_POSSIBLE_PRIMES + i;
printf("%d\n", prime_numbers[qt_calculated_primes]);
qt_calculated_primes++;
else if (!(qt_calculated_primes < MAX_PRIME_NUMBERS))
break;
iteration++;
return 0;
它设置要找到的最大素数，然后用已知的素数（如 2、3、5...29）初始化一个数组。所以我们创建一个缓冲区来存储可能的素数段，这个缓冲区不能大于最大初始素数的幂，在这种情况下是 29。
我确信可以进行大量优化来提高性能，例如并行化段分析过程和跳过 2、3 和 5 的倍数的数字，但它可以作为低内存消耗的示例。
【讨论】：
【参考方案6】：一个数是素数，如果没有一个较小的素数能整除它。由于我们按顺序遍历素数，我们已经将所有可以被至少一个素数整除的数字标记为可整除。因此，如果我们到达一个单元格并且它没有被标记，那么它就不能被任何更小的素数整除，因此必须是素数。
记住以下几点：-
// Generating all prime number up to R
// creating an array of size (R-L-1) set all elements to be true: prime && false: composite
#include<bits/stdc++.h>
using namespace std;
#define MAX 100001
vector<int>* sieve()
bool isPrime[MAX];
for(int i=0;i<MAX;i++)
isPrime[i]=true;
for(int i=2;i*i<MAX;i++)
if(isPrime[i])
for(int j=i*i;j<MAX;j+=i)
isPrime[j]=false;
vector<int>* primes = new vector<int>();
primes->push_back(2);
for(int i=3;i<MAX;i+=2)
if(isPrime[i])
primes->push_back(i);
return primes;
void printPrimes(long long l, long long r, vector<int>*&primes)
bool isprimes[r-l+1];
for(int i=0;i<=r-l;i++)
isprimes[i]=true;
for(int i=0;primes->at(i)*(long long)primes->at(i)<=r;i++)
int currPrimes=primes->at(i);
//just smaller or equal value to l
long long base =(l/(currPrimes))*(currPrimes);
if(base<l)
base=base+currPrimes;
//mark all multiplies within L to R as false
for(long long j=base;j<=r;j+=currPrimes)
isprimes[j-l]=false;
//there may be a case where base is itself a prime number
if(base==currPrimes)
isprimes[base-l]= true;
for(int i=0;i<=r-l;i++)
if(isprimes[i]==true)
cout<<i+l<<endl;
int main()
vector<int>* primes=sieve();
int t;
cin>>t;
while(t--)
long long l,r;
cin>>l>>r;
printPrimes(l,r,primes);
return 0;
【讨论】：
以上是关于埃拉托色尼的分段筛？的主要内容，如果未能解决你的问题，请参考以下文章