Solution:SP8099(TABLE - Crash's number table)

· · 题解

Link:SP8099

\text P1829 为此题的翻译版，数据范围、题面、题目背景、标题都一模一样。

以及此题只能交 C++98。

求

\sum_{i=1}^n\sum_{j=1}^m{\rm lcm}(i,j)

套路地转换 {\rm lcm}：

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

枚举 i,j 的 \gcd：

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

枚举 i,j 为 d 的倍数：

\sum_{d=1}^{\min(n,m)}\sum_{i=1}^{\left\lfloor\frac nd\right\rfloor}\sum_{j=1}^{\left\lfloor\frac md\right\rfloor}dij[\gcd(i,j)=1]

下取整想必都知道是怎么来的。这个时候套 \rm Möbius 反演公式：

\sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\left\lfloor\frac nd\right\rfloor}i\sum_{j=1}^{\left\lfloor\frac md\right\rfloor}j\sum_{k|i,k|j}\mu(k)

交换求和顺序，先枚举 k：

\sum_{k=1}^{\min(n,m)}\mu(k)\sum_{d=1}^{\min(n,m)}d\sum_{i=1}^{\left\lfloor\frac n{dk}\right\rfloor}ik\sum_{j=1}^{\left\lfloor\frac m{dk}\right\rfloor}jk=\sum_{k=1}^{\min(n,m)}\mu(k)k^2\sum_{d=1}^{\min(n,m)}{\rm sum}\left(\left\lfloor\frac n{dk}\right\rfloor\right){\rm sum}\left(\left\lfloor\frac m{dk}\right\rfloor\right)d

其中

{\rm sum}(n)=\sum_{i=1}^ni

然后……感觉枚举两层不太可做？别担心，我们有除法分块。众所周知

\left\lfloor\frac n{ab}\right\rfloor=\left\lfloor\frac {\left\lfloor\dfrac na\right\rfloor}{b}\right\rfloor

所以我们最内层的柿子可以先对 k 除法分块，在每一块内又对 d 除法分块。别忘记这一块的和要乘上这一块 \mu(k) 的总和。

看起来复杂度在 {\rm O}(n) 级别，预处理一下 \mu 的前缀和就可以算了。

不过，这里除法分块 d 的复杂度依赖于 k 这一块的块长，所以除法分块的复杂度还是：

{\rm O}\left(\sum_{i=1}^{\sqrt{n}}\sqrt\frac ni\right)={\rm O}\left(\int_{1}^{\sqrt{n}}\sqrt\frac nx{\rm d}x\right)={\rm O}(n^{\frac 34}) 
#include <iostream>
#include <bitset>

using namespace std;

long long   mu_square[10000001];
int         primes[664580];
bitset<10000001>            prime;
const long long             moder=20101009;

inline long long sum(long long num)
{
    return (num*num+num)%moder*10050505%moder;
}

int main(int argc,char *argv[],char *envp[])
{
    ios::sync_with_stdio(false);
    cin.tie(0);
    long long n,m,cntprime=0,lim;
    cin>>n>>m,lim=min(n,m);
    mu_square[1]=1;
    for(int i=2;i<=lim;i++)
    {
        if(!prime[i])
            primes[++cntprime]=i,mu_square[i]=moder-1ll*i*i%moder;
        for(int j=1;j<=cntprime&&i*primes[j]<=lim;j++)
        {
            prime[i*primes[j]]=true;
            if(!(i%primes[j]))
            {
                mu_square[i*primes[j]]=0;
                break;
            }
            mu_square[i*primes[j]]=mu_square[i]*mu_square[primes[j]]%moder;
        }
    }
    for(int i=2;i<=lim;i++)
        mu_square[i]=(mu_square[i]+mu_square[i-1])%moder;
    long long ans=0;
    for(long long left_r=1,right_r;left_r<=lim;left_r=right_r+1)
    {
        right_r=min(n/(n/left_r),m/(m/left_r));
        long long blc_szn=n/left_r,blc_szm=m/left_r,inner_ans=0;
        for(long long left_d=1,right_d;left_d<=lim/left_r;left_d=right_d+1)
        {
            right_d=min(blc_szn/(blc_szn/left_d),blc_szm/(blc_szm/left_d));
            inner_ans=(inner_ans+(mu_square[right_d]-mu_square[left_d-1]+moder)%moder*sum(n/left_r/left_d)%moder*sum(m/left_r/left_d)%moder)%moder;
        }
        ans=(ans+inner_ans*(sum(right_r)-sum(left_r-1)+moder)%moder)%moder;
    }
    cout<<ans;
    return 0;
}