Back to index

glibc  2.9
e_exp2f.c
Go to the documentation of this file.
00001 /* Single-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, and for
00026    single-precision.
00027    */
00028 #ifndef _GNU_SOURCE
00029 # define _GNU_SOURCE
00030 #endif
00031 #include <stdlib.h>
00032 #include <float.h>
00033 #include <ieee754.h>
00034 #include <math.h>
00035 #include <fenv.h>
00036 #include <inttypes.h>
00037 #include <math_private.h>
00038 
00039 #include "t_exp2f.h"
00040 
00041 static const volatile float TWOM100 = 7.88860905e-31;
00042 static const volatile float TWO127 = 1.7014118346e+38;
00043 
00044 float
00045 __ieee754_exp2f (float x)
00046 {
00047   static const float himark = (float) FLT_MAX_EXP;
00048   static const float lomark = (float) (FLT_MIN_EXP - FLT_MANT_DIG - 1);
00049 
00050   /* Check for usual case.  */
00051   if (isless (x, himark) && isgreaterequal (x, lomark))
00052     {
00053       static const float THREEp14 = 49152.0;
00054       int tval, unsafe;
00055       float rx, x22, result;
00056       union ieee754_float ex2_u, scale_u;
00057       fenv_t oldenv;
00058 
00059       feholdexcept (&oldenv);
00060 #ifdef FE_TONEAREST
00061       /* If we don't have this, it's too bad.  */
00062       fesetround (FE_TONEAREST);
00063 #endif
00064 
00065       /* 1. Argument reduction.
00066         Choose integers ex, -128 <= t < 128, and some real
00067         -1/512 <= x1 <= 1/512 so that
00068         x = ex + t/512 + x1.
00069 
00070         First, calculate rx = ex + t/256.  */
00071       rx = x + THREEp14;
00072       rx -= THREEp14;
00073       x -= rx;  /* Compute x=x1. */
00074       /* Compute tval = (ex*256 + t)+128.
00075         Now, t = (tval mod 256)-128 and ex=tval/256  [that's mod, NOT %; and
00076         /-round-to-nearest not the usual c integer /].  */
00077       tval = (int) (rx * 256.0f + 128.0f);
00078 
00079       /* 2. Adjust for accurate table entry.
00080         Find e so that
00081         x = ex + t/256 + e + x2
00082         where -7e-4 < e < 7e-4, and
00083         (float)(2^(t/256+e))
00084         is accurate to one part in 2^-64.  */
00085 
00086       /* 'tval & 255' is the same as 'tval%256' except that it's always
00087         positive.
00088         Compute x = x2.  */
00089       x -= __exp2f_deltatable[tval & 255];
00090 
00091       /* 3. Compute ex2 = 2^(t/255+e+ex).  */
00092       ex2_u.f = __exp2f_atable[tval & 255];
00093       tval >>= 8;
00094       unsafe = abs(tval) >= -FLT_MIN_EXP - 1;
00095       ex2_u.ieee.exponent += tval >> unsafe;
00096       scale_u.f = 1.0;
00097       scale_u.ieee.exponent += tval - (tval >> unsafe);
00098 
00099       /* 4. Approximate 2^x2 - 1, using a second-degree polynomial,
00100         with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14]
00101         less than 1.3e-10.  */
00102 
00103       x22 = (.24022656679f * x + .69314736128f) * ex2_u.f;
00104 
00105       /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
00106       fesetenv (&oldenv);
00107 
00108       result = x22 * x + ex2_u.f;
00109 
00110       if (!unsafe)
00111        return result;
00112       else
00113        return result * scale_u.f;
00114     }
00115   /* Exceptional cases:  */
00116   else if (isless (x, himark))
00117     {
00118       if (__isinff (x))
00119        /* e^-inf == 0, with no error.  */
00120        return 0;
00121       else
00122        /* Underflow */
00123        return TWOM100 * TWOM100;
00124     }
00125   else
00126     /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
00127     return TWO127*x;
00128 }