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

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

题目描述

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

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

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

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

数据范围

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

解题思路

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

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

[large Ans = sum_{i=1}^nsum_{j=1}^mF_{gcd(i,j)}\ large = sum_{k=1}^{n}F_ksum_{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} floorlfloor frac m{kd} floor\ large = prod_{k=1}^nF_k^{sum_{d=1}^{frac nk}mu(d) lfloor frac n{kd} floorlfloor frac m{kd} floor}\ large = prod_{T=1}^n(prod_{k|T} F_k^{mu(frac Tk)})^{lfloor frac n{T} floorlfloor frac m{T} floor}\ ]

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

数论题都好难啊 /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('
');
}

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;
}