CodeForces - 393E Yet Another Number Sequence

Discription

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).

Examples

Input

1 1

Output

1

Input

4 1

Output

34

Input

5 2

Output

316

Input

7 4

Output

73825

之前做过一道 需要求 f[i] * i 生成函数形式的题，在那道题的题解里(http://www.cnblogs.com/JYYHH/p/8822572.html)已经证明过了这个玩意的生成函数的分母是 (1-x-x^2)^2 。。

当然，更普遍的，我们可以证明 f[i] * i^k 的生成函数表示的分母是 (1-x-x^2)^(k+1) ,这个用二项式定理解一下多项式闭形式就ojbk了。    于是我们可以得到这个函数的递推式，于是就可以直接预处理出前若干项然后直接用矩阵做了。

你问我它还要求前缀和？？？？  这个是矩阵的常规操作啊23333

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int maxn=93;
const int ha=1000000007;
int K,A[maxn],n,F[maxn];
inline int add(int x,int y){ x+=y; return x>=ha?x-ha:x;}
inline int ksm(int x,int y){ int an=1; for(;y;y>>=1,x=x*(ll)x%ha) if(y&1) an=an*(ll)x%ha; return an;}
struct node{
	int a[maxn][maxn];
	inline void clear(){ memset(a,0,sizeof(a));}
	inline void BASE(){ clear(); for(int i=n+1;i;i--) a[i][i]=1;}
	node operator *(const node &u)const{
		node r; r.clear();
		for(int k=n+1;k;k--)
		    for(int i=n+1;i;i--)
		        for(int j=n+1;j;j--) r.a[i][j]=add(r.a[i][j],a[i][k]*(ll)u.a[k][j]%ha);
		return r;
	}
}X,ANS;
ll N;

inline void init(){
	n=2,A[0]=1,A[1]=A[2]=ha-1;
	for(int i=1;i<=K;i++){
		for(int j=n;j>=0;j--){
			A[j+2]=add(A[j+2],ha-A[j]);
			A[j+1]=add(A[j+1],ha-A[j]);
		}
		n+=2;
	}
	for(int i=1;i<=n;i++) A[i]=ha-A[i];

	F[0]=F[1]=1;
	for(int i=2;i<=n;i++) F[i]=add(F[i-1],F[i-2]);
	for(int i=1;i<=n;i++) F[i]=F[i]*(ll)ksm(i,K)%ha;
}

inline void build(){
	X.clear(),ANS.BASE();
	for(int i=1;i<n;i++) X.a[i][i+1]=1;
	for(int i=1;i<=n;i++) X.a[i][1]=X.a[i][n+1]=A[i];
	X.a[n+1][n+1]=1;
}

inline int calc(){
	int ans=0,cnt=0;
	if(N<=n) for(int i=1;i<=N;i++) ans=add(ans,F[i]);
	else{
		N-=n,memset(A,0,sizeof(A));
		for(int i=1;i<=n;i++){
			A[i]=F[n-i+1];
			A[n+1]=add(A[n+1],F[i]);
		}

		for(;N;N>>=1,X=X*X) if(N&1) ANS=ANS*X;

		for(int i=n+1;i;i--) ans=add(ans,A[i]*(ll)ANS.a[i][n+1]%ha);
	}
	return ans;
}

inline void solve(){
    cin>>N>>K;
	init(),build();
	printf("%d\n",calc());
}

int main(){
	solve();
	return 0;
}

原文地址：https://www.cnblogs.com/JYYHH/p/8869738.html