P3704 [SDOI2017]数字表格

\(\color{#0066ff}{ 题目描述 }\)

Doris刚刚学习了fibonacci数列。用\(f[i]\)表示数列的第\(i\)项，那么

\(f[0]=0\),\(f[1]=1\),

\(f[n]=f[n-1]+f[n-2]\),\(n\geq 2\)

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

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

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

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

\(\color{#0066ff}{输入格式}\)

有多组测试数据。

第一个一个数\(T\)，表示数据组数。

接下来\(T\)行，每行两个数\(n,m\)

\(\color{#0066ff}{输出格式}\)

输出\(T\)行，第\(i\)行的数是第\(i\)组数据的结果

\(\color{#0066ff}{输入样例}\)

3
2 3
4 5
6 7

\(\color{#0066ff}{输出样例}\)

1
6
960

\(\color{#0066ff}{数据范围与提示}\)

对\(10\%\)的数据，\(1\leq n,m\leq 100\)

对\(30\%\)的数据，\(1\leq n,m\leq 1000\)

另外存在\(30\%\)的数据，\(T\leq 3\)

对\(100\%\)的数据，\(T\leq1000,1\leq n,m\leq 10^6\)

时间限制：5s

内存限制：128MB

\(\color{#0066ff}{ 题解 }\)

题目要求这东西

\[\prod_{i=1}^n \prod_{j=1}^m f[gcd(i,j)] \]

枚举gcd

\[\prod_{d=1}^{min(n,m)}\prod_{i=1}^n \prod_{j=1}^m [gcd(i,j)==d]f(d) \]

改成枚举gcd然后变成\(\sum\),并放到指数位置

\[\prod_{d=1}^{min(n,m)}f(d)^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} [gcd(i,j)==1]} \]

拿\(\mu\)换一下

\[\prod_{d=1}^{min(n,m)}f(d)^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{k|gcd(i,j)} \mu(k)} \]

\[\prod_{d=1}^{min(n,m)}f(d)^{\sum_{k=1}^{min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor} \mu(k) * \lfloor\frac{n}{kd}\rfloor*\lfloor\frac{m}{kd}\rfloor} \]

kd换q

\[\prod_{q=1}^{min(n,m)}\prod_{d|q}f(d)^{ \mu(\frac q d) * \lfloor\frac{n}{q}\rfloor*\lfloor\frac{m}{q}\rfloor} \]

\[\prod_{q=1}^{min(n,m)}\prod_{d|q}(f(d)^{ \mu(\frac q d) }) ^ { \lfloor\frac{n}{q}\rfloor*\lfloor\frac{m}{q}\rfloor} \]

发现括号的部分可以\(O(nlogn)\)处理出来，然后数列分块就行了

#include<bits/stdc++.h>
#define LL long long
LL in() {
	char ch; LL x = 0, f = 1;
	while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
	for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
	return x * f;
}
const int maxn = 1e6 + 10;
const int mod = 1e9 + 7;
LL f[maxn], mu[maxn], pri[maxn], tot, h[maxn];
bool vis[maxn];
LL ksm(LL x, LL y) {
	LL re = 1LL;
	while(y) {
		if(y & 1) re = re * x % mod;
		x = x * x % mod;
		y >>= 1;
	}
	return re;
}
int mi(LL a, LL b) {
	if(b == 1) return a;
	if(b == 0) return 1;
	if(b == -1) return ksm(a, mod - 2);
	assert(0);
	return 2333333;
}
void predoit() {
	mu[1] = 1;
	for(int i = 2; i < maxn; i++) {
		if(!vis[i]) pri[++tot] = i, mu[i] = -1;
		for(int j = 1; j <= tot && (LL)i * pri[j] < maxn; j++) {
			vis[i * pri[j]] = true;
			if(i % pri[j] == 0) break;
			else mu[i * pri[j]] = -mu[i];
		}
	}
	f[0] = 0, f[1] = 1;
	for(int i = 2; i < maxn; i++) f[i] = (f[i - 1] + f[i - 2]) % mod;
	for(int i = 0; i < maxn; i++) h[i] = 1;
	for(int i = 1; i < maxn; i++)
		for(int j = i; j < maxn; j += i) 
			h[j] = h[j] * mi(f[i], mu[j / i]) % mod;
	for(int i = 2; i < maxn; i++) h[i] = (h[i] * h[i - 1]) % mod;
}
LL work(LL n, LL m) {
	LL ans = 1;
	for(LL l = 1, r; l <= std::min(n, m); l = r + 1) {
		r = std::min(n / (n / l), m / (m / l));
		ans = ans * ksm(h[r] * ksm(h[l - 1], mod - 2) % mod, (n / l) * (m / l)) % mod;
	}
	return ans;
}
int main() {
	predoit();
	for(int T = in(); T --> 0;) printf("%lld\n", work(in(), in()));
	return 0;
}