【莫比乌斯反演】BZOJ2154 Crash的数字表格

Description

　　求sigma lcm(x,y)，x<=n,y<=m。n,m<=1e7。

 

Solution

　　lcm没有什么直接做的好方法，用lcm=x*y/gcd转成gcd来做

　　就是要求sigma d*f(x/d,y/d)

　　f(x,y)为x和y以内gcd正好为1的对数

　　F为所有对数，于是有F(x,y)=x*(x+1)/2*y*(y+1)/2

　　f(x,y)=sigma (1<=i<=x) i*i*mu(i)*F(x/i,y/i) 

　　f用莫比乌斯反演解决，这两个式子都套上分块优化到sqrt，于是总复杂度sqrt*sqrt=n

　　分块优化具体可以见上一篇blog

 

Code

　　一开始全开LL MLE了一发

　　然后又WA了两发，第一次是有一地方算的时候溢出

　　一开始为了避免MLE把prime数组/50，但其实只能/20的样子

#include<cstdio>
#include<algorithm>
#include<cstring>
#define ll long long
using namespace std;
const int maxn=1e7+5,mod=20101009;

bool flag[maxn];
int prime[maxn],mu[maxn],cnt;
int sum[maxn],s[maxn];
int n,m;

void getmu(){
    mu[1]=1;
    for(int i=2;i<=n;i++){
        if(!flag[i]){
            mu[i]=-1;
            prime[++cnt]=i;
        }
        for(int j=1;i*prime[j]<=n&&j<=cnt;j++){
            flag[i*prime[j]]=1;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++)
        sum[i]=(sum[i-1]+(ll)i*i*mu[i]%mod)%mod;
}

ll F(int x,int y){
    return (ll)((ll)x*(x+1)/2%mod)*((ll)y*(y+1)/2%mod)%mod;
}

ll f(int x,int y){
    ll ret=0;
    int p;
    for(int i=1;i<=x;i=p+1){
        p=min(x/(x/i),y/(y/i));
        ret=(ret+(ll)(sum[p]-sum[i-1])*F(x/i,y/i))%mod;
    }
    return ret;
}

int main(){
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    for(int i=1;i<=n;i++) s[i]=(s[i-1]+i)%mod;
    getmu();
    
    int pos;
    ll ans=0;
    for(int i=1;i<=n;i=pos+1){
        pos=min(n/(n/i),m/(m/i));
        ans=(ans+(ll)(s[pos]-s[i-1])*f(n/i,m/i))%mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}