HDU 5895 矩阵快速幂+高次幂取模

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
LL n,x,y,s1,S;
LL euler(LL n){
    LL res=n, a=n;
    for(LL i=2;i*i<=n;i++){
        if(a%i==0){
            res=res/i*(i-1);
            while(a%i==0) a/=i;
        }
    }
    if(a>1) res=res/a*(a-1);
    return res;
}
LL pow_mod(LL a,LL n){
    LL t=a, res=1;
    while(n){
        if(n&1) res=(res*t)%(S+1);
        n/=2;
        t=(t*t)%(S+1);
    }
    return res;
}
struct Mat{
    LL a[2][2];
    void init(){
        a[0][0]=2;
        a[1][0]=a[0][1]=1;
        a[1][1]=0;
    }
};
Mat operator *(Mat a,Mat b){
    Mat c;
    for(LL i=0;i<2;i++)
    for(LL j=0;j<2;j++){
        c.a[i][j]=0;
        for(LL k=0;k<2;k++)
            c.a[i][j]+=a.a[i][k]*b.a[k][j];
        c.a[i][j]=c.a[i][j]%s1;
    }
    return c;
}
Mat operator ^(Mat p,LL k){
    Mat ans; ans.init();
    while(k){
        if(k&1)
            ans=ans*p;
        k/=2;
        p=p*p;
    }
    return ans;
}
int main(){
    LL t; cin>>t;
    while(t--){
        cin>>n>>y>>x>>S;
        s1=2*euler(S+1);
        if(n==-1) break;
        if(n==0){
            cout<<1<<endl;
            continue;
        }
        Mat s; s.init();
        s=s^(n*y-1);
        LL u=(s.a[0][0]*s.a[1][0]);
        u=(u%s1+s1)/2;
        LL ans=pow_mod(x,u);
        cout<<ans<<endl;
    }

    return 0;
}