筛法
概述
- 筛法原本是筛取质数的一种算法。在 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)\)。还没写...
