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 log2 = __log2
31
32 /* INDENT OFF */
33 /*
34 * log2(x) = log(x)/log2
35 *
36 * Base on Table look-up algorithm with product polynomial
37 * approximation for log(x).
38 *
39 * By K.C. Ng, Nov 29, 2004
40 *
41 * (a). For x in [1-0.125, 1+0.125], from log.c we have
42 * log(x) = f + ((a1*f^2) *
43 * ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
44 * (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
45 * ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
46 * where f = x - 1.
47 * (i) modify a1 <- a1 / log2
48 * (ii) 1/log2 = 1.4426950408889634...
49 * = 1.5 - 0.057304959... (4 bit shift)
50 * Let lv = 1.5 - 1/log2, then
71 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
72 * significant bits. Then
73 * log2(x) = n + log2(z).
74 *
75 * Special cases:
76 * log2(x) is NaN with signal if x < 0 (including -INF) ;
77 * log2(+INF) is +INF; log2(0) is -INF with signal;
78 * log2(NaN) is that NaN with no signal.
79 *
80 * Maximum error observed: less than 0.84 ulp
81 *
82 * Constants:
83 * The hexadecimal values are the intended ones for the following constants.
84 * The decimal values may be used, provided that the compiler will convert
85 * from decimal to binary accurately enough to produce the hexadecimal values
86 * shown.
87 */
88 /* INDENT ON */
89
90 #include "libm.h"
91 #include "libm_synonyms.h"
92 #include "libm_protos.h"
93
94 extern const double _TBL_log[];
95
96 static const double P[] = {
97 /* ONE */ 1.0,
98 /* TWO52 */ 4503599627370496.0,
99 /* LN10V */ 1.4426950408889634073599246810018920433347, /* 1/log10 */
100 /* ZERO */ 0.0,
101 /* A1 */ -9.6809362455249638217841932228967194640116e-02,
102 /* A2 */ 1.99628461483039965074226529395673424005508422852e+0000,
103 /* A3 */ 2.26812367662950720159642514772713184356689453125e+0000,
104 /* A4 */ -9.05030639084976384900471657601883634924888610840e-0001,
105 /* A5 */ -1.48275767132434044270894446526654064655303955078e+0000,
106 /* A6 */ 1.88158320939722756293122074566781520843505859375e+0000,
107 /* A7 */ 1.83309386046986411145098827546462416648864746094e+0000,
108 /* A8 */ 1.24847063988317086291601754055591300129890441895e+0000,
109 /* A9 */ 1.98372421445537705508854742220137268304824829102e+0000,
110 /* A10 */ -3.94711735767898475035764249696512706577777862549e-0001,
111 /* A11 */ 3.07890395362954372160402272129431366920471191406e+0000,
|
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 __log2 = log2
31
32 /* INDENT OFF */
33 /*
34 * log2(x) = log(x)/log2
35 *
36 * Base on Table look-up algorithm with product polynomial
37 * approximation for log(x).
38 *
39 * By K.C. Ng, Nov 29, 2004
40 *
41 * (a). For x in [1-0.125, 1+0.125], from log.c we have
42 * log(x) = f + ((a1*f^2) *
43 * ((a2 + (a3*f)*(a4+f)) + (f^3)*(a5+f))) *
44 * (((a6 + f*(a7+f)) + (f^3)*(a8+f)) *
45 * ((a9 + (a10*f)*(a11+f)) + (f^3)*(a12+f)))
46 * where f = x - 1.
47 * (i) modify a1 <- a1 / log2
48 * (ii) 1/log2 = 1.4426950408889634...
49 * = 1.5 - 0.057304959... (4 bit shift)
50 * Let lv = 1.5 - 1/log2, then
71 * z in [1,2). Then similar to (b) find a Y[i] that matches z to 5.5
72 * significant bits. Then
73 * log2(x) = n + log2(z).
74 *
75 * Special cases:
76 * log2(x) is NaN with signal if x < 0 (including -INF) ;
77 * log2(+INF) is +INF; log2(0) is -INF with signal;
78 * log2(NaN) is that NaN with no signal.
79 *
80 * Maximum error observed: less than 0.84 ulp
81 *
82 * Constants:
83 * The hexadecimal values are the intended ones for the following constants.
84 * The decimal values may be used, provided that the compiler will convert
85 * from decimal to binary accurately enough to produce the hexadecimal values
86 * shown.
87 */
88 /* INDENT ON */
89
90 #include "libm.h"
91 #include "libm_protos.h"
92
93 extern const double _TBL_log[];
94
95 static const double P[] = {
96 /* ONE */ 1.0,
97 /* TWO52 */ 4503599627370496.0,
98 /* LN10V */ 1.4426950408889634073599246810018920433347, /* 1/log10 */
99 /* ZERO */ 0.0,
100 /* A1 */ -9.6809362455249638217841932228967194640116e-02,
101 /* A2 */ 1.99628461483039965074226529395673424005508422852e+0000,
102 /* A3 */ 2.26812367662950720159642514772713184356689453125e+0000,
103 /* A4 */ -9.05030639084976384900471657601883634924888610840e-0001,
104 /* A5 */ -1.48275767132434044270894446526654064655303955078e+0000,
105 /* A6 */ 1.88158320939722756293122074566781520843505859375e+0000,
106 /* A7 */ 1.83309386046986411145098827546462416648864746094e+0000,
107 /* A8 */ 1.24847063988317086291601754055591300129890441895e+0000,
108 /* A9 */ 1.98372421445537705508854742220137268304824829102e+0000,
109 /* A10 */ -3.94711735767898475035764249696512706577777862549e-0001,
110 /* A11 */ 3.07890395362954372160402272129431366920471191406e+0000,
|