[LGOJ]P1939模板矩阵加速（数列）[矩阵加速递推]

Posted 2021-02-06 cgazn

题面

矩阵加速递推的原理:

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

以斐波拉契数列为例,

要从矩阵A

[ egin{bmatrix} f[n-1] & f[n]  end{bmatrix} ]

得到矩阵B

[ egin{bmatrix} f[n] & f[n+1]  end{bmatrix} ]

显然可以[egin{bmatrix} f[n-1] & f[n] end{bmatrix} imes egin{bmatrix}  0 & 1\ 1 & 1 end{bmatrix}=egin{bmatrix} f[n] & f[n-1]+f[n] end{bmatrix}=egin{bmatrix} f[n] & f[n+1] end{bmatrix}]

那么要求数列第n项,

用初始矩阵
[ egin{bmatrix}  1 & 1 end{bmatrix}]
乘上
[{ egin{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=egin{bmatrix}  f[n-2] & f[n-1] & f[n] end{bmatrix}]

[C=egin{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=egin{bmatrix}  0 & 0 & 1\ 1 & 0 & 0 \ 0 & 1 & 1end{bmatrix}]

[A imes B = C]

[egin{bmatrix}  f[n-2] & f[n-1] & f[n] end{bmatrix} imesegin{bmatrix}  0 & 0 & 1\ 1 & 0 & 0 \ 0 & 1 & 1end{bmatrix} ]

[= egin{bmatrix}  f[n-1] & f[n] & f[n-2]+f[n] end{bmatrix}]

[ = egin{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
",ans.v[1][1]%mod);
    }
    return 0;
}