【SDOI2017】数字表格

题面

Description

Doris刚刚学习了(fibonacci)数列。用(f[i])表示数列的第(i)项，那么(f[0]=0,f[1]=1,f[n]=f[n-1]+f[n-2],ngeq2)。
Doris用老师的超级计算机生成了一个(ncdot m)的表格，第(i)行第(j)列的格子中的数是(f[gcd(i,j)])，其中(gcd(i,j))表示(i,j)的最大公约数。Doris的表格中共有(ncdot m)个数，她想知道这些数的乘积是多少。答案对(10^9+7)取模。

Input

有多组测试数据。

第一个一个数(T)，表示数据组数。

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

(（Tleq 1000,1leq n,mleq 10^6）)

Output

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

Sample Input

3

2

3

4

5

6

7

Sample Output

1

6

960

Hint

题目分析

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

[egin{split} ans&=prodlimits_{i=1}^nprodlimits_{j=1}^mF(gcd(i,j))\ &=prodlimits_{d=1}^nprodlimits_{i=1}^nprodlimits_{j=1}^mF(d)[gcd(i,j)==d]\ &=prodlimits_{d=1}^nF(d)^{sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)==d]}\ &=prodlimits_{d=1}^nF(d)^{sumlimits_{i=1}^{lfloorfrac nd floor}sumlimits_{j=1}^{lfloorfrac md floor}[gcd(i,j)==1]}\ end{split} ]

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

设(f(x)=sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)==x] Rightarrow g[x]=sumlimits_{i=1}^nsumlimits_{j=1}^m[x|gcd(i,j)]=lfloorfrac nx floorlfloorfrac mx floor)；

所以(f(x)=sumlimits_{x|d}^nmu(frac dx)lfloorfrac nd floorlfloorfrac md floorRightarrow ans=prodlimits_{d=1}^nF(d)^{sumlimits_{i=1}^{lfloorfrac nd floor}mu(i)lfloorfrac n{id} floorlfloorfrac m{id} floor})

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

我们设(T=id)

所以有

[egin{split} ans&=prodlimits_{d=1}^nF(d)^{sumlimits_{i=1}^{lfloorfrac nd floor}mu(i)lfloorfrac n{id} floorlfloorfrac m{id} floor}\ &=prodlimits_{d=1}^nF(d)^{sumlimits_{d|T}^{n}mu(frac Td)lfloorfrac n{T} floorlfloorfrac m{T} floor}\ &=prodlimits_{d=1}^nprodlimits_{d|T}^nF(d)^{mu(frac Td)lfloorfrac n{T} floorlfloorfrac m{T} floor}\ &=prodlimits_{T=1}^nprodlimits_{d|T}F(d)^{mu(frac Td)lfloorfrac n{T} floorlfloorfrac m{T} floor}\ &=prodlimits_{T=1}^n(prodlimits_{d|T}F(d)^{mu(frac Td)})^{lfloorfrac n{T} floorlfloorfrac m{T} floor}\ 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<<'
';
	}
	return 0;
}