Luogu 1829 - [国家集训队]Crash的数字表格 / JZPTAB （莫比乌斯反演、分块）

Luogu 1829 - [国家集训队]Crash的数字表格 / JZPTAB

题意

给定\(n,m\)，求解下式

\[\sum_{i=1}^n\sum_{j=1}^m lcm(i,j)\ (mod\ 20101009) \]

限制

\(1\leq n,m\leq 10^7\)

思路

令\(n\leq m\)

已知

\[lcm(i,j)=\frac{i\times j}{gcd(i,j)} \]

故原式可以化成

\[\sum_{i=1}^n\sum_{j=1}^m \frac{i\times j}{\gcd(i,j)}\\ =\sum_{i=1}^n\sum_{j=1}^m\sum_{d|i,d|j,gcd(\frac i d,\frac j d)=1}\frac{i\times j}d\\ =\sum_{i=1}^n\sum_{j=1}^m[\gcd(\frac i d,\frac j d)=1]\frac{i\times j}d \]

转为枚举\(d\)，可得

\[\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\frac{i\times j}d\times d^2[\gcd(i,j)=1]\\ =\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}i\times j[\gcd(i,j)=1] \]

将单位函数\([\gcd(i,j)=1]\)进行莫比乌斯反演，可得

\[\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{x|\gcd(i,j)}\mu(x)\times i\times j \]

转为枚举\(x\)，可得

\[\sum_{d=1}^nd\sum_{x=1}^{\lfloor\frac nd\rfloor}\mu(x)\sum_{x|i}^{\lfloor\frac nd\rfloor}\sum_{x|j}^{\lfloor\frac md\rfloor} i\times j \]

由于枚举的\(i,j\)为\(x\)的倍数，我们可以转为枚举倍数，得到

\[\sum_{d=1}^nd\sum_{x=1}^{\lfloor\frac nd\rfloor}\mu(x)\times x^2\sum_{i=1}^{\lfloor\frac n {dx}\rfloor}\sum_{j=1}^{\lfloor\frac m {dx}\rfloor} i\times j \]

发现后半部分可以直接根据等差数列来求解

\[\sum_{i=1}^{\lfloor\frac n {dx}\rfloor}\sum_{j=1}^{\lfloor\frac m {dx}\rfloor} i\times j\\ =\frac{1+\lfloor\frac n {dx}\rfloor}2\times \lfloor\frac n {dx}\rfloor\times\frac{1+\lfloor\frac m {dx}\rfloor}2\times \lfloor\frac m {dx}\rfloor \]

令其为一函数来表示，得

\[g(n,m)=\sum_{i=1}^n\sum_{j=1}^mi\times j\\ =\frac{(1+n)\times n}2\frac{(1+m)\times m}2 \]

再令

\[f(n,m)=\sum_{x=1}^n\mu(x)\times x^2\sum_{i=1}^{\lfloor\frac n x\rfloor}\sum_{j=1}^{\lfloor\frac m x\rfloor} i\times j\\ =\sum_{x=1}^n\mu(x)\times x^2\times g(\lfloor\frac nx\rfloor,\lfloor\frac mx\rfloor) \]

可以发现该函数可以运用分块求解，并且\(\sum_{x=1}^n\mu(x)\times x^2\)可以通过前缀和来预处理

将\(f(n,m)\)代回原式得到

\[\sum_{d=1}^nd\times f(\lfloor\frac n d\rfloor,\lfloor\frac m d\rfloor) \]

也可运用分块求解

时间复杂度为\(O(n+\sqrt n\sqrt n)\)

代码

Case Max (317ms/2000ms)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=20101009;
const int N=10000000;

int primcnt;
ll mu[N+50],prim[N+50];
bool vis[N+50];
ll sumpre[N+50];

void init(int n)
{
    memset(vis,false,sizeof vis);
    mu[1]=1;
    primcnt=0;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])
        {
            prim[primcnt++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<primcnt;j++)
        {
            int k=i*prim[j];
            if(k>n)
                break;
            vis[k]=true;
            if(i%prim[j]==0)
            {
                mu[k]=0;
                break;
            }
            else
                mu[k]=-mu[i];
        }
    }
    for(ll i=1;i<=n;i++)
        sumpre[i]=(sumpre[i-1]+mu[i]*i*i)%mod;
}

ll g(ll n,ll m)
{
    return ((n+1)*n/2)%mod*((m+1)*m/2%mod)%mod;
}

ll func(ll n,ll m)
{
    ll res=0;
    for(ll x=1;x<=n;)
    {
        ll nxt=min(n/(n/x),m/(m/x));
        res=(res+(sumpre[nxt]-sumpre[x-1])*g(n/x,m/x))%mod;
        x=nxt+1;
    }
    return res;
}

int main()
{
    ll n,m,ans=0;
    scanf("%lld%lld",&n,&m);
    if(n>m)
        swap(n,m);
    init(n);
    for(ll d=1;d<=n;)
    {
        ll nxt=min(n/(n/d),m/(m/d));
        ans=(ans+func(n/d,m/d)*((d+nxt)*(nxt-d+1)/2%mod)%mod)%mod;
        d=nxt+1;
    }
    printf("%lld\n",(ans+mod)%mod);
    
    return 0;
}