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

https://www.luogu.com.cn/problem/P1829

莫比乌斯反演

\[令n\le m \\ \sum_{i=1}^n\sum_{j=1}^m lcm(i,j)\\ =\sum_{k=1}^n \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{k} \rfloor} ijk[\gcd(i,j)=1]\\ =\sum_{k=1}^n k \sum_{i=1}^{\lfloor \frac{n}{k} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{k} \rfloor} ij \sum_{d|i,d|j}\mu(d)\\ =\sum_{k=1}^n k \sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d)d^2 \sum_{i=1}^{\lfloor \frac{n}{kd} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{kd} \rfloor}ij\\ 令g(\lfloor \frac{n}{k} \rfloor,\lfloor \frac{m}{k} \rfloor)=\sum_{d=1}^{\lfloor \frac{n}{k} \rfloor} \mu(d)d^2 \sum_{i=1}^{\lfloor \frac{n}{kd} \rfloor}\sum_{j=1}^{\lfloor \frac{m}{kd} \rfloor}ij，直接数论分块处理\\ 原式化为\sum_{k=1}^n k\times g(\lfloor \frac{n}{k} \rfloor,\lfloor \frac{m}{k} \rfloor)\\ 再次数论分块即可 \]

注意点：数据范围较小，能够支持两次数论分块

\(C++ Code:\)

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cmath>
#define N 10000005
#define mod 20101009
#define ll long long
using namespace std;
ll n,m;
bool pri[N];
int cnt,prime[N],mu[N];
ll smu[N];
ll ans=0;
void Pre()
{
    mu[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!pri[i])
        {
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=cnt;j++)
        {
            ll g=(ll)i*prime[j];
            if (g>n)
                break;
            pri[g]=true;
            if (i%prime[j]==0)
            {
                mu[g]=0;
                break;
            }
            mu[g]=mu[i]*mu[prime[j]];
        }
    }
    for (int i=1;i<=n;i++)
        smu[i]=(smu[i-1]+(ll)mu[i]*i%mod*i%mod)%mod;
}
ll g(ll n,ll m)
{
    if (n>m)
        swap(n,m);
    ll ans=0;
    ll l,r;
    for (l=1;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));
        ll w1=n/l,w2=m/l;
        ll e=(w1+1)*w1/2;
        e%=mod;
        ll o=(w2+1)*w2/2;
        o%=mod;
        ll d=(smu[r]-smu[l-1])%mod;
        ans=(ans+e*o%mod*d%mod)%mod;
    }
    ans=(ans%mod+mod)%mod;
    return ans;
}
int main()
{
    scanf("%lld%lld",&n,&m);
    if (n>m)
        swap(n,m);
    Pre();
    ll l,r;
    for (l=1;l<=n;l=r+1)
    {
        r=min(n/(n/l),m/(m/l));//注意：n,m都要分块，虽然这一步式子中没有，但会对g函数产生影响
        ll e=(l+r)*(r-l+1)/2;
        e%=mod;
        ll o=g(n/l,m/l);
        ans=(ans+e*o%mod)%mod;
    }
    ans=(ans%mod+mod)%mod;
    cout << ans << endl;
    return 0;
}