更快的质数筛: Wheel Factorization
能在短时间内筛到 \(10^9\)。运用了类似分块的思想
如何做
首先,我们取一个合数 \(k\)。先把 \([1,k]\) 之间筛一下。
然后我们顺便得到哪些数和 \(k\) 互质。
有这样一个引理:
对于 \(k\),设 \(Q(i)\) 表示 \(i\) 是否和 \(k\) 互质,是为 \(1\),否为 \(0\)。则:
\(Q(i)=Q(i+k),\forall i\in \N^{*}\)。即 \(Q\) 的周期是 \(k\)。
接下来,我们连续枚举长度为 \(k\) 的区间。即,\([k+1,2k],[2k+1,3k],...\)。对于这些区间中的每个数,它们是否和 \(k\) 互质的情况,根据那个引理,我们可以知道。而我们从 \([k+1,2k]\) 开始枚举,它们都 \(>k\)。显然,一个 \(>k\) 的数,如果和 \(k\) 不互质,一定不可能是质数,因此就被筛掉了。
对于每个区间,我们只需要考虑剩下的 \(\varphi(k)\) 个和 \(k\) 互质的就行了。那它们一定没有 \(k\) 中的质因子。也就是说,我们在筛质数的时候,跳过 \(k\) 的质因子,就可以完美遍历到这些和 \(k\) 互质的数了。
因此复杂度是 \(O(n\times \dfrac{\varphi(k)}{k})\)。写一个类似埃氏筛的东西,复杂度就是对的了。
如何让 \(\dfrac{\varphi(k)}{k}\) 尽可能小呢?首先,\(k\) 中的质因子次数超过 \(1\) 肯定是浪费的,无效;然后 \(k\) 的质因子越多越好
对于 \(1e9\) 的范围,取 \(k=2\times 3\times 5\times 7\times 11\times 13\times 17=510510\),事实证明,它就够快了。 此时
\(\dfrac{\varphi(k)}{k}\approx 0.1805\)
模板题:spoj某题
代码
#include <bits/stdc++.h>
using namespace std;
namespace Flandre_Scarlet
{
#define N 510515
#define M 181000000
#define ull unsigned long long
#define F(i,l,r) for(int i=l;i<=r;++i)
#define D(i,r,l) for(int i=r;i>=l;--i)
#define Tra(i,u) for(int i=G.h[u],v=G.to(i);~i;i=G.nx(i),v=G.to(i)) if (i>=0)
#define MEM(a,x) memset(a,x,sizeof(a))
#define FK(a) MEM(a,0)
#define sz(x) ((int)x.size())
#define all(x) x.begin(),x.end()
#define p_b push_back
#define pii pair<int,int>
#define fir first
#define sec second
int K=510510,phiK=92160; // 2*3*5*7*11*13*17, 前7个
int pr[60000000]; bool notp[N];
int pos[N],id[N]; // 互质位置, 互质数的编号
bool tmp[N];
#define pid(x) ((x/K)*phiK+id[x%K])
void Init()
{
int n=1000000000; // 1e9
int&c=pr[0],&cp=pos[0];
notp[1]=1; c=cp=0;
F(i,2,K) // 先筛一下 [1,K]
{
if (!notp[i]) pr[++c]=i;
for(int j=1;j<=c and pr[j]<=K/i;++j)
{
notp[i*pr[j]]=1;
if (i%pr[j]==0) break;
}
}
F(i,1,K) if (i%2 and i%3 and i%5 and i%7 and i%11 and i%13 and i%17) // 互质
{
++cp;
pos[cp]=i;
id[i]=cp;
}
F(i,2,1959)
{
int r=i*K,l=r-K+1;
FK(tmp);
// 每次的tmp[i]表示l+i-1是否不是质数
// tmp[i]=1表示被筛掉, 最后tmp[i]=0说明是质数
for(int j=8;j<=c and pr[j]*pr[j]<=r;++j) // pr[j]筛到的最小数是pr[j]*pr[j], 比它小的会更早被筛到
{
int fir=pr[j]*max(pr[j],(l-1)/pr[j]+1);
if (!(fir&1)) fir+=pr[j];
fir-=(l-1);
while(fir<=K) tmp[fir]=1,fir+=(pr[j]<<1); // 只需要筛奇数倍, 因为偶数倍含有2, 肯定不互质
}
F(j,1,cp) if (!tmp[pos[j]]) pr[++c]=l+pos[j]-1;
}
for(int i=1;i<=c and pr[i]<=n;i+=500)
{
printf("%d\n",pr[i]);
}
}
void Input()
{
}
void Sakuya()
{
}
void IsMyWife()
{
Init();
Input();
Sakuya();
}
}
#undef int //long long
int main()
{
Flandre_Scarlet::IsMyWife();
return 0;
}