[P3768]简单的数学题

Description：

求出((sum_{i=1}^n sum_{j=1}^n ij gcd (i,j)) mod p)

Hint:

(n<=10^{10}​)

Solution:

(Ans=sum_{d=1}^nd^3 sum_{i=1}^{lfloor frac{n}{d} floor} sum_{j=1}^{lfloor frac{n}{d} floor} ij [gcd(i,j)==1]​)

(Ans=sum_{d=1}^nd^3sum_{k=1}^{lfloorfrac{n}{d} floor}mu(k) k^2 sum_{i=1}^{lfloor frac{n}{kd} floor} sum_{j=1}^{lfloor frac{n}{kd} floor} ij​)

(Ans=sum_{T=1}^n sum_{k=1}^{T} mu(k) k^2 (frac{T}{k})^3 Sum(lfloor frac{n}{T} floor)^2 ​)

(Ans=sum_{T=1}^n T^2 phi(T) Sum(lfloor frac{n}{T} floor)^2 ​)

杜教筛出 (T^2 phi(T)) 的前缀和

(令g(x)=x^2)

(sum_{d=1}^n f(d)*g(frac{n}{d}) = sum phi(d) d^2(frac{n}{d})^2=n^3)

至此可求

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mxn=8e6;
ll mod,tot,y,z;
int p[mxn+5],vis[mxn+5];
ll ph[mxn+5];
map<ll ,ll > sph;

ll qpow(ll a,ll b) 
{
    ll ans=1,base=a;
    while(b) {
        if(b&1) ans=1ll*ans*base%mod;
        base=1ll*base*base%mod;
        b>>=1;
    }
    return ans;
}

void sieve() 
{
    vis[1]=ph[1]=1;
    for(int i=2;i<=mxn;++i) {
        if(!vis[i]) ph[i]=i-1,p[++tot]=i;
        for(int j=1;j<=tot&&p[j]*i<=mxn;++j) {
            vis[p[j]*i]=1;
            if(i%p[j]) ph[p[j]*i]=ph[p[j]]*ph[i]%mod;
            else {ph[p[j]*i]=ph[i]*p[j]%mod;break;}
        }
    }
    for(int i=1;i<=mxn;++i) ph[i]=(ph[i]*i%mod*i%mod+ph[i-1])%mod;
    y=qpow(2,mod-2),z=qpow(6,mod-2);
}

inline ll cal1(ll x) {
    x%=mod;
    return (1ll*x*(x+1)%mod*y%mod)*(1ll*x*(x+1)%mod*y%mod)%mod;
}

inline ll cal2(ll x) {
    x%=mod;
    return 1ll*x*(x+1)%mod*(2*x%mod+1)%mod*z%mod;
}

ll get(ll n)
{
    if(n<=mxn) return ph[n]; 
    if(sph[n]) return sph[n]; ll ans=0;
    for(ll l=2,r;l<=n;l=r+1) {
        r=n/(n/l);
        ans=(ans+(cal2(r)-cal2(l-1)+mod)%mod*get(n/l)%mod)%mod;
    }
    return sph[n]=((cal1(n)-ans)%mod+mod)%mod;
}

int main()
{
    ll n; ll ans=0;
    scanf("%d %lld",&mod,&n); sieve();
    for(ll l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        ans=(ans+1ll*cal1(n/l)%mod*(get(r)-get(l-1)+mod)%mod)%mod;
    }
    printf("%lld",ans);
    return 0;
}