HDU5950-Recursive sequence(矩阵快速幂）

题目链接：Recursive sequence

题意：给出n头母牛，第一头报a，第二头报b，第i头报f[i-2]*2+f[i-1]+i^4,问第n头母牛报数多少

题目大意：

　　Fi=Fi-1+2Fi-2+i4。给定F1和F2求Fn。

题目思路：

　　【递推+矩阵快速幂】

　　现场用算了1个多小时的公式过了。

　　主要还是我太菜。递推写的太少。

　　先考虑f(i)=f(i-1)+2f(i-2)，很容易写出递推矩阵

　　　　0 2

　　　　1 1

　　(i+1)4=i4+4i3+6i2+4i+1。

　　所以需要在递推矩阵种存下i的4 3 2 1 0次幂，以便推出(i+1)4，矩阵为

　　　　1 0 0 0 0

　　　　4 1 0 0 0

　　　　6 3 1 0 0

　　　　4 3 2 1 0

　　　　1 1 1 1 1

　　于是f={fi-1,fi,i4,i3,i2,i1,i0}，将以上两个矩阵合并，即可推出{fi,fi+1,(i+1)4,(i+1)3,(i+1)2,(i+1)1,(i+1)0}.矩阵如下

　　　　0 2 0 0 0 0 0

　　　　1 1 0 0 0 0 0

　　　　0 1 1 0 0 0 0

　　　　0 4 4 1 0 0 0

　　　　0 6 6 3 1 0 0

　　　　0 4 4 3 2 1 0

　　　　0 1 1 1 1 1 1

　　推出转移矩阵后只需要根据n求矩阵快速幂即可。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;
const long long mod = 2147493647;
struct prog
{
    long long a[8][8];
};
prog s,B;
prog matrixmul(prog a,prog b)
{
    prog c;
    for(int i=1;i<8;++i)for(int j=1;j<8;++j)
    {
        c.a[i][j]=0;
        for(int k=1;k<8;k++)
            c.a[i][j]+=(a.a[i][k]*b.a[k][j])%mod;
        c.a[i][j]%=mod;
    }
    return c;
}
prog mul(prog s,int k)
{
    prog ans;
    for(int i=1;i<8;++i)for(int j=1;j<8;++j) ans.a[i][j]=(i==j)?1:0;
    while(k){
        if(k&1)
            ans=matrixmul(ans,s);
        k>>=1;
        s=matrixmul(s,s);
    }
    return ans;
}
int main()
{
    int n,t,a,b;
    for(scanf("%d",&t);t--;){
        scanf("%d %d %d",&n,&a,&b);
        if(n==1){printf("%lld\n",a%mod);continue;}
        if(n==2){printf("%lld\n",b%mod);continue;}
        if(n==3){printf("%lld\n",(81+2*a%mod+b%mod)%mod);continue;} 
        n-=2;
        for(int i=1;i<=7;++i)for(int j=1;j<=7;++j) s.a[i][j]=0,B.a[i][j]=0;
        for(int i=1; i<=5; i++)s.a[i][1]=1;
        for(int i=2; i<=5; i++)s.a[i][2]=i-1;
        s.a[3][3]=1;s.a[4][3]=3;s.a[5][3]=6;
        s.a[4][4]=1;s.a[5][4]=4;
        s.a[5][5]=1;s.a[6][5]=1;
        s.a[6][6]=1;s.a[7][6]=1;
        s.a[6][7]=2;
        B.a[1][1]=1;B.a[2][1]=3;B.a[3][1]=9;B.a[4][1]=27;B.a[5][1]=81;B.a[6][1]=b;B.a[7][1]=a;
        s=mul(s,n);
        s=matrixmul(s,B);
        printf("%lld\n",s.a[6][1]%mod);
    }
    return 0;
}