首先有一定经验都应该能预见到两个线性递推数列的乘积还是个线性递推数列,然而用 BM 求解显得有些粗暴……
初级理解:注意当
an 是
N 阶递推,
bn 是
M 阶递推,那么记
an 的递推多项式为
∏i(x−λi)si,
bn 的递推式为
∏j(x−μj)tj,则可得
anbnanbn=i∑k=0∑siuiknkλin=j∑k=0∑tjvjknkμjn=i∑j∑k=0∑sil=0∑tjuikvjlnk+lλinμjn
因此可知新的递推式的根是多重集
{λiμj} 组成的,其中包含的次数有
si+tj+1≤(si+1)(tj+1) 所以为
≤NM 阶多项式。
但是,这个推导只在多项式总是可以完全分解的域上才可以。同时这个方法在实际得到递推式方面并不好用。
接下来的推导并不是严谨的,并没有考虑在一般域上的推导……但是足够显示出已有的结论。
我们记
f(x),g(x) 的根的多重集
{αi},{βj},我们在其上定义结式(resultant):
Resx(f(x),g(x))=fNgMi=1∏Nj=1∏M(αi−βj)
我们一会将指出通过等价的修改,可以在不一定可以多项式求根的环上定义结式。
现在让我们看看类似之处,我们要求一个
NM 阶多项式,它是
i=1∏Nj=1∏M(x−αiβj)
它必然是一个线性递推,但不一定是最小的,但是对于任意多项式来说显然达到了一个上界。我们对其变形为
i=1∏NαiMj=1∏M(x/αi−βj)=i=1∏NαiMg(x/αi)
注意对于结式,我们有变形:
i=1∏Nj=1∏M(αi−βj)=gM1i=1∏Ng(αi)
可以得到对于
∏i=1NαiMg(x/αi),令
t=αi,注意到
tMg(x/t)=∑j=0MtM−jgjxj,将
x 看做系数,则是关于
t 的多项式。 因此原式就是
Rest(f(t),tMg(x/t))。
结式的计算可以用 Sylvester 矩阵的行列式表示:
Rest(f(t),tMg(x/t))=∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣10⋮0βM0⋮0α11⋯βM−1xβM⋯α2α1⋱0⋯⋯⋱0⋯⋯1β1xM−1β2xM−1βMαNαN−1α1xMβ1xM−1βM−1x0αN0xM000⋯⋯⋯⋯⋯⋯00⋮αN00⋮xM∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣∣
行列式可以在
Θ(Nω) 时间内完成,即矩阵乘法复杂度。
鈤,这个东西插值的复杂度是
Θ(NM(N+M)ω)。
这个好像复杂度都不如 BM 把?
解 2. 多项式欧几里得
我们考虑研究结式的对称性。
考察式子
Resx(f(x),g(x))=i=1∏Ng(αi),那么如果
αi 又不巧是
g 的一个根,那么
g(αi)=0。
因此我们可以发现,对于
Resx(f(x),g(x)) 且那么任取
q(x),有
Resx(f(x),g(x))=Resx(f(x),g(x)+q(x)f(x))。
因此我们可以得到
Resx(f(x),g(x))=(−1)NMResx(g,fmodg)。至此我们将结式计算归结与多项式欧几里得。多项式欧几里得的朴素复杂度为
Θ(NM),因此朴素做法的复杂度为
Θ(N2M2)。
优化:Half-GCD 算法
事实上,我们可以通过分治方法快速进行多项式欧几里得,这一方法进行的多项式欧几里得可以做到
Θ((N+M)log2(N+M))。因此复杂度为
Θ(NM(N+M)log2(N+M))。
这看起来还是不见得跑得过 BM……
解 3. 另一个变形 //idea:djq_cpp
我们尝试翻转多项式,取对数:
ln[i=1∏Nj=1∏M(1−αiβjx)]=i=1∑Nj=1∑Mk≥1∑−k(αiβjx)k=−k≥1∑kxk(i=1∑Nαik)(j=1∑Mβjk)=−k≥1∑kxk([xk−1]f′/f)([xk−1]g′/g)
复杂度瓶颈为做一次长
NM 的多项式
exp,复杂度为
Θ(NMlog(NM))。
参考资料
https://math.stackexchange.com/questions/1348838/sum-and-product-of-linear-recurrences
https://cp-algorithms.com/algebra/polynomial.html