Problem
Solution
\[\sum_{d=1}^n \sum_{i=1}^n\sum_{j=1}^m [\gcd(i,j)=d](\frac {ij} d)^d
\]
\[\sum_{d=1}^n d^d\sum_{i=1}^{n/d}\sum_{j=1}^{m/d} [\gcd(i,j)=1](ij)^d
\]
考虑后面的式子
\[F(k)=\sum_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=k](ij)^d
\]
\[G(k)=k^{2d}\sum_{i=1}^{x/k}\sum_{j=1}^{y/k}(ij)^d
\]
\[F(1)=\sum_{i=1}^x\mu(i)G(i)
\]
\[\sum_{d=1}^n d^d\sum_{k=1}^{n/d}\mu(k)k^{2d}\sum_{i=1}^{n/kd}i^d\sum_{j=1}^{m/kd}j^d
\]
考虑怎么暴力。枚举d再枚举k是调和级数，那么就要快速计算后面的式子，也可以通过调和级数预处理。
时间复杂度\(O(n\log n+n\ln n)\)
Code
#include <algorithm>
#include <cstdio>
using namespace std;
typedef long long ll;
const int maxn=500010,N=500000,mod=1e9+7;
template <typename Tp> inline int getmin(Tp &x,Tp y){return y<x?x=y,1:0;}
template <typename Tp> inline int getmax(Tp &x,Tp y){return y>x?x=y,1:0;}
template <typename Tp> inline void read(Tp &x)
{
x=0;int f=0;char ch=getchar();
while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
if(ch=='-') f=1,ch=getchar();
while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=getchar();
if(f) x=-x;
}
int n,m,lim,tot,res,ans,pri[maxn],vis[maxn],mu[maxn],md[maxn],sum[maxn];
int pls(int x,int y){return x+y>=mod?x+y-mod:x+y;}
int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
int sqr(int x){return (ll)x*x%mod;}
int power(int x,int y)
{
int res=1;
for(;y;y>>=1,x=(ll)x*x%mod)
if(y&1)
res=(ll)res*x%mod;
return res;
}
void init()
{
mu[1]=md[1]=sum[1]=1;
for(int i=2;i<=N;i++)
{
md[i]=1;
if(!vis[i]) pri[++tot]=i,mu[i]=mod-1;
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
vis[i*pri[j]]=1;
if(i%pri[j]==0){mu[i*pri[j]]=0;break;}
mu[i*pri[j]]=dec(0,mu[i]);
}
}
}
int main()
{
init();
read(n);read(m);
if(n>m) swap(n,m);
for(int i=1;i<=n;i++)
{
lim=n/i;res=0;
for(int j=2;i*j<=m;j++)
{
md[j]=(ll)md[j]*j%mod;
sum[j]=pls(sum[j-1],md[j]);
}
for(int j=1;j<=lim;j++)
res=pls(res,(ll)mu[j]*sqr(md[j])%mod*sum[n/i/j]%mod*sum[m/i/j]%mod);
res=(ll)res*power(i,i)%mod;
ans=pls(ans,res);
}
printf("%d\n",ans);
return 0;
}