标题：fibonacci 的进化
出处：Felix021
时间：Mon, 18 Mar 2013 01:06:15 +0000
作者：felix021
地址：https://www.felix021.com/blog/read.php?2111

内容：
最近在看SICP，抛弃旧的世界观和方法论压力很大，不过还是很有收获的，比如学习了个O(log(n))的fibonacci算法，大涨姿势啊。

想起前几天看到的某个笔试题，说是用最快的办法计算fibonacci数列的第n项。虽然我知道它是有个通项公式的，但是不适用于精确计算，因此写了个迭代的算法，自以为已经很好了，现在看了真是too simple sometimes naive了。

众所皆知 fibonacci 的定义是f(n) = f(n - 1) + f(n - 2); f(1) = f(2) = 1（从f(1)=1开始算起）
#Python版：
fibonacci = lambda n: 1 if n <= 2 else fibonacci(n - 1) + fibonacci(n - 2)

/* C版 */
int fibonacci(int n) { 
    return n <= 2 ? 1 : fibonacci(n - 1) + fibonacci(n - 2); 
}

;Scheme版
(define (fibonacci n) (if (<= n 2) 1 (+ (fibonacci (- n 1)) (fibonacci (- n 2)))))

不幸的是这样树形展开效率太低了，当n=42的时候，C语言需要的时间已经超过1s了。

因此需要改成迭代版，使用 (a, b) <- (a + b, a) 这样的方式。
#Python版
def fibonacci(n):
    a = 1; b = 0;
    for i in range(n):
        a, b = a + b, a
    return b

/* C版 */
int fibonacci(int n) {
    int a = 1, b = 0, c;
    while (n--) {
        c = a;
        a = a + b;
        b = c;
    }
    return b;
}

;Scheme版
(define (fibonacci n)
  (define (iter n a b)
    (if (= n 0) b (iter (- n 1) (+ a b) a)))
  (iter n 1 0)) 

虽然O(n)的效率已经有显著的提升，但是由于这个数列的增长超过了2^n，所以当对于稍大的n，就需要使用大整数的运算，效率很低。python版的代码，当n=300,000的时候，需要超过1s才能得出结果；Scheme版则需要大约9s。

SICP的练习1-19里面则提到一个O(log(n))的巧妙算法：将计算fibonacci的每次迭代 (a, b) <- (a + b, a) 表示为一个变换T[p=0, q=1]，具体表示为（似乎是用矩阵乘法倒推过来的）
引用
T[pq](a, b) = (a(p+q) + bq, aq + bp)



通过计算 T[pq](T[pq](a, b)) （即对(a, b)进行两次T[pq]变换），可以得到
引用
(a((pp+qq) + (2pq+qq)) + b(2pq+qq), a(2pq+qq) + b(pp+qq))

记 p' = pp + qq, q' = 2pq+qq 则有
T[p'q'](a, b) = T[pq](T[pq](a, b)) = (a(p' + q') + bq', aq' + bp')

于是
f(n+1) = T[pq]n(f(1))



也就是说——可以通过类似分治计算 《a的n次方模b》 的算法来计算fibonacci了！

给出了算法以后，代码就容易写了：
#Python版
def fibonacci(n):
    def iter(a, b, p, q, n): 
        if n == 0:
            return b
        elif n % 2 == 0:
            return iter(a, b, p * p + q * q, 2 * p * q + q * q, n / 2)
        else:
            return iter(a * (p + q) + b * q, a * q + b * p, p, q, n - 1)
    return iter(1, 0, 0, 1, n)

/* C版 */
typedef unsigned long long ull;
ull fibo_iter(ull a, ull b, ull p, ull q, int n) {
    if (n == 0)
        return b;
    else if (n % 2 == 0)
        return fibo_iter(a, b, p *p + q * q, 2 * p * q + q * q, n / 2);
    else 
        return fibo_iter(a * (p + q) + b * q, a * q + b * p, p, q, n - 1);
}

ull fibonacci(int n) {
    return fibo_iter(1, 0, 0, 1, n);
}

/* Scheme版 */
(define (even? n) (= (remainder n 2) 0))
(define (fibonacci n)
    (define (iter a b p q n)
        (cond ((= n 0) b)
              ((even? n)
                (iter 
                    a
                    b
                    (+ (* p p) (* q q))
                    (+ (* 2 p q) (* q q))
                    (/ n 2)))
              (else
                (iter
                    (+ (* a (+ p q)) (* b q))
                    (+ (* a q) (* b p))
                    p
                    q
                    (- n 1)))))
    (iter 1 0 0 1 n))

经过这样的优化以后，效率显著，对于n=300,000，python版仅需要不到0.04s，而Scheme版也仅需要0.12s即可得出结果。

p.s. 以上代码均经过测试。


Generated by Bo-blog 2.1.0