思路与算法

类型1二进制快速幂。类型2先逐次除去gcd(a,P)，检查可除性并累计已消去的指数，再把剩余互质方程用BSGS找最小非负指数。类型3分解P为互质素数幂，对小素数幂求去p阶乘及阶乘中p的指数，再CRT合并；较大素模使用Lucas逐位组合数，乘法用精确Barrett约减。

正确性说明

快速幂保持已选位乘积。非互质约分保留原方程可解性及已确定幂次，互质后BSGS按指数商块递增，块内保存最小余数指数，首次匹配就是最小解。组合数在素数幂下由p指数与互质阶乘部分唯一决定，互质部分可逆；各素数幂模互质，CRT唯一恢复模P的值。Lucas只用于素模，未用于合数。

复杂度

类型1 O(logz)；类型2 O(sqrtP+logP)时间及O(sqrtP)空间。小素数幂组合各项预处理O(p^a)，递归阶乘和幂约O(log²z)，CRT及试除另计O(sqrtP)。大素模Lucas枚举各位较小侧，总O(Σmin(k_i,n_i-k_i)+logP)，只用常数辅助空间；初始化按实际表长计。

易错点

P=1全部余数为0，离散对数最小指数0。目标余数1优先指数0，不能把BSGS的第一个碰撞当任意指数。组合y>z输出0，模合数不能直接对阶乘套费马逆元；所有模乘控制在64位/128位已证明范围。

参考代码（C++11）

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll pw(ll a, ll b, ll p) {
    ll z = 1 % p;
    for (; b; b /= 2, a = a * a % p)
        if (b & 1)
            z = z * a % p;
    return z;
}
ll eg(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1;
        y = 0;
        return a;
    }
    ll u, v, g = eg(b, a % b, u, v);
    x = v;
    y = u - a / b * v;
    return g;
}
ll inv(ll a, ll p) {
    ll x, y;
    eg(a, p, x, y);
    return (x % p + p) % p;
}
ll gcd0(ll a, ll b) {
    while (b) {
        ll t = a % b;
        a = b;
        b = t;
    }
    return a;
}
ll logx(ll a, ll b, ll p) {
    a %= p;
    b %= p;
    if (p == 1 || b == 1)
        return 0;
    ll cnt = 0, t = 1;
    for (ll g; (g = gcd0(a, p)) > 1;) {
        if (b % g)
            return -1;
        b /= g;
        p /= g;
        t = t * (a / g) % p;
        cnt++;
        if (t == b)
            return cnt;
    }
    b = b * inv(t, p) % p;
    ll m = sqrt((long double)p) + 1, cur = 1;
    unordered_map<ll, int> baby;
    baby.reserve(2 * m);
    for (int j = 0; j < m; j++) {
        if (!baby.count(cur))
            baby[cur] = j;
        cur = cur * a % p;
    }
    ll shift = inv(cur, p);
    cur = b;
    for (int i = 0; i <= m; i++) {
        auto it = baby.find(cur);
        if (it != baby.end())
            return cnt + i * m + it->second;
        cur = cur * shift % p;
    }
    return -1;
}
struct Fast {
    unsigned long long p, mu;
    Fast(ll p) : p(p), mu(ULLONG_MAX / p + 1) {
    }
    ll mul(ll a, ll b) {
        unsigned long long x = (unsigned long long)a * b, q = ((__uint128_t)x * mu) >> 64,
                           t = q * p;
        return x >= t ? x - t : x + p - t;
    }
    ll power(ll a, ll n) {
        ll z = 1;
        for (; n; n /= 2, a = mul(a, a))
            if (n & 1)
                z = mul(z, a);
        return z;
    }
};
ll primeChoose(ll n, ll k, ll p) {
    Fast f(p);
    ll ans = 1;
    while (n || k) {
        ll a = n % p, b = k % p;
        if (b > a)
            return 0;
        b = min(b, a - b);
        ll num = 1, den = 1;
        for (ll i = 1; i <= b; i++) {
            num = f.mul(num, a - b + i);
            den = f.mul(den, i);
        }
        ans = f.mul(ans, f.mul(num, f.power(den, p - 2)));
        n /= p;
        k /= p;
    }
    return ans;
}
ll ppChoose(ll n, ll k, ll p, ll pk, int power) {
    vector<ll> fact(pk + 1, 1);
    for (ll i = 1; i <= pk; i++)
        fact[i] = fact[i - 1] * (i % p ? i : 1) % pk;
    function<ll(ll)> fac = [&](ll x) -> ll {
        if (!x)
            return 1;
        return pw(fact[pk], x / pk, pk) * fact[x % pk] % pk * fac(x / p) % pk;
    };
    auto exp = [&](ll x) {
        ll z = 0;
        while (x) {
            x /= p;
            z += x;
        }
        return z;
    };
    ll e = exp(n) - exp(k) - exp(n - k);
    if (e >= power)
        return 0;
    return fac(n) * inv(fac(k), pk) % pk * inv(fac(n - k), pk) % pk * pw(p, e, pk) % pk;
}
ll choose(ll n, ll k, ll P) {
    if (k > n || P == 1)
        return 0;
    ll rem = P, ans = 0;
    for (ll p = 2; p * p <= rem; p++)
        if (rem % p == 0) {
            ll pk = 1;
            int e = 0;
            while (rem % p == 0) {
                rem /= p;
                pk *= p;
                e++;
            }
            ll z = ppChoose(n, k, p, pk, e), M = P / pk;
            ans = (ans + z * M % P * inv(M % pk, pk)) % P;
        }
    if (rem > 1) {
        ll pk = rem, z = pk > 100000 ? primeChoose(n, k, pk) : ppChoose(n, k, pk, pk, 1),
           M = P / pk;
        ans = (ans + z * M % P * inv(M % pk, pk)) % P;
    }
    return ans;
}
int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int T;
    cin >> T;
    while (T--) {
        int type;
        ll y, z, P;
        cin >> type >> y >> z >> P;
        if (type == 1)
            cout << pw(y % P, z, P) << '\n';
        else if (type == 2) {
            ll x = logx(y, z, P);
            if (x < 0)
                cout << "Math Error\n";
            else
                cout << x << '\n';
        } else
            cout << choose(z, y, P) << '\n';
    }
}