[Shootout-list] mandelbrot errors (OCaml)

[Christophe TROESTLER]

----Next_Part(Mon_Feb_28_10_37_12_2005_436)--
Content-Type: Text/Plain; charset=us-ascii
Content-Transfer-Encoding: 7bit

On Sun, 27 Feb 2005, Christophe TROESTLER <debian00@tiscali.be> wrote:
> 
> (OCaml) After investigation, it seems that the pixel #300 difference
> may be due to a bug in the compiler

Well, after more investigation, it turns out we come back to the idea
of my original post about the instability of the computations.
Indeed, the problem is that some part of the computation was done in
extended precision (80 bits) while an intermediate result was stored
in double precision (64 bits) -- thanks to Jon Harrop for pointing
this to me.

Code that gives the same result as the C one (with both copilers) is
attached.

ChriS

----Next_Part(Mon_Feb_28_10_37_12_2005_436)--
Content-Type: Text/Plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
Content-Disposition: inline; filename="mandelbrot.ml"

(* http://shootout.alioth.debian.org/benchmark.php?test=mandelbrot&lang=all *)
(* Plot the Mandelbrot set [-1.5-i,0.5+i] on an N-by-N bitmap
   (http://www-info2.informatik.uni-wuerzburg.de/mitarbeiter/wolfram/lehre/bildformate.html#pbm). *)

let niter = 50
let limit = 2.

let limit2 = limit *. limit

let add_bit0 cr ci byte =
  let rec loop i zr zi =
    if zr *. zr +. zi *. zi > limit2 then (byte lsl 1) lor 0x00
    else if i > niter then (byte lsl 1) lor 0x01
    else loop (i + 1) (zr *. zr -. zi *. zi +. cr) (2. *. zr *. zi +. ci) in
  loop 1 0. 0.

let () =
  let w = int_of_string Sys.argv.(1) in
  let h = w in
  let fw = float w
  and fh = float h
  and cplmt8 = 8 - w mod 8 in
  Printf.printf "P4\n%i %i\n" w h;
  let byte = ref 0
  and bit = ref 0 in
  for y = 0 to h - 1 do
    let ci = 2. *. float y /. fh -. 1. in
    for x = 0 to w - 1 do
      byte := add_bit0 (2. *. float x /. fw -. 1.5) ci !byte;
      incr bit;
      if !bit = 8 then (
	output_byte stdout !byte;
	byte := 0;
	bit := 0;
      )
    done;
    if !bit <> 0 then (
      output_byte stdout (!byte lsl cplmt8);
      byte := 0;
      bit := 0;
    )
  done

----Next_Part(Mon_Feb_28_10_37_12_2005_436)----