@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\)。

于是：
\[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)}*\lfloor\frac{n}{p}\rfloor*\lfloor\frac{m}{p}\rfloor\]

根据调和级数，枚举 (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) 过的代码。
好嘛，欺负我常数大。

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

@hdu - 6390@ GuGuFishtion

原文：https://www.cnblogs.com/Tiw-Air-OAO/p/11299994.html