[Bug math/13824] New: exp2l(small_integer) produces rounding errors

[Bug math/13824] New: exp2l(small_integer) produces rounding errors

[bruno at clisp dot org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13824

             Bug #: 13824
           Summary: exp2l(small_integer) produces rounding errors
           Product: glibc
           Version: 2.8
            Status: NEW
          Severity: minor
          Priority: P2
         Component: math
        AssignedTo: unassigned@sourceware.org
        ReportedBy: bruno@clisp.org
    Classification: Unclassified


On Linux/SPARC, with glibc 2.7, exp2l(3.0L) needlessly introduces rounding
errors:
The exact value is 8.0L, the actual result is slightly smaller than 8.0L.

But the main purpose of the function exp2l() is to produce exact values for
small integers. Otherwise I could use
expl(x*0.693147180559945309417232121458176568075L).

In the other glibc ports, exp2l(3.0L) is actually an exact 8.0L.

How to reproduce:
====================== foo.c =================
#define _GNU_SOURCE 1
#include <math.h>
#include <stdio.h>
long double x = 3.0L;
int main ()
{
  long double y = exp2l (x);
  printf ("y = %.38Lf = %LA\n", y, y);
  return 0;
}
==============================================
$ gcc -m32 -Wall foo.c -lm

$ ./a.out 
y = 7.99999999999999999999999999999999845926 =
0X1.FFFFFFFFFFFFFFFFFFFFFFFFFFFEP+2

$ gcc -m64 -Wall foo.c -lm

$ ./a.out 
y = 7.99999999999999999999999999999999845926 =
0X1.FFFFFFFFFFFFFFFFFFFFFFFFFFFEP+2

Expected output:

y = 8.00000000000000000000000000000000000000 = 0X8P+0

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13824] exp2l(small_integer) produces rounding errors

[bruno at clisp dot org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13824

Bruno Haible <bruno at clisp dot org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Target|                            |sparc64-unknown-linux-gnu
               Host|                            |sparc64-unknown-linux-gnu
              Build|                            |sparc64-unknown-linux-gnu

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13824] exp2l(small_integer) produces rounding errors

[joseph at codesourcery dot com]

http://sourceware.org/bugzilla/show_bug.cgi?id=13824

--- Comment #1 from joseph at codesourcery dot com <joseph at codesourcery dot com> 2012-03-09 00:19:42 UTC ---
As far as I can see, there are only "real" exp2l implementations for x86, 
x86_64 and m68k.  Any platform using ldbl-128 or ldbl-128ibm gets the 
fallback __ieee754_expl (M_LN2l * x).  That will have large errors for 
some inputs where the result is near LDBL_MAX, as well as the problem 
noted here for small integers.

The approach used for exp2 could be implemented for the two affected long 
double formats - but a simpler interim fix to the fallback implementation 
would be to make it separate the integer and fractional parts of the input 
and only use __ieee754_expl on M_LN2l * (fractional part); that should 
keep errors down to a few ulp.  Obviously testcases should be added to the 
testsuite that before such a fix did show large errors for large inputs.

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13824] exp2l(small_integer) produces rounding errors

[bruno at clisp dot org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13824

--- Comment #2 from Bruno Haible <bruno at clisp dot org> 2012-03-09 01:24:49 UTC ---
(In reply to comment #1)
> a simpler interim fix to the fallback implementation 
> would be to make it separate the integer and fractional parts of the input 
> and only use __ieee754_expl on M_LN2l * (fractional part)

Exactly. That's also the approach used by the portable gnulib implementation of
exp2l:
<http://git.savannah.gnu.org/gitweb/?p=gnulib.git;a=blob;f=lib/exp2l.c;hb=HEAD>

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13824] exp2l(small_integer) produces rounding errors

[bruno at clisp dot org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13824

--- Comment #3 from Bruno Haible <bruno at clisp dot org> 2012-03-11 14:41:23 UTC ---
The bug is also seen on Linux/PowerPC:
====================== foo.c =================
#define _GNU_SOURCE 1
#include <math.h>
#include <stdio.h>
long double x = 3.0L;
int main ()
{
  long double y = exp2l (x);
  printf ("y = %LA, y - 8 = %LA\n", y, y - 8.0L);
  return 0;
}
==============================================
$ gcc -m32 -Wall foo.c -lm

$ ./a.out
y = 0X1.0000000000000000000000000027P+3, y - 8 = 0X1.38DF91E931B46P-104

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13824] exp2l(small_integer) produces rounding errors

[jsm28 at gcc dot gnu.org]

http://sourceware.org/bugzilla/show_bug.cgi?id=13824

Joseph Myers <jsm28 at gcc dot gnu.org> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|NEW                         |RESOLVED
         Resolution|                            |FIXED

--- Comment #4 from Joseph Myers <jsm28 at gcc dot gnu.org> 2012-03-22 12:56:28 UTC ---
Fixed by:

commit 48e44791e4d4d755bf7a7dd083d87584dc4779e4
Author: Joseph Myers <joseph@codesourcery.com>
Date:   Thu Mar 22 12:55:19 2012 +0000

    Fix exp2l inaccuracy (bug 13824).

-- 
Configure bugmail: http://sourceware.org/bugzilla/userprefs.cgi?tab=email
------- You are receiving this mail because: -------
You are on the CC list for the bug.

[Bug math/13824] exp2l(small_integer) produces rounding errors

[fweimer at redhat dot com]

https://sourceware.org/bugzilla/show_bug.cgi?id=13824

Florian Weimer <fweimer at redhat dot com> changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
              Flags|                            |security-

-- 
You are receiving this mail because:
You are on the CC list for the bug.