题解 P7890 「MCOI-06」Lost Desire

· · 题解

本题解旨在提供一种简洁的推导方法，减少各位的脑细胞的损耗。

在此题解中，我们约定：

F^ * (x)=\prod_{i=1}^xF(i) F^ {* * } (x)=\prod_{i=1}^xF^*(i)

即 F^ * (x) 相当于对 F(x) 做一次前缀积，F^{* * }(x) 就是对 F(x) 做两次前缀积。

那么我们就有 定理 1：

\prod_{i=1}^n\prod_{j=1}^mF(i+j)=\prod_{i=1}^n\dfrac{F^* (i+m)}{F^* (i)}=\dfrac{F^{* * }(n + m)}{F^{* * }(n)F^{* * }(m)}

同样地，在 \sum 中也有类似结论。

注：以下推导省略向下取整符号 \lfloor\rfloor。

接下来开始推式子，假设 n\le m，先看分子部分：

\begin{aligned} \prod_{i=1}^n\prod_{j=1}^m(i+j-1)!^{[i\perp j]}=\prod_{d=1}^n\prod_{i=1}^{n/d}\prod_{j=1}^{m/d}[d(i+j)-1]!^{\mu(d)} \end{aligned}

设 F_d(x)=(dx-1)!^{\mu(d)}，那么由定理 1 得：

\begin{aligned}\prod_{d=1}^n\prod_{i=1}^{n/d}\prod_{j=1}^{m/d}[d(i+j)-1]!^{\mu(d)}=\prod_{d=1}^n\dfrac{F^{* * }_ d (\frac nd+\frac md)}{F^{* * }_ d(\frac nd)F^{* * }_ d(\frac md)} \end{aligned}

这样分子就算推完了，接下来是分母：

\prod_{i=1}^n\prod_{j=1}^m(i!j!)^{[i\perp j]}=\prod_{d=1}^n\prod_{i=1}^{n/d}\prod_{j=1}^{m/d}[(di)!(dj)!]^{\mu(d)}

记 G_d(x)=(dx)!^{\mu(d)}，则：

\prod_{d=1}^n\prod_{i=1}^{n/d}\prod_{j=1}^{m/d}[(di)!(dj)!]^{\mu(d)}=\prod_{d=1}^nG_d^* (\frac nd)^{\frac md}G_d^ * (\frac md)^{\frac nd}

把分子分母结合起来，就可以得到：

\prod_{d=1}^n\dfrac{F^{* * }_ d (\frac nd+\frac md)}{F^{* * }_ d(\frac nd)F^{* * }_ d(\frac md)G_d^* (\frac nd)^{\frac md}G_d^ * (\frac md)^{\frac nd}}

为了能够整除分块，我们设：

H_d(x)=\prod_{i=1}^dF_i^{* * }(x) L_d(x)=\prod_{i=1}^dG_i^* (x)

那么对于整除分块中的区间 [l,r]，其计算式就是：

\dfrac{H_ r(\frac nl+\frac ml)H_ {l-1}(\frac nl)H_ {l-1}(\frac ml)L_{l-1} (\frac nl)^{\frac ml}L_{l-1}(\frac ml)^{\frac nl}}{H_ {l-1}(\frac nl+\frac ml)H_ r(\frac nl)H_ r(\frac ml)L_r (\frac nl)^{\frac ml}L_r (\frac ml)^{\frac nl}}

那么我们只需计算出所有需要的 H_d(x),L_d(x)，再加上整除分块，就能在 O(n\ln n+T\sqrt n\log n) 的时间复杂度内完成计算。

如果再加上 Index Calculus 或其他求离散对数的算法，那就可以通过此题了。

以下给出一个只能过前四个 Subtask 的代码，仅供参考 or 对拍：

#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define ll long long
using namespace std;

int Read(){
    int res = 0; char c = getchar();
    while(c < '0' || c > '9') c = getchar();
    while(c >= '0' && c <= '9') res = res * 10 + (c ^ 48), c = getchar();
    return res;
}
void Write(int x){ if(x > 9) Write(x / 10); putchar(48 ^ (x % 10));}

const int Mx = 1e6 + 1, Nx = 13970040;
int Mod;

int FastPow(ll a, int b){
    int res = 1;
    while(b){
        if(b & 1) res = res * a % Mod;
        b >>= 1, a = a * a % Mod;
    }
    return res;
}

bool Vis[Mx];
int Prime[78500], tot, Miu[Mx];
int Frac[Mx], Infs[Mx];
int St[Mx];
int MF[Nx], MG[Nx];
void TestifyInit(){
    Miu[1] = Frac[0] = Frac[1] = 1;
    for(int i = 2; i < Mx; ++i){
        if(!Vis[i]) Prime[++tot] = i, Miu[i] = -1;
        for(int j = 1; j <= tot && Prime[j] * i < Mx; ++j){
            Vis[i * Prime[j]] = 1;
            if(i % Prime[j] == 0) break;
            Miu[i * Prime[j]] = -Miu[i];
        }
        Frac[i] = 1ll * Frac[i - 1] * i % Mod;
    }
    St[1] = -1, Infs[Mx - 1] = FastPow(Frac[Mx - 1], Mod - 2);
    for(int i = 1; i < Mx - 1; ++i) St[i + 1] = St[i] + (Mx - 1) / i;
    for(int i = Mx - 1; i; --i) Infs[i - 1] = 1ll * Infs[i] * i % Mod;
    for(int d = 1; d < Mx; ++d){
        int lim = St[d] + (Mx - 1) / d;
        if(Miu[d] & 2){
            MF[St[d] + 1] = Infs[d], MG[St[d] + 1] = Infs[d - 1];
            for(int x = St[d] + 2, t = d; x <= lim; ++x){
                MF[x] = 1ll * Infs[t += d] * MF[x - 1] % Mod;
                MG[x] = 1ll * Infs[t - 1] * MG[x - 1] % Mod;
            }
            for(int x = St[d] + 2; x <= lim; ++x) MG[x] = 1ll * MG[x] * MG[x - 1] % Mod;
        }
        else if(Miu[d]){
            MF[St[d] + 1] = Frac[d], MG[St[d] + 1] = Frac[d - 1];
            for(int x = St[d] + 2, t = d; x <= lim; ++x){
                MF[x] = 1ll * Frac[t += d] * MF[x - 1] % Mod;
                MG[x] = 1ll * Frac[t - 1] * MG[x - 1] % Mod;
            }
            for(int x = St[d] + 2; x <= lim; ++x) MG[x] = 1ll * MG[x] * MG[x - 1] % Mod;
        }
        else for(int x = St[d] + 1; x <= lim; ++x) MF[x] = MG[x] = 1;
    }
    for(int d = 2; d < Mx; ++d){
        int lim = St[d] + (Mx - 1) / d;
        for(int x = St[d] + 1, y = St[d - 1] + 1; x <= lim; ++x, ++y)
            MF[x] = 1ll * MF[x] * MF[y] % Mod, MG[x] = 1ll * MG[x] * MG[y] % Mod;
    }
}

int main(){
    int T = Read(); Mod = Read();
    TestifyInit();
    while(T--){
        int n = Read(), m = Read(), k = Read();
        if(n > m) swap(n, m);
        int A = MG[n + m - 1];
        int B = 1ll * MG[n - 1] * MG[m - 1] % Mod * FastPow(MF[n - 1], m) % Mod * FastPow(MF[m - 1], n) % Mod; // Ans = A / B
        for(int l = 2, r; l <= n; l = r + 1){
            int N = n / l, M = m / l;
            r = min(n / N, m / M);
            A = 1ll * A * MG[St[r] + N + M] % Mod
                    * MG[St[l - 1] + N] % Mod
                    * MG[St[l - 1] + M] % Mod
                    * FastPow(MF[St[l - 1] + N], M) % Mod
                    * FastPow(MF[St[l - 1] + M], N) % Mod;
            B = 1ll * B * MG[St[l - 1] + N + M] % Mod
                    * MG[St[r] + N] % Mod
                    * MG[St[r] + M] % Mod
                    * FastPow(MF[St[r] + N], M) % Mod
                    * FastPow(MF[St[r] + M], N) % Mod;
        }
        printf("%d\n", FastPow(1ll * A * FastPow(B, Mod - 2) % Mod, k));
    }
    return 0;
}