(f)为fbi数，即(f[0]=0,f[1]=1,f[i]=f[i-1]+f[i-2])，现在给出T组询问，询问(prod_{i=1}^nprod_{j=1}^mf(gcd(i,j)),T≤1000,1≤n,m≤10^6)。
解
通常我们所做的约数计数问题都是与累加有关，而此处与累乘有关，所以最自然的想法，要么是类比，要么必然转为累加，设(nleq m)，所以
[ans=prod_{i=1}^nprod_{j=1}^mf(gcd(i,j))=prod_{d=1}^nprod_{i=1}^nprod_{j=1}^mf(d)(gcd(i,j)==d)
]
(注：其实这样的表达不够严谨，因为当不满足条件，对应的值为0,所以含义不应该是简单的相乘，应该是(gcd(i,j)==d?f(d):1))
[=prod_{d=1}^nf(d)^{sum_{i=1}^nsum_{j=1}^m(gcd(i,j)==d)}
]
于是设
[f(d)=sum_{i=1}^nsum_{j=1}^m(gcd(i,j)==d)
]
[F(d)=[n/d][m/d]
]
由Mobius反演定理有
[f(d)=sum_{d|x}[n/x][m/x]mu(x/d)
]
代入即
[ans=prod_{d=1}^nf(d)^{sum_{d|x}[n/x][m/x]mu(x/d)}=
]
[prod_{d=1}^nprod_{d|x}f(d)^{[n/x][m/x]mu(x/d)}=
]
[prod_{d=1}^nprod_{d|x}f(d)^{[n/x][m/x]mu(x/d)}
]
[prod_{x=1}^n(prod_{d|x}f(d)^{mu(x/d)})^{[n/x][m/x]}
]
显然里面可以维护，采取倍数型维护方式即可，而对于次方自然可以整除分块，时间复杂度应为(O(nlog(n)+Tsqrt{n}log(n)))。
参考代码：
#include <iostream>
#include <cstdio>
#include <cstring>
#define il inline
#define ri register
#define ll long long
#define lim 1000000
#define yyb 1000000007
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[lim+1];
int mu[lim+1],prime[100000],pt;
ll f[lim+1],pf[lim+1],pfi[lim+1],ans;
il ll pow(ll,ll);
il void prepare(int);
template<class free>il free Min(free,free);
int main(){
int sxr,n,m,i,j;
scanf("%d",&sxr),prepare(lim);
while(sxr--){
scanf("%d%d",&n,&m);
if(n>m)swap(n,m);ans=1;
for(i=1;i<=n;i=j+1)
j=Min(n/(n/i),m/(m/i)),
(ans*=pow(pf[j]*pfi[i-1]%yyb,(ll)(n/i)*(m/i)))%=yyb;
printf("%lld
",ans);
}
return 0;
}
template<class free>
il free Min(free a,free b){
return a<b?a:b;
}
il ll pow(ll x,ll y){
ll ans(1);while(y){
if(y&1)(ans*=x)%=yyb;
(x*=x)%=yyb,y>>=1;
}return ans;
}
il void prepare(int n){
mu[1]|=check[1]|=true;
f[1]|=pf[0]|=pf[1]|=pfi[0]|=true;
int i,j;
for(i=2;i<=n;++i){
if(!check[i])prime[++pt]=i,mu[i]=-1;
for(j=1;j<=pt&&prime[j]*i<=n;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j]))break;
mu[i*prime[j]]=~mu[i]+1;
}f[i]=(f[i-1]+f[i-2])%yyb,pf[i]=pf[i-1]*f[i]%yyb;
}pfi[n]=pow(pf[n],yyb-2);
for(i=n;i>1;--i)
pfi[i-1]=pfi[i]*f[i]%yyb;
for(i=0;i<=n;++i)f[i]=1;
for(i=1;i<=n;++i)
for(j=i;j<=n;j+=i)
if(mu[j/i]==1)(f[j]*=pf[i]*pfi[i-1]%yyb)%=yyb;
else if(mu[j/i]==-1)(f[j]*=pfi[i]*pf[i-1]%yyb)%=yyb;
for(i=1;i<=n;++i)pf[i]=f[i]*pf[i-1]%yyb;pfi[n]=pow(pf[n],yyb-2);
for(i=n;i>1;--i)pfi[i-1]=pfi[i]*f[i]%yyb;
}