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 * Copyright 2011 Nexenta Systems, Inc. All rights reserved.
23 */
24 /*
25 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
26 * Use is subject to license terms.
27 */
28
29 #pragma weak __logf = logf
30
31 /*
32 * Algorithm:
33 *
34 * Let y = x rounded to six significant bits. Then for any choice
35 * of e and z such that y = 2^e z, we have
36 *
37 * log(x) = e log(2) + log(z) + log(1+(x-y)/y)
38 *
39 * Note that (x-y)/y = (x'-y')/y' for any scaled x' = sx, y' = sy;
40 * in particular, we can take s to be the power of two that makes
41 * ulp(x') = 1.
42 *
43 * From a table, obtain l = log(z) and r = 1/y'. For |s| <= 2^-6,
104 3.33368809981254554946e-01,
105 -5.00000008402474976565e-01
106 };
107
108 #define ln2 C[0]
109 #define K3 C[1]
110 #define K2 C[2]
111 #define K1 C[3]
112
113 float
114 logf(float x)
115 {
116 double v, t;
117 float f;
118 int hx, ix, i, exp, iy;
119
120 hx = *(int *)&x;
121 ix = hx & ~0x80000000;
122
123 if (ix >= 0x7f800000) /* nan or inf */
124 return ((hx < 0)? x * 0.0f : x * x);
125
126 exp = 0;
127 if (hx < 0x00800000) { /* negative, zero, or subnormal */
128 if (hx <= 0) {
129 f = 0.0f;
130 return ((ix == 0)? -1.0f / f : f / f);
131 }
132
133 /* subnormal; scale by 2^149 */
134 f = (float)ix;
135 ix = *(int *)&f;
136 exp = -149;
137 }
138
139 exp += (ix - 0x3f320000) >> 23;
140 ix &= 0x007fffff;
141 iy = (ix + 0x20000) & 0xfffc0000;
142 i = iy >> 17;
143 t = ln2 * (double)exp + TBL[i];
144 v = (double)(ix - iy) * TBL[i + 1];
145 v += (v * v) * (K1 + v * (K2 + v * K3));
146 f = (float)(t + v);
147 return (f);
148 }
|
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 /*
27 * Copyright 2005 Sun Microsystems, Inc. All rights reserved.
28 * Use is subject to license terms.
29 */
30
31 #pragma weak __logf = logf
32
33 /*
34 * Algorithm:
35 *
36 * Let y = x rounded to six significant bits. Then for any choice
37 * of e and z such that y = 2^e z, we have
38 *
39 * log(x) = e log(2) + log(z) + log(1+(x-y)/y)
40 *
41 * Note that (x-y)/y = (x'-y')/y' for any scaled x' = sx, y' = sy;
42 * in particular, we can take s to be the power of two that makes
43 * ulp(x') = 1.
44 *
45 * From a table, obtain l = log(z) and r = 1/y'. For |s| <= 2^-6,
106 3.33368809981254554946e-01,
107 -5.00000008402474976565e-01
108 };
109
110 #define ln2 C[0]
111 #define K3 C[1]
112 #define K2 C[2]
113 #define K1 C[3]
114
115 float
116 logf(float x)
117 {
118 double v, t;
119 float f;
120 int hx, ix, i, exp, iy;
121
122 hx = *(int *)&x;
123 ix = hx & ~0x80000000;
124
125 if (ix >= 0x7f800000) /* nan or inf */
126 return ((hx < 0) ? x * 0.0f : x *x);
127
128 exp = 0;
129
130 if (hx < 0x00800000) { /* negative, zero, or subnormal */
131 if (hx <= 0) {
132 f = 0.0f;
133 return ((ix == 0) ? -1.0f / f : f / f);
134 }
135
136 /* subnormal; scale by 2^149 */
137 f = (float)ix;
138 ix = *(int *)&f;
139 exp = -149;
140 }
141
142 exp += (ix - 0x3f320000) >> 23;
143 ix &= 0x007fffff;
144 iy = (ix + 0x20000) & 0xfffc0000;
145 i = iy >> 17;
146 t = ln2 * (double)exp + TBL[i];
147 v = (double)(ix - iy) * TBL[i + 1];
148 v += (v * v) * (K1 + v * (K2 + v * K3));
149 f = (float)(t + v);
150 return (f);
151 }
|