题目链接：LOJ 2000

Doris 刚刚学习了 $\text{Fibnacci}$ 数列，用 $f[i]$ 表示数列的第 $i$ 项，那么：

$$ \begin{align*} f[0] &= 0 \\ f[1] &= 1 \\ f[n] &= f[n - 1] + f[n - 2], n \ge 2 \end{align*} $$

Doris 用老师的超级计算机生成了一个 $n \times m$ 的表格，第 $ i $ 行第 $ j $ 列的格子中的数是 $f[\gcd(i, j)]$，其中 $\gcd(i, j)$ 表示 $i$ 与 $j$ 的最大公约数。

Doris 的表格中共有 $ n \times m $ 个数，她想知道这些数的乘积是多少。

这些数的乘积实在是太大了，所以 Doris 只想知道乘积对 $10 ^ 9 + 7$ 取模后的结果。

本题有 $T$ 组数组。

数据范围：$1 \le T \le 1000$，$1 \le n, m \le 10 ^ 6$。

Solution

直接列出答案的式子：

$$ \prod_{i = 1} ^ n \prod_{j = 1} ^ m f[\gcd(i, j)] $$

按照套路，我们枚举 $\gcd$ 化为：

$$ \prod_{d = 1} ^ {\min(n, m)} f[d] ^ {\sum_{i = 1} ^ n \sum_{j = 1} ^ m [\gcd(i, j) = d]} \\ \prod_{d = 1} ^ {\min(n, m)} f[d] ^ {\sum_{i = 1} ^ {\min\left(\left\lfloor\frac{n}{d}\right\rfloor, \left\lfloor\frac{m}{d}\right\rfloor\right)} \mu(i) \cdot \left\lfloor\frac{n}{id}\right\rfloor \cdot \left\lfloor\frac{m}{id}\right\rfloor} \\ \prod_{d = 1} ^ {\min(n, m)} f[d] ^ {\sum_{d \mid T} \mu\left(\frac{T}{d}\right) \cdot \left\lfloor\frac{n}{T}\right\rfloor \cdot \left\lfloor\frac{m}{T}\right\rfloor} \\ \prod_{T = 1} ^ {\min(n, m)} \prod_{d \mid T} f[d] ^ {\mu\left(\frac{T}{d}\right) \cdot \left\lfloor\frac{n}{T}\right\rfloor \cdot \left\lfloor\frac{m}{T}\right\rfloor} \\ \prod_{i = 1} ^ {\min(n, m)} \left(\prod_{d \mid i} f[d] ^ {\mu\left(\frac{i}{d}\right)}\right) ^ {\left\lfloor\frac{n}{i}\right\rfloor \cdot \left\lfloor\frac{m}{i}\right\rfloor} $$

我们记 $F(n) = \prod_{d \mid n} f[d] ^ {\mu\left(\frac{n}{d}\right)}$，那么 $F(1) \dots F(n)$ 可以在 $O(n \log n)$ 的时间内预处理得到。

原式化为：

$$ \prod_{i = 1} ^ {\min(n, m)} F(i) ^ {\left\lfloor\frac{n}{i}\right\rfloor \cdot \left\lfloor\frac{m}{i}\right\rfloor} $$

可以直接数论分块计算。

时间复杂度：$\mathcal O(n \log n + T \sqrt n \log n)$。

Code

#include <cstdio>
#include <algorithm>

const int N = 1e6 + 5;
const int P = 1e9 + 7, phiP = 1e9 + 6;

int T, n, m, tot, p[N], mu[N], f[N], iF[N], g[N], pre[N];
bool flg[N];

void add(int &x, int y) {
    (x += y) >= P && (x -= P);
}
int pow(int x, int p) {
    if (p < 0) p += phiP;
    int ans = 1;
    for (; p; p >>= 1, x = 1LL * x * x % P) {
        if (p & 1) ans = 1LL * ans * x % P;
    }
    return ans;
}
int inv(int x) {
    return pow(x, P - 2);
}
void sieve(int n) {
    std::fill(flg + 1, flg + n + 1, true);
    mu[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (flg[i]) {
            p[++tot] = i, mu[i] = -1;
        }
        for (int j = 1; j <= tot && i * p[j] <= n; j++) {
            flg[i * p[j]] = false;
            if (i % p[j] == 0) {
                mu[i * p[j]] = 0;
                break;
            } else {
                mu[i * p[j]] = -mu[i];
            }
        }
    }
    f[0] = 0, f[1] = 1;
    for (int i = 2; i <= n; i++) f[i] = (f[i - 1] + f[i - 2]) % P;
    for (int i = 1; i <= n; i++) iF[i] = inv(f[i]);
    std::fill(g + 1, g + n + 1, 1);
    for (int i = 1; i <= n; i++) {
        for (int j = i; j <= n; j += i) {
            int m = mu[j / i];
            if (!m) continue;
            g[j] = 1LL * g[j] * (m > 0 ? f[i] : iF[i]) % P;
        }
    }
    pre[0] = 1;
    for (int i = 1; i <= n; i++) {
        pre[i] = 1LL * pre[i - 1] * g[i] % P;
    }
}
int main() {
    sieve(N - 5);
    scanf("%d", &T);
    for (int _ = 1; _ <= T; _++) {
        scanf("%d%d", &n, &m);
        int ans = 1;
        for (int l = 1, r; l <= std::min(n, m); l = r + 1) {
            r = std::min(n / (n / l), m / (m / l));
            int pro = 1LL * pre[r] * inv(pre[l - 1]) % P;
            ans = 1LL * ans * pow(pro, 1LL * (n / l) * (m / l) % phiP) % P;
        }
        printf("%d\n", ans);
    }
    return 0;
}