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

http://cjjlsdy.blog.163.com/blog/static/18037056320146198140985/

#include<cstdio>
#include<cstring>
#include<iostream>
#define maxn 10000021
#define LL long long
#define Mod 20101009ll
using namespace std;
bool vis[maxn];
LL up,n,m,u[maxn],sum[maxn],pri[maxn],cnt;
void pre(){
	u[1]=1;sum[1]=1;
	for(LL i=2;i<=up;i++){
		if(!vis[i])u[i]=-1,pri[++cnt]=i;
		for(int j=1;j<=cnt&&pri[j]*i<=up;j++){
			vis[pri[j]*i]=1;
			u[pri[j]*i]=-u[i];
			if(i%pri[j]==0){u[pri[j]*i]=0;break;}
		}
	}
	for(LL i=1;i<=up;i++)sum[i]=(sum[i-1]+u[i]*i%Mod*i%Mod)%Mod;
}
LL Q(LL x,LL y){return (x*(x+1)/2ll%Mod)*(y*(y+1)/2ll%Mod)%Mod;}
LL f(LL x,LL y){
	LL p=min(x,y),ans=0,r;
	for(LL l,r,d=1;d<=p;d=r+1){
		l=d,r=min(x/(x/d),y/(y/d));
		ans=(ans+(sum[r]-sum[l-1])*Q(x/d,y/d)%Mod+Mod)%Mod;
	}return ans;
}
int main(){
	scanf("%lld%lld",&n,&m);up=min(n,m);
	pre();LL r,ans=0;
	for(LL l=1;l<=up;l=r+1){
		r=min(n/(n/l),m/(m/l));
		ans=(ans+(r+l)*(r-l+1)/2ll%Mod*f(n/l,m/l)%Mod)%Mod;
	}printf("%lld",ans);
	return 0;
}