主要思路

先写出答案的式子：

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

老套路，枚举公约数：

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

后面那个部分拿下来进行反演：

\[ \begin{aligned} \sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = d] &= \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{d} \rfloor} [\gcd(i, j) = 1] \\ &= \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{d} \rfloor} \sum_{x|i, x|j} \mu(x) \\ &= \sum_{x = 1}^n \mu(x) \sum_{i = 1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{d} \rfloor} [x|i] [x|j] \\ &= \sum_{x = 1}^n \mu(x) \sum_{i = 1}^{\lfloor \frac{n}{dx} \rfloor} \sum_{j = 1}^{\lfloor \frac{m}{dx} \rfloor} 1 \\ &= \sum_{x = 1}^n \mu(x) {\lfloor \frac{n}{dx} \rfloor} {\lfloor \frac{m}{dx} \rfloor} \\ &= \sum_{T = 1}^n {\lfloor \frac{n}{T} \rfloor} {\lfloor \frac{m}{T} \rfloor} \mu(\frac{T}{d}) \end{aligned} \]

放到原式：

\[ \begin{aligned} \prod_{d = 1} d^{\sum_{T = 1}^d {\lfloor \frac{n}{T} \rfloor} {\lfloor \frac{m}{T} \rfloor} \mu(\frac{T}{d})} &= \prod_{d = 1} \prod_{T = 1}^d d^{{\lfloor \frac{n}{T} \rfloor} {\lfloor \frac{m}{T} \rfloor} \mu(\frac{T}{d})} \\ &= \prod_{T = 1}^n (\prod_{d|T} d^{\mu(\frac{T}{d})})^{{\lfloor \frac{n}{T} \rfloor} {\lfloor \frac{m}{T} \rfloor}} \end{aligned} \]

我们可以考虑把\(\prod_{d|T} d^{\mu(\frac{T}{d})}\)作为整体来求，然后再整除分块来搞（当然前面那个部分要从前缀和改成前缀积）。且我们发现\(d\)一定是 Fibonacci 数列的一项，所以这个部分我们可以用埃筛的方式来进行求解：把贡献一个个乘上去。

代码

// P3704.cpp
#include <bits/stdc++.h>

using namespace std;

const int MAX_N = 1e6 + 200, mod = 1e9 + 7;

int T, n, m, mu[MAX_N], primes[MAX_N], tot, f[MAX_N], g[MAX_N];
bool vis[MAX_N];

int quick_pow(int bas, int tim)
{
    int ret = 1;
    while (tim)
    {
        if (tim & 1)
            ret = 1LL * ret * bas % mod;
        bas = 1LL * bas * bas % mod;
        tim >>= 1;
    }
    return ret;
}

void init()
{
    mu[1] = 1;
    for (int i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            primes[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && 1LL * i * primes[j] < MAX_N; j++)
        {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0)
            {
                mu[i * primes[j]] = 0;
                break;
            }
            mu[i * primes[j]] = -mu[i];
        }
    }
    for (int i = 1; i < MAX_N; i++)
        f[i] = 1, g[i] = 1;
    int A = 1, B = 0;
    for (int i = 1; i < MAX_N; i++)
    {
        B = (A + B) % mod, A = (B - A + mod) % mod;
        int bucket[3] = {quick_pow(B, mod - 2), 1, B};
        for (int j = i, cnt = 1; j < MAX_N; j += i, cnt++)
        {
            f[j] = 1LL * f[j] * bucket[1 + mu[cnt]] % mod;
            g[j] = 1LL * g[j] * bucket[1 - mu[cnt]] % mod;
        }
    }
    f[0] = g[0] = 1;
    for (int i = 1; i < MAX_N; i++)
        f[i] = (1LL * f[i] * f[i - 1]) % mod, g[i] = (1LL * g[i] * g[i - 1]) % mod;
}

int main()
{
    scanf("%d", &T), init();
    while (T--)
    {
        scanf("%d%d", &n, &m);
        if (n > m)
            swap(n, m);
        int ans = 1;
        for (int x = 1, gx; x <= n; x = gx + 1)
        {
            gx = min(n / (n / x), m / (m / x));
            ans = 1LL * ans * quick_pow(1LL * f[gx] * g[x - 1] % mod, 1LL * (n / x) * (m / x) % (mod - 1)) % mod;
        }
        printf("%d\n", ans);
    }
    return 0;
}