无需大内存与强力CPU，普通电脑也能用的一万亿以内质数计算器！

最后的完整C++代码在文末，需要请自取。

最近学习了质数筛法，想研究一下计算质数的最快速方法，于是便有了这篇文章。

Part 1 试除法

根据质数的定义，质数是除了1和它本身没有任何其他因数的数，由此得出了我们的第一版算法，枚举所有大于 \(1\) 小于 \(x\) 的数，判断是否能被它整除：

bool isPrime1(int x){
    if(x<2) return false;//质数不能小于2
    if(x==2||x==3) return true;//2和3都是质数
    for(int i=2;i<x;i++)
        if(x%i==0) return false;
    return true;
}

在计算 \(n\) 以内的质数时，这个算法的复杂度是 \(\operatorname{O}(n^2)\)，在普通电脑上大约每秒能计算10万以内的质数，约1万个。

我们很容易发现，我们并没有必要枚举这么多数，因为一个数的一对因数一定有一个小于等于 \(\sqrt{x}\)，所以我们只要枚举 \(\sqrt x\) 以内的数就可以了。由此得出第二版代码：

bool isPrime1(int x){
    if(x<2) return false;//质数不能小于2
    if(x==2||x==3) return true;//2和3都是质数
	int limit=sqrt(x);//最大根号x
    for(int i=2;i<limit;i++)
        if(x%i==0) return false;
    return true;
}

这个代码计算 \(n\) 以内的质数复杂度为 \(\operatorname{O}(n\sqrt{n})\)，大约能够计算一百万以内的质数。

这看起来很快……但是还不够。于是，筛法便应运而生。

Part 2 筛法

埃拉托斯特尼筛法（Eratosthenes 筛法，简称埃氏筛），是质数筛法中最简单的一种。

它的思想很简单：考虑这样一件事情：对于任意一个大于 \(1\) 的正整数 \(x\)，他的 \(k(k>=2)\) 倍数一定是合数。利用这个结论，我们可以避免很多次不必要的检测。

如果我们从小到大考虑每个数，然后同时把当前这个数二倍以上的倍数记为合数，那么运行结束的时候没有被标记的数就是质数了。

为什么能保证正确性？我们在枚举到 \(n\) 时，\(n\) 以前的数都已经被筛过，也就是小于 \(n\) 大于 \(1\) 的数都不是 \(n\) 的因数，即 \(n\) 是质数。由此也保证了复杂度。

具体实现见代码：

vector<int> prime;
bool is_prime[N];

void Eratosthenes(int n) {
  is_prime[0] = is_prime[1] = false;
  for (int i = 2; i <= n; ++i) is_prime[i] = true;
  for (int i = 2; i <= n; ++i) {
    if (is_prime[i]) {
      prime.push_back(i);
      if ((long long)i * i > n) continue;
      for (int j = i * i; j <= n; j += i)
        // 因为从 2 到 i - 1 的倍数我们之前筛过了，这里直接从 i 的倍数开始，提高了运行速度
        is_prime[j] = false;  // 是 i 的倍数的均不是素数
    }
  }
  //此时，prime[i] 若为 1，则说明 i 是质数
}

该算法复杂度为 \(\operatorname{O}(n \log \log n)\)，具体证明较为复杂，详见OI Wiki上的证明部分，这里不再赘述。

注意，受限于内存大小，普通埃氏筛法在 \(8\) GB 内存的加持下大约只能计算 \(8.5 \times 10^8\)（八亿五千万）以内的质数。

这就是极限了吗？

不。在埃氏筛之上，还有一种线性筛法——欧拉筛，这种筛法的理论时间复杂度为 \(\operatorname{O}(n)\)，也即线性时间复杂度。

我们注意到在埃氏筛中一个合数可能会被标记多次，如果能让每个合数都只被标记一次，那么时间复杂度就可以降到线性了。

具体解释请看代码：

vector<int> pri;
bool not_prime[N];

void pre(int n) {
  for (int i = 2; i <= n; ++i) {
    if (!not_prime[i]) {
      pri.push_back(i);
    }
    for (int pri_j : pri) {
      if (i * pri_j > n) break;
      not_prime[i * pri_j] = true;
      if (i % pri_j == 0) {
        // i % pri_j == 0
        // 换言之，i 之前被 pri_j 筛过了
        // 由于 pri 里面质数是从小到大的，所以 i 乘上其他的质数的结果一定会被
        // pri_j 的倍数筛掉，就不需要在这里先筛一次，所以这里直接 break
        // 掉就好了
        break;
      }
    }
  }
}

Part 3 内存不够了怎么办？

理论来说，埃氏筛一秒钟可以计算大约一亿个质数，但是在欧拉筛和埃氏筛中，保存 \(n\) 以内的质数在值域很大时会占用很多内存，普通电脑无法承受。

在 Part 1 中我们就有过一个结论：只需要筛 \(\sqrt n\) 以内的因子，就能判断 \(n\) 是否为质数。

我们采用分块思想，不需要一直保留整个 is_prime[] 数组，而是只保留 \(\sqrt{n}\) 以内的质数，用这些质数去筛剩下的质数，而这就是今天的终极算法：分块埃氏筛。

我们把 \(n\) 以内的数分成 \(\sqrt{n}\) 块，于是每块内也含有 \(\sqrt{n}\) 个数。先计算出第一块内（1 到根号 n）的所有质数，然后对于每个块依次用这些质数去筛选，这样空间复杂度就变成了 \(\operatorname{O}(\sqrt{n})\)。

以下是分块埃氏筛的实现：

int count_primes(int n) {
  constexpr static int S = 10000;
  vector<int> primes;
  int nsqrt = sqrt(n);
  vector<char> is_prime(nsqrt + 1, true);
  for (int i = 2; i <= nsqrt; i++) {
    if (is_prime[i]) {
      primes.push_back(i);
      for (int j = i * i; j <= nsqrt; j += i) is_prime[j] = false;
    }
  }
  int result = 0;
  vector<char> block(S);
  for (int k = 0; k * S <= n; k++) {
    fill(block.begin(), block.end(), true);
    int start = k * S;
    for (int p : primes) {
      int start_idx = (start + p - 1) / p;
      int j = max(start_idx, p) * p - start;
      for (; j < S; j += p) block[j] = false;
    }
    if (k == 0) block[0] = block[1] = false;
    for (int i = 0; i < S && start + i <= n; i++) {
      if (block[i]) result++;
    }
  }
  return result;
}

为什么不能用分块欧拉筛？

经典的欧拉筛（线性筛）依赖一个全局数组记录每个数的最小质因子（或一个 low 数组），同时需要在筛的过程中知道当前质数能筛到的所有合数，并且这些合数是连续的、递推增长的。

当范围很大时，不可能一次性分配这么大的数组。

分段后，当前块内某个合数的最小质因子可能不在块内的质数表中，而是更小的质数，而这些质数在上一个块中已经处理完毕。如果不保留跨块信息，就无法知道该合数是否已经被更小的质数筛过。

Part 4 质数计算器

说了这么多，我们放一个（理论上）能计算 \(10^{16}\)（万万亿）级别的质数计算器。

实际使用中，受限于 CPU 性能，在可接受的时间（约几个小时到一天不等）的时间里，能计算出 \(10^{12}\) （一万亿）以内的质数并保存到文件中。

#include<iostream>
#include<fstream>
#include<vector>
#include<cmath>
#include<conio.h>
#include<windows.h>
#include<time.h>
#include<bitset>
#include<sstream>
#define int long long
using namespace std;
int stt,yyy;
void v(){};
int calc(int n){
    ofstream out(to_string(n)+"underprime.txt");
    const int sz=sqrt(n);
    vector<int> prime;
    vector<bool> ispri(sz+1,1);
    int res=0;
    for(int i=2;i<=sz;i++){
        if(ispri[i]){
            prime.push_back(i);
            for(int j=i*i;j<=sz;j+=i){
                ispri[j]=0;
            }
        }
    }
    vector<int> block(sz);
    for(int k=0;k*sz<=n;k++){
        if(k%(sz/(max((int)((int)log10(n)*(int)log10(n)*0.25-20),1ll)*100))==0||k==sz) printf("(%.2fs)Compelete  %.2f%%  (calculating block%lld/%lld)\n", (clock()-stt)/1000.0, k*100.0/sz, k, sz);
        for(register int i=0;i<sz;i++) block[i]=1;
        int st=k*sz;
        for(int i:prime){
            int x=(st+i-1)/i;
            int j=max(x,i)*i-st;
            for(;j<sz;j+=i) block[j]=0;
        }
        if(k==0) block[0]=block[1]=0;
        stringstream buffer;
        int bufcnt=0;
        for(int i=0;i<sz&&st+i<=n;i++){
            if(block[i]){
                res++;
                if(yyy){
                    buffer<<i+k*sz<<",";
                    bufcnt++;
                    if(bufcnt>=10000){
                        out<<buffer.rdbuf();
                        buffer.str("");
                        bufcnt=0;
                    }
                }
            } 
        }
        if(bufcnt>0) out<<buffer.rdbuf();
    }
    return res;
}
int bl(int n){
    ofstream out(to_string(n)+"underprime.txt");
    int res=0;
    for(int i=2;i<=n;i++){
        int yn=1;
        for(int j=2;j*j<=i;j++){
            if(i%j==0){
                yn=0;
                break;
            }
        }
        if(yn) res++;
        if(yn&&yyy) out<<i<<",";
    }
    return res;
}
signed main(){
    int n;
    printf("需要计算质数的范围上界：");
    cin>>n;
    yyy=0;
    while(n<2) cout<<"输入不合法！\n：",cin>>n;
    double cost=((n/log(n))*(log10(n)+1)/1024/1024)*1.04;
    printf("将自动写入至%s，文件大小约%.2lfMB，是否继续？（按Y/y写入文件，按N/n不写入文件，不继续请直接关闭程序）\n", (to_string(n)+"underprime.txt").c_str(),cost);
    char ch=getch();
    while(ch!='y'&&ch!='n') ch=getch();
    if(ch=='y') yyy=1;
    stt=clock();
    int m=(n<=100000?bl(n):calc(n));
    printf("共%lld个质数，计算用时%lldms(合%.2f秒)%s\n", m, clock()-stt, (clock()-stt)/1000.0,yyy==1?("，存入文件"+(to_string(n)+"underprime.txt")).c_str():"");
    printf("按任意键退出程序......");
    getch();
}

代码风格奇特，勿喷qwq

这个程序可以计算很大范围的质数，在时间花费较长时还会显示当前进度。

运行环境：Windows 7+

程序可选择是否写入文件并帮你预估一下文件的大小，这里给出一些文件大小精确值（亲自验证）：

范围	文件大小
\(10^8\)（一亿）	49MB
\(10^9\)（十亿）	479MB
\(10^{10}\)（百亿）	4.6GB
\(10^{12}\)（万亿）	451GB

在搭载 R5 5600X（超频至4.9GHz）和 16GB DDR4 内存的平台上，计算一万亿以内的质数共用时约四小时。