【洛谷P3704】数字表格

题目

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

思路

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

[=prod_{d} ext{fib}(d)^{sum^{lfloorfrac{n}{d} floor}_{i=1} sum^{lfloorfrac{m}{d} floor}_{j=1}[gcd(i,j)=1]} ]

[=prod_{d} ext{fib}(d)^{sum_{i}mu(i)lfloorfrac{n}{id} floorlfloorfrac{m}{id} floor} ]

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

[=prod^{n}_{T=1}prod_{d|T} ext{fib}(d)^{mu(frac{T}{d})lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor} ]

[=prod^{n}_{T=1}left (prod_{d|T} ext{fib}(d)^{mu(frac{T}{d})} ight )^{lfloorfrac{n}{T} floorlfloorfrac{m}{T} floor} ]

括号外面整除分块，括号里面预处理即可。
时间复杂度 (O(n(log p+log n)+Qsqrt{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
",(ans%MOD+MOD)%MOD);
	}
	return 0;
}