bug in floating point implementation ?

From: Shyjan Mahamud (mahamud@marr.ius.cs.cmu.edu)
Date: Thu Feb 25 1999 - 01:43:05 MET


Message-Id: <199902250043.BAA21786@nez-perce.inria.fr>
Date: Wed, 24 Feb 99 19:43:05 EST
From: Shyjan Mahamud <mahamud@marr.ius.cs.cmu.edu>
To: caml-list@inria.fr
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



This archive was generated by hypermail 2b29 : Sun Jan 02 2000 - 11:58:20 MET