# [luoguP3768]简单的数学题

### 传送门

$f(d)=d\sum_{i=1}^{n}\sum_{j=1}^{n}ij[gcd(i,j)==d]\\$

$f(d)=d^3\sum_{T=1}^{n/d}\mu(T)T^2\sum_{i=1}^{n/Td}\sum_{j=1}^{n/Td}ij\\$

$ans=\sum_{d=1}^{n}f(d)=\sum_{d=1}^{n}d^3\sum_{T=1}^{n/d}\mu(T)T^2\sum_{i=1}^{n/Td}\sum_{j=1}^{n/Td}ij\\$

$ans=\sum_{k=1}^{n}\sum_{i=1}^{n/k}\sum_{j=1}^{n/k}ij\sum_{d|k}\mu(\frac{k}{d})(\frac{k}{d})^2d^3\\ ans=\sum_{k=1}^{n}\sum_{i=1}^{n/k}\sum_{j=1}^{n/k}ijk^2\sum_{d|k}\mu(\frac{k}{d})d\\$

$\mu*id=\varphi\\ g(k)=k^2\varphi(k)$

$S(n)h(1)=\sum_{i=1}^{n}(h*g)-\sum_{i=2}^{n}S(\lfloor\frac{n}{i}\rfloor)h(i)$

$h*g=\sum_{d|n}h(\lfloor\frac{n}{d}\rfloor)g(d)=\sum_{d|n}h(\lfloor\frac{n}{d}\rfloor)d^2\varphi(d)\\$

$\sum_{d|n}\varphi(d)=n$

$h*g=\sum_{d|n}h(\lfloor\frac{n}{d}\rfloor)d^2n\\$

$h*g=n^3\\ S(n)=\sum_{i=1}^{n}i^3-\sum_{i=2}^{n}i^2S(\lfloor\frac{n}{i}\rfloor)$

$ans=\sum_{k=1}^{n}\sum_{i=1}^{n/k}\sum_{j=1}^{n/k}ijk^2\sum_{d|k}\mu(\frac{k}{d})d\\$

$ans=\sum_{k=1}^{n}s(\lfloor\frac{n}{k}\rfloor)g(k)$

#include<cstdio>
#include<algorithm>
#include<iostream>
#include<cstring>
#include<tr1/unordered_map>
using namespace std;
{
char ch;bool ok;
for(ok=0,ch=getchar();!isdigit(ch);ch=getchar())if(ch=='-')ok=1;
for(x=0;isdigit(ch);x=x*10+ch-'0',ch=getchar());if(ok)x=-x;
}
#define rg register
#define ll long long
const int maxn=3e6+10;
int mod,ans,phi[maxn],pri[maxn],tot,inv2,inv6,f[maxn];
ll n;bool vis[maxn];
tr1::unordered_map<ll,int>mp;
int mul(int x,int y){return 1ll*x*y-1ll*x*y/mod*mod;}
int del(int x,int y){return x-y<0?x-y+mod:x-y;}
int mi(int a,int b){
int ans=1;
while(b){
if(b&1)ans=mul(ans,a);
b>>=1,a=mul(a,a);
}
return ans;
}
void prepare()
{
phi[1]=1;
for(rg int i=2;i<=3e6;i++){
if(!vis[i])pri[++tot]=i,phi[i]=i-1;
for(rg int j=1;j<=tot&&pri[j]*i<=3e6;j++){
vis[pri[j]*i]=1;
if(!(i%pri[j])){phi[pri[j]*i]=mul(pri[j],phi[i]);break;}
else phi[pri[j]*i]=mul(phi[pri[j]],phi[i]);
}
}
for(rg int i=1;i<=3e6;i++)f[i]=mul(mul(i,i),phi[i]);
}
int s(int n) {return mul(mul(mul(n*2+1,n+1),n),inv6);}
int square(ll l,ll r){return del(s(r%mod),s((l-1)%mod));}
int sum(ll n){return mul(mul((n+1)%mod,n%mod),inv2);}
int sum1(ll n){return mul(mul(n%mod,(n+1)%mod),inv2);}
int get_phi(ll n){
if(n<=3e6)return f[n];
if(mp[n])return mp[n];
int ans=mul(sum1(n),sum1(n));
for(rg ll i=2,j;i<=n;i=j+1)j=n/(n/i),ans=del(ans,mul(square(i,j),get_phi(n/i)));
return mp[n]=ans;
}
int main()
{