hdoj 5950 Recursive sequence

题目链接：Recursive sequence

题目大意：有一个递推公式$f(i) = 2*f(i-2)+f(i-1)+i^4,f(1) = a,f(2) = b$,现在给出a,b,和n，求f(n)，a,b,n < 2^31

题目思路：很容易想到矩阵快速幂，然后推一个公式就好：
原矩阵：

$$\begin{Bmatrix} b & a & 16 & 8 & 4 & 2 & 1 \end{Bmatrix}$$ 递推矩阵： $$\begin{Bmatrix} 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 2 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 1 & 0 & 0 & 0 & 0\\4 & 0 & 4 & 1 & 0 & 0 & 0\\6 & 0 & 6 & 3 & 1 & 0 & 0\\4 & 0 & 4 & 3 & 2 & 1 & 0\\1 & 0 & 1 & 1 & 1 & 1 & 1 \end{Bmatrix}$$

实际上就是用
$\begin{Bmatrix} f(i) & f(i-1) & i^4 & i^3 & i^2 & i & 1 \end{Bmatrix}$
去推出
$\begin{Bmatrix} f(i+1) & f(i) & (i+1)^4 & (i+1)^3 & (i+1)^2 & (i+1) & 1 \end{Bmatrix}$这个罢了，学过线性代数的很轻松能做出来这个罢，然后就直接矩阵快速幂就好了

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;
const int maxn =  7;
const ll MOD = 2147493647;
#define mod(x) ((x)%MOD)

ll n;

struct mat{
    ll m[maxn][maxn];
}unit;

mat operator * (mat a,mat b){
    mat ret;
    ll x;
    for(ll i = 0;i < 7;i++){
        for(ll j = 0;j < 7;j++){
            x = 0;
            for(ll k = 0;k < 7;k++)
                x += mod((ll)a.m[i][k]*b.m[k][j]);
            ret.m[i][j] = mod(x);
        }
    }
    return ret;
}

void init_unit(){
    for(ll i = 0;i < maxn;i++)
        unit.m[i][i] = 1;
    return ;
}

mat pow_mat(mat a,ll n){
    mat ret = unit;
    while(n){
        if(n&1) ret = ret*a;
        a = a*a;
        n >>= 1;
    }
    return ret;
}

int main(){
    ll n,t,A,B;
    init_unit();
    scanf("%lld",&t);
    while(t--){
        scanf("%lld%lld%lld",&n,&A,&B);
        if(n == 1) printf("%lld\n",A);
        else if(n == 2) printf("%lld\n",B);
        else{
            mat a,b;
            b.m[0][0] = 1,b.m[0][1] = 1,b.m[0][2] = 0,b.m[0][3] = 0,b.m[0][4] = 0,b.m[0][5] = 0,b.m[0][6] = 0;
            b.m[1][0] = 2,b.m[1][1] = 0,b.m[1][2] = 0,b.m[1][3] = 0,b.m[1][4] = 0,b.m[1][5] = 0,b.m[1][6] = 0;
            b.m[2][0] = 1,b.m[2][1] = 0,b.m[2][2] = 1,b.m[2][3] = 0,b.m[2][4] = 0,b.m[2][5] = 0,b.m[2][6] = 0;
            b.m[3][0] = 4,b.m[3][1] = 0,b.m[3][2] = 4,b.m[3][3] = 1,b.m[3][4] = 0,b.m[3][5] = 0,b.m[3][6] = 0;
            b.m[4][0] = 6,b.m[4][1] = 0,b.m[4][2] = 6,b.m[4][3] = 3,b.m[4][4] = 1,b.m[4][5] = 0,b.m[4][6] = 0;
            b.m[5][0] = 4,b.m[5][1] = 0,b.m[5][2] = 4,b.m[5][3] = 3,b.m[5][4] = 2,b.m[5][5] = 1,b.m[5][6] = 0;
            b.m[6][0] = 1,b.m[6][1] = 0,b.m[6][2] = 1,b.m[6][3] = 1,b.m[6][4] = 1,b.m[6][5] = 1,b.m[6][6] = 1;

            a.m[0][0] = B,a.m[0][1] = A,a.m[0][2] = 16,a.m[0][3] = 8,a.m[0][4] = 4,a.m[0][5] = 2,a.m[0][6] = 1;
            b = pow_mat(b,n-2);
            a = a*b;
            printf("%lld\n",a.m[0][0]%MOD);
        }
    }
    return 0;
}