# 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.

• 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 – JBernardo Jul 16 '11 at 18:14
• Can we output the number in hex? – Keith Randall Jul 18 '11 at 0:25
• @Keith Nope. That's part of the spec. – FUZxxl Jul 18 '11 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 '11 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 '11 at 17:07

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>

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
DBL
DBL
DBL
DBL
DBL
DBL
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 '11 at 20:37
• @J B: FIRST! :D – boothby Jul 19 '11 at 20:42
• Keep it ;) The next trick up my sleeve will benefit from Haskell more than from C. – J B Jul 19 '11 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 '11 at 22:10
• 3.6 secs on my machine. – FUZxxl Jul 19 '11 at 22:11

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. – FUZxxl Jul 19 '11 at 14:55
• For the record, I just put this up to be a smartass; you didn't say not to use builtins. – boothby Jul 19 '11 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 '11 at 19:47 • Gotcha. Should be enumFromStepN (s-1) instead of enumFromStepN s – J B Jul 19 '11 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. – FUZxxl Jul 19 '11 at 19:53 • @boothby No, I did not change the desired fibonacci number, but rather the reference output, which was wrong. – FUZxxl Jul 19 '11 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 '11 at 20:44 ## 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... – Timtech Jan 13 '14 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. – boothby Jul 18 '11 at 19:36 • @boothby Wrong. It writes to standard output if you use the command line interface. See for example stackoverflow.com/questions/6542537/… – Dr. belisarius Jul 18 '11 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. – boothby Jul 18 '11 at 20:36 • Sorry, can't reproduce since I don't have mathematica. – FUZxxl Jul 19 '11 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. ;-) – ESultanik Jul 25 '11 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 – ReyCharles Oct 12 '12 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"

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. – Arlen Aug 5 '11 at 3:39
• I ran this in ghci. I was pretty impressed. Haskell is great for these types of mathematical code problems. – Undreren Oct 12 '12 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);
System.out.println(f.toString());
}


# 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)
return res
}

func main() {
fmt.Println(fib_log(20000000))
}

• I tried parallelizing it (before adding the cache) using go routines and it started using 19 GiB of memory :/ – ReyCharles Oct 12 '12 at 10:32

pseudo code (I don't know what you guys are using)

product = 1
multiplier = 3 // 3 is fibonacci sequence, but this can be any number,
// generating an infinite amount of sequences
y = 28 // the 2^x-1 term, so 2^28-1=1,284,455,535th term
for (int i = 1; int < y; i++) {
product= sum*multiplier-1
multiplier= multiplier^2-2
}
multiplier=multiplier-product // 2^28+1 1,284,455,537th


It took my computer 56 hours to do those two terms. My computer is kind of crappy. I'll have the number in a text file on October 22nd. 1.2 gigs is a bit big to be sharing on my connection.

• I'm confused by your answer. Pseudocode? And yet you have timings? Post the code! Language doesn't matter! – boothby Oct 12 '12 at 5:06
• That, and the output is only supposed to be 4 million or so digits... – Wug Oct 12 '12 at 21:12