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;
}