题目

求:

\[ \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) \]

故有:

\[ \begin{align*} & \sum_{p \in \mathbb P, p\le \min(n, m)}\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = p]\\ = & \sum_{p \in \mathbb P, p\le \min(n, m)} f(p) \\ = & \sum_{p \in \mathbb P, p\le \min(n, m)} \sum_{p|d,d\le \min(n, m)} \mu(\frac dp) \left\lfloor \frac nd\right\rfloor\left\lfloor \frac md\right\rfloor \\ = & \sum_{d = 1}^{\min(n, m)} \sum_{p \in \mathbb P,p|d, p\le \min(n, m)} \mu(\frac dp) \left\lfloor \frac nd\right\rfloor\left\lfloor \frac md\right\rfloor \\ = & \sum_{d = 1}^{\min(n, m)} \left\lfloor \frac nd\right\rfloor\left\lfloor \frac md\right\rfloor \sum_{p \in \mathbb P,p|d, p\le \min(n, m)} \mu(\frac dp) \end{align*} \]

设 $g(x) = \sum_{p \in \mathbb P,p|x, p\le \min(n, m)} \mu(\frac xp) $

求法(暴力出奇迹, 经测试下面算法耗时不到 \(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 
    
     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;}