[MLton] Re: [Haskell-cafe] fastest Fibonacci numbers in the West

William Lee Irwin III wli@holomorphy.com
Thu, 27 Jan 2005 07:29:16 -0800

Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline

On Thu, Jan 27, 2005 at 04:20:50PM +0100, Wesley W. Terpstra wrote:
> I can confirm this bug independently; attached is fib.sml.
> ./fib 99999999 segfaults as he mentioned.
> Since I doubt he used exactly the same algorithm, this is probably just a
> problem with the sheer size of the computed numbers. This program completes
> for 10^7-1 in 11 seconds (faster machine).
> In fact, my suspicion (unverified) is a bug in libgmp's conversion to base
> 10 during the output phase. This is a very hard thing to do well. My
> reasoning is that the program takes longer as you grow the number beyond 
> the crash point, but always segfaults before output.
> The program works for 
> 	./fib 19999999
> and creates a 4MB file, but crashes for
> 	./fib 25999999

My testcase is attached. It's equivalent to the matrix exponentiation
algorithm, except it uses bitreversal for the exponentiation, and does
the recurrence on A^n (0,1) for doubling steps explicitly to cut down
on the number of state variables (when these are large integers, they
may consume a great deal of space). It wasn't derived this way, but
rather, via Fibonacci number identities.

-- wli

Content-Type: text/plain; charset=us-ascii
Content-Description: fib.sml
Content-Disposition: attachment; filename="fib.sml"

fun	  unfold 0 = []
	| unfold n = let val ((q:int), (r:int)) = (n div 2, n mod 2)
		     in (unfold q) @ [if r = 1 then true else false]

fun crunch (p, (f, g)) = if p then (f*(f+2*g), f*f+g*g) else (f*f+g*g, g*(2*f-g))

fun myfib n = #2 (foldl crunch (1, 0) (unfold n))

val n = valOf(Int.fromString(hd(CommandLine.arguments())))
  handle Option =>
    (print ("usage: tailfib n\n");

val _ = print (IntInf.toString (myfib n) ^ "\n")