「BZOJ」2154 Crash的数字表格-莫比乌斯反演

Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122

【数据规模和约定】

100%的数据满足N, M ≤ 10^7。

Solution

题目中如此之大的 $n$ 和 $m$ 已经将解法赤裸裸地暴露了出来……
最近学习了莫(meng)比(bi)乌(yi)斯(dun)反(luan)演(gao)，于是还是要写一篇题解麻痹自己。
回到这道题，这道题要求出一个很奇葩简单的值——
$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) \% mod$
其中 $mod=\dots$ 反正是一个很大的数就对了
乍看之下——毫无头绪，莫比乌斯反演似乎不太喜欢 $lcm$ 的样子……
我们有什么快速求出 $lcm$ 的式子吗？（先尝试 $O(n^2)$ 左右的暴力）
有的—— $lcm(a,b)=\cfrac{a·b}{gcd(a,b)}$
诶，这样一变，就出现了莫比乌斯反演很喜欢的 $gcd$
好，我们写下来
$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\cfrac {i·j}{gcd(i,j)}$
同时按照惯例我们要设 $n<m$
这样就是一个 $O(n^2\log n)$ 的暴力了
当然， $n≤10^7$ 是要干嘛的？跑一年都跑不完吧？！！
我们必须要优化。
考虑到很多数的 $gcd$ 很有可能相同，我们可以考虑把它们一起做
我们枚举他们的 $gcd$ —— $d$

$ans=\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}[gcd(i,j)=d]\cfrac {i·j}{d}$

这样似乎还把时间变长了，不急，我们不必枚举$i、j$，只需要枚举他们是$d$的多少倍就行了，不是么？以下的$i、j$代表倍数

$ans=\sum_{d=1}^{n}\sum_{i=1}^{\frac nd}\sum_{j=1}^{\frac md}[gcd(i,j)=1]i·j·d$

诶，我们又看到了一个很好玩的式子$gcd(i,j)=1$，我们知道，一个数$a=1$等价于$\sum_{d|a}\mu(d)=1$，所以我们再换一换（注意 $d$ 被提前了）

$ans=\sum_{d=1}^{n}d\sum_{i=1}^{\frac nd}\sum_{j=1}^{\frac md}i·j\sum_{k|gcd(i,j)}\mu(k)$

我们知道，一个数$k|gcd(i,j)$ 相当于$k|i$ 且 $k|j$，我们再次转化

$ans=\sum_{d=1}^{n}d\sum_{i=1}^{\frac nd}\sum_{j=1}^{\frac md}i·j\sum_{k=1}^{min(i,j)}[k|i][k|j]\mu(k)$

同样的，我们也可以枚举 $k$ 来代替枚举 $i、j$，然后枚举倍数即可

$ans=\sum_{d=1}^{n}d\sum_{k=1}^{\frac nd}\mu(k)\sum_{i=1}^{\frac {n}{kd}}\sum_{j=1}^{\frac {m}{kd}}i·j·k^2$

接下来，神奇的事情发生了，这样的$i、j$只有 $[\frac {n}{kd}]·[\frac {m}{kd}]$ 组，并且
$\sum i×j =(1+2+\dots+[\frac {n}{kd}])(1+2+\dots+[\frac {m}{kd}])=\cfrac{[\frac {n}{kd}][\frac {n}{kd}+1][\frac {m}{kd}][\frac {m}{kd}+1]}{4}$
这是可以 $O(1)$ 计算的，那么答案就是

$ans=4^{-1}\sum_{d=1}^{n}d\sum_{k=1}^{\frac nd}\mu(k)k^2[\frac {n}{dk}][\frac {m}{dk}]([\frac {n}{dk}]+1)([\frac {m}{dk}]+1)$

好了，现在降低成 $O(n\log n)$ （$n+\frac n2+\frac n3+\dots+\frac nn$ ≈ $n\log n$）了，然而依然无法接受

但是我们可以发现，这个式子可以使用分块优化

我们记 $f(n,m)=\sum_{k=1}^{n}\mu(k)k^2[\frac {n}{k}][\frac {m}{k}]([\frac {n}{k}]+1)([\frac {m}{k}]+1)$
我们再记 $g(n,m)=4^{-1}\sum_{d=1}^{n}d·f(\frac nd,\frac md)$
很容易就会发现，我们只需要求出 $ans=g(n,m)$，通过分块枚举 $\cfrac nd$ 和 $\cfrac md$，可以使得外层循环（$g$）降低为 $O(\sqrt n)$，同理，内层循环（$f$）也可以降低为 $O(\sqrt n)$，这样总时间复杂度就是 $O(n)$ 了
$Code:$

#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
#define N 10000000
#define M 5000000
#define mod 20101009
#define inv 15075757
int mu[N+5],h[N+5],p[M+5],cnt;
int sum[N+5],num[N+5],n,m;
int ans2(int a,int b)
{
    if(a>b)swap(a,b);
    int last;
    int re=0;
    for(int i=1;i<=a;i=last+1)
    {
        last=min(a/(a/i),b/(b/i));
        re=(re+1ll*(sum[last]-sum[i-1])*(a/i)%mod*(b/i)%mod*(a/i+1)%mod*(b/i+1)%mod)%mod;
    }
    return re;
}
int ans(int a,int b)
{
    if(a>b)swap(a,b);
    int last;
    int re=0;
    for(int i=1;i<=a;i=last+1)
    {
        last=min(a/(a/i),b/(b/i));
        re=(re+1ll*(num[last]-num[i-1])*ans2(a/i,b/i)%mod)%mod;
    }
    return (1ll*inv*re%mod+mod)%mod;
}
int main()
{
    mu[1]=1;
    scanf("%d%d",&n,&m);
    int mx=max(n,m);
    for(int i=2;i<=mx;i++)
    {
        if(!h[i])
        {
            mu[i]=-1;
            p[++cnt]=i;
        }
        for(int j=1;j<=cnt;j++)
        {
            if(p[j]*i>N)break;
            h[p[j]*i]=1;
            if(i%p[j]==0)break;
            mu[p[j]*i]=-mu[i];
        }
    }
    for(int i=1;i<=mx;i++)
    {
        sum[i]=(1ll*i*i*mu[i]%mod+sum[i-1])%mod;
        num[i]=(i+num[i-1])%mod;
    }
    printf("%d\n",ans(n,m));
}