洛咕题解 P1939【模板】矩阵加速(数列)[矩阵加速递推]

  • A+
所属分类:代码笔记 信息技术

题面

矩阵加速递推的原理:

首先你得会矩阵乘法与快速幂.


以斐波拉契数列为例,

要从矩阵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次方即可


对于这道题,

\(f[1]=f[2]=f[3]=1\)

\(f[x]=f[x-3]+f[x-1]\)

显然从一个状态推至下一状态要保存3个连续的元素, 因为关系式\(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的第一个元素, \(f[n-1]\) ,直接从A的第二个转移来即可,所以B的第一列为0 1 0

考虑C的第2个元素,\(f[n]\) ,直接从A的第3个转移来即可,所以B的第2列为0 0 1

考虑C的第3个元素,\(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即可

  1. #include<cstdio>
  2. #include<cstring>
  3. #define re register
  4. #define in inline
  5. #define int long long
  6. inline int read()
  7. {
  8.     int s=0,b=1;
  9.     char ch;
  10.     do{
  11.         ch=getchar();
  12.         if(ch=='-') b=-1;
  13.     }while(ch<'0' || ch>'9');
  14.     while(ch>='0' && ch<='9')
  15.     {
  16.         s=s*10+ch-'0';
  17.         ch=getchar();
  18.     }
  19.     return b*s;
  20. }//快读是不可能出锅的
  21. int size=3,n,T;
  22. const int mod=1000000007;
  23. struct jz{
  24.     int v[5][5];
  25.     void ql()
  26.         {
  27.             memset(v,0,sizeof(v));
  28.         }//矩阵清零
  29.     void csh()
  30.         {
  31.             ql();
  32.             for(re int i=1;i<=size;++i)
  33.                 v[i][i]=1;
  34.         }//赋值为单位矩阵
  35.     jz operator *(const jz &t)const
  36.         {
  37.             jz c;
  38.             c.ql();
  39.             for(re int i=1;i<=size;++i)
  40.                 for(re int j=1;j<=size;++j)
  41.                     for(re int k=1;k<=size;++k)
  42.                         c.v[i][j]=c.v[i][j]%mod+v[i][k]*t.v[k][j]%mod;
  43.             return c;
  44.         }//矩阵乘法
  45. };
  46. jz ksm(jz a,int p)//矩阵快速幂
  47. {
  48.     /*
  49.       非递归快速幂思想:
  50.       假设p=11,即二进制1011
  51.       则p=2^3+2^2+2^0
  52.       那么a^p转化为a^2^3*a^2^2*a^2^0
  53.       a^2^x即base,a^2^(x+1)可以很快由a^2^(x)乘a得到
  54.       只要p的二进值末尾为1, 就乘入ans即可
  55.     */
  56.     jz ans,base=a;//base一开始为a^2^0即a
  57.     ans.csh();//注意ans要初始化为单位矩阵
  58.     while(p)
  59.     {
  60.         if(p&1) ans=ans*base; //如果p的二进制末尾为1,ans就要乘上a^2^x. 这里别去想什么奇数偶数的, 直接看p的二进制
  61.         base=base*base; //base由a^2^x转为a^2^(x+1)
  62.         p>>=1; //p的末尾已处理完(a^2^x已乘入ans)
  63.     }
  64.     return ans;
  65. }
  66. signed main()
  67. {
  68.     jz st,zy;
  69.     st.ql();
  70.     zy.ql();
  71.     st.v[1][1]=1;st.v[1][2]=1;st.v[1][3]=1;
  72.     //初始矩阵st为|1 1 1|
  73.     zy.v[1][1]=0;zy.v[1][2]=0;zy.v[1][3]=1;
  74.     zy.v[2][1]=1;zy.v[2][2]=0;zy.v[2][3]=0;
  75.     zy.v[3][1]=0;zy.v[3][2]=1;zy.v[3][3]=1;
  76.     /*
  77.       转移矩阵为
  78.       |0 0 1|
  79.       |1 0 0|
  80.       |0 1 1|
  81.       f[n]~f[n+1]转移过程:
  82.       |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]|
  83.                              |1 0 0|
  84.                              |0 1 1|
  85.     */
  86.     T=read();
  87.     for(re int i=1;i<=T;++i)
  88.     {
  89.         n=read();
  90.         jz ans=st*ksm(zy,n-1);//转移次数玄学调一下就好.这里st*zy^(n-1)得到[1][1]为答案
  91.         printf("%lld\n",ans.v[1][1]%mod);
  92.     }
  93.     return 0;
  94. }
xcc

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: