【洛谷P3704】数字表格

题目

题目链接：https://www.luogu.com.cn/problem/P3704
Doris 用老师的超级计算机生成了一个 \(n\times m\) 的表格，
第 \(i\) 行第 \(j\) 列的格子中的数是 \(fib_{\gcd(i,j)}\)，其中 \(\gcd(i,j)\) 表示 \(i,j\) 的最大公约数。
Doris 的表格中共有 \(n\times m\) 个数，她想知道这些数的乘积是多少。
答案对 \(10^9+7\) 取模。
\(Q\leq 1000;n,m\leq 10^6\)。

思路

\[\prod_{d}\prod^{n}_{i=1}\prod^{m}_{j=1}[\gcd(i,j)=d]\text{fib}(d) \]

\[=\prod_{d}\text{fib}(d)^{\sum^{\lfloor\frac{n}{d}\rfloor}_{i=1} \sum^{\lfloor\frac{m}{d}\rfloor}_{j=1}[\gcd(i,j)=1]} \]

\[=\prod_{d}\text{fib}(d)^{\sum_{i}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor} \]

看到 \(id\) 果断令 \(T=id\)，原式

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

\[=\prod^{n}_{T=1}\left (\prod_{d|T}\text{fib}(d)^{\mu(\frac{T}{d})}\right )^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor} \]

括号外面整除分块，括号里面预处理即可。
时间复杂度 \(O(n(\log p+\log n)+Q\sqrt{n})\)。\(n,m\) 同阶。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

const int N=1000010,MOD=1e9+7;
int Q,n,m,fib[N],inv[N],prm[N],mu[N];
ll ans,g[N],h[N];
bool v[N];

void findprm(int n)
{
	mu[1]=1;
	for (int i=2;i<=n;i++)
	{
		if (!v[i]) prm[++m]=i,mu[i]=-1;
		for (int j=1;j<=m;j++)
		{
			if (i>n/prm[j]) break;
			v[i*prm[j]]=1; mu[i*prm[j]]=-mu[i];
			if (!(i%prm[j])) { mu[i*prm[j]]=0; break; }
		}
	}
}

ll fpow(ll x,ll k)
{
	ll ans=1;
	for (;k;k>>=1,x=x*x%MOD)
		if (k&1) ans=ans*x%MOD;
	return ans;
}

int main()
{
	findprm(N-1);
	fib[1]=inv[1]=1;
	for (int i=0;i<N;i++) g[i]=1;
	for (int i=2;i<N;i++)
	{
		fib[i]=(fib[i-1]+fib[i-2])%MOD;
		inv[i]=fpow(fib[i],MOD-2);
		for (int j=i;j<N;j+=i)
			g[j]=g[j]*(mu[j/i]==1 ? fib[i] : (mu[j/i]==0 ? 1 : inv[i]))%MOD;
	}
	h[0]=1;
	for (int i=1;i<N;i++)
	{
		g[i]=g[i]*g[i-1]%MOD;
		h[i]=fpow(g[i],MOD-2);
	}
	scanf("%d",&Q);
	while (Q--)
	{
		scanf("%d%d",&n,&m);
		ans=1;
		for (int i=1,j;i<=min(n,m);i=j+1)
		{
			j=min(n/(n/i),m/(m/i));
			ans=ans*fpow(g[j]*h[i-1]%MOD,1LL*(n/i)*(m/i))%MOD;
		}
		printf("%lld\n",(ans%MOD+MOD)%MOD);
	}
	return 0;
}