Ocaml sums the harmonic series -- four ways, four benchmarks: floating point performance

["Will M. Farr" <farr@MIT.EDU>]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I've been looking at using ocaml to implement a gravitational n-body  
code, and therefore have quite a bit of interest in its floating-point  
performance.  Also, I'm learning the language by playing around with  
simple programs.  Here's an implementation (really 4) along with timing  
information of the "harmonic" benchmark (essentially summing the  
harmonic series), which can be found here:

http://shootout.alioth.debian.org/sandbox/benchmark.php? 
test=harmonic&lang=all&sort=cpu

After testing different ways of implementing the ocaml harmonic  
benchmark, I have settled on the following program.  For sizes of 1 000  
000 000 terms, it takes about 25% longer than the corresponding  
algorithm in c (note that I have replaced an int->float conversion for  
each term with a single floating point operation: ifloat := ifloat +.  
1.0).  Since int->float conversions are slow on my machine (PowerBook  
G4), this is a big win (about a factor of 2 in time for the C program).  
  Alas, when the number of terms approaches 16 digits, this method will  
lose accuracy, since <~16-digit number> +. 1.0 = <16-digit number +  
difference in last bit of mantissa>.  However, for sizes like the  
shootout seems to be using, this algorithm works fine (and the usual  
int type won't hold anything close to 16 digits anyway!).  I'm cc-ing  
this to the caml list because there may be people there interested in  
the floating point performance of Caml

Here's the code for the fastest implementation:

let sum_harmonic4 n =
   let sum = ref 1.0 in
   let ifloat = ref 2.0 in
     for i = 2 to n do
       sum := !sum +. 1.0/.(!ifloat);
       ifloat := !ifloat +. 1.0
     done;
     !sum;;

let _ =
   let n = int_of_string (Sys.argv.(1)) in
     Printf.printf "%g\n" (sum_harmonic4 n);;

And here's all the implementations I tried (for those interested in  
such things with ocaml):

let sum_harmonic n =
   let rec loop i sum =
     if i <= n then
       loop (i + 1) (sum +. 1.0/.(float_of_int i))
     else
       sum in
     loop 2 1.0;;

let sum_harmonic2 n =
   let sum = ref 1.0 in
   for i = 2 to n do
     sum := !sum +. 1.0/.(float_of_int i)
   done;
     !sum;;

let sum_harmonic3 n =
   let rec loop i ifloat sum =
     if i <= n then
       loop (i + 1) (ifloat +. 1.0) (sum +. 1.0/.ifloat)
     else
       sum in
     loop 2 2.0 1.0;;

let sum_harmonic4 n =
   let sum = ref 1.0 in
   let ifloat = ref 2.0 in
     for i = 2 to n do
       sum := !sum +. 1.0/.(!ifloat);
       ifloat := !ifloat +. 1.0
     done;
     !sum;;

let _ =
   let n = int_of_string (Sys.argv.(1)) in
     Printf.printf "%g\n" (sum_harmonic4 n);;

The timings for my machine (PowerBook G4, 800 Mhz) are as follows:

time ./harmonic 1000000000:
harmonic: 	user    2m1.590s
			sys     0m0.790s

harmonic2: 	user    2m0.340s
			sys     0m0.440s

harmonic3: 	user    1m44.350s
			sys     0m0.740s

harmonic4: 	user    1m12.680s
			sys     0m0.430s

Each invocation was compiled with "ocamlopt -unsafe -noassert -o  
harmonic harmonic.ml".  It looks like using references and loops is *by  
far* the fastest (and also that my PowerBook is pretty slow to convert  
int->float, but I don't think this is related to ocaml, since the C  
version does the same thing).

Hope you all find this interesting.

Will
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.0 (Darwin)

iD8DBQFB5pl3jFCrhUweU3MRApDzAJ9Ysln/KTQcq4WzxT9060GcDAgKQwCfTsb0
mDm4UyyghIz7m7r4ZpGcI3o=
=dLDI
-----END PGP SIGNATURE-----