[Sdoi2017]数字表格 莫比乌斯反演

Description
$Doris$刚刚学习了$Fibonacci$数列。用$f[i]$表示数列的第i项，那么
$f[0]=0,f[1]=1,f[n]=f[n-1]+f[n-2],n>=2$
$Doris$用老师的超级计算机生成了一个$n\times m$的表格，第$i$行第$j$列的格子中的数是$f[gcd(i,j)]$，其中$gcd(i,j)$表示$i,j$的最大公约数。$Doris$的表格中共有$n\times m$个数，她想知道这些数的乘积是多少。答案对$10^9+7$取模。

---

Sample Input
3
2 3
4 5
6 7

---

Sample Output
1
6
960

---

莫反一下：
$$\prod _{d=1}^{n}f(d{)}^{{\sum }_{i=1}^{\frac{n}{d}}{\sum }_{j=1}^{\frac{m}{d}}{\sum }_{k|i,k|j}\mu (k)}$$
$=>$
$$\prod _{d=1}^{n}f(d{)}^{{\sum }_{k=1}^{\frac{n}{d}}\mu (k)\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor }$$
然后用一个老套路，设$T=kd$
$$\prod _{T=1}^{n}f(d)(\prod _{d|T}f(d{)}^{\mu (\frac{T}{d})}{)}^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor }$$
中间的式子可以$O(n\ast logn)$预处理，然后数论分块搞一下时间复杂度就是$O(T\ast sqrtn\ast logmod)$

---

#include <ctime>
#include <cmath>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>

using namespace std;
typedef long long LL;
int _max(int x, int y) {return x > y ? x : y;}
int _min(int x, int y) {return x < y ? x : y;}
const int mod = 1e9 + 7;
const int N = 1000001;
int read() {
	int s = 0, f = 1; char ch = getchar();
	while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
	while(ch >= '0' && ch <= '9') s = s * 10 + ch - '0', ch = getchar();
	return s * f;
}
void put(int x) {
	if(x >= 10) put(x / 10);
	putchar(x % 10 + '0');
}

bool v[N];
int plen, mu[N], p[N], f[N], inv[N];
LL s[N];

int add(int x, int y) {
	x += y;
	return x >= mod ? x - mod : x;
}

LL pow_mod(LL a, LL k) {
	LL ans = 1;
	while(k) {
		if(k & 1) (ans *= a) %= mod;
		(a *= a) %= mod; k /= 2;
	} return ans;
}

void get_p() {
	mu[1] = s[1] = s[0] = 1;
	for(int i = 2; i < N; i++) {
		if(!v[i]) p[++plen] = i, mu[i] = -1;
		for(int j = 1; j <= plen && (LL)i * p[j] < N; j++) {
			v[i * p[j]] = 1;
			if(i % p[j] == 0) {mu[i * p[j]] = 0; break;}
			else mu[i * p[j]] = -mu[i];
		} s[i] = 1;
	}
}

int main() {
	get_p();
	f[0] = 0, f[1] = 1, inv[1] = 1;
	for(int i = 2; i < N; i++) f[i] = add(f[i - 1], f[i - 2]), inv[i] = pow_mod(f[i], mod - 2);
	for(int i = 1; i < N; i++) {
		for(int j = i, o = 1; j < N; j += i, ++o) {
			(s[j] *= mu[o] == -1 ? inv[i] : pow_mod(f[i], mu[o])) %= mod;
		}
	} for(int i = 1; i < N; i++) (s[i] *= s[i - 1]) %= mod;
	int tt = read();
	while(tt--) {
		int n = read(), m = read();
		if(n > m) swap(n, m); LL ans = 1;
		for(int l = 1, r; l <= n; l = r + 1) {
			r = _min(n / (n / l), m / (m / l));
			LL hh = s[r] * pow_mod(s[l - 1], mod - 2) % mod;
			(ans *= pow_mod(hh, (LL)(n / l) * (m / l) % (mod - 1))) %= mod;
		} put(ans), puts("");
	} return 0;
}