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?
