BZOJ-2154: Crash的数字表格

Description

N, M ≤ 10^7。求的值。

推式子！！！

=
=
=
=
设
=
对于不同的d，只有√n个值，对于不同的n/d，也只有√n个值。所以分块两次，总时间复杂度O(n)

#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 10001000 
using namespace std;
typedef long long ll;
const ll mod=20101009;
int pri[701000],cnt;
bool v[N];
ll n,m,ans,mu[N];
void quick_mu()
{
    mu[1]=1;
    for(int i=2;i<=n;++i)
    {
        if(!v[i]) pri[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*pri[j]<=n;++j)
        {
            v[i*pri[j]]=1;
            if(i%pri[j]==0){mu[i*pri[j]]=0; break;}
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;++i) mu[i]=(mu[i]*i%mod*i%mod+mu[i-1])%mod;
}
inline ll getsum(ll x,ll y)
{
    x=(x*(x+1)/2)%mod%mod;
    y=(y*(y+1)/2)%mod%mod;
    return x*y%mod;
}
ll work(ll x)
{
    ll nn=n/x,mm=m/x,b=1,last,re=0;
    while(b<=nn)
    {
        last=min(nn,min(nn/(nn/b),mm/(mm/b)));
        re=(re+(mu[last]-mu[b-1]+mod)%mod*getsum(nn/b,mm/b)%mod)%mod;
        b=last+1;
    }
    return re;
}
int main()
{
    scanf("%lld%lld",&n,&m); if(n>m) swap(n,m);
    quick_mu();
    ll b=1,last;
    while(b<=n)
    {
        last=min(n,min(n/(n/b),m/(m/b)));
        ans=(ans+(last+b)*(last-b+1)/2%mod*work(b)%mod)%mod;
        b=last+1;
    }
    printf("%lld",ans);
    return 0;
}