以斐波拉契数列为例,
要从矩阵A
$$ \begin{bmatrix} f[n-1] & f[n] \end{bmatrix} $$
得到矩阵B
$$ \begin{bmatrix} f[n] & f[n+1] \end{bmatrix} $$
显然可以$$\begin{bmatrix} f[n-1] & f[n] \end{bmatrix}\times \begin{bmatrix} 0 & 1\\ 1 & 1 \end{bmatrix}=\begin{bmatrix} f[n] & f[n-1]+f[n] \end{bmatrix}=\begin{bmatrix} f[n] & f[n+1] \end{bmatrix}$$
那么要求数列第n项,
用初始矩阵
$$ \begin{bmatrix} 1 & 1 \end{bmatrix}$$
乘上
$${ \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} }^{n-1}$$
对于这道题,
$latex f[1]=f[2]=f[3]=1$
$latex f[x]=f[x-3]+f[x-1]$
显然从一个状态推至下一状态要保存3个连续的元素, 因为关系式$latex f[x]=f[x-3]+f[x-1]$一共跨过了3个元素, 于是两个矩阵就出来了:
$$A=\begin{bmatrix} f[n-2] & f[n-1] & f[n] \end{bmatrix}$$
$$C=\begin{bmatrix} f[n-1] & f[n] & f[n+1] \end{bmatrix}$$
那么怎样从A转移到下一状态C呢
我们来构造一个矩阵B
先考虑C的第一个元素, $latex f[n-1]$ ,直接从A的第二个转移来即可,所以B的第一列为0 1 0
考虑C的第2个元素,$latex f[n]$ ,直接从A的第3个转移来即可,所以B的第2列为0 0 1
考虑C的第3个元素,$latex f[n+1]$ ,跟据递推式,从A的第1个和第3个相加转移来即可, 所以B的第3列为1 0 1
综上,
$$B=\begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 0 \\ 0 & 1 & 1\end{bmatrix}$$
$$A\times B = C$$
$$\begin{bmatrix} f[n-2] & f[n-1] & f[n] \end{bmatrix}\times\begin{bmatrix} 0 & 0 & 1\\ 1 & 0 & 0 \\ 0 & 1 & 1\end{bmatrix} $$
$$= \begin{bmatrix} f[n-1] & f[n] & f[n-2]+f[n] \end{bmatrix}$$
$$ = \begin{bmatrix} f[n-1] & f[n] & f[n+1] \end{bmatrix}$$
通过矩阵乘法实现了状态的转移,接下来通过快速幂加速即可.
为了方便处理, 统一把矩阵变为正方形的, 原来没有元素的地方赋值为0即可.
#include<cstdio> #include<cstring> #define re register #define in inline #define int long long inline int read() { int s=0,b=1; char ch; do{ ch=getchar(); if(ch=='-') b=-1; }while(ch<'0' || ch>'9'); while(ch>='0' && ch<='9') { s=s*10+ch-'0'; ch=getchar(); } return b*s; } int size=3,n,T; const int mod=1000000007; struct jz{ int v[5][5]; void ql() { memset(v,0,sizeof(v)); }//矩阵清零 void csh() { ql(); for(re int i=1;i<=size;++i) v[i][i]=1; }//赋值为单位矩阵 jz operator *(const jz &t)const { jz c; c.ql(); for(re int i=1;i<=size;++i) for(re int j=1;j<=size;++j) for(re int k=1;k<=size;++k) c.v[i][j]=c.v[i][j]%mod+v[i][k]*t.v[k][j]%mod; return c; }//矩阵乘法 }; jz ksm(jz a,int p)//矩阵快速幂 { /* 非递归快速幂思想: 假设p=11,即二进制1011 则p=2^3+2^2+2^0 那么a^p转化为a^2^3*a^2^2*a^2^0 a^2^x即base,a^2^(x+1)可以很快由a^2^(x)乘a得到 只要p的二进值末尾为1, 就乘入ans即可 */ jz ans,base=a;//base一开始为a^2^0即a ans.csh();//注意ans要初始化为单位矩阵 while(p) { if(p&1) ans=ans*base; //如果p的二进制末尾为1,ans就要乘上a^2^x. 这里别去想什么奇数偶数的, 直接看p的二进制 base=base*base; //base由a^2^x转为a^2^(x+1) p>>=1; //p的末尾已处理完(a^2^x已乘入ans) } return ans; } signed main() { jz st,zy; st.ql(); zy.ql(); st.v[1][1]=1;st.v[1][2]=1;st.v[1][3]=1; //初始矩阵st为|1 1 1| zy.v[1][1]=0;zy.v[1][2]=0;zy.v[1][3]=1; zy.v[2][1]=1;zy.v[2][2]=0;zy.v[2][3]=0; zy.v[3][1]=0;zy.v[3][2]=1;zy.v[3][3]=1; /* 转移矩阵为 |0 0 1| |1 0 0| |0 1 1| f[n]~f[n+1]转移过程: |f[n-2] f[n-1] f[n]| * |0 0 1| = |f[n-1] f[n] f[n-2]+f[n]| =|f[n-1] f[n] f[n+1]| |1 0 0| |0 1 1| */ T=read(); for(re int i=1;i<=T;++i) { n=read(); jz ans=st*ksm(zy,n-1);//这里st*zy^(n-1)得到[1][1]为答案 printf("%lld\n",ans.v[1][1]%mod); } return 0; }
最新评论