参考 https://blog.csdn.net/queuelovestack/article/details/52577212

题目链接

https://acm.hdu.edu.cn/showproblem.php?pid=5895

思路

用f(n-1)乘上f(n)=f(n-2)+2*f(n-1),再通过移项、累加后得 g[n]=f[n] * f[n+1]/2

那么就可以首先通过矩阵快速幂计算出g(n*y)的值

关于除以2的处理：
因为2不一定和p互质，不一定能用费马小定理求出2的逆元

要求 $ \frac{a}{b} \ mod \ p $,不妨先设其值为d,
=> \(\frac{a}{b}=k*p+d\)
=> \(a=k*b*p+b*d\)
=> \(\frac{a \ mod (b*p)}{b}=d\)
⭐因此模数变成了 b*p !!!⭐

因为 \(x^{g(n*y)}\) 肯定大于 $\phi(s+1) $，所以这个地方可以用欧拉定理来降幂

最后再用快速幂求出最终答案

code

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N = 2;
int t;
ll y,n,x,s;
struct Matrix{
	long long m[2][2];
};

Matrix mul(Matrix a,Matrix b,ll mod){
	Matrix c;
	for(int i=0;i<N;i++){  //NxN的矩阵 
		for(int j=0;j<N;j++){
			c.m[i][j]=0;
		}
	}
	for(int i=0;i<2;i++){
		for(int j=0;j<2;j++){
			for(int k=0;k<2;k++){
				c.m[i][j]+=a.m[i][k]*b.m[k][j];
				c.m[i][j]%=mod;
			}
		}
	}
	return c;
}

Matrix ksm(Matrix a,ll t,ll mod){   //t可能要开longlong 
	Matrix ans;
	for(int i=0;i<N;i++){  //NxN的矩阵 
		for(int j=0;j<N;j++){
			if(i==j) ans.m[i][j]=1; 
			else ans.m[i][j]=0;
		}
	}
	while(t){
		if(t&1){
			ans=mul(ans,a,mod);
		}
		a=mul(a,a,mod);
		t/=2;
	}
	return ans;
}

ll qpow(ll a,ll b,ll mod){
	ll ans=1;
    a%=mod;
	while(b){
		if(b%2){
			ans*=a;
			ans%=mod; 
		}
		a*=a;
		a%=mod;  //注意这里也需要取模 
		b/=2;
	}
	return ans;
}

ll Euler(ll n){   
	ll nn=n;
	ll ans=n;
	for(ll i=2;i*i<=nn;i++){
		if(n%i==0){
			ans=ans-ans/i;
			while(n%i==0){
				n/=i;
			}
		}
	}
	if(n>1){
		ans=ans-ans/n;
	}
	return ans;
} 

ll f(ll id,ll mod){
    Matrix init;
    for(int i=0;i<2;i++){
        for(int j=0;j<2;j++){
            if(i==0&&j==0) init.m[i][j]=1;
            else init.m[i][j]=0;
        }
    }
    Matrix a;
    a.m[0][0]=2; a.m[0][1]=1;
    a.m[1][0]=1; a.m[1][1]=0;
    Matrix ans=mul(init,ksm(a,id-1,mod),mod);
    return ans.m[0][0]*ans.m[0][1];
}

int main(){
    cin>>t;
    while(t--){
        scanf("%lld%lld%lld%lld",&n,&y,&x,&s);
        s+=1;
        ll X=1ll*n*y;
        ll p=2*Euler(s);
        X=f(X+1,p)%p/2+p/2;
        printf("%lld\n",qpow(x,X,s));
    }
    system("pause");
    return 0;
}