E - GuGuFishtion HDU - 6390（欧拉函数 / 莫比乌斯反演）

Posted 2020-11-10 jvruodejrLS

GuGuFishtion (HDU - 6390)

题意：

定义\(G_u (a,b)=\frac{\phi(ab)}{\phi(a)\phi(b)}\)。

求\((\sum\limits_{a=1}^m\sum\limits_{b=1}^nG_u (a,b))\pmod p\)。

题解：

考虑\(\phi(x) = x*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})...*(1-\frac{1}{p_n})\)。

将\(G_u (a,b)\)的分子与分母按上述分解、约分后，得\(\frac{1}{(1-\frac{1}{p_1})(1-\frac{1}{p_2})(1-\frac{1}{p_3})...(1-\frac{1}{p_k})}\)，其中\(p_1*p_2*p_3...p_k=gcd(a,b)\)。

上下同时乘以\(gcd(a,b)\)，则\(G_u (a,b)=\frac{gcd(a,b)}{\phi(gcd(a,b))}\)。

问题变成了求\((\sum\limits_{a=1}^m\sum\limits_{b=1}^n\frac{gcd(a,b)}{\phi(gcd(a,b))})\pmod p\)。

上式 \(=(\sum_d^{min(n, m)}\frac{d}{\phi(d)}\sum\limits_{a=1}^m\sum\limits_{b=1}^n[gcd(a,b)==d])\pmod p\)。

对于\(\sum\limits_{a=1}^m\sum\limits_{b=1}^n[gcd(a,b)==d]\)，即\(\sum\limits_{a=1}^{m/d}\sum\limits_{b=1}^{n/d}[gcd(a,b)==1]\)，可以用莫比乌斯反演求。

代码：

#include <bits/stdc++.h>
#define fopi freopen("in.txt", "r", stdin)
#define fopo freopen("out.txt", "w", stdout)
using namespace std;
typedef long long LL;
typedef long double ld;
const int maxn = 1e6 + 10;

int check[maxn], phi[maxn], mu[maxn], prime[maxn];
LL inv[maxn], sum[maxn];

void init(int N) {
    memset(check, false, sizeof(check));
    mu[1] = 1;
    phi[1] = 1;
    int tot = 0;
    for (int i = 2; i <= N; i++) {
        if (!check[i]) {
            prime[tot++] = i;
            mu[i] = -1;
            phi[i] = i-1;
        }
        for (int j = 0; j < tot; j++) {
            if (i * prime[j] > N) break;
            check[i * prime[j]] = true;
            if (i % prime[j] == 0) {
                mu[i * prime[j]] = 0;
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else {
                mu[i * prime[j]] = -mu[i];
                phi[i * prime[j]] = phi[i] * (prime[j]-1);
            }
        }
    }
}

LL solve(int n, int m, int p) {
    LL res = 0;
    for (int i = 1, la = 0; i <= m; i = la+1) {
        la = min(n/(n/i), m/(m/i));
        res = (res + 1ll * (sum[la] - sum[i-1]) * (n/i) % p * (m/i) % p) % p;
    }
    return res;
}

int T, n, m, p;
int main() {
    init(maxn - 10);
    
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d%d", &n, &m, &p);
        if (n < m) swap(n, m);

        inv[1] = 1;
        for (int i = 2; i <= m; i++)
            inv[i] = 1ll * inv[p % i] * (p - p/i) % p;
        for (int i = 1; i <= m; i++)
            sum[i] = (sum[i-1] + mu[i] + p) % p;

        LL ans = 0;
        for (int i = 1; i <= m; i++) {
            ans = (ans + 1ll * i * inv[phi[i]] % p * solve(n/i, m/i, p) % p) % p;
        }

        printf("%lld\n", ans);
    }
}