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

容易看出，此题用莫比乌斯反演求解。
首先，要把∏ni=1∏mj=1f[gcd(i,j)]换一个方向去思考，即不枚举i,j，而是枚举d=gcd(i,j)。可以得到，在这个表格里，f[d]出现的次数就是∑ni=1∑mj=1[gcd(i,j)=d]。
所以：
Ans=∏df[d]∑ni=1∑mj=1[gcd(i,j)=d]
=∏df[d]∑⌊nd⌋i=1∑⌊md⌋j=1[gcd(i,j)=1]
而∑xi=1∑yj=1[gcd(i,j)=1]是一个经典问题，它等于∑d⌊xd⌋⌊yd⌋μ(d)。
即原式等于∏df[d]∑k⌊ndk⌋⌊mdk⌋μ(k)。
继续考虑。可以发现，⌊ndk⌋和⌊mdk⌋这两个式子不容易直接分块。所以这里令u=dk，原式化为：
∏u∏k|uf[uk]μ(k)⌊nu⌋⌊mu⌋
=∏u(∏k|uf[uk]μ(k))⌊nu⌋⌊mu⌋。
这样就可以将⌊nu⌋和⌊mu⌋的值分块了。令F[u]=∏k|uf[uk]μ(k)，那么就可以预处理F的前缀积，计算区间的积时使用逆元计算。
代码：

#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;}