Problem

https://www.lydsy.com/JudgeOnline/problem.php?id=3283

Solution

第一问快速幂。

第二问扩展 BSGS。使用递归进行子任务处理。记得特判 $z=1$ 的情况，否则 bsgs 会出无解。

第三问扩展 Lucas。使用玄妙的找规律方式快速处理阶乘，将因子单独剥离出来以方便求得逆元，最后再单独算因子的贡献，随后用 CRT 进行答案合并。

网上的代码写的都 emm 一言难尽

参考资料：
https://blog.csdn.net/clove_unique/article/details/54571216 （玄妙规律阶乘讲解）
https://www.cnblogs.com/DreamlessDreams/p/8711503.html （Lucas 好板子）
https://blog.csdn.net/moon1125666900/article/details/80555577 （EXBSGS 好板子核心）

PS：资料中 Lucas 板子可以不用事先求前缀积，如果你第二问效率不错的话。建议 Hash 挂链处理冲突，速度快。

下次我一定记得开 namespace

Code

// Code by ajcxsu
// Problem: caculator

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

typedef long long ll;

const int NN=1e6, M=1e6+10;
int h[NN], nexp[M], sp=1;
ll a[M], b[M];
inline void ins(int x, ll a, ll b) { nexp[sp]=h[x], h[x]=sp, ::a[sp]=a, ::b[sp]=b, sp++; }
void clr() { memset(h, 0, sizeof(h)), sp=1; } // Hash表

ll qpow(ll x, ll y, ll mo) {
    ll ret=1;
    while(y) {
        if(y&1ll) ret=ret*x%mo;
        x=x*x%mo, y>>=1ll;
    }
    return ret;
}

ll f(ll x, ll mo, ll p) {
    if(x<=1) return 1;
    ll ret=1;
    for(int i=2; i<=p; i++) if(i%mo) ret=ret*i%p;
    ret=qpow(ret, x/p, p);
    for(int i=2; i<=x%p; i++) if(i%mo) ret=ret*i%p;
    return ret*f(x/mo, mo, p)%p;
} // 阶乘规律玄学

void exgcd(ll a, ll b, ll &x, ll &y) {
    if(b==0) x=1, y=0;
    else exgcd(b, a%b, y, x), y-=a/b*x;
}

ll inv(ll x, ll mo) {
    ll nx, ny;
    exgcd(x, mo, nx, ny);
    return (nx%mo+mo)%mo;
}

ll exbsgs(ll y, ll z, ll p) {
    if(p==1 || z==1) return 0; // 特判
    ll d=__gcd(y, p);
    if(d>1) {
        if(z%d) return -0x3f3f3f3f;
        z/=d, p/=d, z=z*inv(y/d, p)%p;
        return exbsgs(y, z, p)+1; // exbsgs
    }
    int unit=ceil(sqrt(p)); // bsgs
    clr();
    ll yz=1;
    for(int i=0;i<unit;i++, yz=yz*y%p)
        ins(yz*z%p%NN, yz*z%p, i);
    ll nx=1;
    for(int i=0;i<=unit;i++, nx=nx*yz%p)
        for(int u=h[nx%NN];u;u=nexp[u])
            if(a[u]==nx && i*unit-b[u]>=0) return i*unit-b[u];
    return -0x3f3f3f3f;
}

const int N=30;
int pri[M], p;
char npri[M];
ll aa[N], mm[N];
ll exlucas(ll n, ll m, ll mo) {
    int t=0;
    ll rmo=mo;
    for(int i=0;i<p && mo>1;i++)
        if(mo%pri[i]==0) {
            mm[++t]=1;
            while(mo%pri[i]==0) mm[t]*=pri[i], mo/=pri[i];
            aa[t]=f(n, pri[i], mm[t])*inv(f(m, pri[i], mm[t]), mm[t])%mm[t]*inv(f(n-m, pri[i], mm[t]), mm[t])%mm[t];
            int k=0;
            for(int j=n/pri[i];j;j/=pri[i]) k+=j;
            for(int j=m/pri[i];j;j/=pri[i]) k-=j;
            for(int j=(n-m)/pri[i];j;j/=pri[i]) k-=j; // 单独处理因子
            aa[t]=aa[t]*qpow(pri[i], k, mm[t])%mm[t];
        }
    mo=rmo;
    ll ret=0;
    for(int i=1;i<=t;i++)
        (ret+=aa[i]*mo/mm[i]%mo*inv(mo/mm[i], mm[i]))%=mo; // CRT合并答案
    return ret;
}

int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    npri[1]=1;
    for(int i=2;i<M;i++) {
        if(!npri[i]) pri[p++]=i;
        for(int j=0; j<p && i*pri[j]<M; j++) {
            npri[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
    int n;
    cin>>n;
    int x;
    ll y, z, p, na;
    while(n--) {
        cin>>x>>y>>z>>p;
        if(x==1) cout<<qpow(y, z, p)<<endl;
        else if(x==2) {
            na=exbsgs(y, z, p);
            if(na<0) cout<<"Math Error"<<endl;
            else cout<<na<<endl;
        }
        else if(x==3) cout<<exlucas(z, y, p)<<endl;
    }
    return 0;
}