[模板]杜教筛

Posted 2021-02-01 ressed

用途

 比线性更快($O(n^{frac{2}{3}})$)地求积性函数的前缀和

前置知识：狄利克雷卷积

形如$h(n)=sumlimits_{d|n}f(d)g(frac{n}{d})$，则称$h(n)=f(x)*g(x)$

如果f和g都是积性函数，则卷出的h也是积性函数

可以证明，狄利克雷卷积满足交换律、结合律、分配律

比较重要的卷积式子（抄的..）：

$$mu*1=varepsilon , varepsilon(n)=[n=1]$$

$$varphi*1=id , id(n)=n $$

$$id*mu=varphi$$

$$id*1=sigma , sigma表示因数和$$

 

做法

设要求的是$sum{f(i)}$

构造积性函数$h,g$，使得$h=g*f$,而且需要容易求出$h$和$g$的前缀和

推式子：设S(n)是f的前缀和（以下不是整除的都向下取整）

$$sumlimits_{i=1}^{n}h(i)=sumlimits_{i=1}^{n}sumlimits_{d|i}g(d)f(frac{n}{d})$$

$$sumlimits_{i=1}^{n}h(i)=sumlimits_{d=1}^{n}g(d)sumlimits_{i=1}^{frac{n}{d}}f(i)$$

$$sumlimits_{i=1}^{n}h(i)=sumlimits_{d=1}^{n}g(d)S(frac{n}{d})$$

$$把后面的第一项提出来，sumlimits_{i=1}^{n}h(i)=g(1)S(n)+sumlimits_{d=2}^{n}g(d)S(frac{n}{d})$$

$$g(1)S(n)=sumlimits_{i=1}^{n}h(i)-sumlimits_{d=2}^{n}g(d)S(frac{n}{d})$$

于是就可以整除分块，然后递归地求S了（需要先手动求出比较小范围内的S）

需要记忆化，不想带log的话可以用unordered_map（c++11或者tr1/unordered_map）或者手写hash

而且一开始求的那些S不要放到map里，单开一个数组存，不然常数会很大

 

所以主要是构造比较难构造，看着上面的卷积式瞎怼吧...

例题

bzoj4916 神犇和蒟蒻

注意到$varphi (i^2) = ivarphi(i)$

然后构造$g=id , h=id^2$即可

 1 #include<bits/stdc++.h>
 2 #include<tr1/unordered_map>
 3 #define CLR(a,x) memset(a,x,sizeof(a))
 4 #define MP make_pair
 5 using namespace std;
 6 typedef long long ll;
 7 typedef unsigned long long ull;
 8 typedef pair<int,int> pa;
 9 const int maxn=1e5+10,P=1e9+7,inv6=166666668;
10 
11 inline ll rd(){
12     ll x=0;char c=getchar();int neg=1;
13     while(c<‘0‘||c>‘9‘){if(c==‘-‘) neg=-1;c=getchar();}
14     while(c>=‘0‘&&c<=‘9‘) x=x*10+c-‘0‘,c=getchar();
15     return x*neg;
16 }
17 
18 int N,phi[maxn],pri[maxn],cnt;
19 bool np[maxn];
20 tr1::unordered_map<int,int> sum;
21 
22 inline int calc(int n){
23     if(sum.count(n)) return sum[n];
24     ll re=0;
25     re=1ll*n*(n+1)%P*(2*n+1)%P*inv6%P;
26     for(int i=2,j;i<=n;i=j+1){
27         j=n/(n/i);
28         re=(re-1ll*(j-i+1)*(j+i)/2%P*calc(n/i))%P;
29     }
30     sum[n]=re;
31     return re;
32 }
33 
34 int main(){
35     //freopen("","r",stdin);
36     int i,j,k;
37     phi[1]=1;
38     for(i=2;i<=1e5;i++){
39         if(!np[i]){
40             pri[++cnt]=i;
41             phi[i]=i-1;
42         }for(j=1;j<=cnt&&pri[j]*i<=1e5;j++){
43             np[i*pri[j]]=1;
44             if(i%pri[j]==0){
45                 phi[i*pri[j]]=phi[i]*pri[j];
46                 break;
47             }
48             phi[i*pri[j]]=phi[i]*phi[pri[j]];
49         }
50     }
51     sum[0]=0;
52     for(i=1;i<=1e5;i++){
53         sum[i]=(sum[i-1]+1ll*i*phi[i])%P;
54     }
55     printf("1
%d
",(calc(rd())+P)%P);
56     return 0;
57 }