【HDOJ】3221 Brute-force Algorithm

归来吧很好推导。T(n) = a^f(n-1)*b^f(n)%p。
主要难点在于求mod和fibo。引用如下公式
A^B%C = A^(B%phi(C) + phi(C))%C, 满足B>=phi(C)。
  1 /*  */
  2 #include <iostream>
  3 #include <string>
  4 #include <map>
  5 #include <queue>
  6 #include <set>
  7 #include <stack>
  8 #include <vector>
  9 #include <deque>
 10 #include <algorithm>
 11 #include <cstdio>
 12 #include <cmath>
 13 #include <ctime>
 14 #include <cstring>
 15 #include <climits>
 16 #include <cctype>
 17 #include <cassert>
 18 using namespace std;
 19 
 20 #define rep(i,a,n)     for (int i=a;i<n;i++)
 21 #define per(i,a,n)     for (int i=n-1;i>=a;i--)
 22 #define pb             push_back
 23 #define mp             make_pair
 24 #define all(x)         (x).begin(),(x).end()
 25 #define SZ(x)         ((int)(x).size())
 26 #define lson        l, mid, rt<<1
 27 #define rson        mid+1, r, rt<<1|1
 28 
 29 const int maxn = 1000001;
 30 __int64 a, b, p, n;
 31 __int64 fibo[maxn];
 32 __int64 phi;
 33 
 34 typedef struct {
 35     __int64 m[2][2];
 36 } mat_t;
 37 
 38 mat_t res;
 39 
 40 mat_t matMult(mat_t a, mat_t &b) {
 41     mat_t c;
 42     
 43     memset(c.m, 0, sizeof(c.m));
 44     rep(i, 0, 2) {
 45         rep(j, 0, 2) {
 46             rep(k, 0, 2) {
 47                 c.m[i][j] += a.m[i][k] * b.m[k][j];
 48                 if (c.m[i][j] > phi) {
 49                     c.m[i][j] = c.m[i][j]%phi+phi;
 50                 }
 51             }
 52         }
 53     }
 54     return c;
 55 }
 56 
 57 __int64 euler(__int64 x) {
 58     __int64 ret = x, i;
 59     
 60     for (i=2; i*i<=x; ++i) {
 61         if (x%i == 0) {
 62             ret = ret / i * (i-1);
 63             while (x%i == 0)
 64                 x /= i;
 65         }
 66     }
 67     if (x > 1)
 68         ret = ret / x * (x-1);
 69     return ret;
 70 }
 71 
 72 __int64 myPow(__int64 base, __int64 n) {
 73     __int64 ret = 1;
 74     
 75     while (n) {
 76         if (n & 1)
 77             ret = ret * base %p;
 78         base = base * base %p;
 79         n >>= 1;
 80     }
 81     return ret;
 82 }
 83 
 84 void matFibo(__int64 n) {
 85     mat_t base;
 86     
 87     res.m[0][0] = base.m[0][0] = 0;
 88     res.m[0][1] = res.m[1][0] = res.m[1][1] =
 89     base.m[0][1] = base.m[1][0] = base.m[1][1] = 1;
 90     
 91     while (n) {
 92         if (n & 1)
 93             res = matMult(base, res);
 94         base = matMult(base, base);
 95         n >>= 1;
 96     }
 97 }
 98 
 99 void getFibo(__int64 &an, __int64 &bn) {
100     int i;
101     
102     fibo[1] = 0;
103     fibo[2] = 1;
104     for (i=3; i<=n; ++i) {
105         fibo[i] = fibo[i-1] + fibo[i-2];
106         if (fibo[i] > phi)
107             break;
108     }
109     if (i > n) {
110         an = fibo[n-1];
111         bn = fibo[n];
112         return ;
113     }
114     matFibo(n-3);
115     an = res.m[1][0];
116     bn = res.m[1][1];
117 }
118 
119 int main() {
120     int i, j, k;
121     int t;
122     __int64 an, bn, ans;
123     
124     #ifndef ONLINE_JUDGE
125         freopen("data.in", "r", stdin);
126         freopen("data.out", "w", stdout);
127     #endif
128     
129     scanf("%d", &t);
130     rep(tt, 1, t+1) {
131         scanf("%I64d %I64d %I64d %I64d", &a, &b, &p, &n);
132         printf("Case #%d: ", tt);
133         if (a==0 || b==0 || p==1) {
134             printf("0
");
135         } else if (n == 1) {
136             printf("%I64d
", a%p);
137         } else if (n == 2) {
138             printf("%I64d
", b%p);
139         } else {
140             phi = euler(p);
141             getFibo(an, bn);
142             ans = myPow(a, an)*myPow(b, bn)%p;
143             printf("%I64d
", ans);
144         }
145     }
146     
147     #ifndef ONLINE_JUDGE
148         printf("time = %d.
", (int)clock());
149     #endif
150     
151     return 0;
152 }