习题：Crash的数字表格（莫比乌斯反演）

题目

传送门

思路

貌似可以直接推式子来做

\[\begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)&=\sum_{i=1}^n\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}\\&=\sum_{t=1}^{min(n,m)}\frac{1}{t}\sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)==t]\\&=\sum_{t=1}^{min(n,m)}\frac{1}{t}\sum_{i=1}^{it\le n}\sum_{j=1}^{jt\le m}itjt[gcd(i,j)==1]\\&=\sum_{t=1}^{min(n,m)}\frac{1}{t}\sum_{i=1}^{it\le n}\sum_{j=1}^{jt\le m}itjt\sum_{d|gcd(i,j)}\mu (d)\\&=\sum_{t=1}^{min(n,m)}t\sum_{i=1}^{it\le n}\sum_{j=1}^{jt\le m}ij\sum_{d|gcd(i,j)}\mu (d)\\&=\sum_{t=1}^{min(n,m)}t\sum_{d=1}^{\frac{min(n,m)}{t}}\mu(d)d^2\sum_{i=1}^{itd\le n}\sum_{j=1}^{jtd\le m}ij\\&=\sum_{t=1}^{min(n,m)}t\sum_{d=1}^{\frac{min(n,m)}{t}}\mu(d)d^2(\sum_{i=1}^{itd\le n}i)(\sum_{j=1}^{jtd\le m}j)\\&=\sum_{t=1}^{min(n,m)}t\sum_{d=1}^{\frac{min(n,m)}{t}}\mu(d)d^2(\sum_{i=1}^{\frac{n}{td}}i)(\sum_{j=1}^{\frac{m}{td}}j)\end{aligned} \]

注意到这里虽然是\(\frac{n}{td}\)，但是实际上应该是\(\lfloor\frac{\lfloor\frac{n}{t}\rfloor}{d}\rfloor\)

即可以进行两次数论分块，具体而言，

\[设g(n,m)=\sum_{d=1}^{min(n,m)}\mu(d)d^2\sum_{i=1}^{\frac{n}{d}}i\sum_{i=1}^{\frac{m}{d}}i\\那么原式即为\sum_{t=1}^{min(n,m)}t*g(\frac{n}{t},\frac{m}{t}) \]

代码

#include<iostream>
#include<cmath>
using namespace std;
const int mod=20101009;
const int inv2=10050505;
int n,m;
bool vis[10000005];
int pri[1000005],lenp;
long long mu[10000005];
long long f[10000005];
long long ans;
void sieve(int n)
{
    mu[1]=1;
    f[1]=1;
    for(int i=2;i<=n;i++)
    {
        f[i]=(f[i-1]+i)%mod;
        if(vis[i]==0)
        {
            pri[++lenp]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=lenp&&1ll*i*pri[j]<=n;j++)
        {
            vis[pri[j]*i]=1;
            mu[pri[j]*i]=-mu[i];
            if(i%pri[j]==0)
            {
                mu[i*pri[j]]=0;
                break;
            }
        }
        mu[i]=mu[i]*i*i%mod;
		mu[i]=(mu[i]+mu[i-1])%mod;
    }
}
long long getsum(long long x)
{
    return x*(x+1)%mod*inv2%mod;
}
long long solve(int n,int m)
{
    long long temp=0;
    for(int l=1,r;l<=min(n,m);l=r+1)
    {
        r=min(min(n,m),min(n/(n/l),m/(m/l)));
        temp=(temp+(mu[r]-mu[l-1])%mod*getsum(n/r)%mod*getsum(m/r)%mod)%mod;
    }
    temp=(temp%mod+mod)%mod;
    return temp;
}
int main()
{
    cin>>n>>m;
    sieve(min(n,m));
    for(int l=1,r;l<=min(n,m);l=r+1)
    {
        r=min(min(n,m),min(n/(n/l),m/(m/l)));
        //cout<<"debug:"<<solve(n/r,m/r)<<'\n';
        ans=(ans+((f[r]-f[l-1])%mod+mod)%mod*solve(n/r,m/r)%mod)%mod;
    }
    cout<<ans;
    return 0;
}
/*
n/td
n/i
*/