【BZOJ】2154: Crash的数字表格 莫比乌斯反演

【题意】给定n,m，求Σlcm(i,j),1<=i<=n,1<=j<=m，n,m<=10^7。

【算法】数论（莫比乌斯反演）

【题解】

$$ans=sum_{ileq n}sum_{jleq m}frac{i*j}{gcd(i,j)}$$

$$ans=sum_{dleq min(n,m)}1/dsum_{ileq n}sum_{jleq m}[gcd(i,j)=d]i*j$$

$$ans=sum_{dleq min(n,m)}dsum_{ileq n/d}sum_{jleq m/d}[gcd(i,j)=1]i*j$$

发现后面部分只和n/d，m/d有关，于是封装后分块取值优化。

★$$ans=sum_{dleq min(n,m)}d*F(n/d,m/d)$$

$$F(n,m)=sum_{ileq n}sum_{jleq m}[gcd(i,j)=1]i*j$$

运用e=i*μ反演易得

★$$F(n,m)=sum_{dleq min(n,m)}mu (d)*d^2*sum(n/d,m/d)$$

$$sum(n,m)=sum_{ileq n}sum_{jleq m}i*j=sum_{ileq n}isum_{jleq m}j$$

★$$sum(n,m)=frac{n(n+1)}{2}*frac{m(m+1)}{2}$$

（这步由一般分配律）

最后就是两次分块取值优化，√n*√n，复杂度O(n)。

#include<cstdio>
#include<algorithm>
using namespace std;
const int maxn=10000010,MOD=20101009;
int sum[maxn],n,m,miu[maxn],prime[maxn],tot;
bool mark[maxn];
int M(int x){return x>=MOD?x-MOD:x;}
int SUM(int x,int y){return 1ll*x*(x+1)/2%MOD*(1ll*y*(y+1)/2%MOD)%MOD;}
int solve(int x,int y){
    int z=min(x,y),ans=0,pos=1;
    for(int i=1;i<=z;i=pos+1){
        pos=min(x/(x/i),y/(y/i));
        ans=(ans+1ll*(sum[pos]-sum[i-1])*SUM(x/i,y/i)%MOD)%MOD;
    }
    return ans;
}
int main(){
    scanf("%d%d",&n,&m);
    int z=min(n,m);
    miu[1]=sum[1]=1;
    for(int i=2;i<=z;i++){
        if(!mark[i])miu[prime[++tot]=i]=-1;
        for(int j=1;j<=tot&&i*prime[j]<=z;j++){
            mark[i*prime[j]]=1;
            if(i%prime[j]==0)break;
            miu[i*prime[j]]=-miu[i];
        }
        sum[i]=(sum[i-1]+1ll*i*i*miu[i]%MOD)%MOD;
    }
    int pos=1,ans=0;
    for(int i=1;i<=z;i=pos+1){
        pos=min(n/(n/i),m/(m/i));
        ans=(ans+1ll*(pos+i)*(pos-i+1)/2%MOD*solve(n/i,m/i))%MOD;
    }
    printf("%d",(ans+MOD)%MOD);
    return 0;
}