Page 201 : 专栏 更快地计算递归式

事实上,要求m项递推式的第n项的值可以不使用矩阵,而是使用初项的线性表示,通过快速幂在O(m ^ 2 * log n)的时间内求出答案。有兴趣的读者可以试着思考看看。

看到这里,我首先想到的是《SICP》里有类似的内容,在Page31,练习1.19。

下文第一段是一次“失败”的思考,我最后发现自己并没有给出O(m ^ 2 * log n)的算法,而是自己重新建立了对矩阵幂的理解。

经过一天的思考(在缺觉的头疼中),我终于可以给出问题的一个解答,见第二段

从问题本身推演到矩阵乘积(O(m ^ 3 * log n)

以斐波那契为例,首先有递推式

Fn = Fn-1 + Fn-2

注意,这里所讨论的递推式的一般结构不仅仅局限于

fn = a1 * fn-1 + a2 * fn-2 ...

可以是更广义的

Ai = a1 * Ai-1 + b1 * Bi-1 ...
Bi = a2 * Ai-1 + b2 * Bi-1 ...
...

回来,讨论斐波那契。

为了便于理解,变换问题的形式,对于Fn = Fn-1 + Fn-2,对应有变换f1接受一个元组(a, b)得到(c, d),使得c = a + b, d = a

例如,输入是(Fn-1, Fn-2),结果就是(Fn, Fn-1)

我们有了给出第i项值给出第i+1项值的变换f1。

f1跟原递推式不同的地方在于它是一个变换,有统一的输入和输出。

现在,我们可以定义类同于函数复合的操作,例如

f3(x) = f1(f1(f1(x)))

f3 = f1 · f1 · f1 = f1 ^ 3

于是,可以理解

fn = f1 ^ n

我们现在定义了变换f1的幂乘,但是还不知道如何实现快速幂。

关键是如何定义fi(i > 1)上的·(乘法)操作,即

fi+j = fi · fj

回头观察斐波那契的例子,定义f1的,实际上是

c = a + b
d = a

更明白一点,是

c = 1 * a + 1 * b
d = 1 * a + 0 * b

可以观察到一个固定的结构,而所有的fi都基于这一结构,构成一族变换。

对于任意fi,有唯一一组p,q,s,t,使

c = p * a + q * b
d = s * a + t * b

给定fi{pi, qi, si, ti}fj{pj, qj, sj, tj},类比函数复合,定义fi+j = fi · fj

c = pi * (pj * a + qj * b) + qi * (sj * a + tj * b)
d = si * (pj * a + qj * b) + ti * (sj * a + tj * b)

调整一下,

c = a * (pi * pj + qi * sj) + b * (pi * qj + qi * tj)
d = a * (si * pj + ti * sj) + b * (si * qj + ti * tj)

pi+j = pi * pj + qi * sj
qi+j = pi * qj + qi * tj
si+j = si * pj + ti * sj
ti+j = si * qj + ti * tj

这就是m=2时变换的一般结构。

在这里我注意到这个部分的复杂度是O(m ^ 3)(简单地说,此处m取2,有8个乘)。

实际上,上面的等式就是矩阵相乘的展开。

[pi+j, qi+j]        [pi, qi]        [pj, qj]
               =                *   
[si+j, ti+j]        [si, ti]        [sj, tj]

这样,可以理解矩阵其实是变换的一种非常自然的表现形式。

注意,因为我们的例子“斐波那契”的递推式中没有表现常数项,所以这里的结构是不包括常数项的。

O(m ^ 2)

经过思考,我发现自己扩大了问题域,即上段一开始给出的广义结构。

如果我们仅着眼于形如

Fn = a * Fn-1 + b * Fn-2 ...

这样的递推式,问题可以进一步简化。

上段所提出的模型的问题在于它维护了多条“轨道”,具体说,是m条轨道。

m=2时,p,q维护了一条轨道,s,t维护一条。

对于我们现在所关注的递推式,这种维护其实是不必要的。

m条轨道的意义在于我们必须维护m个值以支持变换操作。

但是

Fn = a * Fn-1 + b * Fn-2 ...

可以共用一条轨道:初始有m项初值和f1

Fm+1 = f1(F1 ... Fm)

继续直到

F2m-1 = f1(Fm-1 ... F2m-2)

现在,我们有F1 ... F2m-1

对于任意fi,都可以通过带入这m个m连续区间,得到

Fx+j = fi(Fj ... Fj+m)

即fi变换后又有m个值,可以支持下一步变换(fj),因为我们只关心变换的结构,我们这次不再扩展,只做一个值的变换

Fy = fj(Fx ... Fx+m)

接下来就可以像上一段一样计算参数了。

以斐波那契为例

对于斐波那契,有

//任意的初值组
a
b

//初值拓展所使用的常数,就是f1的结构
p1 = 1
q1 = 1

做一次(m-1次)拓展

c = p1 * a + q1 * b

fi{pi, qi}变换

x = pi * a + qi * b
y = pi * b + qi * c

fj{pj, qj}变换

z = pj * x + qj * y

化简

pj * x  = pj * (pi * a + qi * b)
        = pi * pj * a + qi * pj * b

qi * c  = qi * (p1 * a + q1 * b)
        = p1 * qi * a + q1 * qi * b

qj * y  = qj * (pi * b + qi * c)
        = qj * (pi * b + p1 * qi * a + q1 * qi * b)
        = p1 * qi * qj * a + pi * qj * b + q1 * qi * qj * b

z       = pi * pj * a + qi * pj * b
            + p1 * qi * qj * a + pi * qj * b + q1 * qi * qj * b
        = (pi * pj + p1 * qi * qj) * a
            + (2 * qi * pj + q1 * qi * qj) * b

即有

pi+j = pi * pj + p1 * qi * qj
qi+j = 2 * qi * pj + q1 * qi * qj

这就是m=2时

Fn = a * Fn-1 + b * Fn-2 ...

的一般变换组合结构

其中,p1和q1是常数,唯一对应于f1{p1, q1}

复杂度分析

这个变换显然比之前的矩阵形式有更少的乘积项,不过每一个乘积项乘的次数增多了(上例为7)

经过观察可以发现,2, p1, q1都是常数项,共有m ^ 2个乘积项,即只有不超过m ^ 2个常数项,不影响复杂度。

而每个乘积项中之多有2个变量出现(观察上式,计算x,y时引入一个,计算z时引入另一个,这里也解释为什么只有m ^ 2个乘积项),

所以单次变换复合的复杂度为O(m ^ 2)

注意,这里的变换都是以m的倍数为跨度的,比如f1通过1...m个初值得到第m+1项,f2则得到第2m + 1