学步园

题意：定义Fi为fibonacci数列的第i项，定义Ai(k)=Fi×i^k(i>=1)，求sigma(j=1 to n){Aj(k)}。

思路：终于把这道题A了，好开心，矩阵乘法什么的好神奇。。。首先推导一下Ai(k)，可以得到：

Ai(k) = Fi * i^k
         = ( Fi-1 + Fi-2 ) * i^k
         = Fi-1 * [ (i-1) + 1 ]^k + Fi-2 * [ (i-2) + 2 ]^k
         = sigma(j=0 to k){ C(k,j) * Fi-1 * (i-1)^j }  +  sigma(j=0 to k){ C(k,j) * Fi-2 * (i-2)^j * 2^(k-j)}
         = sigma(j=0 to k){Ai-1(j)*C(k,j)} + sigma(j=0 to k){Ai-1(j)*C(k,j)* 2^(k-j)}

可以看出，Ai(k)可以由Ai-1(k)和Ai-2(k)推出。因此我们可以根据这个来构造矩阵，矩阵中的元素分别为 sum(Ai(k))，Ai(0),Ai(1)……，Ai(k),Ai-1(0),Ai-1(1)……,Ai-1(k)。构造好矩阵后用快速幂就可以算出来了。

代码：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<string>
#include<algorithm>
#include<map>
#include<queue>
#include<set>
#include<stack>
#include<cmath>
#include<vector>
#define inf 0x3f3f3f3f
#define Inf 0x3FFFFFFFFFFFFFFFLL
#define eps 1e-9
#define pi acos(-1.0)
using namespace std;
typedef long long ll;
const int maxn=85;
const int mod=1000000007;
struct Matrix
{
    ll mat[maxn][maxn];
    int n;
    Matrix(){};
    Matrix(int nn){n=nn;}
    void Init()
    {
        for(int i=0;i<n;++i)
            for(int j=0;j<n;++j)
                mat[i][j]=(i==j);
    }
    void clear(){memset(mat,0,sizeof(mat));}
};
Matrix operator *(const Matrix& a,const Matrix & b)
{
    Matrix c(a.n);
    c.clear();
    for(int k=0;k<c.n;++k)
        for(int i=0;i<c.n;++i)
            for(int j=0;j<c.n;++j)
                c.mat[i][j]=(c.mat[i][j]+a.mat[i][k]*b.mat[k][j])%mod;
    return c;
}
ll C[55][55],p[55],num[maxn];
void Init()
{
    memset(C,0,sizeof(C));
    p[0]=1;
    for(int i=1;i<=50;++i)
        p[i]=(p[i-1]*2)%mod;
    for(int i=0;i<=50;++i)
        C[i][0]=1,C[i][1]=i;
    for(int i=1;i<=50;++i)
    {
        for(int j=2;j<=i;++j)
            C[i][j]=(C[i-1][j-1]+C[i-1][j])%mod;
    }
}
ll solve(ll n,int k)
{
    if(n==1) return 1;
    if(n==2) return (p[k+1]+1)%mod;
    ll sum=0;
    int nn=2*k+3;
    num[0]=p[k+1]+1;
    num[1]=2;num[k+2]=1;
    for(int i=1;i<=k;++i)
    {
        num[i+1]=2*p[i];
        num[i+k+2]=1;
    }
    Matrix x(nn),y(nn);
    x.Init(); y.clear();
    y.mat[0][0]=1;
    y.mat[1][1]=y.mat[k+2][1]=1;
    for(int i=0;i<=k;++i)
    {
        y.mat[i+1][0]=C[k][i];
        y.mat[i+k+2][0]=(C[k][i]*p[k-i])%mod;
        y.mat[i+1][i+k+2]=1;
    }
    for(int i=0;i<=k;++i)
    {
        for(int j=0;j<=i;++j)
        {
            y.mat[j+1][i+1]=C[i][j];
            y.mat[j+k+2][i+1]=(C[i][j]*p[i-j])%mod;
        }
    }
    n-=2;
    while(n)
    {
        if(n&1) x=x*y;
        y=y*y;
        n>>=1;
    }
    for(int i=0;i<nn;++i)
        sum=(sum+num[i]*x.mat[i][0])%mod;
    return sum;
}
int main()
{
    //freopen("in.txt","r",stdin);
    //freopen("out.txt","w",stdout);
    ll n;
    int k;
    Init();
    while(cin>>n>>k)
    {
        cout<<solve(n,k)<<endl;
    }
    return 0;
}