bzoj2154 (洛谷P1829) ： Crash的数字表格 （莫比乌斯反演）

bzoj2154 (洛谷P1829) ： Crash的数字表格

题意

求 $\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} lcm(i,j)$

题解

显然

$$lcm(i,j)=\frac{ij}{gcd(i,j)}$$

原式变成

$$\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} \frac{ij}{gcd(i,j)}$$

我们换成去枚举 $gcd$ ，可以得到

$$\sum\limits_{k=1}^{min(n,m)}\frac{1}{k}\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}ij[gcd(i,j)==k]$$
$$\sum\limits_{k=1}^{min(n,m)}\frac{1}{k}\sum\limits_{i=1}^{ \frac{n}{k}}\sum\limits_{j=1}^{\frac{m}{k}}ijk^2 \varepsilon(gcd(i,j)==1)$$
$$\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{i=1}^{ \frac{n}{k}}\sum\limits_{j=1}^{\frac{m}{k}}ij \varepsilon(gcd(i,j)==1)$$
$$\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{i=1}^{ \frac{n}{k}}\sum\limits_{j=1}^{\frac{m}{k}}\sum_{d|gcd(i,j)}\mu(d)ij$$

再去枚举gcd，得到

$$\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\frac{n}{k},\frac{m}{k})}\mu(d)\sum\limits_{i=1}^{\frac{n}{kd}}id\sum\limits_{j=1}^{\frac{m}{kd}}jd$$
$$\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\frac{n}{k},\frac{m}{k})}\mu(d)\cdot d^2\sum\limits_{i=1}^{\frac{n}{kd}}i\sum\limits_{j=1}^{\frac{m}{kd}}j$$
$$\sum\limits_{k=1}^{min(n,m)}k\sum\limits_{d=1}^{min(\frac{n}{k},\frac{m}{k})}\mu(d)\cdot d^2\cdot sum(1,\frac{n}{kd})\cdot sum(1,\frac{m}{kd})$$

---

代码

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
#define mod 20101009
#define N 10000010
#define ll long long
int prime[N],n,m,tot=0;
ll mu[N],ans=0;
bool notprime[N];
void mobius(){
    mu[1]=1;
    memset(notprime,false,sizeof(notprime));
    for(int i=2;i<=n;i++){
        if(!notprime[i]){
            mu[i]=-1;
            prime[++tot]=i;
        }for(int j=1;j<=tot,prime[j]*i<=n;j++){
            notprime[prime[j]*i]=true;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }mu[prime[j]*i]=-mu[i];
        }
    }for(int i=2;i<=n;i++) mu[i]=(mu[i-1]+mu[i]*i%mod*i%mod)%mod;
}
ll sum(int x,int y){return ((ll)(x+y)*(y-x+1)/2)%mod;}
ll f(int a,int b){
    if(a>b) swap(a,b);
    ll ans1=0;
    for(int i=1,last=1;i<=a;i=last+1){
        last=min(a/(a/i),b/(b/i));
        ans1=(ans1+(mu[last]-mu[i-1]+mod)%mod*sum(1,a/i)%mod*sum(1,b/i)%mod)%mod;
    }return ans1;
}
int main(){
    scanf("%d%d",&n,&m);
    if(n>m) swap(n,m);
    mobius();
    for(int i=1,last=1;i<=n;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ans=(ans+f(n/i,m/i)%mod*sum(i,last))%mod; 
    }printf("%lld",ans);
    return 0;
}