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

Code:

#include<bits/stdc++.h>
#define maxn 10010006
#define ll long long  
using namespace std;
const long long mod = 1ll*20101009;  
void setIO(string s)
{
	string in=s+".in"; 
	freopen(in.c_str(),"r",stdin); 
}
int cnt; 
int prime[maxn], vis[maxn], mu[maxn]; 
ll sumv[maxn]; 
inline void init(int hh) 
{
	int i,j;
	mu[1]=sumv[1]=1; 
	for(i=2;i<hh;++i)
	{
		if(!vis[i]) prime[++cnt]=i, mu[i]=-1; 
		for(j=1;j<=cnt&&1ll*prime[j]*i<hh;++j) 
		{
			vis[prime[j]*i]=1; 
			if(i%prime[j]!=0) mu[i*prime[j]]=-mu[i]; 
			else 
			{
				mu[i*prime[j]]=0; 
				break; 
			}
		}
		sumv[i]=(sumv[i-1]+1ll*i*i%mod*(mu[i]+mod))%mod;    
		//  sumv2[i]=(i+sumv2[i-1])%mod; 
	} 
}
inline ll getarr(ll n,ll m)
{
	ll tmp=(1ll*n*(n+1)/2%mod) * (1ll*m*(m+1)/2%mod)%mod;  
    return tmp; 
}
inline ll calc(ll n,ll m)
{
	ll i,j; 
	ll ans=0; 
	for(i=1;i<=n;i=j+1)
	{
		j=min(n/(n/i), m/(m/i));  
		ans=(ans+1ll*(sumv[j]-sumv[i-1]+mod)*getarr(n/i,m/i)%mod)%mod;  
	}
	return ans; 
}
int main()
{
// 	setIO("input");  
	ll n,m,i,j; 
	ll ans=0; 
	scanf("%lld%lld",&n,&m); 
	if(n>m)swap(n,m);        
	init(m + 1); 
	for(i=1;i<=n;i=j+1)
	{
		j=min(n/(n/i), m/(m/i));  
		ans=(ans+1ll*((j-i+1)*(i+j)/2)%mod*calc(n/i,m/i)%mod)%mod;   
		ans=(ans+mod)%mod;  
	}  
	printf("%lld\n",ans); 
	return 0; 
}