Yet Another Number Sequence（矩阵快速幂、构造矩阵典型）

Yet Another Number Sequence

求：${\sum }_1^n{f}_i\ast i^k$
观察到k很小，对于k很小的时候可以通过二项式展开构造大小与k有关的矩阵进行递推。
$f_{i+1}\ast {(i+1)}^k=(f_i+f_{i-1})\ast {(i+1)}^k={\sum }_{j=0}^k(C_k^j\ast i^j\ast f_i+C_k^j\ast i^j\ast f_{i-1})$
令$G_{i,j}=f_i\ast i^j,H_{i,j}=f_{i-1}\ast i^j$
则$G_{i+1,k}={\sum }_{j=0}^k(C_k^j\ast G_{i,j}+C_k^j\ast H_{i,j})$
那么得到了G的递推式，但是式子里有H是没有递推的，我们尝试将H纳入递推范畴
$H_{i+1,k}=f_i\ast (i+1{)}^k={\sum }_{j=0}^k{C}_k^j\ast f_i\ast i^j={\sum }_{j=0}^k{C}_k^j\ast G_{i,j}$
那么所有需要的递推式都出来了，这个时候就可以愉快的构造矩阵了。
观察到首先有两列数是要进行递推的分别是$G_{i,0}{G}_{i,k}和H_{i,0}{H}_{i,k}$，分别将他们放进矩阵
那么左边
$\left(\begin{array}{c}{G}_{i+1,0}\\ G_{i+1,1}\\ .\\ .\\ .\\ G_{i+1,k}\\ H_{i+1,0}\\ H_{i+1,1}\\ .\\ .\\ .\\ H_{i+1,k}\end{array}\right)$
构造递推的时候的系数矩阵，就是将推出来的系数填进去就好了
$\left(\begin{array}{ccccccccccc}1& & & & & 1& & & & & \\ 1& 1& & & & 1& 1& & & & \\ 1& 2& 1& & & 1& 2& 1& & & \\ .& .& .& .& .& .& .& .& .& .& \\ C_k^0& C_k^1& .& .& C_k^k& C_k^0& C_k^1& .& .& C_k^k& \\ 1& & & & & & & & & & \\ 1& 1& & & & & & & & & \\ 1& 2& 1& & & & & & & & \\ .& .& .& .& .& .& .& .& .& .& \\ C_k^0& C_k^1& .& .& C_k^k& & & & & & \\ & & & & & & & & & & \end{array}\right)$
那么
$\left(\begin{array}{c}{G}_{i+1,0}\\ G_{i+1,1}\\ .\\ .\\ .\\ G_{i+1,k}\\ H_{i+1,0}\\ H_{i+1,1}\\ .\\ .\\ .\\ H_{i+1,k}\end{array}\right)==\left(\begin{array}{ccccccccccc}1& & & & & 1& & & & & \\ 1& 1& & & & 1& 1& & & & \\ 1& 2& 1& & & 1& 2& 1& & & \\ .& .& .& .& .& .& .& .& .& .& \\ C_k^0& C_k^1& .& .& C_k^k& C_k^0& C_k^1& .& .& C_k^k& \\ 1& & & & & & & & & & \\ 1& 1& & & & & & & & & \\ 1& 2& 1& & & & & & & & \\ .& .& .& .& .& .& .& .& .& .& \\ C_k^0& C_k^1& .& .& C_k^k& & & & & & \\ & & & & & & & & & & \end{array}\right)\ast \left(\begin{array}{c}{G}_{i,0}\\ G_{i,1}\\ .\\ .\\ .\\ G_{i,k}\\ H_{i,0}\\ H_{i,1}\\ .\\ .\\ .\\ H_{i,k}\end{array}\right)$
这样就得到了每一项的递推式，问题问我们前缀和所以我们加一维sum求前缀和
得到了式子
$\left(\begin{array}{c}{G}_{i+1,0}\\ G_{i+1,1}\\ .\\ .\\ .\\ G_{i+1,k}\\ H_{i+1,0}\\ H_{i+1,1}\\ .\\ .\\ .\\ H_{i+1,k}\\ su{m}_i\end{array}\right)==\left(\begin{array}{ccccccccccc}1& & & & & 1& & & & & 0\\ 1& 1& & & & 1& 1& & & & 0\\ 1& 2& 1& & & 1& 2& 1& & & 0\\ .& .& .& .& .& .& .& .& .& .& 0\\ C_k^0& C_k^1& .& .& C_k^k& C_k^0& C_k^1& .& .& C_k^k& 0\\ 1& & & & & & & & & & 0\\ 1& 1& & & & & & & & & 0\\ 1& 2& 1& & & & & & & & 0\\ .& .& .& .& .& .& .& .& .& .& 0\\ C_k^0& C_k^1& .& .& C_k^k& & & & & & 0\\ & & & & 1& & & & & & 1\end{array}\right)\ast \left(\begin{array}{c}{G}_{i,0}\\ G_{i,1}\\ .\\ .\\ .\\ G_{i,k}\\ H_{i,0}\\ H_{i,1}\\ .\\ .\\ .\\ H_{i,k}\\ su{m}_{i-1}\end{array}\right)$
.代表延续前面的规则，空的代表0
剩下的就是典型的矩阵快速幂，就不用多说了吧。
求答案的时候只要访问当前的前缀和$su{m}_{n-1}$在加上当前的$G_n^k$就可以了

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<vector>
#include<queue>
using namespace std;
typedef long long LL;
const int maxn=305;
const LL INF=1e18+7;
const LL mod=1e9+7;
const double PI=acos(-1);
const double eps=1e-6;
#define pii pair<LL,int>
#define mp(x,y) make_pair(x,y)
#define sfi(a) scanf("%d",&a)
#define sfl(a) scanf("%lld",&a)
#define sff(a) scanf("%lf",&a)
#define sfs(a) scanf("%s",a)
LL C[110][110];
void init()
{
	memset(C,0,sizeof(C));
	C[0][0]=1;
	for(int i=1;i<100;++i)
	{
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;++j)
		{
			C[i][j]=(C[i-1][j]+C[i-1][j-1])%mod;
		}
	}
}
struct matrix
{
	LL jie[110][110];
	LL n,m;
	void init()
	{
		for(int i=1;i<=n;++i)
		for(int j=1;j<=m;++j)jie[i][j]=0;
	}
	matrix operator * (const matrix &a)const
	{
		matrix b;b.n=n;b.m=a.m;
		b.init();
		for(int i=1;i<=n;++i)
		for(int j=1;j<=a.m;++j)
		{
			b.jie[i][j]=0;
			for(int k=1;k<=m;++k)
			{
				b.jie[i][j]=((jie[i][k]*a.jie[k][j])%mod+b.jie[i][j])%mod;
			}
		}
		return b;
	}
}t1,t2,t3;
matrix mpow(matrix a,LL b)
{
	matrix kk;kk.n=kk.m=a.n;
	kk=a;b--;
	while(b)
	{
		if(b&1)kk=kk*a;
		a=a*a;
		b>>=1;
	}
	return kk;
}
LL n,m;
void gao()
{
	t1.n=2*m+3;t1.m=1;
	t1.init();
	for(int i=1;i<=t1.n-1;++i)t1.jie[i][1]=1;
	t1.jie[t1.n][1]=0;
	t2.n=t2.m=2*m+3;
	t2.init();
	for(int i=1;i<=m+1;++i)
	{
		for(int j=1;j<=i;++j)
		{
			t2.jie[i][j]=t2.jie[i+1+m][j]=t2.jie[i][1+m+j]=C[i-1][j-1];
		}
	}
	t2.jie[2*m+3][1+m]=1;t2.jie[2*m+3][2*m+3]=1;
}
int main()
{
	init();
	while(sfl(n)!=EOF)
	{
		sfl(m);
		if(n==1)
		{
			printf("1\n");return 0;
		}
		gao();
		t3=mpow(t2,n-1)*t1;
		LL ans=0;
		ans=(t3.jie[m+1][1]%mod+t3.jie[m*2+3][1]%mod+mod)%mod;
		printf("%lld\n",ans);
	}
		
}