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 /*
23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
24*25c28e83SPiotr Jasiukajtis */
25*25c28e83SPiotr Jasiukajtis /*
26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms.
28*25c28e83SPiotr Jasiukajtis */
29*25c28e83SPiotr Jasiukajtis
30*25c28e83SPiotr Jasiukajtis #include "libm.h"
31*25c28e83SPiotr Jasiukajtis #include "longdouble.h"
32*25c28e83SPiotr Jasiukajtis
33*25c28e83SPiotr Jasiukajtis /*
34*25c28e83SPiotr Jasiukajtis * exp10l(x)
35*25c28e83SPiotr Jasiukajtis * n = nint(x*(log10/log2)) ;
36*25c28e83SPiotr Jasiukajtis * exp10(x) = 10**x = exp(x*ln(10)) = exp(n*ln2+(x*ln10-n*ln2))
37*25c28e83SPiotr Jasiukajtis * = 2**n*exp(ln10*(x-n*log2/log10)))
38*25c28e83SPiotr Jasiukajtis * If x is an integer <= M then use repeat multiplication. For
39*25c28e83SPiotr Jasiukajtis * 10**M is the largest representable integer, where
40*25c28e83SPiotr Jasiukajtis * M = 10 single precision (24 bits)
41*25c28e83SPiotr Jasiukajtis * M = 22 double precision (53 bits)
42*25c28e83SPiotr Jasiukajtis * M = 48 quadruple precision (113 bits)
43*25c28e83SPiotr Jasiukajtis */
44*25c28e83SPiotr Jasiukajtis
45*25c28e83SPiotr Jasiukajtis #define TINY 1.0e-20L /* single: 1e-5, double: 1e-10, quad: 1e-20 */
46*25c28e83SPiotr Jasiukajtis #define LG10OVT 4933.L /* single: 39, double: 309, quad: 4933 */
47*25c28e83SPiotr Jasiukajtis #define LG10UFT -4966.L /* single: -45, double: -323, quad: -4966 */
48*25c28e83SPiotr Jasiukajtis #define M 48
49*25c28e83SPiotr Jasiukajtis /* logt2hi : last 32 bits is zero for quad prec */
50*25c28e83SPiotr Jasiukajtis #define LOGT2HI 0.30102999566398119521373889472420986034688L
51*25c28e83SPiotr Jasiukajtis #define LOGT2LO 2.831664213089468167896664371953e-31L
52*25c28e83SPiotr Jasiukajtis
53*25c28e83SPiotr Jasiukajtis static const long double
54*25c28e83SPiotr Jasiukajtis zero = 0.0L,
55*25c28e83SPiotr Jasiukajtis tiny = TINY * TINY,
56*25c28e83SPiotr Jasiukajtis one = 1.0L,
57*25c28e83SPiotr Jasiukajtis lg10 = 3.321928094887362347870319429489390175865e+0000L,
58*25c28e83SPiotr Jasiukajtis ln10 = 2.302585092994045684017991454684364207601e+0000L,
59*25c28e83SPiotr Jasiukajtis logt2hi = LOGT2HI,
60*25c28e83SPiotr Jasiukajtis logt2lo = LOGT2LO,
61*25c28e83SPiotr Jasiukajtis lg10ovt = LG10OVT,
62*25c28e83SPiotr Jasiukajtis lg10uft = LG10UFT;
63*25c28e83SPiotr Jasiukajtis
64*25c28e83SPiotr Jasiukajtis long double
exp10l(long double x)65*25c28e83SPiotr Jasiukajtis exp10l(long double x) {
66*25c28e83SPiotr Jasiukajtis long double t, tenp;
67*25c28e83SPiotr Jasiukajtis int k;
68*25c28e83SPiotr Jasiukajtis
69*25c28e83SPiotr Jasiukajtis if (!finitel(x)) {
70*25c28e83SPiotr Jasiukajtis if (isnanl(x) || x > zero)
71*25c28e83SPiotr Jasiukajtis return (x + x);
72*25c28e83SPiotr Jasiukajtis else
73*25c28e83SPiotr Jasiukajtis return (zero);
74*25c28e83SPiotr Jasiukajtis }
75*25c28e83SPiotr Jasiukajtis if (fabsl(x) < tiny)
76*25c28e83SPiotr Jasiukajtis return (one + x);
77*25c28e83SPiotr Jasiukajtis if (x <= lg10ovt)
78*25c28e83SPiotr Jasiukajtis if (x >= lg10uft) {
79*25c28e83SPiotr Jasiukajtis k = (int) x;
80*25c28e83SPiotr Jasiukajtis tenp = 10.0L;
81*25c28e83SPiotr Jasiukajtis /* x is a small +integer */
82*25c28e83SPiotr Jasiukajtis if (0 <= k && k <= M && (long double) k == x) {
83*25c28e83SPiotr Jasiukajtis t = one;
84*25c28e83SPiotr Jasiukajtis if (k & 1)
85*25c28e83SPiotr Jasiukajtis t *= tenp;
86*25c28e83SPiotr Jasiukajtis k >>= 1;
87*25c28e83SPiotr Jasiukajtis while (k) {
88*25c28e83SPiotr Jasiukajtis tenp *= tenp;
89*25c28e83SPiotr Jasiukajtis if (k & 1)
90*25c28e83SPiotr Jasiukajtis t *= tenp;
91*25c28e83SPiotr Jasiukajtis k >>= 1;
92*25c28e83SPiotr Jasiukajtis }
93*25c28e83SPiotr Jasiukajtis return (t);
94*25c28e83SPiotr Jasiukajtis }
95*25c28e83SPiotr Jasiukajtis t = anintl(x * lg10);
96*25c28e83SPiotr Jasiukajtis return (scalbnl(expl(ln10 * ((x - t * logt2hi) -
97*25c28e83SPiotr Jasiukajtis t * logt2lo)), (int) t));
98*25c28e83SPiotr Jasiukajtis } else
99*25c28e83SPiotr Jasiukajtis return (scalbnl(one, -50000)); /* underflow */
100*25c28e83SPiotr Jasiukajtis else
101*25c28e83SPiotr Jasiukajtis return (scalbnl(one, 50000)); /* overflow */
102*25c28e83SPiotr Jasiukajtis }
103