ShinriiTin's Blog - 博主是sb

suansuansuan

已知数列f:f0=f1=1,fn+1 = fn+fn-1+1(n>=1).

求Sn = ∑i∈[0,n] fik(n<264,k∈[1,13])

 

k=1时:

构造转移矩阵A:{A0,A1,A2,A3},列向量B:(Sn,fn,fn-1,1),使得AxB=(Sn+1,fn+1,fn,1).

        A= (1,1,1,1)

        A= (0,1,1,1)

        A= (0,1,0,0)

        A= (0,0,0,1)

 

k=2时:

构造转移矩阵A:{A0,A1,A2,A3,A4,A5,A6},列向量B:(Sn,fn2,fnfn-1,fn-12,fn,fn-1,1),使得AxB=(Sn+1,fn+12,fn+1fn,fn2,fn+1,fn,1).

       A0 = (1,1,2,1,2,2,1)

       A1 = (0,1,2,1,2,2,1)

       A2 = (0,1,1,0,1,0,0)

       A3 = (0,1,0,0,0,0,0)

       A4 = (0,0,0,0,1,1,1)

       A5 = (0,0,0,0,1,0,0)

       A6 = (0,0,0,0,0,0,1)

 

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

并且fn+1afn-1b(a+b<=k)都能用fnifn-1j(i+j<=k)来线性表示.

 

fn+1afnb = (fn+fn-1+1)afnb

              = ∑i∈[0,a] C(a,i)*(fn+fn-1)afnb

                  = ∑i∈[0,a] C(a,i)*∑j∈[0,i] C(i,j)*fnjfn-1i-jfnb

               = i∈[0,a] ∑j∈[0,i] C(a,i)*C(i,j)*fnb-jfn-1i-j 

           

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

然后就可以用快速幂加速矩阵乘法O(k6logn)求解Sn

 

 

#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;
}