P3768 简单的数学题 题解

Statement

\(\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j),n\leq 1e10\)

Solution

\[\begin{align*} &\sum_{i=1}^n\sum_{j=1}^nij\gcd(i,j)\\ =&\sum_{x=1}^nx^3\sum_{i=1}^{\frac nx}\sum_{j=1}^{\frac nx}ij[i\perp j]\\ =&\sum_{x=1}^nx^3\sum_{i=1}^{\frac nx}\sum_{j=1}^{\frac nx}ij\sum_{d|gcd(i,j)}\mu(d)\\ =&\sum_{x=1}^nx^3\sum_{d=1}^{\frac nx}\mu(d)\sum_{i=1}^{\frac n{dx}}\sum_{j=1}^{\frac n{dx}}ij\\ =&\sum_{x=1}^nx^3\sum_{d=1}^{\frac nx}\mu(d)g(\frac n{dx})\\ =&\sum_{s=1}^ns^2g(\frac ns)\sum_{d|s}\frac sd\mu(d)\\ =&\sum_{s=1}^ns^2g(\frac ns)\varphi(s)\\ \end{align*} \]

考虑杜教筛

设 \(f(x)=x^2\varphi(x),g(x)=x^2\)

则 \((f*g)(n)=\sum_{d|n}f(d)g(\frac nd)=n^2\sum_{d|n}\varphi(d)=n^3\)

有 \(S(n)g(1)=\sum_{i=1}^n i^3-\sum_{i=2}^nS(\frac ni)g(i)\)

Code

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N = 5e6+5;

int read(){
    int s=0,w=1;char ch=getchar();
    while(ch<'0'||ch>'9'){if(ch=='-')w=-1;ch=getchar();}
    while(ch>='0'&&ch<='9')s=s*10+ch-'0',ch=getchar();
    return s*w;
}

int prime[N],phi[N];
int n,p,cnt,ans,inv2,inv6;
map<int,int>res;
bool vis[N];

void getprime(int n){
    phi[1]=1;
    for(int i=2;i<n;++i){
        if(!vis[i])prime[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*prime[j]<n;++j){
            vis[i*prime[j]]=true;
            if(i%prime[j]==0){phi[i*prime[j]]=phi[i]*prime[j];break;}
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(int i=1;i<n;++i)phi[i]=(phi[i-1]+phi[i]*i%p*i%p)%p;
}
int ksm(int a,int b){
    int res=1;
    while(b){
        if(b&1)res=res*a%p;
        a=a*a%p,b>>=1;
    }
    return res;
}
int pw(int n){return n*n%p;}
int csum(int n){return n%=p,pw(n*(n+1)%p*inv2%p);}
int ssum(int n){return n%=p,n*(n+1)%p*(2*n+1)%p*inv6%p;}
int calc(int n){
    if(n<N)return phi[n];
    if(res[n])return res[n];
    int ans=csum(n);
    for(int i=2,j;i<=n;i=j+1)
        j=n/(n/i),ans=(ans-(ssum(j)-ssum(i-1)+p)*calc(n/i)%p+p)%p;
    return res[n]=ans;
}

signed main(){
    p=read(),n=read(),getprime(N);
    inv2=ksm(2,p-2),inv6=ksm(6,p-2);
    for(int i=1,j;i<=n;i=j+1)
        j=n/(n/i),(ans+=(calc(j)-calc(i-1)+p)%p*csum(n/i)%p)%=p;
    printf("%lld\n",ans);
    return 0;
}