[BZOJ4816][SDOI2017]数字表格（莫比乌斯反演）

容易看出，此题用莫比乌斯反演求解。
首先，要把$\prod_{i=1}^n\prod_{j=1}^mf[\gcd(i,j)]$换一个方向去思考，即不枚举$i,j$，而是枚举$d=\gcd(i,j)$。可以得到，在这个表格里，$f[d]$出现的次数就是$\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]$。
所以：
$Ans=\prod_df[d]^{\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=d]}$
$=\prod_df[d]^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[\gcd(i,j)=1]}$
而$\sum_{i=1}^x\sum_{j=1}^y[\gcd(i,j)=1]$是一个经典问题，它等于$\sum_d\lfloor\frac{x}{d}\rfloor\lfloor\frac{y}{d}\rfloor\mu(d)$。
即原式等于$\prod_df[d]^{\sum_k\lfloor\frac{n}{dk}\rfloor\lfloor\frac{m}{dk}\rfloor\mu(k)}$。
继续考虑。可以发现，$\lfloor\frac{n}{dk}\rfloor$和$\lfloor\frac{m}{dk}\rfloor$这两个式子不容易直接分块。所以这里令$u=dk$，原式化为：
$\prod_u\prod_{k|u}f[\frac{u}{k}]^{\mu(k)\lfloor\frac{n}{u}\rfloor\lfloor\frac{m}{u}\rfloor}$
$=\prod_u(\prod_{k|u}f[\frac{u}{k}]^{\mu(k)})^{\lfloor\frac{n}{u}\rfloor\lfloor\frac{m}{u}\rfloor}$。
这样就可以将$\lfloor\frac{n}{u}\rfloor$和$\lfloor\frac{m}{u}\rfloor$的值分块了。令$F[u]=\prod_{k|u}f[\frac{u}{k}]^{\mu(k)}$，那么就可以预处理 F <script type="math/tex" id="MathJax-Element-2861">F</script>的前缀积，计算区间的积时使用逆元计算。
代码：

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;
inline int read() {
    int res = 0; bool bo = 0; char c;
    while (((c = getchar()) < '0' || c > '9') && c != '-');
    if (c == '-') bo = 1; else res = c - 48;
    while ((c = getchar()) >= '0' && c <= '9')
        res = (res << 3) + (res << 1) + (c - 48);
    return bo ? ~res + 1 : res;
}
const int MaxN = 1e6, PYZ = 1e9 + 7, N = MaxN + 5;
int tot, f[N], g[N], pri[N], miu[N], F[N], prod[N];
bool mark[N];
int qpow(int a, int b) {
    int res = 1;
    while (b) {
        if (b & 1) res = 1ll * res * a % PYZ;
        a = 1ll * a * a % PYZ;
        b >>= 1;
    }
    return res;
}
void sieve() {
    int i, j; f[0] = 0; f[1] = miu[1] = g[1] = F[1] = prod[0] = 1;
    mark[0] = mark[1] = 1;
    for (i = 2; i <= MaxN; i++) {
        f[i] = (f[i - 1] + f[i - 2]) % PYZ;
        g[i] = qpow(f[i], PYZ - 2); F[i] = 1;
        if (!mark[i]) pri[++tot] = i, miu[i] = -1;
        for (j = 1; j <= tot; j++) {
            if (1ll * i * pri[j] > MaxN) break;
            mark[i * pri[j]] = 1;
            if (i % pri[j] == 0) break;
            else miu[i * pri[j]] = -miu[i];
        }
    }
    for (i = 1; i <= MaxN; i++) {
        if (miu[i] != 0) for (j = i; j <= MaxN; j += i)
            F[j] = 1ll * F[j] * (miu[i] == 1 ? f[j / i] : g[j / i]) % PYZ;
        prod[i] = 1ll * prod[i - 1] * F[i] % PYZ;
    }
}
int solve(int a, int b) {
    int i, n = min(a, b), ans = 1;
    for (i = 1; i <= n;) {
        int nxt = min(a / (a / i), b / (b / i)), pro;
        pro = 1ll * prod[nxt] * qpow(prod[i - 1], PYZ - 2) % PYZ;
        ans = 1ll * ans * qpow(pro, 1ll * (a / i) * (b / i) % (PYZ - 1)) % PYZ;
        i = nxt + 1;
    }
    return ans;
}
int main() {
    int a, b, T = read(); sieve();
    while (T--) a = read(), b = read(),
        printf("%d\n", solve(a, b));
    return 0;
}