【国家集训队】Crash 的数字表格 / JZPTAB

前言

以前一直都懒得写莫比乌斯反演题的题解。。。（因为字数太多而且自己是蒟蒻）

不过正因为是蒟蒻才要写题解不是吗啊啊啊希望写完后可以熟练掌握markdown滑稽

Solution

首先我们写出答案的式子（我们令 (mathtt{n <= m})）：

[mathtt{sum_{i=1}^nsum_{j=1}^m lcm(i,j)} ]

然后我们很自然地就想到了将 lcm 转化为 gcd。（毕竟公式是这么来的）

[mathtt{sum_{i=1}^nsum_{j=1}^m frac{i*j}{gcd(i,j)}} ]

枚举 i 与 j 的 gcd。（P.S.(mathtt{[x=1]}) 指的是当 x 取 1 时表达式的值为 1，其余为 0）

[mathtt{sum_{k=1}^nsum_{i=1}^nsum_{j=1}^m [gcd(i,j)=k]*frac{i*j}{k}} ]

把 k 提出来。

[mathtt{sum_{k=1}^nk*sum_{i=1}^{frac nk}sum_{j=1}^{frac mk} [gcd(i,j)=1]*i*j} ]

我们讲一讲这个式子是怎么化过去的。我们知道，如果 (mathtt{gcd(i*k,j*k)==k})，那么其实可以转化成 (mathtt{gcd(i,j)==1})。可不可以这样理解：在 (mathtt{[1,n]}) 与 (mathtt{[1,m]}) 里分别选取 x，y，使其 gcd 为 k。在 (mathtt{[1,n/k]}) 与 (mathtt{[1,m/k]}) 中选择 i，j 使其互质。那么 (mathtt{i*k})，(mathtt{j*k}) 和 (mathtt{x})，(mathtt{y}) 其实是一一对应的关系。因为 (mathtt{n/k}) 就是在 n 里面选取 k 的倍数，我们这样改写不会错过任何一个 gcd 为 k 的解。

至于后面的 (mathtt{i*j}) 可以这样理解：通过前面的式子，我们知道后半部分其实就是 lcm。我们尝试列出此时的 lcm：(mathtt{frac{(i*k)*(j*k)}{k}})。

即 (mathtt{i*j*k})，将 k 提前，就是 (mathtt{i*j})。

好了接下来我们利用莫比乌斯函数本身性质可以推出：

[mathtt{sum_{k=1}^nk*sum_{i=1}^{frac nk}sum_{j=1}^{frac mk}sum_{d|gcd(i,j)} μ(d)*i*j} ]

枚举 d，就是转化成 (mathtt{μ(d)}) 可以匹配什么。。。至于为什么 d 的上界是 (mathtt{frac nk})，我们将 (mathtt{gcd(i,j)}) 都除了个 k 嘛。

[mathtt{sum_{k=1}^nk*sum_{d=1}^{frac nk}μ(d)*sum_{i=1}^{frac nk}sum_{j=1}^{frac mk}[d|gcd(i,j)]*i*j} ]

我们枚举 i 或 j 倍的 d，把 (mathtt{i*k}) 和 (mathtt{j*k}) 看作原来的 i，j。因为肯定满足条件，所以直接乘起来。然后可以推出这样一个鬼式子（i，j 均为整）：

[mathtt{sum_{k=1}^nk*sum_{d=1}^{frac nk}μ(d)*sum_{i*d=1}^{frac nk}sum_{j*d=1}^{frac mk}(i*d)*(j*d)} ]

显然下一步可以化成

[mathtt{sum_{k=1}^nk*sum_{d=1}^{frac nk}μ(d)*sum_{i=1}^{frac {n}{k*d}}sum_{j=1}^{frac {m}{k*d}}(i*d)*(j*d)} ]

继续，

[mathtt{sum_{k=1}^nk*sum_{d=1}^{frac nk}μ(d)*d^2sum_{i=1}^{frac {n}{k*d}}sum_{j=1}^{frac {m}{k*d}}i*j} ]

用等差数列求和公式：

[mathtt{sum_{k=1}^nk*sum_{d=1}^{frac nk}μ(d)*d^2*frac{frac{n}{k*d}*(frac{n}{k*d}+1)}{2}*frac{frac{m}{k*d}*(frac{m}{k*d}+1)}{2}} ]

最后我们筛出 (mathtt{μ(d)*d^2})，套两层数论分块就行了QwQ。

终于写完了。。。

emm 感觉证了半天自己对莫比乌斯反演一点感觉都没有我果然是蒟蒻。

Code

#include <cstdio>
#include <cmath>
#include <iostream>
using namespace std;
#define rep(i,_l,_r) for(signed i=(_l),_end=(_r);i<=_end;++i)
#define print(x,y) write(x),putchar(y)

#define int long long

int read() {
	int x=0,f=1; char s;
	while((s=getchar())>'9'||s<'0') if(s=='-') f=-1;
	while(s>='0'&&s<='9') x=(x<<1)+(x<<3)+(s^48),s=getchar();
	return x*f;
}

void write(int x) {
	if(x<0) return (void)(putchar('-'),write(-x));
	if(x>9) write(x/10);
	putchar(x%10^48);
}

const int maxn=1e7+5,mod=20101009;

int n,m,mu[maxn],p[(int)1e6],pc,ans,inv;
bool is[maxn];

int get(int x) {
	return 1ll*(x+1)*x%mod*inv%mod;
}

int Duck(int N,int M) {
	int r,ret=0,lim=min(N,M);
	for(int l=1;l<=lim;l=r+1) {
		r=min(N/(N/l),M/(M/l));
		ret=(ret+1ll*(mu[r]-mu[l-1]+mod)%mod*get(N/l)%mod*get(M/l)%mod)%mod;
	}
	return ret;
}

void init() {
	mu[1]=1;
	rep(i,2,maxn-5) {
		if(!is[i]) p[++pc]=i,mu[i]=-1;
		rep(j,1,pc) {
			if(p[j]*i>maxn-5) break;
			is[p[j]*i]=1;
			if(i%p[j]==0) break;
			mu[p[j]*i]=-mu[i];
		}
	}
	rep(i,1,maxn-5) mu[i]=(0ll+mu[i-1]+1ll*(mu[i]+mod)%mod*i%mod*i%mod)%mod;
}

int Goose() {
	int lim=min(n,m),r,ret=0;
	for(int l=1;l<=lim;l=r+1) {
		r=min(n/(n/l),m/(m/l));
		ret=(ret+1ll*(l+r)*(r-l+1)%mod*inv%mod*Duck(n/l,m/l)%mod)%mod;
	}
	return ret;
}

int qkpow(int x,int y) {
	int r=1;
	while(y) {
		if(y&1) r=1ll*r*x%mod;
		x=1ll*x*x%mod; y>>=1;
	}
	return r;
}

signed main() {
	inv=qkpow(2,mod-2);
	n=read(),m=read();
	init();
	printf("%lld
",Goose());
	return 0;
}