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

#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>#include<cmath>#define P 20101009using namespace std;typedef long long ll;inline char nc(){static char buf[100000],*p1=buf,*p2=buf;if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }return *p1++;}inline void read(ll &x){char c=nc(),b=1;for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;}inline int write(ll x,char c){if (x==0) return putchar('0'),putchar(c),0;if (x<0) putchar('-'),x=-x;char s[50]={0},len=0;while (x) s[++len]=x%10+'0',x/=10;for (int i=len;i;i--) putchar(s[i]); putchar(c);}const int maxn=10000000;int mobius[maxn+5];int vst[maxn+5],prime[maxn+5],num=0;ll f[maxn+5];inline void Pre(){mobius[1]=1;for (int i=2;i<=maxn;i++){if (!vst[i])mobius[i]=-1,prime[++num]=i;for (int j=1;j<=num && prime[j]*i<=maxn;j++){vst[i*prime[j]]=1;if (i%prime[j]==0){mobius[i*prime[j]]=0; break;}elsemobius[i*prime[j]]=-mobius[i];}}for (ll i=1;i<=maxn;i++)f[i]=i*i*mobius[i];for (ll i=1;i<=maxn;i++)(f[i]+=f[i-1])%=P;}ll n,m;ll ans;inline ll sum(ll x,ll y){return ((x*(x+1)/2)%P)*((y*(y+1)/2)%P)%P;}inline ll F(ll x,ll y){if (x>y) swap(x,y);ll ret=0,i,j;for (j,i=1;i<=x;i=j+1){j=min(x/(x/i),y/(y/i));(ret+=(f[j]-f[i-1])*sum(x/i,y/i)%P)%=P;}return ret;}int main(){read(n); read(m);Pre();if (n>m) swap(n,m);for (ll j,i=1;i<=n;i=j+1){j=min(n/(n/i),m/(m/i));(ans+=(j+i)*(j-i+1)/2%P*F(n/i,m/i)%P)%=P;}write((ans+P)%P,'\n');return 0;}