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
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.
(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>
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
Fixed by: commit 48e44791e4d4d755bf7a7dd083d87584dc4779e4 Author: Joseph Myers <joseph@codesourcery.com> Date: Thu Mar 22 12:55:19 2012 +0000 Fix exp2l inaccuracy (bug 13824).