[模板]杜教筛

用途

 比线性更快($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$即可

#include<bits/stdc++.h>
#include<tr1/unordered_map>
#define CLR(a,x) memset(a,x,sizeof(a))
#define MP make_pair
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pa;
const int maxn=1e5+10,P=1e9+7,inv6=166666668;

inline ll rd(){
    ll x=0;char c=getchar();int neg=1;
    while(c<'0'||c>'9'){if(c=='-') neg=-1;c=getchar();}
    while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
    return x*neg;
}

int N,phi[maxn],pri[maxn],cnt;
bool np[maxn];
tr1::unordered_map<int,int> sum;

inline int calc(int n){
    if(sum.count(n)) return sum[n];
    ll re=0;
    re=1ll*n*(n+1)%P*(2*n+1)%P*inv6%P;
    for(int i=2,j;i<=n;i=j+1){
        j=n/(n/i);
        re=(re-1ll*(j-i+1)*(j+i)/2%P*calc(n/i))%P;
    }
    sum[n]=re;
    return re;
}

int main(){
    //freopen("","r",stdin);
    int i,j,k;
    phi[1]=1;
    for(i=2;i<=1e5;i++){
        if(!np[i]){
            pri[++cnt]=i;
            phi[i]=i-1;
        }for(j=1;j<=cnt&&pri[j]*i<=1e5;j++){
            np[i*pri[j]]=1;
            if(i%pri[j]==0){
                phi[i*pri[j]]=phi[i]*pri[j];
                break;
            }
            phi[i*pri[j]]=phi[i]*phi[pri[j]];
        }
    }
    sum[0]=0;
    for(i=1;i<=1e5;i++){
        sum[i]=(sum[i-1]+1ll*i*phi[i])%P;
    }
    printf("1
%d
",(calc(rd())+P)%P);
    return 0;
}