前两块可以看成是不是二次剩余，快速幂计算即可。

后半部分可以看成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
     
       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
      
       >=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\n",tmp1%p*tmp2%p*calc(a,b,n,p)%p);71     }72     return 0;73 }