已知数列f：f₀=f₁=1,f_n+1 = f_n+f_n-1+1(n>=1).

求S_n = ∑_i∈[0,n] f_i^k(n<2⁶⁴,k∈[1,13])

k=1时：

构造转移矩阵A：{A₀,A₁,A₂,A₃},列向量B：(S_n,f_n,f_n-1,1),使得AxB=(S_n+1,f_n+1,f_n,1).

        A₀ = (1,1,1,1)

        A₁ = (0,1,1,1)

        A₂ = (0,1,0,0)

        A₃ = (0,0,0,1)

k=2时：

构造转移矩阵A：{A₀,A₁,A₂,A₃,A₄,A₅,A₆},列向量B：(S_n,f_n²,f_nf_n-1,f_n-1²,f_n,f_n-1,1),使得AxB=(S_n+1,f_n+1²,f_n+1f_n,f_n²,f_n+1,f_n,1).

       A₀ = (1,1,2,1,2,2,1)

       A₁ = (0,1,2,1,2,2,1)

       A₂ = (0,1,1,0,1,0,0)

       A₃ = (0,1,0,0,0,0,0)

       A₄ = (0,0,0,0,1,1,1)

       A₅ = (0,0,0,0,1,0,0)

       A₆ = (0,0,0,0,0,0,1)

容易发现我们所用到的项除了Sn之外就只有f_nⁱf_n-1^j(i+j<=k)

并且f_n+1^af_n-1^b(a+b<=k)都能用f_nⁱf_n-1^j(i+j<=k)来线性表示.

f_n+1^af_n^b = (f_n+f_n-1+1)^af_n^b

              = ∑_i∈[0,a] C(a,i)*(f_n+f_n-1)^af_n^b

                  = ∑_i∈[0,a] C(a,i)*∑_j∈[0,i] C(i,j)*f_n^jf_n-1^i-jf_n^b

               = ∑_i∈[0,a] ∑_j∈[0,i] C(a,i)*C(i,j)*f_n^b-jf_n-1^i-j 

这样就能O(k⁴)预处理出转移矩阵A

然后就可以用快速幂加速矩阵乘法O(k⁶logn)求解S_n了

#include <map>
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <cstdio>
#define MAXN 110
#define Mod 1000000000
typedef long long ll;
typedef unsigned long long ull;
using namespace std;

typedef pair<int,int>cp;

inline void set_IO(){
	freopen("suansuansuan.in","r",stdin);
	freopen("suansuansuan.out","w",stdout);
}

ull n;

int k,C[20][20];

int A[MAXN][MAXN],B[MAXN];

inline void mul(int A[MAXN][MAXN],int B[MAXN][MAXN],int n){
	static int C[MAXN][MAXN];
	for(int i=0;i<=n;++i)
		for(int j=0;j<=n;++j){
			C[i][j]=0;
		}
	for(int i=0;i<=n;++i)
		for(int k=0;k<=n;++k){
			if(!A[i][k])continue;
			for(int j=0;j<=n;++j){
				C[i][j]+=(ll)A[i][k]*B[k][j]%Mod;
				C[i][j]%=Mod;
			}
		}
	for(int i=0;i<=n;++i)
		for(int j=0;j<=n;++j){
			A[i][j]=C[i][j];
		}
}

inline void qpow(int x[MAXN][MAXN],int n,ull k){
	static int ans[MAXN][MAXN];
	for(int i=0;i<=n;++i)
		for(int j=0;j<=n;++j){
			ans[i][j]=bool(i==j);
		}
	for(;k;mul(x,x,n),k>>=1)
		if(k&1)mul(ans,x,n);
	for(int i=0;i<=n;++i)
		for(int j=0;j<=n;++j){
			x[i][j]=ans[i][j];
		}
}

map<cp,int>id;

int Id;

int main(){
	set_IO();
	cin>>n>>k;
	if(n<2){
		printf("%09d",(int)n+1);
		return 0;
	}
	for(int i=0;i<=k;++i){
		C[i][0]=1;
		for(int j=1;j<=i;++j){
			C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mod;
		}
	}
	for(int i=k;~i;--i)
		for(int j=i;~j;--j){
			id[cp(j,i-j)]=++Id;
		}
	for(int a=0;a<=k;++a)
		for(int b=0;a+b<=k;++b){
			int x=id[cp(a,b)];
			for(int i=0;i<=a;++i)
				for(int j=0;j<=i;++j){
					int y=id[cp(b+j,i-j)];
					A[x][y]=(ll)C[a][i]*C[i][j]%Mod;
				}
		}
	A[0][0]=1;
	for(int i=1,x=id[cp(k,0)];i<=Id;++i){
		A[0][i]=A[x][i];
	}
	qpow(A,Id,n-1);
	B[0]=2;
	for(int i=1;i<=Id;++i){
		B[i]=1;
	}
	int ans=0;
	for(int i=0;i<=Id;++i){
		ans=(ans+(ll)A[0][i]*B[i]%Mod)%Mod;
	}
	printf("%09d\n",ans);	
	return 0;
}