这个题目，是个数论题，里面涉及到了很多数论的性质
首先，根据费马小定理，a^(p-1) = 1 (mod p) 其中p为质数
由于对于p有这样的性质 so当指数大于p-1的时候 就可以将指数对p-1取余即可~~这样就可以将指数控制在p-1的范围内。。。。
其次对于前面的a^(p-1)/2 和b^(p-1)/2...使用二分快速幂即可
另外有一个很特别的结论（由uncle luo提出） 对于a^(p-1)/2 (mod p) 一定会是-1或+1 ...
对于后面那部分就是重头戏了，首先求Fn mod p-1...其次根号里面那玩意，是一个数列，（首先既然能够mod p，说明这数列是整数列啦）
其次，这是个通过n就能求出通项的通项公式，这是在高中用递推式最渴望达到的结果 （= =） ，可是呢，你要知道有根号，所以直接代入是一个很可怕的事。。。
实际上两项一起的时候达到相除根号的目的，而这个时候应该想办法通项公式求递推式。。。相当于以前求特征方程然后求通项的逆运算 - -
然后对于z=c1*x^n+c2*y^n ，你要构造以x，y为解的一元二次方程...然后把这个特征方程(x^2=px+q)转成递归式(Zn=p*Zn-1+q*Zn-2)就可以了
在这题，我们令Zn=x^2n+y^2n 。。通过代入n=0,1可以得到初始值（这跟求特征方程然后求通项的时候通过初始值得到c1,c2同一个道理） = =
最后注意负数取余的问题就行了~~
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include<iostream>
6 using namespace std;
7 #define MSIZE 2
8 #define bigt __int64
9 bigt mod;
10 bigt mod2;
11 typedef struct matr{
12 bigt a[MSIZE][MSIZE];
13 int top;
14 }mat;
15 mat m1,m2,unit,m3;
16 bigt powermod(bigt a,bigt b,bigt c) //a^b %c
17 {
18 bigt d=1,e;
19 a%=c;
20 e=a;
21 while(b)
22 {
23 if(b&1) d=(d*e)%c;
24 b>>=1;
25 e=(e*e)%c;
26 }
27 return d;
28 }
29
30 bigt modgod(bigt kk,bigt modd)
31 {
32 while(kk<0)
33 {
34 kk+=modd*1000000;
35 }
36 kk%=modd;
37 return kk;
38 }
39
40 void init(int SIZE)
41 {
42 int i,j;
43 unit.top=m1.top=m2.top=m3.top=SIZE;
44 for(i=0;i<SIZE;i++)
45 {
46 for(j=0;j<SIZE;j++)
47 {
48 unit.a[i][j]=(i==j);
49 m1.a[i][j]=1;
50 m3.a[i][j]=0;
51 }
52 }
53 m1.a[0][0]=0;
54 }
55
56 mat matadd(mat p,mat q)
57 {
58 int i,j;
59 for (i=0;i<p.top;i++)
60 {
61 for (j=0;j<p.top;j++)
62 {
63 p.a[i][j]+=q.a[i][j];
64 p.a[i][j]%=mod;
65 }
66 }
67 return p;
68 }
69
70 mat matmul(mat p,mat q,bigt modd)
71 {
72 mat c;
73
74 int i,j,k;
75 c.top=p.top;
76 for (i=0;i<c.top;i++)
77 {
78 for (j=0;j<c.top;j++)
79 {
80 c.a[i][j] = 0;
81 for (k=0;k<c.top;k++)
82 {
83 c.a[i][j]+=(p.a[i][k]*q.a[k][j]);
84 c.a[i][j]=modgod(c.a[i][j],modd);
85 }
86 }
87 }
88 return c;
89 }
90
91 mat matmop(int t,mat m,bigt modd)
92 {
93 mat mm=unit;
94 if(t>0)
95 while(1)
96 {
97 if(t&1)
98 {
99 mm=matmul(mm,m,modd);
100 }
101 t>>=1;
102 if(!t) break;
103 m=matmul(m,m,modd);
104 }
105 return mm;
106 }
107
108
109 int main()
110 {
111 int n,i,j,k,t,m;
112 bigt x,y;
113 mat fib;
114 mat ans;
115 bigt zz;
116 bigt pp,qq;
117 init(2);
118 scanf("%d",&t);
119 while(t--)
120 {
121 cin>>x>>y>>n>>mod2;
122 mod=mod2-1;
123 fib=matmop(n,m1,mod);
124 pp=powermod(x,mod2/2,mod2);
125 pp%=mod2;
126 qq=powermod(y,mod2/2,mod2);
127 qq%=mod2;
128 if(pp==-1 || qq==-1) {puts("0");continue;}
129 pp++;qq++;
130 pp=(pp*qq)%mod2;
131 m2.a[0][0]=0;m2.a[0][1]=1;
132 m2.a[1][0]=modgod(-(x-y)*(x-y),mod2);
133 m2.a[1][1]=(2*(x+y))%mod2;
134 m3.a[0][0]=2;
135 m3.a[1][0]=(2*(x+y))%mod2;
136 fib=matmop(fib.a[1][1]-1,m2,mod2);
137 ans=matmul(fib,m3,mod2);
138 zz=ans.a[1][0]%mod2;
139 zz=(zz*pp)%mod2;
140 cout<<zz<<endl;
141 }
142 return 0;
143 }
