题目
题目链接：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;
}