1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis /* 22*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 23*25c28e83SPiotr Jasiukajtis */ 24*25c28e83SPiotr Jasiukajtis /* 25*25c28e83SPiotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc. All rights reserved. 26*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 27*25c28e83SPiotr Jasiukajtis */ 28*25c28e83SPiotr Jasiukajtis 29*25c28e83SPiotr Jasiukajtis #pragma weak log10 = __log10 30*25c28e83SPiotr Jasiukajtis 31*25c28e83SPiotr Jasiukajtis /* INDENT OFF */ 32*25c28e83SPiotr Jasiukajtis /* 33*25c28e83SPiotr Jasiukajtis * log10(x) = log(x)/log10 34*25c28e83SPiotr Jasiukajtis * 35*25c28e83SPiotr Jasiukajtis * Base on Table look-up algorithm with product polynomial 36*25c28e83SPiotr Jasiukajtis * approximation for log(x). 37*25c28e83SPiotr Jasiukajtis * 38*25c28e83SPiotr Jasiukajtis * By K.C. Ng, Nov 29, 2004 39*25c28e83SPiotr Jasiukajtis * 40*25c28e83SPiotr Jasiukajtis * (a). For x in [1-0.125, 1+0.125], from log.c we have 41*25c28e83SPiotr Jasiukajtis * log(x) = f + ((a1*f^2) * 42*25c28e83SPiotr Jasiukajtis * ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) * 43*25c28e83SPiotr Jasiukajtis * (((a6 + f*(a7+f)) + (f^3)*(a8+f)) * 44*25c28e83SPiotr Jasiukajtis * ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f))) 45*25c28e83SPiotr Jasiukajtis * where f = x - 1. 46*25c28e83SPiotr Jasiukajtis * (i) modify a1 <- a1 / log10 47*25c28e83SPiotr Jasiukajtis * (ii) 1/log10 = 0.4342944819... 48*25c28e83SPiotr Jasiukajtis * = 0.4375 - 0.003205518... (7 bit shift) 49*25c28e83SPiotr Jasiukajtis * Let lgv = 0.4375 - 1/log10, then 50*25c28e83SPiotr Jasiukajtis * lgv = 0.003205518096748172348871081083395..., 51*25c28e83SPiotr Jasiukajtis * (iii) f*0.4375 is exact because f has 3 trailing zero. 52*25c28e83SPiotr Jasiukajtis * (iv) Thus, log10(x) = f*0.4375 - (lgv*f - PPoly) 53*25c28e83SPiotr Jasiukajtis * 54*25c28e83SPiotr Jasiukajtis * (b). For 0.09375 <= x < 24 55*25c28e83SPiotr Jasiukajtis * Let j = (ix - 0x3fb80000) >> 15. Look up Y[j], 1/Y[j], and log(Y[j]) 56*25c28e83SPiotr Jasiukajtis * from _TBL_log.c. Then 57*25c28e83SPiotr Jasiukajtis * log10(x) = log10(Y[j]) + log10(1 + (x-Y[j])*(1/Y[j])) 58*25c28e83SPiotr Jasiukajtis * = log(Y[j])(1/log10) + log10(1 + s) 59*25c28e83SPiotr Jasiukajtis * where 60*25c28e83SPiotr Jasiukajtis * s = (x-Y[j])*(1/Y[j]) 61*25c28e83SPiotr Jasiukajtis * From log.c, we have log(1+s) = 62*25c28e83SPiotr Jasiukajtis * 2 2 2 63*25c28e83SPiotr Jasiukajtis * (b s) (b + b s + s ) [b + b s + s (b + s)] (b + b s + s ) 64*25c28e83SPiotr Jasiukajtis * 1 2 3 4 5 6 7 8 65*25c28e83SPiotr Jasiukajtis * 66*25c28e83SPiotr Jasiukajtis * By setting b1 <- b1/log10, we have 67*25c28e83SPiotr Jasiukajtis * log10(x) = 0.4375 * T - (lgv * T - POLY(s)) 68*25c28e83SPiotr Jasiukajtis * 69*25c28e83SPiotr Jasiukajtis * (c). Otherwise, get "n", the exponent of x, and then normalize x to 70*25c28e83SPiotr Jasiukajtis * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5 71*25c28e83SPiotr Jasiukajtis * significant bits. Then 72*25c28e83SPiotr Jasiukajtis * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]). 73*25c28e83SPiotr Jasiukajtis * log10(x) = n*(ln2/ln10) + log10(z). 74*25c28e83SPiotr Jasiukajtis * 75*25c28e83SPiotr Jasiukajtis * Special cases: 76*25c28e83SPiotr Jasiukajtis * log10(x) is NaN with signal if x < 0 (including -INF) ; 77*25c28e83SPiotr Jasiukajtis * log10(+INF) is +INF; log10(0) is -INF with signal; 78*25c28e83SPiotr Jasiukajtis * log10(NaN) is that NaN with no signal. 79*25c28e83SPiotr Jasiukajtis * 80*25c28e83SPiotr Jasiukajtis * Maximum error observed: less than 0.89 ulp 81*25c28e83SPiotr Jasiukajtis * 82*25c28e83SPiotr Jasiukajtis * Constants: 83*25c28e83SPiotr Jasiukajtis * The hexadecimal values are the intended ones for the following constants. 84*25c28e83SPiotr Jasiukajtis * The decimal values may be used, provided that the compiler will convert 85*25c28e83SPiotr Jasiukajtis * from decimal to binary accurately enough to produce the hexadecimal values 86*25c28e83SPiotr Jasiukajtis * shown. 87*25c28e83SPiotr Jasiukajtis */ 88*25c28e83SPiotr Jasiukajtis /* INDENT ON */ 89*25c28e83SPiotr Jasiukajtis 90*25c28e83SPiotr Jasiukajtis #include "libm.h" 91*25c28e83SPiotr Jasiukajtis 92*25c28e83SPiotr Jasiukajtis extern const double _TBL_log[]; 93*25c28e83SPiotr Jasiukajtis 94*25c28e83SPiotr Jasiukajtis static const double P[] = { 95*25c28e83SPiotr Jasiukajtis /* ONE */ 1.0, 96*25c28e83SPiotr Jasiukajtis /* TWO52 */ 4503599627370496.0, 97*25c28e83SPiotr Jasiukajtis /* LNAHI */ 3.01029995607677847147e-01, /* 3FD34413 50900000 */ 98*25c28e83SPiotr Jasiukajtis /* LNALO */ 5.63033480667509769841e-11, /* 3DCEF3FD E623E256 */ 99*25c28e83SPiotr Jasiukajtis /* A1 */ -2.9142521960136582507385480707044582802184e-02, 100*25c28e83SPiotr Jasiukajtis /* A2 */ 1.99628461483039965074226529395673424005508422852e+0000, 101*25c28e83SPiotr Jasiukajtis /* A3 */ 2.26812367662950720159642514772713184356689453125e+0000, 102*25c28e83SPiotr Jasiukajtis /* A4 */ -9.05030639084976384900471657601883634924888610840e-0001, 103*25c28e83SPiotr Jasiukajtis /* A5 */ -1.48275767132434044270894446526654064655303955078e+0000, 104*25c28e83SPiotr Jasiukajtis /* A6 */ 1.88158320939722756293122074566781520843505859375e+0000, 105*25c28e83SPiotr Jasiukajtis /* A7 */ 1.83309386046986411145098827546462416648864746094e+0000, 106*25c28e83SPiotr Jasiukajtis /* A8 */ 1.24847063988317086291601754055591300129890441895e+0000, 107*25c28e83SPiotr Jasiukajtis /* A9 */ 1.98372421445537705508854742220137268304824829102e+0000, 108*25c28e83SPiotr Jasiukajtis /* A10 */ -3.94711735767898475035764249696512706577777862549e-0001, 109*25c28e83SPiotr Jasiukajtis /* A11 */ 3.07890395362954372160402272129431366920471191406e+0000, 110*25c28e83SPiotr Jasiukajtis /* A12 */ -9.60099585275022149311041630426188930869102478027e-0001, 111*25c28e83SPiotr Jasiukajtis /* B1 */ -5.4304894950350052960838096752491540286689e-02, 112*25c28e83SPiotr Jasiukajtis /* B2 */ 1.87161713283355151891381127914642725337613123482e+0000, 113*25c28e83SPiotr Jasiukajtis /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000, 114*25c28e83SPiotr Jasiukajtis /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000, 115*25c28e83SPiotr Jasiukajtis /* B5 */ 1.64822828085258366037635369139024987816810607910e+0000, 116*25c28e83SPiotr Jasiukajtis /* B6 */ -1.24409107065868340669112512841820716857910156250e+0000, 117*25c28e83SPiotr Jasiukajtis /* B7 */ 1.70534231658220414296067701798165217041969299316e+0000, 118*25c28e83SPiotr Jasiukajtis /* B8 */ 1.99196833784655646937267192697618156671524047852e+0000, 119*25c28e83SPiotr Jasiukajtis /* LGH */ 0.4375, 120*25c28e83SPiotr Jasiukajtis /* LGL */ 0.003205518096748172348871081083395, 121*25c28e83SPiotr Jasiukajtis /* LG10V */ 0.43429448190325182765112891891660509576226, 122*25c28e83SPiotr Jasiukajtis }; 123*25c28e83SPiotr Jasiukajtis 124*25c28e83SPiotr Jasiukajtis #define ONE P[0] 125*25c28e83SPiotr Jasiukajtis #define TWO52 P[1] 126*25c28e83SPiotr Jasiukajtis #define LNAHI P[2] 127*25c28e83SPiotr Jasiukajtis #define LNALO P[3] 128*25c28e83SPiotr Jasiukajtis #define A1 P[4] 129*25c28e83SPiotr Jasiukajtis #define A2 P[5] 130*25c28e83SPiotr Jasiukajtis #define A3 P[6] 131*25c28e83SPiotr Jasiukajtis #define A4 P[7] 132*25c28e83SPiotr Jasiukajtis #define A5 P[8] 133*25c28e83SPiotr Jasiukajtis #define A6 P[9] 134*25c28e83SPiotr Jasiukajtis #define A7 P[10] 135*25c28e83SPiotr Jasiukajtis #define A8 P[11] 136*25c28e83SPiotr Jasiukajtis #define A9 P[12] 137*25c28e83SPiotr Jasiukajtis #define A10 P[13] 138*25c28e83SPiotr Jasiukajtis #define A11 P[14] 139*25c28e83SPiotr Jasiukajtis #define A12 P[15] 140*25c28e83SPiotr Jasiukajtis #define B1 P[16] 141*25c28e83SPiotr Jasiukajtis #define B2 P[17] 142*25c28e83SPiotr Jasiukajtis #define B3 P[18] 143*25c28e83SPiotr Jasiukajtis #define B4 P[19] 144*25c28e83SPiotr Jasiukajtis #define B5 P[20] 145*25c28e83SPiotr Jasiukajtis #define B6 P[21] 146*25c28e83SPiotr Jasiukajtis #define B7 P[22] 147*25c28e83SPiotr Jasiukajtis #define B8 P[23] 148*25c28e83SPiotr Jasiukajtis #define LGH P[24] 149*25c28e83SPiotr Jasiukajtis #define LGL P[25] 150*25c28e83SPiotr Jasiukajtis #define LG10V P[26] 151*25c28e83SPiotr Jasiukajtis 152*25c28e83SPiotr Jasiukajtis double 153*25c28e83SPiotr Jasiukajtis log10(double x) { 154*25c28e83SPiotr Jasiukajtis double *tb, dn, dn1, s, z, r, w; 155*25c28e83SPiotr Jasiukajtis int i, hx, ix, n, lx; 156*25c28e83SPiotr Jasiukajtis 157*25c28e83SPiotr Jasiukajtis n = 0; 158*25c28e83SPiotr Jasiukajtis hx = ((int *)&x)[HIWORD]; 159*25c28e83SPiotr Jasiukajtis ix = hx & 0x7fffffff; 160*25c28e83SPiotr Jasiukajtis lx = ((int *)&x)[LOWORD]; 161*25c28e83SPiotr Jasiukajtis 162*25c28e83SPiotr Jasiukajtis /* subnormal,0,negative,inf,nan */ 163*25c28e83SPiotr Jasiukajtis if ((hx + 0x100000) < 0x200000) { 164*25c28e83SPiotr Jasiukajtis if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */ 165*25c28e83SPiotr Jasiukajtis return (x * x); 166*25c28e83SPiotr Jasiukajtis if (((hx << 1) | lx) == 0) /* zero */ 167*25c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, x, 18)); 168*25c28e83SPiotr Jasiukajtis if (hx < 0) /* negative */ 169*25c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, x, 19)); 170*25c28e83SPiotr Jasiukajtis if (((hx - 0x7ff00000) | lx) == 0) /* +inf */ 171*25c28e83SPiotr Jasiukajtis return (x); 172*25c28e83SPiotr Jasiukajtis 173*25c28e83SPiotr Jasiukajtis /* x must be positive and subnormal */ 174*25c28e83SPiotr Jasiukajtis x *= TWO52; 175*25c28e83SPiotr Jasiukajtis n = -52; 176*25c28e83SPiotr Jasiukajtis ix = ((int *)&x)[HIWORD]; 177*25c28e83SPiotr Jasiukajtis lx = ((int *)&x)[LOWORD]; 178*25c28e83SPiotr Jasiukajtis } 179*25c28e83SPiotr Jasiukajtis 180*25c28e83SPiotr Jasiukajtis i = ix >> 19; 181*25c28e83SPiotr Jasiukajtis if (i >= 0x7f7 && i <= 0x806) { 182*25c28e83SPiotr Jasiukajtis /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */ 183*25c28e83SPiotr Jasiukajtis if (ix >= 0x3fec0000 && ix < 0x3ff20000) { 184*25c28e83SPiotr Jasiukajtis /* 0.875 <= x < 1.125 */ 185*25c28e83SPiotr Jasiukajtis s = x - ONE; 186*25c28e83SPiotr Jasiukajtis z = s * s; 187*25c28e83SPiotr Jasiukajtis if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */ 188*25c28e83SPiotr Jasiukajtis return (z); 189*25c28e83SPiotr Jasiukajtis r = (A10 * s) * (A11 + s); 190*25c28e83SPiotr Jasiukajtis w = z * s; 191*25c28e83SPiotr Jasiukajtis return (LGH * s - (LGL * s - ((A1 * z) * 192*25c28e83SPiotr Jasiukajtis ((A2 + (A3 * s) * (A4 + s)) + w * (A5 + s))) * 193*25c28e83SPiotr Jasiukajtis (((A6 + s * (A7 + s)) + w * (A8 + s)) * 194*25c28e83SPiotr Jasiukajtis ((A9 + r) + w * (A12 + s))))); 195*25c28e83SPiotr Jasiukajtis } else { 196*25c28e83SPiotr Jasiukajtis i = (ix - 0x3fb80000) >> 15; 197*25c28e83SPiotr Jasiukajtis tb = (double *)_TBL_log + (i + i + i); 198*25c28e83SPiotr Jasiukajtis s = (x - tb[0]) * tb[1]; 199*25c28e83SPiotr Jasiukajtis return (LGH * tb[2] - (LGL * tb[2] - ((B1 * s) * 200*25c28e83SPiotr Jasiukajtis (B2 + s * (B3 + s))) * 201*25c28e83SPiotr Jasiukajtis (((B4 + s * B5) + (s * s) * (B6 + s)) * 202*25c28e83SPiotr Jasiukajtis (B7 + s * (B8 + s))))); 203*25c28e83SPiotr Jasiukajtis } 204*25c28e83SPiotr Jasiukajtis } else { 205*25c28e83SPiotr Jasiukajtis dn = (double)(n + ((ix >> 20) - 0x3ff)); 206*25c28e83SPiotr Jasiukajtis dn1 = dn * LNAHI; 207*25c28e83SPiotr Jasiukajtis i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */ 208*25c28e83SPiotr Jasiukajtis ((int *)&x)[HIWORD] = i; 209*25c28e83SPiotr Jasiukajtis i = (i - 0x3fb80000) >> 15; 210*25c28e83SPiotr Jasiukajtis tb = (double *)_TBL_log + (i + i + i); 211*25c28e83SPiotr Jasiukajtis s = (x - tb[0]) * tb[1]; 212*25c28e83SPiotr Jasiukajtis dn = dn * LNALO + tb[2] * LG10V; 213*25c28e83SPiotr Jasiukajtis return (dn1 + (dn + ((B1 * s) * 214*25c28e83SPiotr Jasiukajtis (B2 + s * (B3 + s))) * 215*25c28e83SPiotr Jasiukajtis (((B4 + s * B5) + (s * s) * (B6 + s)) * 216*25c28e83SPiotr Jasiukajtis (B7 + s * (B8 + s))))); 217*25c28e83SPiotr Jasiukajtis } 218*25c28e83SPiotr Jasiukajtis } 219