【SDOI2017】数字表格

题面

Description

Doris刚刚学习了\(fibonacci\)数列。用\(f[i]\)表示数列的第\(i\)项，那么\(f[0]=0,f[1]=1,f[n]=f[n-1]+f[n-2],n\geq2\)。
Doris用老师的超级计算机生成了一个\(n\cdot m\)的表格，第\(i\)行第\(j\)列的格子中的数是\(f[gcd(i,j)]\)，其中\(gcd(i,j)\)表示\(i,j\)的最大公约数。Doris的表格中共有\(n\cdot 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

Hint

题目分析

设\(F(i)\)表示斐波那契数列的第\(i\)项，

\[\begin{split} ans&=\prod\limits_{i=1}^n\prod\limits_{j=1}^mF(gcd(i,j))\\ &=\prod\limits_{d=1}^n\prod\limits_{i=1}^n\prod\limits_{j=1}^mF(d)[gcd(i,j)==d]\\ &=\prod\limits_{d=1}^nF(d)^{\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)==d]}\\ &=\prod\limits_{d=1}^nF(d)^{\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}\sum\limits_{j=1}^{\lfloor\frac md\rfloor}[gcd(i,j)==1]}\\ \end{split} \]

\(F(d)\)的指数是一个比较模板的莫比乌斯反演。

设\(f(x)=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[gcd(i,j)==x] \Rightarrow g[x]=\sum\limits_{i=1}^n\sum\limits_{j=1}^m[x|gcd(i,j)]=\lfloor\frac nx\rfloor\lfloor\frac mx\rfloor\)；

所以\(f(x)=\sum\limits_{x|d}^n\mu(\frac dx)\lfloor\frac nd\rfloor\lfloor\frac md\rfloor\Rightarrow ans=\prod\limits_{d=1}^nF(d)^{\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}\mu(i)\lfloor\frac n{id}\rfloor\lfloor\frac m{id}\rfloor}\)

此时已经可以用整除分块\(O(n)\)求解，但是答案需要一个\(O(\sqrt n)\)的做法。

我们设\(T=id\)

所以有

\[\begin{split} ans&=\prod\limits_{d=1}^nF(d)^{\sum\limits_{i=1}^{\lfloor\frac nd\rfloor}\mu(i)\lfloor\frac n{id}\rfloor\lfloor\frac m{id}\rfloor}\\ &=\prod\limits_{d=1}^nF(d)^{\sum\limits_{d|T}^{n}\mu(\frac Td)\lfloor\frac n{T}\rfloor\lfloor\frac m{T}\rfloor}\\ &=\prod\limits_{d=1}^n\prod\limits_{d|T}^nF(d)^{\mu(\frac Td)\lfloor\frac n{T}\rfloor\lfloor\frac m{T}\rfloor}\\ &=\prod\limits_{T=1}^n\prod\limits_{d|T}F(d)^{\mu(\frac Td)\lfloor\frac n{T}\rfloor\lfloor\frac m{T}\rfloor}\\ &=\prod\limits_{T=1}^n(\prod\limits_{d|T}F(d)^{\mu(\frac Td)})^{\lfloor\frac n{T}\rfloor\lfloor\frac m{T}\rfloor}\\ \end{split} \]

此时，括号内的内容可以在调和级数的时间复杂度下预处理，而括号外的内容可以用整除分块\(O(\sqrt n)\)计算。

此时，时间复杂度可以满足要求。

代码实现

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdio>
#include<iomanip>
#include<cstdlib>
#define MAXN 0x7fffffff
#define int LL
typedef long long LL;
const int N=1e6+5,mod=1000000007;
using namespace std;
inline int Getint(){register int x=0,F=1;register char ch=getchar();while(!isdigit(ch)){if(ch=='-')F=-1;ch=getchar();}while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}return x*F;}
LL ksm(LL x,LL k){
	LL ret=1;x%=mod;
	while(k){
		if(k&1)ret=ret*x%mod;
		x=x*x%mod;
		k>>=1; 
	}
	return ret;
}
int mu[N],prime[N];
bool vis[N];
int F[N],G[N];
int g[N];
signed main(){
	mu[1]=F[1]=G[1]=g[1]=1;g[0]=1;
	for(int i=2;i<=1e6;i++){
		if(!vis[i])prime[++prime[0]]=i,mu[i]=-1;
		for(int j=1;j<=prime[0]&&i*prime[j]<=1e6;j++){
			vis[i*prime[j]]=1;
			if(i%prime[j]==0)break;
			mu[i*prime[j]]=-mu[i];
		}
		F[i]=(F[i-1]+F[i-2])%mod;
		G[i]=ksm(F[i],mod-2);
		g[i]=1;
	}
	
	for(int i=1;i<=1e6;i++){
		g[i]=1ll*(g[i]*g[i-1])%mod; 
		if(!mu[i])continue;
		for(int j=1;i*j<=1e6;j++)
			g[i*j]=g[i*j]*(mu[i]==-1?G[j]:F[j])%mod;
	} 
	
	int T=Getint();
	while(T--){
		int n=Getint(),m=Getint();
		if(n>m)swap(n,m);
		int ans=1;
		for(int l=1,r;l<=n;l=r+1){
			r=min(n/(n/l),m/(m/l));
			ans=1ll*ans*ksm(g[r]*ksm(g[l-1],mod-2)%mod,1ll*(n/l)*(m/l))%mod;
		} 
		cout<<(ans+mod)%mod<<'\n';
	}
	return 0;
}