P3768 【简单的数学题】

P3768 【简单的数学题】

\(Ans=\sum ^{n}_{i=1}\sum ^{n}_{j=1}ijgcd(i,j)\)

\(=\sum ^{n}_{i=1}\sum ^{n}_{j=1}ij\sum _{k|i,k|j} φ(k)\)

\(=\sum ^{n}_{k=1} φ(k) \sum _{k|i}\sum _{k|j}ij\)

\(=\sum ^{n}_{k=1}\varphi (k) k^{2} (\sum ^{n/k}_{i=1}i)^{2}\)

\(=\sum ^{n}_{k=1}\varphi (k) k^{2} \sum ^{n/k}_{i=1}i^{3}\)

\(n<=1e10\)杜教筛筛 \(\varphi (k) k^{2}\)

#include<bits/stdc++.h>
#include<tr1/unordered_map>
using namespace std;
typedef long long LL;
const int maxn=5000000+9;
inline LL Read(){
    LL x=0,f=1; char c=getchar();
    while(c<'0'||c>'9'){
        if(c=='-') f=-1; c=getchar();
    }
    while(c>='0'&&c<='9')
        x=(x<<3)+(x<<1)+c-'0',c=getchar();
    return x*f;
}
LL n,p,inv;
int phi[maxn],prime[maxn];
LL sum[maxn];
bool visit[maxn];
inline LL Pow(LL base,LL b){
    LL a=1;
    while(b){
        if(b&1)
            a=(a*base)%p;
        base=(base*base)%p;
        b>>=1;
    }
    return a;
}
inline void F_phi(int N){
    inv=Pow(6,p-2);
    phi[1]=1; int tot(0);
    for(int i=2;i<=N;++i){
        if(!visit[i]){
            phi[i]=i-1;
            prime[++tot]=i;
        }
        for(int j=1;j<=tot&&i*prime[j]<=N;++j){
            visit[i*prime[j]]=true;
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            else
                phi[i*prime[j]]=phi[i]*phi[prime[j]];
        }
    }
    for(int i=1;i<=N;++i)
        sum[i]=(sum[i-1]+1ll*i*i%p*phi[i]%p)%p;
    //for(int i=1;i<=2000;++i)
    //    printf("%lld ",sum[i]);printf("\n");
}
tr1::unordered_map<LL,LL> w;
inline LL S2(LL x){
    x%=p;
    return x*(x+1)%p*(2*x+1)%p*inv%p;
}
LL S3(LL x){
    x%=p;
    return (x*(x+1)/2)%p*((x*(x+1)/2)%p)%p;
}
LL Calc(LL now){
    if(now<=5000000)
        return sum[now];
    if(w[now])
        return w[now];
    LL num=S3(now);
    for(LL l=2,r;l<=now;l=r+1){
        r=now/(now/l);
        num=(num-Calc(now/l)*(S2(r)-S2(l-1)+p)%p+p)%p;
    }
    return w[now]=num;
}
int main(){
    p=Read(),n=Read();
    F_phi(5000000);

    LL ans(0);
    for(LL l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=(ans+S3(n/l)*((Calc(r)-Calc(l-1)+p)%p))%p;
    }
    printf("%lld",ans);
    return 0;
}/*
998244353 2000
883968974
*/

原文地址：https://www.cnblogs.com/y2823774827y/p/10229183.html