P3704 [SDOI2017]数字表格

P3704 [SDOI2017]数字表格

链接

分析：

$ prodlimits_{i = 1}^{n} prodlimits_{j = 1}^{m} f[gcd(i, j)]$

$=prodlimits_{d = 1}^{min(n, m)} prodlimits_{i = 1}^{n} prodlimits_{j = 1}^{m} [gcd(i, j) = d] imes f[d]$

$=prodlimits_{d = 1}^{min(n, m)} f[d] ^ {sumlimits_{i = 1}^{n} sumlimits_{j = 1}^{m} [gcd(i, j) = d]}$

$=prodlimits_{d = 1}^{min(n, m)} f[d] ^ {sumlimits_{k = 1}^{min( frac{n}{d} , frac{m}{d} )} mu(k) frac{n}{kd} frac{m}{kd}}$

设$T=kd$

$prodlimits_{T = 1} ^ {min(n, m)} (prodlimits_{d | T} f[d] ^ {mu(frac{T}{d}) } ) ^ {frac{n}{T} frac{m}{T} }$

对中间的部分$nlogn$预处理，$O(sqrt n)$处理每个询问。

代码：

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#include<cmath>
#include<cctype>
#include<set>
#include<queue>
#include<vector>
#include<map>
using namespace std;
typedef long long LL;

inline int read() {
    int x=0,f=1;char ch=getchar();for(;!isdigit(ch);ch=getchar())if(ch=='-')f=-1;
    for(;isdigit(ch);ch=getchar())x=x*10+ch-'0';return x*f;
}

const int N = 2000005, mod = 1e9 + 7;
int mu[N], pri[N], f[N], g1[N], g2[N], inv1[N], inv2[N];
bool nopri[N];

int ksm(int a,LL b) {
    register int res = 1;
    while (b) {
        if (b & 1) res = 1ll * res * a % mod;
        a = 1ll * a * a % mod;
        b >>= 1;
    }
    return res % mod;
}
void init(int n) {
    nopri[1] = true; mu[1] = 1;
    int tot = 0;
    for (int i = 2; i <= n; ++i) {
        if (!nopri[i]) pri[++tot] = i, mu[i] = -1;
        for (register int j = 1; j <= tot && pri[j] * i <= n; ++j) {
            nopri[pri[j] * i] = 1; 
            if (i % pri[j] == 0) { mu[i * pri[j]] = 0; break; }
            mu[pri[j] * i] = -mu[i];
        }
    }
    f[0] = 0, f[1] = 1; g1[1] = g2[1] = 1;
    for (register int i = 2; i <= n; ++i) f[i] = (f[i - 1] + f[i - 2]) % mod, g1[i] = g2[i] = 1;
    for (int i = 1; i <= n; ++i) 
        for (int j = i; j <= n; j += i) {
            if (mu[j / i] == 0) continue;
            else if (mu[j / i] == 1) g1[j] = 1ll * g1[j] * f[i] % mod;
            else g2[j] = 1ll * g2[j] * f[i] % mod;
        }
    g1[0] = g2[0] = inv1[0] = inv2[0] = 1;
    for (int i = 1; i <= n; ++i) {
        g1[i] = 1ll * g1[i] * g1[i - 1] % mod,
        g2[i] = 1ll * g2[i] * g2[i - 1] % mod;
        inv1[i] = ksm(g1[i], mod - 2);
        inv2[i] = ksm(g2[i], mod - 2);
    }
}
void solve() {
    int n = read(), m = read(), nm = min(n, m), pos = 0, ans = 1;
    for (int t1, t2, i = 1; i <= nm; i = pos + 1) {
        pos = min(n / (n / i), m / (m / i));
        LL t = 1ll * (n / i) * (m / i); // !!!
        t1 = 1ll * g1[pos] * inv1[i - 1] % mod; 
        t2 = 1ll * g2[pos] * inv2[i - 1] % mod;
        ans = 1ll * ans * ksm(t1, t) % mod * ksm(ksm(t2, t), mod - 2) % mod;
    }
    cout << ans << "
";
}
int main() {
    init(1000000);
    for (int T = read(); T --; solve());
    return 0;
}