BZOJ 4816. [Sdoi2017]数字表格
$$\prod_{i=1}^{n}\prod_{j=1}^{m}f(\gcd(i,j))$$
$$=\prod_{d=1}^{\min\{n,m\}}f(d)^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}[(i,j)=1]}$$
$$=\prod_{d=1}^{\min\{n,m\}}f(d)^{\sum_{d'}\mu(d')\lfloor\frac{n}{dd'}\rfloor\lfloor\frac{m}{dd'}\rfloor}$$
$$=\prod_{T}(\prod_{d|T}f(d)^{\mu(\frac{T}{d})})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$$
求出 $g(n) =\prod_{d|n}f(d)^{\mu(\frac{n}{d})}$ 的前缀积即可。
#include <bits/stdc++.h> const int MOD = 1e9 + 7; int qp(int a, int b = MOD - 2) { int ans = 1; while (b) { if (b & 1) ans = 1LL * ans * a % MOD; b >>= 1; a = 1LL * a * a % MOD; } return ans; } const int N = 1e6; int prime[N + 7], prin, mu[N + 7], g[N + 7], f[N + 7], fi[N + 7]; bool vis[N + 7]; inline int add(int x, int y) { return x + y >= MOD ? x + y - MOD : x + y; } void init() { mu[1] = 1; f[0] = 0; f[1] = 1; fi[1] = 1; g[0] = g[1] = 1; for (int i = 2; i <= N; i++) { g[i] = 1; f[i] = add(f[i - 1], f[i - 2]); fi[i] = qp(f[i]); if (!vis[i]) { prime[++prin] = i; mu[i] = -1; } for (int j = 1; j <= prin && i * prime[j] <= N; j++) { vis[i * prime[j]] = 1; if (i % prime[j] == 0) break; mu[i * prime[j]] = -mu[i]; } } for (int i = 1; i <= N; i++) if (mu[i]) for (int j = 1; 1LL * i * j <= N; j++) g[i * j] = 1LL * g[i * j] * (mu[i] == 1 ? f[j] : fi[j]) % MOD; for (int i = 2; i <= N; i++) g[i] = 1LL * g[i] * g[i - 1] % MOD; } int main() { init(); int T; scanf("%d", &T); while (T--) { int n, m; scanf("%d%d", &n, &m); if (n > m) std::swap(n, m); int ans = 1; for (int i = 1, j; i <= n; i = j + 1) { j = std::min(n / (n / i), m / (m / i)); ans = 1LL * ans * qp(1LL * g[j] * qp(g[i - 1]) % MOD, 1LL * (n / i) * (m / i) % (MOD - 1)) % MOD; } printf("%d\n", ans); } return 0; }