Luogu3704 SDOI2017数字表格

Description

link

题目背景

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

(f_0=0,f_1=1)

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

题目描述

(Doris) 用老师的超级计算机生成了一个(n imes m) 的表格，

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

(Doris) 的表格中共有 (n imes m) 个数，她想知道这些数的乘积是多少。

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

Solution

然而还是不行

表示一下式子吧

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

把这题换一换，然后转成求每个 (fib) 数列中的项会出现多少次

用快速幂解决就好了，经过处理最优可以做到 (O(Tnsqrt n))

式子写出来不复杂

[prod_{i=1}^{min(n,m)} f_i^{sumlimits_{d=1}^{lfloorfrac{n}{d} floor}mu(d)lfloorfrac{n}{id} floorlfloorfrac{m}{i d} floor}]

然而依旧爆炸……

接着来

提出来 (T=id)

[prod_{T=1}^{n}prod _ {d|T} (f_d)^{mu(frac T d) lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor} ]

有可以算的固定值了

[lfloorfrac{n}{d} floorlfloorfrac{m}{d} floor ]

里面的算出来就完事了

调和级数助力过掉此题

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
",ans);
		return ;
	}
	signed main()
	{
		prework();
		int T=read(); while(T--) work();
		return 0;
	}
}
signed main(){return yspm::main();}
```**