【学习笔记】常系数线性齐次递推

应用(LGP4723)

求一个满足(k)阶线性齐次递推数列(a_i)的第(n)项；

及：(a_n = sum_{i=1}^{k} f_i imes a_{n-i}) ;

$N le 10^9 , K le 32000 $ ；

矩阵

• 矩阵的秩(rk(A))：
  • (rk(A) = 0 Leftrightarrow A=0) ;
  • (rk(A+B) le rk(A) + rk(B)) ;
  • (rk(AB) le min(rk(A),rk(B))) ;
• 行列式按行展开：
  • (det(A) = sum_{j=1}^{n} a_{i,j}A_{i,j})
  • (A_{i,j})为行列式的余子式；
  • 利用行列式的分拆性质可以证明；

特征多项式

• 对于(n)阶方阵(A)；
• 若存在列向量(v)和常数(lambda) ，满足 (lambda vec{v} = A vec{v});
• 则称(v)为(A)的特征向量;
• 若(k = rk(A))，则存在(k)个线性无关的(v)；
• 变换形式：((lambda I - A) vec{v} = 0) ;
• 有解的充要条件是(det(lambda I - A) = 0);
• 所以可以得到一个关于 $ lambda $ 的 $ n $ 阶方程 $ F(lambda)=0 $ ， $ F(x) $ 称为 $ A $ 的特征多项式；
• Cayley-Hamilton定理
• 矩阵的特征多项式为矩阵的零化多项式：(F(A) =0) ；
• 证明：https://blog.csdn.net/qq_35649707/article/details/78688235 ；

线性递推

• 设递推数列 (a_{0} - a_{k-1}) 的初始列向量为(G)(下标从下到上)，(A)为(k)阶递推矩阵，实际上就是要求([A^nG]_{k,1}) ；

• 考虑如何求(A^n)，注意到(F(A)=0)；

• 设(A^n = P(A) F(A)+ Q(A));

• 那么相当于(A^n)减去了很多个0矩阵，(Q(A))和它意义相同，只需要算多项式(Q(A))即可；

• 这部分可以快速幂+多项式取模完成，复杂度(O(k log k log m))或者$k ^ 2 log n $；

• 考虑如何求(A)的特征多项式，利用定义$det(lambda I - A) $的形式为:

  [egin{pmatrix} lambda - f_{1} & -f_2 & -f_3 & cdots & -f_{k-1} & -f_{k} \ 1 & lambda & 0 & cdots & 0 & 0 \ 0 & 1 & lambda & cdots & 0 & 0 \ vdots & vdots & vdots & ddots & vdots & vdots \ 0 & 0 & 0 & cdots & lambda & 0 \ 0 & 0 & 0 & cdots & 1 & lambda \ end{pmatrix}* ]

• 按第一行展开，会剩下一个下三角行列式，值为对角线乘积；

• 得$F(lambda) = lambda^k - sum_{i=1}^{k}f_i lambda^{k-i} $ ;

• [g_{n} = [A^n G]_{k,1} = sum_{i=0}^{k-1}[Q_{i} A^i imes G]_{k,1} =sum_{i=0}^{k-1}Q_{i}[A^{i}G]_{k,1} = sum_{i=0}^{k-1}Q_{i} a_{i} ]

• 总复杂度就是求(Q)的复杂度，另外如果前(k)项没有给出的话似乎可以利用生成函数求；

  #include<bits/stdc++.h>
  #define mod 998244353 
  #define ll long long
  using namespace std;
  const int N=1000010;
  int n,k,m,f[N],h[N],a[N],b[N],ib[N];
  char gc(){
      static char*p1,*p2,s[1000000];
      if(p1==p2)p2=(p1=s)+fread(s,1,1000000,stdin);
      return(p1==p2)?EOF:*p1++;
  }
  int rd(){
      int x=0,f=1;char c=gc();
      while(c<'0'||c>'9'){if(c=='-')f=-1;c=gc();}
      while(c>='0'&&c<='9'){x=(x<<1)+(x<<3)+c-'0',c=gc();}
      return x*f;
  }
  int pw(int x,int y){
      int re=1;
      if(y<0)y+=mod-1;
      while(y){
          if(y&1)re=(ll)re*x%mod;
          y>>=1;x=(ll)x*x%mod;
      }
      return re;
  }
  void inc(int&x,int y){x+=y;if(x>=mod)x-=mod;}
  namespace poly{
      const int G=3;
      int rev[N],L;
      void ntt(int*A,int len,int f){
          for(L=0;(1<<L)<len;++L);
          for(int i=0;i<len;++i){
              rev[i]=(rev[i>>1]>>1)|((i&1)<<(L-1));
              if(i<rev[i])swap(A[i],A[rev[i]]);
          }
          for(int i=1;i<len;i<<=1){
              int wn=pw(G,f*(mod-1)/(i<<1));
              for(int j=0;j<len;j+=i<<1){
                  int w=1;
                  for(int k=0;k<i;++k,w=(ll)w*wn%mod){
                      int x=A[j+k],y=(ll)w*A[j+k+i]%mod;
                      A[j+k]=(x+y)%mod,A[j+k+i]=(x-y+mod)%mod;
                  }
              }
          }
          if(!~f){
              int iv=pw(len,mod-2);
              for(int i=0;i<len;++i)A[i]=(ll)A[i]*iv%mod;
          }
      }
      void cls(int*A,int l,int r){for(int i=l;i<r;++i)A[i]=0;}
      void cpy(int*A,int*B,int l){for(int i=0;i<l;++i)A[i]=B[i];}
      void inv(int*A,int*B,int l){
          if(l==1){B[0]=pw(A[0],mod-2);return;}
          static int t[N];
          int len=l<<1;
          inv(A,B,l>>1);
          cpy(t,A,l);cls(t,l,len);
          ntt(t,len,1);ntt(B,len,1);
          for(int i=0;i<len;++i)B[i]=(ll)B[i]*(2-(ll)t[i]*B[i]%mod+mod)%mod;
          ntt(B,len,-1);cls(B,l,len);
      }
      void pmod(int*A){
          static int t[N];
          int l=k+1,len=1;while(len<=(k<<1))len<<=1;
          cpy(t,A,(k<<1)+1);
          reverse(t,t+(k<<1)+1);
          cls(t,l,len);
          ntt(t,len,1);
          for(int i=0;i<len;++i)t[i]=(ll)t[i]*ib[i]%mod;
          ntt(t,len,-1);
          cls(t,l,len);
          reverse(t,t+l);
          ntt(t,len,1);
          for(int i=0;i<len;++i)t[i]=(ll)t[i]*b[i]%mod;
          ntt(t,len,-1);
          cls(t,l,len);
          for(int i=0;i<k;++i)A[i]=(A[i]-t[i]+mod)%mod;
          cls(A,k,len);
      }
      void pow(int*A,int n){
          if(n==1){cls(A,0,k+1);A[1]=1;return;}
          pow(A,n>>1);
          int len=1;while(len<=(k<<1))len<<=1;
          ntt(A,len,1);
          for(int i=0;i<len;++i)A[i]=(ll)A[i]*A[i]%mod;
          ntt(A,len,-1);
          pmod(A);
          if(n&1){
              for(int i=k;i;--i)A[i]=A[i-1];A[0]=0;
              pmod(A);
          }
      }
  }
  int main(){
  //	freopen("P4723.in","r",stdin);
  //	freopen("P4723.out","w",stdout);
      n=rd();k=rd();
      for(int i=1;i<=k;++i)f[i]=(mod+rd())%mod;
      for(int i=0;i<k;++i)h[i]=(mod+rd())%mod;
      for(int i=a[k]=b[k]=1;i<=k;++i)a[k-i]=b[k-i]=(mod-f[i])%mod;
      int len=1;while(len<=(k<<1))len<<=1;
      reverse(a,a+k+1);
      poly::inv(a,ib,len);
      poly::cls(ib,k+1,len);
      /*{
          for(int i=0;i<len;++i)printf("%d ",b[i]);puts("");
          for(int i=0;i<len;++i)printf("%d ",ib[i]);puts("");
      }*/
      poly::ntt(b,len,1);
      poly::ntt(ib,len,1);
      poly::pow(a,n);
      int ans=0;
      for(int i=0;i<k;++i)inc(ans,(ll)a[i]*h[i]%mod);
      printf("%d
  ",ans);
      return 0;
  }