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

最新推荐文章于 2024-03-09 11:39:09 发布

原创最新推荐文章于 2024-03-09 11:39:09 发布 · 276 阅读

1 ·

CC 4.0 BY-SA版权

数学专栏收录该内容

3 篇文章

订阅专栏

该博客介绍了如何运用矩阵快速幂的方法解决一类涉及二项式展开的数列递推求和问题。通过构造特定的矩阵并进行矩阵运算，可以高效地计算出给定指数k下的前缀和。博主详细阐述了从递推关系建立矩阵到矩阵快速幂应用的全过程，并给出了C++实现的代码示例。

Yet Another Number Sequence

求： $∑1nfi∗ik\sum_1^n{f_i*i^k}$
观察到k很小，对于k很小的时候可以通过二项式展开构造大小与k有关的矩阵进行递推。
$fi+1∗(i+1)k=(fi+fi−1)∗(i+1)k=∑j=0k(Ckj∗ij∗fi+Ckj∗ij∗fi−1){f_{i+1}*{(i+1)}^k}=(f_{i}+f_{i-1})*{(i+1)}^k=\sum_{j=0}^k{(C_k^j*i^j*f_{i}+C_k^j*i^j*f_{i-1})}$
令 $G_{i,j}=f_{i}*i^j,H_{i,j}=f_{i-1}*i^j$
则 $Gi+1,k=∑j=0k(Ckj∗Gi,j+Ckj∗Hi,j)G_{i+1,k}=\sum_{j=0}^k{(C_k^j*G_{i,j}+C_k^j*H_{i,j})}$
那么得到了G的递推式，但是式子里有H是没有递推的，我们尝试将H纳入递推范畴
$Hi+1,k=fi∗(i+1)k=∑j=0kCkj∗fi∗ij=∑j=0kCkj∗Gi,jH_{i+1,k}=f_{i}*(i+1)^k=\sum_{j=0}^k{C_k^j*f_{i}*i^j}=\sum_{j=0}^k{C_k^j*G_{i,j}}$
那么所有需要的递推式都出来了，这个时候就可以愉快的构造矩阵了。
观察到首先有两列数是要进行递推的分别是 $G_{i,0}~G_{i,k}和H_{i,0}~H_{i,k}$ ，分别将他们放进矩阵
那么左边
$(Gi+1,0Gi+1,1...Gi+1,kHi+1,0Hi+1,1...Hi+1,k)\begin{pmatrix} 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{pmatrix}$
构造递推的时候的系数矩阵，就是将推出来的系数填进去就好了
$(111111121121..........Ck0Ck1..CkkCk0Ck1..Ckk111121..........Ck0Ck1..Ckk)\begin{pmatrix} 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{pmatrix}$
那么
$(Gi+1,0Gi+1,1...Gi+1,kHi+1,0Hi+1,1...Hi+1,k)==(111111121121..........Ck0Ck1..CkkCk0Ck1..Ckk111121..........Ck0Ck1..Ckk)∗(Gi,0Gi,1...Gi,kHi,0Hi,1...Hi,k)\begin{pmatrix} 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{pmatrix}==\begin{pmatrix} 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{pmatrix}*\begin{pmatrix} G_{i,0} \\ G_{i,1} \\ . \\ . \\ . \\ G_{i,k} \\ H_{i,0} \\ H_{i,1} \\ . \\ . \\ . \\ H_{i,k} \\ \end{pmatrix}$
这样就得到了每一项的递推式，问题问我们前缀和所以我们加一维sum求前缀和
得到了式子
$(Gi+1,0Gi+1,1...Gi+1,kHi+1,0Hi+1,1...Hi+1,ksumi)==(110111101211210..........0Ck0Ck1..CkkCk0Ck1..Ckk0101101210..........0Ck0Ck1..Ckk011)∗(Gi,0Gi,1...Gi,kHi,0Hi,1...Hi,ksumi−1)\begin{pmatrix} G_{i+1,0} \\ G_{i+1,1} \\ . \\ . \\ . \\ G_{i+1,k} \\ H_{i+1,0} \\ H_{i+1,1} \\ . \\ . \\ . \\ H_{i+1,k} \\ sum_{i} \end{pmatrix}==\begin{pmatrix} 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{pmatrix}*\begin{pmatrix} G_{i,0} \\ G_{i,1} \\ . \\ . \\ . \\ G_{i,k} \\ H_{i,0} \\ H_{i,1} \\ . \\ . \\ . \\ H_{i,k} \\ sum_{i-1} \end{pmatrix}$
.代表延续前面的规则，空的代表0
剩下的就是典型的矩阵快速幂，就不用多说了吧。
求答案的时候只要访问当前的前缀和 $sum_{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);
	}
		
}