矩阵快速幂的多项式取模优化

原问题

给定一个数列前k项,并给出其k阶递推关系 h n = ∑ i = 1 k a i h n − i h_n=\sum_{i=1}^k a_ih_{n-i} hn=i=1kaihni,求 h n h_n hn

矩阵快速幂

大家都会矩阵快速幂的方法。

构造一个转移矩阵 A A A

A = [ a 1 a 2 a 3 ⋯ a k 1 0 0 ⋯ 0 0 1 0 ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ 0 ] A=\left[ \begin{matrix} a_1 & a_2 & a_3 &\cdots & a_k\\ 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 0 \end{matrix} \right] A=a1100a2010a3000ak000

用前k项构造初始矩阵,也是一个k维向量 H 1 H_1 H1

H 1 = ( h k h k − 1 h k − 2 ⋮ h 1 ) H_1=\left( \begin{matrix} h_k\\ h_{k-1}\\ h_{k-2}\\ \vdots \\ h_1 \end{matrix} \right) H1=hkhk1hk2h1

即可得 A n − 1 H 1 = H n A^{n-1}H_1=H_n An1H1=Hn

时间复杂度 O ( k 3 log ⁡ n ) O(k^3\log n) O(k3logn)

k如果比较大怎么办?我会卡常,jki玄学卡常,wys优化

多项式取模优化

前置知识

1.矩阵的特征值

对于矩阵 A A A和一个n维向量 x x x,如果有一个值 λ \lambda λ满足 A x = λ x Ax=\lambda x Ax=λx,则我们称 λ \lambda λ为矩阵 A A A的特征值。

2.矩阵的特征多项式

简单的缩好像就是特征方程再移项?

我们把之前的那个方程移项得 ( A − λ E ) x = 0 (A-\lambda E)x=0 (AλE)x=0,显然若x有非零解,那么需要满足

∣ A − λ E ∣ = ∣ a 11 − λ a 12 a 13 ⋯ a 1 k a 21 a 22 − λ a 23 ⋯ a 2 k a 31 a 32 a 33 − λ ⋯ a 3 k ⋮ ⋮ ⋮ ⋱ ⋮ a k 1 a k 2 a k 3 ⋯ a k k − λ ∣ = 0 |A-\lambda E|= \left| \begin{matrix} a_{11}-\lambda & a_{12} & a_{13} & \cdots & a_{1k}\\ a_{21} & a_{22}-\lambda & a_{23} & \cdots & a_{2k}\\ a_{31} & a_{32} & a_{33}-\lambda & \cdots & a_{3k} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{k1} & a_{k2} & a_{k3} & \cdots & a_{kk}-\lambda \end{matrix} \right| =0 AλE=a11λa21a31ak1a12a22λa32ak2a13a23a33λak3a1ka2ka3kakkλ=0

这其实是一个一元k阶方程,即特征方程。我们记 f A ( x ) = ∣ A − x E ∣ f_A(x)=|A-xE| fA(x)=AxE,称矩阵A的特征多项式

3.Cayley-Hamilton定理

f A ( A ) = 0 f_A(A)=0 fA(A)=0。具体的证明可以看wikipedia,因为我也不会

~~为了方便记忆,~~有一个著名的伪证: f A ( A ) = ∣ A − A E ∣ = 0 f_A(A)=|A-AE|=0 fA(A)=AAE=0

4.拉普拉斯展开

对于行列式A,我们定义它的关于第i行第j列的余子式 M i j M_{ij} Mij为原行列式删除第i行和第j列之后剩下的k-1阶行列式。

拉普拉斯展开是一个有关行列式的值的等式,按第i行展开得到如下等式

∣ A ∣ = ∑ j = 1 k ( − 1 ) j − 1 a i j ∣ M i j ∣ |A|=\sum_{j=1}^k(-1)^{j-1}a_{ij}|M_{ij}| A=j=1k(1)j1aijMij

类似的,我们也可以对第j列展开

∣ A ∣ = ∑ i = 1 k ( − 1 ) i a i j ∣ M i j ∣ |A|=\sum_{i=1}^k(-1)^ia_{ij}|M_{ij}| A=i=1k(1)iaijMij

5.多项式取模

f ( x ) = g ( x ) p ( x ) + r ( x ) f(x)=g(x)p(x)+r(x) f(x)=g(x)p(x)+r(x)

它们满足 d e g ( r ) &lt; d e g ( g ) &lt; d e g ( f ) deg(r)&lt;deg(g)&lt;deg(f) deg(r)<deg(g)<deg(f)

优化

我们知道如果给矩阵的某一行全部元素全部变为相反数,那么矩阵的行列式也会变成相反数,那么我们可以做这样一个操作,并记作 g g g

g ( x ) = ( − 1 ) k f A ( x ) = ∣ x − a 1 − a 2 − a 3 ⋯ − a k − 1 x 0 ⋯ 0 0 − 1 x ⋯ 0 ⋮ ⋮ ⋮ ⋱ ⋮ 0 0 0 ⋯ x ∣ g(x)=(-1)^kf_A(x)=\left| \begin{matrix} x-a_1 &amp; -a_2 &amp; -a_3 &amp;\cdots &amp; -a_k\\ -1 &amp; x &amp; 0 &amp; \cdots &amp; 0\\ 0 &amp; -1 &amp; x &amp; \cdots &amp; 0 \\ \vdots &amp; \vdots &amp; \vdots &amp; \ddots &amp; \vdots\\ 0 &amp; 0 &amp; 0 &amp; \cdots &amp; x \end{matrix} \right| g(x)=(1)kfA(x)=xa1100a2x10a30x0ak00x

为了求 g ( x ) g(x) g(x)的表达式,我们可以按照第一行对其进行拉普拉斯展开:

g ( x ) = x k − ∑ i = 1 k a i x k − i g(x)=x^k-\sum_{i=1}^k a_ix^{k-i} g(x)=xki=1kaixki

由于我们要求的东西就是 A n A^n An,那不妨构造一个这样的式子:

x n = g ( x ) p ( x ) + r ( x ) x^n=g(x)p(x)+r(x) xn=g(x)p(x)+r(x)

x = A x=A x=A则有 A n = g ( A ) p ( A ) + r ( A ) A^n=g(A)p(A)+r(A) An=g(A)p(A)+r(A)

由之前的Cayley-Hamilton定理我们知道 g ( A ) = ( − 1 ) k f A ( A ) = 0 g(A)=(-1)^kf_A(A)=0 g(A)=(1)kfA(A)=0。那么我们就得到了 A n = r ( A ) A^n=r(A) An=r(A),而 d e g ( r ) &lt; k deg(r)&lt;k deg(r)<k。为了方便,我们不妨看作 d e g ( r ) = k − 1 deg(r)=k-1 deg(r)=k1

那么我们怎么求 r ( A ) r(A) r(A)呢?原来的取模板子显然不可取,因为n可能很大,数组都开不下。考虑用倍增,即

x 2 &VeryThinSpace; m o d &VeryThinSpace; g ( x ) = ( x 1 &VeryThinSpace; m o d &VeryThinSpace; g ( x ) ) 2 &VeryThinSpace; m o d &VeryThinSpace; g ( x ) , x 4 &VeryThinSpace; m o d &VeryThinSpace; g ( x ) = ( x 2 &VeryThinSpace; m o d &VeryThinSpace; g ( x ) ) 2 &VeryThinSpace; m o d &VeryThinSpace; g ( x ) x^2\bmod g(x)=(x^1\bmod g(x))^2\bmod g(x),x^4\bmod g(x)=(x^2\bmod g(x))^2\bmod g(x) x2modg(x)=(x1modg(x))2modg(x),x4modg(x)=(x2modg(x))2modg(x)

考虑用快速幂实现,这需要做两个操作:

  • 1.多项式乘法
    • 暴力 O ( k 2 ) O(k^2) O(k2)
    • FFT O ( k log ⁡ k ) O(k\log k) O(klogk)
  • 2.多项式取模
    • 暴力 O ( k 2 ) O(k^2) O(k2),跳过中间的零项 O ( k t ) O(kt) O(kt)
    • FFT O ( k log ⁡ k ) O(k\log k) O(klogk),限制多,有时用不了?

你问我怎么暴力取模?当然是模拟大除法啊

至此,我们已经得到了 A n = ∑ i = 0 k − 1 r i A i A^n=\sum_{i=0}^{k-1} r_iA^i An=i=0k1riAi,再在两边乘上初始矩阵 H 1 H_1 H1即得

H n + 1 = ∑ i = 0 k − 1 r i H i + 1 H_{n+1}=\sum_{i=0}^{k-1} r_iH_{i+1} Hn+1=i=0k1riHi+1

由于 H H H是一个n维向量,而且中间涉及到的数乘及加法的运算,对于每一行是独立的,直接把第k行上的运算提出来就得到

h n + 1 = ∑ i = 0 k − 1 r i h i + 1 h_{n+1}=\sum_{i=0}^{k-1}r_ih_{i+1} hn+1=i=0k1rihi+1

时间复杂度常数较大的 O ( k 2 log ⁡ n ) O(k^2\log n) O(k2logn),或者常数更大的 O ( k log ⁡ k log ⁡ n ) O(k\log k\log n) O(klogklogn),fft的话要注意精度问题,mod较大就可能有问题。

相关题目

看了这么久,然而貌似应用不多?

bzoj4161

hdu4914

NOI2017泳池

  • 1
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值