[矩阵快速幂] CodeForces 392C Yet Another Number Sequence

题意：

题目意思很明朗~

思路：

A(n+1)=F(n+1)*(n+1)^k

A(n)=F(n)*(n)^k

A(n-1)=F(n-1)*(n-1)^k

这里拿k=2来举例

A(n+1)=F(n)*(n+1)^2+F(n-1)*(n+1)^2

对于A(n+1)发现可以由A(n)和A(n-1)得到

实际上就是多了 2*n+1个F(n) 和4*n个F(n-1)

其实就是n^k -> (n+1)^k  以及 (n-1)^k ->(n+1)^k

大家算一算就发现是和杨辉三角有关的，就是系数

这样我们就需要构造 n^(k-1)*F(n),n^(k-2)F(n), ... ,F(n)，以及 n^(k-1)*F(n-1),n^(k-2)F(n-1), ... ,F(n-1)

这样我们通过乘上系数，就能够完成转换了。

接着就是  (n+1)^(k-1)*F(n+1),(n+1)^(k-2)F(n+1), ... ,F(n+1) 和 (n+1)^(k-1)*F(n),(n+1)^(k-2)F(n), 

大家也可以算一算 都是可以通过之前的得到的。这样就能够构造矩阵转移了。

矩阵的大小可能比较大，所以需要细心一点。

代码：

#include"cstdlib"
#include"cstdio"
#include"cstring"
#include"cmath"
#include"queue"
#include"algorithm"
#include"iostream"
#include"map"
#define ll long long
#define mod 1000000007
using namespace std;
struct matrix
{
    ll mat[85][85];
    matrix()
    {
        memset(mat,0,sizeof(mat));
    }
};
matrix matadd(matrix a,matrix b,int n,int m)
{
    int i,j;
    matrix c;
    memset(c.mat,0,sizeof(c.mat));
    for(i=0; i<n; i++)  for(j=0; j<n; j++)  c.mat[i][j]=(a.mat[i][j]+b.mat[i][j])%m;
    return c;
}
matrix matmul(matrix a,matrix b,int n,int m)
{
    int i,j,k;
    matrix c;
    memset(c.mat,0,sizeof(c.mat));
    for(i=0; i<n; i++)
    {
        for(j=0; j<n; j++)
        {
            for(k=0; k<n; k++)
            {
                c.mat[i][j]+=a.mat[i][k]*b.mat[k][j];
                if(c.mat[i][j]>=m) c.mat[i][j]%=m;
            }
        }
    }
    return c;
}
matrix matpow(matrix a,ll k,int n,int m)
{
    matrix b;
    int i;
    memset(b.mat,0,sizeof(b.mat));
    for(i=0; i<n; i++) b.mat[i][i]=1;
    while(k)
    {
        if(k&1) b=matmul(a,b,n,m);
        a=matmul(a,a,n,m);
        k>>=1;
    }
    return b;
}
matrix matsum(matrix a,ll k,int n,int m)
{
    if(k<=0) return matpow(a,0,n,m);
    ll x=(k+1)/2;
    matrix b,ans;
    b=matsum(a,x-1,n,m);
    ans=matadd(b,matmul(b,matpow(a,x,n,m),n,m),n,m);
    if(k%2==0)  ans=matadd(ans,matpow(a,k,n,m),n,m);
    return ans;
}
ll bit[44];
ll c[44][44];

int main()
{
    ll n;
    int k;
    bit[0]=1;
    for(int i=1; i<42; i++) bit[i]=(bit[i-1]*2)%mod;
    memset(c,0,sizeof(c));
    for(int i=0; i<=41; i++)
    {
        c[i][0]=c[i][i]=1;
        for(int j=1; j<i; j++) c[i][j]=(c[i-1][j-1]+c[i-1][j])%mod;
    }
    while(scanf("%I64d%d",&n,&k)!=-1)
    {
        if(n==1)
        {
            puts("1");
            continue;
        }
        matrix a,b,ans;
        for(int i=0; i<k; i++)
        {
            a.mat[0][i]=bit[k-i];
            a.mat[0][i+k]=bit[k-1-i];
        }
        a.mat[0][2*k]=bit[k+1];
        a.mat[0][2*k+1]=1;
        a.mat[0][2*k+2]=(bit[k+1]+1)%mod;
        for(int i=0;i<k;i++)
        {
            for(int j=0;j<=i;j++)
            {
                b.mat[i][j]=c[k-1-j][i-j];
                b.mat[i+k][j]=b.mat[i][j+k]=c[k-1-j][i-j];
            }
        }
        for(int i=1;i<=k;i++)
        {
            b.mat[i-1][2*k]=b.mat[i-1][2*k+2]=c[k][i];
            if(i%2!=0) b.mat[i-1+k][2*k]=b.mat[i-1+k][2*k+2]=(c[k][i]*2)%mod;
        }
        b.mat[2*k][2*k]=b.mat[2*k+1][2*k]=b.mat[2*k][2*k+1]=1;
        b.mat[2*k][2*k+2]=b.mat[2*k+1][2*k+2]=b.mat[2*k+2][2*k+2]=1;

        ans=matmul(a,matpow(b,n-2,2*k+3,mod),2*k+3,mod);
        printf("%I64d\n",ans.mat[0][2*k+2]);
    }
    return 0;
}