[莫反] P1829 Crash的数字表格

posted on 2024-03-09 05:21:06 | under | source

前言：自己推的，和部分题解（包括当前第一篇）不太一样。

为了方便令 \(n\le m\)，下文的除号表示向下取整，并且和式都默认从 \(1\) 开始枚举。

问题是求 \(\sum\limits_i^n\sum\limits_j^mlcm(i,j)\)。

转成熟悉的公因数形式：

\(\sum\limits_i^n\sum\limits_j^m\frac {ij}{gcd(i,j)}\)

枚举公因数，并将 \(i,j\) 同时除以 \(d\)：

\(\sum\limits_d^n\frac1d\sum\limits_i^{\frac nd}\sum\limits_j^{\frac md}id\times jd [gcd(i,j)==1]\)

\(\sum\limits_d^nd\sum\limits_i^{\frac nd}\sum\limits_j^{\frac md}ij [gcd(i,j)==1]\)

莫反，这里跳了一步：

\(\sum\limits_d^nd\sum\limits_i^{\frac nd}\sum\limits_j^{\frac md}ij\sum\limits_{p\mid i,p\mid j} \mu(p)\)

将后半部分的 \(\mu\) 提上来，并将 \(i,j\) 同时除以 \(p\)：

\(\sum\limits_d^nd\sum\limits_p^n\mu(p)\sum\limits_i^{\frac n{dp}}\sum\limits_j^{\frac m{dp}}ip\times jp\)

把能提出来变量的都提:

\(\sum\limits_d^nd\sum\limits_p^np^2\mu(p)\sum\limits_i^{\frac n{dp}}i\sum\limits_j^{\frac m{dp}}j\)

令 \(G(n)=\sum\limits_{i}^n i\)，化简原式：

\(\sum\limits_d^nd\sum\limits_p^np^2\mu(p)G(\frac n{dp})G(\frac m{dp})\)

枚举 \(dp=T\)：

\(\sum\limits_T^n G(\frac nT)G(\frac mT) \sum\limits_{d\mid T} d^2\mu(d)\frac Td\)

最终式子长这样：

\(\sum\limits_T^n G(\frac nT)G(\frac mT) \sum\limits_{d\mid T} dT\mu(d)\)

后半部分的和式显然是积性函数，线性筛预处理前缀和，然后整除分块即可。复杂度 \(O(n+\sqrt n)\)

#include<bits/stdc++.h>
using namespace std;

#define int long long
const int N = 1e7 + 5, mod = 20101009;
int T, n, m, k, ans;
int mi[N], f[N], ts[N], mu[N], prime[N], pc, tot;
bool flag[N];

inline int qstp(int a, int k) {int res = 1; for(; k; k >>= 1, a = a * a % mod) if(k & 1) res = res * a % mod; return res;}
inline void init(){
	f[1] = 1, flag[1] = 1, mu[1] = 1;
	for(int i = 2; i < N; ++i){
		if(!flag[i]) prime[++pc] = i, mi[i] = i, ts[i] = i, mu[i] = -1, f[i] = (i * (1 + mu[i] * i) + mod) % mod;
		for(int j = 1; j <= pc && i * prime[j] < N; ++j){
			int q = i * prime[j];
			flag[q] = 1, mi[q] = prime[j], ts[q] = prime[j] * (mi[i] == prime[j] ? ts[i] : 1);
			mu[q] = i % prime[j] == 0 ? 0 : -mu[i];
			int A = ts[q], B = q / A;
			if(B == 1)
				f[q] = (q * (1 + mu[prime[j]] * prime[j]) % mod + mod) % mod;
			else
				f[q] = f[A] * f[B] % mod;
			if(i % prime[j] == 0) break; 
		} 
	}
	for(int i = 2; i < N; ++i) f[i] = (f[i] + f[i - 1]) % mod;
}
inline int G(int a) {return (a + 1) * a / 2 % mod;}
inline int getF(int a){
	int res = 0;
	for(int d = 1; d <= a; ++d){
		if(a % d) continue;
		res += mu[d] * d * a;
	}
	return ((res % mod) + mod) % mod;
}
signed main(){
	init();
	cin >> n >> m;
	if(n > m) swap(n, m);
	for(int l = 1; l <= n;){
		int r = min(n / (n / l), m / (m / l));
		ans = (ans + G(n / l) * G(m / l) % mod * ((f[r] - f[l - 1]) % mod + mod) % mod) % mod;
		l = r + 1;
	}
	cout << ans;
	return 0;
}