[HDOJ5667]Sequence（矩阵快速幂，费马小定理）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=5667

费马小定理：

假如p是质数，且gcd(a,p)=1，那么 a^(p-1)≡1（mod p）。

即：假如a是整数，p是质数，且a,p互质(即两者只有一个公约数1)，那么a的(p-1)次方除以p的余数恒等于1。

注意这里使用快速幂的时候要根据费马小定理对p-1取模。还有注意a%p=0的情况。

递推式：f(n)=f(n-1)*c+f(n-2)+1 非齐次。

构造矩阵：

|c 1 1|
|1 0 0|
|0 0 1|

初始的矩阵：

|1|
|0|
|1|

  1 #include <algorithm>
  2 #include <iostream>
  3 #include <iomanip>
  4 #include <cstring>
  5 #include <climits>
  6 #include <complex>
  7 #include <fstream>
  8 #include <cassert>
  9 #include <cstdio>
 10 #include <bitset>
 11 #include <vector>
 12 #include <deque>
 13 #include <queue>
 14 #include <stack>
 15 #include <ctime>
 16 #include <set>
 17 #include <map>
 18 #include <cmath>
 19
 20 using namespace std;
 21
 22 typedef long long ll;
 23 const ll maxn = 10;
 24 ll n, a, b, c, p;
 25
 26 typedef struct Matrix {
 27     ll m[maxn][maxn];
 28     ll r;
 29     ll c;
 30     Matrix(){
 31         r = c = 0;
 32         memset(m, 0, sizeof(m));
 33     }
 34 } Matrix;
 35
 36 Matrix mul(Matrix m1, Matrix m2, ll mod) {
 37     Matrix ans = Matrix();
 38     ans.r = m1.r;
 39     ans.c = m2.c;
 40     for(ll i = 1; i <= m1.r; i++) {
 41         for(ll j = 1; j <= m2.r; j++) {
 42                for(ll k = 1; k <= m2.c; k++) {
 43                 if(m2.m[j][k] == 0) continue;
 44                 ans.m[i][k] = ((ans.m[i][k] + m1.m[i][j] * m2.m[j][k] % mod) % mod) % mod;
 45             }
 46         }
 47     }
 48     return ans;
 49 }
 50
 51 Matrix quickmul(Matrix m, ll n, ll mod) {
 52     Matrix ans = Matrix();
 53     for(ll i = 1; i <= m.r; i++) {
 54         ans.m[i][i]  = 1;
 55     }
 56     ans.r = m.r;
 57     ans.c = m.c;
 58     while(n) {
 59         if(n & 1) {
 60             ans = mul(m, ans, mod);
 61         }
 62         m = mul(m, m, mod);
 63         n >>= 1;
 64     }
 65     return ans;
 66 }
 67
 68 ll qm(ll x, ll n, ll mod) {
 69     ll ans = 1, t = x;
 70     while(n) {
 71         if(n & 1) ans = (ans * t) % mod;
 72         t = (t * t) % mod;
 73         n >>= 1;
 74     }
 75     return ans;
 76 }
 77 int main() {
 78     // freopen("in", "r", stdin);
 79     int T;
 80     scanf("%d", &T);
 81     while(T--) {
 82         cin >> n >> a >> b >> c >> p;
 83         Matrix r;
 84         r.r = 3, r.c = 1;
 85         r.m[1][1] = 1;
 86         r.m[2][1] = 0;
 87         r.m[3][1] = 1;
 88         if(n == 1) {
 89             printf("1\n");
 90             continue;
 91         }
 92         if(n == 2) {
 93             printf("%I64d\n", qm(a, b, p));
 94             continue;
 95         }
 96         if(a % p == 0) {
 97             printf("0\n");
 98             continue;
 99         }
100         Matrix s;
101         s.r = s.c = 3;
102         s.m[1][1] = c, s.m[1][2] = 1, s.m[1][3] = 1;
103         s.m[2][1] = 1, s.m[2][2] = 0, s.m[2][3] = 0;
104         s.m[3][1] = 0, s.m[3][2] = 0, s.m[3][3] = 1;
105         s = quickmul(s, n-2, p-1);
106         ll ans = 0;
107         for(int i = 1; i <= r.r; i++) {
108             ans = (ans + (s.m[1][i] * r.m[i][1]) % (p - 1)) % (p - 1);
109         }
110         printf("%I64d\n", qm(a, (ans*b)%(p-1), p));
111     }
112     return 0;
113 }