bug-gnulib
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Note about gnulib's roundl/expl under NetBSD (limited precision leading


From: Assaf Gordon
Subject: Note about gnulib's roundl/expl under NetBSD (limited precision leading to segfault)
Date: Fri, 03 Jul 2015 00:37:08 -0400
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:31.0) Gecko/20100101 Thunderbird/31.7.0

Hello,

This is not a bug, but it took me a long time to troubleshoot, so I figured 
I'll share it (perhaps save someone else some time).

In short: under NetBSD, because of defaulting to 'double' precision (instead of 'long 
double', as nicely documented in gunlib's fpucw.h file), the "roundl" 
implementation might return incorrect values (mostly zero), and as a result, the 'expl' 
implementation can segfault.
The solution is to use  DECL_LONG_DOUBLE_ROUNDING/BEGIN_LONG_DOUBLE_ROUNDING 
from fpucw.h .

How it happens:
1. Using NetBSD 6.1.4 x86_64 (under QEMU VM).

2. 'roundl' and 'expl' are not provided, using gnulib's implementation instead.

3. The detected implementation for round is  FLOOR_FREE_ROUND (without ceil() 
or floor()),
    from lib/round.c lines 112-174.

4. In round.c line 115, the value of TWO_MANT_DIG is very large
    (  9223372036854775808.000000 ).

5. In round.c lines 146,147 and 165,166, the variable z contains
    the input value to round (+/-0.5), and the following code is used:

   /* Round to the next integer (nearest or up or down, doesn't
      matter).  */
   z -= TWO_MANT_DIG;
   z += TWO_MANT_DIG;

If z is negative and small (e.g. -31), and rounding mode is only 'double' 
precision,
the subtraction will make Z equal exactly TWO_MANT_DIG.
Then the addition will make Z equal zero.


6. If roundl() was called through expl(), the following code is executed 
(abbreviated from expl.c, lines 113-144):

     long double nm = roundl (x * LOG2_BY_256_INVERSE); /* = 256 * n + m */
     ...
     int n = (int) roundl (nm * (1.0L / 256.0L));
     int m = (int) nm - 256 * n;
     return ldexpl (gl_expl_table[128 + m] * exp_y, n);

Assuming expl was called with 'x=-20',
it incorrectly sets 'nm = -8191' (due to rounding errors in roundl),
and incorrectly sets 'n = 0',
then sets m to '-8191', which is outside the implicit assumed range of +/-128,
then it segfaults (access to invalid index in the array of 257 elements).

===

There is nothing to fix in the code, but perhaps adding a comment mentioning 
'fpucw.h' in 'roundl' or 'expl' would be helpful.

regards,
 -assaf



reply via email to

[Prev in Thread] Current Thread [Next in Thread]