【莫比乌斯反演】BZOJ2154[Crash的数字表格]题解

题目概述

求 ${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$ 。

解题报告

第一道反演题……首先先要知道 $e=\mu \ast 1$ ，证明戳这里。

莫比乌斯反演： $f=g\ast 1\iff g=f\ast \mu$ ，证明：
$$f=g\ast 1\iff g\ast (\mu \ast 1)=f\ast \mu \iff g\ast e=f\ast \mu \iff g=f\ast \mu$$
也就是说若 $F(n)={\sum }_{d|n}f(d)$ ，则 $f(n)={\sum }_{d|n}F(d)\mu (\frac{n}{d})$ 。

---

这道题先列出答案式子：
$$\sum _{i=1}^n\sum _{j=1}^m[i,j]=\sum _{i=1}^n\sum _{j=1}^m\frac{ij}{gcd(i,j)}$$
令 $f(n)=\frac{1}{n},f(n)={\sum }_{d|n}F(d)$ ，则：
$$\begin{gathered} =\sum _{i=1}^n\sum _{j=1}^{m}ij\times f[gcd(i,j)]=\sum _{i=1}^n\sum _{j=1}^{m}ij\sum _{d|i,d|j}F(d) \\ =\sum _{d=1}^{min(n,m)}F(d)\sum _{i=1}^n\sum _{j=1}^{m}ij[d|i,d|j]=\sum _{d=1}^{min(n,m)}F(d)d^2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }ij \end{gathered}$$
后面的 ${\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }i{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }j$ 就是两个等差数列前 $x$ 项和的乘积，我们记为 $S(d)$ ，接下来就用到了反演：
$$\begin{gathered} =\sum _{d=1}^{min(n,m)}F(d)d^{2}S(d)=\sum _{d=1}^{min(n,m)}S(d)\sum _{k|d}f(\frac{d}{k})\mu (k)d^2 \\ =\sum _{d=1}^{min(n,m)}S(d)\sum _{k|d}\frac{k}{d}\times d^2\mu (k)=\sum _{d=1}^{min(n,m)}S(d)\sum _{k|d}\mu (k)kd \end{gathered}$$
我们会发现后面的 $H(d)={\sum }_{k|d}\mu (k)kd$ 是积性函数，所以可以用线性筛：
$$\begin{gathered} H(p)=p-p^2 \\ H(ap)=H(a)H(p),gcd(a,p)=1 \\ H(ap)=H(a)p,gcd(a,p)=p \end{gathered}$$

简单证明： $gcd(a,p)=p$ 时， $H(ap)$ 比 $H(a)$ 多出的 $k$ 一定含有 $p^2$ 而使 $\mu (k)=0$ ，相同 $k$ 则 $\times p$ ，所以 $H(ap)=H(a)p$ 。

那么答案就是 ${\sum }_{i=1}^{min(n,m)}S(d)H(d)$ ，其实可以分块加快速度，不过我太懒就没写:P。

示例程序

#include<cstdio>
using namespace std;
typedef long long LL;
const int maxn=1e7,MOD=20101009;

int n,m,H[maxn+5],ans;
int p[maxn+5];bool pri[maxn+5];

inline void AMOD(int &x,int tem) {if ((x+=tem)>=MOD) x-=MOD;}
void Make(int n)
{
	pri[1]=true;H[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!pri[i]) p[++p[0]]=i,AMOD(H[i]=i,MOD-(LL)i*i%MOD);
		for (int j=1,t;j<=p[0]&&(t=i*p[j])<=n;j++)
		{
			pri[t]=true;H[t]=(LL)H[i]*H[p[j]]%MOD;
			if (!(i%p[j])) {H[t]=(LL)H[i]*p[j]%MOD;continue;}
		}
	}
}
inline LL S(int d) {LL A=n/d,B=m/d;return ((A+1)*A/2)%MOD*((B+1)*B/2%MOD)%MOD;}
int main()
{
	freopen("program.in","r",stdin);
	freopen("program.out","w",stdout);
	scanf("%d%d",&n,&m);Make(n);
	for (int i=1;i<=n&&i<=m;i++) AMOD(ans,S(i)*H[i]%MOD);
	return printf("%d\n",ans),0;
}