( ext{Problem})
求
[prod_{i=1}^n prod_{j=1}^m f({gcd(i,j)})
]
其中
[f(n)=
egin{cases}
0 & n=0 \
1 & n=1 \
f(n-1)+f(n-2) & n > 1
end{cases}
]
(T) 组数据，(1 le T le 10^3,1 le n,m le 10^6)
( ext{Analysis})
[egin{aligned}
prod_{i=1}^n prod_{j=1}^m f_{gcd(i,j)}
&= prod_{d=1} f(d)^{sum_{d|i} sum_{d|j} [gcd(i,j)=d]}
end{aligned}
]
我们把指数抽出来处理
[egin{aligned}
sum_{d|i} sum_{d|j} [gcd(i,j)=d]
&= sum_{i=1} sum_{j=1} [gcd(i,j)=1] \
&= sum_{g=1} mu(g) lfloor frac{n}{dg}
floor lfloor frac{m}{dg}
floor
end{aligned}
]
令 (T=dg)，那么把外面的积式也换乘 (T)
那么就得到了
[prod_{T=1}^n prod_{d|T} f(d)^{mu(frac{T}{d})lfloor frac{n}{T}
floor lfloor frac{m}{T}
floor}
]
(prod_{d|T} f(d)^{mu(frac{T}{d})}) 这一部分直接 (O(n log n)) 预处理即可
然后数论分快解决
( ext{Code})
#include<cstdio>
#include<iostream>
#define re register
using namespace std;
typedef long long LL;
const int N = 1e6, P = 1e9 + 7;
int n, m, k, totp, pr[N], vis[N + 5], sum[N + 5], mu[N + 5], f[N + 5], g[N + 5];
LL F[N + 5];
inline int fpow(LL x, LL y)
{
LL res = 1;
for(; y; y >>= 1)
{
if (y & 1) res = res * x % P;
x = x * x % P;
}
return res;
}
inline void Euler()
{
vis[1] = mu[1] = f[1] = g[1] = F[0] = F[1] = 1;
for(re int i = 2; i <= N; i++)
{
f[i] = (f[i - 1] + f[i - 2]) % P, g[i] = fpow(f[i], P - 2), F[i] = 1;
if (!vis[i]) pr[++totp] = i, mu[i] = -1;
for(re int j = 1; j <= totp && i * pr[j] <= N; j++)
{
vis[i * pr[j]] = 1;
if (!(i % pr[j])) break;
mu[i * pr[j]] = -mu[i];
}
}
for(re int i = 1; i <= N; i++)
{
if (!mu[i]) continue;
for(re int j = i; j <= N; j += i)
F[j] = F[j] * (mu[i] == 1 ? f[j / i] : g[j / i]) % P;
}
for(re int i = 2; i <= N; i++) F[i] = F[i] * F[i - 1] % P;
}
int main()
{
freopen("product.in", "r", stdin);
freopen("product.out", "w", stdout);
Euler();
int T; scanf("%d", &T);
for(; T; --T)
{
scanf("%d%d", &n, &m);
LL ans = 1;
for(re int l = 1, r, v; l <= min(n, m); l = r + 1)
{
r = min(n / (n / l), m / (m / l));
v = F[r] * fpow(F[l - 1], P - 2) % P;
ans = ans * fpow(v, (LL)(n / l) * (m / l) % (P - 1)) % P;
}
printf("%lld
", ans);
}
}