Brute-force Algorithm(hdu3221)

Brute-force Algorithm

Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 2740 Accepted Submission(s): 728

Problem Description

Professor Brute is not good at algorithm design. Once he was asked to solve a path finding problem. He worked on it for several days and finally came up with the following algorithm:

Any fool but Brute knows that the function “funny” will be called too many times. Brute wants to investigate the number of times the function will be called, but he is too lazy to do it.

Now your task is to calculate how many times the function “funny” will be called, for the given a, b and n. Because the answer may be too large, you should output the answer module by P.

Input

There are multiple test cases. The first line of the input contains an integer T, meaning the number of the test cases.

For each test cases, there are four integers a, b, P and n in a single line.
You can assume that 1≤n≤1000000000, 1≤P≤1000000, 0≤a, b<1000000.

Output

For each test case, output the answer with case number in a single line.

Sample Input

3

3 4 10 3

4 5 13 5

3 2 19 100

Sample Output

Case #1: 2

Case #2: 11

Case #3: 12

思路：a的指数和b的指数是斐波数列，所以按照指数来，然后当指数小于p时直接快速幂，否则矩阵快速幂，加上oula降幂，前面小于p直接快速幂，也保证了后面使用欧拉降幂的条件。

  1  #include<stdio.h>
  2 #include<algorithm>
  3 #include<iostream>
  4 #include<string.h>
  5 #include<queue>
  6 #include<set>
  7 #include<math.h>
  8 using namespace std;
  9 typedef long long LL;
 10 typedef struct node
 11 {
 12         LL m[4][4];
 13         node()
 14         {
 15                 memset(m,0,sizeof(m));
 16         }
 17 } maxtr;
 18 int f1[1000];
 19 int f2[1000];
 20 maxtr E();
 21 void Init(maxtr *q);
 22 maxtr quick_m(maxtr ans,LL n,LL mod);
 23 LL quick(LL n,LL m,LL mod);
 24 bool prime[1000005];
 25 int aa[1000005];
 26 int oula[1000005];
 27 
 28 int main(void)
 29 {
 30         int T;
 31         int __ca = 0;
 32         int i,j;
 33         int k1,k2;
 34         f1[2] = 0;
 35         f1[3] = 1;
 36         f2[2] = 1;
 37         f2[3] = 1;
 38         for(i = 2; i < 1000; i++)
 39         {
 40                 for(j = i; (i*j) <= 1000000; j++)
 41                 {
 42                         prime[i*j] = true;
 43                 }
 44         }
 45         int cn = 0;
 46         for(i = 2; i <= 1000000; i++)
 47         {
 48                 oula[i] =i;
 49                 if(!prime[i])
 50                         aa[cn++] = i;
 51         }
 52         for(i = 0; i < cn; i++)
 53         {
 54                 for(j = 1; (j*aa[i])<=1000000; j++)
 55                 {
 56                         oula[j*aa[i]]/=aa[i];
 57                         oula[j*aa[i]]*=(aa[i]-1);
 58                 }
 59         }
 60         for(i = 4;; i++)
 61         {
 62                 f1[i] = f1[i-1]+f1[i-2];
 63                 if(f1[i] > 1000000)
 64                 {
 65                         k1 = i;
 66                         break;
 67                 }
 68         }
 69         for(i = 4;; i++)
 70         {
 71                 f2[i] = f2[i-1] + f2[i-2];
 72                 if(f2[i] > 1000000)
 73                 {
 74                         k2 = i;
 75                         break;
 76                 }
 77         }
 78         k1=max(k1,k2);
 79         scanf("%d",&T);
 80         while(T--)
 81         {
 82                 __ca++;
 83                 LL n,p,a,b;
 84                 scanf("%lld %lld %lld %lld",&a,&b,&p,&n);
 85                 printf("Case #%d: ",__ca);
 86                 if(n==1)
 87                 {
 88                         printf("%lld
",a%p);
 89                 }
 90                 else if(n == 2)
 91                 {
 92                         printf("%lld
",b%p);
 93                 }
 94                 else if(n == 3)
 95                 {
 96                         printf("%lld
",a*b%p);
 97                 }
 98                 else  if(n<=k1)
 99                 {
100                         int x1 = f1[n];
101                         int x2 = f2[n];
102                         LL ask = quick(a,x1,p);
103                         LL ack = quick(b,x2,p);
104                         printf("%lld
",ask*ack%p);
105                 }
106                 else
107                 {
108                         if(p==1)printf("0
");
109                         else
110                         {
111                                 maxtr cc ;
112                                 Init(&cc);
113                                 cc = quick_m(cc,n-3,oula[p]);
114                                 LL xx = cc.m[0][0] +cc.m[0][1];
115                                 LL yy = cc.m[1][0] + cc.m[1][1];
116                                 LL ask = quick(a,yy+oula[p],p)*quick(b,xx+oula[p],p)%p;
117                                 printf("%lld
",ask);
118                         }
119                 }
120         }
121         return 0;
122 }
123 maxtr E()
124 {
125         maxtr e;
126         for(int i = 0; i < 4; i++)
127         {
128                 for(int j = 0; j < 4; j++)
129                 {
130                         if(i == j)
131                                 e.m[i][j] = 1;
132                 }
133         }
134         return e;
135 }
136 void Init(maxtr *q)
137 {
138         q->m[0][0] = 1;
139         q->m[0][1] = 1;
140         q->m[1][0] = 1;
141         q->m[1][1] = 0;
142 }
143 maxtr quick_m(maxtr ans,LL n,LL mod)
144 {
145         int i,j,s;
146         maxtr ak = E();
147         while(n)
148         {
149                 if(n&1)
150                 {
151                         maxtr c;
152                         for(i = 0; i < 2; i++)
153                         {
154                                 for(j = 0; j < 2; j++)
155                                 {
156                                         for(s = 0; s < 2; s++)
157                                         {
158                                                 c.m[i][j] = c.m[i][j] + ans.m[i][s]*ak.m[s][j]%mod;
159                                                 c.m[i][j]%=mod;
160                                         }
161                                 }
162                         }
163                         ak = c;
164                 }
165                 maxtr d;
166                 for(i = 0; i < 2; i++)
167                 {
168                         for(j = 0; j < 2; j++)
169                         {
170                                 for(s= 0; s < 2; s++)
171                                 {
172                                         d.m[i][j] = d.m[i][j] + ans.m[i][s]*ans.m[s][j]%mod;
173                                         d.m[i][j]%=mod;
174                                 }
175                         }
176                 }
177                 ans = d;
178                 n>>=1;
179         }
180         return ak;
181 }
182 LL quick(LL n,LL m,LL mod)
183 {
184         LL ak = 1;
185         while(m)
186         {
187                 if(m&1)
188                         ak = ak*n%mod;
189                 n = n*n%mod;
190                 m>>=1;
191         }
192         return ak;
193 }

油！油！you@