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 /* INDENT OFF */ 31 /* 32 * double __k_cexp(double x, int *n); 33 * Returns the exponential of x in the form of 2**n * y, y=__k_cexp(x,&n). 34 * 35 * Method 36 * 1. Argument reduction: 37 * Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658. 38 * Given x, find r and integer k such that 39 * 40 * x = k*ln2 + r, |r| <= 0.5*ln2. 41 * 42 * Here r will be represented as r = hi-lo for better 43 * accuracy. 44 * 45 * 2. Approximation of exp(r) by a special rational function on 46 * the interval [0,0.34658]: 47 * Write 48 * R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ... 49 * We use a special Remez algorithm on [0,0.34658] to generate 50 * a polynomial of degree 5 to approximate R. The maximum error 51 * of this polynomial approximation is bounded by 2**-59. In 52 * other words, 53 * R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5 54 * (where z=r*r, and the values of P1 to P5 are listed below) 55 * and 56 * | 5 | -59 57 * | 2.0+P1*z+...+P5*z - R(z) | <= 2 58 * | | 59 * The computation of exp(r) thus becomes 60 * 2*r 61 * exp(r) = 1 + ------- 62 * R - r 63 * r*R1(r) 64 * = 1 + r + ----------- (for better accuracy) 65 * 2 - R1(r) 66 * where 67 * 2 4 10 68 * R1(r) = r - (P1*r + P2*r + ... + P5*r ). 69 * 70 * 3. Return n = k and __k_cexp = exp(r). 71 * 72 * Special cases: 73 * exp(INF) is INF, exp(NaN) is NaN; 74 * exp(-INF) is 0, and 75 * for finite argument, only exp(0)=1 is exact. 76 * 77 * Range and Accuracy: 78 * When |x| is really big, say |x| > 50000, the accuracy 79 * is not important because the ultimate result will over or under 80 * flow. So we will simply replace n = 50000 and r = 0.0. For 81 * moderate size x, according to an error analysis, the error is 82 * always less than 1 ulp (unit in the last place). 83 * 84 * Constants: 85 * The hexadecimal values are the intended ones for the following 86 * constants. The decimal values may be used, provided that the 87 * compiler will convert from decimal to binary accurately enough 88 * to produce the hexadecimal values shown. 89 */ 90 /* INDENT ON */ 91 92 #include "libm.h" /* __k_cexp */ 93 #include "complex_wrapper.h" /* HI_WORD/LO_WORD */ 94 95 /* INDENT OFF */ 96 static const double 97 one = 1.0, 98 two128 = 3.40282366920938463463e+38, 99 halF[2] = { 100 0.5, -0.5, 101 }, 102 ln2HI[2] = { 103 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */ 104 -6.93147180369123816490e-01, /* 0xbfe62e42, 0xfee00000 */ 105 }, 106 ln2LO[2] = { 107 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */ 108 -1.90821492927058770002e-10, /* 0xbdea39ef, 0x35793c76 */ 109 }, 110 invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */ 111 P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */ 112 P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */ 113 P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */ 114 P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */ 115 P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */ 116 /* INDENT ON */ 117 118 double 119 __k_cexp(double x, int *n) { 120 double hi = 0.0L, lo = 0.0L, c, t; 121 int k, xsb; 122 unsigned hx, lx; 123 124 hx = HI_WORD(x); /* high word of x */ 125 lx = LO_WORD(x); /* low word of x */ 126 xsb = (hx >> 31) & 1; /* sign bit of x */ 127 hx &= 0x7fffffff; /* high word of |x| */ 128 129 /* filter out non-finite argument */ 130 if (hx >= 0x40e86a00) { /* if |x| > 50000 */ 131 if (hx >= 0x7ff00000) { 132 *n = 1; 133 if (((hx & 0xfffff) | lx) != 0) 134 return (x + x); /* NaN */ 135 else 136 return ((xsb == 0) ? x : 0.0); 137 /* exp(+-inf)={inf,0} */ 138 } 139 *n = (xsb == 0) ? 50000 : -50000; 140 return (one + ln2LO[1] * ln2LO[1]); /* generate inexact */ 141 } 142 143 *n = 0; 144 /* argument reduction */ 145 if (hx > 0x3fd62e42) { /* if |x| > 0.5 ln2 */ 146 if (hx < 0x3FF0A2B2) { /* and |x| < 1.5 ln2 */ 147 hi = x - ln2HI[xsb]; 148 lo = ln2LO[xsb]; 149 k = 1 - xsb - xsb; 150 } else { 151 k = (int) (invln2 * x + halF[xsb]); 152 t = k; 153 hi = x - t * ln2HI[0]; 154 /* t*ln2HI is exact for t<2**20 */ 155 lo = t * ln2LO[0]; 156 } 157 x = hi - lo; 158 *n = k; 159 } else if (hx < 0x3e300000) { /* when |x|<2**-28 */ 160 return (one + x); 161 } else 162 k = 0; 163 164 /* x is now in primary range */ 165 t = x * x; 166 c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5)))); 167 if (k == 0) 168 return (one - ((x * c) / (c - 2.0) - x)); 169 else { 170 t = one - ((lo - (x * c) / (2.0 - c)) - hi); 171 if (k > 128) { 172 t *= two128; 173 *n = k - 128; 174 } else if (k > 0) { 175 HI_WORD(t) += (k << 20); 176 *n = 0; 177 } 178 return (t); 179 } 180 } 181