P3704 [SDOI2017]数字表格——（莫比乌斯反演）

题目链接：https://www.luogu.org/problem/P3704

题目大意

Doris 刚刚学习了 fibnacci 数列，用 f[i] f[i]f[i] 表示数列的第 i ii 项，那么：

f[0]=0f[1]=1f[n]=f[n−1]+f[n−2],n≥2 \begin{aligned} f[0] &= 0 \\ f[1] &= 1 \\ f[n] &= f[n - 1] + f[n - 2], n \geq 2 \end{aligned}​f[0]​f[1]​f[n]​​​=0​=1​=f[n−1]+f[n−2],n≥2​​

Doris 用老师的超级计算机生成了一个 n×m n \times mn×m 的表格，第 i ii 行第 j jj 列的格子中的数是 f[gcd(i,j)] f[\gcd(i, j)]f[gcd(i,j)]，其中 gcd(i,j) \gcd(i, j)gcd(i,j) 表示 i ii 与 j jj 的最大公约数。

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

这些数的乘积实在是太大了，所以 Doris 只想知道乘积对 1000000007 10000000071000000007 取模后的结果。

输入格式

有多组测试数据。
第一行一个数 T TT，表示数据组数。
接下来 T TT 行，每行两个数 n nn 和 m mm。

输出格式

输出 T TT 行，第 i ii 行的数是第 i ii 组数据的结果。

样例

样例输入

3
2 3
4 5
6 7

样例输出

1
6
960

数据范围与提示

对于 10% 10\%10% 的数据，1≤n,m≤100 1 \leq n, m \leq 1001≤n,m≤100；
对于 30% 30\%30% 的数据，1≤n,m≤1000 1 \leq n, m \leq 10001≤n,m≤1000；
另外存在 30% 30\%30% 的数据，T≤3 T \leq 3T≤3；
对于 100% 100\%100% 的数据，1≤n,m≤106,1≤T≤1000 1 \leq n, m \leq 10 ^ 6, 1 \leq T \leq 10001≤n,m≤10​6​​,1≤T≤1000。

思路

推式子，个人觉得还是比较好推的（虽然当初没推出来

反正就尽量往找gcd个数的模板上面靠吧

我就不误人子弟了，贴个大神博客

点我看大神博客

#include<cstdio>
#include<cstring>
#include<iostream>
#include <algorithm>
#define int long long 

typedef long long ll;
using namespace std;

const int N = 1000000;
const ll p = 1e9 + 7;
int sshu[N + 5], tot = 0, mu[N + 5], n, m;
int f[N + 5], sum[N + 5], inv[N + 5];
bool no[N + 5];

inline int KSM(int a, int b) {
	int t = 1;
	a %= p;
	while (b) {
		if (b & 1) t = (1LL * t * a) % p;
		b >>= 1;
		a = (1LL * a * a) % p;
	}
	return t % p;
}

void prepare() {
	mu[1] = 1;
	for (int i = 2; i <= N; i++) {
		if (!no[i]) {
			sshu[++tot] = i;
			mu[i] = -1;
		}
		for (int j = 1; j <= tot && sshu[j] * i <= N; j++) {
			no[sshu[j] * i] = 1;
			if (i % sshu[j] == 0) {
				mu[sshu[j] * i] = 0;
				break;
			}
			mu[sshu[j] * i] = -mu[i];
		}
	}

	f[0] = 0; f[1] = 1;
	sum[0] = 1; sum[1] = 1;                               //前缀积
	inv[1] = 1;                                         //f的逆元 
	for (int i = 2; i <= N; i++) f[i] = (f[i - 1] + f[i - 2]) % p, inv[i] = KSM(f[i], p - 2), sum[i] = 1;
	for (int i = 1; i <= N; i++)
		if (mu[i] != 0) {
			for (int j = 1; j * i <= N; j++)                  // j*i=T i=T/d j=d
				if (mu[i] > 0) sum[j * i] = (1LL * sum[j * i] * f[j]) % p;
				else sum[j * i] = (1LL * sum[j * i] * inv[j]) % p;
		}
	for (int i = 2; i <= N; i++) sum[i] = (1LL * sum[i] * sum[i - 1]) % p;
}

signed main()
{
	prepare();
	int T;
	cin >> T;

	while (T--) {
		cin >> n >> m;
		int last, ans = 1;
		for (int i = 1; i <= min(n, m); i = last + 1) {
			last = min(n / (n / i), m / (m / i));
			ans = 1LL * ans % p * KSM(1LL * sum[last] * KSM(sum[i - 1], p - 2) % p, 1LL * (n / i) * (m / i) % (p - 1)) % p;
		}
		cout << ans % p << endl;
	}
	return 0;
}