[BZOJ4816]数字表格_bzoj4816数字表格-优快云博客

本文链接：https://blog.youkuaiyun.com/qq_31776811/article/details/70148197

题目链接：BZOJ4816

题目大意
求 $\prod_{1\le i \le n, 1 \le j \le m}F(\gcd(i,j))$ ，答案对 $1e9+7$ 取模，其中 $F$ 为 $fibonacci$ 数列， $F(0) = 0, F(1) = 1$ 。

分析
1. 枚举gcd，转化为求 $\prod_{d = 1}^nF(d)^{h(d)}$ ，其中 $h(d) = \sum_{i = 1}^n \sum_{j = 1} ^ m [\gcd(i,j) == d]$ 。
2. 而求 $h(d)$ 作为莫比乌斯反演的经典问题，我们知道结论，于是变成了 $\prod_{d = 1}^n F(d)^{(\sum_{k = 1}^{\lfloor \frac{n}{d} \rfloor}\mu(k)\lfloor \frac{n}{dk} \rfloor \lfloor \frac{m}{dk} \rfloor)}$ 。
3. 设 $T = dk$ ，同时外移T，就变成了 $\prod_{T = 1}^{n}( \prod_{d\mid T}F(d)^{\mu(\frac{T}{d})})^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor}$ 。
4. 设 $g(T) = \prod_{d\mid T}F(d)^{\mu(\frac{T}{d})}$ ，那么预处理出g(x)的前缀积以及前缀积的逆元即可 $\sqrt{N}\log N$ 求出每次询问。
5. 那么如何求 $g(x)$ 呢？首先初始化 $g(x) = 1$ ，我们枚举倍数 $d$ ，对所有 $d\mid T$ ， $g(T) = g(T)*F(d)^{\mu(\frac{T}{d})}$ ；因此，我们在求出 $g$ 之前，需要先筛出 $\mu$ ，求出 $F$ ，以及 $F$ 的逆元。
6. 然而因为模数过大，不能 $O(P)$ 预处理逆元，于是我们需要另一种方法求数列的逆元：对于数列 $a_1, a_2, a_3 ... a_n$ ，我们设 $S_i$ 为前i项积；那么我们可以 $O(logP)$ 求出 $S_n$ 的逆元 $invs[n] = (a_1a_2a_3...a_n)^{-1}$ ，那么 $a_n$ 的逆元 $inva[n] = S_{n-1} \cdot invs[n] = a_n^{-1}$ ，而且 $invs[n-1] = invs[n] \cdot a_n$ ，这样就可以递推出数列元素的逆元以及元素前缀积的逆元了。
7. 总结：预处理筛出 $\mu$ ，求出F以及其前缀积，然后求逆元；枚举因数求出g；然后求前缀积，求逆元；最后就可以 $O(T\sqrt N\log N)$ 过掉了。虽然步骤多，但是代码不长。

上代码

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
using namespace std;

typedef long long LL;
const int N = 1e6 + 10;
const int MOD = 1e9 + 7;

int n, m;
inline int read() {
    register char ch = getchar();
    register int ans = 0, neg = 1;
    for (; !isdigit(ch); ch = getchar())
        if (ch == '-') neg = -1;
    for (; isdigit(ch); ch = getchar())
        ans = ans * 10 + ch - '0';
    return ans * neg;
}

inline int pls(int a, int b) {
    LL ans = (LL)a + b;
    return ans >= (LL)MOD ? ans - MOD : ans;
}
inline int mul(int a, int b) {
    return (LL)a * b % MOD;
}
inline int power(int a, int b) {
    int ans = 1;
    while (b) {
        if (b & 1) ans = mul(ans, a);
        b >>= 1, a = mul(a, a);
    }
    return ans;
}

bool notPrime[N];
int prime[N], mu[N];
int fi[N], pf[N], invfi[N], invpf[N];
int gi[N], pg[N], invpg[N];
void preCalc(int a) {
    int cnt = 0, mut = 1;
    mu[1] = 1, fi[0] = 0, fi[1] = pf[1] = 1;
    for (int i = 2; i <= a; ++i) {
        if (!notPrime[i]) prime[++cnt] = i, mu[i] = -1;
        for (int j = 1, k; j <= cnt && (k = i * prime[j]) <= a; ++j) {
            notPrime[k] = true;
            if (i % prime[j] == 0) {
                mu[k] = 0; break;
            }
            mu[k] = -mu[i];
        }
        pf[i] = mul(pf[i - 1], fi[i] = pls(fi[i - 1], fi[i - 2]));
    }
    invpf[a] = power(pf[a], MOD - 2);
    for (int i = a; i >= 1; --i) {
        gi[i] = 1;
        invpf[i - 1] = mul(invpf[i], fi[i]);
        invfi[i] = mul(pf[i - 1], invpf[i]);
    }
    for (int i = 2; i <= a; ++i)
        for (int j = 1, k; (k = i * j) <= a; ++j)
            if (mu[j] == 1) gi[k] = mul(gi[k], fi[i]);
            else if (mu[j] == -1) gi[k] = mul(gi[k], invfi[i]);
    pg[1] = 1;
    for (int i = 2; i <= a; ++i)
        pg[i] = mul(pg[i - 1], gi[i]);
    invpg[a] = power(pg[a], MOD - 2);
    for (int i = a; i >= 1; --i)
        invpg[i - 1] = mul(invpg[i], gi[i]);
}

int main() {
    int T = read();
    preCalc(1e6);
    while (T--) {
        n = read(), m = read();
        if (n > m) swap(n, m);
        int ans = 1;
        for (int i = 1, nxt; i <= n; i = ++nxt) {
            nxt = min(n / (n / i), m / (m / i));
            ans = mul(ans, power(mul(pg[nxt], invpg[i - 1]), (LL)(n / i) * (m / i) % (MOD - 1)));
        }
        printf("%d\n", ans);
    }
    return 0;
}

以上