前两块可以看成是不是二次剩余，快速幂计算即可。
后半部分可以看成x1=a+b+2ab,x2=a+b-2ab为特征方程x^2-px-qx=0的两根
然后可以通过韦达定理求出p和q，因此递推式为A(n+2)=pA(n+1)+qA(n)
还要用费马小定理化简一下斐波那契数。
矩阵快速幂即可求。
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 ll mod; 5 struct Max 6 { 7 ll a[4][4];int n,m; 8 Max(){n=m=0;memset(a,0,sizeof(a));}; 9 Max operator *(const Max &b)const{ 10 Max c;c.n=n;c.m=b.m; 11 for(int k=0;k<m;++k) 12 for(int i=0;i<n;++i) 13 for(int j=0;j<b.m;++j) 14 c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%mod)%mod; 15 return c; 16 } 17 }; 18 ll qmod(ll n,ll m,ll p) 19 { 20 ll ans=1; 21 while(m) 22 { 23 if(m&1)ans=ans*n%p; 24 n=n*n%p;m>>=1; 25 } 26 return ans; 27 } 28 Max Qmod(Max a,ll b) 29 { 30 Max ans;ans.n=a.n;ans.m=a.m; 31 ans.a[0][0]=ans.a[1][1]=1; 32 while(b) 33 { 34 if(b&1)ans=ans*a; 35 a=a*a;b>>=1; 36 } 37 return ans; 38 } 39 ll calc(ll a,ll b,ll n,ll p) 40 { 41 mod=p-1; 42 Max f;f.n=f.m=2; 43 f.a[0][0]=1;f.a[0][1]=1; 44 f.a[1][0]=1;f.a[1][1]=0; 45 f=Qmod(f,n); 46 mod=p; 47 int k=f.a[0][0]; 48 Max z;z.n=z.m=2; 49 z.a[0][0]=2*(a+b)%mod;z.a[1][0]=-(a-b)*(a-b)%mod; 50 z.a[0][1]=1;z.a[1][1]=0; 51 z=Qmod(z,k-1); 52 Max ans;ans.n=1,ans.m=2; 53 ans.a[0][1]=2,ans.a[0][0]=2*(a+b)%mod; 54 ans=ans*z; 55 return (ans.a[0][0]+mod)%mod; 56 } 57 int main() 58 { 59 int T;ll a,b,n,p; 60 scanf("%d",&T); 61 while(T--) 62 { 63 scanf("%lld%lld%lld%lld",&a,&b,&n,&p); 64 int tmp1=qmod(a,(p-1)/2,p); 65 int tmp2=qmod(b,(p-1)/2,p); 66 if(tmp1==-1||tmp2==-1) 67 { 68 puts("0");continue; 69 }tmp1++;tmp2++; 70 printf("%lld ",tmp1%p*tmp2%p*calc(a,b,n,p)%p); 71 } 72 return 0; 73 }