筛法
概述
- 筛法原本是筛取质数的一种算法。在 OI 中,它被推广到了筛积性函数值。
埃氏筛
-
枚举每个质数,筛去它们的整数倍。
-
复杂度为 \(O(n\log\log n)\),证明非常困难,我们从心一下。
-
埃氏筛有很多 trick,例如只取 \(\leqslant \sqrt{n}\) 的质数来筛,分块筛(这个是时间换空间),etc.,但大都不实用。
欧拉筛(线性筛)
-
枚举所有数 \(i\),筛去他们的质数 \(p_i\) 倍。
-
证明:
-
显然所有合数都可以被写成 \(i\times p\)(\(p\) 是质数)。
-
当 \(i\bmod p\neq 0,i\times p\) 在枚举至 \(i\) 时筛掉;
-
当遇到第一个 \(i\bmod p=0\) ,显然这个 \(p\) 是 \(i\) 的最小质因子(也就是最小因子),记为 \(p_0\) ,把 \(i\times p_0\) 筛掉然后 break;
-
对于 \(p>p_0\) 的 \(i\times p\),\(i\times p=p_0\times k\times p=p_0\times (k\times p)\)。
-
其中 \(p_0<p,k\times p>i\) ,从而会在枚举到更大的 \(i\) 时删掉。
-
综上,每个数只会被筛掉一次,所以复杂度为线性。
-
-
这一证明有如下等价形式:将每个数写成质数-指数向量,则 \(n\) 由 \(\dfrac{n}{lowbit(n)}\) 筛出,这里 \(lowbit(n)\) 表示最低非零位的位权。
-
小技巧:可以通过记录每个数是被谁筛的来递归求解每个数的质因数分解。单次复杂度 \(O(\sum\limits_{p_i\mid x} k_i)\)。
-
欧拉筛也可以使用根号加速的优化,但通常很少用:
-
如果不根号加速会 T,那么所筛的范围根本存不下。
-
如果是在卡常...好吧,这还是有一点意义的。
-
不过事实上一般欧拉筛不是简单的筛质数,求某些积性函数的前缀和就已经 \(O(n)\) 了,别老惦记卡这点不在瓶颈上的常数了。
-
-
给出示范代码:
bool np[maxv];
int p[maxv],ptot;
il void init(){
np[1]=1; int res;
For(i,2,usev){
if(!np[i]) p[++ptot]=i;
for(int j=1;j<=ptot && i*(ll)p[j]<=usev;++j){
res=i*p[j],np[res]=1;
if(!(i%p[j])) break;
}
}
return;
}
杜教筛
-
杜教筛可以在亚线性时间内求出某些特征积性函数的前缀和。
-
所谓某些特征积性函数,指的是存在一个对应的积性函数 \(g\),使得 \(f\ast g=h\),且 \(g,h\) 的前缀和可以 \(O(1)\) 得知的函数 \(f\)。
-
对于这种函数,我们考虑化式子,记 \(F,G,H\) 分别为 \(f,g,h\) 的前缀和:
- 我们这里用到的数论函数肯定都是 \(f(1)=1\) 的啦(积性嘛),所以事实上上式代表着
-
注意到式右的 \(F\) 值我们是已知的(如果我们从小到大来筛的话),\(G,H\) 可以 \(O(1)\) 得知,则对 \(\lfloor\frac{n}{d}\rfloor\) 整除分块,单轮复杂度 \(O(\sqrt{n})\)。
-
当然事实上我们只关心一个点的 \(F(n)\),故该式一般递归求解。对复杂度积分,发现在线筛预处理 \(n^\frac{2}{3}\) 内的 \(f\) 值时取得最优复杂度为 \(\Theta(n^\frac{2}{3})\)。
对 \(\mu\)
- 注意到 \(\mu\ast I=\epsilon\),且 \(I,\epsilon\) 的前缀和都可以 \(O(1)\) 得知,于是是板子。
对 \(\varphi\)
-
注意到 \(\varphi\ast I=id\),且 \(I,id\) 的前缀和都可以 \(O(1)\) 得知,于是还是板子。
-
不过,事实上我们有更好的方法。将 \(\varphi(i)\) 展开,得到如下式子:
-
那么我们对这个东西整除分块,直接杜来查 \(\mu\) 的前缀和,就可以了。当然,常数未必更小,无非是杜 \(\mu\) 和杜 \(\varphi\) 的区别。另外也可以把这个式子拆一拆,把 \(2,1\) 拆出来什么的。需要指出的是这只在单点求值时可以接受,多点求值或两个函数的 \(n\) 不同时,用杜求 \(S\mu\) 的复杂度显然会炸。
-
给出对这两个函数的示范代码(这里为了卡常使用了自定义哈希函数,另外因为是早期作品,实现不太规范):
namespace du{
bool np[maxn23];
int p[maxn23],ptot;
int mu[maxn23],smu[maxn23];
il void init(){
np[1]=1,mu[1]=smu[1]=1; int res;
For(i,2,usen23){
if(!np[i]) p[++ptot]=i,mu[i]=-1;
smu[i]=smu[i-1]+mu[i];
for(int j=1;j<=ptot && i*(ll)p[j]<=usen23ll;++j){
res=i*p[j],np[res]=1;
if(i%p[j]) mu[res]=-mu[i];
else break;
}
}
return;
}
struct my_hash{
static ull splitmix64(ull x){
x+=0x9e3779b97f4a7c15;
x=(x^(x>>30))*0xbf58476d1ce4e5b9;
x=(x^(x>>27))*0x94d049bb133111eb;
return x^(x>>31);
}
size_t operator()(ull x) const {
static const ull FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x + FIXED_RANDOM);
}
};
u_map<int,int,my_hash> smu2;
il int getmu(int n){
if(n<=usen23) return smu[n];
if(smu2.find(n)!=smu2.end()) return smu2[n];
int ret=1,res;
for(int l=2,r;;l=r+1){
res=n/l,r=n/res;
ret=ret-(r-l+1)*getmu(res);
if(r==n) break;
}
return smu2[n]=ret;
}
il ll getphi(int n){
static ll ret; ret=0;
static int res;
for(int l=1,r;;l=r+1){
res=n/l,r=n/res;
ret=ret+(getmu(r)-getmu(l-1))*(res*(res+1ll)/2);
if(r==n) break;
}
return ret;
}
}
对 \(id\mu\)
- 直接瞪不太容易瞪出来合适的 \(g\),我们先假设 \(g\) 已经找到,把卷积展开:
- 意识到 \(f=id\mu\),而
- 考虑类似的手法,令 \(g=id\),于是
- 额...有点出乎意料。不过确实 \(g\) 和 \(h=\epsilon\) 都可以 \(O(1)\) 求前缀和,直接杜就好了。
对 \(id\varphi\)
- 故技重施,快进到
- 那么可以直接杜了。
P3768 简单的数学题
-
题意:求 \(\sum\limits_{i=1}^n \sum\limits_{j=1}^n ij(i,j)\pmod{p}\)。
-
数据范围:\(n\leqslant 10^{10}\),\(4s\)。
-
考虑莫反:
-
...真(炎国粗口)长。
-
事实上,我们可以用欧拉反演:
-
就是这么简洁。
-
问题变成求 \(f(n)=n^2\varphi(n)\) 的前缀和 \(F\),考虑怎么杜:取 \(g=id_2\),有
- 那么 \(g=id_2,h=id_3\),问题解决。应当指出 \(G(n)=\frac{n(n+1)(2n+1)}{6},H(n)=ID^2(n)\),建议背过。推一下杜的式子:
- 所以有
- 终于结束了。呼~
23.1.12 T3 Product
-
来源:某 ACM。
-
题意:求 \(\prod\limits_{i=1}^n \prod\limits_{j=1}^n \prod\limits_{t=1}^n a^{(i,j)[t\mid (i,j)]} \pmod p\)。不保证 \(p\) 是质数。
-
数据范围:\(n\leqslant 10^9,m,p\) 在 int 范围内。
-
首先我们来一波令人窒息的化式子:
- 为了美观和能看清,我们记 \(g(n)=id(n)d(n),h(n)=\sum\limits_{i=1}^n \sum\limits_{j=1}^n [(i,j)=1]\),于是上式变成 \(f(n)=\prod\limits_{t=1}^n a^{g(t)h(\lfloor\frac{n}{t}\rfloor)}\)。我们继续推这个 \(h\):
-
好的...所以现在如果我们能快速求得 \(g\) 的前缀和,那么这就是一个整除分块套整除分块。
-
打表容易发现整除分块套整除分块,或者说 \(\sum\limits_{i=1}^n \sum\limits_{j=1}^{\lfloor\frac{n}{i}\rfloor} \lfloor\frac{n}{ij}\rfloor\) 这个式子在不合并 \(ij\)(莫反常用技巧,可惜这里用不了因为...因为什么都不允许这么合并啊,求和号和求积号显然不能这样简单合并,而且还是人家的指数)的时候,直接做是 \(O(n^{\frac{3}{4}})\) 的复杂度,可以接受。
-
好像有积分证明,但注意这里的 \(O\) 真的很 \(O\),实际的量大概要多个四五倍常数,更不要提整除分块本身的常数,所以 \(10^{10}\) 大概没机会。。
-
考虑求 \(G\)——不用试了!哈哈。\(G\) 不可杜,至少我没找到合适的配凑函数;也不可 PN,因为其质数拟合函数找不到(\(2id\) 显然不是积性函数)。难道没有机会了吗?不——只要运用一个另辟蹊径的小技巧。事实上,我们有
- 类似地,可以发现
-
没想到吧!Kono 整除分块的深刻应用哒!那么 \(g,h\) 的求解都是一个内层的整除分块。复杂度为 \(O(n^{\frac{3}{4}})\),足够通过本题。哦当然 \(\mu\) 是要杜的。
-
然而实现之后发现 T 了。脑测+机测发现算 \(h\) 的复杂度太高,仔细思考,意识到我们认为杜后的查询复杂度为 \(O(1)\)。这在单层整除分块中是成立的,因为在杜一次 \(S\mu(n)\) 的过程中已经求得,总计 \(O(n^\frac{2}{3})\)。然而这里我们是双层整除分块,所求的点有很多是不曾求出的,炸了。
-
究其原因...式子太丑。重新化一遍 \(h\) 试试:
-
那么这个式子符合杜的形式,筛出 \(S\varphi(n)\) 的过程中我们已经求出所有需求的点值。
-
另外按道理讲有一波令人窒息的扩展欧拉定理的分讨,需要记录“是否模过”。不过这里出题人偷懒只造了 \(p=10^9+7\) 的数据,所以我也偷个懒。
darkbzOJ #4916. 神犇和蒟蒻
-
题意:求 \(\sum\limits_{i=1}^n \mu(i^2),\varphi(i^2)\)。
-
被诈骗了!被狠狠地诈骗了!
-
第一问的诈骗很明显,由 \(\mu\) 的定义就可以看出。答案总为 \(1\)。
-
第二问的话...其实只要瞪一下 \(\varphi\) 的积性展开式就可以看出,\(\varphi(i^2)=i\varphi(i)\)。
-
典,这不就 \(id\varphi\) 嘛。直接杜之,\(O(n^\frac{2}{3})\)。
Powerful Number 筛
-
下简称 PN 筛。PN 筛可以在亚线性时间内求出某些特征积性函数的前缀和。
-
所谓某些特征积性函数,指的是存在一个对应的积性函数 \(g\),使得 \(g(p)=f(p)\),且求 \(g\) 的前缀和不比杜更复杂。称 \(g\) 为 \(f\) 的质数拟合函数。
-
对于此种函数,我们构造一个函数 \(h\) 使得 \(f=g\ast h\)。\(h\) 的存在性显然,因为 \(h\) 实际上是 \(g^{-1}\ast f\),而逆元存在性我们在狄利克雷卷积处已经证过了。从这里也能得出 \(h\) 的积性。
-
将上面的卷积式展开,得到:
- 这里 \(g\) 的前缀和我们可以 \(O(1)\) 得知,于是只要解决 \(h\) 的问题,就可以使用和杜教筛一样的整除分块办法来求解。关键是 \(h\) 要满足一个性质:其只在 Powerful Nubmer 处有值。
Powerful Number
-
定义:记正整数 \(n\) 的质数-指数向量为 \(n=\prod\limits_{i=1}^{m} p_i^{k_i}\)。\(n\) 是 Powerful Number(下简称 PN),当且仅当 \(\forall i\in [1,m],k_i>1\)。
-
性质 1:所有 PN 都可以写成 \(a^2b^3\) 的形式。
-
对于 \(k_i\bmod 2=0\),将 \(p_i^{k_i}\) 合并入 \(a^2\);
-
对于 \(k_i\bmod 2=1\),先将 \(p_i^3\) 合并入 \(b^3\),再将剩余的 \(p_i^{k_i-3}\) 合并入 \(a^2\)。
-
-
性质 2:\([1,n]\) 中的 PN 只有 \(O(\sqrt{n})\) 个。
- 考虑枚举 \(a\),于是 PN 的数量为
\[\sum\limits_{a=1}^{\sqrt{n}} \sqrt[3]{\frac{n}{a^2}}=O(\sqrt{n}) \]-
这个等号是积分求得的,我暂时不会。
-
怎么求出?线筛出 \(\sqrt{n}\) 内的所有质数,然后搜所有质数的指数即可。
-
性质 3:对任何的 \(f=g\ast h\And f(p)=g(p)\),\(h\) 函数仅在 PN 处可能有值。
- 考虑将卷积式在 \(p\) 处展开:
\[\begin{aligned} f(p) & =g(1)h(p)+g(p)h(1) \\ & =h(p)+g(p) \end{aligned} \]-
故显然 \(h(p)=0\)。因为 \(h\) 是积性函数,则 h 在非 PN 处一定无值。
-
注意 \(1\) 是 PN,因为 \(1\) 的质因数分解形式中 \(m=0\)。
-
于是我们可以在筛 PN 的过程中直接计算 \(F\),唯一的问题是 \(h\) 在 PN 处的值到底是多少。两种求法:
-
推出 \(h(p^k)\) 的公式,或者说,\(h(p^k)\) 仅关于 \(p\) 和 \(k\) 的表达式。比较难,但有一个很优的复杂度为 \(O(\dfrac{\sqrt{n}}{\ln \sqrt{n}})\)。
-
由 \(f=g\ast h\),展开得
\[\begin{aligned} f(p^k) & =\sum\limits_{i=0}^kg (p^i)h(p^{k-i}) \\ & =h(p^k)+\sum\limits_{i=1}^k g(p^i)h(p^{k-i}) \end{aligned} \]- 于是有
\[h(p^k)=f(p^k)-\sum\limits_{i=1}^k g(p^i)h(p^{k-i}) \]- 枚举 \(p,k\) 计算之,若已知 \(f(p^k)\),则复杂度为
\[O(\frac{\sqrt{n}}{\ln \sqrt{n}})\times \log^2 n\approx O(\sqrt{n}\log) \]-
分别是质数个数、枚举 \(k\)、枚举 \(i\)。事实上,这是一个很宽的上界,远远不满(较大质数的次数很小,远不及一般意义上的 \(\log_2 n\))。
-
后来在 Min_25 筛的板题中的实验也表明,纯暴力求(理论复杂度 \(O(\sqrt{n}\log)\))和公式法(理论复杂度 \(O(\dfrac{\sqrt{n}}{\log})\) 的用时一样,可见瓶颈似乎不在这部分(该题需要杜)。
-
是的,按道理讲这里我们还没有筛出 \(f(p^k)\),但积性函数就是通过给出 \(f(p^k)\) 的取值定义的,故我们在题面中就知道了。注意这里的 \(g\) 也应当用这种方式求,如果杜的话...感觉上显然复杂度不对,试了一下也没跑出来。
-
-
然后搜索所有 PN 直接计算 \(F\),式子为
-
若 \(G\) 可以 \(O(1)\) 得知,则复杂度为搜索的 \(O(\sqrt{n})\)。
-
若 \(G\) 需要杜,则复杂度为 \(O(n^\frac{2}{3})\),因若先杜一次 \(G(n)\),则访问到的所有 \(G\) 都是已经筛出的值。
-
理由:对于 \(i=1\),\(G\) 是 \(G(n)\);否则 \(G\) 在杜中是求 \(G(n)\) 时被减去的项,整除分块过了。
-
综上,两部分复杂度加合,可以认为 PN 筛的复杂度上界是 \(O(n^\frac{2}{3})\)。
-
最后,再度给出便于实现的一般过程:
-
构造 \(f\) 的质数拟合函数 \(g\)。
-
构造快速计算 \(G\) 的方法。
-
计算 \(h(p^k)\),通常是预处理把表打出来,这样第 4 步直接查表就好。
-
搜索所有 PN,过程中累加答案。
-
-
给出示范代码,这里 \(f(p^k)=p^k(p^k-1)\),即 P5325 【模板】Min_25 筛(\(h\) 部分使用纯暴力计算):
namespace du{
//f=iphi(i),g=id,h=id_2
const int maxn23=4641593,maxpn23=325167;
const ll usen23=4641590;
ll n23;
bool np[maxn23];
ll p[maxpn23]; int ptot;
ll f[maxn23],sf[maxn23];
il void init(){
n23=ceil(pow(n,2/(long double)3));
np[1]=1,f[1]=sf[1]=1; ll res;
for(ll i=2;i<=n23;++i){
if(!np[i]) p[++ptot]=i,f[i]=i*(i-1)%mod;
sf[i]=addr(sf[i-1],f[i]);//printf("phi[%lld]:%lld f[%lld]:%lld sf[%lld]:%lld\n",i,phi[i],i,f[i],i,sf[i]);
for(int j=1;i*p[j]<=n23 && j<=ptot;++j){
res=i*p[j],np[res]=1;
if(i%p[j]) f[res]=f[i]*p[j]%mod*(p[j]-1)%mod;
else{
f[res]=f[i]*p[j]%mod*p[j]%mod;
break;
}
}
}
return;
}
il ll getid(ll nnow){nnow%=mod; return nnow*(nnow+1)%mod*inv[2]%mod;}
il ll getid2(ll nnow){nnow%=mod; return nnow*(nnow+1)%mod*(2*nnow+1)%mod*inv[6]%mod;}
struct my_hash{
static ull splitmix64(ull x){
x+=0x9e3779b97f4a7c15;
x=(x^(x>>30))*0xbf58476d1ce4e5b9;
x=(x^(x>>27))*0x94d049bb133111eb;
return x^(x>>31);
}
size_t operator()(ull x) const {
static const ull FIXED_RANDOM = chrono::steady_clock::now().time_since_epoch().count();
return splitmix64(x + FIXED_RANDOM);
}
};
u_map<ll,ll,my_hash> sf2;
il ll getf(ll nnow){
if(nnow<=n23) return sf[nnow];
if(sf2.find(nnow)!=sf2.end()) return sf2[nnow];
ll ret=getid2(nnow); static ll res;
for(ll l=2,r;;l=r+1){
res=nnow/l,r=nnow/res;
minu(ret,minur(getid(r),getid(l-1))*getf(res)%mod);
if(r==nnow) break;
}
return sf2[nnow]=ret;
}
}
namespace pn{
//f(p^k)=p^k(p^k-1),g=idphi,h=? 这里 g 是 du::f
const int maxpsq=9597;
ll nsq; int uptot;
vector<ll> h[maxpsq];//h[i][j] 表示 h(p[i]^j) 的值
il void init(){//计算 h(p^k)。
//注意到至少 h(p^2) 处才有值,故只需要考虑 sqrt(n) 以内的质数
nsq=ceil(sqrt(n)); if(n==1) return;//这是因为 n=1 时 p[1]=0,但我不想在循环里加判断
ll powp[36]={1}; int top=0;//2^35>10^10,足够了
for(int i=1;du::p[i]<=nsq;++i){
while(powp[top]*du::p[i]<=n) ++top,powp[top]=powp[top-1]*du::p[i];
h[i].resize(top+1),h[i][0]=1,h[i][1]=0;
For(k,2,top){
powp[k]%=mod;
h[i][k]=powp[k]*(powp[k]-1)%mod;
For(j,1,k) minu(h[i][k],(powp[j]*minur(powp[j],powp[j-1])%mod)*h[i][k-j]%mod);
}
top=0,uptot=i;
}
return;
}
ll nnow,lim[maxpsq],sum;
il void dfs(int now,ll vn,ll hn){//考虑到第几个质数,当前的值是多少,当前的 h 是多少
if(now>uptot || vn*du::p[now]>=lim[now]){//计算贡献。这个剪枝非常重要。
sum=(sum+hn*du::getf(nnow/vn))%mod;
return;
}
dfs(now+1,vn,hn);//k 为 0
int k=2;
for(ll res=vn*du::p[now];res<lim[now];res=res*du::p[now],++k) dfs(now+1,res*du::p[now],hn*h[now][k]%mod);
return;
}
il ll getf(ll nnowin){
nnow=nnowin;
For(i,1,uptot) lim[i]=nnow/du::p[i]+1;//这是一个防止 10^10 相乘爆 ll 的 trick,记 v*p=x,则 x*p<=nnow->x<=nnow/p->x<floor(nnow/p)+1
sum=0;
dfs(1,1,1);
return sum;
}
}
P5325 【模板】Min_25 筛
-
题意:求 \(F(n)\),这里积性函数 \(f(p^k)=p^k(p^k-1)\)。
-
数据范围:\(n\leqslant 10^{10}\)。
-
首先我们令 \(k=1\) 看看 \(f\) 的表现,发现 \(f(p)=p(p-1)\),注意到 \(\varphi(p)=p-1\),故考虑构造 \(f\) 的质数拟合函数 \(g=id\varphi\)。
-
嗯?这个函数好像似曾相识。是的,这个函数可以杜,参看杜教筛的最后一个例子。
-
于是是板子,\(O(\sqrt{n}\log+n^\frac{2}{3})\)。
-
但本题有特殊的一点在于 \(h\) 是可以求公式的,思路类似杜,过程如下:
-
(未完成)
-
又 \(h(p)=0\),于是通过累加法可以得到 \(h(p^k)=(k-1)(p-1)p^k\)。
HDU7186 Jo loves counting
-
题意:
-
对一个数 \(n\),我们定义 \(good_n\) 为所有质数-指数向量的质数部分和 \(n\) 的质数-指数向量的质数部分相同的数构成的集合。
-
形式化地,记 \(n=\prod\limits_{i=1}^m p_i^{k_i}\),则 \(i=\prod\limits_{i=1}^m p_i^{k_i}\And k_i\geqslant 0 \leftrightarrow i\in good_n\)。
-
定义函数 \(f(n)=\frac{n}{|good_n|}\),求 \(\frac{1}{m}F(m)\),对 \(4179340454199820289\)(一个大质数,\(=29\times 2^{57}+1\))取模。
-
-
数据范围:\(T\leqslant 12,m\leqslant 10^{12}\),但 \(m>10^6\) 的 \(T\leqslant 6\)。
-
我们先考察一下 \(|good_n|\),发现其其实是 \(\prod\limits_{i=1}^m k_i\),因为 \(good\) 中元素的质数-指数向量中,\(p_i\) 的次数可以取 \(1\sim k_i\)。
-
于是发现 \(f\) 在 \(p\) 处的表现为 \(f(p)=\frac{p}{1}=p\),考虑取 \(id\) 为其质数拟合函数。
-
后面就是板子咯,而且还不用杜,\(O(T\sqrt{m}\log)\)。还没写...