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 2006 Sun Microsystems, Inc. All rights reserved.
26*25c28e83SPiotr Jasiukajtis * Use is subject to license terms.
27*25c28e83SPiotr Jasiukajtis */
28*25c28e83SPiotr Jasiukajtis
29*25c28e83SPiotr Jasiukajtis /* INDENT OFF */
30*25c28e83SPiotr Jasiukajtis /*
31*25c28e83SPiotr Jasiukajtis * exp10(x)
32*25c28e83SPiotr Jasiukajtis * Code by K.C. Ng for SUN 4.0 libm.
33*25c28e83SPiotr Jasiukajtis * Method :
34*25c28e83SPiotr Jasiukajtis * n = nint(x*(log10/log2));
35*25c28e83SPiotr Jasiukajtis * exp10(x) = 10**x = exp(x*ln(10)) = exp(n*ln2+(x*ln10-n*ln2))
36*25c28e83SPiotr Jasiukajtis * = 2**n*exp(ln10*(x-n*log2/log10)))
37*25c28e83SPiotr Jasiukajtis * If x is an integer < 23 then use repeat multiplication. For
38*25c28e83SPiotr Jasiukajtis * 10**22 is the largest representable integer.
39*25c28e83SPiotr Jasiukajtis */
40*25c28e83SPiotr Jasiukajtis /* INDENT ON */
41*25c28e83SPiotr Jasiukajtis
42*25c28e83SPiotr Jasiukajtis #include "libm.h"
43*25c28e83SPiotr Jasiukajtis
44*25c28e83SPiotr Jasiukajtis static const double C[] = {
45*25c28e83SPiotr Jasiukajtis 3.3219280948736234787, /* log(10)/log(2) */
46*25c28e83SPiotr Jasiukajtis 2.3025850929940456840, /* log(10) */
47*25c28e83SPiotr Jasiukajtis 3.0102999565860955045E-1, /* log(2)/log(10) high */
48*25c28e83SPiotr Jasiukajtis 5.3716447674669983622E-12, /* log(2)/log(10) low */
49*25c28e83SPiotr Jasiukajtis 0.0,
50*25c28e83SPiotr Jasiukajtis 0.5,
51*25c28e83SPiotr Jasiukajtis 1.0,
52*25c28e83SPiotr Jasiukajtis 10.0,
53*25c28e83SPiotr Jasiukajtis 1.0e300,
54*25c28e83SPiotr Jasiukajtis 1.0e-300,
55*25c28e83SPiotr Jasiukajtis };
56*25c28e83SPiotr Jasiukajtis
57*25c28e83SPiotr Jasiukajtis #define lg10 C[0]
58*25c28e83SPiotr Jasiukajtis #define ln10 C[1]
59*25c28e83SPiotr Jasiukajtis #define logt2hi C[2]
60*25c28e83SPiotr Jasiukajtis #define logt2lo C[3]
61*25c28e83SPiotr Jasiukajtis #define zero C[4]
62*25c28e83SPiotr Jasiukajtis #define half C[5]
63*25c28e83SPiotr Jasiukajtis #define one C[6]
64*25c28e83SPiotr Jasiukajtis #define ten C[7]
65*25c28e83SPiotr Jasiukajtis #define huge C[8]
66*25c28e83SPiotr Jasiukajtis #define tiny C[9]
67*25c28e83SPiotr Jasiukajtis
68*25c28e83SPiotr Jasiukajtis double
exp10(double x)69*25c28e83SPiotr Jasiukajtis exp10(double x) {
70*25c28e83SPiotr Jasiukajtis double t, pt;
71*25c28e83SPiotr Jasiukajtis int ix, hx, k;
72*25c28e83SPiotr Jasiukajtis
73*25c28e83SPiotr Jasiukajtis ix = ((int *)&x)[HIWORD];
74*25c28e83SPiotr Jasiukajtis hx = ix & ~0x80000000;
75*25c28e83SPiotr Jasiukajtis
76*25c28e83SPiotr Jasiukajtis if (hx >= 0x4074a000) { /* |x| >= 330 or x is nan */
77*25c28e83SPiotr Jasiukajtis if (hx >= 0x7ff00000) { /* x is inf or nan */
78*25c28e83SPiotr Jasiukajtis if (ix == 0xfff00000 && ((int *)&x)[LOWORD] == 0)
79*25c28e83SPiotr Jasiukajtis return (zero);
80*25c28e83SPiotr Jasiukajtis return (x * x);
81*25c28e83SPiotr Jasiukajtis }
82*25c28e83SPiotr Jasiukajtis t = (ix < 0)? tiny : huge;
83*25c28e83SPiotr Jasiukajtis return (t * t);
84*25c28e83SPiotr Jasiukajtis }
85*25c28e83SPiotr Jasiukajtis
86*25c28e83SPiotr Jasiukajtis if (hx < 0x3c000000)
87*25c28e83SPiotr Jasiukajtis return (one + x);
88*25c28e83SPiotr Jasiukajtis
89*25c28e83SPiotr Jasiukajtis k = (int)x;
90*25c28e83SPiotr Jasiukajtis if (0 <= k && k < 23 && (double)k == x) {
91*25c28e83SPiotr Jasiukajtis /* x is a small positive integer */
92*25c28e83SPiotr Jasiukajtis t = one;
93*25c28e83SPiotr Jasiukajtis pt = ten;
94*25c28e83SPiotr Jasiukajtis if (k & 1)
95*25c28e83SPiotr Jasiukajtis t = ten;
96*25c28e83SPiotr Jasiukajtis k >>= 1;
97*25c28e83SPiotr Jasiukajtis while (k) {
98*25c28e83SPiotr Jasiukajtis pt *= pt;
99*25c28e83SPiotr Jasiukajtis if (k & 1)
100*25c28e83SPiotr Jasiukajtis t *= pt;
101*25c28e83SPiotr Jasiukajtis k >>= 1;
102*25c28e83SPiotr Jasiukajtis }
103*25c28e83SPiotr Jasiukajtis return (t);
104*25c28e83SPiotr Jasiukajtis }
105*25c28e83SPiotr Jasiukajtis t = x * lg10;
106*25c28e83SPiotr Jasiukajtis k = (int)((ix < 0)? t - half : t + half);
107*25c28e83SPiotr Jasiukajtis return (scalbn(exp(ln10 * ((x - k * logt2hi) - k * logt2lo)), k));
108*25c28e83SPiotr Jasiukajtis }
109