[SDOI2017]数字表格

\(\text{Problem}\)

求

\[\prod_{i=1}^n \prod_{j=1}^m f({\gcd(i,j)}) \]

其中

\[f(n)= \begin{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\)

\(\text{Analysis}\)

\[\begin{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} \]

我们把指数抽出来处理

\[\begin{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} \rfloor \lfloor \frac{m}{dg} \rfloor \end{aligned} \]

令 \(T=dg\)，那么把外面的积式也换乘 \(T\)
那么就得到了

\[\prod_{T=1}^n \prod_{d|T} f(d)^{\mu(\frac{T}{d})\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor} \]

\(\prod_{d|T} f(d)^{\mu(\frac{T}{d})}\) 这一部分直接 \(O(n \log n)\) 预处理即可
然后数论分快解决

\(\text{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\n", ans);
	}
}