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} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{d} \rfloor} ij\ \ [gcd(i,j)==1]\)
\(Ans=\sum_{d=1}^nd^3\sum_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\ k^2 \sum_{i=1}^{\lfloor \frac{n}{kd} \rfloor} \sum_{j=1}^{\lfloor \frac{n}{kd} \rfloor} ij\)
\(Ans=\sum_{T=1}^n \sum_{k=1}^{T} \mu(k) \ k^2\ (\frac{T}{k})^3 \ Sum(\lfloor \frac{n}{T} \rfloor)^2 \)
\(Ans=\sum_{T=1}^n T^2 \phi(T) \ Sum(\lfloor \frac{n}{T} \rfloor)^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;
}