前置
幂转和
(prodlimits_{i} g^{d_i}=g^{sumlimits_{i}d_i})，(g) 是常数，或者至少在这次 (prod) 中不会变化
gcd 的套路卷积
[egin{aligned}
ecause[x==1]&=sumlimits_{d|x}mu(d)\
herefore[gcd(x, y)==1]&=sumlimits_{d|gcd(x, y)}mu(d)
\&=sumlimits_{d}mu(d)[d|x][d|y]
end{aligned}
]
此时常常可以将 (sumlimits_{d}mu(d)) 提到式子最前面，而 ([d|x][d|y]) 常常可以与枚举 (x, y) 的上界进行抵消，就像这样：
(sumlimits_{x=1}^{n}[x|d]mathbf{f}(x)=sumlimits_{x=1}^{lfloorfrac{n}{d} floor}mathbf{f}(xd))
而当 (mathbf{f}(x)=mathbf{id_0}(x)=mathbf{1}(x)=1) 时，原式 (=lfloorfrac{n}{d} floor)
推公式
[egin{aligned}
prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}f_{gcd(i,j)}
&=prodlimits_{i=1}^{n}prodlimits_{j=1}^{m}prodlimits_{g=1}^{min(i, j)}f_g^{[gcd(i, j)==g]}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}[gcd(i, j)==g]}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{i=1}^{lfloorfrac{n}{g}
floor}sumlimits_{j=1}^{lfloorfrac{m}{g}
floor}[gcd(i, j)==1]}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{i=1}^{lfloorfrac{n}{g}
floor}sumlimits_{j=1}^{lfloorfrac{m}{g}
floor}sumlimits_{d|gcd(i, j)}mu(d)}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{i=1}^{lfloorfrac{n}{g}
floor}sumlimits_{j=1}^{lfloorfrac{m}{g}
floor}sumlimits_{d}[d|i][d|j]mu(d)}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{d}mu(d)sumlimits_{i=1}^{lfloorfrac{n}{g}
floor}[d|i]sumlimits_{j=1}^{lfloorfrac{m}{g}
floor}[d|j]}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{d}mu(d)lfloorfrac{n}{gd}
floorlfloorfrac{m}{gd}
floor}
\&=prodlimits_{g=1}^{min(n, m)}f_g^{sumlimits_{g|T}mu(frac{T}{g})lfloorfrac{n}{T}
floorlfloorfrac{m}{T}
floor}&( exttt{代换 }T=gd)
\&=prodlimits_{T}prodlimits_{g|T}^{min(n, m)}f_g^{mu(frac{T}{g})lfloorfrac{n}{T}
floorlfloorfrac{m}{T}
floor}
end{aligned}
]
设 (mathbf{F}(T)=prodlimits_{g|T}f_g^{mu(frac{T}{g})})
( herefore exttt{ 原式}=prodlimits_{T=1}^{min(n, m)}mathbf{F}(T)^{lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor})
(mathbf{F}(T)) 进行 (Theta(nlog n)) 预处理前缀积
其余部分用数论分块 (Theta(sqrt{n})) 计算
代码
const int N = 1e6 + 7, P = 1e9 + 7;
inline int64 fpow(int64 a, int64 b) {
int64 ret = 1;
while (b) {
if (b & 1) ret = ret * a % P;
a = a * a % P;
b >>= 1;
}
return ret;
}
#define inv(x) fpow(x, P - 2)
namespace euler {
bool IsntPrime[N];
int prs[N], pcnt, mu[N];
unsigned a[N];
void init_prs(int n) {
mu[1] = 1;
for (int i = 2; i <= n; ++i) {
if (!IsntPrime[i]) {
prs[++pcnt] = i;
mu[i] = -1;
}
for (int j = 1; j <= pcnt && prs[j] * i <= n; ++j) {
IsntPrime[prs[j] * i] = true;
if (i % prs[j] == 0) break;
mu[prs[j] * i] = -mu[i];
}
}
}
}
namespace fib {
int f[N], invf[N];
void init(int n) {
invf[1] = f[1] = 1;
for (int i = 2; i <= n; ++i) {
f[i] = (f[i - 1] + f[i - 2]) % P;
invf[i] = inv(f[i]);
}
}
}
namespace F {
int64 f[N], sum[N];
void init(int n) {
for (int g = 0; g <= n; ++g) f[g] = 1;
sum[0] = 1;
for (int g = 1; g <= n; ++g) {
for (int d = 1; d * g <= n; ++d) // T = d * g
switch (euler::mu[d]) {
case 1: f[g * d] = f[g * d] * fib::f[g] % P; break;
case -1: f[g * d] = f[g * d] * fib::invf[g] % P; break;
}
sum[g] = sum[g - 1] * f[g] % P;
// printf("f[%d] = %lld
", g, f[g]);
}
}
int64 prod(int l, int r) {
return sum[r] * inv(sum[l - 1]) % P;
}
}
signed main() {
euler::init_prs(1e6);
fib::init(1e6);
F::init(1e6);
int T, n, m;
for (cin >> T; T--;) {
cin >> n >> m;
int64 ans = 1;
for (int l = 1, r, bound = min(n, m); l <= bound; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = ans * fpow(F::prod(l, r), int64(n / l) * (m / l) % (P - 1)) % P;
}
cout << ans << endl;
}
return 0;
}