杜教筛学习笔记

杜教筛是用于在低于线性时间内求解一些积性函数前缀和的工具。

原理&模板

假设现在要求积性函数(f)的前缀和，设(sum_{i=1}^{n} f(i)=S(n))

再找一个积性函数，计算它们卷积的前缀和：

[egin{alignat}{1} &sum_{i=1}^{n}(f * g)(i) \ =&sum_{i=1}^{n} sum_{d mid i} f(d) gleft(frac{i}{d} ight) \ =&sum_{d=1}^{n} g(d) sum_{i=1}^{leftlfloorfrac{n}{d} ight floor} f(i) \ =&sum_{d=1}^{n} g(d) Sleft(leftlfloorfrac{n}{d} ight floor ight) end{alignat} ]

将（4）的第一项和其他项分离，可以得到

[egin{alignat}{1} g(1) S(n)&=sum_{i=1}^{n}(f * g)(i)-sum_{i=2}^{n} g(i) Sleft(leftlfloorfrac{n}{i} ight floor ight) \ S(n)&=frac{sum_{i=1}^{n}(f * g)(i)-sum_{i=2}^{n} g(i) Sleft(leftlfloorfrac{n}{i} ight floor ight)}{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
",getphi(n),getmu(n));
    }
}

P3768杜教筛+整除分块

最后的式子(d^2varphi(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
",getsum(n));
}