[BZOJ4816][Sdoi2017]数字表格 数学

考虑每个$f_i$对答案的贡献，就能得到式子

$$\sum_{k=1}^nf_k^{\sum_{d|k}\mu (d)\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor}$$
转化成枚举 $kd$，令 $T=kd$
$$\prod_{T=1}^n(\prod_{d|T}f_d^{\mu(\frac{T}{d})})^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor}$$
预处理小括号内的部分，令 $g_T=\prod_{d|T}f_d^{\mu(\frac{T}{d})}$
单次询问只需要回答
$$\prod_{T=1}^ng_T^{\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor}$$
易证 $\lfloor \frac{n}{T}\rfloor$和 $\lfloor \frac{m}{T}\rfloor$只有 $O(\sqrt{N})$种取值，分段计数即可
预处理 $O(N\sqrt{N})$，询问 $O(T\sqrt{N}log_2N)$，log来自快速幂

#include <bits/stdc++.h>
#define mod 1000000007
#define N 1000050
#define tp 1000000
using namespace std;
typedef long long LL;
inline int rd() {int r;scanf("%d",&r);return r;}

bool np[N];
int eu[N],pr[N],g[N],fib[N],sum[N],siv[N],cnt,n,m,ans;
inline int qp(int a,long long b) {
    int ret = 1;
    while (b) {
        if (b&1) ret = 1LL * ret * a % mod;
        b >>= 1, a = 1LL * a * a % mod;
    }
    return ret;
}

void preset() {
    eu[1] = 1;
    for (int i=2;i<=tp;i++) {
        if (!np[i]) pr[++cnt] = i, eu[i] = -1;
        for (int j=1;j<=cnt && pr[j]*i<=tp;j++) {
            np[i*pr[j]] = 1;
            if (i%pr[j] == 0) {
                eu[i*pr[j]] = 0;
                break;
            }
            eu[i*pr[j]] = -eu[i];
        }
    }   


    fib[0] = 0, fib[1] = 1;
    for (int i=2;i<=tp;i++) fib[i] = (fib[i-1] + fib[i-2]) % mod;
    for (int i=1;i<=tp;i++) g[i] = 1;   

    sum[0] = siv[0] = 1;
    for (int i=1;i<=tp;++i) {
        int cur = fib[i];
        int inv = qp(cur, mod-2);
        for (int j=1;i*j<=tp;++j) {
            if (eu[j] == 0) continue;
            g[i*j] = eu[j] == 1 ?
                1LL * g[i*j] * cur % mod:
                1LL * g[i*j] * inv % mod;
        }
        sum[i] = 1LL * sum[i-1] * g[i] % mod;
        siv[i] = 1LL * siv[i-1] * qp(g[i], mod-2) % mod;
    }

}

void solve() {
    n = rd(), m = rd(), ans = 1;
    if (n>m) swap(n,m);
    for (int i=1,lst,cur;i<=n;i=lst+1) {
        lst = min(n/(n/i), m/(m/i));
        cur = 1LL * sum[lst] * siv[i-1] % mod;
        ans = 1LL * ans * qp(cur, 1LL*(n/i)*(m/i)) % mod;
    }
    printf("%d\n",ans);
}

int main() {
    freopen("product.in","r",stdin);
    freopen("product.out","w",stdout);
    preset();
    for (int T=rd();T;T--) solve();
    return 0;
}