[SDOI2017]数字表格
\(f\)为fbi数，即\(f[0]=0,f[1]=1,f[i]=f[i-1]+f[i-2]\)，现在给出T组询问，询问\(\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j)),T≤1000,1≤n,m≤10^6\)。
解
通常我们所做的约数计数问题都是与累加有关，而此处与累乘有关，所以最自然的想法，要么是类比，要么必然转为累加，设\(n\leq m\)，所以
\[ans=\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j))=\prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^mf(d)(gcd(i,j)==d)
\]
(注：其实这样的表达不够严谨，因为当不满足条件，对应的值为0,所以含义不应该是简单的相乘，应该是\(gcd(i,j)==d?f(d):1\))
\[=\prod_{d=1}^nf(d)^{\sum_{i=1}^n\sum_{j=1}^m(gcd(i,j)==d)}
\]
于是设
\[f(d)=\sum_{i=1}^n\sum_{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}^n\prod_{d|x}f(d)^{[n/x][m/x]\mu(x/d)}=
\]
\[\prod_{d=1}^n\prod_{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)+T\sqrt{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\n",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;
}
