1/* Copyright JS Foundation and other contributors, http://js.foundation 2 * 3 * Licensed under the Apache License, Version 2.0 (the "License"); 4 * you may not use this file except in compliance with the License. 5 * You may obtain a copy of the License at 6 * 7 * http://www.apache.org/licenses/LICENSE-2.0 8 * 9 * Unless required by applicable law or agreed to in writing, software 10 * distributed under the License is distributed on an "AS IS" BASIS 11 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 12 * See the License for the specific language governing permissions and 13 * limitations under the License. 14 * 15 * This file is based on work under the following copyright and permission 16 * notice: 17 * 18 * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved. 19 * 20 * Permission to use, copy, modify, and distribute this 21 * software is freely granted, provided that this notice 22 * is preserved. 23 * 24 * @(#)e_exp.c 1.6 04/04/22 25 */ 26 27#include "jerry-libm-internal.h" 28 29/* exp(x) 30 * Returns the exponential of x. 31 * 32 * Method: 33 * 1. Argument reduction: 34 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 35 * Given x, find r and integer k such that 36 * 37 * x = k*ln2 + r, |r| <= 0.5*ln2. 38 * 39 * Here r will be represented as r = hi-lo for better 40 * accuracy. 41 * 42 * 2. Approximation of exp(r) by a special rational function on 43 * the interval [0,0.34658]: 44 * Write 45 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 46 * We use a special Remes algorithm on [0,0.34658] to generate 47 * a polynomial of degree 5 to approximate R. The maximum error 48 * of this polynomial approximation is bounded by 2**-59. In 49 * other words, 50 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 51 * (where z=r*r, and the values of P1 to P5 are listed below) 52 * and 53 * | 5 | -59 54 * | 2.0+P1*z+...+P5*z - R(z) | <= 2 55 * | | 56 * The computation of exp(r) thus becomes 57 * 2*r 58 * exp(r) = 1 + ------- 59 * R - r 60 * r*R1(r) 61 * = 1 + r + ----------- (for better accuracy) 62 * 2 - R1(r) 63 * where 64 * 2 4 10 65 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). 66 * 67 * 3. Scale back to obtain exp(x): 68 * From step 1, we have 69 * exp(x) = 2^k * exp(r) 70 * 71 * Special cases: 72 * exp(INF) is INF, exp(NaN) is NaN; 73 * exp(-INF) is 0, and 74 * for finite argument, only exp(0)=1 is exact. 75 * 76 * Accuracy: 77 * according to an error analysis, the error is always less than 78 * 1 ulp (unit in the last place). 79 * 80 * Misc. info: 81 * For IEEE double 82 * if x > 7.09782712893383973096e+02 then exp(x) overflow 83 * if x < -7.45133219101941108420e+02 then exp(x) underflow 84 * 85 * Constants: 86 * The hexadecimal values are the intended ones for the following 87 * constants. The decimal values may be used, provided that the 88 * compiler will convert from decimal to binary accurately enough 89 * to produce the hexadecimal values shown. 90 */ 91 92static const double halF[2] = 93{ 94 0.5, 95 -0.5, 96}; 97static const double ln2HI[2] = 98{ 99 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 100 -6.93147180369123816490e-01, /* 0xbfe62e42, 0xfee00000 */ 101}; 102static const double ln2LO[2] = 103{ 104 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 105 -1.90821492927058770002e-10, /* 0xbdea39ef, 0x35793c76 */ 106}; 107 108#define one 1.0 109#define huge 1.0e+300 110#define twom1000 9.33263618503218878990e-302 /* 2**-1000=0x01700000,0 */ 111#define o_threshold 7.09782712893383973096e+02 /* 0x40862E42, 0xFEFA39EF */ 112#define u_threshold -7.45133219101941108420e+02 /* 0xc0874910, 0xD52D3051 */ 113#define invln2 1.44269504088896338700e+00 /* 0x3ff71547, 0x652b82fe */ 114#define P1 1.66666666666666019037e-01 /* 0x3FC55555, 0x5555553E */ 115#define P2 -2.77777777770155933842e-03 /* 0xBF66C16C, 0x16BEBD93 */ 116#define P3 6.61375632143793436117e-05 /* 0x3F11566A, 0xAF25DE2C */ 117#define P4 -1.65339022054652515390e-06 /* 0xBEBBBD41, 0xC5D26BF1 */ 118#define P5 4.13813679705723846039e-08 /* 0x3E663769, 0x72BEA4D0 */ 119 120double 121exp (double x) /* default IEEE double exp */ 122{ 123 double hi, lo, c, t; 124 int k = 0, xsb; 125 unsigned hx; 126 127 hx = __HI (x); /* high word of x */ 128 xsb = (hx >> 31) & 1; /* sign bit of x */ 129 hx &= 0x7fffffff; /* high word of |x| */ 130 131 /* filter out non-finite argument */ 132 if (hx >= 0x40862E42) /* if |x| >= 709.78... */ 133 { 134 if (hx >= 0x7ff00000) 135 { 136 if (((hx & 0xfffff) | __LO (x)) != 0) /* NaN */ 137 { 138 return x + x; 139 } 140 else /* exp(+-inf) = {inf,0} */ 141 { 142 return (xsb == 0) ? x : 0.0; 143 } 144 } 145 if (x > o_threshold) /* overflow */ 146 { 147 return huge * huge; 148 } 149 if (x < u_threshold) /* underflow */ 150 { 151 return twom1000 * twom1000; 152 } 153 } 154 155 /* argument reduction */ 156 if (hx > 0x3fd62e42) /* if |x| > 0.5 ln2 */ 157 { 158 if (hx < 0x3FF0A2B2) /* and |x| < 1.5 ln2 */ 159 { 160 hi = x - ln2HI[xsb]; 161 lo = ln2LO[xsb]; 162 k = 1 - xsb - xsb; 163 } 164 else 165 { 166 k = (int) (invln2 * x + halF[xsb]); 167 t = k; 168 hi = x - t * ln2HI[0]; /* t * ln2HI is exact here */ 169 lo = t * ln2LO[0]; 170 } 171 x = hi - lo; 172 } 173 else if (hx < 0x3e300000) /* when |x| < 2**-28 */ 174 { 175 if (huge + x > one) /* trigger inexact */ 176 { 177 return one + x; 178 } 179 } 180 else 181 { 182 k = 0; 183 } 184 185 double_accessor ret; 186 187 /* x is now in primary range */ 188 t = x * x; 189 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); 190 if (k == 0) 191 { 192 return one - ((x * c) / (c - 2.0) - x); 193 } 194 else 195 { 196 ret.dbl = one - ((lo - (x * c) / (2.0 - c)) - hi); 197 } 198 if (k >= -1021) 199 { 200 ret.as_int.hi += (((unsigned int) k) << 20); /* add k to y's exponent */ 201 return ret.dbl; 202 } 203 else 204 { 205 ret.as_int.hi += ((k + 1000) << 20); /* add k to y's exponent */ 206 return ret.dbl * twom1000; 207 } 208} /* exp */ 209 210#undef one 211#undef huge 212#undef twom1000 213#undef o_threshold 214#undef u_threshold 215#undef invln2 216#undef P1 217#undef P2 218#undef P3 219#undef P4 220#undef P5 221