hdu 3221 Bruteforce Algorithm

http://acm.hdu.edu.cn/showproblem.php?pid=3221

　　矩阵快速幂加上x^n取模公式的一道题。

　　题意是计算函数在题目图片中的暴力程序中被调用了多少次。只要稍微模拟一下就可以知道是求a^f(n)+b^f(n+1) (mod P)是多少，其中f(n)为a[1]=1、a[2]=0且a[i]=a[i-1]+a[i-2]中的a[n]。这题有一条公式可以直接利用，x^n (mod m) = x^(n%phi(m)+phi(m)) (mod m)，其中n>m，phi(m)为m的欧拉函数。这种方法的复杂度为O(sqrt(m))。

　　其实这题也可以用之前一题的cycle detection的技术来求出循环节的长度，这样做的复杂度是O(m)。

直接套公式的代码如下：

View Code

  1 #include <cstdio>
  2 #include <iostream>
  3 #include <cstring>
  4 #include <algorithm>
  5 #include <cstdlib>
  6 
  7 using namespace std;
  8 
  9 typedef long long ll;
 10 const int maxSize = 2;
 11 const int initMod = 1E9 + 7;
 12 int curMod = initMod;
 13 
 14 struct Mat {
 15     ll v[maxSize][maxSize];
 16 
 17     Mat (ll ONE = 0) {
 18         for (int i = 0; i < maxSize; i++) {
 19             for (int j = 0; j < maxSize; j++) {
 20                 v[i][j] = 0;
 21             }
 22             v[i][i] = ONE;
 23         }
 24     }
 25 } Base, op;
 26 
 27 Mat operator * (Mat &_a, Mat &_b) {
 28     Mat ret = Mat ();
 29 
 30     for (int i = 0; i < maxSize; i++) {
 31         for (int k = 0; k < maxSize; k++) {
 32             if (_a.v[i][k]) {
 33                 for (int j = 0; j < maxSize; j++) {
 34                     ret.v[i][j] += _a.v[i][k] * _b.v[k][j];
 35                     ret.v[i][j] %= curMod;
 36                 }
 37             }
 38         }
 39     }
 40 
 41     return ret;
 42 }
 43 
 44 Mat operator ^ (Mat &__a, int p) {
 45     Mat ret = Mat (1);
 46     Mat _a = __a;
 47 
 48     while (p) {
 49         if (p & 1) {
 50             ret = ret * _a;
 51         }
 52         _a = _a * _a;
 53         p >>= 1;
 54     }
 55 
 56     return ret;
 57 }
 58 
 59 const int maxn = 1000005;
 60 int prm[maxn / 5], prmCnt;
 61 bool notPrm[maxn];
 62 
 63 void getPrime() {
 64     notPrm[0] = notPrm[1] = true;
 65     prmCnt = 0;
 66     for (int i = 2; i < maxn; i++) {
 67         if (!notPrm[i]) prm[prmCnt++] = i;
 68         for (int j = 0; j < prmCnt && prm[j] * i < maxn; j++) {
 69             notPrm[prm[j] * i] = true;
 70             if (i % prm[j] == 0) break;
 71         }
 72     }
 73 }
 74 
 75 int phi(int n) {
 76     int ret = 1;
 77 
 78     for (int i = 0; i < prmCnt; i++) {
 79         if (n <= 1) break;
 80         if (!notPrm[n]) {
 81             ret *= n - 1;
 82             break;
 83         }
 84         if (n % prm[i]) continue;
 85         n /= prm[i];
 86         ret *= prm[i] - 1;
 87         while (n % prm[i] == 0) {
 88             n /= prm[i];
 89             ret *= prm[i];
 90         }
 91     }
 92 
 93     return ret;
 94 }
 95 
 96 ll powerMod(ll a, int p, int m) {
 97     ll ret = 1;
 98 
 99     while (p) {
100         if (p & 1) {
101             ret *= a;
102             ret %= m;
103         }
104         a *= a;
105         a %= m;
106         p >>= 1;
107     }
108 
109     return ret;
110 }
111 
112 void makeMat() {
113     Base = Mat();
114     Base.v[0][0] = 1;
115     op.v[0][0] = 0;
116     op.v[1][0] = op.v[0][1] = op.v[1][1] = 1;
117 }
118 
119 int deal(int a, int b, int P, int n) {
120     if (n <= 30) {
121         ll f1 = -1, f2 = 1;
122 
123         while (n--) {
124             f2 += f1;
125             f1 = f2 - f1;
126         }
127         return powerMod(a, f1, P) * powerMod(b, f2, P) % P;
128     }
129     curMod = phi(P);
130 
131     makeMat();
132     op = op ^ (n - 1);
133     Base = Base * op;
134 //    cout << Base.v[0][0] << ' ' << Base.v[0][1] << ' ' << curMod << endl;
135 
136     return powerMod(a, Base.v[0][0] + curMod, P) * powerMod(b, Base.v[0][1] + curMod, P) % P;
137 }
138 
139 int main() {
140     int a, b, P, n;
141     int T;
142 
143     freopen("in", "r", stdin);
144     getPrime();
145     cin >> T;
146     for (int c = 1; c <= T; c++) {
147         cin >> a >> b >> P >> n;
148         cout << "Case #" << c << ": " << deal(a, b, P, n) << endl;
149     }
150 
151     return 0;
152 }

——written by Lyon