P3704 [SDOI2017]数字表格

题目描述

洛谷

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

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

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

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

数据范围： \(T\leq 10^3, 1\leq n,m\leq 10^6\)

solution

莫比乌斯反演好题。

题目让我们求的是这个式子： \(\displaystyle\prod _{i=1}^{n}\prod_{j=1}^{m} f_{\gcd(i,j)}\) 。

我们可以考虑 \(f_{i}\) 被算了多少次，则有：

\(\large \displaystyle\prod_{d=1}^{n} f_{d} ^{\sum_{i=1}^{n}\sum_{j=1}^{m} [gcd(i,j) == d]}\)

我们尝试把指数反演一下，则有：

\(\ \ \displaystyle\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j) = d]\)

\(=\displaystyle\sum_{i=1}^{n\over d}\sum_{j=1}^{m\over d} [\gcd(i,j) = 1]\)

\(= \displaystyle\sum_{i=1}^{n\over d}\sum_{j=1}^{m\over d} \sum_{p\mid \gcd(i,j)} \mu(p)\)

\(= \displaystyle\sum_{p=1}^{n\over d} {n\over dp} {m\over dp} \mu(p)\)

把我们化简后的式子带入原式可得:

\(\large \displaystyle\prod_{d=1}^{n} f_d ^{\sum_{p=1}^{n\over d} {n\over dp} {m\over dp} \mu(p)}\)

设 \(Q = dp\), 则原式等于：

\(\large \displaystyle\prod_{Q=1} \left(\prod_{d\mid Q} f_d^{\mu({Q\over d})}\right)^{{n\over Q} {m\over Q}}\)

设 \(g(n) = \displaystyle\prod _{d\mid n} f_d^{\mu({n\over d})}\) （也就是上式中括号里面的式子）, 我们可以通过枚举倍数，预处理出 \(g(1)-g(n)\)。

根据调和级数可知，预处理的复杂度为 \(O(nlnn)\) 的。

维护一下 \(g(i)\) 的前缀积，我们就可以对 \({n\over Q}, {m\over Q}\) 整除分块了。

复杂度： \(O(T\sqrt{n} + nlnn)\) 。

Code

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
#define int long long
const int N = 1e6+10;
const int p = 1e9+7;
int T,n,m,tot;
int f[N],mu[N],sum[N],prime[N];
bool check[N]; 
inline int read()
{
    int s = 0, w = 1; char ch = getchar();
    while(ch < '0' || ch > '9'){if(ch == '-') w = -1; ch = getchar();}
    while(ch >= '0' && ch <= '9'){s = s * 10 + ch - '0'; ch = getchar();}
    return s * w;
}
int ksm(int a,int b)
{
	int res = 1;
	for(; b; b >>= 1)
	{
		if(b & 1) res = res * a % p;
		a = a * a % p;
	}
	return res; 
}
void YYCH()
{
	mu[1] = 1; f[0] = 0, f[1] = 1; sum[0] = 1;
	for(int i = 2; i <= N-5; i++) f[i] = (f[i-1] + f[i-2]) % p;
	for(int i = 2; i <= N-5; i++)
	{
		if(!check[i])
		{
			prime[++tot] = i;
			mu[i] = -1;
		}
		for(int j = 1; j <= tot && i * prime[j] <= N-5; j++)
		{
			check[i*prime[j]] = 1;
			if(i % prime[j] == 0)
			{
				mu[i * prime[j]] = 0;
				break;
			} 
			else mu[i * prime[j]] = -mu[i];
		}
	}
	for(int i = 1; i <= N-5; i++) if(mu[i] == -1) mu[i] = p-2;
	for(int i = 1; i <= N-5; i++) sum[i] = 1;
	for(int i = 1; i <= N-5; i++)//调和级数算g(i)
	{
		for(int j = 1; i * j <= N-5; j++)
		{
			sum[i * j] = sum[i * j] * ksm(f[i],mu[j]) % p;
		}
	}
	for(int i = 1; i <= N-5; i++) sum[i] = sum[i-1] * sum[i] % p;//sum数组为前缀积
}
signed main()
{
	T = read(); YYCH();
	while(T--)
	{
		n = read(); m = read();
		if(n > m) swap(n,m);
		int ans = 1;
		for(int l = 1, r; l <= n; l = r+1)
		{
			r = min(n/(n/l),m/(m/l));
			ans = ans * ksm(ksm(sum[r] * ksm(sum[l-1],p-2) % p,n/l),m/l) % p; 
		}
		printf("%lld\n",ans);
	}
	return 0;
}