【莫比乌斯反演】[国家集训队]Crash的数字表格 / JZPTAB

题目链接

P1829 [国家集训队]Crash的数字表格 / JZPTAB

题目大意

求：

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

题解

首先将 $\mathrm{lcm}$ 转成 $\gcd$：

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

考虑对于每一个 $\gcd (i,j)$，如何求出他的所有贡献：

$$\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{d}\times [\gcd (i,j)=d]$$

由于本质上一定有 $i=i^{\prime }d,j=j^{\prime }d$，所以我们换一种枚举方式，把 $d$ 全部提出来，则有 $\frac{ij}{d}=d{i}^{\prime }{j}^{\prime }$：

$$\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }ij\times [\gcd (i,j)=1]$$

根据反演结论：$[n=1]=\sum _{d\mid n}\mu (d)$ 可得：

$$\sum _{d=1}^{n}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }ij\times \sum _{d^{\prime }\mid \gcd (i,j)}\mu (d^{\prime })$$

令 $s(n,m)=\sum _{i=1}^n\sum _{j=1}^{m}ij\sum _{d\mid \gcd (i,j)}\mu (d)$，则原式化为：

$$\sum _{d=1}^{n}ds\left(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor \right)\text{\hspace{1em}(1)}$$

现在考虑 $s(n,m)$ 怎么求，我们把 $\gcd (i,j)$ 的所有约数都枚举一下，改变一下式子：

$$\sum _{d=1}^n\sum _{i=1}^n\sum _{j=1}^{m}ij\times [d\mid i\text{and}d\mid j]\mu (d)$$

将所有式子归位：

$$\sum _{d=1}^n\mu (d)\sum _{i=1}^n(i\times [d\mid i])\sum _{j=1}^m(j\times [d\mid j])$$

发现此时还是可以令 $i=i^{\prime }d,j=j^{\prime }d$，改变枚举关系：

$$\sum _{d=1}^n\mu (d)d^2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }(i\times [1\mid i])\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }(j\times [1\mid j])$$

显然 $[1\mid d]=1$，所以这个式子等价于：

$$\sum _{d=1}^n\mu (d)d^2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }j\text{\hspace{1em}(2)}$$

令 $\mathrm{sum}(n,m)=\sum _{i=1}^n\sum _{j=1}^{m}ij$，显然这个式子就等于：

$$\frac{n(n+1)}{2}\times \frac{m(m+1)}{2}$$

可以 $\mathcal{O}(1)$ 处理。考虑式 $(2)$，我们可以提前预处理出 $\mu (d)d^2$ 的前缀和，观察 $\lfloor \frac{n}{d}\rfloor$，显然是可以对他数论分块，考虑一大块 $\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor$ 均相同的 $d$，可以直接利用前缀和 $\mathcal{O}(1)$ 求出该块的和，所以这里运用数论分块可以达到 $\mathcal{O}(\sqrt{n})$ 的复杂度。每一块的答案是：

$$(\mu (r)r^2-\mu (l)l^2)\times \mathrm{sum}\left(\lfloor \frac{n}{l}\rfloor ,\lfloor \frac{m}{l}\rfloor \right)$$

回到式 $(1)$，显然对于一大块 $\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor$ 均相同的 $d$，$s_1$ 的值仍然不变，所以我们再利用一次数论分块，这一层也是 $\mathcal{O}(\sqrt{n})$ 的复杂度。每一块的答案是：

$$\frac{l+r}{2}\times (r-l+1)\times s\left(\lfloor \frac{n}{l}\rfloor ,\lfloor \frac{m}{l}\rfloor \right)$$

总复杂度的瓶颈在于线性筛，总时间复杂度 $\mathcal{O}(n)$。

代码

//你 cnt=1 了吗？
#include<bits/stdc++.h>
#define int long long
#define mid ((l+r)>>1)
#define fir first
#define sec second
#define lowbit(i) (i&(-i))
using namespace std;
const int N=2e7+5;
const int inf=1e18;
struct edge{int to,nxt,l;};
inline int read(){
    char op=getchar();
    int w=0,s=1;
    while(op<'0'||op>'9'){
        if(op=='-') s=-1;
        op=getchar();
    }
    while(op>='0'&&op<='9'){
        w=(w<<1)+(w<<3)+op-'0';
        op=getchar();
    }
    return w*s;
}
//数论部分
const double pi=acos(-1);
const int mod=20101009;
int Mul(int a,int b){return (a%mod*b%mod)%mod;}
int Add(int a,int b){return (a+b)%mod;}
int Dec(int a,int b){return (a-b+mod)%mod;}
int Pow(int a,int k){
    int ans=1;
    while(k){
        if(k&1) ans=Mul(ans,a);
        a=Mul(a,a);
        k>>=1;
    }
    return ans;
}
int gcd(int x,int y){return y==0?x:gcd(y,x%y);}
int lcm(int x,int y){return x/gcd(x,y)*y;}
int inv(int x){return Pow(x,mod-2);}
void exgcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
}
const int inv2=inv(2);
int n,m,mu[N],vis[N],p[N],sum[N],cnt=0;
void mobius(){
    mu[1]=1;
    for(register int i=2;i<=min(n,m);i++){
        if(!vis[i]){
            mu[i]=-1;
            p[++cnt]=i;
        }
        for(register int j=1;j<=cnt&&p[j]*i<=min(n,m);j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=min(n,m);i++) sum[i]=Add(sum[i-1],Mul(Mul(i,i),mu[i]));
}
int solve(int x,int y){
    return Mul(Mul(Mul(x,x+1),inv2),Mul(Mul(y,y+1),inv2));
}
int sqrt2(int n,int m){
    /* 第二层式子
     * ans=\sum\limits_{i=1}^n mu(d)\times d^2 \sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{m/d}ij
     */
    int l=1,r,ans=0;
    for(l=1;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=Add(ans,Mul(solve(n/l,m/l),sum[r]-sum[l-1]));
    }
    return ans;
}
int sqrt1(int n,int m){
    /* 第一层式子
     * ans=\sum\limits_{d=1}^n d\times f(n/d,m/d)
     */
    int l=1,r,ans=0;
    for(l=1;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans=Add(ans,Mul(Mul(Mul(l+r,inv2),r-l+1),sqrt2(n/l,m/l)));
    }
    return ans;
}
signed main(){
    n=read(),m=read();
    mobius();
    printf("%lld\n",(sqrt1(n,m)+mod)%mod);
}