xref: /illumos-gate/usr/src/lib/libm/common/Q/exp10l.c (revision 25c28e83beb90e7c80452a7c818c5e6f73a07dc8)
1 /*
2  * CDDL HEADER START
3  *
4  * The contents of this file are subject to the terms of the
5  * Common Development and Distribution License (the "License").
6  * You may not use this file except in compliance with the License.
7  *
8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
9  * or http://www.opensolaris.org/os/licensing.
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 exp10l = __exp10l
31 
32 #include "libm.h"
33 #include "longdouble.h"
34 
35 /*
36  * exp10l(x)
37  *	n = nint(x*(log10/log2)) ;
38  *	exp10(x) = 10**x = exp(x*ln(10)) = exp(n*ln2+(x*ln10-n*ln2))
39  *		 = 2**n*exp(ln10*(x-n*log2/log10)))
40  *	If x is an integer <= M then use repeat multiplication. For
41  *	10**M is the largest representable integer, where
42  *		M = 10		single precision (24 bits)
43  *		M = 22		double precision (53 bits)
44  *		M = 48		quadruple precision (113 bits)
45  */
46 
47 #define	TINY 	1.0e-20L	/* single: 1e-5, double: 1e-10, quad: 1e-20 */
48 #define	LG10OVT	4933.L		/* single:  39, double:  309, quad:  4933 */
49 #define	LG10UFT	-4966.L		/* single: -45, double: -323, quad: -4966 */
50 #define	M	48
51 			/* logt2hi : last 32 bits is zero for quad prec */
52 #define	LOGT2HI	0.30102999566398119521373889472420986034688L
53 #define	LOGT2LO	2.831664213089468167896664371953e-31L
54 
55 static const long double
56 	zero	= 0.0L,
57 	tiny	= TINY * TINY,
58 	one	= 1.0L,
59 	lg10	= 3.321928094887362347870319429489390175865e+0000L,
60 	ln10	= 2.302585092994045684017991454684364207601e+0000L,
61 	logt2hi	= LOGT2HI,
62 	logt2lo	= LOGT2LO,
63 	lg10ovt	= LG10OVT,
64 	lg10uft	= LG10UFT;
65 
66 long double
67 exp10l(long double x) {
68 	long double t, tenp;
69 	int k;
70 
71 	if (!finitel(x)) {
72 		if (isnanl(x) || x > zero)
73 			return (x + x);
74 		else
75 			return (zero);
76 	}
77 	if (fabsl(x) < tiny)
78 		return (one + x);
79 	if (x <= lg10ovt)
80 		if (x >= lg10uft) {
81 			k = (int) x;
82 			tenp = 10.0L;
83 					/* x is a small +integer */
84 			if (0 <= k && k <= M && (long double) k == x) {
85 				t = one;
86 				if (k & 1)
87 					t *= tenp;
88 				k >>= 1;
89 				while (k) {
90 					tenp *= tenp;
91 					if (k & 1)
92 						t *= tenp;
93 					k >>= 1;
94 				}
95 				return (t);
96 			}
97 			t = anintl(x * lg10);
98 			return (scalbnl(expl(ln10 * ((x - t * logt2hi) -
99 				t * logt2lo)), (int) t));
100 		} else
101 			return (scalbnl(one, -50000));	/* underflow */
102 	else
103 			return (scalbnl(one, 50000));	/* overflow  */
104 }
105