杜教筛
参考资料：
https://blog.csdn.net/Ike940067893/article/details/84781307
https://www.luogu.com.cn/problemnew/solution/P4213
前置知识：
主要内容：
杜教筛
杜教筛能用低于线性的时间复杂度求积性函数前缀和。
比如要求积性函数 (f) 的前缀和 (sum(n)=sumlimits_{i=1}^nf_i)。
需要找一个辅助积性函数 (g)，考虑 ((f*g))（(*) 是狄利克雷卷积，不懂自学）的前缀和 (sum')：
[egin{split}
sum'(n)=&sumlimits_{i=1}^n(f*g)(i)\
=&sumlimits_{i=1}^nsumlimits_{d|i}f(frac id)g(d)\
=&sumlimits_{d=1}^ng(d)sumlimits_{i=1}^{lfloorfrac nd
floor}f(i)\
=&sumlimits_{d=1}^ng(d)color{#77cc33}{sum(lfloorfrac nd
floor)}\
end{split}
]
[egin{split}
herefore g(1)sum(n)=&sumlimits_{d=1}^ng(d)sum(lfloorfrac nd
floor)-sumlimits_{d=2}^ng(d)sum(lfloorfrac nd
floor)\
=&sumlimits_{i=1}^n(f*g)(i)-sumlimits_{d=2}^ng(d)sum(lfloorfrac nd
floor)\
end{split}
]
于是得到了最主要的式子，单独提出来：
[g(1)sum(n)=sumlimits_{i=1}^n(f*g)(i)-sumlimits_{d=2}^ng(d)sum(lfloorfrac nd
floor)
]
如果 (g) 找得合适，大部分时候 (g(1)=1)，而 (sumlimits_{i=1}^n(f*g)(i)) 和 (sumlimits_{i=1}^ng(i)) 可以 (Theta(1)) 算，而 (sumlimits_{d=2}^ng(d)sum(lfloorfrac nd floor)) 可以整除分块算。
于是直接递归求 (sum(n)) 的时间复杂度为 (Theta(n^{frac{3}{4}}))。
证明
设求出 $sum(n)$ 的时间复杂度为 $Theta(T(n))$：[egin{split}
T(n)=&sqrt n+sumlimits_{i=2}^{sqrt n}left(T(i)cdot T(lfloorfrac{n}{i}
floor)
ight)\
ge&sqrt n+sumlimits_{i=2}^{sqrt n}left(sqrt icdot sqrt{lfloorfrac{n}{i}
floor)}
ight)\
ge&sqrt n+sumlimits_{i=2}^{sqrt n}2sqrt{sqrt n}\
=&n^{frac{3}{4}}\
end{split}
]
虽然等式步步是大于，但是如果记忆化 (sum(n))，时间复杂度就是 (Theta(n^{frac{3}{4}}))。
优化：
线性筛求出前 (n^{frac{2}{3}}) 个 (sum(i))，然后再杜教筛剩下的，优化后的时间复杂度为 (Theta(n^{frac{2}{3}}))。
证明
设求出 $sum(n)$ 的时间复杂度为 $Theta(T(n))$：[egin{split}
T(n)=&m+sumlimits_{i=2}^{lfloorfrac nm
floor}sqrt{lfloorfrac ni
floor}\
=&m+frac{n}{sqrt m}\
ge&2n^{frac{2}{3}}(=: operatorname{while} m=n^{frac{2}{3}})\
end{split}
]
经典例题
【模板】杜教筛（Sum）
(T) 组测试数据，给定 (n)，求 (sumlimits_{i=1}^nvarphi(i)) 和 (sumlimits_{i=1}^nmu(i))。
数据范围：(1le Tle 10)，(1le nle 2^{31}-1)。
求 (mu) 的前缀和：
此时 (f=mu)，如果取 (g=1)，那么 (f*g=mu*1=epsilon)。
同时满足：
- (g(1)=1)。
- (sumlimits_{i=1}^n(f*g)(i)=1) 可以 (Theta(1)) 算出。
- (sumlimits_{i=1}^ng(i)=n) 可以 (Theta(1)) 算出。
所以有：
[sum(n)=1-sumlimits_{d=2}^nsum(lfloorfrac nd
floor)
]
右边分块递归求解，(sum) 记忆化，时间复杂度 (Theta(n^{frac{2}{3}}))。
( exttt{Code})
hash<int,lng> Mu;
il lng DuMu(re int n){
if(n<=N) return mu[n];
if(Mu[n]) return Mu[n];
re lng res=1ll;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuMu(n/l);
return Mu[n]=res;
}
求 (varphi) 的前缀和：
此时 (f=varphi)，如果取 (g=1)，那么 (f*g=varphi*1=ID)。
同时满足：
- (g(1)=1)。
- (sumlimits_{i=1}^n(f*g)(i)=frac{n(n+1)}{2}) 可以 (Theta(1)) 算出。
- (sumlimits_{i=1}^ng(i)=n) 可以 (Theta(1)) 算出。
所以有：
[sum(n)=frac{n(n+1)}{2}-sumlimits_{d=2}^nsum(lfloorfrac nd
floor)
]
右边分块递归求解，(sum) 记忆化，时间复杂度 (Theta(n^{frac{2}{3}}))。
( exttt{Code})
hash<int,lng> Phi;
il lng DuPhi(re int n){
if(n<=N) return phi[n];
if(Phi[n]) return Phi[n];
re lng res=1ll*n*(n+1)/2;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuPhi(n/l);
return Phi[n]=res;
}
$ exttt{Code}$
#include <bits/stdc++.h>
using namespace std;
//&Start
#define inf 0x3f3f3f3f
#define re register
#define il inline
#define hash unordered_map
typedef long long lng;
typedef unsigned long long ulng;
typedef vector<int> veci;
//&Data
#define N 5000000
//&Sieve
bitset<N+10> np;
int p[N+10];
lng mu[N+10],phi[N+10];
il void Sieve(){
np[1]=true;
mu[1]=1,phi[1]=1;
re int cnt=0;
for(re int i=2;i<=N;i++){
if(!np[i]) p[++cnt]=i,mu[i]=-1,phi[i]=i-1;
for(re int j=1;j<=cnt&&i*p[j]<=N;j++){
np[i*p[j]]=1;
if(i%p[j]==0){mu[i*p[j]]=0,phi[i*p[j]]=phi[i]*p[j];break;}
mu[i*p[j]]=-mu[i],phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
for(re int i=2;i<=N;i++) mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
//&Du-Sieve
hash<int,lng> Mu,Phi; //记忆化sum
il lng DuPhi(re int n){
if(n<=N) return phi[n];
if(Phi[n]) return Phi[n];
re lng res=1ll*n*(n+1)/2;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuPhi(n/l);
return Phi[n]=res;
}
il lng DuMu(re int n){
if(n<=N) return mu[n];
if(Mu[n]) return Mu[n];
re lng res=1ll;
for(re int l=2,r;l<=n;l=r+1)
r=n/(n/l),res-=(r-l+1)*DuMu(n/l);
return Mu[n]=res;
}
//&Main
int main(){
Sieve();
re int n,t;
scanf("%d",&t);
for(re int i=1;i<=t;i++){
scanf("%d",&n);
printf("%lld %lld
",DuPhi(n),DuMu(n));
}
return 0;
}
毒瘤的数学题
给定 (n) 和 (p)，求
[left(sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j) ight)mod p ]
数据范围：(1le nle color{red}{10^{10}})，(5 imes10^{8}le ple 1.1 imes 10^{9})。
稍微改了一下题目抬头，因为对这题的毒瘤深有体会。如果你看不到这道题的前半段推导，说明你没有学好莫比乌斯反演。
第一部分：普通の推导
[egin{split}
F(n)=&sumlimits_{i=1}^nsumlimits_{j=1}^nijgcd(i,j)\
=&sumlimits_{d=1}^ndsumlimits_{i=1}^nisumlimits_{j=1}^nj[gcd(i,j)=d]\
=&sumlimits_{d=1}^ndsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}idsumlimits_{j=1}^{lfloorfrac{n}{d}
floor}jd[gcd(i,j)=1]\
=&sumlimits_{d=1}^ndsumlimits_{i=1}^{lfloorfrac{n}{d}
floor}idsumlimits_{j=1}^{lfloorfrac{n}{d}
floor}jdsumlimits_{k|gcd(i,j)}mu(k)\
=&sumlimits_{d=1}^nd^3sumlimits_{k=1}^{lfloorfrac{n}{d}
floor}mu(k)sumlimits_{i=1}^{lfloorfrac{n}{d}
floor}[k|i]isumlimits_{j=1}^{lfloorfrac{n}{d}
floor}[k|j]j\
=&sumlimits_{d=1}^nd^3sumlimits_{k=1}^{lfloorfrac{n}{d}
floor}mu(k)sumlimits_{i=1}^{lfloorfrac{n}{dk}
floor}iksumlimits_{j=1}^{lfloorfrac{n}{dk}
floor}jk\
=&sumlimits_{d=1}^nd^3sumlimits_{k=1}^{lfloorfrac{n}{d}
floor}mu(k)k^2sumlimits_{i=1}^{lfloorfrac{n}{dk}
floor}isumlimits_{j=1}^{lfloorfrac{n}{dk}
floor}j\
end{split}
]
设 (S(n)=sumlimits_{i=1}^ni=frac{n(n+1)}{2})，所以
[F(n)=sumlimits_{d=1}^nd^3sumlimits_{k=1}^{lfloorfrac{n}{d}
floor}mu(k)k^2S(lfloorfrac{n}{dk}
floor)^2
]
第二部分：代换
将 (T=dk) 带入：
[egin{split}
F(n)=&sumlimits_{d=1}^nd^3sumlimits_{k=1}^{lfloorfrac{n}{d}
floor}mu(k)k^2S(lfloorfrac{n}{dk}
floor)^2\
=&sumlimits_{d=1}^nd^3sumlimits_{k=1}^{lfloorfrac{n}{d}
floor}mu(k)k^2S(lfloorfrac{n}{T}
floor)^2\
=&sumlimits_{T=1}^nS(lfloorfrac{n}{T}
floor)^2sumlimits_{d|T}d^3mu(frac Td)left(frac Td
ight)^2\
=&sumlimits_{T=1}^nS(lfloorfrac{n}{T}
floor)^2T^2sumlimits_{d|T}dmu(frac Td)\
=&sumlimits_{T=1}^nS(lfloorfrac{n}{T}
floor)^2T^2(mu*ID)(T)\
end{split}
]
将公式 (mu*ID=varphi) 带入得：
[F(n)=sumlimits_{T=1}^nS(lfloorfrac{n}{T}
floor)^2T^2varphi(T)
]
第三部分：杜教卷积
这里给出其中一种方法，可能不是最优的：
将整个 (T^2varphi(T)) 提出来杜教，即令 (f(n)=n^2varphi(n))。
为了防止遗忘，把最重要的套路式再拿出来遛一下：
[g(1)sum(n)=sumlimits_{i=1}^n(f*g)(i)-sumlimits_{d=2}^ng(d)sum(lfloorfrac nd
floor)
]
为了抵消掉 (f) 中的 (n^2)，(g) 应该等于 (n^2)。
则有：
[egin{split}
(f*g)(n)=&sumlimits_{d|n}f(d)g(frac nd)\
=&sumlimits_{d|n}d^2varphi(d)cdot left(frac nd
ight)^2\
=&n^2sumlimits_{d|n}varphi(d)\
end{split}
]
将公式 (varphi*1=ID) 带入得：
[(f*g)(n)=n^3
]
于是便满足：
- (g(1)=1)。
- (sumlimits_{i=1}^n(f*g)(i)=S(n)^2) 可以 (Theta(1)) 算出。
- (sumlimits_{i=1}^ng(i)=frac{n(n+1)(2n+1)}{6}) 可以 (Theta(1)) 算出。
所以可以杜教筛得：
[sum(n)=sumlimits_{i=1}^ni^2varphi(i)=S(n)^2-sumlimits_{d=2}^nd^2sum(lfloorfrac{n}{d}
floor)
]
然后再看原式：
[F(n)=sumlimits_{T=1}^nS(lfloorfrac{n}{T}
floor)^2T^2varphi(T)
]
就可以分块加杜教解决了，复杂度 (Theta(n^{frac {2}{3}}))。
代码：要取模，要 ( exttt{long long})。
$ exttt{Code}$
#include <bits/stdc++.h>
using namespace std;
//&Start
#define inf 0x3f3f3f3f
#define re register
#define il inline
#define hash unordered_map
typedef long long lng;
typedef unsigned long long ulng;
typedef vector<int> veci;
//&Data
#define N 5000000
int m;
//&Pow
il int Pow(re int a,re int x){ //用于求逆元
re int res=1;
for(;x;a=1ll*a*a%m,x>>=1)if(x&1) res=1ll*res*a%m;
return res;
}
int inv;
//&Sum
il int sum(re lng x){x%=m;return 1ll*x*(x+1)/2%m;} //即S
il int sum2(re lng x){x%=m;return 1ll*x*(x+1)%m*(2*x+1)%m*inv%m;} //即n(n+1)(2n+1)/6
//&Sieve
bitset<N+10> np;
int p[N+10],phi[N+10],f[N+10];
il void Sieve(){ //优化筛
np[1]=true,f[1]=phi[1]=1;
re int cnt=0;
for(re int i=2;i<=N;i++){
if(!np[i]) p[++cnt]=i,phi[i]=i-1;
f[i]=1ll*i*i%m*phi[i]%m;
for(re int j=1;j<=cnt&&i*p[j]<=N;j++){
np[i*p[j]]=1;
if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
phi[i*p[j]]=phi[i]*phi[p[j]];
}
}
for(re int i=2;i<=N;i++) (f[i]+=f[i-1])%=m;
}
//&Du-Sieve0-------------->杜教筛
hash<lng,int> F;
il int DuF(re lng n){
if(n<=N) return f[n];
if(F[n]) return F[n];
re int res=sum(n); res=1ll*res*res%m;
for(re lng l=2,r;l<=n;l=r+1)
r=n/(n/l),(res-=1ll*(sum2(r)-sum2(l-1))%m*DuF(n/l)%m)%=m;
return F[n]=(res+m)%m;
}
//&Main
int main(){
re lng n;
scanf("%d%lld",&m,&n);
inv=Pow(6,m-2),Sieve();
re int ans=0,tmp;
for(re lng l=1,r;l<=n;l=r+1){
r=n/(n/l),tmp=sum(n/l);
(ans+=1ll*(DuF(r)-DuF(l-1))%m*tmp%m*tmp%m)%=m; //分块加杜教
}
printf("%d
",(ans+m)%m);
return 0;
}
如果有感悟，会更新。
祝大家学习愉快！