【SDOI2017】数字表格

题面

题解

这道题目还有一种比较有意思的解法。

定义一种运算\((\mathbf f\oplus\mathbf g)(x) = \prod\limits_{d\mid x}\mathbf f(d)^{\mathbf g(\frac xd)}\)

研究一下这种运算的性质：

虽然这个运算没有交换律也没有结合律，但是它有一个比较奇特的性质：

设运算\(*\)是狄利克雷卷积，那么可以证明\((\mathbf f \oplus \mathbf g) \oplus \mathbf h = \mathbf f \oplus (\mathbf g * \mathbf h)\)。

于是就有一种基于\(\prod\)的莫比乌斯反演：

\(\mathbf f = \mathbf g \oplus \mathbf 1 \Rightarrow \mathbf g = \mathbf f \oplus \mu\)

也就是\(\mathbf f(x) = \prod_{d|x} \mathbf g(d) \Rightarrow \mathbf g(x) = \prod_{d|x} \mathbf f(d)^{\mu(\frac xd)}\)

那么这道题目就很好推了。

\[\begin{aligned} &\prod_{i=1}^n\prod_{j=1}^m f[\gcd(i, j)] \\ =&\prod_{i=1}^n\prod_{j=1}^m\prod_{d|i, d|j} \mathbf g(d) \quad (\mathbf g = \mathbf f \oplus \mu) \\ =&\prod_{d=1}^n \mathbf g(d)^{\sum_{d|i}\sum_{d|j}1} \\ =&\prod_{d=1}^n \mathbf g(d)^{\left\lfloor \frac nd\right\rfloor \left\lfloor \frac md\right\rfloor} \end{aligned} \]

我们发现\(\mathbf g\)可以\(\mathrm{O}(n\log n)\)算，于是就做完了。

代码

这个代码貌似很古老了QAQ

#include<bits/stdc++.h>
#define RG register
#define clear(x, y) memset(x, y, sizeof(x));
using namespace std;

inline int read()
{
	int data=0, w=1;
	char ch=getchar();
	while(ch!='-'&&(ch<'0'||ch>'9')) ch=getchar();
	if(ch=='-') w=-1, ch=getchar();
	while(ch>='0'&&ch<='9') data=(data<<3)+(data<<1)+(ch^48), ch=getchar();
	return data*w;
}

const int mod(1e9+7), maxn(1e6+10), lim(1e6);
inline int fastpow(int x, int y)
{
	int ans=1;
	while(y)
	{
		if(y&1) ans=1ll*ans*x%mod;
		x=1ll*x*x%mod; y>>=1;
	}
	return ans;
}

int f[maxn], prime[maxn], g[maxn], cnt, inv, sum[maxn], mu[maxn], n, m, T;
bool not_prime[maxn];

inline void init()
{
	not_prime[1]=f[1]=g[1]=sum[0]=sum[1]=mu[1]=1;
	for(RG int i=2;i<=lim;i++)
	{
		f[i]=(f[i-1]+f[i-2])%mod;
		g[i]=fastpow(f[i], mod-2);
		sum[i]=1;
		if(!not_prime[i]) prime[++cnt]=i, mu[i]=-1;
		for(RG int j=1;j<=cnt && i*prime[j]<=lim;j++)
		{
			not_prime[i*prime[j]]=true;
			if(i%prime[j]) mu[i*prime[j]]=-mu[i];
			else break;
		}
	}
	for(RG int i=1;i<=lim;i++)
	{
		if(!mu[i]) continue;
		for(RG int j=i;j<=lim;j+=i)
			sum[j]=1ll*sum[j]*((~mu[i])?f[j/i]:g[j/i])%mod;
	}
	for(RG int i=1;i<=lim;i++) sum[i]=1ll*sum[i]*sum[i-1]%mod;
}

int main()
{
	init();
	T=read();
	while(T--)
	{
		n=read(); m=read();
		if(n>m) swap(n, m);
		RG int i=1, j, k, l, tmp, ans=1;
		while(i<=n)
		{
			k=n/i; l=m/i;
			j=min(n/k, m/l);
			tmp=1ll*sum[j]*fastpow(sum[i-1], mod-2)%mod;
			ans=1ll*ans*fastpow(tmp, 1ll*k*l%(mod-1))%mod;
			i=j+1;
		}
		printf("%d\n", (ans+mod)%mod);
	}
	return 0;
}