Codeforces 392-C Yet Another Number Sequence (矩阵快速幂)

C. Yet Another Number Sequence

time limit per test

1 second

memory limit per test

256 megabytes

input

standard input

output

standard output

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

F1 = 1, F2 = 2, Fi = Fi - 1 + Fi - 2 (i > 2).

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

Ai(k) = Fi × ik (i ≥ 1).

In this problem, your task is to calculate the following sum: A1(k) + A2(k) + ... + An(k). The answer can be very large, so print it modulo1000000007 (109 + 7).

Input

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

Output

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

Examples

input

1 1

output

1

input

4 1

output

34

input

5 2

output

316

input

7 4

output

73825

求 前A[i]和  工作量巨大,   真是煞人;

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

另 U(n+1,k)= Fn+1*(n+1)^k; V(n+1,k)=Fn*n^k

(n+1)^k 二项展开   Fn+1 = Fn+Fn-1   ==  (n+1)^k*U(n,i)+(n+1)^k*V(n,i)

构造矩阵;   二项展开与 杨辉三角有关

2k+3 的矩阵; 三个杨辉三角;

注意 n 和 k  都用long long

#include <iostream>
#include <stdio.h>
#include <algorithm>
#include <cmath>
#include <cstring>
#include <string>
#include <queue>
#include <stack>
#include <set>
#include <bitset>
#include <vector>
using namespace std;
typedef long long ll;
const double PI=acos(-1);
const int INF=0x3f3f3f3f;
const double esp=1e-6;
const int MAXN=100;
const int N=120;
const int MOD=1e9+7;
#define mem(a,b) memset(a,b,sizeof(a))
ll gcd(ll a,ll b){ return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){ return b/gcd(a,b)*a;}
ll qpow(ll x,ll n){ll res=1;for(;n;n>>=1){if(n&1)res=(res*x)%MOD;x=(x*x)%MOD;}return res;}
ll C[202][205];
void tab_init()
{
    C[0][0]=1;
    for(int i=1;i<=50;i++)
        for(int j=0;j<=i;j++){
        	C[i][j]=(j==0)?1:(C[i-1][j]+C[i-1][j-1])%MOD;
		}
}
ll n,k;
struct Matrix{
	ll arr[N][N];
	void init()
	{
		memset(arr,0,sizeof(arr));
		for(int i=0;i<MAXN;i++)
			arr[i][i]=1;//初始化
	}
	void iinit()
	{
		memset(arr,0,sizeof(arr));
		for(int i=0;i<=k;i++)
            for(int j=0;j<=i;j++)
            {
                arr[i][j]=arr[i+k+1][j]=arr[i][j+k+1]=C[i][j];
            }
        arr[2*k+2][k]=arr[2*k+2][2*k+2]=1;
    }
};
Matrix mul(Matrix X,Matrix Y)// 矩阵乘法
{
	Matrix ans;
	for(int i=0;i<MAXN;i++)
		for(int j=0;j<MAXN;j++){
			ans.arr[i][j]=0;
			for(int k=0;k<MAXN;k++){
				ans.arr[i][j]+=X.arr[i][k]*Y.arr[k][j];
			ans.arr[i][j]%=MOD;
			}
		}
	return ans;
}
Matrix Q_pow(Matrix B,ll n)// 矩阵快速幂
{
	Matrix ans;
	ans.init();
	while(n)
	{
		if(n&1)
			ans=mul(ans,B);
		n>>=1;
		B=mul(B,B);
	}
	return ans;
}

int main()
{
    tab_init();
    while(~scanf("%I64d %I64d",&n,&k))
    {
        Matrix ans;
        ans.iinit();
        ans=Q_pow(ans,n);
        ll res=0;
        for(int i=0;i<2*k+2;i++)
            res+=(ans.arr[2*k+2][i])%MOD;
        printf("%I64d\n",res%MOD);
    }
    return 0;
}

123