@hdu - 6390@ GuGuFishtion
@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\)。
于是：
可以发现这个式子与 A, B 的公共质因子有关，不妨考虑给它添几项，使它与 gcd(A, B) 相关。
问题很愉快地转为求 \((\sum_{i=1}^{m}\sum_{j=1}^{n}\frac{gcd(i, j)}{\phi(gcd(i, j))})\)。
出现了 gcd，考虑对式子进行反演。
根据调和级数，枚举 (d, p) 的复杂度为 O(nlog n)，但我常数大写不过。
一种方法是整除分块：考虑 \(\lfloor\frac{n}{p}\rfloor*\lfloor\frac{m}{p}\rfloor\) 的不同取值一共有 \(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\rfloor}\mu(a) - \sum_{a=1}^{\lfloor l/d\rfloor}\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\n", (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) 过的代码。
好嘛，欺负我常数大。
关于那个新的积性函数的递推式子，打一下表就可以发现了。
