@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),但我常数大写不过。
一种方法是整除分块:考虑 (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) 过的代码。
好嘛,欺负我常数大。
关于那个新的积性函数的递推式子,打一下表就可以发现了。