BZOJ3561: DZY Loves Math VI

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cstdlib>
using namespace std;
#define ll long long
const
ll mod=1000000007;
const
ll maxn=500001;
ll mu[maxn],mul[maxn],pre[maxn];
ll prime[maxn],tot;
bool check[maxn];
char c;
inline void read(ll &a)
{
a=0;do c=getchar();while(c<'0'||c>'9');
while(c<='9'&&c>='0')a=(a<<3)+(a<<1)+c-'0',c=getchar();
}
inline ll D(ll x)
{
ll y=x,base=1,res=1;
while(x)
{
if(x&base)x^=base,res=res*y;
base<<=1;
y*=y;
if(y>=mod)y%=mod;
if(res>=mod)res%=mod;
}
return res;
}
int main()
{
ll d,t,i,j,k;
mu[1]=1;
for(i=2;i<maxn;i++)
{
if(!check[i])prime[++tot]=i,mu[i]=-1;
for(j=1;j<=tot;j++)
{
k=i*prime[j];
if(k>=maxn)break;
if(i%prime[j]==0){check[k]=true,mu[k]=0;break; }
check[k]=true,mu[k]=-mu[i];
}
}
for(i=1;i<maxn;i++)mul[i]=1;
ll n,m;
read(n),read(m);
if(m<n)swap(n,m);
ll ans=0,base,ed,ed2,base2;
for(d=1;d<=n;d++)
{
base=D(d),ed=n/d,ed2=m/d;
for(i=1;i<=ed2;i++)
{
mul[i]*=i;
if(mul[i]>=mod)
mul[i]%=mod;
pre[i]=mul[i]+pre[i-1];
if(pre[i]>=mod)pre[i]-=mod;
}
t=1;
for(t=1;t<=ed;t++)
if(mu[t])             ans=(ans+((base*mul[t]%mod)*mul[t]*mu[t]%mod)*(pre[n/(t*d)]*pre[m/(t*d)]%mod))%mod;
}
printf("%lld\n",ans);
return 0;
}