[ExBsgs]垃圾计算机

第一问 快速幂。
第二问 模数是质数的话用Bsgs，因为不是质数所以用ExBsgs。
第三问 模数是质数可以Lucas，不是质数所以就在crt的条件下乱搞，考虑一种特殊的求组合数的方法，可以发现在p[i]^c[i]的意义下，阶乘是有循环节的，就考虑拆成若干个循环节，然后那些p[i]次幂的东西就用分配律拆出来，就成为了p[i]的次幂*一个阶乘的形式，然后把拆出来的阶乘递归下去做，就可以在接近log的效率下完成在p[i]^c[i]的意义下求组合数。然后crt合并。
详细的看Bsgs的文章。

#include<cstdio>
#include<algorithm>
#include<string>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<map>
#include<iostream>
#define ll long long
#define Pair pair<int,int>
using namespace std;
const int N=1e6+7;
int n,y,z;
int c[N],p[N],m[N],H[N],M[N]; 
ll ans[N];
map<int,int> hash;
inline int read()
{
    char c;
    bool flag=false;
    while((c=getchar())>'9'||c<'0')
    if(c=='-')flag=true;
    int res=c-'0';
    while((c=getchar())>='0'&&c<='9')
    res=(res<<3)+(res<<1)+c-'0';
    return flag?-res:res;
}
inline int ksm(int s,int t)
{
    int res=1;
    while(t)
    {
        if(t&1) res=(ll)res*s;
        s=(ll)s*s;
        t>>=1;
    }
    return res;
}
inline int ksm(int s,int t,int p)
{
    int res=1;
    while(t)
    {
        if(t&1) res=(ll)res*s%p;
        s=(ll)s*s%p;
        t>>=1;
    }
    return res;
}
inline ll inv(ll a,int x)
{
    return ksm(a,(ll)(p[x]-1)*(m[x]/p[x])-1,m[x]);
}
inline Pair mul(Pair a,Pair b,int x)
{
    return Pair((ll)a.first*b.first%m[x],a.second+b.second);
}
inline Pair div(Pair a,Pair b,int x)
{
    return Pair((ll)a.first*inv(b.first,x)%m[x],a.second-b.second);
}
inline Pair trans(int a,int x)
{
    Pair res=Pair(0,0);
    while(a%p[x]==0&&a)
    {
        a/=p[x];
        res.second++;
    }
    res.first=a;
    return res;
}
Pair fac[N];
inline int get_num(Pair a,int x)
{
    return (ll)ksm(p[x],a.second,m[x])*a.first%m[x];
}
inline Pair fac_ksm(Pair s,int t,int x)
{
    Pair res=Pair(1,0);
    while(t)
    {
        if(t&1) res=mul(res,s,x);
        s=mul(s,s,x);
        t>>=1;
    }
    return res;
}
inline Pair get_fac(int n,int x)
{
    if(n==0) return Pair(1,0);
    Pair res=fac_ksm(fac[m[x]-1],n/m[x],x);
    res.second+=n/p[x];
    res=mul(res,get_fac(n/p[x],x),x);
//  for(int i=n%m[x];i>=1;--i)
//  if(i%p[x]!=0) res=mul(res,trans(i,x),x);
    res=mul(res,fac[n%m[x]],x);
    return res;
}
inline Pair C(int nn,int mm,int x)
{
//  return div(fac[n],mul(fac[m],fac[n-m],x),x);
    fac[0]=Pair(1,0);
    for(int i=1;i<m[x];++i)
    {
        if(i%p[x]!=0) fac[i]=mul(fac[i-1],trans(i,x),x);
        else fac[i]=fac[i-1];
    }
    return div(get_fac(nn,x),mul(get_fac(mm,x),get_fac(nn-mm,x),x),x);
}
inline Pair solve(int x)
{
//  fac[0]=Pair(1,0);
//  for(int i=1;i<=z;++i)
//  fac[i]=mul(fac[i-1],trans(i,x),x);
    ans[x]=get_num(C(z,y,x),x);
} 
int Fac[N],Inv[N];
inline int lucas(int n,int m,int P)
{
    if(n<m) return 0;
    if(n<P&&m<P) return (ll)Fac[n]*Inv[m]%P*Inv[n-m]%P;
    return (ll)lucas(n/P,m/P,P)*lucas(n%P,m%P,P)%P;
} 
int tot;
inline int bsgs(int y, int z, int p)
{
    hash.clear();
    int q=ceil(sqrt(p));
    int tmp=z;
    for(int j=1;j<=q;++j)
    {
        hash[tmp]=j;
        tmp=(ll)tmp*y%p;
    }
    int tm=ksm(y,q,p);
    tmp=1;
    for(int i=1;i<=q;++i)
    {
        tmp=(ll)tmp*tm%p;
        if(hash[tmp]) return i*q-hash[tmp];
    }
    return -1;
}
inline int gcd(int m,int n)
{
    return m%n?gcd(n,m%n):n;
}
int d;
inline int exgcd(int y,int &z,int &p,int &n)
{
    int m=gcd(y,p);
    if(m==1) return n;
    if(z%m!=0) return -1;
    z/=m; p/=m;
    d=(ll)d*(y/m)%p;
    ++n;
    if(z==d) return n;
    return exgcd(y,z,p,n);
}
inline int exbsgs(int y,int z,int P)
{
    int n=0,flag=0;
    d=1;
    if(z==1) return 0;
    flag=exgcd(y,z,P,n);
    if(flag==-1) return -1;
    int phi=1,e=sqrt(P),p=P;
    for(int i=2;i<=e;++i)
    if(p%i==0)
    {
        phi=(ll)phi*(i-1)%P;
        p/=i;
        while(p%i==0)
        {
            p/=i;
            phi=(ll)phi*i%P;
        }
    }
    if(p!=1) phi=(ll)phi*(p-1)%P;
    d=ksm(d,phi-1,P);
    return bsgs(y,(ll)z*d%P,P)+n;
}
int main()
{
    freopen("calculator.in","r",stdin);
    freopen("calculator.out","w",stdout);
    int P,type;
    ll tmp;
    n=read();
    while(n--)
    {
        type=read();
        y=read();
        z=read();
        P=read();
        if(type==1) printf("%d\n",ksm(y,z,P));
        if(type==2)
        {
            int ans = exbsgs(y, z, P);
            if(ans != -1) printf("%d\n", ans);
            else printf("Math Error\n");
        }
        if(type==3)
        {
            if(z<y)
            {
                printf("0\n");
                continue;
            }
            tot=0;tmp=P;
            int l=sqrt(P);
            for(int i=2;i<=l;++i)
            if(tmp%i==0)
            {
                p[++tot]=i;
                while(tmp%i==0)
                {
                    tmp/=i;
                    ++c[tot];
                }
            }
            if(tmp!=1)
            {
                p[++tot]=tmp;
                c[tot]=1;
            }
/*          if(tot==1&&c[1]==1)
            {
                Fac[0]=Inv[0]=1;
                for(int i=1;i<P;++i)
                Fac[i]=(ll)Fac[i-1]*i%P;
                Inv[P-1]=ksm(Fac[P-1],P-2,P);
                for(int i=P-2;i>=1;--i)
                Inv[i]=(ll)Inv[i+1]*(i+1)%P;
                printf("%d\n",lucas(z,y,P));
            }
            else
            {*/
                for(int i=1;i<=tot;++i)
                {
                    m[i]=ksm(p[i],c[i]);
                    M[i]=P/m[i];
                    H[i]=inv(M[i],i);
                    solve(i);
                } 
                ll Ans=0;
                for(int i=1;i<=tot;++i)
                    Ans=(Ans+(ll)ans[i]*M[i]%P*H[i]%P)%P;
                cout<<Ans<<endl;
//          }
            for(int i=1;i<=tot;++i) c[i]=0;
        }
    }
}