題目

求:

\[ \sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) \in \mathbb P] \]

有 \(T\) 組數據, \(T\le 10^4, n, m\le 10^7\)

分析

莫比烏斯反演:

\[ \begin{align*} & \sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) \in \mathbb P]\ = & \sum_{p \in \mathbb P, p\le \min(n, m)}\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = p]\\end{align*} \]

設 \(f(x) = \sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = x]\), $F(x) = \sum_{i = 1}^n \sum_{j = 1}^m [x?|\gcd(i, j)]=\left\lfloor \frac nx\right\rfloor\left\lfloor \frac mx\right\rfloor $

則有:
\[ F(x) = \sum_{x|d, d\le \min(n, m)} f(d) \Rightarrow f(x) = \sum_{x|d, d\le\min(n, m)} \mu(\frac dx)F(d) \]

求法(暴力出奇跡, 經測試下面算法耗時不到 \(0.5 s\)):

void get_g(int n) {
    for(int i = 1; i <= n; i++) {
        int tmp = i;
        while(tmp != 1) {
            g[i] += mobius[i / s_factor[tmp]];
            int p = s_factor[tmp];
            if(tmp % (p * p) == 0) break;
            for(; tmp % p == 0; tmp /= p);
        }
        g_prefix[i] = g_prefix[i - 1] + g[i];
    }
}

有上式等於:
\[ \sum_{d = 1}^{\min(n, m)} \left\lfloor \frac nd\right\rfloor\left\lfloor \frac md\right\rfloor g(d) \]
對於 \(\left\lfloor \frac nd\right\rfloor\left\lfloor \frac md\right\rfloor\) 相同的值整除分塊即可.

代碼

#include <bits/stdc++.h>

typedef long long Int64;

const int kMaxSize = 1e7 + 5;

int s_factor[kMaxSize], prime[kMaxSize], mobius[kMaxSize], g[kMaxSize],
    g_prefix[kMaxSize], prime_tot = 0;
bool isn_prime[kMaxSize];

void euler_sieve(int n) {
    mobius[1] = 1;
    isn_prime[0] = isn_prime[1] = true;
    for(int i = 2; i <= n; i++) {
        if(!isn_prime[i]) {
            prime[prime_tot++] = i;
            s_factor[i] = i;
            mobius[i] = -1;
        }
        for(int j = 0; j < prime_tot && i * prime[j] <= n; j++) {
            isn_prime[i * prime[j]] = true;
            s_factor[i * prime[j]] = prime[j];
            mobius[i * prime[j]] = -mobius[i];
            if(i % prime[j] == 0) {
                mobius[i * prime[j]] = 0;
                break;
            }
        }
    }
}

void get_g(int n) {
    for(int i = 1; i <= n; i++) {
        int tmp = i;
        while(tmp != 1) {
            g[i] += mobius[i / s_factor[tmp]];
            int p = s_factor[tmp];
            if(tmp % (p * p) == 0) break;
            for(; tmp % p == 0; tmp /= p);
        }
        g_prefix[i] = g_prefix[i - 1] + g[i];
    }
}

int main() {
    euler_sieve(1e7);
    get_g(1e7);
    int t;
    scanf("%d", &t);
    while(t--) {
        int n, m;
        Int64 ans = 0;
        scanf("%d%d", &n, &m);
        if(!n) break;
        for(int l = 1, r; l <= n && l <= m; l = r + 1) {
            r = std::min(n / (n / l), m / (m / l));
            ans += (Int64)(n / l) * (m / l) * (g_prefix[r] - g_prefix[l - 1]);
        }
        printf("%lld\n", ans);
    }
    return 0;
}

【Luogu P2257】YY 的 GCD