[洛谷P1939]【模板】矩阵加速（数列）[矩阵加速递推]

题面

以斐波拉契数列为例,

要从矩阵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;
}