Mar 18

fibonacci 的进化 不指定

felix021 @ 2013-3-18 01:06 [IT » 程序设计] 评论(4) , 引用(0) , 阅读(17978) | Via 本站原创 | |
最近在看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. 以上代码均经过测试。

转载请注明出自 ,如是转载文则注明原出处,谢谢:)
RSS订阅地址: http://www.felix021.com/blog/feed.php
xiami Email
2013-4-27 21:58
原来矩阵乘法不慢的 0////0faint\n上个暑假学线代什么的学了这个~当时直接用矩阵乘法低效把算法给毙了……\ntoo young too simple...unhappy  喂喂换行符怎么被吃了……
ctqmumu
2013-3-21 11:48
说起来真是惭愧... 虽然大四就买了SICP,但是到现在也只看完第一章,我看书还是太慢了 cry
felix021 回复于 2013-3-22 00:45
这个跟慢没关系吧。。。
snoopy
2013-3-18 15:36
跟着网易公开课看, 至少前 10 节还都是挺简单的 :P
felix021 回复于 2013-3-18 16:31
看视频好慢的感觉。。。
snoopy
2013-3-18 08:32
算法导论公开课的第一节还是第二节就讲了这个方法, 用 2*2 矩阵的 N 阶方来实现, 然后矩阵 N 阶方是可以做到 log(n) 的. 似乎以前集训队队内训练时有过这样的题
felix021 回复于 2013-3-18 09:33
我的算法一直很半吊子的来着。。。算导一直没有鼓起勇气去看。
分页: 1/1 第一页 1 最后页
发表评论
表情
emotemotemotemotemot
emotemotemotemotemot
emotemotemotemotemot
emotemotemotemotemot
emotemotemotemotemot
打开HTML
打开UBB
打开表情
隐藏
记住我
昵称   密码   *非必须
网址   电邮   [注册]