1/* Double-precision floating point 2^x. 2 Copyright (C) 1997, 1998, 2000, 2001 Free Software Foundation, Inc. 3 This file is part of the GNU C Library. 4 Contributed by Geoffrey Keating <geoffk@ozemail.com.au> 5 6 The GNU C Library is free software; you can redistribute it and/or 7 modify it under the terms of the GNU Lesser General Public 8 License as published by the Free Software Foundation; either 9 version 2.1 of the License, or (at your option) any later version. 10 11 The GNU C Library is distributed in the hope that it will be useful, 12 but WITHOUT ANY WARRANTY; without even the implied warranty of 13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14 Lesser General Public License for more details. 15 16 You should have received a copy of the GNU Lesser General Public 17 License along with the GNU C Library; if not, write to the Free 18 Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 19 02111-1307 USA. */ 20 21/* The basic design here is from 22 Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical 23 Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft., 24 17 (1), March 1991, pp. 26-45. 25 It has been slightly modified to compute 2^x instead of e^x. 26 */ 27#ifndef _GNU_SOURCE 28#define _GNU_SOURCE 29#endif 30#include <stdlib.h> 31#include <float.h> 32#include <ieee754.h> 33#include <math.h> 34#include <fenv.h> 35#include <inttypes.h> 36#include <math_private.h> 37 38#include "t_exp2.h" 39 40static const volatile double TWO1023 = 8.988465674311579539e+307; 41static const volatile double TWOM1000 = 9.3326361850321887899e-302; 42 43double 44__ieee754_exp2 (double x) 45{ 46 static const double himark = (double) DBL_MAX_EXP; 47 static const double lomark = (double) (DBL_MIN_EXP - DBL_MANT_DIG - 1); 48 49 /* Check for usual case. */ 50 if (isless (x, himark) && isgreaterequal (x, lomark)) 51 { 52 static const double THREEp42 = 13194139533312.0; 53 int tval, unsafe; 54 double rx, x22, result; 55 union ieee754_double ex2_u, scale_u; 56 fenv_t oldenv; 57 58 feholdexcept (&oldenv); 59#ifdef FE_TONEAREST 60 /* If we don't have this, it's too bad. */ 61 fesetround (FE_TONEAREST); 62#endif 63 64 /* 1. Argument reduction. 65 Choose integers ex, -256 <= t < 256, and some real 66 -1/1024 <= x1 <= 1024 so that 67 x = ex + t/512 + x1. 68 69 First, calculate rx = ex + t/512. */ 70 rx = x + THREEp42; 71 rx -= THREEp42; 72 x -= rx; /* Compute x=x1. */ 73 /* Compute tval = (ex*512 + t)+256. 74 Now, t = (tval mod 512)-256 and ex=tval/512 [that's mod, NOT %; and 75 /-round-to-nearest not the usual c integer /]. */ 76 tval = (int) (rx * 512.0 + 256.0); 77 78 /* 2. Adjust for accurate table entry. 79 Find e so that 80 x = ex + t/512 + e + x2 81 where -1e6 < e < 1e6, and 82 (double)(2^(t/512+e)) 83 is accurate to one part in 2^-64. */ 84 85 /* 'tval & 511' is the same as 'tval%512' except that it's always 86 positive. 87 Compute x = x2. */ 88 x -= exp2_deltatable[tval & 511]; 89 90 /* 3. Compute ex2 = 2^(t/512+e+ex). */ 91 ex2_u.d = exp2_accuratetable[tval & 511]; 92 tval >>= 9; 93 unsafe = abs(tval) >= -DBL_MIN_EXP - 1; 94 ex2_u.ieee.exponent += tval >> unsafe; 95 scale_u.d = 1.0; 96 scale_u.ieee.exponent += tval - (tval >> unsafe); 97 98 /* 4. Approximate 2^x2 - 1, using a fourth-degree polynomial, 99 with maximum error in [-2^-10-2^-30,2^-10+2^-30] 100 less than 10^-19. */ 101 102 x22 = (((.0096181293647031180 103 * x + .055504110254308625) 104 * x + .240226506959100583) 105 * x + .69314718055994495) * ex2_u.d; 106 107 /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex). */ 108 fesetenv (&oldenv); 109 110 result = x22 * x + ex2_u.d; 111 112 if (!unsafe) 113 return result; 114 else 115 return result * scale_u.d; 116 } 117 /* Exceptional cases: */ 118 else if (isless (x, himark)) 119 { 120 if (__isinf (x)) 121 /* e^-inf == 0, with no error. */ 122 return 0; 123 else 124 /* Underflow */ 125 return TWOM1000 * TWOM1000; 126 } 127 else 128 /* Return x, if x is a NaN or Inf; or overflow, otherwise. */ 129 return TWO1023*x; 130} 131