线性递推乘积

首先有一定经验都应该能预见到两个线性递推数列的乘积还是个线性递推数列,然而用 BM 求解显得有些粗暴……

初级理解:注意当 a n a_n an N N N 阶递推, b n b_n bn M M M 阶递推,那么记 a n a_n an 的递推多项式为 ∏ i ( x − λ i ) s i \prod_i (x-\lambda_i)^{s_i} i(xλi)si b n b_n bn 的递推式为 ∏ j ( x − μ j ) t j \prod_j(x-\mu_j)^{t_j} j(xμj)tj,则可得
a n = ∑ i ∑ k = 0 s i u i k n k λ i n b n = ∑ j ∑ k = 0 t j v j k n k μ j n a n b n = ∑ i ∑ j ∑ k = 0 s i ∑ l = 0 t j u i k v j l n k + l λ i n μ j n \begin{aligned} a_n &= \sum_i \sum_{k=0}^{s_i} u_{ik} n^k \lambda_i^n\\ b_n &= \sum_j\sum_{k=0}^{t_j} v_{jk}n^k \mu_j^n \\ a_nb_n & = \sum_i \sum_j \sum_{k=0}^{s_i} \sum_{l=0}^{t_j}u_{ik} v_{j l} n^{k+l} \lambda_i^n \mu_j^n \end{aligned} anbnanbn=ik=0siuiknkλin=jk=0tjvjknkμjn=ijk=0sil=0tjuikvjlnk+lλinμjn
因此可知新的递推式的根是多重集 { λ i μ j } \{ \lambda_i \mu_j \} {λiμj} 组成的,其中包含的次数有 s i + t j + 1 ≤ ( s i + 1 ) ( t j + 1 ) s_i + t_j + 1\le (s_i + 1)(t_j + 1) si+tj+1(si+1)(tj+1) 所以为 ≤ N M \le NM NM 阶多项式。
但是,这个推导只在多项式总是可以完全分解的域上才可以。同时这个方法在实际得到递推式方面并不好用。

接下来的推导并不是严谨的,并没有考虑在一般域上的推导……但是足够显示出已有的结论。

我们记 f ( x ) , g ( x ) f(x), g(x) f(x),g(x) 的根的多重集 { α i } , { β j } \{\alpha_i\}, \{\beta_j\} {αi},{βj},我们在其上定义结式(resultant):
Res ⁡ x ( f ( x ) , g ( x ) ) = f N g M ∏ i = 1 N ∏ j = 1 M ( α i − β j ) \operatorname{Res}_x (f(x),g(x)) =f_N g_M\prod_{i=1}^N \prod_{j=1}^M (\alpha_i - \beta_j) Resx(f(x),g(x))=fNgMi=1Nj=1M(αiβj)
我们一会将指出通过等价的修改,可以在不一定可以多项式求根的环上定义结式。

现在让我们看看类似之处,我们要求一个 N M NM NM 阶多项式,它是
∏ i = 1 N ∏ j = 1 M ( x − α i β j ) \prod_{i=1}^N \prod_{j=1}^M (x - \alpha_i\beta_j) i=1Nj=1M(xαiβj)
它必然是一个线性递推,但不一定是最小的,但是对于任意多项式来说显然达到了一个上界。我们对其变形为
∏ i = 1 N α i M ∏ j = 1 M ( x / α i − β j ) = ∏ i = 1 N α i M g ( x / α i ) \begin{aligned} &\quad \prod_{i=1}^N \alpha_i^M \prod_{j=1}^M (x/\alpha_i - \beta_j)\\ &=\prod_{i=1}^N \alpha_i^M g(x/\alpha_i) \end{aligned} i=1NαiMj=1M(x/αiβj)=i=1NαiMg(x/αi)

注意对于结式,我们有变形:
∏ i = 1 N ∏ j = 1 M ( α i − β j ) = 1 g M ∏ i = 1 N g ( α i ) \prod_{i=1}^N \prod_{j=1}^M (\alpha_i - \beta_j) = \frac1{g_M}\prod_{i=1}^N g(\alpha_i) i=1Nj=1M(αiβj)=gM1i=1Ng(αi)
可以得到对于 ∏ i = 1 N α i M g ( x / α i ) \prod_{i=1}^N \alpha_i^M g(x/\alpha_i) i=1NαiMg(x/αi),令 t = α i t = \alpha_i t=αi,注意到 t M g ( x / t ) = ∑ j = 0 M t M − j g j x j t^M g(x/t) = \sum_{j=0}^M t^{M-j} g_j x^j tMg(x/t)=j=0MtMjgjxj,将 x x x 看做系数,则是关于 t t t 的多项式。 因此原式就是 R e s t ⁡ ( f ( t ) , t M g ( x / t ) ) \operatorname{Res_t} (f(t), t^Mg (x/t)) Rest(f(t),tMg(x/t))

解 1. Sylvester 矩阵

结式的计算可以用 Sylvester 矩阵的行列式表示:
Res ⁡ t ( f ( t ) , t M g ( x / t ) ) = ∣ 1 α 1 α 2 ⋯ α N 0 0 ⋯ 0 0 1 α 1 ⋯ α N − 1 α N 0 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 0 1 α 1 ⋯ α N β M β M − 1 x ⋯ β 1 x M − 1 x M 0 ⋯ 0 0 β M ⋯ β 2 x M − 1 β 1 x M − 1 x M 0 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 0 β M β M − 1 x ⋯ x M ∣ \operatorname{Res}_t(f(t), t^Mg(x/t)) = \left|\begin{array}{cccc} 1&\alpha_1&\alpha_2&\cdots&\alpha_N&0&0&\cdots&0\\ 0&1&\alpha_1&\cdots&\alpha_{N-1}&\alpha_N&0&\cdots&0\\ \vdots&&\ddots&&&&&&\vdots\\ 0&\cdots&0&1&\alpha_1&&&\cdots&\alpha_N\\ \beta_M&\beta_{M-1}x&\cdots&\beta_1 x^{M-1}&x^M&0&&\cdots&0\\ 0&\beta_M&\cdots&\beta_2 x^{M-1}&\beta_1 x^{M-1}&x^M&0&\cdots&0\\ \vdots&&\ddots&&&&&&\vdots\\ 0&\cdots&0&\beta_M&\beta_{M-1}x&&&\cdots&x^{M} \end{array}\right| Rest(f(t),tMg(x/t))=100βM00α11βM1xβMα2α1001β1xM1β2xM1βMαNαN1α1xMβ1xM1βM1x0αN0xM00000αN00xM
行列式可以在 Θ ( N ω ) \Theta(N^{\omega}) Θ(Nω) 时间内完成,即矩阵乘法复杂度。

鈤,这个东西插值的复杂度是 Θ ( N M ( N + M ) ω ) \Theta(NM(N+M)^{\omega}) Θ(NM(N+M)ω)

这个好像复杂度都不如 BM 把?

解 2. 多项式欧几里得

我们考虑研究结式的对称性。

考察式子 Res ⁡ x ( f ( x ) , g ( x ) ) = ∏ i = 1 N g ( α i ) \operatorname{Res}_x(f(x), g(x)) = \prod_{i=1}^N g(\alpha_i) Resx(f(x),g(x))=i=1Ng(αi),那么如果 α i \alpha_i αi 又不巧是 g g g 的一个根,那么 g ( α i ) = 0 g(\alpha_i) = 0 g(αi)=0

因此我们可以发现,对于 Res ⁡ x ( f ( x ) , g ( x ) ) \operatorname{Res}_x(f(x), g(x)) Resx(f(x),g(x)) 且那么任取 q ( x ) q(x) q(x),有 Res ⁡ x ( f ( x ) , g ( x ) ) = Res ⁡ x ( f ( x ) , g ( x ) + q ( x ) f ( x ) ) \operatorname{Res}_x (f(x), g(x)) = \operatorname{Res}_x (f(x), g(x) + q(x)f(x)) Resx(f(x),g(x))=Resx(f(x),g(x)+q(x)f(x))

因此我们可以得到 Res ⁡ x ( f ( x ) , g ( x ) ) = ( − 1 ) N M Res ⁡ x ( g , f   m o d   g ) \operatorname{Res}_x (f(x), g(x)) = (-1)^{NM}\operatorname{Res}_x (g, f\bmod g) Resx(f(x),g(x))=(1)NMResx(g,fmodg)。至此我们将结式计算归结与多项式欧几里得。多项式欧几里得的朴素复杂度为 Θ ( N M ) \Theta(NM) Θ(NM),因此朴素做法的复杂度为 Θ ( N 2 M 2 ) \Theta(N^2M^2) Θ(N2M2)

优化:Half-GCD 算法

事实上,我们可以通过分治方法快速进行多项式欧几里得,这一方法进行的多项式欧几里得可以做到 Θ ( ( N + M ) log ⁡ 2 ( N + M ) ) \Theta((N+M)\log^2 (N+M)) Θ((N+M)log2(N+M))。因此复杂度为 Θ ( N M ( N + M ) log ⁡ 2 ( N + M ) ) \Theta(NM(N+M)\log^2 (N+M)) Θ(NM(N+M)log2(N+M))

这看起来还是不见得跑得过 BM……

解 3. 另一个变形 //idea:djq_cpp

我们尝试翻转多项式,取对数:
ln ⁡ [ ∏ i = 1 N ∏ j = 1 M ( 1 − α i β j x ) ] = ∑ i = 1 N ∑ j = 1 M ∑ k ≥ 1 − ( α i β j x ) k k = − ∑ k ≥ 1 x k k ( ∑ i = 1 N α i k ) ( ∑ j = 1 M β j k ) = − ∑ k ≥ 1 x k k ( [ x k − 1 ] f ′ / f ) ( [ x k − 1 ] g ′ / g ) \begin{aligned} &\quad \ln \left[\prod _{i=1}^N \prod_{j=1}^M (1 - \alpha_i\beta_j x)\right]\\ &= \sum_{i=1}^N \sum_{j=1}^M \sum_{k\ge 1} - \frac{(\alpha_i\beta_jx)^k}k\\ &= -\sum_{k\ge 1}\frac{x^k}k \left(\sum_{i=1}^N \alpha_i^k \right) \left(\sum_{j=1}^M \beta_j^k\right) \\ &= -\sum_{k\ge 1} \frac{x^k}k ([x^{k-1}] f'/f) ([x^{k-1}]g'/g) \end{aligned} ln[i=1Nj=1M(1αiβjx)]=i=1Nj=1Mk1k(αiβjx)k=k1kxk(i=1Nαik)(j=1Mβjk)=k1kxk([xk1]f/f)([xk1]g/g)
复杂度瓶颈为做一次长 N M NM NM 的多项式 exp ⁡ \exp exp,复杂度为 Θ ( N M log ⁡ ( N M ) ) \Theta(NM\log (NM)) Θ(NMlog(NM))

参考资料

https://math.stackexchange.com/questions/1348838/sum-and-product-of-linear-recurrences

https://cp-algorithms.com/algebra/polynomial.html

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值