P3704 [SDOI2017]数字表格(莫比乌斯反演)

P3704 [SDOI2017]数字表格(莫比乌斯反演)

题目描述

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

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

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

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

数据范围

• 对于 \(10\%\) 的数据，保证 \(n,m\leq 10^2\)。
• 对于 \(30\%\) 的数据，保证 \(n,m\leq 10^3\)。
• 另有 \(30\%\) 的数据，保证 \(T\leq 3\)。
• 对于 \(100\%\) 的数据，保证 \(1 \leq T\leq 10^3\)，\(1\leq n,m\leq 10^6\)。

解题思路

推一波式子，设 \(n \geq m\)

淦，推了半天发现都写成 \(\sum\) 了。。。不过没啥影响就是了

\[\large Ans = \sum_{i=1}^n\sum_{j=1}^mF_{\gcd(i,j)}\\ \large = \sum_{k=1}^{n}F_k\sum_{i=1}^{\frac nk}\sum_{j=1}^{\frac mk}[gcd(i,j)=1]\\ \large = \sum_{k=1}^{n}F_k \sum_{d=1}^{\frac nk}\mu(d)\sum_{i=1}^{\frac n{kd}}\sum_{j=1}^{\frac m{kd}}1\\ \large = \sum_{k=1}^{n}F_k \sum_{d=1}^{\frac nk}\mu(d) \lfloor \frac n{kd} \rfloor\lfloor \frac m{kd} \rfloor\\ \large = \prod_{k=1}^nF_k^{\sum_{d=1}^{\frac nk}\mu(d) \lfloor \frac n{kd} \rfloor\lfloor \frac m{kd} \rfloor}\\ \large = \prod_{T=1}^n(\prod_{k|T} F_k^{\mu(\frac Tk)})^{\lfloor \frac n{T} \rfloor\lfloor \frac m{T} \rfloor}\\ \]

括号里的东西可以预处理，然后做整除分块即可

数论题都好难啊 /kk

#include <queue>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define MP make_pair
#define ll long long
#define fi first
#define se second
using namespace std;

template <typename T>
void read(T &x) {
    x = 0; bool f = 0;
    char c = getchar();
    for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
    for (;isdigit(c);c=getchar()) x=x*10+(c^48);
    if (f) x=-x;
}

template<typename F>
inline void write(F x)
{
	static short st[30];short tp=0;
	if(x<0) putchar('-'),x=-x;
	do st[++tp]=x%10,x/=10; while(x);
	while(tp) putchar('0'|st[tp--]);
	putchar('\n');
}

template <typename T>
inline void Mx(T &x, T y) { x < y && (x = y); }

template <typename T>
inline void Mn(T &x, T y) { x > y && (x = y); }

const int N = 2005000;
const int P = 1e9+7;
int e[N], prime[N], mu[N], tot;
ll f[N], g[N], gv[N], T;

ll fpw(ll x, ll mi) {
	ll res = 1;
	while (mi) {
		if (mi & 1) res = res * x % P;
		x = x * x % P, mi >>= 1;
	}
	return res;
}

void prework(void) {
	const int N = 1000000;
	f[1] = mu[1] = g[1] = 1;
	for (int i = 2;i <= N; i++) {
		if (!e[i]) prime[++tot] = e[i] = i, mu[i] = -1;
		for (int j = 1;j <= tot && prime[j] * i <= N; j++) {
			int t = prime[j] * i; e[t] = prime[j];
			if (prime[j] == e[i]) break; mu[t] = -mu[i];
		}
		f[i] = (f[i-1] + f[i-2]) % P, g[i] = 1;
	}
	g[0] = gv[0] = 1;
	for (int i = 1;i <= N; i++) {
		int G[3] = { fpw(f[i], P - 2), 1, f[i] };
		for (int j = 1; j * i <= N; j++)
			g[j * i] = g[j * i] * G[mu[j] + 1] % P;
		gv[i] = fpw(g[i], P - 2) * gv[i-1] % P, g[i] = g[i] * g[i-1] % P;
	}
}

int main() {
	for (prework(), read(T); T; T--) {
		int m, n; read(n), read(m);
		ll ans = 1; if (n > m) swap(n, m);
//		for (int i = 1;i <= 5; i++) write(g[i]);
		for (int l = 1, r; l <= n; l = r + 1) {
			r = min(n / (n / l), m / (m / l));
			ans = ans * fpw(g[r] * gv[l-1] % P, (ll)(n/l) * (m/l) % (P - 1)) % P;
		}
		write(ans);
	}
	return 0;
}