貼兩個連結吧反正我也說不清楚。

留個程式碼紀念一下。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<map>
using namespace std;
#define REP(I,ST,ED) for(int I=ST;I<=ED;++I)
#define DREP(I,ST,ED) for(int I=ST;I>=ED;--I)
const int maxn=1005,maxp=100005,maxm=10000005;
namespace Math{
	int ksm(int x,int y,int mod){
		int res=1;
		while(y){
			if(y&1)res=1ll*res*x%mod;
			x=1ll*x*x%mod,y>>=1;
		}return res;
	}
	int gcd(int a,int b){
		if(b==0)return a;
		else return gcd(b,a%b);
	}
	int Ex_gcd(int a,int b,int &x,int &y){
		if(b==0){
			x=1,y=0;
			return a;
		}else{
			int d=Ex_gcd(b,a%b,y,x);
			y-=a/b*x;
			return d;
		}
	}
	int inv(int x,int mod){
		int ans,tmp;
		Ex_gcd(x,mod,ans,tmp);
		return (ans+mod)%mod;
	}
	map<int,int>vis;
	int Ex_bsgs(int x,int y,int mod,int k){
		x%=mod,y%=mod;
		if(x==y)return k+1;
		int d=gcd(x,mod);
		if(d==1){
			if(!x)return y?-1:1;
			vis.clear();
			int m=ceil(sqrt(mod)),tmp=1;
			int inv_x=inv(ksm(x,m,mod),mod);
			y=1ll*y*inv(ksm(x,k,mod),mod)%mod;
			REP(i,0,m){
				if(!vis.count(tmp))vis[tmp]=i;
				tmp=1ll*tmp*x%mod;
			}REP(i,0,m){
				if(vis.count(y))return vis[y]+i*m+k;
				y=1ll*y*inv_x%mod;
			}return -1;
		}else if(y%d!=0)return y==1?0:-1;
		else return Ex_bsgs(x,y,mod/d,k+1);
	}
	int pi[35],pk[35],a[35],tot;
	int crt(int M){
		int ans=0;
		REP(i,1,tot){
			int mi=M/pk[i],x,y;
			Ex_gcd(mi,pk[i],x,y);
			(ans+=1ll*(x+M)%M*mi%M*a[i]%M)%=M;
		}return ans;
	}
	int calc(int n,int pi,int pk,int mod){
		if(!n)return 1;
		int tmp=1;
		REP(i,1,pk)
			if(i%pi)
				tmp=1ll*tmp*i%mod;
		tmp=ksm(tmp,n/pk,mod);
		REP(i,n/pk*pk+1,n)
			if(i%pi)
				tmp=1ll*tmp*i%mod;
		return 1ll*tmp*calc(n/pi,pi,pk,mod)%mod;
	}
	int multi_lucas(int x,int y,int pi,int pk){
		int tmp=0;
		for(int i=y;i;i/=pi)tmp+=i/pi;
		for(int i=x;i;i/=pi)tmp-=i/pi;
		for(int i=y-x;i;i/=pi)tmp-=i/pi;
		tmp=ksm(pi,tmp,pk);
		int a=calc(y,pi,pk,pk);
		int b=calc(x,pi,pk,pk);
		int c=calc(y-x,pi,pk,pk);
		return 1ll*tmp*a%pk*inv(b,pk)%pk*inv(c,pk)%pk;
	}
	int Ex_lucas(int x,int y,int p){
		if(x>y)return 0;
		int tmp=p;
		tot=0;
		for(int i=2;i*i<=tmp;i++)
			if(tmp%i==0){
				pi[++tot]=i;
				pk[tot]=1;
				while(tmp%i==0)
					tmp/=i,pk[tot]*=i;
			}
		if(tmp>1)
			pi[++tot]=pk[tot]=tmp;
		REP(i,1,tot){
			a[i]=multi_lucas(x,y,pi[i],pk[i]);
			(a[i]+=pk[i])%=pk[i];
		}return crt(p);
	}
	void solve(int type,int x,int y,int p){
		if(type==2){
			int ans=Ex_bsgs(x,y,p,0);
			if(ans==-1)puts("Math Error");
			else printf("%d\n",ans);
		}else printf("%d\n",type==1?ksm(x,y,p):Ex_lucas(x,y,p));
	}
}
int main(){
#ifndef ONLINE_JUDGE
	freopen("calculator.in","r",stdin);
	freopen("calculator.out","w",stdout);
#endif
	int T,type,x,y,p;
	scanf("%d",&T);
	while(T--){
		scanf("%d%d%d%d",&type,&x,&y,&p);
		Math::solve(type,x,y,p);
	}return 0;
}