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

题目传送门
题意
求 (∑ni=1∑nj=1ijgcd(i,j)) mod p$(\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j))modp$，n<=1010$n<={10}^{10}$，p$p$为质数。
题解
之前一直想着用μ$\mu$函数搞，结果死都只能搞到O(n)$O(n)$。后来，我fa♂现了bzoj3518点组计数，这道题用了欧拉函数的一个性质n=∑i|nφ(i)$n=\sum _{i|n}\phi (i)$，好像叫做欧拉反演？
于是新日暮里的大门打开了。
于是开始推式子。
∑ni=1∑nj=1ijgcd(i,j)$\sum _{i=1}^n\sum _{j=1}^{n}ijgcd(i,j)$
=∑ni=1∑nj=1ij∑k|gcd(i,j)φ(k)$=\sum _{i=1}^n\sum _{j=1}^{n}ij\sum _{k|gcd(i,j)}\phi (k)$
=∑nk=1k2φ(k)∑⌊nk⌋i=1∑⌊nk⌋j=1ij$=\sum _{k=1}^n{k}^2\phi (k)\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{k}\rfloor }ij$
我们令S(n)=∑ni=1i=n(n+1)2$S(n)=\sum _{i=1}^{n}i=\frac{n(n+1)}{2}$
则原式
=∑nk=1k2φ(k)S(⌊nk⌋)2$=\sum _{k=1}^n{k}^2\phi (k)S(\lfloor \frac{n}{k}\rfloor {)}^2$
于是我们只需要能够快速求出∑nk=1k2φ(k)$\sum _{k=1}^n{k}^2\phi (k)$就可以分块求和计算了。
这不就是杜教筛干的活吗？
∑ni=1i2∑j|iφ(j)$\sum _{i=1}^n{i}^2\sum _{j|i}\phi (j)$
=∑ni=1i3$=\sum _{i=1}^n{i}^3$
=S(n)2$=S(n{)}^2$
=∑nj=1∑⌊nj⌋i=1i2j2φ(j)$=\sum _{j=1}^n\sum _{i=1}^{\lfloor \frac{n}{j}\rfloor }{i}^2{j}^2\phi (j)$
因为i≤nj$i\le \frac{n}{j}$
所以j≤ni$j\le \frac{n}{i}$
因此原式
=∑ni=1i2∑⌊ni⌋j=1j2φ(j)$=\sum _{i=1}^n{i}^2\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }{j}^2\phi (j)$
把i=1$i=1$带人
得∑ni=1i2φ(i)=S(n)2−∑ni=2i2∑⌊ni⌋j=1j2φ(j)$\sum _{i=1}^n{i}^2\phi (i)=S(n{)}^2-\sum _{i=2}^n{i}^2\sum _{j=1}^{\lfloor \frac{n}{i}\rfloor }{j}^2\phi (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
",(ans+mod)%mod);
    return 0;
}

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