Bug 14803 - Different ULPs depending on size of long int in GCC
Summary: Different ULPs depending on size of long int in GCC
Status: RESOLVED FIXED
Alias: None
Product: glibc
Classification: Unclassified
Component: math (show other bugs)
Version: 2.17
: P2 normal
Target Milestone: ---
Assignee: Not yet assigned to anyone
URL:
Keywords:
Depends on:
Blocks:
 
Reported: 2012-11-04 22:24 UTC by H.J. Lu
Modified: 2014-06-14 11:16 UTC (History)
2 users (show)

See Also:
Host:
Target:
Build:
Last reconfirmed:
fweimer: security-


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description H.J. Lu 2012-11-04 22:24:14 UTC
Since GCC encodes float-pointing constants slightly different,
depend on size of long int:

http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55145

x86-64 glibc built by x32 GCC has

[hjl@gnu-tools-1 glibc-test]$ cat /export/build/gnu/glibc-test/build-x86_64-linx/math/test-ildoubl.out
testing long double (inline functions)
Failure: Test: asin (0x0.ffffffp0) == 1.5704510598101804156437184421571127056013
Result:
 is:          1.57045105981018041572e+00   0xc.9048a52ed37c86000000p-3
 should be:   1.57045105981018041561e+00   0xc.9048a52ed37c85f00000p-3
 difference:  1.08420217248550443401e-19   0x8.00000000000000000000p-66
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: asin (-0x0.ffffffp0) == -1.5704510598101804156437184421571127056013
Result:
 is:         -1.57045105981018041572e+00  -0xc.9048a52ed37c86000000p-3
 should be:  -1.57045105981018041561e+00  -0xc.9048a52ed37c85f00000p-3
 difference:  1.08420217248550443401e-19   0x8.00000000000000000000p-66
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: asin (0x0.ffffffff8p0) == 1.5707810680058339712015850710748035974710
Result:
 is:          1.57078106800583397126e+00   0xc.90f5aa22168bce000000p-3
 should be:   1.57078106800583397115e+00   0xc.90f5aa22168bcdf00000p-3
 difference:  1.08420217248550443401e-19   0x8.00000000000000000000p-66
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: asin (-0x0.ffffffff8p0) == -1.5707810680058339712015850710748035974710
Result:
 is:         -1.57078106800583397126e+00  -0xc.90f5aa22168bce000000p-3
 should be:  -1.57078106800583397115e+00  -0xc.90f5aa22168bcdf00000p-3
 difference:  1.08420217248550443401e-19   0x8.00000000000000000000p-66
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: asin_downward (-1.0) == -pi/2
Result:
 is:         -1.57079632679489661937e+00  -0xc.90fdaa22168c23600000p-3
 should be:  -1.57079632679489661926e+00  -0xc.90fdaa22168c23500000p-3
 difference:  1.08420217248550443400e-19   0x8.00000000000000000000p-66
 ulp       :  1.0000
 max.ulp   :  0.0000
Failure: Test: asin_upward (1.0) == pi/2
Result:
 is:          1.57079632679489661937e+00   0xc.90fdaa22168c23600000p-3
 should be:   1.57079632679489661926e+00   0xc.90fdaa22168c23500000p-3
 difference:  1.08420217248550443401e-19   0x8.00000000000000000000p-66
 ulp       :  1.0000
 max.ulp   :  0.0000

Test suite completed:
  5855 test cases plus 5008 tests for exception flags executed.
  6 errors occurred.
[hjl@gnu-tools-1 glibc-test]$ 

comparing against x86-64 glibc built with x86-64 GCC.
Comment 1 joseph@codesourcery.com 2012-11-04 22:42:15 UTC
On Sun, 4 Nov 2012, hjl.tools at gmail dot com wrote:

> Since GCC encodes float-pointing constants slightly different,
> depend on size of long int:

Which constants are involved here?  Ones in the asinl source code, or ones 
in the expected test results?  In either case we should regenerate the 
affected constant in hex float or with more decimal places, though I don't 
know offhand the right software for recomputing the rational function 
approximation for asinl.

> http://gcc.gnu.org/bugzilla/show_bug.cgi?id=55145

Really I'd consider this just a variant on bug 21718 (real.c rounding not 
perfect).  That would ideally be fixed by using MPFR for this in GCC ... 
except that for any MPFR versions before 3.1.1p2, the bug I found with the 
ternary value from mpfr_strtofr could cause problems for subnormal 
results.

But since GCC uses 160 bits internally, getting different results for a 
64-bit mantissa is a bit surprising unless one of the constants is 
extremely close to a half-way value.

> /export/build/gnu/glibc-test/build-x86_64-linx/math/test-ildoubl.out
> testing long double (inline functions)

We can of course also update libm-test-ulps to include these ulps, though 
preferably all computed constants in libm sources would be hex float to 
make unambiguous exactly what floating-point number is intended.
Comment 2 H.J. Lu 2012-11-04 22:54:25 UTC
(In reply to comment #1)

> But since GCC uses 160 bits internally, getting different results for a 
> 64-bit mantissa is a bit surprising unless one of the constants is 
> extremely close to a half-way value.

I believe it is the case here.

> > /export/build/gnu/glibc-test/build-x86_64-linx/math/test-ildoubl.out
> > testing long double (inline functions)
> 
> We can of course also update libm-test-ulps to include these ulps, though

Since libm-test-ulps can be auto-generated, we need to collect x86-64 ULPs
from both 32-bit long GCC and 64-bit long GCC.

> preferably all computed constants in libm sources would be hex float to 
> make unambiguous exactly what floating-point number is intended.

This would be ideal.
Comment 3 H.J. Lu 2012-11-04 23:02:52 UTC
(In reply to comment #1)
> On Sun, 4 Nov 2012, hjl.tools at gmail dot com wrote:
> 
> > Since GCC encodes float-pointing constants slightly different,
> > depend on size of long int:
> 
> Which constants are involved here?  Ones in the asinl source code, or ones 
> in the expected test results?  In either case we should regenerate the 
> affected constant in hex float or with more decimal places, though I don't 
> know offhand the right software for recomputing the rational function 
> approximation for asinl.

It is:

ieee754/ldbl-96/e_asinl.c: pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L,

32-bit long GCC encodes it as

1.57079632679489661925640447970309310221637133509

and 64-bit long GCC encodes it as

1.57079632679489661914798426245454265881562605500221252441

This difference causes the extra ULP.
Comment 4 H.J. Lu 2012-11-04 23:05:35 UTC
(In reply to comment #1)
>
> Really I'd consider this just a variant on bug 21718 (real.c rounding not 
> perfect).  That would ideally be fixed by using MPFR for this in GCC ... 
> except that for any MPFR versions before 3.1.1p2, the bug I found with the 
> ternary value from mpfr_strtofr could cause problems for subnormal 
> results.
> 

Fedora 17 has mpfr-3.1.0. So most of system MPFR won't be suitable for
this.
Comment 5 Vincent Lefèvre 2012-11-04 23:32:38 UTC
(In reply to comment #1)
> But since GCC uses 160 bits internally, getting different results for a 
> 64-bit mantissa is a bit surprising unless one of the constants is 
> extremely close to a half-way value.

The constant is *exactly* a halfway value (see below).

(In reply to comment #3)
> It is:
> 
> ieee754/ldbl-96/e_asinl.c: pio2_hi =
> 1.5707963267948966192021943710788178805159986950457096099853515625L,

In binary, this is the 65-bit number: 1.1001001000011111101101010100010001000010110100011000010001101001

So, for a long double with 64-bit precision, this is a halfway value. BTW, I suspect that this is not wanted, i.e. I would see this as a bug in the e_asinl.c code.

> 32-bit long GCC encodes it as
> 
> 1.57079632679489661925640447970309310221637133509

Exactly 1.570796326794896619256404479703093102216371335089206695556640625, or in binary: 1.100100100001111110110101010001000100001011010001100001000110101 (i.e. the constant rounded upward).

> and 64-bit long GCC encodes it as
> 
> 1.57079632679489661914798426245454265881562605500221252441

Exactly 1.5707963267948966191479842624545426588156260550022125244140625, or in binary: 1.100100100001111110110101010001000100001011010001100001000110100 (i.e. the constant rounded downward).
Comment 6 Florian Weimer 2012-11-05 09:33:53 UTC
(In reply to comment #1)

> We can of course also update libm-test-ulps to include these ulps, though 
> preferably all computed constants in libm sources would be hex float to 
> make unambiguous exactly what floating-point number is intended.

Using unsigned long long non-floating-point constants would be even safer, but more difficult to maintain.
Comment 7 Joseph Myers 2012-11-28 21:47:50 UTC
Fixed for 2.17 by:

commit 9984dd012623633f6e48d6d189a8df7d6c308535
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Wed Nov 28 21:46:16 2012 +0000

    Use hex float 64-bit values in ldbl-96 asinl (bug 14803).
Comment 8 Jackie Rosen 2014-02-16 19:31:10 UTC Comment hidden (spam)