BZOJ 2154 Crash的数字表格 莫比乌斯反演

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<ctime>
#include<cmath>
#include<algorithm>
#include<queue>
#include<set>
#include<map>
#include<iomanip>
using namespace std;
#define LL long long
#define up(i,j,n) for(int i=j;i<=n;i++)
#define pii pair<LL,LL>
#define db double
#define eps 1e-4
#define FILE "dealing"
int read(){
int x=0,f=1,ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0',ch=getchar();}
return x*f;
}
const LL maxn=(LL)(1e7+2),inf=1000000000000000LL,limit=(LL)(1e7+2),mod=20101009;
bool cmin(LL& a,LL b){return a>b?a=b,true:false;}
bool cmax(LL& a,LL b){return a<b?a=b,true:false;}
LL n,m;
LL mu[maxn],b[maxn],prime[1000000],tail=0;
LL sum(LL x,LL y){
LL c=x%2?(x+1)/2*x%mod:x/2*(x+1)%mod;
LL d=y%2?(y+1)/2*y%mod:y/2*(y+1)%mod;
return c*d%mod;
}
void getmu(){
mu[1]=1;
up(i,2,n){
if(!b[i])prime[++tail]=i,mu[i]=-1;
for(LL j=1;prime[j]*i<=n;j++){
b[i*prime[j]]=1;
if(i%prime[j]){mu[i*prime[j]]=-mu[i];}
else {mu[i*prime[j]]=0;break;}
}
}
up(i,1,n)mu[i]=((LL)i*i%mod*mu[i])%mod;
up(i,1,n)mu[i]=(mu[i]+mu[i-1])%mod;
}
LL f[1000][1000];
LL getf(LL x,LL y){
if(x<1000&&y<1000&&f[x][y])return f[x][y];
LL ans=0,pos;
for(int i=1;i<=x;i++){
pos=min(x/(x/i),y/(y/i));
ans=(ans+(mu[pos]-mu[i-1]+mod)%mod*sum(x/i,y/i)%mod)%mod;
i=pos;
}
if(x<1000&&y<1000)f[x][y]=ans;
return ans;
}
int main(){
freopen(FILE".in","r",stdin);
freopen(FILE".out","w",stdout);
n=read(),m=read();
if(n>m)swap(n,m);
getmu();
LL ans=0;
up(i,1,n)ans=(ans+i*getf(n/i,m/i))%mod;
printf("%lld\n",ans);
//cout<<clock()<<endl;
return 0;
}