数论-矩阵快速幂及其简单构造方法

[用途]

n较大时，求简单递推式$F_n=\sum _{i=1}^{n-1}{k}_i{F}_{n-i}+something,(k_i\in Z)$

[构造形式]

构建矩阵X，A，B,使矩阵X乘以矩阵A等于矩阵B ， (注意一般情况下矩阵不满足交换律！！！）

[构造方法]

• 含$F_{n-i}$的项依次按照顺序填入矩阵A中
• $something$中常数项与$(n-1{)}^k$项也依次填入矩阵A中
• $something$含$n^k$的项根据二项式定理
  $(1+x{)}^n=\sum _{i=0}^n{C}_n^i{x}^n$
  则$n^k=(1+(n-1){)}^k=\sum _{i=0}^k{C}_k^i(n-1{)}^i$
  这样就转化为了含$(n-1{)}^k$的项

$构建后的矩阵A类似于：\left[\begin{array}{ccccc}{F}_{n-1}& 0& 0& ...& 0\\ F_{n-2}& 0& 0& ...& 0\\ ...& 0& 0& ...& 0\\ (n-1{)}^k& 0& 0& ...& 0\\ (n-1{)}^{k-1}& 0& 0& ...& 0\\ ...& 0& 0& ...& 0\\ (n-1{)}^0& 0& 0& ...& 0\end{array}\right]$

$构建后的矩阵B类似于：\left[\begin{array}{ccccc}{F}_n& 0& 0& ...& 0\\ F_{n-1}& 0& 0& ...& 0\\ ...& 0& 0& ...& 0\\ n^k& 0& 0& ...& 0\\ n^{k-1}& 0& 0& ...& 0\\ ...& 0& 0& ...& 0\\ n^0& 0& 0& ...& 0\end{array}\right]$
然后根据矩阵乘法求X
$B[i][j]为矩阵X的第i行与矩阵A的第j列对应相乘所得的和,由已知递推式求系数即可$

[代码]

以hdu6470为例

$F_n=2F_{n-2}+F_{n-1}+n^3$
根据上述有
$A=\left[\begin{array}{cccccc}{F}_{n-1}& 0& 0& 0& 0& 0\\ F_{n-2}& 0& 0& 0& 0& 0\\ (n-1{)}^3& 0& 0& 0& 0& 0\\ (n-1{)}^2& 0& 0& 0& 0& 0\\ (n-1{)}^1& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\end{array}\right]$ , $B=\left[\begin{array}{cccccc}{F}_n& 0& 0& 0& 0& 0\\ F_{n-1}& 0& 0& 0& 0& 0\\ n^3& 0& 0& 0& 0& 0\\ n^2& 0& 0& 0& 0& 0\\ n^1& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\end{array}\right]$

则 $X=\left[\begin{array}{cccccc}1& 2& 1& 3& 3& 1\\ 1& 0& 0& 0& 0& 0\\ 0& 0& 1& 3& 3& 1\\ 0& 0& 0& 1& 2& 1\\ 0& 0& 0& 0& 1& 1\\ 0& 0& 0& 0& 0& 1\end{array}\right]$

设基准点是F[2],矩阵$A=\left[\begin{array}{cccccc}2(即F[2])& 0& 0& 0& 0& 0\\ 1即(F[1])& 0& 0& 0& 0& 0\\ 8& 0& 0& 0& 0& 0\\ 4& 0& 0& 0& 0& 0\\ 2& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0\end{array}\right]$

求F[n],矩阵$B=X^{n-2}\ast A$，即F[n]=B[0][0]

#include<cstdio>
#include<cstring>
using namespace std;
typedef long long LL;
const LL N=6;
const LL mod=123456789;
struct jz{
    LL A[N][N];
    inline jz friend operator*(jz x,jz y){
        jz ans;
        memset(ans.A,0,sizeof(ans.A));
        for(LL i=0;i<N;++i)
            for(LL j=0;j<N;++j)
                for(LL k=0;k<N;++k)
                    ans.A[i][j]=(ans.A[i][j]+x.A[i][k]*y.A[k][j])%mod;
        return ans;
    }
};
jz jzksm(jz x,LL n)
{
    jz ans;
    memset(ans.A,0,sizeof(ans.A));
    for(LL i=0;i<N;++i)ans.A[i][i]=1;
    while(n){
        if(n&1)ans=ans*x;
        n>>=1;
        x=x*x;
    }
    return ans;
}
int main()
{
    jz first={
        2,0,0,0,0,0,
        1,0,0,0,0,0,
        8,0,0,0,0,0,
        4,0,0,0,0,0,
        2,0,0,0,0,0,
        1,0,0,0,0,0,
    };
    jz x={
        1,2,1,3,3,1,
        1,0,0,0,0,0,
        0,0,1,3,3,1,
        0,0,0,1,2,1,
        0,0,0,0,1,1,
        0,0,0,0,0,1,
    };
    LL n,T;
    scanf("%lld",&T);
    while(T--){
        scanf("%lld",&n);
        jz ans=jzksm(x,n-2)*first;
        printf("%lld\n",ans.A[0][0]);
    }
    return 0;
}