1 /* 2 * CDDL HEADER START 3 * 4 * The contents of this file are subject to the terms of the 5 * Common Development and Distribution License (the "License"). 6 * You may not use this file except in compliance with the License. 7 * 8 * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9 * or http://www.opensolaris.org/os/licensing. 10 * See the License for the specific language governing permissions 11 * and limitations under the License. 12 * 13 * When distributing Covered Code, include this CDDL HEADER in each 14 * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15 * If applicable, add the following below this CDDL HEADER, with the 16 * fields enclosed by brackets "[]" replaced with your own identifying 17 * information: Portions Copyright [yyyy] [name of copyright owner] 18 * 19 * CDDL HEADER END 20 */ 21 22 /* 23 * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24 */ 25 /* 26 * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27 * Use is subject to license terms. 28 */ 29 30 #pragma weak expm1 = __expm1 31 32 /* INDENT OFF */ 33 /* 34 * expm1(x) 35 * Returns exp(x)-1, the exponential of x minus 1. 36 * 37 * Method 38 * 1. Arugment reduction: 39 * Given x, find r and integer k such that 40 * 41 * x = k*ln2 + r, |r| <= 0.5*ln2 ~ 0.34658 42 * 43 * Here a correction term c will be computed to compensate 44 * the error in r when rounded to a floating-point number. 45 * 46 * 2. Approximating expm1(r) by a special rational function on 47 * the interval [0,0.34658]: 48 * Since 49 * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 - r^4/360 + ... 50 * we define R1(r*r) by 51 * r*(exp(r)+1)/(exp(r)-1) = 2+ r^2/6 * R1(r*r) 52 * That is, 53 * R1(r**2) = 6/r *((exp(r)+1)/(exp(r)-1) - 2/r) 54 * = 6/r * ( 1 + 2.0*(1/(exp(r)-1) - 1/r)) 55 * = 1 - r^2/60 + r^4/2520 - r^6/100800 + ... 56 * We use a special Reme algorithm on [0,0.347] to generate 57 * a polynomial of degree 5 in r*r to approximate R1. The 58 * maximum error of this polynomial approximation is bounded 59 * by 2**-61. In other words, 60 * R1(z) ~ 1.0 + Q1*z + Q2*z**2 + Q3*z**3 + Q4*z**4 + Q5*z**5 61 * where Q1 = -1.6666666666666567384E-2, 62 * Q2 = 3.9682539681370365873E-4, 63 * Q3 = -9.9206344733435987357E-6, 64 * Q4 = 2.5051361420808517002E-7, 65 * Q5 = -6.2843505682382617102E-9; 66 * (where z=r*r, and the values of Q1 to Q5 are listed below) 67 * with error bounded by 68 * | 5 | -61 69 * | 1.0+Q1*z+...+Q5*z - R1(z) | <= 2 70 * | | 71 * 72 * expm1(r) = exp(r)-1 is then computed by the following 73 * specific way which minimize the accumulation rounding error: 74 * 2 3 75 * r r [ 3 - (R1 + R1*r/2) ] 76 * expm1(r) = r + --- + --- * [--------------------] 77 * 2 2 [ 6 - r*(3 - R1*r/2) ] 78 * 79 * To compensate the error in the argument reduction, we use 80 * expm1(r+c) = expm1(r) + c + expm1(r)*c 81 * ~ expm1(r) + c + r*c 82 * Thus c+r*c will be added in as the correction terms for 83 * expm1(r+c). Now rearrange the term to avoid optimization 84 * screw up: 85 * ( 2 2 ) 86 * ({ ( r [ R1 - (3 - R1*r/2) ] ) } r ) 87 * expm1(r+c)~r - ({r*(--- * [--------------------]-c)-c} - --- ) 88 * ({ ( 2 [ 6 - r*(3 - R1*r/2) ] ) } 2 ) 89 * ( ) 90 * 91 * = r - E 92 * 3. Scale back to obtain expm1(x): 93 * From step 1, we have 94 * expm1(x) = either 2^k*[expm1(r)+1] - 1 95 * = or 2^k*[expm1(r) + (1-2^-k)] 96 * 4. Implementation notes: 97 * (A). To save one multiplication, we scale the coefficient Qi 98 * to Qi*2^i, and replace z by (x^2)/2. 99 * (B). To achieve maximum accuracy, we compute expm1(x) by 100 * (i) if x < -56*ln2, return -1.0, (raise inexact if x != inf) 101 * (ii) if k=0, return r-E 102 * (iii) if k=-1, return 0.5*(r-E)-0.5 103 * (iv) if k=1 if r < -0.25, return 2*((r+0.5)- E) 104 * else return 1.0+2.0*(r-E); 105 * (v) if (k<-2||k>56) return 2^k(1-(E-r)) - 1 (or exp(x)-1) 106 * (vi) if k <= 20, return 2^k((1-2^-k)-(E-r)), else 107 * (vii) return 2^k(1-((E+2^-k)-r)) 108 * 109 * Special cases: 110 * expm1(INF) is INF, expm1(NaN) is NaN; 111 * expm1(-INF) is -1, and 112 * for finite argument, only expm1(0)=0 is exact. 113 * 114 * Accuracy: 115 * according to an error analysis, the error is always less than 116 * 1 ulp (unit in the last place). 117 * 118 * Misc. info. 119 * For IEEE double 120 * if x > 7.09782712893383973096e+02 then expm1(x) overflow 121 * 122 * Constants: 123 * The hexadecimal values are the intended ones for the following 124 * constants. The decimal values may be used, provided that the 125 * compiler will convert from decimal to binary accurately enough 126 * to produce the hexadecimal values shown. 127 */ 128 /* INDENT ON */ 129 130 #include "libm_synonyms.h" /* __expm1 */ 131 #include "libm_macros.h" 132 #include <math.h> 133 134 static const double xxx[] = { 135 /* one */ 1.0, 136 /* huge */ 1.0e+300, 137 /* tiny */ 1.0e-300, 138 /* o_threshold */ 7.09782712893383973096e+02, /* 40862E42 FEFA39EF */ 139 /* ln2_hi */ 6.93147180369123816490e-01, /* 3FE62E42 FEE00000 */ 140 /* ln2_lo */ 1.90821492927058770002e-10, /* 3DEA39EF 35793C76 */ 141 /* invln2 */ 1.44269504088896338700e+00, /* 3FF71547 652B82FE */ 142 /* scaled coefficients related to expm1 */ 143 /* Q1 */ -3.33333333333331316428e-02, /* BFA11111 111110F4 */ 144 /* Q2 */ 1.58730158725481460165e-03, /* 3F5A01A0 19FE5585 */ 145 /* Q3 */ -7.93650757867487942473e-05, /* BF14CE19 9EAADBB7 */ 146 /* Q4 */ 4.00821782732936239552e-06, /* 3ED0CFCA 86E65239 */ 147 /* Q5 */ -2.01099218183624371326e-07 /* BE8AFDB7 6E09C32D */ 148 }; 149 #define one xxx[0] 150 #define huge xxx[1] 151 #define tiny xxx[2] 152 #define o_threshold xxx[3] 153 #define ln2_hi xxx[4] 154 #define ln2_lo xxx[5] 155 #define invln2 xxx[6] 156 #define Q1 xxx[7] 157 #define Q2 xxx[8] 158 #define Q3 xxx[9] 159 #define Q4 xxx[10] 160 #define Q5 xxx[11] 161 162 double 163 expm1(double x) { 164 double y, hi, lo, c = 0.0L, t, e, hxs, hfx, r1; 165 int k, xsb; 166 unsigned hx; 167 168 hx = ((unsigned *) &x)[HIWORD]; /* high word of x */ 169 xsb = hx & 0x80000000; /* sign bit of x */ 170 if (xsb == 0) 171 y = x; 172 else 173 y = -x; /* y = |x| */ 174 hx &= 0x7fffffff; /* high word of |x| */ 175 176 /* filter out huge and non-finite argument */ 177 /* for example exp(38)-1 is approximately 3.1855932e+16 */ 178 if (hx >= 0x4043687A) { 179 /* if |x|>=56*ln2 (~38.8162...) */ 180 if (hx >= 0x40862E42) { /* if |x|>=709.78... -> inf */ 181 if (hx >= 0x7ff00000) { 182 if (((hx & 0xfffff) | ((int *) &x)[LOWORD]) 183 != 0) 184 return (x * x); /* + -> * for Cheetah */ 185 else 186 /* exp(+-inf)={inf,-1} */ 187 return (xsb == 0 ? x : -1.0); 188 } 189 if (x > o_threshold) 190 return (huge * huge); /* overflow */ 191 } 192 if (xsb != 0) { /* x < -56*ln2, return -1.0 w/inexact */ 193 if (x + tiny < 0.0) /* raise inexact */ 194 return (tiny - one); /* return -1 */ 195 } 196 } 197 198 /* argument reduction */ 199 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 200 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 201 if (xsb == 0) { /* positive number */ 202 hi = x - ln2_hi; 203 lo = ln2_lo; 204 k = 1; 205 } else { 206 /* negative number */ 207 hi = x + ln2_hi; 208 lo = -ln2_lo; 209 k = -1; 210 } 211 } else { 212 /* |x| > 1.5 ln2 */ 213 k = (int) (invln2 * x + (xsb == 0 ? 0.5 : -0.5)); 214 t = k; 215 hi = x - t * ln2_hi; /* t*ln2_hi is exact here */ 216 lo = t * ln2_lo; 217 } 218 x = hi - lo; 219 c = (hi - x) - lo; /* still at |x| > 0.5 ln2 */ 220 } else if (hx < 0x3c900000) { 221 /* when |x|<2**-54, return x */ 222 t = huge + x; /* return x w/inexact when x != 0 */ 223 return (x - (t - (huge + x))); 224 } else 225 /* |x| <= 0.5 ln2 */ 226 k = 0; 227 228 /* x is now in primary range */ 229 hfx = 0.5 * x; 230 hxs = x * hfx; 231 r1 = one + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5)))); 232 t = 3.0 - r1 * hfx; 233 e = hxs * ((r1 - t) / (6.0 - x * t)); 234 if (k == 0) /* |x| <= 0.5 ln2 */ 235 return (x - (x * e - hxs)); 236 else { /* |x| > 0.5 ln2 */ 237 e = (x * (e - c) - c); 238 e -= hxs; 239 if (k == -1) 240 return (0.5 * (x - e) - 0.5); 241 if (k == 1) { 242 if (x < -0.25) 243 return (-2.0 * (e - (x + 0.5))); 244 else 245 return (one + 2.0 * (x - e)); 246 } 247 if (k <= -2 || k > 56) { /* suffice to return exp(x)-1 */ 248 y = one - (e - x); 249 ((int *) &y)[HIWORD] += k << 20; 250 return (y - one); 251 } 252 t = one; 253 if (k < 20) { 254 ((int *) &t)[HIWORD] = 0x3ff00000 - (0x200000 >> k); 255 /* t = 1 - 2^-k */ 256 y = t - (e - x); 257 ((int *) &y)[HIWORD] += k << 20; 258 } else { 259 ((int *) &t)[HIWORD] = (0x3ff - k) << 20; /* 2^-k */ 260 y = x - (e + t); 261 y += one; 262 ((int *) &y)[HIWORD] += k << 20; 263 } 264 } 265 return (y); 266 } 267