杜教筛是用于在低于线性时间内求解一些积性函数前缀和的工具。
原理&模板
假设现在要求积性函数(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));
}