FIOC is SAF (Slow as Fuck)

While reading my SICP, I found a neat fibonacci function that has log n operations. I looks something like this.


(define (fib n)
  (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   (+ (* p p) (* q q))     ; compute p'
                   (+ (* q q) (* 2 p q))   ; compute q'
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

So I implemented it in Haskell, just for fun.


fib n = fib' 1 0 0 1 n
  where fib' a b p q 0 = b
        fib' a b p q n | even n =
          fib' a b (p * p + q * q) (q * q + 2 * p * q) (n `div` 2)
        fib' a b p q n =
          fib' (b * q + a * q + a * p) (b * p + a * q) p q (n - 1)

To calculate and print fib(1000000) and fib(10000000) (i.e. the millionth and ten-millionth fibonacci numbers) takes about 3.26 seconds with the Haskell version (when compiled with GHC), and 43.24 seconds with the Scheme version (when interpreted by plt-r5rs, cause I can't figure out how to get it to be r6rs compliant, nor can I figure out how to make the compiler accept r5rs code :( )

Just to have some more fun, I also implemented it iteratively in Python.


def fib(n):
    a = 1
    b = 0
    p = 0
    q = 1

    while n > 0:
        if n % 2 == 0:
            new_p = p * p + q * q
            new_q = q * q + 2 * p * q
            n = n / 2
            p = new_p
            q = new_q
        else:
            new_a = b * q + a * q + a * p
            new_b = b * p + a * q
            n = n - 1
            a = new_a
            b = new_b

    return b

This takes CPython a whopping 1926.49 seconds to calculate the same numbers. Holy fuck, /prog/, maybe you're right about FIOC after all.

Here are the other implementations I've written and their times.

Ruby (995.35 seconds)

def fib(n)
    a = 1
    b = 0
    p = 0
    q = 1

    while n > 0
        if n % 2 == 0:
            new_p = p * p + q * q
            new_q = q * q + 2 * p * q
            n = n / 2
            p = new_p
            q = new_q
        else
            new_a = b * q + a * q + a * p
            new_b = b * p + a * q
            n = n - 1
            a = new_a
            b = new_b
        end
    end
    b
end

C using GMP (3.24 seconds)

void fib(unsigned long n)
{
    mpz_t a, b, p, q, tmp1, tmp2, tmp3;

    mpz_init(a);
    mpz_init(b);
    mpz_init(p);
    mpz_init(q);
    mpz_init(tmp1);
    mpz_init(tmp2);
    mpz_init(tmp3);

    mpz_set_ui(a, 1);
    mpz_set_ui(b, 0);
    mpz_set_ui(p, 0);
    mpz_set_ui(q, 1);

    while(n > 0)
    {
        if(n % 2 == 0)
        {
            // p * p + q * q
            mpz_mul(tmp1, p, p);
            mpz_mul(tmp2, q, q);
            mpz_add(tmp3, tmp1, tmp2);

            // q * q + 2 * p * q;
            mpz_mul_ui(tmp1, p, 2);
            mpz_add(tmp1, tmp1, q);
            mpz_mul(tmp1, tmp1, q);

            mpz_set(p, tmp3);
            mpz_set(q, tmp1);
            n >>= 1;
        }
        else
        {
            // b * q + a * q + a * p;
            mpz_mul(tmp1, b, q);
            mpz_mul(tmp2, a, q);
            mpz_mul(tmp3, a, p);
            mpz_add(tmp3, tmp3, tmp2);
            mpz_add(tmp3, tmp3, tmp1);

            // b * p + a * q;
            mpz_mul(tmp1, b, p);
            mpz_mul(b, a, q);
            mpz_add(b, tmp1, b);

            mpz_set(a, tmp3);
            n--;
        }
    }

    gmp_printf("%Zd\n", b);

    mpz_clear(a);
    mpz_clear(b);
    mpz_clear(p);
    mpz_clear(q);
    mpz_clear(tmp1);
    mpz_clear(tmp2);
    mpz_clear(tmp3);
}

>>43
[u][o][b][i]^BBCODE[sub]FAILURE[/i][/b][/o][/u]