hdu 1588 矩阵运算的应用

这题弄了两天才做出来，还是去请教了竹教主。

贴个别人的解说吧，自己懒得写了

把斐波那契数列转化为矩阵:
A={1,1}
    {1,0};

{f[n+1],f[n]}
{f[n],f[n-1]} = A^n ;最后输出右上角那项
或者用
{f[n+2],f[n+1]}
{f[n+1], f[n] } = A^(n+1); 最后输出右下角那项

我们用第一个公式
所求即为A^b + A^(k+b) + A^(2*k+b) + ... + A^((n-1)*k+b)
=A^b * ( A^0 + A^k + A^(2*k) + ... + A^((n-1)*k) )
=A^b * ( (A^k)^0 + (A^k)^1 + (A^k)^2 + ...+ (A^k)^(n-1) );
B=A^k;
上式=A^b * ( B^0 + B^1 + B^2 + ... + B^(n-1) );
B^0 + B^1 + B^2 + ... + B^(n-1)用上篇介绍到的递归二分 方法求解
最后输出矩阵的第二项(右上角)即可；

对于求解 B^0 + B^1 + B^2 + ... + B^(n-1)
我们也可以构造矩阵的方法。

我们来设置这样一个矩阵
B I
O I
其中O是零矩阵，I是单位矩阵

将它乘方，得到
B^2 I+B
O I
乘三方，得到
B^3 I+B+B^2
O I
乘四方，得到
B^4 I+B+B^2+B^3
O I
用快速幂求出n方，让第二项再与A^b相乘即可。

/*
 * hdu1588/win.cpp
 * Created on: 2011-11-30
 * Author    : ben
 */
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include <iostream>
#include <algorithm>
#include <queue>
#include <set>
#include <map>
#include <stack>
#include <string>
#include <vector>
#include <deque>
#include <list>
#include <functional>
#include <numeric>
#include <cctype>
using namespace std;
typedef long long typec;
int MOD;
// 矩阵的最大阶
const int MAX_ORDER = 4;
struct MyMatrix {
    int order;
    typec num[MAX_ORDER][MAX_ORDER];
    MyMatrix(int ord) {
        order = ord;
    }
    void init() {
        for (int i = 0; i < order; i++) {
            for (int j = 0; j < order; j++) {
                num[i][j] = 0;
            }
        }
    }
};
MyMatrix operator*(MyMatrix ma, MyMatrix mb) {
    int ord = ma.order;
    MyMatrix numc(ord);
    numc.init();
    int i, j, k;
    for (i = 0; i < ord; i++) {
        for (j = 0; j < ord; j++) {
            for (k = 0; k < ord; k++) {
                numc.num[i][j] += ma.num[i][k] * mb.num[k][j];
                numc.num[i][j] %= MOD;
            }
        }
    }
    return numc;
}
MyMatrix mpow(MyMatrix ma, typec x) {
    int ord = ma.order;
    MyMatrix numc(ord);
    numc.init();
    for (int i = 0; i < ord; i++) {
        numc.num[i][i] = 1;
    }
    for (; x; x >>= 1) {
        if (x & 1) {
            numc = numc * ma;
        }
        ma = ma * ma;
    }
    return numc;
}

inline void initA(MyMatrix &m) {
    m.num[0][0] = 1;
    m.num[0][1] = 1;
    m.num[1][0] = 1;
    m.num[1][1] = 0;
}

inline void initC(MyMatrix &C, MyMatrix B) {
    C.init();
    for(int i = 0; i < 2; i++) {
        for(int j = 0; j < 2; j++) {
            C.num[i][j] = B.num[i][j];
        }
    }
    C.num[0][2] = 1;
    C.num[1][3] = 1;
    C.num[2][2] = 1;
    C.num[3][3] = 1;
}

MyMatrix getBrFromD(MyMatrix &D) {
    MyMatrix res(2);
    for(int i = 0; i < 2; i++) {
        for(int j = 0; j < 2; j++) {
            res.num[i][j] = D.num[i][j + 2];
        }
    }
    return res;
}

int main() {
#ifndef ONLINE_JUDGE
    freopen("data.in", "r", stdin);
#endif
    int k, b, n;
    while(scanf("%d%d%d%d", &k, &b, &n, &MOD) == 4 ){
        MyMatrix A(2);
        initA(A);
        MyMatrix B = mpow(A, k);
        MyMatrix C(4);
        initC(C, B);
        MyMatrix D = mpow(C, n);
        MyMatrix E = getBrFromD(D);
        MyMatrix ans = mpow(A, b) * E;
        printf("%d\n", (int)ans.num[1][0]);
    }
    return 0;
}