Description

link

題目背景

\(Doris\) 剛剛學習了 \(fibonacci\) 數列。用 \(f_i\) 表示數列的第 \(i\) 項，那麼

\(f_0=0,f_1=1\)

\(f_n=f_{n-1}+f_{n-2},n\geq 2\)

題目描述

\(Doris\) 用老師的超級計算機生成了一個\(n\times m\) 的表格，

第 \(i\) 行第 \(j\) 列的格子中的數是 \(f_{\gcd(i,j)}\)，其中 \(\gcd(i,j)\) 表示 \(i,j\) 的最大公約數。

\(Doris\) 的表格中共有 \(n\times m\) 個數，她想知道這些數的乘積是多少。

答案對 \(10^9+7\) 取模。

Solution

然而還是不行

表示一下式子吧

\[\prod_{i=1}^n\prod _ {j=1}^m f_{gcd(i,j)} \]

把這題換一換，然後轉成求每個 \(fib\) 數列中的項會出現多少次

用快速冪解決就好了，經過處理最優可以做到 \(O(Tn\sqrt n)\)

式子寫出來不復雜

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

然而依舊爆炸……

接著來

提出來 \(T=id\)

\[\prod_{T=1}^{n}\prod _ {d|T} (f_d)^{\mu(\frac T d) \lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor} \]

有可以算的固定值了

\[\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor \]

裡面的算出來就完事了

調和級數助力過掉此題

Code

#include<bits/stdc++.h>
using namespace std;
#define int long long
namespace yspm{
	inline int read()
	{
		int res=0,f=1; char k;
		while(!isdigit(k=getchar())) if(k=='-') f=-1;
		while(isdigit(k)) res=res*10+k-'0',k=getchar();
		return res*f;
	}
	const int N=1e6+10,mod=1e9+7;
	int pri[N],mu[N],cnt,fl[N],res[N],fib[N],inv[N];
	inline int add(int x,int y){return x+y>=mod? x+y-mod:x+y;}
	inline int ksm(int x,int y)
	{
		int res=1; for(;y;y>>=1,x=x*x%mod) if(y&1) res=res*x%mod;
		return res;
	}
	inline void prework()
	{
		fib[1]=1; res[0]=res[1]=1;  for(int i=2;i<N;++i) fib[i]=add(fib[i-1],fib[i-2]),inv[i]=ksm(fib[i],mod-2),res[i]=1;
		mu[1]=1; inv[1]=1;
		for(int i=2;i<N;++i)
		{
			if(!fl[i]) pri[++cnt]=i,mu[i]=-1;
			for(int j=1;j<=cnt&i*pri[j]<N;++j)
			{
				fl[i*pri[j]]=1;
				if(i%pri[j]==0) break;
				else mu[i*pri[j]]=-mu[i];
			}
		}
		for(int i=1;i<N;++i)
		{
			if(!mu[i]) continue;
			for(int j=i;j<N;j+=i)
			{
				res[j]=res[j]*(mu[i]==1?fib[j/i]:inv[j/i])%mod;
			}  
		}
		for(int i=2;i<N;++i) res[i]=res[i]*res[i-1]%mod;
		return ;
	}
	inline void work()
	{
		int inv=1,ans=1,n=read(),m=read(),l,r; if(n>m) swap(n,m);
		for(l=1;l<=n;l=r+1)
		{
			r=min(n/(n/l),m/(m/l));
			ans=ans*ksm(res[r]*ksm(res[l-1],mod-2)%mod,(n/l)%(mod-1)*(m/l)%(mod-1))%mod; 
		}
		printf("%lld\n",ans);
		return ;
	}
	signed main()
	{
		prework();
		int T=read(); while(T--) work();
		return 0;
	}
}
signed main(){return yspm::main();}
```**