P3704 [SDOI2017]数字表格

题目描述

洛谷

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

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

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

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

数据范围： (Tleq 10^3, 1leq n,mleq 10^6)

solution

莫比乌斯反演好题。

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

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

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

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

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

(=displaystylesum_{i=1}^{nover d}sum_{j=1}^{mover d} [gcd(i,j) = 1])

(= displaystylesum_{i=1}^{nover d}sum_{j=1}^{mover d} sum_{pmid gcd(i,j)} mu(p))

(= displaystylesum_{p=1}^{nover d} {nover dp} {mover dp} mu(p))

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

(large displaystyleprod_{d=1}^{n} f_d ^{sum_{p=1}^{nover d} {nover dp} {mover dp} mu(p)})

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

(large displaystyleprod_{Q=1} left(prod_{dmid Q} f_d^{mu({Qover d})} ight)^{{nover Q} {mover Q}})

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

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

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

复杂度： (O(Tsqrt{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
",ans);
	}
	return 0;
}