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 /* 27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 28 * Use is subject to license terms. 29 */ 30 31 #pragma weak __log = log 32 33 34 /* 35 * log(x) 36 * Table look-up algorithm with product polynomial approximation. 37 * By K.C. Ng, Oct 23, 2004. Updated Oct 18, 2005. 38 * 39 * (a). For x in [1-0.125, 1+0.1328125], using a special approximation: 40 * Let f = x - 1 and z = f*f. 41 * return f + ((a1*z) * 42 * ((a2 + (a3*f)*(a4+f)) + (f*z)*(a5+f))) * 43 * (((a6 + f*(a7+f)) + (f*z)*(a8+f)) * 44 * ((a9 + (a10*f)*(a11+f)) + (f*z)*(a12+f))) 45 * a1 -6.88821452420390473170286327331268694251775741577e-0002, 46 * a2 1.97493380704769294631262255279580131173133850098e+0000, 47 * a3 2.24963218866067560242072431719861924648284912109e+0000, 48 * a4 -9.02975906958474405783476868236903101205825805664e-0001, 49 * a5 -1.47391630715542865104339398385491222143173217773e+0000, 50 * a6 1.86846544648220058704168877738993614912033081055e+0000, 51 * a7 1.82277370459347465292410106485476717352867126465e+0000, 52 * a8 1.25295479915214102994980294170090928673744201660e+0000, 53 * a9 1.96709676945198275177517643896862864494323730469e+0000, 54 * a10 -4.00127989749189894030934055990655906498432159424e-0001, 55 * a11 3.01675528558798333733648178167641162872314453125e+0000, 56 * a12 -9.52325445049240770778453679668018594384193420410e-0001, 57 * 58 * with remez error |(log(1+f) - P(f))/f| <= 2**-56.81 and 59 * 60 * (b). For 0.09375 <= x < 24 61 * Use an 8-bit table look-up (3-bit for exponent and 5 bit for 62 * significand): 63 * Let ix stands for the high part of x in IEEE double format. 64 * Since 0.09375 <= x < 24, we have 65 * 0x3fb80000 <= ix < 0x40380000. 66 * Let j = (ix - 0x3fb80000) >> 15. Then 0 <= j < 256. Choose 67 * a Y[j] such that HIWORD(Y[j]) ~ 0x3fb8400 + (j<<15) (the middle 68 * number between 0x3fb80000 + (j<<15) and 3fb80000 + ((j+1)<<15)), 69 * and at the same time 1/Y[j] as well as log(Y[j]) are very close 70 * to 53-bits floating point numbers. 71 * A table of Y[j], 1/Y[j], and log(Y[j]) are pre-computed and thus 72 * log(x) = log(Y[j]) + log(1 + (x-Y[j])*(1/Y[j])) 73 * = log(Y[j]) + log(1 + s) 74 * where 75 * s = (x-Y[j])*(1/Y[j]) 76 * We compute max (x-Y[j])*(1/Y[j]) for the chosen Y[j] and obtain 77 * |s| < 0.0154. By applying remez algorithm with Product Polynomial 78 * Approximiation, we find the following approximated of log(1+s) 79 * (b1*s)*(b2+s*(b3+s))*((b4+s*b5)+(s*s)*(b6+s))*(b7+s*(b8+s)) 80 * with remez error |log(1+s) - P(s)| <= 2**-63.5 81 * 82 * (c). Otherwise, get "n", the exponent of x, and then normalize x to 83 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5 84 * significant bits. Then 85 * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]). 86 * 87 * Special cases: 88 * log(x) is NaN with signal if x < 0 (including -INF) ; 89 * log(+INF) is +INF; log(0) is -INF with signal; 90 * log(NaN) is that NaN with no signal. 91 * 92 * Maximum error observed: less than 0.90 ulp 93 * 94 * Constants: 95 * The hexadecimal values are the intended ones for the following constants. 96 * The decimal values may be used, provided that the compiler will convert 97 * from decimal to binary accurately enough to produce the hexadecimal values 98 * shown. 99 */ 100 101 #include "libm.h" 102 103 extern const double _TBL_log[]; 104 105 static const double P[] = { 106 /* ONE */ 107 1.0, 108 /* TWO52 */ 4503599627370496.0, 109 /* LN2HI */ 6.93147180369123816490e-01, /* 3fe62e42, fee00000 */ 110 /* LN2LO */ 1.90821492927058770002e-10, /* 3dea39ef, 35793c76 */ 111 /* A1 */ -6.88821452420390473170286327331268694251775741577e-0002, 112 /* A2 */ 1.97493380704769294631262255279580131173133850098e+0000, 113 /* A3 */ 2.24963218866067560242072431719861924648284912109e+0000, 114 /* A4 */ -9.02975906958474405783476868236903101205825805664e-0001, 115 /* A5 */ -1.47391630715542865104339398385491222143173217773e+0000, 116 /* A6 */ 1.86846544648220058704168877738993614912033081055e+0000, 117 /* A7 */ 1.82277370459347465292410106485476717352867126465e+0000, 118 /* A8 */ 1.25295479915214102994980294170090928673744201660e+0000, 119 /* A9 */ 1.96709676945198275177517643896862864494323730469e+0000, 120 /* A10 */ -4.00127989749189894030934055990655906498432159424e-0001, 121 /* A11 */ 3.01675528558798333733648178167641162872314453125e+0000, 122 /* A12 */ -9.52325445049240770778453679668018594384193420410e-0001, 123 /* B1 */ -1.25041641589283658575482149899471551179885864258e-0001, 124 /* B2 */ 1.87161713283355151891381127914642725337613123482e+0000, 125 /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000, 126 /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000, 127 /* B5 */ 1.64822828085258366037635369139024987816810607910e+0000, 128 /* B6 */ -1.24409107065868340669112512841820716857910156250e+0000, 129 /* B7 */ 1.70534231658220414296067701798165217041969299316e+0000, 130 /* B8 */ 1.99196833784655646937267192697618156671524047852e+0000, 131 }; 132 133 #define ONE P[0] 134 #define TWO52 P[1] 135 #define LN2HI P[2] 136 #define LN2LO P[3] 137 #define A1 P[4] 138 #define A2 P[5] 139 #define A3 P[6] 140 #define A4 P[7] 141 #define A5 P[8] 142 #define A6 P[9] 143 #define A7 P[10] 144 #define A8 P[11] 145 #define A9 P[12] 146 #define A10 P[13] 147 #define A11 P[14] 148 #define A12 P[15] 149 #define B1 P[16] 150 #define B2 P[17] 151 #define B3 P[18] 152 #define B4 P[19] 153 #define B5 P[20] 154 #define B6 P[21] 155 #define B7 P[22] 156 #define B8 P[23] 157 158 double 159 log(double x) 160 { 161 double *tb, dn, dn1, s, z, r, w; 162 int i, hx, ix, n, lx; 163 164 n = 0; 165 hx = ((int *)&x)[HIWORD]; 166 ix = hx & 0x7fffffff; 167 lx = ((int *)&x)[LOWORD]; 168 169 /* subnormal,0,negative,inf,nan */ 170 if ((hx + 0x100000) < 0x200000) { 171 if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */ 172 return (x * x); 173 174 if (((hx << 1) | lx) == 0) /* zero */ 175 return (_SVID_libm_err(x, x, 16)); 176 177 if (hx < 0) /* negative */ 178 return (_SVID_libm_err(x, x, 17)); 179 180 if (((hx - 0x7ff00000) | lx) == 0) /* +inf */ 181 return (x); 182 183 /* x must be positive and subnormal */ 184 x *= TWO52; 185 n = -52; 186 ix = ((int *)&x)[HIWORD]; 187 lx = ((int *)&x)[LOWORD]; 188 } 189 190 i = ix >> 19; 191 192 if (i >= 0x7f7 && i <= 0x806) { 193 /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */ 194 if (ix >= 0x3fec0000 && ix < 0x3ff22000) { 195 /* 0.875 <= x < 1.125 */ 196 s = x - ONE; 197 z = s * s; 198 199 if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */ 200 return (z); 201 202 r = (A10 * s) * (A11 + s); 203 w = z * s; 204 return (s + ((A1 * z) * (A2 + ((A3 * s) * (A4 + s) + w * 205 (A5 + s)))) * ((A6 + (s * (A7 + s) + w * 206 (A8 + s))) * (A9 + (r + w * (A12 + s))))); 207 } else { 208 i = (ix - 0x3fb80000) >> 15; 209 tb = (double *)_TBL_log + (i + i + i); 210 s = (x - tb[0]) * tb[1]; 211 return (tb[2] + ((B1 * s) * (B2 + s * (B3 + s))) * 212 (((B4 + s * B5) + (s * s) * (B6 + s)) * 213 (B7 + s * (B8 + s)))); 214 } 215 } else { 216 dn = (double)(n + ((ix >> 20) - 0x3ff)); 217 dn1 = dn * LN2HI; 218 i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */ 219 ((int *)&x)[HIWORD] = i; 220 i = (i - 0x3fb80000) >> 15; 221 tb = (double *)_TBL_log + (i + i + i); 222 s = (x - tb[0]) * tb[1]; 223 dn = dn * LN2LO + tb[2]; 224 return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) * 225 (((B4 + s * B5) + (s * s) * (B6 + s)) * 226 (B7 + s * (B8 + s))))); 227 } 228 }