【BZOJ 2154】Crash的数字表格 （莫比乌斯+分块）

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
#include<queue>
#include<cmath>
using namespace std;
#define Mod 20101009
#define Maxn 10000010
#define LL long long

LL mu[Maxn],pri[Maxn],h[Maxn],pl;
bool q[Maxn];

LL mymin(LL x,LL y) {return x<y?x:y;}

void get_mu(LL mx)
{
pl=0;
memset(q,1,sizeof(q));
mu[1]=1;
for(LL i=2;i<=mx;i++)
{
if(q[i])
{
pri[++pl]=i;
mu[i]=-1;
}
for(LL j=1;j<=pl;j++)
{
if(i*pri[j]>mx) break;
q[i*pri[j]]=0;
if(i%pri[j]==0) mu[i*pri[j]]=0;
else mu[i*pri[j]]=-mu[i];
if(i%pri[j]==0) break;
}
}
h[1]=(mu[1]*1*1)%Mod;
for(LL i=2;i<=mx;i++) h[i]=(h[i-1]+mu[i]*i*i)%Mod;
}

LL get_g(LL x,LL y)
{
return ( ( ((x+1)*x/2)%Mod )*( ((y+1)*y/2)%Mod ) )%Mod;
}

LL get_f(LL n,LL m)
{
LL t;
if(n>m) t=n,n=m,m=t;

LL sq=(LL)ceil(sqrt((double)m));

LL ans=0;
for(LL i=1;i<=mymin(n,sq);i++)
{
ans=( ans+((mu[i]*i*i)%Mod)*get_g(n/i,m/i) )%Mod;
}
for(int i=sq+1;i<=n;)
{
int x=n/i,y=m/i;
int r1=n/x+1,r2=m/y+1;
if(r1>n+1) r1=n+1;
if(r2>n+1) r2=n+1;
int r=mymin(r1,r2);
ans=(ans+ ((h[r-1]+Mod-h[i-1])%Mod)*get_g(x,y) )%Mod;
i=r;
}
return ans;
}

LL get_ans(int n,int m)
{
LL ans=0;

LL sq=(LL)ceil(sqrt((double)m));
for(LL i=1;i<=mymin(sq,n);i++)
{
ans=(ans+i*get_f(n/i,m/i) )%Mod;
}

for(LL i=sq+1;i<=n;)
{
LL x=n/i,y=m/i;
LL r1=n/x+1,r2=m/y+1;
LL r=mymin(r1,r2);
if(r>m+1) r=m+1;
ans=( ans+(((r-i)*(i+r-1)/2)%Mod)*get_f(x,y) )%Mod;
i=r;
}
return ans;
}

int main()
{
int T;
T=1;

while(T--)
{
LL n,m,t;
scanf("%lld%lld",&n,&m);
if(n>m) t=n,n=m,m=t;
get_mu(m);

LL ans=get_ans(n,m);

printf("%lld\n",ans);
}
return 0;
}