杜教筛学习笔记
杜教筛是用于在低于线性时间内求解一些积性函数前缀和的工具。
原理&模板
假设现在要求积性函数\(f\)的前缀和，设\(\sum_{i=1}^{n} f(i)=S(n)\)
再找一个积性函数，计算它们卷积的前缀和：
\[\begin{alignat}{1}
&\sum_{i=1}^{n}(f * g)(i) \\
=&\sum_{i=1}^{n} \sum_{d \mid i} f(d) g\left(\frac{i}{d}\right) \\
=&\sum_{d=1}^{n} g(d) \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} f(i) \\
=&\sum_{d=1}^{n} g(d) S\left(\left\lfloor\frac{n}{d}\right\rfloor\right)
\end{alignat}
\]
将（4）的第一项和其他项分离，可以得到
\[\begin{alignat}{1}
g(1) S(n)&=\sum_{i=1}^{n}(f * g)(i)-\sum_{i=2}^{n} g(i) S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\
S(n)&=\frac{\sum_{i=1}^{n}(f * g)(i)-\sum_{i=2}^{n} g(i) S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)}{g(1)}
\end{alignat}
\]
此时若\(f*g\)的前缀和可以快速求出，则\(S(n)\)就可以利用这个公式进行递归+整除分块来求解。
一般的代码框架如下：
ll GetSum(int n) { // 算 f 前缀和的函数
ll ans = f_g_sum(n); // 算 f * g 的前缀和
// 以下这个 for 循环是数论分块
for(ll l = 2, r; l <= n; l = r + 1) { // 注意从 2 开始
r = (n / (n / l));
ans -= (g_sum(r) - g_sum(l - 1)) * GetSum(n / l);
// g_sum 是 g 的前缀和
// 递归 GetSum 求解
} return ans; //实际求解时需要利用unordered_map进行记忆化
}
这个代码的时间复杂度时\(O(n^\frac{3}{4})\)
但是如果用线性筛先将一部分前\(n^\frac{2}{3}\)的答案筛出，之后再用杜教筛，此时的时间复杂度就可以达到\(O(n^\frac{2}{3})\)
例题
P4213求欧拉函数和莫比乌斯函数的前缀和
//给定一个数n，求出1~n的欧拉函数前缀和，莫比乌斯函数前缀和
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=5e6+5;
bool isPrime[maxn];
int Prime[maxn], primecnt = 0;
ll phi[maxn],mu[maxn];
void getPrime(int n){
memset(isPrime, 1, sizeof(isPrime));
isPrime[1] = 0;
phi[1] = mu[1] = 1;
for(int i = 2; i <= n; i++){
if(isPrime[i]){
Prime[++primecnt] = i;
phi[i]=i-1;
mu[i]=-1;
}
for(int j = 1; j <= primecnt && i*Prime[j] <= n; j++) {
isPrime[i*Prime[j]] = 0;
if(i % Prime[j] == 0){
phi[i*Prime[j]]=phi[i]*Prime[j];
mu[i*Prime[j]]=0;
break;
}
else{
phi[i*Prime[j]]=phi[i]*(Prime[j]-1);
mu[i*Prime[j]]=-1*mu[i];
}
}
}
}
void init(int n){
for(int i=1;i<=n;i++){
phi[i]+=phi[i-1];
mu[i]+=mu[i-1];
}
}
unordered_map<int,int>ansmu;
unordered_map<int,ll>ansphi;
ll getphi(ll x){
if(x<maxn-5)return phi[x];
if(ansphi[x])return ansphi[x];
ll ans=(1+x)*x/2;
for(ll l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-(l-1))*getphi(x/l);
}
return ansphi[x]=ans;
}
ll getmu(ll x){
if(x<maxn-5)return mu[x];
if(ansmu[x])return ansmu[x];
ll ans=1;
for(ll l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-(l-1))*getmu(x/l);
}
return ansmu[x]=ans;
}
int main () {
getPrime(maxn-5);
init(maxn-5);
int T;
scanf("%d",&T);
while(T--){
ll n;
scanf("%lld",&n);
printf("%lld %lld\n",getphi(n),getmu(n));
}
}
P3768杜教筛+整除分块
最后的式子\(d^2\varphi(d)\)可以用杜教筛求，后面部分的可以用整除分块，杜教筛套一层数论分块，复杂度还是\(O(n^\frac{2}{3})\)
#include <bits/stdc++.h>
#define stdd(x) (x>=mod?x-mod:x)
using namespace std;
typedef long long ll;
const int maxn=5e6+5;
bool ntp[maxn];
int Prime[maxn], primecnt = 0;
ll phi[maxn],pre[maxn];
ll mod;
ll _2,_6;
void getPrime(int n){
ntp[1] = 1;
phi[1] = 1;
for(int i = 2; i <= n; i++){
if(!ntp[i]){
Prime[++primecnt] = i;
phi[i]=i-1;
}
for(int j = 1; j <= primecnt && i*Prime[j] <= n; j++) {
ntp[i*Prime[j]] = 1;
if(i % Prime[j] == 0){
phi[i*Prime[j]]=phi[i]*Prime[j];
break;
}
else{
phi[i*Prime[j]]=phi[i]*(Prime[j]-1);
}
}
}
pre[0]=0;
for(int i=1;i<=n;i++){
pre[i]=((ll)pre[i-1]+(ll)phi[i]*i%mod*i%mod)%mod;
}
}
ll fp(ll b,ll p){
ll ans=1;
while(p){
if(p&1){
ans=ans*b%mod;
}
p>>=1;
b=b*b%mod;
}
return ans;
}
ll S2(ll n){
n%=mod;
return n*(n+1)%mod*(2*n+1)%mod*_6%mod;
}
ll S3(ll n){
n%=mod;
return (n*(n+1)/2%mod)*(n*(n+1)/2%mod)%mod;
}
unordered_map<ll,ll>ans;
ll getans(ll x){
if(x<maxn-5)return pre[x];
if(ans.count(x))return ans[x];
ll res=S3(x);
for(ll l=2,r;l<=x;l=r+1){
r=x/(x/l);
res-=(S2(r)-S2(l-1)+mod)%mod*getans(x/l)%mod;
res=stdd(res+mod);
}
return ans[x]=res;
}
ll getsum(ll x){
ll res=0;
for(ll l=1,r;l<=x;l=r+1){
r=x/(x/l);
res+=S3(x/l)*(getans(r)-getans(l-1)+mod)%mod;
res=stdd(res);
}
return res;
}
int main () {
ll n;
scanf("%lld%lld",&mod,&n);
getPrime(min(n,(ll)maxn-5));
_2=fp(2,mod-2);
_6=fp(6,mod-2);
printf("%lld\n",getsum(n));
}
