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 __log10l = log10l
31
32 /*
33 * log10l(X)
34 *
35 * Method :
36 * Let log10_2hi = leading 98(SPARC)/49(x86) bits of log10(2) and
37 * log10_2lo = log10(2) - log10_2hi,
38 * ivln10 = 1/log(10) rounded.
39 * Then
40 * n = ilogb(x),
41 * if (n<0) n = n+1;
42 * x = scalbn(x,-n);
43 * LOG10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
44 *
48 *
49 * Special cases:
50 * log10(x) is NaN with signal if x < 0;
51 * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
52 * log10(NaN) is that NaN with no signal;
53 * log10(10**N) = N for N=0,1,...,22.
54 *
55 * Constants:
56 * The hexadecimal values are the intended ones for the following constants.
57 * The decimal values may be used, provided that the compiler will convert
58 * from decimal to binary accurately enough to produce the hexadecimal values
59 * shown.
60 */
61
62 #include "libm.h"
63 #include "longdouble.h"
64
65 #if defined(__x86)
66 #define __swapRD __swap87RD
67 #endif
68 extern enum fp_direction_type __swapRD(enum fp_direction_type);
69
70 static const long double
71 zero = 0.0L,
72 ivln10 = 4.342944819032518276511289189166050822944e-0001L,
73 one = 1.0L,
74 #if defined(__x86)
75 log10_2hi = 3.010299956639803653501985536422580480576e-01L,
76 log10_2lo = 8.298635403410822349787106337291183585413e-16L;
77 #elif defined(__sparc)
78 log10_2hi = 3.010299956639811952137388947242098603469e-01L,
79 log10_2lo = 2.831664213089468167896664371953210945664e-31L;
80 #else
81 #error Unknown Architecture!
82 #endif
83
84 long double
85 log10l(long double x) {
86 long double y, z;
87 enum fp_direction_type rd;
88 int n;
89
90 if (!finitel(x))
91 return (x + fabsl(x)); /* x is +-INF or NaN */
92 else if (x > zero) {
93 n = ilogbl(x);
94 if (n < 0)
95 n += 1;
96 rd = __swapRD(fp_nearest);
97 y = n;
98 x = scalbnl(x, -n);
99 z = y * log10_2lo + ivln10 * logl(x);
100 z += y * log10_2hi;
101 if (rd != fp_nearest)
102 (void) __swapRD(rd);
103 return (z);
104 } else if (x == zero) /* -INF */
105 return (-one / zero);
106 else /* x <0, return NaN */
107 return (zero / zero);
108 }
|
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 /*
27 * Copyright 2006 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __log10l = log10l
32
33 /*
34 * log10l(X)
35 *
36 * Method :
37 * Let log10_2hi = leading 98(SPARC)/49(x86) bits of log10(2) and
38 * log10_2lo = log10(2) - log10_2hi,
39 * ivln10 = 1/log(10) rounded.
40 * Then
41 * n = ilogb(x),
42 * if (n<0) n = n+1;
43 * x = scalbn(x,-n);
44 * LOG10(x) := n*log10_2hi + (n*log10_2lo + ivln10*log(x))
45 *
49 *
50 * Special cases:
51 * log10(x) is NaN with signal if x < 0;
52 * log10(+INF) is +INF with no signal; log10(0) is -INF with signal;
53 * log10(NaN) is that NaN with no signal;
54 * log10(10**N) = N for N=0,1,...,22.
55 *
56 * Constants:
57 * The hexadecimal values are the intended ones for the following constants.
58 * The decimal values may be used, provided that the compiler will convert
59 * from decimal to binary accurately enough to produce the hexadecimal values
60 * shown.
61 */
62
63 #include "libm.h"
64 #include "longdouble.h"
65
66 #if defined(__x86)
67 #define __swapRD __swap87RD
68 #endif
69
70 extern enum fp_direction_type __swapRD(enum fp_direction_type);
71 static const long double zero = 0.0L,
72 ivln10 = 4.342944819032518276511289189166050822944e-0001L,
73 one = 1.0L,
74 #if defined(__x86)
75 log10_2hi = 3.010299956639803653501985536422580480576e-01L,
76 log10_2lo = 8.298635403410822349787106337291183585413e-16L;
77 #elif defined(__sparc)
78 log10_2hi = 3.010299956639811952137388947242098603469e-01L,
79 log10_2lo = 2.831664213089468167896664371953210945664e-31L;
80 #else
81 #error Unknown Architecture!
82 #endif
83
84 long double
85 log10l(long double x)
86 {
87 long double y, z;
88 enum fp_direction_type rd;
89 int n;
90
91 if (!finitel(x)) {
92 return (x + fabsl(x)); /* x is +-INF or NaN */
93 } else if (x > zero) {
94 n = ilogbl(x);
95
96 if (n < 0)
97 n += 1;
98
99 rd = __swapRD(fp_nearest);
100 y = n;
101 x = scalbnl(x, -n);
102 z = y * log10_2lo + ivln10 * logl(x);
103 z += y * log10_2hi;
104
105 if (rd != fp_nearest)
106 (void) __swapRD(rd);
107
108 return (z);
109 } else if (x == zero) { /* -INF */
110 return (-one / zero);
111 } else { /* x <0, return NaN */
112 return (zero / zero);
113 }
114 }
|