题意
求如下表达式的值
\[\prod_{i=1}^{n} \prod_{j=1}^{m} f_{gcd(i,j)} \pmod{10^9 + 7}
\]
其中，\(f_i\)为 fibonacci 数列的第\(i\)项，\(n, m \leqslant 10^6\)
Solution
\[\prod_{i=1}^{n} \prod_{j=1}^{m} f_{gcd(i,j)}
\]
改变枚举顺序，优先枚举\(d = gcd(i,j)\)，
\[=\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]}
\]
取出指数，记：
\[g(d) = \sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{d} \rfloor} [gcd(i,j)=1]
\]
\[=\sum_{i=1}^{\lfloor \frac{n}{d} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{d} \rfloor} \sum_{x|gcd(i,j)} \mu(x)
\]
莫比乌斯反演经典处理手段，优先枚举因子，
\[= \sum_{x=1}^{min\{\lfloor \frac{n}{d} \rfloor , \lfloor \frac{m}{d} \rfloor \}} \mu(d) \lfloor \frac{n}{xd} \rfloor \lfloor \frac{m}{xd} \rfloor
\]
带回原表达式，则有，
\[\prod_{d=1}^{min\{n, m\}} \ \prod_{x=1}^{min\{\lfloor \frac{n}{d} \rfloor, \lfloor \frac{m}{d}\ \rfloor\}}f_d^{\mu(x) \lfloor \frac{n}{x\ d} \rfloor \lfloor \frac{m}{x\ d} \rfloor}
\]
可以比较容易发现，该式子可以通过两次数论分块解决，时间复杂度为\(O(T \ n \log n)\)，还需要继续优化
考虑将此式子转化为更加便于预处理\(\mu(x)\)的表达式，则记\(T = xd\)，
\[= \prod_{T=1}^{min\{n, m\}} \ \prod_{d|T} f_d^{\mu(\lfloor \frac{T}{d} \rfloor) \lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}
\]
改变一下幂指数的顺序，
\[= \prod_{T=1}^{min\{n, m\}} \left( \prod_{d|T} f_d^{\mu(\lfloor \frac{T}{d} \rfloor)} \right)^{\lfloor \frac{n}{T} \rfloor \lfloor \frac{m}{T} \rfloor}
\]
外面是我们熟悉的数论分块的形式，里面考虑到\(\mu(\lfloor \frac{T}{d} \rfloor)\)的取值只能为\(-1, 1, 0\)，可以\(O(n \log {n})\)筛出\(\mu\)的值并且计算，然后在\(O(n \sqrt{n})\)的时间内完成。
\(O(n \sqrt{n})\)？感觉在\(10^6\)下不能容纳这个时间复杂度，但是考虑到\(\mu(i) \ne 0\)的个数只有\(607926\)个，时间复杂度应该重新表达为\(O(n + 607926\sqrt{n})\)。
费马小定理与欧拉定理
费马小定理
若\(p\)为质数，\(gcd(a,p)=1\)，则\(a^{p-1} \equiv 1 \pmod{p}\)
欧拉定理
若\(gcd(a,m)=1\)则\(a^{\varphi(m)} \equiv 1 \pmod{m}\)
拓展欧拉定理
\[a^b \equiv \begin{cases} a^{b \bmod \varphi(m)}, & \text {$gcd(a,m)=1,$} \\
a^b, &\text {$gcd(a,m) \ne 1, b < \varphi(m),$} \\ a^{(b \bmod \varphi(m)) + \varphi(m)} , &\text {$gcd(a,m) \ne 1, b \geqslant \varphi(m).$} \end{cases} \pmod{m}\]
以上证明略。
因为本题要求对幂指数取模求解，并且本题的模数$m = 10^9 + 7 $，可以通过拓展欧拉定理加以求解。
Code
点击查看代码
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <iostream>
using namespace std;
const int N = 1e6 + 50;
const int mod = 1e9 + 7;
inline int read() {
int res = 0, f = 1;
char c = getchar();
while (c > '9' || c < '0') {
if (c == '-') f = -1;
c = getchar();
}
while (c >= '0' && c <= '9') {
res = (res << 3) + (res << 1) + c - '0';
c = getchar();
}
return res * f;
}
inline int add_mod(int a, int b) {
a += b;
if (a < 0) a += mod;
if (a >= mod) a -= mod;
return a;
}
inline int del_mod(int a, int b) {
a -= b;
if (a < 0) a += mod;
if (a >= mod) a -= mod;
return a;
}
inline int mul_mod(int a, int b) {
a = 1ll * a * b % mod;
return a;
}
int num;
int prime[N], mu[N];
int v[N];
int f[N], F[N], inv[N];
inline int qpow(int a, int b) {
int base = 1;
while (b) {
if (b & 1) base = mul_mod(base, a);
a = mul_mod(a, a);
b >>= 1;
}
return base;
}
void init() {
int MAX = 1e6;
F[0] = F[1] = f[1] = 1;
inv[0] = inv[1] = 1;
mu[1] = 1;
for (int i = 2; i <= MAX; ++i) {
f[i] = add_mod(f[i - 1], f[i - 2]);
inv[i] = qpow(f[i], mod - 2);
F[i] = 1;
if (!v[i]) {
prime[++num] = i;
v[i] = i;
mu[i] = -1;
}
for (int j = 1; j <= num; ++j) {
if (prime[j] > v[i] || prime[j] > MAX / i) break;
v[i * prime[j]] = prime[j];
if (i % prime[j] == 0) {
mu[i * prime[j]] = 0;
break;
}
mu[i * prime[j]] = -mu[i];
}
}
for (int i = 1; i <= MAX; ++i) {
if (!mu[i]) continue;
for (int j = i; j <= MAX; j += i) {
F[j] = mul_mod(F[j], (mu[i] == 1 ? f[j / i] : inv[j / i]));
}
}
for (int i = 2; i <= MAX; ++i) F[i] = mul_mod(F[i], F[i - 1]);
}
int n, m;
int ans;
int main () {
init();
int T = read();
while (T--) {
ans = 1;
n = read(), m = read();
for (int i = 1, j; i <= min(n, m); i = j + 1) {
j = min(n / (n / i), m / (m / i));
int tmp = mul_mod(F[j], qpow(F[i - 1], mod - 2));
ans = mul_mod(ans, qpow(tmp, 1ll * (n / i) * (m / i) % (mod - 1)));
}
printf("%d\n", ans);
}
return 0;
}