# Write the fastest Fibonacci

This is yet another challenge about the Fibonacci numbers.

The goal is to compute the 20'000'000th Fibonacii number as fast as possible. The decimal output is about 4 MiB large; it starts with:

28543982899108793710435526490684533031144309848579

The MD5 sum of the output is

fa831ff5dd57a830792d8ded4c24c2cb

You have to submit a program that calculates the number while running and puts the result to stdout. The fastest program, as measured on my own machine, wins.

Here are some additional rules:

• You have to submit the source code and a binary runnable on an x64 Linux
• The source code must be shorter than 1 MiB, in case of assembly it is also acceptable if only the binary is that small.
• You must not include the number to be computed in your binary, even in a disguised fashion. The number has to be calculated at runtime.
• My computer has two cores; you are allowed to use parallelism

I took a small implementation from the Internet which runs in about 4.5 seconds. It should not be very difficult to beat this, assuming that you have a good algorithm.

• Dude, anything like Sage that has indeterminate float precision will run that thing in less thant 1/10th of second. It's just a simple expression as phi = (1+sqrt(5))/2 Jul 16, 2011 at 18:14
• Can we output the number in hex? Jul 18, 2011 at 0:25
• @Keith Nope. That's part of the spec. Jul 18, 2011 at 10:16
• Since it's to be measured on your CPU, we might as well have some more information about it, couldn't we? Intel or AMD? Size of the L1 and instruction caches? Instruction set extensions?
– J B
Jul 19, 2011 at 13:37
• As I compute it, your start string and MD5 are for the 20'000'000th number, not the mere 2'000'000th.
– J B
Jul 19, 2011 at 17:07

Sage

Hmm, you seem to assume that the fastest is going to be a compiled program. No binary for you!

print fibonacci(2000000)


On my machine, it takes 0.10 cpu seconds, 0.15 wall seconds.

edit: timed on the console, instead of the notebook

• My idea was not to know, how fast your CAS can do this, but rather how fast you can code this by yourself. Jul 19, 2011 at 14:55
• For the record, I just put this up to be a smartass; you didn't say not to use builtins. Jul 19, 2011 at 20:19

This is my own try, although I did not wrote the algorithm by myself. I rather copied it from haskell.org and adapted it to use Data.Vector with its famous stream fusion:

import Data.Vector as V
import Data.Bits

main :: IO ()
main = print $fib 20000000 fib :: Int -> Integer fib n = snd . V.foldl' fib' (1,0) . V.dropWhile not$ V.map (testBit n) $V.enumFromStepN (s-1) (-1) s where s = bitSize n fib' (f,g) p | p = (f*(f+2*g),ss) | otherwise = (ss,g*(2*f-g)) where ss = f*f+g*g  This takes around 4.5 seconds when compiled with GHC 7.0.3 and the following flags: ghc -O3 -fllvm fib.hs  • Weird... I needed to change 20000000 to 40000000 to get it to print the expected number. – J B Jul 19, 2011 at 19:47 • Gotcha. Should be enumFromStepN (s-1) instead of enumFromStepN s – J B Jul 19, 2011 at 19:53 • @J B Sorry for all this confusion. I initially tested the program with different values to get a reasonably big number and saved the output into different files. But some how I confused them. I have updated the number to match the desired result. Jul 19, 2011 at 19:53 • @boothby No, I did not change the desired fibonacci number, but rather the reference output, which was wrong. Jul 19, 2011 at 20:07 • Side note: it's about 1.5s on my machine, but neither LLVM not Data.Vector seem to bring any significant advantage. – J B Jul 19, 2011 at 20:44 C with GMP, 3.6s Gods, but GMP makes code ugly. With a Karatsuba-style trick, I managed to cut down to 2 multiplies per doubling step. Now that I'm reading FUZxxl's solution, I'm not the first to have the idea. I've got a couple more tricks up my sleeve... maybe I'll try 'em out later on. #include <gmp.h> #include <stdio.h> #define DBL mpz_mul_2exp(u,a,1);mpz_mul_2exp(v,b,1);mpz_add(u,u,b);mpz_sub(v,a,v);mpz_mul(b,u,b);mpz_mul(a,v,a);mpz_add(a,b,a); #define ADD mpz_add(a,a,b);mpz_swap(a,b); int main(){ mpz_t a,b,u,v; mpz_init(a);mpz_set_ui(a,0); mpz_init(b);mpz_set_ui(b,1); mpz_init(u); mpz_init(v); DBL DBL DBL ADD DBL ADD DBL DBL DBL DBL ADD DBL DBL DBL ADD DBL DBL ADD DBL ADD DBL DBL ADD DBL DBL DBL DBL DBL DBL DBL DBL /*Comment this line out for F(10M)*/ mpz_out_str(stdout,10,b); printf("\n"); }  Built with gcc -O3 m.c -o m -lgmp. • LOL. Apart from a identifier naming, that's exactly my solution :) – J B Jul 19, 2011 at 20:37 • @J B: FIRST! :D Jul 19, 2011 at 20:42 • Keep it ;) The next trick up my sleeve will benefit from Haskell more than from C. – J B Jul 19, 2011 at 20:46 • First trick up my sleeve bumped in a GHC bug. Drat. I'll have to fall back to the second one, which isn't remotely as fun to implement, so it'll take time and motivation. – J B Jul 19, 2011 at 22:10 • 3.6 secs on my machine. Jul 19, 2011 at 22:11 ## COW  MoO moO MoO mOo MOO OOM MMM moO moO MMM mOo mOo moO MMM mOo MMM moO moO MOO MOo mOo MoO moO moo mOo mOo moo  Moo! (Takes a while. Drink some milk...) • Note: Although this really does work, it will probably never reach 20,000,000... Jan 13, 2014 at 22:19 Mathematica, interpreted: First@Timing[Fibonacci[2 10^6]]  Timed: 0.032 secs on my poor man's laptop.  And of course, no binary. • Doesn't print to stdout. Jul 18, 2011 at 19:36 • @boothby Wrong. It writes to standard output if you use the command line interface. See for example stackoverflow.com/questions/6542537/… Jul 18, 2011 at 19:44 • Nope, I'm using the commandline interface, version 6.0. Even using -batchoutput, it only prints timing info and not the Fibonacci number. Jul 18, 2011 at 20:36 • Sorry, can't reproduce since I don't have mathematica. Jul 19, 2011 at 14:49 • curl 'http://www.wolframalpha.com/input/?i=Fibonacci%5B2+10^6%5D' | grep 'Decimal approximation:' | sed ... It runs in constant time with respect to the speed of your Internet connection. ;-) Jul 25, 2011 at 17:06 # Ocaml, 0.856s on my laptop Requires the zarith library. I used Big_int but it's dog slow compared to zarith. It took 10 minutes with the same code! Most of the time was spent printing the damn number (9½ minutes or so)! module M = Map.Make (struct type t = int let compare = compare end) let double b = Z.shift_left b 1 let ( +. ) b1 b2 = Z.add b1 b2 let ( *. ) b1 b2 = Z.mul b1 b2 let cache = ref M.empty let rec fib_log n = if n = 0 then Z.zero else if n = 1 then Z.one else if n mod 2 = 0 then let f_n_half = fib_log_cached (n/2) and f_n_half_minus_one = fib_log_cached (n/2-1) in f_n_half *. (f_n_half +. double f_n_half_minus_one) else let f_n_half = fib_log_cached (n/2) and f_n_half_plus_one = fib_log_cached (n/2+1) in (f_n_half *. f_n_half) +. (f_n_half_plus_one *. f_n_half_plus_one) and fib_log_cached n = try M.find n !cache with Not_found -> let res = fib_log n in cache := M.add n res !cache; res let () = let res = fib_log 20_000_000 in Z.print res; print_newline ()  I can't believe how much a difference the library made! • For comparison @boothby's solution takes 0.875s to run on my laptop. It seems the difference is neglible. Also, apparently my laptop is fast :o Oct 12, 2012 at 13:28 ## Haskell On my system, this runs almost as fast as FUZxxl's answer (~18 seconds instead of ~17 seconds). main = print$ fst \$ fib2 20000000

-- | fib2: Compute (fib n, fib (n+1)).
--
-- Having two adjacent Fibonacci numbers lets us
-- traverse up or down the series efficiently.
fib2 :: Int -> (Integer, Integer)

-- Guard against negative n.
fib2 n | n < 0 = error "fib2: negative index"

-- Start with a few base cases.
fib2 0 = (0, 1)
fib2 1 = (1, 1)
fib2 2 = (1, 2)
fib2 3 = (2, 3)

-- For larger numbers, derive fib2 n from fib2 (n div 2)
-- This takes advantage of the following identity:
--
--    fib(n) = fib(k)*fib(n-k-1) + fib(k+1)*fib(n-k)
--             where n > k
--               and k ≥ 0.
--
fib2 n =
let (a, b) = fib2 (n div 2)
in if even n
then ((b-a)*a + a*b, a*a + b*b)
else (a*a + b*b, a*b + b*(a+b))

• Nice. I love Haskell. Aug 5, 2011 at 3:39
• I ran this in ghci. I was pretty impressed. Haskell is great for these types of mathematical code problems. Oct 12, 2012 at 8:03

## C, naive algorithm

Was curious, and I hadn't used gmp before... so:

#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

int main(int argc, char *argv[]){
int n = (argc>1)?atoi(argv):0;

mpz_t temp,prev,result;
mpz_init(temp);
mpz_init_set_ui(prev, 0);
mpz_init_set_ui(result, 1);

for(int i = 2; i <= n; i++) {
mpz_swap(temp, result);
mpz_swap(temp, prev);
}

printf("fib(%d) = %s\n", n, mpz_get_str (NULL, 10, result));

return 0;
}


fib(1 million) takes about 7secs... so this algorithm won't win the race.

I implemented the matrix multiplication method (from sicp, http://sicp.org.ua/sicp/Exercise1-19) in SBCL but it takes about 30 seconds to finish. I ported it to C using GMP, and it returns the correct result in about 1.36 seconds on my machine. It's about as fast as boothby's answer.

#include <gmp.h>
#include <stdio.h>

int main()
{
int n = 20000000;

mpz_t a, b, p, q, psq, qsq, twopq, bq, aq, ap, bp;
int count = n;

mpz_init_set_si(a, 1);
mpz_init_set_si(b, 0);
mpz_init_set_si(p, 0);
mpz_init_set_si(q, 1);
mpz_init(psq);
mpz_init(qsq);
mpz_init(twopq);
mpz_init(bq);
mpz_init(aq);
mpz_init(ap);
mpz_init(bp);

while(count > 0)
{
if ((count % 2) == 0)
{
mpz_mul(psq, p, p);
mpz_mul(qsq, q, q);
mpz_mul(twopq, p, q);
mpz_mul_si(twopq, twopq, 2);

mpz_add(p, psq, qsq);    // p -> (p * p) + (q * q)
mpz_add(q, twopq, qsq);  // q -> (2 * p * q) + (q * q)
count/=2;
}

else
{
mpz_mul(bq, b, q);
mpz_mul(aq, a, q);
mpz_mul(ap, a, p);
mpz_mul(bp, b, p);

mpz_add(a, bq, aq);      // a -> (b * q) + (a * q)
mpz_add(a, a, ap);       //              + (a * p)

mpz_add(b, bp, aq);      // b -> (b * p) + (a * q)

count--;
}

}

gmp_printf("%Zd\n", b);
return 0;
}


## Java: 8 seconds to compute, 18 seconds to write

public static BigInteger fibonacci1(int n) {
if (n < 0) explode("non-negative please");
short charPos = 32;
boolean[] buf = new boolean;
do {
buf[--charPos] = (n & 1) == 1;
n >>>= 1;
} while (n != 0);
BigInteger a = BigInteger.ZERO;
BigInteger b = BigInteger.ONE;
BigInteger temp;
do {
if (buf[charPos++]) {
a = temp;
} else {
a = a.multiply(b.shiftLeft(1).subtract(a));
b = temp;
}
} while (charPos < 32);
return a;
}

public static void main(String[] args) {
BigInteger f;
f = fibonacci1(20000000);
// about 8 seconds
System.out.println(f.toString());
// about 18 seconds
}


# Go

It's embarrasingly slow. On my computer it takes a little less than 3 minutes. It's only 120 recursive calls, though (after adding the cache). Note that this may use a lot of memory (like 1.4 GiB)!

package main

import (
"math/big"
"fmt"
)

var cache = make(map[int64] *big.Int)

func fib_log_cache(n int64) *big.Int {
if res, ok := cache[n]; ok {
return res
}
res := fib_log(n)
cache[n] = res
return res
}

func fib_log(n int64) *big.Int {
if n <= 1 {
return big.NewInt(n)
}

if n % 2 == 0 {
f_n_half := fib_log_cache(n/2)
f_n_half_minus_one := fib_log_cache(n/2-1)
res := new(big.Int).Lsh(f_n_half_minus_one, 1)
res.Mul(f_n_half, res)
return res
}
f_n_half := fib_log_cache(n/2)
f_n_half_plus_one := fib_log_cache(n/2+1)
res := new(big.Int).Mul(f_n_half_plus_one, f_n_half_plus_one)
tmp := new(big.Int).Mul(f_n_half, f_n_half)