BZOJ 3561 DZY Loves Math VI

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define maxn 500500
#define mod 1000000007
using namespace std;
long long n,m,a[maxn],sum[maxn],ans=0,prime[maxn],tot=0,miu[maxn];
bool vis[maxn];
void get_table()
{
miu[1]=1;
for (long long i=2;i<=maxn-500;i++)
{
if (!vis[i]) {prime[++tot]=i;miu[i]=-1;}
for (long long j=1;i*prime[j]<=maxn-500 && j<=tot;j++)
{
vis[i*prime[j]]=true;
if (i%prime[j]) miu[i*prime[j]]=-miu[i];
else {miu[i*prime[j]]=0;break;}
}
}
}
long long f_pow(long long x,long long y)
{
long long ans=1,base=x;
while (y)
{
if (y&1) ans=(ans*base)%mod;
base=(base*base)%mod;
y>>=1;
}
return ans;
}
int main()
{
scanf("%lld%lld",&n,&m);get_table();
if (n>m) swap(n,m);
for (long long i=1;i<=m;i++) a[i]=1,sum[i]=i;
for (long long i=1;i<=n;i++)
{
for (long long j=1;j*i<=m;j++) {a[j]=a[j]*j%mod;sum[j]=(sum[j-1]+a[j])%mod;}
long long ret=0;
for (long long j=1;j*i<=n;j++) {ret+=miu[j]*sum[n/i/j]*sum[m/i/j]%mod*a[j]%mod*a[j]%mod;ret%=mod;}
ans+=ret*f_pow(i,i);ans%=mod;
}
printf("%lld\n",ans);
return 0;
}