Re: strtod rounding wrt. long double

[Aurelien Jarno <aurelien at aurel32 dot net>]

On Sat, Jan 25, 2014 at 06:03:59PM +0000, Joseph S. Myers wrote:
> On Fri, 24 Jan 2014, David Miller wrote:
> 
> > While regenerating the ULPs for sparc, I noticed that I get failures
> > for tst-strtod-round of the form:
> > 
> > strtold (0x1p-16446) returned 0x0.000000000000000100000001p-16382 not 0x0.0000000000000001p-16382 (default rounding mode)
> > strtold (0x1p-16446) returned 0x0.000000000000000100000001p-16382 not 0x0.0000000000000001p-16382 (FE_DOWNWARD)
> > strtold (0x1p-16446) returned 0x0.000000000000000100000001p-16382 not 0x0.0000000000000001p-16382 (FE_TOWARDZERO)
> > strtold (0x1p-16446) returned 0x0.000000000000000100000001p-16382 not 0x0.0000000000000001p-16382 (FE_UPWARD)
> > strtold (-0x1p-16446) returned -0x0.000000000000000100000001p-16382 not -0x0.0000000000000001p-16382 (default rounding mode)
> > strtold (-0x1p-16446) returned -0x0.000000000000000100000001p-16382 not -0x0.0000000000000001p-16382 (FE_DOWNWARD)
> > strtold (-0x1p-16446) returned -0x0.000000000000000100000001p-16382 not -0x0.0000000000000001p-16382 (FE_TOWARDZERO)
> > strtold (-0x1p-16446) returned -0x0.000000000000000100000001p-16382 not -0x0.0000000000000001p-16382 (FE_UPWARD)
> > 
> > Does this look like a test error or is this more likely some sparc
> > specific bug in the long double routines?
> 
> That sounds like some sparc-specific bug in strtold or one of the 
> functions it calls (you don't say whether this appears for 32-bit, 64-bit 
> or both).  I'd suggest first seeing if the computed number as passed to 
> round_and_return is correct, then whether the rounded (i.e. just shifted 
> right in this case, as an exact subnormal value) value passed to 
> __mpn_construct_long_double is correct.

The problem is at actually in the generic code, or rather at the
interface level. The value is correct up to round_and_return, the broken
value is computed by the first call__mpn_rshift. My first guess was that
the sparc32 specific version of this function was broken, but I also get
the same behavior with the generic version.

In turns out that for this input value, a shift of 64 has to be done,
which calls __mpn_rshift with cnt = 0. Looking at the generic code in 
stdlib/rshift.c, one can find this interesting comment:

|   Argument constraints:
|   1. 0 < CNT < BITS_PER_MP_LIMB
|   2. If the result is to be written over the input, WP must be <= UP.

Even more interesting, the code to handle the case cnt = 0 exists but is
disabled...

The problem is therefore that both the sparc32 and the generic version
assumes that cnt can't be 0. When it is the case, we end up doing a
logical right shift of the same size than the register (32 bit there),
which is undefined in the C standard. On some architectures (I tested
s390 which is 32-bit big endian with 128-bit long double like sparc) a
0 is returned in that case, while on sparc the shift is done using only
the last 5 bits, leaving the value unchanged.

The correct way is probably to fix the call to __mpn_rshift in 
round_and_return, but a temporary solution is to  use the generic
version of rshift with the #if 0 changed into #if 1.

-- 
Aurelien Jarno                          GPG: 4096R/1DDD8C9B
aurelien@aurel32.net                 http://www.aurel32.net