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

题目概述

求 ∑ni=1∑mj=1lcm(i,j) 。

解题报告

第一道反演题……首先先要知道 e=μ∗1 ，证明戳这里。

莫比乌斯反演： f=g∗1⇔g=f∗μ ，证明：

f=g∗1⇔g∗(μ∗1)=f∗μ⇔g∗e=f∗μ⇔g=f∗μ

也就是说若 F(n)=∑d|nf(d) ，则 f(n)=∑d|nF(d)μ(nd) 。

这道题先列出答案式子：

∑i=1n∑j=1m[i,j]=∑i=1n∑j=1mijgcd(i,j)

令 f(n)=1n,F(n)=∑d|nf(d) ，则：

=∑i=1n∑j=1mij×f[gcd(i,j)]=∑i=1n∑j=1mij∑d|i,d|jF(d)=∑d=1min(n,m)F(d)∑i=1n∑j=1mij[d|i,d|j]=∑d=1min(n,m)F(d)d2∑i=1⌊nd⌋∑j=1⌊md⌋ij

后面的 ∑⌊nd⌋i=1i∑⌊md⌋j=1j 就是两个等差数列前 x 项和的乘积，我们记为 S(d) ，接下来就用到了反演：

=∑d=1min(n,m)F(d)d2S(d)=∑d=1min(n,m)S(d)∑k|df(dk)μ(k)d2=∑d=1min(n,m)S(d)∑k|dkd×d2μ(k)=∑d=1min(n,m)S(d)∑k|dμ(k)kd

我们会发现后面的 H(d)=∑k|dμ(k)kd 是积性函数，所以可以用线性筛：

H(p)=p−p2H(ap)=H(a)H(p),gcd(a,p)=1H(ap)=H(a)p,gcd(a,p)=p

简单证明： gcd(a,p)=p 时， H(ap) 比 H(a) 多出的 k 一定含有 p2 而使 μ(k)=0 ，相同 k 则 ×p ，所以 H(ap)=H(a)p 。

那么答案就是 ∑min(n,m)i=1S(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;
}