Double-precision floating point arithmetics behaves differently in Guile and in C

[Panicz Maciej Godek <godek.maciek@gmail.com>]

Hi,

I recently ran into a very weird problem. I am using code for computing
Singular Value Decomposition that's been ported do Scheme by Gerald Sussman
as a part of the scmutils package (more specifically, I use some bits of
code extracted from guile-scmutils
<https://www.cs.rochester.edu/%7Egildea/guile-scmutils/>)

The code is too long to paste here, but I have uploaded a copy to my github
<https://github.com/panicz/slayer/blob/master/guile-modules/extra/scmutils/numerics/linear/svd.scm>
.

I have found an ellipsoid-fitting algorithm that leans on this SVD code:
when I run it on a randomly generated set of points lying on an ellipsoid,
the ellipsoid parameters are reconstructed correctly.

The trouble begins when I try to port that code to C. Initially, I have
used a C code for computing SVD that I've found elsewhere (which was based
on the same Fortran code from LINPACK), but I eventually wrote code to
convert Scheme code to C to make sure that both algorithms are identical
(and at this point I am quite confident that they are).

The problem is, that the results of floating-point operations behave
slightly differently in C than in Guile, and these small numerical
differences accumulate, and eventually start making huge differences -- at
least that's something that I've noticed. Moreover, I can quite confidently
say that the results that I get in C are incorrect -- that the *eigenvalues*
(which translate to the radii of the ellipsoid) are very different, the
calculated center of the ellipsoid is slightly different, and that the
*eigenvectors* (which can be interpreted as a rotation of the ellipsoid)
are very similar.

Either way, the ellipsoid isn't reconstructed correctly, because the points
used for reconstruction do not lie near that ellipsoid (unlike in the case
of the C code).

I am really troubled with this behavior -- according to the documentation,
and to the bits of source code that I managed to analyze, Guile uses
exactly the same representation of floating point numbers as the one I use
in C -- namely, the double type. (I only use single-precision floats to
express the input points, but I ran some tests using small integers on
inputs, to make sure, that they are exactly representable -- and the
results were still very different)

Is there a way to make sure that the floating point operations in C behave
identically as those in Guile?