2016ACM-ICPC亚洲区沈阳站 C-Recursive sequence 矩阵快速幂

递推式为F[n]=2*F[n-2]+F[n-1]+n^4.用矩阵快速幂求解:

构造 { F[n-2],F[n-1],n^4,n^3,n^2,n^1,1 } * X =
{ F[n-1],F[n],(n+1)^4,(n+1)^3,(n+1)^2,(n+1)^1,1 }.

可求得矩阵 X =
{ 0 2 0 0 0 0 0 }
{ 1 1 0 0 0 0 0 }
{ 0 1 1 0 0 0 0 }
{ 0 0 4 1 0 0 0 }
{ 0 0 6 3 1 0 0 }
{ 0 0 4 3 2 1 0 }
{ 0 0 1 1 1 1 1 }

由题意可知F[1]=A,F[2]=B,所以我们递推式的n从3开始.
特殊判断一下当n=1或者n=2时直接输出解,n>=3时构造一个 I = { A,B,3^4,3^3,3^2,3^1,1} .
求Ans = I * X^(n-3),答案即为2*Ans[0][0]+Ans[0][1]+Ans[0][2].

#include<bits/stdc++.h>
using namespace std;
const int INF=0x3f3f3f3f;
const long long MOD=2147493647;
const int SIZE=16;
typedef long long ll;
struct Matrix
{
    int n;
    ll Mat[SIZE][SIZE];
    Matrix(int a):n(a)
    {
        for (int i=0; i<n; i++)
            for (int j=0; j<n; j++)
                Mat[i][j]=0;
    }
    Matrix operator * (Matrix& x)
    {
        Matrix result(n);
        for (int k=0; k<n; k++)
            for (int i=0; i<n; i++)
                for (int j=0; j<n; j++)
                    result.Mat[i][j]=(result.Mat[i][j]+Mat[i][k]%MOD*x.Mat[k][j]%MOD)%MOD;
        return result;
    }
    Matrix operator ^ (ll x)
    {
        Matrix temp(n),a(*this);
        for (int i=0; i<n; i++)
            temp.Mat[i][i]=1;
        while (x)
        {
            if (x&1)
                temp=a*temp;
            a=a*a;
            x>>=1;
        }
        return temp;
    }
    ll* operator [] (const int& a)
    {
        return Mat[a];
    }
};
ll T,N,A,B;
int main()
{
    scanf("%lld",&T);
    Matrix X(7);
    X[0][1]=2;
    X[1][0]=1,X[1][1]=1;
    X[2][1]=1,X[2][2]=1;
    X[3][2]=4,X[3][3]=1;
    X[4][2]=6,X[4][3]=3,X[4][4]=1;
    X[5][2]=4,X[5][3]=3,X[5][4]=2,X[5][5]=1;
    X[6][2]=1,X[6][3]=1,X[6][4]=1,X[6][5]=1,X[6][6]=1;
    while (T--)
    {
        scanf("%lld%lld%lld",&N,&A,&B);
        Matrix I(7),Ans(7);
        I[0][0]=A,I[0][1]=B,I[0][2]=81,I[0][3]=27,I[0][4]=9,I[0][5]=3,I[0][6]=1;
        if (N>=3)
        {
            Ans=X^(N-3);
            Ans=I*Ans;
            cout<<(2*Ans[0][0]+Ans[0][1]+Ans[0][2])%MOD<<endl;
        }
        else if (N==2)
            cout<<B<<endl;
        else if (N==1)
            cout<<A<<endl;
    }
    return 0;
}