125c28e83SPiotr Jasiukajtis /*
225c28e83SPiotr Jasiukajtis * CDDL HEADER START
325c28e83SPiotr Jasiukajtis *
425c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the
525c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License").
625c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License.
725c28e83SPiotr Jasiukajtis *
825c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
925c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing.
1025c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions
1125c28e83SPiotr Jasiukajtis * and limitations under the License.
1225c28e83SPiotr Jasiukajtis *
1325c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each
1425c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
1525c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the
1625c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying
1725c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner]
1825c28e83SPiotr Jasiukajtis *
1925c28e83SPiotr Jasiukajtis * CDDL HEADER END
2025c28e83SPiotr Jasiukajtis */
2125c28e83SPiotr Jasiukajtis /*
2225c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
2325c28e83SPiotr Jasiukajtis */
2425c28e83SPiotr Jasiukajtis /*
2525c28e83SPiotr Jasiukajtis * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
2625c28e83SPiotr Jasiukajtis * Use is subject to license terms.
2725c28e83SPiotr Jasiukajtis */
2825c28e83SPiotr Jasiukajtis
29*ddc0e0b5SRichard Lowe #pragma weak __log = log
3025c28e83SPiotr Jasiukajtis
3125c28e83SPiotr Jasiukajtis /* INDENT OFF */
3225c28e83SPiotr Jasiukajtis /*
3325c28e83SPiotr Jasiukajtis * log(x)
3425c28e83SPiotr Jasiukajtis * Table look-up algorithm with product polynomial approximation.
3525c28e83SPiotr Jasiukajtis * By K.C. Ng, Oct 23, 2004. Updated Oct 18, 2005.
3625c28e83SPiotr Jasiukajtis *
3725c28e83SPiotr Jasiukajtis * (a). For x in [1-0.125, 1+0.1328125], using a special approximation:
3825c28e83SPiotr Jasiukajtis * Let f = x - 1 and z = f*f.
3925c28e83SPiotr Jasiukajtis * return f + ((a1*z) *
4025c28e83SPiotr Jasiukajtis * ((a2 + (a3*f)*(a4+f)) + (f*z)*(a5+f))) *
4125c28e83SPiotr Jasiukajtis * (((a6 + f*(a7+f)) + (f*z)*(a8+f)) *
4225c28e83SPiotr Jasiukajtis * ((a9 + (a10*f)*(a11+f)) + (f*z)*(a12+f)))
4325c28e83SPiotr Jasiukajtis * a1 -6.88821452420390473170286327331268694251775741577e-0002,
4425c28e83SPiotr Jasiukajtis * a2 1.97493380704769294631262255279580131173133850098e+0000,
4525c28e83SPiotr Jasiukajtis * a3 2.24963218866067560242072431719861924648284912109e+0000,
4625c28e83SPiotr Jasiukajtis * a4 -9.02975906958474405783476868236903101205825805664e-0001,
4725c28e83SPiotr Jasiukajtis * a5 -1.47391630715542865104339398385491222143173217773e+0000,
4825c28e83SPiotr Jasiukajtis * a6 1.86846544648220058704168877738993614912033081055e+0000,
4925c28e83SPiotr Jasiukajtis * a7 1.82277370459347465292410106485476717352867126465e+0000,
5025c28e83SPiotr Jasiukajtis * a8 1.25295479915214102994980294170090928673744201660e+0000,
5125c28e83SPiotr Jasiukajtis * a9 1.96709676945198275177517643896862864494323730469e+0000,
5225c28e83SPiotr Jasiukajtis * a10 -4.00127989749189894030934055990655906498432159424e-0001,
5325c28e83SPiotr Jasiukajtis * a11 3.01675528558798333733648178167641162872314453125e+0000,
5425c28e83SPiotr Jasiukajtis * a12 -9.52325445049240770778453679668018594384193420410e-0001,
5525c28e83SPiotr Jasiukajtis *
5625c28e83SPiotr Jasiukajtis * with remez error |(log(1+f) - P(f))/f| <= 2**-56.81 and
5725c28e83SPiotr Jasiukajtis *
5825c28e83SPiotr Jasiukajtis * (b). For 0.09375 <= x < 24
5925c28e83SPiotr Jasiukajtis * Use an 8-bit table look-up (3-bit for exponent and 5 bit for
6025c28e83SPiotr Jasiukajtis * significand):
6125c28e83SPiotr Jasiukajtis * Let ix stands for the high part of x in IEEE double format.
6225c28e83SPiotr Jasiukajtis * Since 0.09375 <= x < 24, we have
6325c28e83SPiotr Jasiukajtis * 0x3fb80000 <= ix < 0x40380000.
6425c28e83SPiotr Jasiukajtis * Let j = (ix - 0x3fb80000) >> 15. Then 0 <= j < 256. Choose
6525c28e83SPiotr Jasiukajtis * a Y[j] such that HIWORD(Y[j]) ~ 0x3fb8400 + (j<<15) (the middle
6625c28e83SPiotr Jasiukajtis * number between 0x3fb80000 + (j<<15) and 3fb80000 + ((j+1)<<15)),
6725c28e83SPiotr Jasiukajtis * and at the same time 1/Y[j] as well as log(Y[j]) are very close
6825c28e83SPiotr Jasiukajtis * to 53-bits floating point numbers.
6925c28e83SPiotr Jasiukajtis * A table of Y[j], 1/Y[j], and log(Y[j]) are pre-computed and thus
7025c28e83SPiotr Jasiukajtis * log(x) = log(Y[j]) + log(1 + (x-Y[j])*(1/Y[j]))
7125c28e83SPiotr Jasiukajtis * = log(Y[j]) + log(1 + s)
7225c28e83SPiotr Jasiukajtis * where
7325c28e83SPiotr Jasiukajtis * s = (x-Y[j])*(1/Y[j])
7425c28e83SPiotr Jasiukajtis * We compute max (x-Y[j])*(1/Y[j]) for the chosen Y[j] and obtain
7525c28e83SPiotr Jasiukajtis * |s| < 0.0154. By applying remez algorithm with Product Polynomial
7625c28e83SPiotr Jasiukajtis * Approximiation, we find the following approximated of log(1+s)
7725c28e83SPiotr Jasiukajtis * (b1*s)*(b2+s*(b3+s))*((b4+s*b5)+(s*s)*(b6+s))*(b7+s*(b8+s))
7825c28e83SPiotr Jasiukajtis * with remez error |log(1+s) - P(s)| <= 2**-63.5
7925c28e83SPiotr Jasiukajtis *
8025c28e83SPiotr Jasiukajtis * (c). Otherwise, get "n", the exponent of x, and then normalize x to
8125c28e83SPiotr Jasiukajtis * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
8225c28e83SPiotr Jasiukajtis * significant bits. Then
8325c28e83SPiotr Jasiukajtis * log(x) = n*ln2 + log(Y[i]) + log(z/Y[i]).
8425c28e83SPiotr Jasiukajtis *
8525c28e83SPiotr Jasiukajtis * Special cases:
8625c28e83SPiotr Jasiukajtis * log(x) is NaN with signal if x < 0 (including -INF) ;
8725c28e83SPiotr Jasiukajtis * log(+INF) is +INF; log(0) is -INF with signal;
8825c28e83SPiotr Jasiukajtis * log(NaN) is that NaN with no signal.
8925c28e83SPiotr Jasiukajtis *
9025c28e83SPiotr Jasiukajtis * Maximum error observed: less than 0.90 ulp
9125c28e83SPiotr Jasiukajtis *
9225c28e83SPiotr Jasiukajtis * Constants:
9325c28e83SPiotr Jasiukajtis * The hexadecimal values are the intended ones for the following constants.
9425c28e83SPiotr Jasiukajtis * The decimal values may be used, provided that the compiler will convert
9525c28e83SPiotr Jasiukajtis * from decimal to binary accurately enough to produce the hexadecimal values
9625c28e83SPiotr Jasiukajtis * shown.
9725c28e83SPiotr Jasiukajtis */
9825c28e83SPiotr Jasiukajtis /* INDENT ON */
9925c28e83SPiotr Jasiukajtis
10025c28e83SPiotr Jasiukajtis #include "libm.h"
10125c28e83SPiotr Jasiukajtis
10225c28e83SPiotr Jasiukajtis extern const double _TBL_log[];
10325c28e83SPiotr Jasiukajtis
10425c28e83SPiotr Jasiukajtis static const double P[] = {
10525c28e83SPiotr Jasiukajtis /* ONE */ 1.0,
10625c28e83SPiotr Jasiukajtis /* TWO52 */ 4503599627370496.0,
10725c28e83SPiotr Jasiukajtis /* LN2HI */ 6.93147180369123816490e-01, /* 3fe62e42, fee00000 */
10825c28e83SPiotr Jasiukajtis /* LN2LO */ 1.90821492927058770002e-10, /* 3dea39ef, 35793c76 */
10925c28e83SPiotr Jasiukajtis /* A1 */ -6.88821452420390473170286327331268694251775741577e-0002,
11025c28e83SPiotr Jasiukajtis /* A2 */ 1.97493380704769294631262255279580131173133850098e+0000,
11125c28e83SPiotr Jasiukajtis /* A3 */ 2.24963218866067560242072431719861924648284912109e+0000,
11225c28e83SPiotr Jasiukajtis /* A4 */ -9.02975906958474405783476868236903101205825805664e-0001,
11325c28e83SPiotr Jasiukajtis /* A5 */ -1.47391630715542865104339398385491222143173217773e+0000,
11425c28e83SPiotr Jasiukajtis /* A6 */ 1.86846544648220058704168877738993614912033081055e+0000,
11525c28e83SPiotr Jasiukajtis /* A7 */ 1.82277370459347465292410106485476717352867126465e+0000,
11625c28e83SPiotr Jasiukajtis /* A8 */ 1.25295479915214102994980294170090928673744201660e+0000,
11725c28e83SPiotr Jasiukajtis /* A9 */ 1.96709676945198275177517643896862864494323730469e+0000,
11825c28e83SPiotr Jasiukajtis /* A10 */ -4.00127989749189894030934055990655906498432159424e-0001,
11925c28e83SPiotr Jasiukajtis /* A11 */ 3.01675528558798333733648178167641162872314453125e+0000,
12025c28e83SPiotr Jasiukajtis /* A12 */ -9.52325445049240770778453679668018594384193420410e-0001,
12125c28e83SPiotr Jasiukajtis /* B1 */ -1.25041641589283658575482149899471551179885864258e-0001,
12225c28e83SPiotr Jasiukajtis /* B2 */ 1.87161713283355151891381127914642725337613123482e+0000,
12325c28e83SPiotr Jasiukajtis /* B3 */ -1.89082956295731507978530316904652863740921020508e+0000,
12425c28e83SPiotr Jasiukajtis /* B4 */ -2.50562891673640253387134180229622870683670043945e+0000,
12525c28e83SPiotr Jasiukajtis /* B5 */ 1.64822828085258366037635369139024987816810607910e+0000,
12625c28e83SPiotr Jasiukajtis /* B6 */ -1.24409107065868340669112512841820716857910156250e+0000,
12725c28e83SPiotr Jasiukajtis /* B7 */ 1.70534231658220414296067701798165217041969299316e+0000,
12825c28e83SPiotr Jasiukajtis /* B8 */ 1.99196833784655646937267192697618156671524047852e+0000,
12925c28e83SPiotr Jasiukajtis };
13025c28e83SPiotr Jasiukajtis
13125c28e83SPiotr Jasiukajtis #define ONE P[0]
13225c28e83SPiotr Jasiukajtis #define TWO52 P[1]
13325c28e83SPiotr Jasiukajtis #define LN2HI P[2]
13425c28e83SPiotr Jasiukajtis #define LN2LO P[3]
13525c28e83SPiotr Jasiukajtis #define A1 P[4]
13625c28e83SPiotr Jasiukajtis #define A2 P[5]
13725c28e83SPiotr Jasiukajtis #define A3 P[6]
13825c28e83SPiotr Jasiukajtis #define A4 P[7]
13925c28e83SPiotr Jasiukajtis #define A5 P[8]
14025c28e83SPiotr Jasiukajtis #define A6 P[9]
14125c28e83SPiotr Jasiukajtis #define A7 P[10]
14225c28e83SPiotr Jasiukajtis #define A8 P[11]
14325c28e83SPiotr Jasiukajtis #define A9 P[12]
14425c28e83SPiotr Jasiukajtis #define A10 P[13]
14525c28e83SPiotr Jasiukajtis #define A11 P[14]
14625c28e83SPiotr Jasiukajtis #define A12 P[15]
14725c28e83SPiotr Jasiukajtis #define B1 P[16]
14825c28e83SPiotr Jasiukajtis #define B2 P[17]
14925c28e83SPiotr Jasiukajtis #define B3 P[18]
15025c28e83SPiotr Jasiukajtis #define B4 P[19]
15125c28e83SPiotr Jasiukajtis #define B5 P[20]
15225c28e83SPiotr Jasiukajtis #define B6 P[21]
15325c28e83SPiotr Jasiukajtis #define B7 P[22]
15425c28e83SPiotr Jasiukajtis #define B8 P[23]
15525c28e83SPiotr Jasiukajtis
15625c28e83SPiotr Jasiukajtis double
log(double x)15725c28e83SPiotr Jasiukajtis log(double x) {
15825c28e83SPiotr Jasiukajtis double *tb, dn, dn1, s, z, r, w;
15925c28e83SPiotr Jasiukajtis int i, hx, ix, n, lx;
16025c28e83SPiotr Jasiukajtis
16125c28e83SPiotr Jasiukajtis n = 0;
16225c28e83SPiotr Jasiukajtis hx = ((int *)&x)[HIWORD];
16325c28e83SPiotr Jasiukajtis ix = hx & 0x7fffffff;
16425c28e83SPiotr Jasiukajtis lx = ((int *)&x)[LOWORD];
16525c28e83SPiotr Jasiukajtis
16625c28e83SPiotr Jasiukajtis /* subnormal,0,negative,inf,nan */
16725c28e83SPiotr Jasiukajtis if ((hx + 0x100000) < 0x200000) {
16825c28e83SPiotr Jasiukajtis if (ix > 0x7ff00000 || (ix == 0x7ff00000 && lx != 0)) /* nan */
16925c28e83SPiotr Jasiukajtis return (x * x);
17025c28e83SPiotr Jasiukajtis if (((hx << 1) | lx) == 0) /* zero */
17125c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, x, 16));
17225c28e83SPiotr Jasiukajtis if (hx < 0) /* negative */
17325c28e83SPiotr Jasiukajtis return (_SVID_libm_err(x, x, 17));
17425c28e83SPiotr Jasiukajtis if (((hx - 0x7ff00000) | lx) == 0) /* +inf */
17525c28e83SPiotr Jasiukajtis return (x);
17625c28e83SPiotr Jasiukajtis
17725c28e83SPiotr Jasiukajtis /* x must be positive and subnormal */
17825c28e83SPiotr Jasiukajtis x *= TWO52;
17925c28e83SPiotr Jasiukajtis n = -52;
18025c28e83SPiotr Jasiukajtis ix = ((int *)&x)[HIWORD];
18125c28e83SPiotr Jasiukajtis lx = ((int *)&x)[LOWORD];
18225c28e83SPiotr Jasiukajtis }
18325c28e83SPiotr Jasiukajtis
18425c28e83SPiotr Jasiukajtis i = ix >> 19;
18525c28e83SPiotr Jasiukajtis if (i >= 0x7f7 && i <= 0x806) {
18625c28e83SPiotr Jasiukajtis /* 0.09375 (0x3fb80000) <= x < 24 (0x40380000) */
18725c28e83SPiotr Jasiukajtis if (ix >= 0x3fec0000 && ix < 0x3ff22000) {
18825c28e83SPiotr Jasiukajtis /* 0.875 <= x < 1.125 */
18925c28e83SPiotr Jasiukajtis s = x - ONE;
19025c28e83SPiotr Jasiukajtis z = s * s;
19125c28e83SPiotr Jasiukajtis if (((ix - 0x3ff00000) | lx) == 0) /* x = 1 */
19225c28e83SPiotr Jasiukajtis return (z);
19325c28e83SPiotr Jasiukajtis r = (A10 * s) * (A11 + s);
19425c28e83SPiotr Jasiukajtis w = z * s;
19525c28e83SPiotr Jasiukajtis return (s + ((A1 * z) *
19625c28e83SPiotr Jasiukajtis (A2 + ((A3 * s) * (A4 + s) + w * (A5 + s)))) *
19725c28e83SPiotr Jasiukajtis ((A6 + (s * (A7 + s) + w * (A8 + s))) *
19825c28e83SPiotr Jasiukajtis (A9 + (r + w * (A12 + s)))));
19925c28e83SPiotr Jasiukajtis } else {
20025c28e83SPiotr Jasiukajtis i = (ix - 0x3fb80000) >> 15;
20125c28e83SPiotr Jasiukajtis tb = (double *)_TBL_log + (i + i + i);
20225c28e83SPiotr Jasiukajtis s = (x - tb[0]) * tb[1];
20325c28e83SPiotr Jasiukajtis return (tb[2] + ((B1 * s) * (B2 + s * (B3 + s))) *
20425c28e83SPiotr Jasiukajtis (((B4 + s * B5) + (s * s) * (B6 + s)) *
20525c28e83SPiotr Jasiukajtis (B7 + s * (B8 + s))));
20625c28e83SPiotr Jasiukajtis }
20725c28e83SPiotr Jasiukajtis } else {
20825c28e83SPiotr Jasiukajtis dn = (double)(n + ((ix >> 20) - 0x3ff));
20925c28e83SPiotr Jasiukajtis dn1 = dn * LN2HI;
21025c28e83SPiotr Jasiukajtis i = (ix & 0x000fffff) | 0x3ff00000; /* scale x to [1,2] */
21125c28e83SPiotr Jasiukajtis ((int *)&x)[HIWORD] = i;
21225c28e83SPiotr Jasiukajtis i = (i - 0x3fb80000) >> 15;
21325c28e83SPiotr Jasiukajtis tb = (double *)_TBL_log + (i + i + i);
21425c28e83SPiotr Jasiukajtis s = (x - tb[0]) * tb[1];
21525c28e83SPiotr Jasiukajtis dn = dn * LN2LO + tb[2];
21625c28e83SPiotr Jasiukajtis return (dn1 + (dn + ((B1 * s) * (B2 + s * (B3 + s))) *
21725c28e83SPiotr Jasiukajtis (((B4 + s * B5) + (s * s) * (B6 + s)) *
21825c28e83SPiotr Jasiukajtis (B7 + s * (B8 + s)))));
21925c28e83SPiotr Jasiukajtis }
22025c28e83SPiotr Jasiukajtis }
221