埃氏筛的优化——Wheel Factorization 学习笔记

· · 算法·理论

Wheel Factorization

埃氏筛可以在 O(n\log \log n) 的时间复杂度内筛出 [1,n] 内的所有质数；线性筛可以进一步做到 O(n)。

Wheel Factorization 就是对埃氏筛的一种优化，可以在低于线性的时间内筛出所有质数。

Algorithm

主要思想：分块筛选。

还是考虑如何减少重复标记的次数。

设立一个阈值 B>1，先用线性筛筛出 [1,B] 的质数。

考虑一个结论：

记 A_k(i)=1/0 表示 i 是否与 k 互质，那么对于正整数 i 有 A_k(i)=A_k(i+k)。

那么对于每个 x\in[1,B]，若 x 不与 B 互质，则形如 kB+x 的数肯定都不与 B 互质，这些数一定是合数，之后筛的过程中不需要考虑。

接下来对值域分块，每 B 个为一段，分成 [1,B],[B+1,2B],[2B+1,3B],\cdots。依次考虑每段（第一段已经筛过），从上面分析知道只需考虑每个区间剩下的 \varphi(B) 个数。类似埃筛，枚举之前区间得到的（不被 B 包含的）质数 p，从 p^2 开始筛这个区间，这样就跳过了所有不与 B 互质的数，且每个被跳到的数只会被最小的那个质因子筛一次。

于是时间复杂度为 O(B+n\frac{\varphi(B)}{B})，空间复杂度取决于质数个数和块长，O(\frac{n}{\ln n}+B)。

考虑让 \frac{\varphi(B)}{B} 尽量小，质因子次数 >1 没有意义，且让 B 中包含质因子的数量更多。对于 10^9 范围的筛质数，取 B=2\times 3\times 5\times 7\times 11\times 13\times 17=510510 较优，此时 \frac{\varphi(B)}{B} 略大于 0.1。

Code

实现的若干小优化：

SP6488

// Problem: PRIMES2 - Printing some primes (Hard)
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/SP6488
// Memory Limit: 1 MB
// Time Limit: 2280 ms
// 
// Powered by CP Editor (https://cpeditor.org)

#include<bits/stdc++.h>
using namespace std;
const int N=5.1e7+5;
const int B=510510; // 2*3*5*7*11*13*17
const int PB=92160;
const int SZ=1959;

inline void write(int x){
    static char buf[20];
    static int len=-1;
    if(x<0)putchar('-'),x=-x;
    do buf[++len]=x%10,x/=10;while(x);
    while(len>=0)putchar(buf[len--]+48);
}

int pr[N],cnt;
bitset<B+5> vs;
int od[B+5],idx;

void init(int n){
    vs.set(1);
    for(int i=2;i<=B;++i){
        if(!vs[i]) pr[++cnt]=i;
        for(int j=1;j<=cnt&&i*pr[j]<=B;++j){
            vs.set(i*pr[j]);
            if(i%pr[j]==0) break;
        }
    }
    for(int i=1;i<=B;++i){
        if((i&1)&&(i%3)&&(i%5)&&(i%7)&&(i%11)&&(i%13)&&(i%17))
            od[++idx]=i;
    }
    vs.reset();
    for(int i=2;i<=SZ;++i){
        int r=i*B,l=r-B+1;
        if(i==SZ) r=n;
        vs.reset();
        for(int j=8;j<=cnt&&pr[j]*pr[j]<=r;++j){
            int u=pr[j]*max(pr[j],(l-1)/pr[j]+1);
            if(!(u&1)) u+=pr[j];
            u-=l-1;
            while(u<=B) vs.set(u),u+=(pr[j]<<1);
        }
        for(int j=1;j<=idx;++j) if(!vs[od[j]]) pr[++cnt]=od[j]+l-1;
    }
}

signed main(){
    int n(1e9); init(n);
    for(int i=1;i<=cnt&&pr[i]<=n;i+=500) write(pr[i]),puts("");
    return 0;
}

手写 bitset：

typedef unsigned long long ull;
ull vs[(B>>6)+5];
#define set(x) vs[x>>6]|=(1ull<<(x&63)) // vs.set(x)
#define ask(x) (vs[x>>6]>>(x&63)&1) // vs[x]
// 清空直接 memset

练习：SP6489（极度卡常）