題面

矩陣加速遞推的原理:

首先你得會矩陣乘法與快速冪.

以斐波拉契數列為例,

要從矩陣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即可

#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;
}

[LGOJ]P1939【模板】矩陣加速（數列）[矩陣加速遞推]