「BZOJ2693」jzptab-莫比乌斯反演

2019.5.2

「BZOJ2693」jzptab-莫比乌斯反演

Description

求

$$\sum _{i=1}^N\sum _{j=1}^{M}lcm(i,j)$$

$n,m\le 10^7$，多组数据。

Solution

设

$$f(n)=\sum _{i=1}^N\sum _{j=1}^M[gcd(i,j)=1]ij$$
$$F(n)=\sum _{i=1}^N\sum _{j=1}^M[n|gcd(i,j)]ij$$

显然

$$F(n)=n^2\frac{\lfloor \frac{N}{n}\rfloor }{2}\frac{\lfloor \frac{M}{n}\rfloor }{2}$$

又因为

$$F(n)=\sum _{n|d}f(d)$$

所以

$$f(n)=\sum _{n|d}\mu (\frac{d}{n})F(d)$$

所以

$$\begin{array}{cc}& {\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)\\ =& {\sum }_{i=1}^n{\sum }_{j=1}^m\frac{ij}{gcd(i,j)}\\ =& {\sum }_{k=1}^n\frac{1}{k}f(k)\\ =& {\sum }_{k=1}^n\frac{1}{k}{\sum }_{k|d}\mu (\frac{d}{k})F(d)\\ =& {\sum }_{d}d\frac{\lfloor \frac{N}{d}\rfloor }{2}\frac{\lfloor \frac{M}{d}\rfloor }{2}{\sum }_{k|d}\mu (\frac{d}{k})\frac{d}{k}\\ =& {\sum }_{d}d\frac{\lfloor \frac{N}{d}\rfloor }{2}\frac{\lfloor \frac{M}{d}\rfloor }{2}{\sum }_{k|d}\mu (d)d\end{array}$$

${\sum }_{k|d}\mu (d)d$可以线性筛（可以发现每次添加一个已有的质因子时，新增加的约数的$\mu$为$0$，当添加一个没有的质因子$p$时，新增加的约数$d\cdot p$的贡献为$d\cdot p\mu (d\cdot p)=d\mu (d)\cdot p\mu (p)$）。

每次询问时数论分块求就好了。

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

inline int gi()
{
	char c = getchar();
	while(c < '0' || c > '9') c = getchar();
	int sum = 0;
	while('0' <= c && c <= '9') sum = sum * 10 + c - 48, c = getchar();
	return sum;
}

typedef long long ll;
const int maxn = 10000005, mod = 100000009;

int T, n, m, vis[maxn], p[maxn], tot, f[maxn], sum[maxn], s[maxn];

void pre()
{
	register int i, j;
	n = 1e7;
	f[1] = 1;
	for (i = 2; i <= n; ++i) {
		if (!vis[i]) p[++tot] = i, f[i] = mod + 1 - i;
		for (j = 1; j <= tot && i * p[j] <= n; ++j) {
			vis[i * p[j]] = 1;
			if (i % p[j]) f[i * p[j]] = (ll)f[i] * f[p[j]] % mod;
			else {f[i * p[j]] = f[i]; break;}
		}
	}
	for (i = 1; i <= n; ++i)
		sum[i] = (sum[i - 1] + (ll)i * f[i]) % mod, s[i] = s[i - 1] + i >= mod ? s[i - 1] + i - mod : s[i - 1] + i;
}

int main()
{
	freopen("jzptab.in", "r", stdin);
	freopen("jzptab.out", "w", stdout);

	int T = gi();
	register int i, j, ans;
	
	pre();
	
	while (T--) {
		n = gi(); m = gi();
		if (n > m) swap(n, m);
		ans = 0;
		for (i = 1; i <= n; i = j + 1) {
			j = min(n / (n / i), m / (m / i));
			ans = (ans + (ll)(sum[j] - sum[i - 1] + mod) * s[n / i] % mod * s[m / i]) % mod;
		}
		cout << ans << endl;
	}
	
	return 0;
}