Version française
Home     About     Download     Resources     Contact us    
Browse thread
bug in floating point implementation ?
[ Home ] [ Index: by date | by threads ]
[ Search: ]

[ Message by date: previous | next ] [ Message in thread: previous | next ] [ Thread: previous | next ]
Date: -- (:)
From: Shyjan Mahamud <mahamud@m...>
Subject: bug in floating point implementation ?
hi folks, 
i'm new to CAML, so i'm not sure if this is a known issue.
is there a known bug or difference in the implementation of floating
point arithmetic between the toplevel and that for the standalone
compiled version. i've run the same exact function from both the
toplevel and a standalone program. the numbers come out slightly
different.

from the toplevel -->

test sample 1 residue = 0.584302
0 : -0.964812
1 : 0.221230
2 : 0.006142
3 : 0.069594
4 : -0.059595
5 : -0.069939
6 : -0.001882
7 : 0.004956
8 : -0.046719
9 : 0.063383
10 : -0.012646
11 : 0.004402

from the stanalone -->

test sample 1 residue = 0.570106
0 : -0.952025
1 : 0.229649
2 : 0.023926
3 : 0.010044
4 : -0.036213
5 : -0.052544
6 : -0.005504
7 : 0.010667
8 : -0.054280
9 : 0.079775
10 : 0.079762
11 : -0.058593

all corresponding numbers are supposed to be exactly the same.
also, as can be seen the errors accumulate as u go from 0 to 11.
both call an SVD routine implemented in C (hence i did have to create my own
toplevel linking in the SVD routine, so this is *not* the standard
toplevel). i'm running a linux box, RH 5.2.

i could provide the code, but it is long and am not sure if it helps
anyway. if it helps, the following routine is what gives the values
above (except for the residue), excuse the imperative style, it was
just a transplant from Numerical Recipes :

let svbksb uwv b =
  let isnan (x : float) = x <> x in
  let (u, w, v) = uwv in
  let m = (Array.length u) in
  let n = (Array.length w) in
  let sum = ref 0.0 in
  let tmp = Array.make n 0.0 in
  let x = Array.make n 0.0 in
  for j = 0 to n - 1 do
    sum := 0.0;
    if ((not (isnan w.(j))) && (w.(j) != 0.0)) then 
      begin
      	for i = 0 to m - 1 do 
	  sum := !sum +. u.(i).(j) *. b.(i)
      	done;
    	sum := !sum /. w.(j)
      end;
    tmp.(j) <- !sum
  done;
  for j = 0 to n - 1 do
    sum := 0.0;
    for jj = 0 to n - 1 do
      sum := !sum +. v.(j).(jj) *. tmp.(jj)
    done;
    x.(j) <- !sum
  done;
  x;;

i'm baffled !!!

- shyjan mahamud
Carnegie Mellon University