@hdu

目录

• @description@
• @solution@
• @accepted code@
• @details@

---

@description@

定义函数 (G_u(a, b) = frac{phi(ab)}{phi(a)phi(b)})。给定 n, m 与质数 p，求解 ((sum_{i=1}^{m}sum_{j=1}^{n}G_u(i, j)) mod p)。

Input
第一行一个整数 T 描述数据组数。
接下来 T 行每行三个整数 m, n, p，含义如上。
1≤T≤3，1≤m,n≤1,000,000，max(m,n)<p≤1,000,000,007 且 p 为质数。

Output
输出 T 行，表示每组数据的答案。

Sample Input
1
5 7 23
Sample Output
2

@solution@

考虑 (phi(x))的定义式，设 (x) 的唯一分解式为 (x = prod_{i=0}^{m}p_i^{x_i})，则 (phi(x) = prod_{i=0}^{m}p_i^{x_i-1}*(p_i - 1))。

将这个定义式代入函数 (G_u) 中，不妨设 (A = A'*prod_{i=0}^{m}p_i^{a_i}, B = B'*prod_{i=0}^{m}p_i^{b_i})，其中要求 (gcd(A', B') = gcd(A', p_i) = gcd(B', p_i) = 1)。

于是：

[G_u(A, B) = frac{phi(AB)}{phi(A)phi(B)} \ = frac{phi(A')phi(B')prod_{i=0}^{m}p_i^{a_i+b_i-1}*(p_i - 1)}{phi(A')phi(B')prod_{i=0}^{m}p_i^{a_i+b_i-2}*(p_i - 1)^2}\ = frac{prod_{i=0}^{m}p_i}{prod_{i=0}^{m}(p_i - 1)}]

可以发现这个式子与 A, B 的公共质因子有关，不妨考虑给它添几项，使它与 gcd(A, B) 相关。

[G_u(A, B) = frac{prod_{i=0}^{m}p_i}{prod_{i=0}^{m}(p_i - 1)}\ = frac{prod_{i=0}^{m}p_i^{min(a_i, b_i)-1}*p_i}{prod_{i=0}^{m}p_i^{min(a_i, b_i)-1}*(p_i - 1)}\ = frac{gcd(A, B)}{phi(gcd(A, B))}]

问题很愉快地转为求 ((sum_{i=1}^{m}sum_{j=1}^{n}frac{gcd(i, j)}{phi(gcd(i, j))}))。

出现了 gcd，考虑对式子进行反演。

[ans = (sum_{i=1}^{m}sum_{j=1}^{n}frac{gcd(i, j)}{phi(gcd(i, j))})\ = sum_{d=1}sum_{d|p}mu(frac{p}{d})*frac{d}{phi(d)}*lfloorfrac{n}{p} floor*lfloorfrac{m}{p} floor]

根据调和级数，枚举 (d, p) 的复杂度为 O(nlog n)，但我常数大写不过。

一种方法是整除分块：考虑 (lfloorfrac{n}{p} floor*lfloorfrac{m}{p} floor) 的不同取值一共有 (O(sqrt{n})) 种，可以枚举出来算出对应的 p 的区间。
则我们相当于是要快速求出 (sum_{p=l}^{r}sum_{d|p}mu(frac{p}{d})*frac{d}{phi(d)})。
再进行变形得到 (sum_{d=1}^{r}frac{d}{phi(d)}*(sum_{a=1}^{lfloor r/d floor}mu(a) - sum_{a=1}^{lfloor l/d floor}mu(a)))。
可以发现又出现了整除，于是我们又可以分块用 (O(sqrt{n})) 的时间计算上式，总时间复杂度降为 (O(n))。

另一种方法，我们考虑使用线性筛。因为 ((sum_{d|p}mu(frac{p}{d})*frac{d}{phi(d)})) 是 p 的积性函数，所以使用线性筛 O(n) 求出后直接枚举 p 计算即可。

代码给出的是线性筛的实现。

@accepted code@

#include<cstdio>
#include<algorithm>
using namespace std;
const int MAXN = 1000000;
int prm[MAXN + 5], pcnt = 0;
int phi[MAXN + 5], miu[MAXN + 5];
bool nprm[MAXN + 5];
void init() {
	miu[1] = phi[1] = 1;
	for(int i=2;i<=MAXN;i++) {
		if( !nprm[i] ) {
			prm[++pcnt] = i;
			miu[i] = -1;
			phi[i] = i-1;
		}
		for(int j=1;1LL*i*prm[j]<=MAXN;j++) {
			nprm[i*prm[j]] = true;
			if( i % prm[j] == 0 ) {
				phi[i*prm[j]] = phi[i]*prm[j];
				miu[i*prm[j]] = 0;
				break;
			}
			else {
				phi[i*prm[j]] = phi[i]*phi[prm[j]];
				miu[i*prm[j]] = miu[i]*miu[prm[j]];
			}
		}
	}
}
int n, m, p;
int inv[MAXN + 5], f[MAXN + 5], g[MAXN + 5];
void solve() {
	int ans = 0;
	scanf("%d%d%d", &m, &n, &p);
	if( n > m ) swap(n, m);
	inv[1] = 1;
	for(int i=1;i<=n;i++) {
		if( i == 1 ) inv[i] = 1;
		else inv[i] = (p - 1LL*inv[p%i]*(p/i)%p)%p;
		f[i] = 1LL*i*inv[phi[i]]%p;
	}
	g[1] = 1;
	for(int i=2;i<=n;i++) {
		if( !nprm[i] ) g[i] = (1LL*miu[i]*f[1] + 1LL*miu[1]*f[i])%p;
		for(int j=1;1LL*i*prm[j]<=n;j++) {
			nprm[i*prm[j]] = true;
			if( i % prm[j] == 0 ) {
				g[i*prm[j]] = 0;
				break;
			}
			else g[i*prm[j]] = 1LL*g[i]*g[prm[j]]%p;
		}
	}
	for(int j=n;j>=1;j--)
		ans = (ans + 1LL*(n/j)*(m/j)%p*g[j]%p)%p;
	printf("%d
", (ans + p)%p);
}
int main() {
	init();
	int T; scanf("%d", &T);
	for(int i=1;i<=T;i++)
		solve();
}

@details@

卡我 O(nlog n) 还行。关键是我上网翻题解竟然翻到 O(nlog n) 过的代码。
好嘛，欺负我常数大。

关于那个新的积性函数的递推式子，打一下表就可以发现了。