【luogu3768】简单的数学题【杜教筛】【欧拉函数】

题目传送门
题意
求 $(\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j))\ mod\ p$，$n<=10^{10}$，$p$为质数。
题解
之前一直想着用$\mu$函数搞，结果死都只能搞到$O(n)$。后来，我fa♂现了bzoj3518点组计数，这道题用了欧拉函数的一个性质$n=\sum_{i|n}\varphi(i)$，好像叫做欧拉反演？
于是新日暮里的大门打开了。
于是开始推式子。
$\sum_{i=1}^n\sum_{j=1}^nijgcd(i,j)$
$=\sum_{i=1}^n\sum_{j=1}^nij\sum_{k|gcd(i,j)}\varphi(k)$
$=\sum_{k=1}^nk^2\varphi(k)\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{k}\rfloor}ij$
我们令$S(n)=\sum_{i=1}^ni=\frac{n(n+1)}{2}$
则原式
$=\sum_{k=1}^nk^2\varphi(k)S(\lfloor\frac{n}{k}\rfloor)^2$
于是我们只需要能够快速求出$\sum_{k=1}^nk^2\varphi(k)$就可以分块求和计算了。
这不就是杜教筛干的活吗？
$\sum_{i=1}^ni^2\sum_{j|i}\varphi(j)$
$=\sum_{i=1}^ni^3$
$=S(n)^2$
$=\sum_{j=1}^n\sum_{i=1}^{\lfloor\frac{n}{j}\rfloor}i^2j^2\varphi(j)$
因为$i≤\frac{n}{j}$
所以$j≤\frac{n}{i}$
因此原式
$=\sum_{i=1}^ni^2\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2\varphi(j)$
把$i=1$带人
得$\sum_{i=1}^ni^2\varphi(i)=S(n)^2-\sum_{i=2}^ni^2\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}j^2\varphi(j)$
然后就递归计算即可。细节详见代码。
这道题就解决了！
代码

#include<cstdio>
typedef long long ll;
const ll N=5000005,M=1000017;
ll mod,n,ans,p[N],phi[N];
bool vis[N];
struct hashmap{
    int cnt,head[M],nxt[M*10];
    ll to[M*10],v[M*10];
    bool count(ll f){
        ll x=f%M;
        for(int i=head[x];i;i=nxt[i]){
            if(to[i]==f){
                return true;
            }
        }
        return false;
    }
    ll get(ll f){
        ll x=f%M;
        for(int i=head[x];i;i=nxt[i]){
            if(to[i]==f){
                return v[i];
            }
        }
        return 0;
    }
    void set(ll f,ll val){
        ll x=f%M;
        to[++cnt]=f;
        nxt[cnt]=head[x];
        v[cnt]=val;
        head[x]=cnt;
    }
}mp;
inline ll sqr(ll x){
    x%=mod;
    return x*x%mod;
}
ll calc(ll x){
    if(!x){
        return 0;
    }
    x%=mod;
    return sqr((1+x)*x/2)%mod;
}
ll get(ll n){
    n%=mod;
    if(n%3!=1){
        return n*(n+1)/6%mod*(2*n+1)%mod;
    }else if(n&1){
        return (n+1)*(2*n+1)/6%mod*n%mod;
    }else{
        return n*(2*n+1)/6%mod*(n+1)%mod;
    }
}
ll solve(ll n){
    if(n<=5000000){
        return phi[n];
    }
    if(mp.count(n)){
        return mp.get(n);
    }
    ll res=calc(n);
    for(ll i=2,last;i<=n;i=last+1){
        last=n/(n/i);
        res-=solve(n/i)*(get(last)-get(i-1))%mod;
        res%=mod;
    }
    mp.set(n,res);
    return res;
}
int main(){
    scanf("%lld%lld",&mod,&n);
    phi[1]=1;
    for(int i=2;i<=5000000;i++){
        if(!vis[i]){
            p[++p[0]]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=p[0]&&i*p[j]<=5000000;j++){
            vis[i*p[j]]=true;
            if(i%p[j]){
                phi[i*p[j]]=phi[i]*(p[j]-1);
            }else{
                phi[i*p[j]]=phi[i]*p[j];
            }
        }
    }
    for(int i=1;i<=5000000;i++){
        phi[i]=phi[i]*i%mod*i%mod;
        phi[i]+=phi[i-1];
        phi[i]%=mod;
    }
    for(ll i=1,last;i<=n;i=last+1){
        last=n/(n/i);
        ans+=(solve(last)-solve(i-1))*calc(n/i)%mod;
        ans%=mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}

最后附上蒟蒻博主龙hou飞yan凤wu舞chi的公式草稿一张2333 鼠标写字真是太方便啦