Yet Another Number Sequence——[矩阵快速幂]

Description

　　Everyone knows what the Fibonacci sequence is. This sequence can be defined by the recurrence relation:

F₁ = 1, F₂ = 2, F_i = F_i - 1 + F_i - 2 (i > 2).

We'll define a new number sequence A_i(k) by the formula:

A_i(k) = F_i × i^k (i ≥ 1).

In this problem, your task is to calculate the following sum: A₁(k) + A₂(k) + ... + A_n(k). The answer can be very large, so print it modulo 1000000007(10⁹ + 7).

Input

　　The first line contains two space-separated integers n, k(1 ≤ n ≤ 10¹⁷; 1 ≤ k ≤ 40).

Output

　　Print a single integer — the sum of the first n elements of the sequence A_i(k) modulo 1000000007(10⁹ + 7).

Sample Input

Input

1 1

Output

1

Input

4 1

Output

34

Input

5 2

Output

316

Input

7 4

Output

73825

解题思路：
　　这是一道矩阵快速幂求多项式的应用，思路很清晰：找到各个量之间的递推关系，用矩阵求幂转移状态。问题的关键在于推导A（n，k），sumA（n），
F（n）之间的关系。过程如下：
　　1.令 U（n＋1，k）＝F（n＋1）＊（n＋1）^k；
　　　　　V（n＋1，k）＝F（n）＊（n＋1）^k；
　　　　Sum（n）＝∑［i＝1～n］A（i，k）；
　　2.利用二项式定理将（1+i）^n展开；
　　则 U（n＋1，k）＝∑［i＝0～k］C（i，k）＊F（n）＊（n）^i＋∑［i＝0～k］C（i，k）＊F（n－1）＊（n）^i
　　　　　　　　　　　　＝∑［i＝0～k］C（i，k）＊U（n，i）＋∑［i＝0～k］C（i，k）＊V（n，i）；
　　同理 V（n＋1，k）＝∑［i＝0～k］C（i，k）＊U（n，i）；
　　Sum（n）＝Sum（n－1）＋U（n，k）；
　　
　　3.状态转移用矩阵向量表示：

sum（n-1）		sum(n)
U(n,k)		U(n+1,k)
...		...
U(n,0)	=>	U(n+1,0)
V(n,k)		V(n+1,k)
...		...
V(n,0)		V(n+1,0)

代码如下：

#include <iostream>
#include <cstdio>
#include <cstring>
#include <time.h>
#include <assert.h>
#define time_ printf("time : %f
",double(clock())/CLOCKS_PER_SEC)
using namespace std;
#define maxk 40
#define MOD 1000000007
#define mod_(x) (x%MOD)

typedef long long LL;
LL n;
int k;

struct Matrix{
    int row,col;
    int m[2*maxk+5][2*maxk+5];
    Matrix(int r=83,int c=83){
        row=r;
        col=c;
        memset(m, 0, sizeof m);
    }
    void operator=(const Matrix& m2){
        memcpy(m, m2.m, sizeof m);
    }
};
Matrix operator*(const Matrix& a,const Matrix& b){
    Matrix c(a.row,b.col);
    for(int i=0;i<c.row;i++)
        for(int j=0;j<c.col;j++){
            for(int p=0;p<a.col;p++){
                c.m[i][j]=(c.m[i][j]+LL(a.m[i][p])*b.m[p][j])%MOD;
                //assert(c.m[i][j]>=0);
            }
        }
    return c;
}
Matrix C;
void set_C(){
    for(int i=k;i>=0;i--)
        for(int j=k;j>=i;j--){
            if(j==k||j==i)
                C.m[i][j]=1;
            else
                C.m[i][j]=(C.m[i+1][j]+C.m[i+1][j+1])%MOD;
            }
}
void cpy_C(Matrix& a,int pi,int pj){
    int len=k+1;
    for(int i=0;i<len;i++)
        for(int j=0;j<len;j++)
            a.m[pi+i][pj+j]=C.m[i][j];
}
void set_M(Matrix& a){
    a.m[0][0]=1,a.m[0][1]=1;
    cpy_C(a, 1, 1);
    cpy_C(a, k+2, 1);
    cpy_C(a, 1, k+2);
}
Matrix pow_mod(Matrix& a,LL n){
    if(n==1) return a;
    Matrix temp=pow_mod(a, n/2);
    temp=temp*temp;
    if(n%2){
        temp=temp*a;
    }
    return temp;
}
int pow_2(int n){
    if(n==0) return 1;
    int temp=pow_2(n/2)%MOD;
    temp=int((LL(temp)*temp)%MOD);
    if(n%2)
        temp=(temp*2)%MOD;
    return temp;
}
int main(int argc, const char * argv[]) {
    while(scanf("%lld%d",&n,&k)==2){
        if(n==1){
            printf("%d
",1);
            continue;
        }
        set_C();
        Matrix I(2*k+3,2*k+3);
        set_M(I);
        Matrix A(2*k+3,1);
        A.m[0][0]=1;
        for(int i=1;i<k+2;i++)
            A.m[i][0]=pow_2(k+1-i+1);
        for(int i=k+2;i<2*k+3;i++)
            A.m[i][0]=pow_2(2*k+2-i);
        Matrix temp=pow_mod(I, n-1);
        A=temp*A;
        printf("%d
",A.m[0][0]);
        //time_;
    }
    return 0;
}