題目連結：

題目

部落格園

題目大意：

快速求：

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

其中 \(f_i\) 表示斐波那契數列第 \(i\) 個數。

正文：

在寫本題之前，建議先拿 【Luogu P2257】 YY的GCD 練練手。

將原式化為我們能夠接受的時限：

\[\begin{aligned}\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)}&=\prod_{d=1}^{\min(n,m)}\prod_{i=1}^{n}\prod_{j=1}^{m}f_d\left[d==\gcd(i,j)\right]\\&=\prod_{d=1}^{\min(n,m)}{f_d}^{\left(\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\left[\gcd(i,j)==1\right]\right)}\\&=\prod_{d=1}^{\min(n,m)}{f_d}^{\left(\sum_{k|d}\mu(k)\left\lfloor\frac{n}{kd}\right\rfloor\left\lfloor\frac{m}{kd}\right\rfloor\right)}\\&=\prod_{T=1}^{\min(n,m)}\left(\prod_{d|T}{f_d}^{\mu(\frac{T}{d})}\right)^{\left\lfloor\frac{n}{T}\right\rfloor\left\lfloor\frac{m}{T}\right\rfloor}\end{aligned} \]

設 \(F_i\) 括號內的數，即 \(F_i=\prod_{d|i}{f_d}^{\mu(\frac{i}{d})}\)，這個函式用 \(O(n\log n)\) 預處理出來，在外面再套一個整除分塊就是答案。

程式碼：

inline ll qpow(ll a, ll b)
{
	if(b < 0) b = p - 2;
	ll ans = 1;
	for (; b; b >>= 1)
	{
		if (b & 1) ans = (ans * a % p + p) % p;
		a = (a * a % p + p) % p;
	}
	return ans;
}

inline void prework()
{
	miu[1] = 1;
	for (ll i = 2; i <= N - 10; i++)
	{
		if(!vis[i]) {pri[++cnt] = i, miu[i] = -1;}
		for (int j = 1; j <= cnt && pri[j] * i <= N - 10; j++)
		{
			vis[pri[j] * i] = 1;
			if (i % pri[j] == 0)
			{
				miu[i * pri[j]] = 0;
				break;
			}
			else
				miu[i * pri[j]] = -miu[i];
		}
	}
	sum[1] = fib[1] = fib_inv[1] = 1LL;			//求斐波那契數列及其逆元（這裡用fib代替）
	for (int i = 2; i <= N - 10; i++)
		fib[i] = (fib[i - 1] + fib[i - 2]) % p,
		fib_inv[i] = qpow(fib[i], p - 2), 
		sum[i] = 1LL;
	for (int i = 1; i <= N - 10; i++)			//求F（n）及其字首和（積）
		for (int j = i; j <= N - 10; j += i)
		{
			int t = j / i;
			if(miu[t] == 1) sum[j] = (sum[j] * fib[i] % p + p) % p;
			else if(miu[t] == -1) sum[j] = (sum[j] * fib_inv[i] % p + p) % p;
		}
	inv[0] = sum[0] = 1;
	for (int i = 1; i <= N - 10; i++)
		sum[i] = (sum[i] * sum[i - 1] % p + p) % p,
		inv[i] = qpow(sum[i], p - 2);
}

int main()
{
	prework();
	for (read(t); t--; )
	{
		ans = 1LL;
		read(n);read(m);
		if(n > m)
		{
			ll c = n; n = m; m = c;
		}
		for (register int l = 1, r; l <= n; l = r + 1)
		{
			r = min (n / (n / l), m / (m / l));
			ans = ans * qpow((sum[r] * inv[l - 1] % p + p) % p, (1ll * (n / r) * (m / r))) % p;
		}
		print(ans);
	}
	return 0;
}