Back to index

glibc  2.9
e_exp2.c
Go to the documentation of this file.
00001 /* Double-precision floating point 2^x.
00002    Copyright (C) 1997,1998,2000,2001,2005,2006 Free Software Foundation, Inc.
00003    This file is part of the GNU C Library.
00004    Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
00005 
00006    The GNU C Library is free software; you can redistribute it and/or
00007    modify it under the terms of the GNU Lesser General Public
00008    License as published by the Free Software Foundation; either
00009    version 2.1 of the License, or (at your option) any later version.
00010 
00011    The GNU C Library is distributed in the hope that it will be useful,
00012    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014    Lesser General Public License for more details.
00015 
00016    You should have received a copy of the GNU Lesser General Public
00017    License along with the GNU C Library; if not, write to the Free
00018    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
00019    02111-1307 USA.  */
00020 
00021 /* The basic design here is from
00022    Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical
00023    Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft.,
00024    17 (1), March 1991, pp. 26-45.
00025    It has been slightly modified to compute 2^x instead of e^x.
00026    */
00027 #ifndef _GNU_SOURCE
00028 #define _GNU_SOURCE
00029 #endif
00030 #include <stdlib.h>
00031 #include <float.h>
00032 #include <ieee754.h>
00033 #include <math.h>
00034 #include <fenv.h>
00035 #include <inttypes.h>
00036 #include <math_private.h>
00037 
00038 #include "t_exp2.h"
00039 
00040 /* XXX I know the assembler generates a warning about incorrect section
00041    attributes. But without the attribute here the compiler places the
00042    constants in the .data section.  Ideally the constant is placed in
00043    .rodata.cst8 so that it can be merged, but gcc sucks, it ICEs when
00044    we try to force this section on it.  --drepper  */
00045 static const volatile double TWO1023 = 8.988465674311579539e+307;
00046 static const volatile double TWOM1000 = 9.3326361850321887899e-302;
00047 
00048 double
00049 __ieee754_exp2 (double x)
00050 {
00051   static const double himark = (double) DBL_MAX_EXP;
00052   static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1);
00053 
00054   /* Check for usual case.  */
00055   if (isless (x, himark) && isgreaterequal (x, lomark))
00056     {
00057       static const double THREEp42 = 13194139533312.0;
00058       int tval, unsafe;
00059       double rx, x22, result;
00060       union ieee754_double ex2_u, scale_u;
00061       fenv_t oldenv;
00062 
00063       feholdexcept (&oldenv);
00064 #ifdef FE_TONEAREST
00065       /* If we don't have this, it's too bad.  */
00066       fesetround (FE_TONEAREST);
00067 #endif
00068 
00069       /* 1. Argument reduction.
00070         Choose integers ex, -256 <= t < 256, and some real
00071         -1/1024 <= x1 <= 1024 so that
00072         x = ex + t/512 + x1.
00073 
00074         First, calculate rx = ex + t/512.  */
00075       rx = x + THREEp42;
00076       rx -= THREEp42;
00077       x -= rx;  /* Compute x=x1. */
00078       /* Compute tval = (ex*512 + t)+256.
00079         Now, t = (tval mod 512)-256 and ex=tval/512  [that's mod, NOT %; and
00080         /-round-to-nearest not the usual c integer /].  */
00081       tval = (int) (rx * 512.0 + 256.0);
00082 
00083       /* 2. Adjust for accurate table entry.
00084         Find e so that
00085         x = ex + t/512 + e + x2
00086         where -1e6 < e < 1e6, and
00087         (double)(2^(t/512+e))
00088         is accurate to one part in 2^-64.  */
00089 
00090       /* 'tval & 511' is the same as 'tval%512' except that it's always
00091         positive.
00092         Compute x = x2.  */
00093       x -= exp2_deltatable[tval & 511];
00094 
00095       /* 3. Compute ex2 = 2^(t/512+e+ex).  */
00096       ex2_u.d = exp2_accuratetable[tval & 511];
00097       tval >>= 9;
00098       unsafe = abs(tval) >= -DBL_MIN_EXP - 1;
00099       ex2_u.ieee.exponent += tval >> unsafe;
00100       scale_u.d = 1.0;
00101       scale_u.ieee.exponent += tval - (tval >> unsafe);
00102 
00103       /* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial,
00104         with maximum error in [-2^-10-2^-30,2^-10+2^-30]
00105         less than 10^-19.  */
00106 
00107       x22 = (((.0096181293647031180
00108               * x + .055504110254308625)
00109              * x + .240226506959100583)
00110             * x + .69314718055994495) * ex2_u.d;
00111 
00112       /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
00113       fesetenv (&oldenv);
00114 
00115       result = x22 * x + ex2_u.d;
00116 
00117       if (!unsafe)
00118        return result;
00119       else
00120        return result * scale_u.d;
00121     }
00122   /* Exceptional cases:  */
00123   else if (isless (x, himark))
00124     {
00125       if (__isinf (x))
00126        /* e^-inf == 0, with no error.  */
00127        return 0;
00128       else
00129        /* Underflow */
00130        return TWOM1000 * TWOM1000;
00131     }
00132   else
00133     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
00134     return TWO1023*x;
00135 }