bzoj2154 Crash的数字表格

#include<bits/stdc++.h>
#define MAXN 10000000
#define mod 20101009
#define inv 10050505
#define LL long long
using namespace std;
bool P[MAXN+1];
int prime[MAXN+1],mu[MAXN+1],size;
int sum_mu_dou[MAXN+1];
int sum_sig[MAXN+1];
LL n,m;
void init(LL n)
{
memset(P,true,sizeof(P));
P[1]=false;size=0;
mu[1]=1;
for(LL i=1;i<=n;i++)
{
if(P[i])
{
prime[++size]=i;
mu[i]=-1;
}
for(LL j=1;j<=size&&i*prime[j]<=n;j++)
{
P[i*prime[j]]=false;
if(!(i%prime[j]))
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(LL i=1;i<=n;i++)
{
sum_mu_dou[i]=(LL)(sum_mu_dou[i-1]+mu[i]*i%mod*i%mod)%mod;
sum_sig[i]=(LL)(sum_sig[i-1]+i)%mod;
}
}
LL S(LL n,LL m)
{
return (((n%mod*(n+1)%mod*m%mod)%mod*(m+1)%mod*inv%mod)%mod*inv%mod)%mod;
}
LL G(LL n,LL m)
{
if(n>m)swap(n,m);
LL ans=0;
LL i=1,j=0;
while(i<=n)
{
j=min(n/(n/i),m/(m/i));
ans=(LL)(ans+(sum_mu_dou[j]-sum_mu_dou[i-1])%mod*S(n/i,m/i)%mod)%mod;
i=j+1;
}
return (ans%mod+mod)%mod;
}
LL ANS(LL n,LL m)
{
if(n>m)swap(n,m);
LL ans=0;
LL i=1,j=0;
while(i<=n)
{
j=min(n/(n/i),m/(m/i));
ans=(LL)(ans+(sum_sig[j]-sum_sig[i-1])%mod*G(n/i,m/i)%mod)%mod;
i=j+1;
}
return (ans%mod+mod)%mod;

}
int main()
{
cin>>n>>m;
init(max(m,n));
cout<<ANS(n,m);
return 0;
}