BZOJ2154: Crash的数字表格

#include <bits/stdc++.h>
#define maxn 10000005
using namespace std;
typedef long long ll;
const ll mod=20101009;
ll mob[maxn];
ll sum[maxn];
int prime[maxn];
bool vis[maxn];
int cnt;
void getmob(int n)
{
    mob[1]=1;
    for(int i=2;i<n;++i)
    {
        if(!vis[i])
        {
            prime[cnt++]=i;
            mob[i]=-1;
        }
        for(int j=0;j<cnt&&prime[j]*i<maxn;++j)
        {
            vis[i*prime[j]]=true;
            if(i%prime[j]==0)
            {
                mob[i*prime[j]]=0;
                break;
            }
            else
            {
                mob[i*prime[j]]=-mob[i];
            }
        }
    }
}
void init()
{
    for(ll i=1;i<maxn;++i)
    {
        sum[i]=(sum[i-1]+(i*i)%mod*mob[i]+mod)%mod;
    }
}
ll solve(ll x,ll y)
{
    ll res=0;
    ll tmp;
    for(ll i=1;i<=min(x,y);i=tmp+1)
    {
        tmp=min(x/(x/i),y/(y/i));
        ll tmp_x=x/i,tmp_y=y/i;
        ll temp=((tmp_x+1)*tmp_x/2%mod)*((tmp_y+1)*tmp_y/2%mod)%mod;
        temp=(temp*(sum[tmp]-sum[i-1]+mod)%mod+mod)%mod;
        res=(res+temp)%mod;
    }
    return res;
}
int main()
{
    int n,m;
    getmob(maxn-2);
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    init();
    ll tmp;
    ll ans=0;
    for(ll i=1;i<=n;i=tmp+1)
    {
        tmp=min(n/(n/i),m/(m/i));
        ll temp=(tmp-i+1)*(tmp+i)/2%mod;
        ans=(ans+temp*solve(n/i,m/i)%mod+mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}