Codeforces 392C Yet Another Number Sequence (矩阵快速幂+二项式展开)

题意：已知斐波那契数列fib(i) , 给你n 和 k , 求∑fib(i)*i^k (1<=i<=n)

思路：不得不说，这道题很有意思，首先我们根据以往得出的一个经验，当我们遇到 X^k 的形式，当 X 很大，k很小时，我们可以利用二项式定理进行展开，然后求出递推式在利用矩阵加速

推导过程：

已知 fib(1) = 1, fib(2) = 1,fib(i) = fib(i-1) + fib(i-2);

       Ai(k) =fib(i)*i^k;

根据数学归纳法，我们可知 fib(i+1)*(i+1)^k = ( fib(i) + fib(i-1) ) * (i+1)^k 必然成立

将上式进行二项式展开 

上式 = ( fib(i) + fib(i-1) ) * {C[k][0]*i^0 + C[k][1]*i^1 +...+C[k][k]*i^k}; 

       = C[k][0]*（i^0* fib(i) + i^0*fib(i-1) ) + C[k][1]*(i^1* fib(i) + i^1*fib(i-1) )....（乘法结合律

我们定义 g1[i][a] = fib(i)*i^a , g2[i][a] = fib(i-1)*i^a

将定义的g1 , g2 带入上式中得到递推公式

g1[i][k] = C[k][0]*(g1[i-1][0]+g2[i-1][0])+C[k][1]*(g1[i-1][1]+g2[i-1][1])+...+C[k][k](g1[i-1][k]+g2[i-1][k]);

g2[i][k]= C[k][0]*g1[i-1][0] +C[k][1]*g1[i-1][1]+..+C[k][k]*g1[i-1][k];

最后再根据递推式写出构造矩阵即可

代码：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#define ll long long
using namespace std;
ll N,M,P;
const ll MOD=1000000007;
ll C[85][85];
struct  Matrix
{
    ll m[85][85];
};

Matrix  A;

Matrix  I;

void get_C()
{
    memset(C,0,sizeof(C));
    C[0][0]=1;
    for(int i=1;i<=82;i++)
    {
        C[0][i]=C[i][i]=1;
        for(int j=1;j<=i;j++)
        {
            C[j][i]=(C[j][i-1]+C[j-1][i-1])%MOD;
        }
    }
}

void make_I(int k)
{
    memset(I.m,0,sizeof(I.m));
    for(int i=0;i<=2*(k+1);i++)
      for(int j=0;j<=2*(k+1)+1;j++)
       if(i==j)
       I.m[i][j]=1;
}

void make_A(int k)
{
    memset(A.m,0,sizeof(A.m));
    for(int i=0;i<=k;i++)
    {
        for(int j=0;j<=i;j++)
        {
            A.m[i][j]=C[j][i];
        }
    }
    for(int i=0;i<=k;i++)
    {
        for(int j=0;j<=i;j++)
        {
            A.m[i+k+1][j]=C[j][i];
        }
    }
    for(int i=0;i<=k;i++)
    {
        for(int j=k+1;j<=2*(k+1)-1;j++)
        {
            A.m[i][j]=C[j-k-1][i];
        }
    }
    A.m[2*(k+1)][k]=1;
    A.m[2*(k+1)][2*(k+1)]=1;
}

Matrix multi(Matrix a,Matrix b)
{
    Matrix ans;
    for(int i=0;i<=N;i++)
    {
       for(int j=0;j<=M;j++)
       {
          ans.m[i][j]=0;
          for(int k=0;k<=P;k++)
          {
             ans.m[i][j]=(ans.m[i][j]+a.m[i][k]*b.m[k][j])%MOD;
          }
          ans.m[i][j]%=MOD;
       }
    }
    return ans;
}

Matrix power(Matrix a,ll k)
{
    Matrix ans=I,p=a;
    while(k)
    {
      if(k&1)
      {
        ans=multi(ans,p);
      }
      k>>=1;
      p=multi(p,p);
    }
    return ans;
}

int main()
{
    ll n;
    int k;
    get_C();
    while(cin>>n>>k)
    {
        if(n==1)
        {
            cout<<"1"<<endl;
            continue;
        }
        else
        {
            make_A(k);
            make_I(k);
            N=M=P=2*(k+1);
            Matrix ans = power(A,n);
            ll num=0;
            for(int i=0;i<M;i++)
            num=(num+ans.m[2*(k+1)][i])%MOD;
            cout<<num<<endl;
        }
    }
    return 0;
}