[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|

#include <algorithm>
#include <iostream>
#include <iomanip>
#include <cstring>
#include <climits>
#include <complex>
#include <fstream>
#include <cassert>
#include <cstdio>
#include <bitset>
#include <vector>
#include <deque>
#include <queue>
#include <stack>
#include <ctime>
#include <set>
#include <map>
#include <cmath>

using namespace std;

typedef long long ll;
const ll maxn = 10;
ll n, a, b, c, p;

typedef struct Matrix {
    ll m[maxn][maxn];
    ll r;
    ll c;
    Matrix(){
        r = c = 0;
        memset(m, 0, sizeof(m));
    } 
} Matrix;

Matrix mul(Matrix m1, Matrix m2, ll mod) {
    Matrix ans = Matrix();
    ans.r = m1.r;
    ans.c = m2.c;
    for(ll i = 1; i <= m1.r; i++) {
        for(ll j = 1; j <= m2.r; j++) {
               for(ll k = 1; k <= m2.c; k++) {
                if(m2.m[j][k] == 0) continue;
                ans.m[i][k] = ((ans.m[i][k] + m1.m[i][j] * m2.m[j][k] % mod) % mod) % mod;
            }
        }
    }
    return ans;
}

Matrix quickmul(Matrix m, ll n, ll mod) {
    Matrix ans = Matrix();
    for(ll i = 1; i <= m.r; i++) {
        ans.m[i][i]  = 1;
    }
    ans.r = m.r;
    ans.c = m.c;
    while(n) {
        if(n & 1) {
            ans = mul(m, ans, mod);
        }
        m = mul(m, m, mod);
        n >>= 1;
    }
    return ans;
}

ll qm(ll x, ll n, ll mod) {
    ll ans = 1, t = x;
    while(n) {
        if(n & 1) ans = (ans * t) % mod;
        t = (t * t) % mod; 
        n >>= 1;
    }
    return ans;
}
int main() {
    // freopen("in", "r", stdin);
    int T;
    scanf("%d", &T);
    while(T--) {
        cin >> n >> a >> b >> c >> p;
        Matrix r;
        r.r = 3, r.c = 1;
        r.m[1][1] = 1;
        r.m[2][1] = 0;
        r.m[3][1] = 1;
        if(n == 1) {
            printf("1
");
            continue;
        }
        if(n == 2) {
            printf("%I64d
", qm(a, b, p));
            continue;
        }
        if(a % p == 0) {
            printf("0
");
            continue;
        }
        Matrix s;
        s.r = s.c = 3;
        s.m[1][1] = c, s.m[1][2] = 1, s.m[1][3] = 1;
        s.m[2][1] = 1, s.m[2][2] = 0, s.m[2][3] = 0;
        s.m[3][1] = 0, s.m[3][2] = 0, s.m[3][3] = 1;
        s = quickmul(s, n-2, p-1);
        ll ans = 0;
        for(int i = 1; i <= r.r; i++) {
            ans = (ans + (s.m[1][i] * r.m[i][1]) % (p - 1)) % (p - 1);
        }
        printf("%I64d
", qm(a, (ans*b)%(p-1), p));
    }
    return 0;
}