bzoj-4816,P3704 [SDOI2017]数字表格

$Hyperlink$

https://www.lydsy.com/JudgeOnline/problem.php?id=4816
https://www.luogu.org/problem/P3704

---

$Description$

求

$$\prod _{i=1}^n\prod _{j=1}^m{f}_{gcd(i,j)}$$

$f$表示斐波那契数列

数据范围：$T\le 1000,1\le n,m\le 10^6$

---

$Solution$

问题转换成
$$\prod _{d=1}^n\prod _{i=1}^n\prod _{j=1}^m[gcd(i,j)==d]f_d$$
把$f_d$拉出来，每个$f_d$乘的次数就是${\sum }_{i=1}^{n/d}{\sum }_{j=1}^{m/d}[gcd(i,j)==1]$，所以

$$\prod _{d=1}^n{f}_d^{{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }{\sum }_{j=1}^{\lfloor \frac{m}{d}\rfloor }[gcd(i,j)==1]}$$

到这里数论分块套数论分块就可以做到$O(n)$了，但是会T飞

再根据之前莫比乌斯反演题的经验啊。。。
${\sum }_{i=1}^{n/d}{\sum }_{j=1}^{m/d}[gcd(i,j)==1]->{\sum }_{i=1}^{n/d}\mu (i)\lfloor \frac{n}{id}\rfloor \lfloor \frac{m}{id}\rfloor$

继续用YY的GCD的方法，设$T=id$，放入整个式子中
$$\prod _{d=1}^n{f}_d^{{\sum }_{i=1}^{\lfloor \frac{n}{d}\rfloor }\mu (i)\lfloor \frac{n}{id}\rfloor \lfloor \frac{m}{id}\rfloor }$$
$$->\prod _{T=1}^n(\prod _{d|T}{f}_d^{\mu (T/d)}{)}^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor }$$

然后暴力求中间那一坨，然后整除分块即可

时间复杂度：$O(A+B+C)$
$A:$线性筛+快速幂，$B:$暴力计算中间那一坨（跑不满），$C:$回答询问
$O(nlogn+nlo{g}^{2}n+T\sqrt{n}logn)$

---

$Code$

#include<cstdio>
#include<cctype>
#include<algorithm>
#define N 1000000
#define LL long long
#define mod 1000000007
using namespace std;
LL mu[N+10],f[N+10],F[N+10],g[N+10];int prime[N>>3],n,m,t;
bool vis[N+10];
inline LL read()
{
	LL f=0,d=1;char c; 
	while(c=getchar(),!isdigit(c)) if(c=='-') d=-1;f=(f<<3)+(f<<1)+c-48;
	while(c=getchar(),isdigit(c)) f=(f<<3)+(f<<1)+c-48;
	return d*f;
}
inline LL ksm(LL x,LL y)
{
	LL ans=1;
	for(;y;y>>=1,(x*=x)%=mod) if(y&1) (ans*=x)%=mod;
	return ans;
}
inline LL solve(int a,int b)//C
{
	if(a>b) a^=b,b=a^b,a^=b;
	LL res=1,inv=0;
	for(register int l=1,r=0;l<=a;l=r+1)
	{
		r=min(a/(a/l),b/(b/l));
		inv=1ll*F[r]*ksm(F[l-1],mod-2)%mod;
		(res*=ksm(inv,1ll*(a/l)*(b/l)%(mod-1)))%=mod;//扩展欧拉定理a^b=a^{b mod phi(p)}
	}
	return res;
}
signed main()
{
	mu[1]=1;f[1]=1;g[1]=1;F[0]=F[1]=1;
	for(register int i=2;i<=N;i++)//A
	{
		f[i]=(f[i-1]+f[i-2])%mod;g[i]=ksm(f[i],mod-2);F[i]=1;
		if(vis[i]==0) {prime[++m]=i;mu[i]=-1;}
		for(register int j=1;j<=m&&prime[j]*i<=N;j++)
		{
			vis[prime[j]*i]=true;
			if(i%prime[j]==0) break;
			mu[prime[j]*i]=-mu[i];
		}
	}
	for(register int i=1;i<=N;i++)//B
	{
		if(mu[i]==0) continue;
		for(register int j=i;j<=N;j+=i)
		 F[j]=1ll*F[j]*(mu[i]==1?f[j/i]:g[j/i])%mod;
	}
	for(register int i=2;i<=N;i++) F[i]=1ll*F[i]*F[i-1]%mod;
	for(t=read();t--;)
	{
		n=read();m=read();
		printf("%lld\n",solve(n,m));
	}
}