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

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}

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

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

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

题解

反演一下就好啦。

设 n n 为 n,m n, m 中的较小值，首先，题目要求求的式子是

∏i=1n∏j=1mf(gcd((i,j)) \prod\limits_{i = 1}^{n} \prod\limits_{j = 1}^{m} f(\gcd((i, j))

考虑枚举一个 d d 使得 d d 作为 gcd(i,j) \gcd(i, j) 的值，原式可化为

∏d=1nf(d)∑i=1n∑j=1m[gcd(i,j)=d] \prod\limits_{d = 1}^{n} f(d) ^ {\sum\limits_{i = 1}^{n} \sum\limits_{j = 1}^{m} [\gcd(i, j) = d]}

把幂指数单独拿下来处理

∑i=1n∑j=1m[gcd(i,j)=d] \sum\limits_{i = 1}^{n} \sum\limits_{j = 1}^{m} [\gcd(i, j) = d]

将 i,j i, j 同时除以 d d 得

∑i=1n∑j=1m[gcd(⌊id⌋⌊jd⌋)=1] \sum\limits_{i = 1}^{n} \sum\limits_{j = 1}^{m} [\gcd(\lfloor \frac i d \rfloor \lfloor \frac j d \rfloor) = 1]

更改求和指标，设 N=⌊nd⌋,M=⌊md⌋ N = \lfloor \frac n d \rfloor, M = \lfloor \frac m d \rfloor

∑i=1N∑j=1M[gcd(i,j)=1] \sum\limits_{i = 1}^{N} \sum\limits_{j = 1}^{M} [\gcd(i, j) = 1]

把 [gcd(i,j)=1] [\gcd(i, j) = 1] 卷一下

∑i=1N∑j=1M∑p[p∣gcd(i,j)]μ(p) \sum\limits_{i = 1}^{N} \sum\limits_{j = 1}^{M} \sum\limits_{p} [p | \gcd(i, j)]\mu(p)

更改求和顺序，改为先枚举 p p

∑p=1N∑i=1N∑j=1M[p∣gcd(i,j)] \sum\limits_{p = 1}^{N} \sum\limits_{i = 1}^{N} \sum\limits_{j = 1}^{M} [p | \gcd(i, j)]

p ∣gcd(i,j) p\ | \gcd(i, j) 的充要条件是 p ∣ i∧p ∣ j p\ |\ i \wedge p\ |\ j

∑p=1Nμ(p)∑i=1N[p ∣ i]∑j=1M[p ∣ j] \sum\limits_{p = 1}^{N} \mu(p) \sum\limits_{i = 1}^{N}[p\ |\ i]\sum\limits_{j = 1}^{M} [p\ |\ j]

可以转化为

∑p=1Nμ(p)⌊Np⌋⌊Mp⌋ \sum\limits_{p = 1}^{N} \mu(p) \lfloor \frac N p \rfloor \lfloor \frac M p \rfloor

预处理 μ \mu 的前缀和，数论分块求。

回到一开始的式子，上式可以看做一个函数 g(N,M) g(N, M) ，其中 N=⌊nd⌋,M=⌊md⌋ N = \lfloor \frac n d \rfloor, M = \lfloor \frac m d \rfloor 。

这样，预处理斐波那契数列前缀积，对于指数上的这个函数，也可以进行数论分块。

代码

#include <cstdio>
#include <cstring>
#include <climits>
#include <algorithm>
#define min(a, b) ((a) < (b) ? (a) : (b))
const int MAX_N = 1e6;
const int MOD = 1e9 + 7;
int primes[MAX_N + 1], mu[MAX_N + 1], mus[MAX_N + 1], cnt;
bool isNotPrime[MAX_N + 1];
inline void sieve() {
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    for (int i = 2; i <= MAX_N; i++) {
        if (!isNotPrime[i]) {
            primes[cnt++] = i;
            mu[i] = -1;
        }
        for (int j = 0; j < cnt && (long long)i * primes[j] <= MAX_N; j++) {
            int p = primes[j];
            isNotPrime[i * p] = true;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            } else {
                mu[i * p] = -mu[i];
            }
        }
    }
    for (int i = 1; i <= MAX_N; i++) mus[i] = mus[i - 1] + mu[i];
}
inline void exgcd(long long a, long long b, long long &g, long long &x, long long &y) {
    if (!b) g = a, x = 1, y = 0;
    else exgcd(b, a % b, g, y, x), y -= x * (a / b);
}
inline long long inv(long long x) {
    long long r, g, tmp;
    exgcd(x, MOD, g, r, tmp);
    return (r % MOD + MOD) % MOD;
}
long long fib[MAX_N + 1], fibProd[MAX_N + 1], fibProdInv[MAX_N + 1];
inline void prepare() {
    fib[0] = 0, fib[1] = 1;
    for (int i = 2; i <= MAX_N; i++) fib[i] = (fib[i - 1] + fib[i - 2]) % MOD;
    fibProd[0] = fibProdInv[0] = 1;
    for (int i = 1; i <= MAX_N; i++) {
        fibProd[i] = fibProd[i - 1] * fib[i] % MOD;
        fibProdInv[i] = inv(fibProd[i]);
    }
}
inline long long g(int N, int M) {
    long long res = 0;
    for (int l = 1, r; l <= N; l = r + 1) {
        r = min(N / (N / l), M / (M / l));
        res += (long long)(mus[r] - mus[l - 1]) * (N / l) * (M / l);
    }
    return res;
}
inline int pow(long long a, long long n) {
    long long res = 1;
    for (; n; n >>= 1, a = a * a % MOD) if (n & 1) res = res * a % MOD;
    return res;
}
inline int calc(int N, int M, int l, int r) {
    return pow(fibProd[r] * fibProdInv[l - 1] % MOD, g(N, M));
}
inline int solve(int n, int m) {
    if (n > m) std::swap(n, m);
    long long res = 1;
    for (int l = 1, r; l <= n; l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        (res *= calc(n / l, m / l, l, r)) %= MOD;
    }
    return res;
}
int main() {
    freopen("product.in", "r", stdin);
    freopen("product.out", "w", stdout);
    int T_T;
    scanf("%d", &T_T);
    sieve();
    prepare();
    while (T_T--) {
        int n, m;
        scanf("%d %d", &n, &m);
        printf("%d\n", solve(n, m));
    }
}
1. 1. 1. 1.