（矩阵快速幂）hdu5950 Recursive sequence

今天的专题是矩阵快速幂，其实就是把快速幂算法应用到矩阵中去，把普通的乘法重载成矩阵相乘，从而解决一些实际问题。基本操作如下：

快速幂取模：应用的是取模运算对乘法的可结合性，以及二进制的原理——一个数可以被拆分成若干个2的幂相加之和。具体代码如下：

int qpow(int a, int b, int p) {
    int res = 1;
    while(b) {//循环到b为0时，所有的因子都已相乘，算法结束
        if(b & 1) {//相当于b%2 == 1
            res = (res * a) % c;//b为奇数，需补乘一个a
        }
        b >>= 1;//相当于b /= 2
        a = a*a % c;//相当于int k = a*a
    }
    return res;
}

应用到矩阵中，则定义一个矩阵乘法运算，并将power运算稍作变化。代码如下：

struct mat {
    ll a[maxn][maxn];
    mat() {//构造方法
        memset(a, 0, sizeof(a));
    }
};

mat init() {//初始化
    mat tt;
    for(int i = 0;i < n;i++) {
        for(int j = 0;j < n;j++) {
            if(i == j)
                tt.a[i][j] = 1;
            else
                tt.a[i][j] = 0;
        }
    }
    return tt;
}

mat mul(mat x, mat y) {//乘法
    mat tt;
    for(int i = 0;i < n;i++) {
        for(int j = 0;j < n;j++) {
            tt.a[i][j] = 0;
            for(int k = 0;k < n;k++) {
                tt.a[i][j] += x.a[i][k] * y.a[k][j];
            }
            tt.a[i][j] %= mod;
        }
    }
    return tt;
}

mat mul(mat x, mat y) {//乘法
    mat tt;
    for(int i = 0;i < n;i++) {
        for(int j = 0;j < n;j++) {
            tt.a[i][j] = 0;
            for(int k = 0;k < n;k++) {
                tt.a[i][j] += x.a[i][k] * y.a[k][j];
            }
            tt.a[i][j] %= mod;
        }
    }
    return tt;
}

而本题就是一个应用的实例。题目要求这样一个以递归定义的表达式第n项的值（将结果mod2147493647）：

F(n) = F(n-1) + 2*F(n-2) + n^4.

由于数据范围很大，直接用循环算或用递归算的话必然会爆long long，用大数运算的话时空效率均很低。故利用线性代数中基的思想，采取矩阵相乘的形式将其拆开运算。

具体过程......这个辣鸡博客不让贴图片好气啊...!

因此仅通过对系数矩阵A做n-2次乘幂运算，再与一个已知的向量相乘，得到的向量第一个元素即为F（n）。具体ac代码如下：

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long ll;

const ll mod = 2147493647;
const ll maxn = 2147483647;

int cas;
ll n, a, b;

struct mat {
    ll a[7][7];
};

void print(mat ans) {//打印一个矩阵，一直wa，调试时方便
    for(int i = 0;i < 7;i++) {
        for(int j = 0;j < 7;j++)
            printf("%d ", ans.a[i][j]);
        printf("\n");
    }
    return;
}

mat init0() {//初始化为0
    mat tt;
    for(int i = 0;i < 7;i++)
        for(int j = 0;j < 7;j++) tt.a[i][j] = 0;
    return tt;
}

mat init1() {//初始化为单位矩阵
    mat tt;
    for(int i = 0;i < 7;i++)
        for(int j = 0;j < 7;j++) tt.a[i][j] = (i==j);
    return tt;
}

mat mul(mat x, mat y) {//矩阵 乘法
    mat tt;
    for(int i = 0;i < 7;i++) {
        for(int j = 0;j < 7;j++) {
            tt.a[i][j] = 0;
            for(int k = 0;k < 7;k++) {
                tt.a[i][j] += x.a[i][k] * y.a[k][j];
            }
            tt.a[i][j] %= mod;
        }
    }
    return tt;
}

mat power(mat tt, ll b) {//快速幂取模
    mat res = init1();
    while(b) {
        if(b & 1) {
            res = mul(res, tt);
        }
        b >>= 1;
        tt = mul(tt, tt);
    }
    return res;
}

int main() {
    while(~scanf("%d", &cas)) {
        mat coef = init0();
        coef.a[0][0] = 1;
        coef.a[0][1] = 2;
        coef.a[0][2] = 1;
        coef.a[0][3] = 4;
        coef.a[0][4] = 6;
        coef.a[0][5] = 4;
        coef.a[0][6] = 1;
        coef.a[1][0] = 1;
        coef.a[2][2] = 1;
        coef.a[2][3] = 4;
        coef.a[2][4] = 6;
        coef.a[2][5] = 4;
        coef.a[2][6] = 1;
        coef.a[3][3] = 1;
        coef.a[3][4] = 3;
        coef.a[3][5] = 3;
        coef.a[3][6] = 1;
        coef.a[4][4] = 1;
        coef.a[4][5] = 2;
        coef.a[4][6] = 1;
        coef.a[5][5] = 1;
        coef.a[5][6] = 1;
        coef.a[6][6] = 1;
        //print(coef);
        while(cas--) {
            scanf("%I64d %I64d %I64d", &n, &a, &b);
            if(n == 1)
                printf("%I64d\n", a);
            else if(n == 2)
                printf("%I64d\n", b);
            else {
                mat aa = power(coef, n-2);
                mat bb = {b, 0, 0, 0, 0, 0, 0,
                          a, 0, 0, 0, 0, 0, 0,
                          16,0, 0, 0, 0, 0, 0,
                          8, 0, 0, 0, 0, 0, 0,
                          4, 0, 0, 0, 0, 0, 0,
                          2, 0, 0, 0, 0, 0, 0,
                          1, 0, 0, 0, 0, 0, 0};
                /*ll ans = 0;
                ans = (ans + aa.a[0][0]*b) % mod;
                ans = (ans + aa.a[0][1]*a) % mod;
                ans = (ans + aa.a[0][2]*16)% mod;
                ans = (ans + aa.a[0][3]*8) % mod;
                ans = (ans + aa.a[0][4]*4) % mod;
                ans = (ans + aa.a[0][5]*2) % mod;
                ans = (ans + aa.a[0][6]*1) % mod;*/
                mat cc = mul(aa, bb);
                //printf("d\n", ans);
                printf("%lld\n", cc.a[0][0]);
            }
        }
    }
    return 0;
}