bzoj4816: [Sdoi2017]数字表格

Description

Doris刚刚学习了\(fibonacci\)数列。用\(f[i]\)表示数列的第\(i\)项，那么

\[\begin{align}f[0]&=0\\f[1]&=1\\f[n]&=f[n-1]+f[n-2],n>=2\end{align} \]

Doris用老师的超级计算机生成了一个\(n\times m\)的表格，第\(i\)行第\(j\)列的格子中的数是\(f[\mathrm{gcd}(i,j)]\)，其中\(\mathrm{gcd}(i,j)\)表示\(i\),\(j\)的最大公约数。Doris的表格中共有\(n\times m\)个数，她想知道这些数的乘积是多少。答案对\(10^9+7\)取模。

Input

有多组测试数据。第一个一个数\(T\)，表示数据组数。

接下来\(T\)行，每行两个数\(n,m\)

\(T \leq 1000\),\(1 \leq n\),\(m \leq 10^6\)

Output

输出\(T\)行，第\(i\)行的数是第\(i\)组数据的结果

Sample Input

3
2 3
4 5
6 7

Sample Output

1
6
960

题解

要求的东西是

\[ans=\prod_{i=1}^n\prod_{j=1}^mf[\mathrm{gcd}(i,j)] \]

假设\(n<m\)，根据常见套路，枚举约数

\[ans=\prod_{d=1}^nf(d)^{\sum_{i=1}^n\sum_{j=1}^m[\mathrm{gcd}(i,j)==d]} \]

把指数拿出来化简

\[\begin{align} &\sum_{i=1}^n\sum_{j=1}^m[\mathrm{gcd}(i,j)==d]\\ =&\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}[\mathrm{gcd}(i,j)==1]\\ =&\sum_{i=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\sum_{j=1}^{\left \lfloor \frac{m}{d} \right \rfloor}\sum_{k|\mathrm{gcd}(i,j)}\mu(k)\\ =&\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{kd} \right \rfloor\left \lfloor \frac{m}{kd} \right \rfloor \end{align} \]

令\(T=kd\)，代入上式

\[\sum_{k=1}^{\left \lfloor \frac{n}{d} \right \rfloor}\mu(k)\left \lfloor \frac{n}{kd} \right \rfloor\left \lfloor \frac{m}{kd} \right \rfloor=\sum_{T=1}^n\mu(\frac{T}{d})\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor\\ \]

把\(T\)放到外面枚举

\[ans=\prod_{T=1}^n\prod_{d|T}f(d)^{\mu(\frac{T}{d})\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor} \]

令函数\(g(x)=\prod_{d|x}f(d)^{\mu(\frac{x}{d})}\)

\[ans=\prod_{T=1}^ng(T)^{\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor} \]

只要我们求出了\(g(x)\)的前缀积以及对应的逆元就可以整除分块了。

只要预处理出莫比乌斯函数和\(f\)以及\(f\)的逆元，\(g(x)\)就可以通过枚举每个数的筛它的倍数求出，复杂度是调和级数。

问题就是如何求所有\(f\)以及对应的逆元\(invf\)。

这里我们需要一个辅助的东西\(mulf\)，表示\(f\)的前缀积，即\(mulf(x)=\prod_{i=1}^x f(i)\).

假设我们知道了\(mulf(i+1)\)的逆元\(invm(i+1)\)，那么有

\[\begin{align} invm(i+1) &\equiv invm(i)\times invf(i+1) (\mathrm{mod}\ P)\\ invm(i)&\equiv invm(i+1)\times f(i+1)(\mathrm{mod}\ P) \end{align} \]

而\(invf(i)\equiv invm(i)*mulf(i-1)(\mathrm{mod} P)\),这样我们就可以求出\(f\)以及\(invf\)了。

用同样的方法，我们也可以求出\(g(x)\)的前缀积以及对应的逆元。

代码

#include<bits/stdc++.h>
#define MAXN 1000010
#define P 1000000007
#define LL long long
namespace IO{
	char buf[1<<15],*fs,*ft;
	inline char gc(){return (fs==ft&&(ft=(fs=buf)+fread(buf,1,1<<15,stdin),fs==ft))?0:*fs++;}
	inline int qr(){
		int x=0,rev=0,ch=gc();
		while(ch<'0'||ch>'9'){if(ch=='-')rev=1;ch=gc();}
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=gc();}
		return rev?-x:x;}
}using namespace IO;
using namespace std;
LL f[MAXN],invm[MAXN],invf[MAXN],g[MAXN],invg[MAXN],mulf[MAXN],mulg[MAXN],ans;
int p[MAXN],cnt,u[MAXN],T,N,M;
bool np[MAXN]; 
inline LL qp(LL x,LL y){
	LL ret=1;
	while(y){
		if(y&1)ret=ret*x%P;
		x=x*x%P;y>>=1; 
	}
	return ret;
}
int main(){
	#ifndef ONLINE_JUDGE
	freopen("product.in","r",stdin);
	freopen("product.out","w",stdout);
	#endif
	f[0]=0;
	f[1]=u[1]=np[1]=g[1]=mulf[1]=mulg[1]=invm[0]=1;
	for(int i=2;i<=1000000;i++){
		f[i]=(f[i-1]+f[i-2])%P;
		mulf[i]=mulf[i-1]*f[i]%P;
		g[i]=1;
		if(!np[i])p[++cnt]=i,u[i]=-1;
		for(int j=1;j<=cnt;j++){
			int t=p[j]*i;
			if(t>1000000)break;
			np[t]=1;
			if(i%p[j]==0){u[t]=0;break;}
			u[t]=-u[i];
		}
	}
	invm[1000000]=qp(mulf[1000000],P-2);
	invf[1000000]=qp(f[1000000],P-2);
	for(int i=999999;i>=1;i--){
		invm[i]=invm[i+1]*f[i+1]%P;
		invf[i]=mulf[i-1]*invm[i]%P;
	} 
	for(int i=2;i<=1000000;i++){
		for(int j=1,k;(k=j*i)<=1000000;j++){
			if(u[j]==1)g[k]=g[k]*f[i]%P;
			else if(u[j]==-1)g[k]=g[k]*invf[i]%P;
		}
		mulg[i]=mulg[i-1]*g[i]%P;
	}
	invg[1000000]=qp(mulg[1000000],P-2);
	for(int i=999999;i>=0;i--)invg[i]=invg[i+1]*g[i+1]%P;
	T=qr();
	while(T--){
		N=qr();M=qr();ans=1;
		if(N>M)swap(N,M);
		for(int i=1,j;i<=N;i=j+1){
			j=min(N/(N/i),M/(M/i));
			ans=ans*qp(mulg[j]*invg[i-1]%P,1LL*(N/i)*(M/i)%(P-1))%P;
		}
		printf("%lld\n",ans);
	}
	return 0;
}