# Description

(1)对任意的正整数a,b，都要满足f(a,b)=f(b,a)；
(2)对任意的正整数a,b，都要满足b×f(a,a+b),(a+b)×f(a,b)。

# Solution

$ans=\sum _{d=1}^{n}f\left(d\right)\sum _{i=1}^{n}\sum _{j=1}^{n}\left[gcd\left(i,j\right)=d\right]i×j$

$ans=\sum _{d=1}^{n}f\left(d\right)×{d}^{2}\sum _{i=1}^{⌊\frac{n}{d}⌋}\sum _{j=1}^{⌊\frac{n}{d}⌋}\left[gcd\left(i,j\right)=1\right]i×j$

$ans=\sum _{d=1}^{n}f\left(d\right)×{d}^{2}\sum _{i=1}^{⌊\frac{n}{d}⌋}\sum _{j=1}^{⌊\frac{n}{d}⌋}ij\sum _{x|gcd\left(i,j}\mu \left(x\right)$

$ans=\sum _{d=1}^{n}f\left(d\right)×{d}^{2}\sum _{x=1}^{⌊\frac{n}{d}⌋}\mu \left(x\right){x}^{2}×S\left(⌊\frac{n}{dx}⌋{\right)}^{2}$

$g\left(n\right)=\sum _{i=1}^{n}\mu \left(i\right){i}^{2}×S\left(⌊\frac{n}{i}⌋{\right)}^{2}$$g(n)=\sum_{i=1}^{n}\mu(i)i^2\times S(\lfloor\frac{n}{i}\rfloor)^2$，可以发现g(n)是能递推的，具体说来就是$g\left(n\right)=g\left(n-1\right)+{n}^{3}\sum _{i|n}\frac{\mu \left(i\right)}{i}$$g(n)=g(n-1)+n^3\sum_{i|n}\frac{\mu(i)}{i}$

# Code

#include <stdio.h>
#include <algorithm>
#include <math.h>
#define rep(i,st,ed) for (register int i=st;i<=ed;++i)
#define drp(i,st,ed) for (register int i=st;i>=ed;--i)

typedef long long LL;
const int MOD=1000000007;
const int N=4000010;
const int Q=2010;

LL inv[N+5],f[N+5],g[N+5],h[N+5];
LL t[Q][Q],sum[Q];
int st[Q],ed[Q],bel[N],prime[N/10],n,m;

bool not_prime[N+5];

void mod(LL &x) {
x-=x>=MOD?MOD:0;
}

LL x=0,v=1; char ch=getchar();
for (;ch<'0'||ch>'9';v=(ch=='-')?(-1):(v),ch=getchar());
for (;ch>='0'&&ch<='9';x=x*10+ch-'0',ch=getchar());
return x*v;
}

inline int gcd(int x,int y) {
return (!y)?(x):gcd(y,x%y);
}

void pre_work() {
inv[1]=f[1]=g[1]=h[1]=1;
rep(i,2,n) {
f[i]=1;
inv[i]=(MOD-MOD/i)*inv[MOD%i]%MOD;
if (!not_prime[i]) {
prime[++prime[0]]=i;
h[i]=(1-inv[i]+MOD)%MOD;
}
for (register int j=1,x=i*prime[j];i*prime[j]<=n&&j<=prime[0];x=i*prime[++j]) {
not_prime[x]=1;
if (i%prime[j]==0) {
h[x]=h[i];
break;
}
h[x]=(h[i]*h[prime[j]])%MOD;
}
mod(g[i]=g[i-1]+((LL)i%MOD*(LL)i%MOD*(LL)i%MOD*h[i]%MOD)%MOD);
}
int size=(int)ceil(sqrt(n));
rep(i,1,n) {
bel[i]=(i-1)/size+1;
if (!st[bel[i]]) st[bel[i]]=i;
ed[bel[i]]=std:: max(ed[bel[i]],i);
}
rep(b,1,bel[n]) {
rep(i,st[b],ed[b]) {
t[b][i-st[b]+1]=(t[b][i-st[b]]+f[i]*i%MOD*i%MOD)%MOD;
}
}
drp(b,bel[n],1) {
sum[b]=(sum[b+1]+t[b][ed[b]-st[b]+1])%MOD;
}
}

inline void modify(int x) {
int bx=bel[x];
rep(i,x,ed[bx]) t[bx][i-st[bx]+1]=(t[bx][i-st[bx]]+f[i]*(LL)i%MOD*(LL)i%MOD)%MOD;
drp(i,bx,1) sum[i]=(sum[i+1]+t[i][ed[i]-st[i]+1])%MOD;
}

inline LL query(int x,int y) {
int bx=bel[x],by=bel[y];
if (bx==by) return (t[bx][y-st[bx]+1]%MOD-t[bx][x-st[bx]]%MOD+MOD)%MOD;
LL ret=0;
mod(ret+=t[bx][ed[bx]-st[bx]+1]-t[bx][x-1-st[bx]+1]+MOD);
mod(ret+=t[by][y-st[by]+1]);
mod(ret+=sum[bx+1]-sum[by]+MOD);
return ret;
}

inline void solve(int n) {
LL ans=0;
for (register int i=1,j;i<=n;i=j+1) {
j=n/(n/i);
mod(ans+=query(i,j)*g[n/i]%MOD);
}
printf("%lld\n", ans);
}

int main(void) {
}