[数论][莫比乌斯反演] BZOJ 4816: 数字表格

$Description$

求

$$\prod_{i=1}^n\prod_{i=1}^mF_{(i,j)}$$
$1\le n,m\le10^6$。

$Solution$

推一下柿子：

$$\begin{eqnarray} \prod_{i=1}^n\prod_{i=1}^mF_{(i,j)}&=&\prod_{d=1}^nF_d^{\sum_{i=1}^n\sum_{j=1}^m[(i,j)=d]}\\ &=&\prod_{d=1}^nF_d^{\sum_{k=1}^{\lfloor{n\over d}\rfloor}\mu(k){\lfloor{n\over{dk}}\rfloor}{\lfloor{m\over{dk}}\rfloor}}\\ &=&\prod_{d=1}^nF_d^{\sum_{D=1}^n\mu({D\over k}){\lfloor{n\over{D}}\rfloor}{\lfloor{m\over{D}}\rfloor}}\\ &=&\prod_{D = 1}^{n}( \prod_{d\mid D}F(d)^{\mu(\frac{T}{d})})^{{\lfloor{n\over{D}}\rfloor}{\lfloor{m\over{D}}\rfloor}} \end{eqnarray}$$
设 $$G(D)=\prod_{d\mid D}F(d)^{\mu(\frac{T}{d})}$$ 所以 $$\prod_{i=1}^n\prod_{i=1}^mF_{(i,j)}=\prod_{D = 1}^{n} G(D)^{{\lfloor{n\over{D}}\rfloor}{\lfloor{m\over{D}}\rfloor}}$$
$G$可以$O(n\ln n)$时间预处理出来。 $F$也可以直接预处理出来。
这里有一个小$trick$：
因为每次累乘计算 $F_i^{-1}$的时候都要乘上一个 $O(\log P)$的时间，是不能跑过去的。所以可以先算出 $F_i$的前缀积 $T_i$。那么就有 $$F_i^{-1}=T_i^{-1}T_{i-1},T_i^{-1}=T_{i+1}^{-1}F_{i+1}$$ 就可以 $O(n)$计算出来啦。

#include <bits/stdc++.h>
using namespace std;

const int N = 1010101;
const int MOD = 1000000007;
typedef long long ll;

inline char get(void) {
    static char buf[100000], *S = buf, *T = buf;
    if (S == T) {
        T = (S = buf) + fread(buf, 1, 100000, stdin);
        if (S == T) return EOF;
    }
    return *S++;
}
inline void read(int &x) {
    static char c; x = 0;
    for (c = get(); c < '0' || c > '9'; c = get());
    for (; c >= '0' && c <= '9'; c = get()) x = x * 10 + c - '0';
}

int f[N], fi[N], pref[N], prefi[N];
int g[N], gi[N], preg[N], pregi[N];
int prime[N], vis[N], mu[N];
int lim = 1000000;
int n, m, Pcnt, x, y, test, pos, ans;

inline void Add(int &x, int t) {
    x += t; while (x >= MOD) x -= MOD;
}
inline int Pow(int a, ll b) {
    int c = 1;
    while (b) {
        if (b & 1) c = (ll)c * a % MOD;
        b >>= 1; a = (ll)a * a % MOD;
    }
    return c;
}
inline int Inv(int x) {
    return Pow(x, MOD - 2);
}

int main(void) {
    freopen("1.in", "r", stdin);
    mu[1] = 1;
    for (int i = 2; i <= lim; i++) {
        if (!vis[i]) {
            prime[++Pcnt] = i; mu[i] = -1;
        }
        for (int j = 1; j <= Pcnt && (x = i * prime[j]) <= lim; j++) {
            vis[x] = 1;
            if (i % prime[j]) {
                mu[x] = -mu[i];
            } else {
                mu[x] = 0; break;
            }
        }
    }
    pref[0] = prefi[0] = 1; f[0] = fi[0] = 0;
    pref[1] = prefi[1] = f[1] = fi[1] = 1;
    for (int i = 2; i <= lim; i++) {
        Add(f[i] = f[i - 1], f[i - 2]);
        pref[i] = (ll)pref[i - 1] * f[i] % MOD;
    }
    prefi[lim] = Inv(pref[lim]);
    for (int i = lim; i >= 2; i--) {
        prefi[i - 1] = (ll)prefi[i] * f[i] % MOD;
        fi[i] = (ll)prefi[i] * pref[i - 1] % MOD;
    }
    for (int i = 1; i <= lim; i++) g[i] = gi[i] = 1;
    for (int i = 1; i <= lim; i++)
        for (int j = i; j <= lim; j += i) {
            if (mu[j / i] == 1) {
                x = f[i]; y = fi[i];
            } else if (mu[j / i] == 0) {
                x = y = 1;
            } else {
                x = fi[i]; y = f[i];
            }
            g[j] = (ll)g[j] * x % MOD;
            gi[j] = (ll)gi[j] * y % MOD;
        }
    preg[0] = pregi[0] = 1;
    for (int i = 1; i <= lim; i++) {
        preg[i] = (ll)preg[i - 1] * g[i] % MOD;
        pregi[i] = (ll)pregi[i - 1] * gi[i] % MOD;
    }
    read(test);
    while (test--) {
        read(n); read(m); ans = 1;
        if (n > m) swap(n, m);
        for (int i = 1; i <= n; i = pos + 1) {
            pos = min(n / (n / i), m / (m / i));
            ans = (ll)ans * Pow((ll)preg[pos] * pregi[i - 1] % MOD, (ll)(n / i) * (m / i)) % MOD;
        }
        printf("%d\n", ans);
    }
    return 0;
}