[BZOJ2693]jzptab

[BZOJ2693]jzptab

试题描述

给出 \(n, m\),求 \(\sum_{i=1}^n \sum_{j=1}^m \mathrm{lcm}(i, j)\)\(100000009\) 取模后的结果。

多组询问。

输入

一个正整数 \(T\) 表示数据组数

接下来 \(T\) 行,每行两个正整数,表示 \(n\)\(m\)

输出

\(T\) 行,每行一个整数,表示第 \(i\) 组数据的结果

输入示例

1
4 5

输出示例

122

数据规模及约定

\(T \le 10000\)

\(n, m \le 10000000\)

题解

这道题是那道题的加强版,强化的地方是多组数据而不是 \(n, m\) 的范围,所以那篇文章后面的优化做法是没有用的。

那篇博文的公式有点丑,所以在这里重推一下……

\[\sum_{i=1}^n \sum_{j=1}^m \mathrm{lcm}(i, j) \\ = \sum_{i=1}^n \sum_{j=1}^m \frac{ij}{\gcd(i, j)} \]

现在考虑枚举最大公约数

\[= \sum_{g=1}^{\min \{n, m\} } \frac{1}{g} \sum_{i=1}^n i \sum_{j=1}^m j \cdot [\gcd(i, j) = g] \\ = \sum_{g=1}^{\min \{n, m\} } \frac{1}{g} \sum_{i=1}^{\lfloor \frac{n}{g} \rfloor} ig \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} jg \sum_{d|i, d|j} \mu(d) \\ = \sum_{g=1}^{\min \{n, m\} } g \sum_{d=1}^{\lfloor \frac{\min \{n, m\}}{g} \rfloor} \mu(d) \sum_{i=1}^{\lfloor \frac{n}{gd} \rfloor} id \sum_{j=1}^{\lfloor \frac{m}{gd} \rfloor} jd \\ = \sum_{g=1}^{\min \{n, m\} } g \sum_{d=1}^{\lfloor \frac{\min \{n, m\}}{g} \rfloor} \mu(d) \cdot d^2 \frac{\lfloor \frac{n}{gd} \rfloor \left( \lfloor \frac{n}{gd} \rfloor + 1 \right)}{2} \cdot \frac{\lfloor \frac{m}{gd} \rfloor \left( \lfloor \frac{m}{gd} \rfloor + 1 \right)}{2} \\ = \frac{1}{4} \sum_{g=1}^{\min \{n, m\} } g \sum_{d=1}^{\lfloor \frac{\min \{n, m\}}{g} \rfloor} \mu(d) \cdot d^2 \lfloor \frac{n}{gd} \rfloor \left( \lfloor \frac{n}{gd} \rfloor + 1 \right) \lfloor \frac{m}{gd} \rfloor \left( \lfloor \frac{m}{gd} \rfloor + 1 \right) \]

从这里开始将和那道题不一样了,现在用 \(T\) 替换 \(gd\)

\[= \frac{1}{4} \sum_{T=1}^{\min \{n, m\} } \lfloor \frac{n}{T} \rfloor \left( \lfloor \frac{n}{T} \rfloor + 1 \right) \lfloor \frac{m}{T} \rfloor \left( \lfloor \frac{m}{T} \rfloor + 1 \right) \sum_{d|T} \mu(d) \cdot d^2 \cdot \frac{T}{d} \\ = \frac{1}{4} \sum_{T=1}^{\min \{n, m\} } \lfloor \frac{n}{T} \rfloor \left( \lfloor \frac{n}{T} \rfloor + 1 \right) \lfloor \frac{m}{T} \rfloor \left( \lfloor \frac{m}{T} \rfloor + 1 \right) T \sum_{d|T} \mu(d) \cdot d \]

\(f(T) = T \sum_{d|T} \mu(d) \cdot d\),则 \(f(T)\) 是积性函数(两个积性函数的乘积是积性函数,两个积性函数的狄利克雷卷积也是积性函数),我们可以线性筛出来。

如何线性筛?假设 \(p\) 是质数,则 \(f(p) = p(1 - p)\),我们要解决的就是已知 \(f(T)\)\(f(Tp)\),分两种情况:

  • \(p|T\),令 \(T = T' p^k\),且 \(T'\) 中不含因子 \(p\),那么 \(f(Tp) = f(T') \cdot f(p^{k+1})\),在线性筛的过程 \(f(T')\) 已知,\(f(p^k)\) 已知,那么 \(f(p^{k+1}) = f(p^k) \cdot p\),这是因为 \(f(p^k) = p^k \sum_{i=0}^k \mu(p^i) p^i\),后面 \(\mu(p^i)(i > 1)\) 都是 \(0\),所以只有前面那个 \(p^k\) 会乘上一个 \(p\),所以最后得到 \(f(Tp) = f(T) \cdot p\)
  • \(p \nmid T\),显然 \(f(Tp) = f(T) \cdot f(p)\)

这样我们可以 \(O(n)\) 的时间求出 \(f(T)\) 的前缀和,后面的询问就变成 \(O(\sqrt n)\) 一次啦。

注意 \(100000009\) 不是素数,但这并不影响我们求 \(4\) 的逆元。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;

#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

int read() {
	int x = 0, f = 1; char c = getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = getchar(); }
	return x * f;
}

#define maxn 10000010
#define MOD 100000009
#define inv4 75000007
#define LL long long

bool vis[maxn];
int prime[maxn], cp, f[maxn], sumf[maxn];

void init() {
	int n = maxn - 1;
	f[1] = sumf[1] = 1;
	rep(i, 2, n) {
		if(!vis[i]) prime[++cp] = i, f[i] = (LL)i * (1 - i + MOD) % MOD;
		for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
			vis[i*prime[j]] = 1;
			if(i % prime[j] == 0) {
				f[i*prime[j]] = (LL)f[i] * prime[j] % MOD;
				break;
			}
			f[i*prime[j]] = (LL)f[i] * f[prime[j]] % MOD;
		}
		sumf[i] = sumf[i-1] + f[i];
		if(sumf[i] >= MOD) sumf[i] -= MOD;
	}
	return ;
}

void work() {
	int n = read(), m = read(), ans = 0;
	for(int T = 1; T <= min(n, m); ) {
		int r = min(min(n / (n / T), m / (m / T)), min(n, m));
		ans += (LL)(n / T) * (n / T + 1) % MOD * (m / T) % MOD * (m / T + 1) % MOD * (sumf[r] - sumf[T-1] + MOD) % MOD;
		if(ans >= MOD) ans -= MOD;
		T = r + 1;
	}
	ans = (LL)ans * inv4 % MOD;
	printf("%d\n", ans);
	return ;
}

int main() {
	init();

	int T = read();
	while(T--) work();

	return 0;
}
posted @ 2018-02-23 15:03  xjr01  阅读(157)  评论(0编辑  收藏  举报