HDOJ 5950 Recursive sequence (矩阵快速幂)

http://acm.hdu.edu.cn/showproblem.php?pid=5950
题意：给定一个数列的前两项a和b，数列的递推公式是$f_i=2f_{i-2}+f_{i-1}+i^4$
求第N个数，其中N,a,b<2^31。

类似矩阵快速幂求斐波那契数列，这里比较难搞的是i^4。
考虑i^4变成(i+1)^4时多出的几项：
$i^4 \to(i+1)^4=i^4+4i^3+6i^2+4i+1$
所以可以构建矩阵B如下：

$$\begin{bmatrix} 1&1&1&1&1&0&0\\ 0&1&2&3&4&0&0\\ 0&0&1&3&6&0&0\\ 0&0&0&1&4&0&0\\ 0&0&0&0&1&0&1\\ 0&0&0&0&0&0&2\\ 0&0&0&0&0&1&1\\ \end{bmatrix}$$
每一列分别对应 $1\ i\ i^2\ i^3\ i^4\ f_{i-1}\ f_{i-2}$
这是一个7*7的矩阵，得到i=3时，相应的A矩阵应构造成一个1*7的矩阵如下：
$$\begin{bmatrix} 1&3&9&27&81&a&b\\ \end{bmatrix}$$
这样A乘B^(n-2)得到的矩阵C的c[1][7]即为所求。

#include<iostream>
#include<cstdio>
#include<cstring>
#define MOD 2147493647
using namespace std;

const int maxn = 10;
int T,A,B,n;

//矩阵结构体
struct Matrix{
    long long a[maxn][maxn];
    void init(){    //初始化为单位矩阵
        memset(a, 0, sizeof(a));
        for(int i=0;i<maxn;++i){
            a[i][i] = 1;
        }
    }
};

//矩阵乘法
Matrix mul(Matrix a, Matrix b){
    Matrix ans;
    for(int i=0;i<maxn;++i){
        for(int j=0;j<maxn;++j){
            ans.a[i][j] = 0;
            for(int k=0;k<maxn;++k){
                ans.a[i][j]+=a.a[i][k]*b.a[k][j];
                ans.a[i][j]%=MOD;
            }
        }
    }
    return ans;
}

//矩阵快速幂
Matrix qpow(Matrix a, long long n){
    Matrix ans;
    ans.init();
    while(n){
        if(n&1) ans = mul(ans, a);
        a = mul(a, a);
        n /= 2;
    }
    return ans;
}

void output(Matrix a){
    printf("%lld\n",a.a[0][6]%MOD);
}

int main()
{
    scanf("%d",&T);
    while(T--)
    {
        scanf("%d%d%d",&n,&A,&B);
        Matrix a;
        a.init();
        int p=1;
        for(int i=0;i<=4;i++)
        {
            a.a[0][i]=p;
            p*=3;
        }
        a.a[0][5] = A;
        a.a[0][6] = B;
        Matrix b;
        b.init();
        for(int i=0;i<7;i++)
            for(int j=0;j<7;j++)
                b.a[i][j]=0;
        b.a[0][0]=b.a[0][1]=b.a[0][2]=b.a[0][3]=b.a[0][4]=1;
        b.a[1][1]=1;b.a[1][2]=2;b.a[1][3]=3;b.a[1][4]=4;
        b.a[2][2]=1;b.a[2][3]=3;b.a[2][4]=6;
        b.a[3][3]=1;b.a[3][4]=4;
        b.a[4][4]=1;b.a[4][6]=1;
        b.a[5][6]=2;
        b.a[6][5]=1;b.a[6][6]=1;
        if(n==1)
        {
            printf("%d\n",A);
        }
        else if(n==2)
        {
            printf("%d\n",B);
        }
        else
        {
            Matrix ans = mul(a, qpow(b, n-2));
            output(ans);
        }
    }
    return 0;
}