[关闭]
@KirinBill 2018-07-13T13:15:00.000000Z 字数 6414 阅读 655

线性递推式的计算

线性代数

目录


特征多项式

对于常系数齐次递推式,虽然我们可以用矩阵快速幂求了,但是如果递推式的阶数比较大,那就很麻烦了,我们这里介绍一种神奇的优化。

特征向量与特征值

算法过程

模板题:[BZOJ 4161] Shlw loves matrixI

朴素版:

  1. #include <cstdio>
  2. #include <cstring>
  3. #define mset(a,x,n) memset(a,x,(n)*sizeof(*(a)))
  4. #define mcpy(a,b,n) memcpy(a,b,(n)*sizeof(*(a)))
  5. const int MAXN=2005,P=1e9+7;
  6. int a[MAXN],b[MAXN];
  7. int f[MAXN],g[MAXN],h[MAXN];
  8. inline void inc(int &x,int y){x+=y,x= x>=P ? x-P:x;}
  9. inline int pow(int x,int y){
  10. int ret=1;
  11. for(;y;x=(long long)x*x%P,y>>=1){
  12. if(y&1) ret=(long long)ret*x%P;
  13. }
  14. return ret;
  15. }
  16. inline int inv(int x){return pow(x,P-2);}
  17. inline void poly_mul(int a[],int b[],int p[],int m){
  18. static int c[MAXN<<1];
  19. mset(c,0,(m<<1)-1);
  20. for(int i=0;i<m;++i){
  21. for(int j=0;j<m;++j)
  22. inc(c[i+j],(long long)a[i]*b[j]%P);
  23. }
  24. for(int i=(m<<1)-2;i>=m;--i){
  25. int k=(long long)c[i]*inv(p[m])%P;
  26. for(int j=0;j<=m;++j)
  27. inc(c[i-j],P-(long long)k*p[m-j]%P);
  28. }
  29. mcpy(a,c,m);
  30. }
  31. inline void poly_pow(int x[],int y,int p[],int m,int ret[]){
  32. ret[0]=1;
  33. for(;y;poly_mul(x,x,p,m),y>>=1){
  34. if(y&1) poly_mul(ret,x,p,m);
  35. }
  36. }
  37. int main(){
  38. int n,k;
  39. scanf("%d%d",&n,&k);
  40. for(int i=1;i<=k;++i){
  41. scanf("%d",&a[i]);
  42. a[i]=(a[i]+P)%P;
  43. }
  44. for(int i=0;i<k;++i){
  45. scanf("%d",&b[i]);
  46. b[i]=(b[i]+P)%P;
  47. }
  48. if(n<k){
  49. printf("%d",b[n]);
  50. return 0;
  51. }
  52. for(int i=1;i<=k;++i) f[k-i]=P-a[i];
  53. f[k]=1;
  54. g[1]=1;
  55. poly_pow(g,n,f,k,h);
  56. int ans=0;
  57. for(int i=0;i<k;++i)
  58. inc(ans,(long long)h[i]*b[i]%P);
  59. printf("%d",ans);
  60. return 0;
  61. }

一个我一直不知道的东西

练习

[NOI 2017] 泳池

代码:

  1. #include <cstdio>
  2. #include <algorithm>
  3. #include <cstring>
  4. #define mset(a,x,n) memset(a,x,(n)*sizeof(*(a)))
  5. #define mcpy(a,b,n) memcpy(a,b,(n)*sizeof(*(a)))
  6. using std::min;
  7. const int MAXN=1005,P=998244353;
  8. int n,q,K;
  9. int powq[MAXN];
  10. int f[MAXN][MAXN],g[MAXN],h[MAXN<<1];
  11. int a[MAXN],b[MAXN];
  12. inline void inc(int &x,int y){x+=y,x= x>=P ? x-P:x;}
  13. inline int pow(int x,int y){
  14. int ret=1;
  15. for(;y;x=(long long)x*x%P,y>>=1){
  16. if(y&1) ret=(long long)ret*x%P;
  17. }
  18. return ret;
  19. }
  20. inline int inv(int x){return pow(x,P-2);}
  21. inline void DP(int K){
  22. mset(g,0,K+1);
  23. g[0]=1;
  24. for(int i=K;i;--i){
  25. int m=min(n,K/i);
  26. mset(f[i],0,m+1);
  27. for(int j=1;j<=m;++j){
  28. for(int k=1;k<=j;++k){
  29. inc(f[i][j],(long long)g[k-1]*g[j-k]%P);
  30. if(k!=K) inc(f[i][j],(long long)g[k-1]*f[i][j-k]%P);
  31. }
  32. f[i][j]=(long long)f[i][j]*(1-q+P)%P;
  33. }
  34. for(int j=1;j<=m;++j)
  35. g[j]=(long long)(g[j]+f[i][j])%P*powq[j]%P;
  36. }
  37. }
  38. inline void poly_mul(int a[],int b[],int ret[],int p[],int m){
  39. static int c[MAXN<<1];
  40. mset(c,0,(m<<1)-1);
  41. for(int i=0;i<m;++i){
  42. for(int j=0;j<m;++j)
  43. inc(c[i+j],(long long)a[i]*b[j]%P);
  44. }
  45. for(int i=(m<<1)-1;i>=m;--i){
  46. int k=(long long)c[i]*p[m]%P;
  47. for(int j=0;j<=m;++j)
  48. inc(c[i-j],P-(long long)k*p[m-j]%P);
  49. }
  50. mcpy(ret,c,m+1);
  51. }
  52. inline void poly_pow(int ret[],int y,int p[],int m){
  53. static int x[MAXN];
  54. mset(x,0,m+1),mset(ret,0,m);
  55. x[1]=1,ret[0]=1;
  56. for(;y;y>>=1){
  57. if(y&1) poly_mul(ret,x,ret,p,m);
  58. poly_mul(x,x,x,p,m);
  59. }
  60. }
  61. inline int cal(int K){
  62. if(K==0) return pow((1-q+P)%P,n);
  63. DP(K);
  64. mcpy(h,g,K+1);
  65. for(int i=K+1;i;--i)
  66. g[i]=(long long)(1-q+P)*g[i-1]%P;
  67. for(int i=1;i<=K;++i){
  68. for(int j=1;j<=i;++j)
  69. inc(h[i],(long long)h[i-j]*g[j]%P);
  70. }
  71. if(n<=K) return h[n];
  72. a[K+1]=1;
  73. for(int i=K+1;i;--i)
  74. a[K+1-i]=P-g[i];
  75. poly_pow(b,n,a,K+1);
  76. int ret=0;
  77. for(int i=0;i<=K;++i)
  78. inc(ret,(long long)b[i]*h[i]%P);
  79. return ret;
  80. }
  81. inline void pre_cal(){
  82. powq[0]=1;
  83. for(int i=1;i<=K;++i)
  84. powq[i]=(long long)powq[i-1]*q%P;
  85. }
  86. int main(){
  87. int x,y;
  88. scanf("%d%d%d%d",&n,&K,&x,&y);
  89. q=(long long)x*inv(y)%P;
  90. pre_cal();
  91. printf("%d",(cal(K)-cal(K-1)+P)%P);
  92. return 0;
  93. }
添加新批注
在作者公开此批注前,只有你和作者可见。
回复批注