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

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

Thanks for the feedback!
It turned out that SVD was fine and I was looking in the wrong place
-- in the C counterpart, I was feeding my ellipsoid algorithm with
uninitialized junk,
and the compiler didn't shout at me.

Sorry for the noise!

niedz., 26 wrz 2021 o 15:16 Eli Zaretskii <eliz@gnu.org> napisał(a):

> > 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
>