(Description)
求$$sum_{i=1}^nsum_{j=1}^ifrac{lcm(i,j)}{gcd(i,j)}$$
答案对(10^9+7)取模。
(n<=10^9)
(Solution)
以前做的反演题都是(j)枚举到(n)，但是现在(j)只枚举到(i)就非常难受，考虑怎么求(sum_{i=1}^nsum_{j=1}^nfrac{lcm(i,j)}{gcd(i,j)})。
可以把它看成是一个(n*n)的网格，第(i)行第(j)列上的数是(frac{lcm(i,j)}{gcd(i,j)})，需要我们求的是包括对角线在内的下三角矩阵的权值和。
所以答案为(所有网格权值之和+对角线上的权值和)/2。
[sum_{i=1}^nsum_{j=1}^nfrac{lcm(i,j)}{gcd(i,j)}
]
[=sum_{d=1}^nsum_{i=1}^nsum_{j=1}^nfrac{ij}{d^2}[gcd(i,j)==d]
]
[=sum_{d=1}^nsum_{i=1}^{n/d}sum_{j=1}^{n/d}ij[gcd(i,j)==1]
]
考虑怎么求后半部分
[sum_{i=1}^nsum_{j=1}^nij[gcd(i,j)==1]
]
[=sum_{i=1}^nsum_{j=1}^nijsum_{d|gcd(i,j)}mu_d
]
枚举(d)
[=sum_{d=1}^nmu_dsum_{d|i}^nsum_{d|j}^nij
]
[=sum_{d=1}^nmu_dd^2sum_{i=1}^{n/d}sum_{j=1}^{n/d}ij
]
令(sum(n)=sum_{i=1}^ni)，
所以原式
[=sum_{i=1}^nmu_ii^2sum(n/i)^2
]
带回到一开始的式子里去
[sum_{d=1}^nsum_{i=1}^{n/d}mu_ii^2sum(frac{n}{id})^2
]
按照套路令(T=id)
[=sum_{T=1}^nsum(n/T)^2sum_{d|T}mu_dd^2
]
令(f(x)=sum_{d|x}mu_dd^2)，现在如果我们可以快速的求出(f(x))的前缀和，那么就可以数论分块算答案了。
可是(f(x))并不是一个熟悉的数论函数，怎么才能用杜教筛呢？
可以把(f(x))写成几个函数的卷积的形式。
令(g(x)=mu_xx^2)。那么(f=g*1)。现在要找一个函数(h)使得(f*h=g*1*h)好算。我们知道(sum_{d|x}mu_d=e)，所以令(h(x)=x^2)来把(g(x)中的乘x^2)消掉。
所以就构造出了(s=f*h=g*1*h=e*1=1)，不难发现(f)是个积性函数，可以线筛。
#include<complex>
#include<cstdio>
#include<map>
using namespace std;
const int mod=1e9+7;
const int N=2e6+7;
int n,tot,inv2=mod+1>>1,inv6=166666668;
int prime[N],mu[N],f[N];
bool check[N];
map<int,int>mp;
int qread()
{
int x=0;
char ch=getchar();
while(ch<'0' || ch>'9')ch=getchar();
while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();}
return x;
}
void Init()
{
int nn=min(n,N-1);
check[1]=f[1]=1;
for(int i=2;i<=nn;i++)
{
if(!check[i])prime[++tot]=i,f[i]=1-1ll*i*i%mod;
for(int j=1;j<=tot && i*prime[j]<=nn;j++)
{
check[i*prime[j]]=1;
if(i%prime[j])f[i*prime[j]]=1ll*f[i]*f[prime[j]]%mod;
else
{
f[i*prime[j]]=f[i];
break;
}
}
}
for(int i=1;i<=nn;i++)
f[i]=(f[i-1]+f[i])%mod;
}
int Calc1(int x)
{
long long res=1ll*x*(x+1)/2%mod;
return res*res%mod;
}
int Calc2(int x)
{
return 1ll*x*(x+1)%mod*(x+x+1)%mod*inv6%mod;
}
int Sum(int x)
{
if(x<N)return f[x];
if(mp[x])return mp[x];
long long res=x;
for(int l=2,r;l<=x;l=r+1)
{
r=x/(x/l);
res=(res-1ll*(Calc2(r)-Calc2(l-1)+mod)*Sum(x/l))%mod;
}
return mp[x]=(res+mod)%mod;
}
int main()
{
scanf("%d",&n);
Init();
long long ans=0;
for(int l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+1ll*Calc1(n/l)*(Sum(r)-Sum(l-1)))%mod;
}
printf("%d
",1ll*(ans+n+mod)*inv2%mod);
return 0;
}