P3704 [SDOI2017]数字表格

神，不会低头！

---

首先，有 出现，说明这道题很阴间；其次，有 出现，说明这道题不可做；最后，多测和取模，说明这道题是给神仙做的。

结束。

---

开玩笑哒，贺肯定还是要贺的。

题目简介

题目来源： 所属不详（可能是集训题之类的）

评测链接 ：https://www.luogu.com.cn/problem/P3704
评测链接 ：https://loj.ac/p/2000

形式化题目：给定一个 的矩阵，权值 ，求 ，多测，模数 。

数据范围：

这显然看起来是一道莫比乌斯反演的题目，我们可以考虑从结果入手：

$$ans=\prod_{i=1}^{n}\prod_{j=1}^{m}fib\Big({\gcd(i,j)} \Big) \nonumber$$

没有任何一个部分是我们能够快速处理的。很难受。

由于从古至今斐波那契都是人类智慧的一大难题，所以暂且不考虑，我们先解决 的问题，并想办法来让 变得像 一样好处理。我们考虑把 提出，通过枚举 的值来解决问题。

这里定义 ，化简答案：

$$\begin{align} ans&=\prod_{d=1}^{N}\Big(fib(d)\Big)^{f(d)} \nonumber \end{align}$$

容易发现这个 我们在之前也经常看见，所以我们可以先把 进行转化，让这个函数能够在 的时间内求出，想看推导过程的读者可以参照 ，这里给出结果：

$$f(d)=\sum_{w=1}^{N}\mu(w)\left\lfloor\frac{n}{wd}\right\rfloor\left\lfloor\frac{m}{wd}\right\rfloor \nonumber$$

然后我们把 带回式子中，得到一个较冗长但是时间复杂度更优的式子：

$$ans=\prod_{d=1}^{N}\Big(fib(d)\Big)^{\Large \Big(\sum\limits_{w=1}^{N}\mu(w)\left\lfloor\frac{n}{wd}\right\rfloor\left\lfloor\frac{m}{wd}\right\rfloor\Big) } \nonumber$$

然后容易发现并不是很好从某一部分入手了，但此时的时间复杂度对于 多测数据高达 的题目而言还是太困难了（虽然已经有 了）。

由于有极其难搞的阴间多测，所以我们不用考虑在 上优化，所以考虑把 提出考虑。因为乘方的性质，我们可以把幂上的那个 给剔下来，变换一下顺序，设 ，有：

$$\begin{align} ans&=\prod_{d=1}^{N}\Big(fib(d)\Big)^{\Large \Big(\sum\limits_{w=1}^{N}\mu(w)\left\lfloor\frac{n}{wd}\right\rfloor\left\lfloor\frac{m}{wd}\right\rfloor\Big) } \nonumber \\ &=\prod_{p=1}^{N}\Big(\prod_{d\mid p}fib(d) \Big)^{\mu(p)\left\lfloor\frac{n}{p}\right\rfloor\left\lfloor\frac{m}{p}\right\rfloor } \nonumber \end{align}$$

然后我们最后来解决 的问题，重设 表示为：

$$h(p)=\prod_{d\mid p}\Big(fib(d)\Big)^{\mu\left(\frac{p}{d}\right) } \nonumber$$

然后通过 来再次表示 ：

$$ans=\prod_{p=1}^{N}h(p)^{\left\lfloor\frac{n}{p}\right\rfloor\left\lfloor\frac{m}{p}\right\rfloor } \nonumber$$

容易发现指数部分可以通过数论分块在根号内结束，那如果我们可以预处理出 的前缀积，询问复杂度就在 ，其中 为数论分块， 为快速幂。

接下来的操作就在于，如何能在有限的时间复杂度内预处理出 的前缀积， 讲了枚举因数和枚举倍数两种方法，这里写出枚举倍数的方法：

$$\begin{align} h(p)&=\prod_{d\mid p}\Big(fib(d)\Big)^{\mu\left(\frac{p}{d}\right) } \nonumber \\ &=\prod_{i=1}^{N}fib\left(\frac{p}{i} \right)^{\mu(i)} \nonumber \end{align}$$

我们用 来替代 得到：

$$h(i\times j)=\prod_{i=1}^{N}\prod_{j=1}^{\lfloor\frac{N}{i}\rfloor} fib(j)^{\mu(i)} \nonumber$$

这样的话，我们可以在 的时间内预处理出所需要的 。

最终我们得到了一个时间复杂度为 的解法，可以通过本题。

注意 std::memset 的用法，不要向我一样被卡了十多分钟才看出来。另外，根据欧拉定理（），我们的整个运算都是建立在模数为 的前提下的，所以涉及到指数的运算都要将模数变为 ，这是必要的。

AC Code

// ----- Eternally question-----
// Problem: P3704 [SDOI2017]数字表格
// Contest: Luogu
// URL: https://www.luogu.com.cn/problem/P3704
// Memory Limit: 125 MB
// Time Limit: 2000 ms
// Written by: Eternity
// Time: 2023-01-03 20:42:24
// ----- Endless solution-------

#include<bits/stdc++.h>
#define re register
typedef long long ll;
template<class T>
inline void read(T &x)
{
	x=0;
	char ch=getchar(),t=0;
	while(ch<'0'||ch>'9') t|=ch=='-',ch=getchar();
	while(ch>='0'&&ch<='9') x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
	if(t) x=-x;
}
template<class T,class ...T1>
inline void read(T &x,T1 &...x1){ read(x),read(x1...); }
template<class T>
inline void write(T x)
{
	if(x<0) putchar('-'),x=-x;
	if(x>9) write(x/10);
	putchar(x%10+'0');
}
template<>
inline void write(bool x){ putchar(x?'1':'0'); }
template<>
inline void write(char c){ putchar(c); }
template<>
inline void write(char *s){ while(*s!='\0') putchar(*s++); }
template<>
inline void write(const char *s){ while(*s!='\0') putchar(*s++); }
template<class T,class ...T1>
inline void write(T x,T1 ...x1){ write(x),write(x1...); }
template<class T>
inline bool checkMax(T &x,T y){ return x<y?x=y,1:0; }
template<class T>
inline bool checkMin(T &x,T y){ return x>y?x=y,1:0; }
const int MAXN=1e6+10;
const ll Mod=1e9+7;
const int S=1e6;
int Test,N,M;
ll H[MAXN],Fib[MAXN];
ll Mu[MAXN],Pri[MAXN],Tot;
ll Inv_fib[MAXN],Mul_H[MAXN];
bool Is[MAXN];
inline ll qPow(ll a,ll b)
{
	// if(b<0) b+=Mod;
	ll res=1;
	while(b)
	{
		if(b&1) res=res*a%Mod;
		a=a*a%Mod;b>>=1;
	}
	return res;
}
inline void sieve(int n)
{
	Mu[1]=1,Is[1]=1;
	for(int i=2;i<=n;++i)
	{
		if(!Is[i]) Pri[++Tot]=i,Mu[i]=-1;
		for(int j=1;j<=Tot&&i*Pri[j]<=n;++j)
		{
			Is[i*Pri[j]]=1;
			if(i%Pri[j]==0)
			{
				Mu[i*Pri[j]]=0;
				break;
			}
			Mu[i*Pri[j]]=-Mu[i];
		}
	}
}
inline void init()
{
	sieve(S);
	// printf("Test:%d\n",qPow(2,5));
	Fib[0]=0,Fib[1]=Inv_fib[1]=1;
	for(int i=2;i<=S;++i) Fib[i]=(Fib[i-1]+Fib[i-2])%Mod;
	// for(int i=1;i<=10;++i) write(Fib[i],' ');puts("");
	for(int i=2;i<=S;++i) Inv_fib[i]=qPow(Fib[i],Mod-2);
	for(int i=1;i<=S;++i) H[i]=1ll;
	for(int i=1;i<=S;++i)
	{
		if(!Mu[i]) continue;
		for(int j=1;i*j<=S;++j)
			H[i*j]=H[i*j]*(Mu[i]==1?Fib[j]:Inv_fib[j])%Mod;
	}
	Mul_H[0]=1;
	for(int i=1;i<=S;++i) Mul_H[i]=Mul_H[i-1]*H[i]%Mod;
}
inline int getEd(int x,int n){ return n/(n/x); }
inline ll solve(int n,int m)
{
	ll res=1;
	int ed=std::min(n,m);
	for(int l=1,r=0;l<=ed;l=r+1)
	{
		r=std::min(std::min(getEd(l,n),getEd(l,m)),ed);
		ll btm=Mul_H[r]*qPow(Mul_H[l-1],Mod-2)%Mod;
		res=res*qPow(btm,1ll*(n/l)*(m/l)%(Mod-1))%Mod;
	}
	return res;
}
int main()
{
	// freopen(".in","r",stdin);
	// freopen(".out","w",stdout);
	read(Test);
	init();
	while(Test--)
	{
		read(N,M);
		write(solve(N,M),'\n');
	}
	return 0;
}
/*

*/