简要题意
令 \(f(i)\) 为斐波那契数列第 \(i\) 项的值。
\(T\) 组数据，对于每一个 \(n,m\)，求出：
\[\prod_{i=1}^{n}\prod_{j=1}^{m}f(\gcd(i,j))\pmod{10^9+7}
\]
\(1 \leq T \leq 10^3,1 \leq n,m\leq 10^6\)
思路
这里将介绍一种自认为比题解更为简便的方法
首先原式有 \(\prod\) 非常难搞，如果是 \(\sum\) 就好了。为了向毒瘤题 P5518 [MtOI2019]幽灵乐团 / 莫比乌斯反演基础练习题 致敬，这里我们使用 \(\ln+\exp\) 大法。
下面令 \(n\leq m\)。
将原式取 \(\ln\) 得：
\[\begin{aligned}
&\sum_{i=1}^{n}\sum_{j=1}^{m}\ln(f(\gcd(i,j)))\\
&=\sum_{d=1}^{n}\ln(f(\gcd(i,j)))\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=d]\\
&=\sum_{d=1}^{n}\ln(f(\gcd(i,j)))\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)\lfloor\frac{n}{di}\rfloor\lfloor\frac{m}{di}\rfloor\\
&=\sum_{d=1}^{n}\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor\sum_{k|d}(\ln(f(k))\cdot\mu(\frac{d}{k}))
\end{aligned}
\]
令 \(g(d)=\sum_{k|d}(\ln(f(k))\cdot\mu(\frac{d}{k}))\)，对其取 \(\exp\) 得 \(G(d)=\prod_{k|d}f(k)^{\mu(\frac{d}{k})}\)。显然这个可以 \(O(n\log n)\) 预处理。
然后对原式取 \(\exp\) 得：
\[\prod_{d=1}^{n}G(d)^{\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor}
\]
对 \(G\) 求预处理前缀积，可以使用数论分块做到 \(O(\sqrt{n}\log n)\)。
代码
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int N = 1e6+5;
int vis[N],mu[N],pri[N],tot,fib[N],f[N],inv[N];
const int mod = 1e9+7;
int pow(int a,int b,int mod) {
int ans=1;
for(; b; b>>=1,a=1ll*a*a%mod) {
if(b&1) {
ans=1ll*ans*a%mod;
}
}
return ans;
}
int M(const int x){
return (x%mod+mod)%mod;
}
void sieve(int ZYB = 1e6){
mu[1]=vis[1]=fib[1]=f[1]=f[0]=inv[1]=1;
for(int i=2;i<=ZYB;i++){
fib[i]=M(fib[i-1]+fib[i-2]);
f[i]=1;
inv[i]=pow(fib[i],mod-2,mod);
if(!vis[i]){
pri[++tot]=i;
mu[i]=-1;
}
for(int j=1;j<=tot&&i*pri[j]<=ZYB;j++){
vis[i*pri[j]]=1;
if(i%pri[j]) mu[i*pri[j]]=-mu[i];
else break;
}
}
for(int i=1;i<=ZYB;i++){
for(int j=i;j<=ZYB;j+=i){
f[j]=1ll*f[j]*(mu[i]==1?fib[j/i]:(mu[i]==0?1:inv[j/i]))%mod;
}
}
// for(int i=2;i<=ZYB;i++) cout<<f[i]<<' ';
for(int i=2;i<=ZYB;i++){
f[i]=M(f[i]*f[i-1]);
// cout<<f[i]<<' ';
}
}
signed main(){
sieve();
int t;cin>>t;
while(t--){
int n,m;
cin>>n>>m;
if(n>m) swap(n,m);
int ans=1;
for(int l=1,r=0;l<=n;l=r+1){
r=min(n/(n/l),m/(m/l));
int ret=M(f[r]*pow(f[l-1],mod-2,mod));
ans=M(ans*M(pow(ret,(n/l)*(m/l)%(mod-1),mod)));
}
cout<<M(ans)<<'\n';
}
}