# 题目链接

https://lydsy.com/JudgeOnline/problem.php?id=4815

# 题解

$f(a,b)=\frac{a}{\gcd(a,b)}\frac{b}{\gcd(a,b)}f(\gcd(a,b),\gcd(a,b))$

$\sum_{g=1}^{k}f(g,g)\sum_{i=1}^{\lfloor k/g\rfloor} \sum_{j=1}^{\lfloor k/g\rfloor}ij[\gcd(i,j)=1]$

$\sum_{g=1}^k f(g,g)\sum_{i=1}^{\lfloor k/g\rfloor}i^2\varphi(i)$

# 代码

#include <cmath>
#include <cstdio>
#include <algorithm>

template<typename T>
{
T x=0;
int f=1;
char ch=getchar();
while((ch<'0')||(ch>'9'))
{
if(ch=='-')
{
f=-f;
}
ch=getchar();
}
while((ch>='0')&&(ch<='9'))
{
x=x*10+ch-'0';
ch=getchar();
}
return x*f;
}

const int maxn=4000000;
const int mod=1000000007;

int p[maxn+10],prime[maxn+10],cnt,phi[maxn+10],inv[maxn+10];

int getprime()
{
p[1]=phi[1]=1;
for(int i=2; i<=maxn; ++i)
{
if(!p[i])
{
prime[++cnt]=i;
phi[i]=i-1;
}
for(int j=1; (j<=cnt)&&(i*prime[j]<=maxn); ++j)
{
int x=i*prime[j];
p[x]=1;
if(i%prime[j]==0)
{
phi[x]=phi[i]*prime[j];
break;
}
phi[x]=phi[i]*(prime[j]-1);
}
}
for(int i=1; i<=maxn; ++i)
{
phi[i]=(phi[i-1]+1ll*phi[i]*i%mod*i)%mod;
}
inv[0]=inv[1]=1;
for(int i=2; i<=maxn; ++i)
{
inv[i]=1ll*(mod-mod/i)*inv[mod%i]%mod;
}
return 0;
}

int block[maxn+10],sum[maxn+10],bs[maxn+10],val[maxn+10],totb;

int modify(int pos,int v)
{
int k=v-val[pos];
if(k<0)
{
k+=mod;
}
val[pos]+=k;
if(val[pos]>=mod)
{
val[pos]-=mod;
}
for(int i=pos; block[i]==block[pos]; ++i)
{
sum[i]+=k;
if(sum[i]>=mod)
{
sum[i]-=mod;
}
}
for(int i=block[pos]; i<=totb; ++i)
{
bs[i]+=k;
if(bs[i]>=mod)
{
bs[i]-=mod;
}
}
return 0;
}

int getsum(int l,int r)
{
(l==0)&&(++l);
int v=bs[block[r]-1]-bs[block[l]-1]+sum[r];
if(v<0)
{
v+=mod;
}
if(v>=mod)
{
v-=mod;
}
if(block[l-1]==block[l])
{
v-=sum[l-1];
if(v<0)
{
v+=mod;
}
}
return v;
}

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

int quickpow(int a,int b)
{
int res=1;
while(b)
{
if(b&1)
{
res=1ll*res*a%mod;
}
a=1ll*a*a%mod;
b>>=1;
}
return res;
}

int m,n,a,b,k;
long long x;

int main()
{
getprime();
totb=sqrt(n);
for(int i=1; i<=n; ++i)
{
block[i]=(i-1)/totb+1;
}
totb=(n-1)/totb+1;
for(int i=1; i<=n; ++i)
{
val[i]=1ll*i*i%mod;
bs[block[i]]+=val[i];
if(bs[block[i]]>=mod)
{
bs[block[i]]-=mod;
}
}
for(int i=1; i<=totb; ++i)
{
bs[i]+=bs[i-1];
if(bs[i]>=mod)
{
bs[i]-=mod;
}
}
for(int i=1; i<=n; ++i)
{
sum[i]=((block[i]==block[i-1])?sum[i-1]:0)+val[i];
if(sum[i]>=mod)
{
sum[i]-=mod;
}
}
for(int i=1; i<=m; ++i)
{