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

[Eli Zaretskii <eliz@gnu.org>]

> From: Panicz Maciej Godek <godek.maciek@gmail.com>
> Date: Sun, 26 Sep 2021 13:55:20 +0200
> 
> I forgot to mention that I run Guile in an Emacs session running in a WSL
> console on Windows 10.
> The tests of the C code that I've been performing so far were executed in
> an MSYS terminal, but I have just tried running them in WSL console, and
> the radii get ridiculous values.
> 
> While the values used to generate the points are the following:
> 
> center '[-65.12 -50.54 88.66]:
> radii: '[83.95 47.13 45.56]
> rotation: '[[9.633871e-001   1.363818e-001   -2.308359e-001]
>                 [-1.734094e-001  9.735887e-001   -1.485064e-001]
>                 [2.044857e-001   1.830983e-001   9.615927e-001]]
> 
> and the values reconstructed in Guile are:
> 
>  (ellipsoid
>   #:center (65.10194623013226 -88.58460232582514 -50.721825868168324)
>   #:radii (83.94677717528019 45.56525864722978 47.12037877216948)
>   #:rotation ((-0.9633227088329351 0.2307406879308479 0.13699669185777974)
>                    (0.173638969876562 0.14674930301995573
> 0.9738142277679884)
>                    (-0.20459439578586436 -0.961885324244193
> 0.18143251146545056)))
> 
> (the signs and the order of radii differ, but it seems that the rotation
> matrix compensates for that)
> 
> in the case of MSYS (MinGW 64-bit) I get the following values:
> 
> approx_e.center = [90.13, -68.76, -42.46]
> approx_e.radius = [178.83, 97.08, 100.43]
> approx_e.rotation(3x3) =
>         -9.634121049254586e-001 1.734020884333972e-001
>  -2.043742444879810e-001
>         -2.314281253565038e-001 -1.535574849495593e-001
> 9.606566096217422e-001
>         -1.351966674037161e-001 -9.728061546592428e-001
> -1.880692600613582e-001
> 
> but in the case of the WSL console, I get
> 
> approx_e.center = [90.13, -68.76, -42.46]
> approx_e.radius = [1347844355973136335704668472606720.00,
> 731678657375689861259088782950400.00, 756961648396782369967223865344000.00]
> approx_e.rotation(3x3) =
>         -9.634121049254428e-01  1.734020884334219e-01
> -2.043742444880351e-01
>         -2.314281253566119e-01  -1.535574849499134e-01
>  9.606566096216593e-01
>         -1.351966674036449e-01  -9.728061546591822e-01
>  -1.880692600617213e-01

Forgive me for saying this, but frankly the most probable reason for
this is some bug in converting the Scheme code to C.  SVD is a
remarkably stable algorithm in the face of numerical roundoff errors,
that's one of its strong points.  So if the issue is with some minor
numerical inaccuracies, SVD results should be insensitive to it, as
long as the condition numbers of the matrices that you decompose are
not preposterously large, like 10^12 or something.  (Did you calculate
the condition numbers, btw? what are they?)

So I would check and recheck your C code, as IMO that's the first
suspect in this story.

Also, please note that on Intel CPUs double-precision FP calculations
usually keep intermediate results in 10-byte (80-bit) wide "long
double" format, unless your program is built with the GCC switch that
disables that (few programs do, so I don't expect you to do that,
eneither while building Guile, nor when compiling your C
implementation of SVD).  But that cannot by itself explain the
problem, because using 80-bit FP values should make the results more
accurate, not less accurate.

So again, I suspect your C program.  Does it work well when you throw
simple problems at it, like systems of linear equations where the
result is known in advance and the condition number of the matrix is
around 1?

HTH