[bug#49659] [PATCH core-updates] gnu: guile: Fix failing tests on i686-linux.

["Ludovic Courtès" <ludo@gnu.org>]

Maxime Devos <maximedevos@telenet.be> skribis:

> Ludovic Courtès schreef op di 20-07-2021 om 15:55 [+0200]:

[...]

>>   1. Should we make it conditional on
>>        (or (string-prefix? "i686-" %host-type)
>>            (string-prefix? "i586-" %host-type))
>
> Rather, (target-x86-32?). target-x86-32? also recognises "i486-linux-gnu"
> even though that's not a ‘supported’ cross-target.

Yes, makes sense.

>>      ?  (I wonder why armhf-linux doesn’t have the same problem.)
>
> AFAIK floats & doubles on arm don't have excess precision.
>
> Floating-point numbers are either 32-bit or 64-bit,
> unlike in x86, where the floating-point registers are 80-bit
> but (sizeof) double==8 (64 bits).
>
> (I'm ignoring MMX, SSE and the like.)
>
> I don't know any architectures beside x86 which have excess precision.
> "-fexcess-precision=standard" should be harmless on architectures
> that don't have excess precision.
>
> I'd make it unconditional, but conditional on x86-target? should work
> for all ‘supported’ targets in Guix.

Alright.

I’d still err on the side of making the change only for target-x86-32?,
because that’s the only case where we know it’s needed.

>>   2. Is there any downside to compiling all of libguile with this flag?
>
> I searched with "git grep -F double" and "git grep -F float".
> Floating-point arithmetic doen't seem to be used much outside numbers.c.
>
> There's vm-engine.c, but the results of the calculations are written
> to some (stack?) memory (not a register), so the excess precision
> would be thrown away anyway, so I don't expect a slow-down.
>
> No code appears to be relying on excess precision.

OK.

>>   3. Do we have a clear explanation of why ‘-fexcess-precision=fast’
>>      (the default) would lead to failures in ‘numbers.test’?
>
> The problem I think is that the rounding choices made in
>   scm_i_inexact_floor_divide
> must be consistent with those made in
>   scm_i_inexact_floor_quotient
> and 
>   scm_i_inexact_floor_remainder
> (There are tests testing whether the results agree.)
>
> "-fexcess-precision=standard" reduces the degrees of freedom GCC has
> in choosing when to round, so it is more likely the rounding choices
> coincide?  It's only an untested hypothesis though.
>
> FWIW, I think this line:
>
>     /* in scm_i_inexact_floor_remainder */
>     return scm_i_from_double (x - y * floor (x / y));
>
> should be (for less fragility in case GCC changes its optimisations and
> register allocation / spilling)
>
>     /* in scm_i_inexact_floor_remainder */
>     return scm_i_from_double (x - y * (double) floor (x / y));
>
> even then, there's no guarantee the rounding choices for x/y
> are the same in scm_i_inexact_floor_divide, scm_i_inexact_floor_remainder
> and scm_i_inexact_floor_quotient.

Makes sense.  Seems to me that this should simply be implemented
differently to avoid the inconsistency in the first place (or one could
ignore IA32 altogether…).

> I dunno if 'floor' is broken.  Let me explain why this output is possible for a
> well-implemented 'floor':
>
> I want to note that printf("%f", floor(x/y))
> can display 16 different strings:
>
>   GCC can choose to round 'x' and/or 'y' by moving it from a register to stack memory.
>   (doesn't apply here I think because SCM values discard the excess precision)
>
>   GCC can choose to round the result of x/y (same reasons)
>
>   GCC can choose to round the result of floor(x/y) (same reasons)
>
> Likewise, printf("%f", x/y) can display 8 different strings, and the rounding
> choices may be different from those made for printf("%f", floor(x/y)).
>
> So this inconsistency (x/y=91.00... > 90.00 = floor(x/y))  doesn't really
> surprise me.  A more concrete scenario: suppose the CPU is configured to round
> upwards, and 'floor' is a mapping between extended-precision floating-point numbers.
>
> Let 'x' and 'y' be two floating-point numbers, such that:
>
>  (1) the mathematical value of 'x/y' is slightly less than 91
>  (2) 'x/y' when calculated in extended precision is slightly less than 91.
>      Denote this approximation as F1.
>  (3) 'x/y' when calculated in double precision is 91 (or slightly larger)
>      (due to the ‘rounding upwards’ mode, in ‘rounding downwards’ it might
>       have been slightly less than 91 as in (2))
>      Denote this approximation as F2.
>
> Then [floor(x/y)=] floor(F1)=floor(90.999...)=90.0,
> and  [x/y=] F2=91.0, thus we seemingly have x/y >= 1 + floor(x/y),
> even though that's mathematically nonsense.
>
> Thus, by choosing when to round (in-)appropriately, it is possible (on x86)
> that printf("x/y=%f, floor(x/y)=%f",x/y,floor(x/y)) will output "x/y=91.0 floor(x/y)=90.0".

I’m no expert but that makes sense to me.

Could you send an updated patch?

If you think of a way to fix the issue in Guile itself, we can also do
that.  :-)

Thanks for the investigation & explanation!

Ludo’.