[SDOI2017]数字表格

Decribe:

求 \(\prod_{i=1}^{n}\prod_{j=1}^{m}f_{\gcd(i,j)}\)，其中 \(f_i\) 代表斐波那契数列的第 \(i\) 项。

Solution：

显然莫反启动！

\[\prod_{i=1}^{\min(n,m)}f_i^{\sum_{j=1}^{n}\sum_{k=1}^{m}[gcd(j,k)==i]} \]

\[\prod_{i=1}^{\min(n,m)}f_i^{\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{i}\rfloor}[\gcd(i,k)==1]} \]

\[\prod_{i=1}^{\min(n,m)}f_i^{\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{i}\rfloor}\epsilon(\gcd(j,k))} \]

\[\prod_{i=1}^{\min(n,m)}f_i^{\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{i}\rfloor}\sum_{d|\gcd(j,k)}\mu(d)} \]

\[\prod_{i=1}^{\min(n,m)}f_i^{\sum_{d=1}^{\min(\lfloor\frac{n}{i}\rfloor,\lfloor\frac{m}{i}\rfloor)}\mu(d)\sum_{j=1}^{\lfloor\frac{n}{id}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{id}\rfloor}1} \]

\[\prod_{i=1}^{\min(n,m)}f_i^{\sum_{d=1}^{\min(\lfloor\frac{n}{i}\rfloor,\lfloor\frac{m}{i}\rfloor)}\mu(d)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor} \]

这里可以两次分块了，但最后 40 过不了，所以请教了一下大佬 @ScatteredHope。以下是他的推导：

令 \(T\) 等于 \(id\)。

\[\prod_{T=1}^{\min(n,m)}\prod_{d|T}f_{d}^{\mu(\frac{T}{d})\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor} \]

然后这一块：

\[\prod_{d|T}f_d^{\mu(\frac{T}{d})} \]

可以 \(O(n \log n)\) 预处理，然后就可以 \(O(\sqrt {n}\log n)\) 单次回答，总时间复杂度 \(O(T \sqrt n \log n)\)。

Code:

bool _Start;
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
namespace IO
{
	#define TP template<typename T>
	#define TP_ template<typename T,typename ... T_>
	#ifdef DEBUG
	#define gc() (getchar())
	#else
	char buf[1<<20],*p1,*p2;
	#define gc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<20,stdin),p1==p2)?EOF:*p1++)
	#endif
	#ifdef DEBUG
	void pc(const char &c)
	{
		putchar(c);
	}
	#else
	char pbuf[1<<20],*pp=pbuf;
	void pc(const char &c)
	{
		if(pp-pbuf==1<<20)
			fwrite(pbuf,1,1<<20,stdout),pp=pbuf;
		*pp++=c;
	}
	struct IO{~IO(){fwrite(pbuf,1,pp-pbuf,stdout);}}_;
	#endif
	TP void read(T &x)
	{
		x=0;static int f;f=0;static char ch;ch=gc();
		for(;ch<'0'||ch>'9';ch=gc())ch=='-'&&(f=1);
		for(;ch>='0'&&ch<='9';ch=gc())x=(x<<1)+(x<<3)+(ch^48);
		f&&(x=-x);
	}
	TP void write(T x)
	{
		if(x<0)
			pc('-'),x=-x;
		static T sta[35],top;top=0;
		do
			sta[++top]=x%10,x/=10;
		while(x);
		while(top)
			pc(sta[top--]^48);
	}
	TP_ void read(T &x,T_&...y){read(x);read(y...);}
	TP void writeln(const T x){write(x);pc('\n');}
	TP void writesp(const T x){write(x);pc(' ');}
	TP_ void writeln(const T x,const T_ ...y){writesp(x);writeln(y...);}
	TP inline T max(const T &a,const T &b){return a>b?a:b;}
	TP_ inline T max(const T &a,const T_&...b){return max(a,max(b...));}
	TP inline T min(const T &a,const T &b){return a<b?a:b;}
	TP_ inline T min(const T &a,const T_&...b){return min(a,min(b...));}
	TP inline void swap(T &a,T &b){static T t;t=a;a=b;b=t;}
	TP inline T abs(const T &a){return a>0?a:-a;}
	#undef TP
	#undef TP_
}
using namespace IO;
using std::cerr;
using LL=long long;
constexpr int N=1e6+5;
constexpr LL mod=1e9+7;
int prime[N],pr;
LL mu[N],f[N],s[N],invs[N];
bool v[N];
LL qpow(LL a,LL b)
{
	LL ans=1ll;a%=mod;
	for(;b;b>>=1,a=a*a%mod)
		if(b&1)
			ans=ans*a%mod;
	return ans;
}
void init()
{
	mu[1]=1;f[1]=1;
	for(int i=2;i<N;i++)
	{
		f[i]=(f[i-1]+f[i-2])%mod;
		if(!v[i])
			prime[++pr]=i,mu[i]=-1;
		for(int j=1;j<=pr&&prime[j]*i<N;j++)
		{
			int k=prime[j]*i;
			v[k]=1;
			if(!(i%prime[j]))
				break;
			mu[k]=-mu[i];
		}
	}
	for(int i=0;i<N;i++)
		s[i]=1;
	for(int i=1;i<N;i++)
		for(int j=1;j*i<N;j++)
			s[j*i]=s[j*i]*(!mu[j]?1:(mu[j]==-1?qpow(f[i],mod-2):f[i]))%mod;
	invs[0]=1;
	for(int i=1;i<N;i++)
		s[i]=(s[i]*s[i-1])%mod,invs[i]=qpow(s[i],mod-2);
}
LL n,m;
bool _End;
int main()
{
//	fprintf(stderr,"%.2 MBlf\n",(&_End-&_Start)/1048576.0);
	int T;read(T);
	init();
	while(T--)
	{
		read(n,m);
		if(n>m)
			swap(n,m);
		LL ans=1ll;
		for(LL l=1,r;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans=ans*qpow(s[r]*invs[l-1],(n/l)*(m/l))%mod;
		}
		writeln(ans);
	}
	return 0;
}