Problem Eight:[Crash的数字表格/BZOJ2693 jzptab（多组）]

http://www.lydsy.com/JudgeOnline/problem.php?id=2154

$$\begin{equation*} ans=\sum_{i=1}^n \sum_{j=1}^m lcm(i,j)=\sum_{i=1}^n \sum_{j=1}^m i*j/gcd(i,j)\\ f(x,y)=\sum_{i=1}^x \sum_{j=1}^y i*j*[gcd(i,j)==1]\\ ans=\sum_{d=1}^{min(n,m)}d*f(n/d,m/d)\\ f(x,y)=\sum_{i=1}^x \sum_{j=1}^y i*j*\sum_{k|i,k|j}\mu(k)\\ =\sum_{k=1}^{min(x,y)}k^2*\mu(k)*\sum_{i=1}^{x/k}\sum_{j=1}^{y/k}i*j\\ 另F(x,y)=\sum_{i=1}^{x}\sum_{j=1}^{y}i*j\\ f(x,y)=\sum_{k=1}^{min(x,y)}k^2*\mu(k)*F(x/k,y/k)\\ ans=\sum_{d=1}^{min(n,m)}d*\sum_{k=1}^{min(n/d,m/d)}k^2*\mu(k)*F(n/d/k,m/d/k)\\ 另D=d*k\\ ans=\sum_{D=1}^{min(n,m)}F(n/D,m/D)*\sum_{i|D}i^2*\mu(i)*(D/i)\\ 另g(D)=\sum_{i|D}i^2*\mu(i)*(D/i)\\ ans=\sum_{D=1}^{min(n,m)}F(n/D,m/D)*g(D) \end{equation*}$$

#include<cstdio>
#include<iostream>
#include<cstring>
using namespace std;
const int mod=20101009;
typedef long long LL;
const int N=1e7+100;
int p[N],mu[N],pri[N],pn;
LL g[N];

void init(){
    int n=1e7;
    //printf("%d\n",n);
    memset(p,0,sizeof(p));
    mu[0]=0;pn=0;mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!p[i]){
            mu[i]=-1;
            pri[pn++]=i;
            g[i]=(LL)(1-i)%mod;
        }
        for(int j=0;j<pn;j++){
            if((LL)i*pri[j]>n)break;
            p[i*pri[j]]=1;
            if(i%pri[j]==0){
                mu[i*pri[j]]=0;
                g[i*pri[j]]=g[i]%mod;
                break;
            }
            else {
                mu[i*pri[j]]=-mu[i];
                g[i*pri[j]]=(g[i]*g[pri[j]])%mod;
            }
        }
    }
    g[0]=0;g[1]=1;
    for(int i=1;i<=n;i++){
        g[i]*=i;g[i]%=mod;
        g[i]+=g[i-1];g[i]%=mod;
    }
}

LL F(int x,int y){
    LL xx=(LL)(1+x)*x/2;xx%=mod;
    LL yy=(LL)(1+y)*y/2;yy%=mod;
    return (xx*yy)%mod;
}

int main(){
    #ifdef DouBi
    freopen("in.cpp","r",stdin);
    //freopen("out1.cpp","w",stdout);
    #endif // DouBi
    init();
    int n,m;
    while(scanf("%d%d",&n,&m)!=EOF){
        int r=min(n,m);
        LL ans=0;
        while(r>0){
            //printf("%d\n",r);
            int k1=n/r;
            int k2=m/r;
            int l1=n/(k1+1)+1;
            int l2=m/(k2+1)+1;
            int l=max(l1,l2);
            ans+=(F(k1,k2)*(g[r]-g[l-1]))%mod;ans%=mod;
            r=l-1;
        }
        printf("%lld\n",(ans+mod)%mod);
    }
    return 0;
}