1 /* $NetBSD: fpu_log.c,v 1.18 2014/01/04 13:23:22 isaki Exp $ */
2
3 /*
4 * Copyright (c) 1995 Ken Nakata
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without
8 * modification, are permitted provided that the following conditions
9 * are met:
10 * 1. Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * 3. Neither the name of the author nor the names of its contributors
16 * may be used to endorse or promote products derived from this software
17 * without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
20 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
23 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
24 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
26 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
27 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
28 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
29 * SUCH DAMAGE.
30 *
31 * @(#)fpu_log.c 10/8/95
32 */
33
34 #include <sys/cdefs.h>
35 __KERNEL_RCSID(0, "$NetBSD: fpu_log.c,v 1.18 2014/01/04 13:23:22 isaki Exp $");
36
37 #include <sys/types.h>
38 #include <sys/systm.h>
39
40 #include "fpu_emulate.h"
41
42 static uint32_t logA6[] = { 0x3FC2499A, 0xB5E4040B };
43 static uint32_t logA5[] = { 0xBFC555B5, 0x848CB7DB };
44 static uint32_t logA4[] = { 0x3FC99999, 0x987D8730 };
45 static uint32_t logA3[] = { 0xBFCFFFFF, 0xFF6F7E97 };
46 static uint32_t logA2[] = { 0x3FD55555, 0x555555A4 };
47 static uint32_t logA1[] = { 0xBFE00000, 0x00000008 };
48
49 static uint32_t logB5[] = { 0x3F175496, 0xADD7DAD6 };
50 static uint32_t logB4[] = { 0x3F3C71C2, 0xFE80C7E0 };
51 static uint32_t logB3[] = { 0x3F624924, 0x928BCCFF };
52 static uint32_t logB2[] = { 0x3F899999, 0x999995EC };
53 static uint32_t logB1[] = { 0x3FB55555, 0x55555555 };
54
55 /* sfpn = shortened fp number; can represent only positive numbers */
56 static struct sfpn {
57 int sp_exp;
58 uint32_t sp_m0, sp_m1;
59 } logtbl[] = {
60 { 0x3FFE - 0x3fff, 0xFE03F80FU, 0xE03F80FEU },
61 { 0x3FF7 - 0x3fff, 0xFF015358U, 0x833C47E2U },
62 { 0x3FFE - 0x3fff, 0xFA232CF2U, 0x52138AC0U },
63 { 0x3FF9 - 0x3fff, 0xBDC8D83EU, 0xAD88D549U },
64 { 0x3FFE - 0x3fff, 0xF6603D98U, 0x0F6603DAU },
65 { 0x3FFA - 0x3fff, 0x9CF43DCFU, 0xF5EAFD48U },
66 { 0x3FFE - 0x3fff, 0xF2B9D648U, 0x0F2B9D65U },
67 { 0x3FFA - 0x3fff, 0xDA16EB88U, 0xCB8DF614U },
68 { 0x3FFE - 0x3fff, 0xEF2EB71FU, 0xC4345238U },
69 { 0x3FFB - 0x3fff, 0x8B29B775U, 0x1BD70743U },
70 { 0x3FFE - 0x3fff, 0xEBBDB2A5U, 0xC1619C8CU },
71 { 0x3FFB - 0x3fff, 0xA8D839F8U, 0x30C1FB49U },
72 { 0x3FFE - 0x3fff, 0xE865AC7BU, 0x7603A197U },
73 { 0x3FFB - 0x3fff, 0xC61A2EB1U, 0x8CD907ADU },
74 { 0x3FFE - 0x3fff, 0xE525982AU, 0xF70C880EU },
75 { 0x3FFB - 0x3fff, 0xE2F2A47AU, 0xDE3A18AFU },
76 { 0x3FFE - 0x3fff, 0xE1FC780EU, 0x1FC780E2U },
77 { 0x3FFB - 0x3fff, 0xFF64898EU, 0xDF55D551U },
78 { 0x3FFE - 0x3fff, 0xDEE95C4CU, 0xA037BA57U },
79 { 0x3FFC - 0x3fff, 0x8DB956A9U, 0x7B3D0148U },
80 { 0x3FFE - 0x3fff, 0xDBEB61EEU, 0xD19C5958U },
81 { 0x3FFC - 0x3fff, 0x9B8FE100U, 0xF47BA1DEU },
82 { 0x3FFE - 0x3fff, 0xD901B203U, 0x6406C80EU },
83 { 0x3FFC - 0x3fff, 0xA9372F1DU, 0x0DA1BD17U },
84 { 0x3FFE - 0x3fff, 0xD62B80D6U, 0x2B80D62CU },
85 { 0x3FFC - 0x3fff, 0xB6B07F38U, 0xCE90E46BU },
86 { 0x3FFE - 0x3fff, 0xD3680D36U, 0x80D3680DU },
87 { 0x3FFC - 0x3fff, 0xC3FD0329U, 0x06488481U },
88 { 0x3FFE - 0x3fff, 0xD0B69FCBU, 0xD2580D0BU },
89 { 0x3FFC - 0x3fff, 0xD11DE0FFU, 0x15AB18CAU },
90 { 0x3FFE - 0x3fff, 0xCE168A77U, 0x25080CE1U },
91 { 0x3FFC - 0x3fff, 0xDE1433A1U, 0x6C66B150U },
92 { 0x3FFE - 0x3fff, 0xCB8727C0U, 0x65C393E0U },
93 { 0x3FFC - 0x3fff, 0xEAE10B5AU, 0x7DDC8ADDU },
94 { 0x3FFE - 0x3fff, 0xC907DA4EU, 0x871146ADU },
95 { 0x3FFC - 0x3fff, 0xF7856E5EU, 0xE2C9B291U },
96 { 0x3FFE - 0x3fff, 0xC6980C69U, 0x80C6980CU },
97 { 0x3FFD - 0x3fff, 0x82012CA5U, 0xA68206D7U },
98 { 0x3FFE - 0x3fff, 0xC4372F85U, 0x5D824CA6U },
99 { 0x3FFD - 0x3fff, 0x882C5FCDU, 0x7256A8C5U },
100 { 0x3FFE - 0x3fff, 0xC1E4BBD5U, 0x95F6E947U },
101 { 0x3FFD - 0x3fff, 0x8E44C60BU, 0x4CCFD7DEU },
102 { 0x3FFE - 0x3fff, 0xBFA02FE8U, 0x0BFA02FFU },
103 { 0x3FFD - 0x3fff, 0x944AD09EU, 0xF4351AF6U },
104 { 0x3FFE - 0x3fff, 0xBD691047U, 0x07661AA3U },
105 { 0x3FFD - 0x3fff, 0x9A3EECD4U, 0xC3EAA6B2U },
106 { 0x3FFE - 0x3fff, 0xBB3EE721U, 0xA54D880CU },
107 { 0x3FFD - 0x3fff, 0xA0218434U, 0x353F1DE8U },
108 { 0x3FFE - 0x3fff, 0xB92143FAU, 0x36F5E02EU },
109 { 0x3FFD - 0x3fff, 0xA5F2FCABU, 0xBBC506DAU },
110 { 0x3FFE - 0x3fff, 0xB70FBB5AU, 0x19BE3659U },
111 { 0x3FFD - 0x3fff, 0xABB3B8BAU, 0x2AD362A5U },
112 { 0x3FFE - 0x3fff, 0xB509E68AU, 0x9B94821FU },
113 { 0x3FFD - 0x3fff, 0xB1641795U, 0xCE3CA97BU },
114 { 0x3FFE - 0x3fff, 0xB30F6352U, 0x8917C80BU },
115 { 0x3FFD - 0x3fff, 0xB7047551U, 0x5D0F1C61U },
116 { 0x3FFE - 0x3fff, 0xB11FD3B8U, 0x0B11FD3CU },
117 { 0x3FFD - 0x3fff, 0xBC952AFEU, 0xEA3D13E1U },
118 { 0x3FFE - 0x3fff, 0xAF3ADDC6U, 0x80AF3ADEU },
119 { 0x3FFD - 0x3fff, 0xC2168ED0U, 0xF458BA4AU },
120 { 0x3FFE - 0x3fff, 0xAD602B58U, 0x0AD602B6U },
121 { 0x3FFD - 0x3fff, 0xC788F439U, 0xB3163BF1U },
122 { 0x3FFE - 0x3fff, 0xAB8F69E2U, 0x8359CD11U },
123 { 0x3FFD - 0x3fff, 0xCCECAC08U, 0xBF04565DU },
124 { 0x3FFE - 0x3fff, 0xA9C84A47U, 0xA07F5638U },
125 { 0x3FFD - 0x3fff, 0xD2420487U, 0x2DD85160U },
126 { 0x3FFE - 0x3fff, 0xA80A80A8U, 0x0A80A80BU },
127 { 0x3FFD - 0x3fff, 0xD7894992U, 0x3BC3588AU },
128 { 0x3FFE - 0x3fff, 0xA655C439U, 0x2D7B73A8U },
129 { 0x3FFD - 0x3fff, 0xDCC2C4B4U, 0x9887DACCU },
130 { 0x3FFE - 0x3fff, 0xA4A9CF1DU, 0x96833751U },
131 { 0x3FFD - 0x3fff, 0xE1EEBD3EU, 0x6D6A6B9EU },
132 { 0x3FFE - 0x3fff, 0xA3065E3FU, 0xAE7CD0E0U },
133 { 0x3FFD - 0x3fff, 0xE70D785CU, 0x2F9F5BDCU },
134 { 0x3FFE - 0x3fff, 0xA16B312EU, 0xA8FC377DU },
135 { 0x3FFD - 0x3fff, 0xEC1F392CU, 0x5179F283U },
136 { 0x3FFE - 0x3fff, 0x9FD809FDU, 0x809FD80AU },
137 { 0x3FFD - 0x3fff, 0xF12440D3U, 0xE36130E6U },
138 { 0x3FFE - 0x3fff, 0x9E4CAD23U, 0xDD5F3A20U },
139 { 0x3FFD - 0x3fff, 0xF61CCE92U, 0x346600BBU },
140 { 0x3FFE - 0x3fff, 0x9CC8E160U, 0xC3FB19B9U },
141 { 0x3FFD - 0x3fff, 0xFB091FD3U, 0x8145630AU },
142 { 0x3FFE - 0x3fff, 0x9B4C6F9EU, 0xF03A3CAAU },
143 { 0x3FFD - 0x3fff, 0xFFE97042U, 0xBFA4C2ADU },
144 { 0x3FFE - 0x3fff, 0x99D722DAU, 0xBDE58F06U },
145 { 0x3FFE - 0x3fff, 0x825EFCEDU, 0x49369330U },
146 { 0x3FFE - 0x3fff, 0x9868C809U, 0x868C8098U },
147 { 0x3FFE - 0x3fff, 0x84C37A7AU, 0xB9A905C9U },
148 { 0x3FFE - 0x3fff, 0x97012E02U, 0x5C04B809U },
149 { 0x3FFE - 0x3fff, 0x87224C2EU, 0x8E645FB7U },
150 { 0x3FFE - 0x3fff, 0x95A02568U, 0x095A0257U },
151 { 0x3FFE - 0x3fff, 0x897B8CACU, 0x9F7DE298U },
152 { 0x3FFE - 0x3fff, 0x94458094U, 0x45809446U },
153 { 0x3FFE - 0x3fff, 0x8BCF55DEU, 0xC4CD05FEU },
154 { 0x3FFE - 0x3fff, 0x92F11384U, 0x0497889CU },
155 { 0x3FFE - 0x3fff, 0x8E1DC0FBU, 0x89E125E5U },
156 { 0x3FFE - 0x3fff, 0x91A2B3C4U, 0xD5E6F809U },
157 { 0x3FFE - 0x3fff, 0x9066E68CU, 0x955B6C9BU },
158 { 0x3FFE - 0x3fff, 0x905A3863U, 0x3E06C43BU },
159 { 0x3FFE - 0x3fff, 0x92AADE74U, 0xC7BE59E0U },
160 { 0x3FFE - 0x3fff, 0x8F1779D9U, 0xFDC3A219U },
161 { 0x3FFE - 0x3fff, 0x94E9BFF6U, 0x15845643U },
162 { 0x3FFE - 0x3fff, 0x8DDA5202U, 0x37694809U },
163 { 0x3FFE - 0x3fff, 0x9723A1B7U, 0x20134203U },
164 { 0x3FFE - 0x3fff, 0x8CA29C04U, 0x6514E023U },
165 { 0x3FFE - 0x3fff, 0x995899C8U, 0x90EB8990U },
166 { 0x3FFE - 0x3fff, 0x8B70344AU, 0x139BC75AU },
167 { 0x3FFE - 0x3fff, 0x9B88BDAAU, 0x3A3DAE2FU },
168 { 0x3FFE - 0x3fff, 0x8A42F870U, 0x5669DB46U },
169 { 0x3FFE - 0x3fff, 0x9DB4224FU, 0xFFE1157CU },
170 { 0x3FFE - 0x3fff, 0x891AC73AU, 0xE9819B50U },
171 { 0x3FFE - 0x3fff, 0x9FDADC26U, 0x8B7A12DAU },
172 { 0x3FFE - 0x3fff, 0x87F78087U, 0xF78087F8U },
173 { 0x3FFE - 0x3fff, 0xA1FCFF17U, 0xCE733BD4U },
174 { 0x3FFE - 0x3fff, 0x86D90544U, 0x7A34ACC6U },
175 { 0x3FFE - 0x3fff, 0xA41A9E8FU, 0x5446FB9FU },
176 { 0x3FFE - 0x3fff, 0x85BF3761U, 0x2CEE3C9BU },
177 { 0x3FFE - 0x3fff, 0xA633CD7EU, 0x6771CD8BU },
178 { 0x3FFE - 0x3fff, 0x84A9F9C8U, 0x084A9F9DU },
179 { 0x3FFE - 0x3fff, 0xA8489E60U, 0x0B435A5EU },
180 { 0x3FFE - 0x3fff, 0x83993052U, 0x3FBE3368U },
181 { 0x3FFE - 0x3fff, 0xAA59233CU, 0xCCA4BD49U },
182 { 0x3FFE - 0x3fff, 0x828CBFBEU, 0xB9A020A3U },
183 { 0x3FFE - 0x3fff, 0xAC656DAEU, 0x6BCC4985U },
184 { 0x3FFE - 0x3fff, 0x81848DA8U, 0xFAF0D277U },
185 { 0x3FFE - 0x3fff, 0xAE6D8EE3U, 0x60BB2468U },
186 { 0x3FFE - 0x3fff, 0x80808080U, 0x80808081U },
187 { 0x3FFE - 0x3fff, 0xB07197A2U, 0x3C46C654U },
188 };
189
190 static struct fpn *__fpu_logn(struct fpemu *fe);
191
192 /*
193 * natural log - algorithm taken from Motorola FPSP,
194 * except this doesn't bother to check for invalid input.
195 */
196 static struct fpn *
__fpu_logn(struct fpemu * fe)197 __fpu_logn(struct fpemu *fe)
198 {
199 static struct fpn X, F, U, V, W, KLOG2;
200 struct fpn *d;
201 int i, k;
202
203 CPYFPN(&X, &fe->fe_f2);
204
205 /* see if |X-1| < 1/16 approx. */
206 if ((-1 == X.fp_exp && (0xf07d0000U >> (31 - FP_LG)) <= X.fp_mant[0]) ||
207 (0 == X.fp_exp && X.fp_mant[0] <= (0x88410000U >> (31 - FP_LG)))) {
208 /* log near 1 */
209 #if FPE_DEBUG
210 printf("__fpu_logn: log near 1\n");
211 #endif
212
213 fpu_const(&fe->fe_f1, FPU_CONST_1);
214 /* X+1 */
215 d = fpu_add(fe);
216 CPYFPN(&V, d);
217
218 CPYFPN(&fe->fe_f1, &X);
219 fpu_const(&fe->fe_f2, FPU_CONST_1);
220 fe->fe_f2.fp_sign = 1; /* -1.0 */
221 /* X-1 */
222 d = fpu_add(fe);
223 CPYFPN(&fe->fe_f1, d);
224 /* 2(X-1) */
225 fe->fe_f1.fp_exp++; /* *= 2 */
226 CPYFPN(&fe->fe_f2, &V);
227 /* U=2(X-1)/(X+1) */
228 d = fpu_div(fe);
229 CPYFPN(&U, d);
230 CPYFPN(&fe->fe_f1, d);
231 CPYFPN(&fe->fe_f2, d);
232 /* V=U*U */
233 d = fpu_mul(fe);
234 CPYFPN(&V, d);
235 CPYFPN(&fe->fe_f1, d);
236 CPYFPN(&fe->fe_f2, d);
237 /* W=V*V */
238 d = fpu_mul(fe);
239 CPYFPN(&W, d);
240
241 /* calculate U+U*V*([B1+W*(B3+W*B5)]+[V*(B2+W*B4)]) */
242
243 /* B1+W*(B3+W*B5) part */
244 CPYFPN(&fe->fe_f1, d);
245 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB5);
246 /* W*B5 */
247 d = fpu_mul(fe);
248 CPYFPN(&fe->fe_f1, d);
249 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB3);
250 /* B3+W*B5 */
251 d = fpu_add(fe);
252 CPYFPN(&fe->fe_f1, d);
253 CPYFPN(&fe->fe_f2, &W);
254 /* W*(B3+W*B5) */
255 d = fpu_mul(fe);
256 CPYFPN(&fe->fe_f1, d);
257 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB1);
258 /* B1+W*(B3+W*B5) */
259 d = fpu_add(fe);
260 CPYFPN(&X, d);
261
262 /* [V*(B2+W*B4)] part */
263 CPYFPN(&fe->fe_f1, &W);
264 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB4);
265 /* W*B4 */
266 d = fpu_mul(fe);
267 CPYFPN(&fe->fe_f1, d);
268 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logB2);
269 /* B2+W*B4 */
270 d = fpu_add(fe);
271 CPYFPN(&fe->fe_f1, d);
272 CPYFPN(&fe->fe_f2, &V);
273 /* V*(B2+W*B4) */
274 d = fpu_mul(fe);
275 CPYFPN(&fe->fe_f1, d);
276 CPYFPN(&fe->fe_f2, &X);
277 /* B1+W*(B3+W*B5)+V*(B2+W*B4) */
278 d = fpu_add(fe);
279 CPYFPN(&fe->fe_f1, d);
280 CPYFPN(&fe->fe_f2, &V);
281 /* V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
282 d = fpu_mul(fe);
283 CPYFPN(&fe->fe_f1, d);
284 CPYFPN(&fe->fe_f2, &U);
285 /* U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
286 d = fpu_mul(fe);
287 CPYFPN(&fe->fe_f1, d);
288 CPYFPN(&fe->fe_f2, &U);
289 /* U+U*V*(B1+W*(B3+W*B5)+V*(B2+W*B4)) */
290 d = fpu_add(fe);
291 } else /* the usual case */ {
292 #if FPE_DEBUG
293 printf("__fpu_logn: the usual case. X=(%d,%08x,%08x...)\n",
294 X.fp_exp, X.fp_mant[0], X.fp_mant[1]);
295 #endif
296
297 k = X.fp_exp;
298 /* X <- Y */
299 X.fp_exp = fe->fe_f2.fp_exp = 0;
300
301 /* get the most significant 7 bits of X */
302 F.fp_class = FPC_NUM;
303 F.fp_sign = 0;
304 F.fp_exp = X.fp_exp;
305 F.fp_mant[0] = X.fp_mant[0] & (0xfe000000U >> (31 - FP_LG));
306 F.fp_mant[0] |= (0x01000000U >> (31 - FP_LG));
307 F.fp_mant[1] = F.fp_mant[2] = 0;
308 F.fp_sticky = 0;
309
310 #if FPE_DEBUG
311 printf("__fpu_logn: X=Y*2^k=(%d,%08x,%08x...)*2^%d\n",
312 fe->fe_f2.fp_exp, fe->fe_f2.fp_mant[0],
313 fe->fe_f2.fp_mant[1], k);
314 printf("__fpu_logn: F=(%d,%08x,%08x...)\n",
315 F.fp_exp, F.fp_mant[0], F.fp_mant[1]);
316 #endif
317
318 /* index to the table */
319 i = (F.fp_mant[0] >> (FP_LG - 7)) & 0x7e;
320
321 #if FPE_DEBUG
322 printf("__fpu_logn: index to logtbl i=%d(%x)\n", i, i);
323 #endif
324
325 CPYFPN(&fe->fe_f1, &F);
326 /* -F */
327 fe->fe_f1.fp_sign = 1;
328 /* Y-F */
329 d = fpu_add(fe);
330 CPYFPN(&fe->fe_f1, d);
331
332 /* fe_f2 = 1/F */
333 fe->fe_f2.fp_class = FPC_NUM;
334 fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[2]
335 = 0;
336 fe->fe_f2.fp_exp = logtbl[i].sp_exp;
337 fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG));
338 fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) |
339 (logtbl[i].sp_m1 >> (31 - FP_LG));
340 fe->fe_f2.fp_mant[2] =
341 (uint32_t)(logtbl[i].sp_m1 << (FP_LG + 1));
342
343 #if FPE_DEBUG
344 printf("__fpu_logn: 1/F=(%d,%08x,%08x...)\n", fe->fe_f2.fp_exp,
345 fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
346 #endif
347
348 /* U = (Y-F) * (1/F) */
349 d = fpu_mul(fe);
350 CPYFPN(&U, d);
351
352 /* KLOG2 = K * ln(2) */
353 /* fe_f1 == (fpn)k */
354 fpu_explode(fe, &fe->fe_f1, FTYPE_LNG, &k);
355 (void)fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
356 #if FPE_DEBUG
357 printf("__fpu_logn: fp(k)=(%d,%08x,%08x...)\n",
358 fe->fe_f1.fp_exp,
359 fe->fe_f1.fp_mant[0], fe->fe_f1.fp_mant[1]);
360 printf("__fpu_logn: ln(2)=(%d,%08x,%08x...)\n",
361 fe->fe_f2.fp_exp,
362 fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
363 #endif
364 /* K * LOGOF2 */
365 d = fpu_mul(fe);
366 CPYFPN(&KLOG2, d);
367
368 /* V=U*U */
369 CPYFPN(&fe->fe_f1, &U);
370 CPYFPN(&fe->fe_f2, &U);
371 d = fpu_mul(fe);
372 CPYFPN(&V, d);
373
374 /*
375 * approximation of LOG(1+U) by
376 * (U+V*(A1+V*(A3+V*A5)))+(U*V*(A2+V*(A4+V*A6)))
377 */
378
379 /* (U+V*(A1+V*(A3+V*A5))) part */
380 CPYFPN(&fe->fe_f1, d);
381 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA5);
382 /* V*A5 */
383 d = fpu_mul(fe);
384
385 CPYFPN(&fe->fe_f1, d);
386 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA3);
387 /* A3+V*A5 */
388 d = fpu_add(fe);
389
390 CPYFPN(&fe->fe_f1, d);
391 CPYFPN(&fe->fe_f2, &V);
392 /* V*(A3+V*A5) */
393 d = fpu_mul(fe);
394
395 CPYFPN(&fe->fe_f1, d);
396 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA1);
397 /* A1+V*(A3+V*A5) */
398 d = fpu_add(fe);
399
400 CPYFPN(&fe->fe_f1, d);
401 CPYFPN(&fe->fe_f2, &V);
402 /* V*(A1+V*(A3+V*A5)) */
403 d = fpu_mul(fe);
404
405 CPYFPN(&fe->fe_f1, d);
406 CPYFPN(&fe->fe_f2, &U);
407 /* U+V*(A1+V*(A3+V*A5)) */
408 d = fpu_add(fe);
409
410 CPYFPN(&X, d);
411
412 /* (U*V*(A2+V*(A4+V*A6))) part */
413 CPYFPN(&fe->fe_f1, &V);
414 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA6);
415 /* V*A6 */
416 d = fpu_mul(fe);
417 CPYFPN(&fe->fe_f1, d);
418 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA4);
419 /* A4+V*A6 */
420 d = fpu_add(fe);
421 CPYFPN(&fe->fe_f1, d);
422 CPYFPN(&fe->fe_f2, &V);
423 /* V*(A4+V*A6) */
424 d = fpu_mul(fe);
425 CPYFPN(&fe->fe_f1, d);
426 fpu_explode(fe, &fe->fe_f2, FTYPE_DBL, logA2);
427 /* A2+V*(A4+V*A6) */
428 d = fpu_add(fe);
429 CPYFPN(&fe->fe_f1, d);
430 CPYFPN(&fe->fe_f2, &V);
431 /* V*(A2+V*(A4+V*A6)) */
432 d = fpu_mul(fe);
433 CPYFPN(&fe->fe_f1, d);
434 CPYFPN(&fe->fe_f2, &U);
435 /* U*V*(A2+V*(A4+V*A6)) */
436 d = fpu_mul(fe);
437 CPYFPN(&fe->fe_f1, d);
438 i++;
439 /* fe_f2 = logtbl[i+1] (== LOG(F)) */
440 fe->fe_f2.fp_class = FPC_NUM;
441 fe->fe_f2.fp_sign = fe->fe_f2.fp_sticky = fe->fe_f2.fp_mant[2]
442 = 0;
443 fe->fe_f2.fp_exp = logtbl[i].sp_exp;
444 fe->fe_f2.fp_mant[0] = (logtbl[i].sp_m0 >> (31 - FP_LG));
445 fe->fe_f2.fp_mant[1] = (logtbl[i].sp_m0 << (FP_LG + 1)) |
446 (logtbl[i].sp_m1 >> (31 - FP_LG));
447 fe->fe_f2.fp_mant[2] = (logtbl[i].sp_m1 << (FP_LG + 1));
448
449 #if FPE_DEBUG
450 printf("__fpu_logn: ln(F)=(%d,%08x,%08x,...)\n",
451 fe->fe_f2.fp_exp,
452 fe->fe_f2.fp_mant[0], fe->fe_f2.fp_mant[1]);
453 #endif
454
455 /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
456 d = fpu_add(fe);
457 CPYFPN(&fe->fe_f1, d);
458 CPYFPN(&fe->fe_f2, &X);
459 /* LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */
460 d = fpu_add(fe);
461
462 #if FPE_DEBUG
463 printf("__fpu_logn: ln(Y)=(%c,%d,%08x,%08x,%08x)\n",
464 d->fp_sign ? '-' : '+', d->fp_exp,
465 d->fp_mant[0], d->fp_mant[1], d->fp_mant[2]);
466 #endif
467
468 CPYFPN(&fe->fe_f1, d);
469 CPYFPN(&fe->fe_f2, &KLOG2);
470 /* K*LOGOF2+LOG(F)+U+V*(A1+V*(A3+V*A5))+U*V*(A2+V*(A4+V*A6)) */
471 d = fpu_add(fe);
472 }
473
474 return d;
475 }
476
477 struct fpn *
fpu_log10(struct fpemu * fe)478 fpu_log10(struct fpemu *fe)
479 {
480 struct fpn *fp = &fe->fe_f2;
481 uint32_t fpsr;
482
483 fpsr = fe->fe_fpsr & ~FPSR_EXCP; /* clear all exceptions */
484
485 if (fp->fp_class >= FPC_NUM) {
486 if (fp->fp_sign) { /* negative number or Inf */
487 fp = fpu_newnan(fe);
488 fpsr |= FPSR_OPERR;
489 } else if (fp->fp_class == FPC_NUM) {
490 /* the real work here */
491 fp = __fpu_logn(fe);
492 if (fp != &fe->fe_f1)
493 CPYFPN(&fe->fe_f1, fp);
494 (void)fpu_const(&fe->fe_f2, FPU_CONST_LN_10);
495 fp = fpu_div(fe);
496 } /* else if fp == +Inf, return +Inf */
497 } else if (fp->fp_class == FPC_ZERO) {
498 /* return -Inf */
499 fp->fp_class = FPC_INF;
500 fp->fp_sign = 1;
501 fpsr |= FPSR_DZ;
502 } else if (fp->fp_class == FPC_SNAN) {
503 fpsr |= FPSR_SNAN;
504 fp = fpu_newnan(fe);
505 } else {
506 fp = fpu_newnan(fe);
507 }
508
509 fe->fe_fpsr = fpsr;
510
511 return fp;
512 }
513
514 struct fpn *
fpu_log2(struct fpemu * fe)515 fpu_log2(struct fpemu *fe)
516 {
517 struct fpn *fp = &fe->fe_f2;
518 uint32_t fpsr;
519
520 fpsr = fe->fe_fpsr & ~FPSR_EXCP; /* clear all exceptions */
521
522 if (fp->fp_class >= FPC_NUM) {
523 if (fp->fp_sign) { /* negative number or Inf */
524 fp = fpu_newnan(fe);
525 fpsr |= FPSR_OPERR;
526 } else if (fp->fp_class == FPC_NUM) {
527 /* the real work here */
528 if (fp->fp_mant[0] == FP_1 && fp->fp_mant[1] == 0 &&
529 fp->fp_mant[2] == 0) {
530 /* fp == 2.0 ^ exp <--> log2(fp) == exp */
531 fpu_explode(fe, &fe->fe_f3, FTYPE_LNG,
532 &fp->fp_exp);
533 fp = &fe->fe_f3;
534 } else {
535 fp = __fpu_logn(fe);
536 if (fp != &fe->fe_f1)
537 CPYFPN(&fe->fe_f1, fp);
538 (void)fpu_const(&fe->fe_f2, FPU_CONST_LN_2);
539 fp = fpu_div(fe);
540 }
541 } /* else if fp == +Inf, return +Inf */
542 } else if (fp->fp_class == FPC_ZERO) {
543 /* return -Inf */
544 fp->fp_class = FPC_INF;
545 fp->fp_sign = 1;
546 fpsr |= FPSR_DZ;
547 } else if (fp->fp_class == FPC_SNAN) {
548 fpsr |= FPSR_SNAN;
549 fp = fpu_newnan(fe);
550 } else {
551 fp = fpu_newnan(fe);
552 }
553
554 fe->fe_fpsr = fpsr;
555 return fp;
556 }
557
558 struct fpn *
fpu_logn(struct fpemu * fe)559 fpu_logn(struct fpemu *fe)
560 {
561 struct fpn *fp = &fe->fe_f2;
562 uint32_t fpsr;
563
564 fpsr = fe->fe_fpsr & ~FPSR_EXCP; /* clear all exceptions */
565
566 if (fp->fp_class >= FPC_NUM) {
567 if (fp->fp_sign) { /* negative number or Inf */
568 fp = fpu_newnan(fe);
569 fpsr |= FPSR_OPERR;
570 } else if (fp->fp_class == FPC_NUM) {
571 /* the real work here */
572 fp = __fpu_logn(fe);
573 } /* else if fp == +Inf, return +Inf */
574 } else if (fp->fp_class == FPC_ZERO) {
575 /* return -Inf */
576 fp->fp_class = FPC_INF;
577 fp->fp_sign = 1;
578 fpsr |= FPSR_DZ;
579 } else if (fp->fp_class == FPC_SNAN) {
580 fpsr |= FPSR_SNAN;
581 fp = fpu_newnan(fe);
582 } else {
583 fp = fpu_newnan(fe);
584 }
585
586 fe->fe_fpsr = fpsr;
587
588 return fp;
589 }
590
591 struct fpn *
fpu_lognp1(struct fpemu * fe)592 fpu_lognp1(struct fpemu *fe)
593 {
594 struct fpn *fp;
595
596 /* if src is +0/-0, return +0/-0 */
597 if (ISZERO(&fe->fe_f2))
598 return &fe->fe_f2;
599
600 /* build a 1.0 */
601 fp = fpu_const(&fe->fe_f1, FPU_CONST_1);
602 /* fp = 1.0 + f2 */
603 fp = fpu_add(fe);
604
605 /* copy the result to the src opr */
606 CPYFPN(&fe->fe_f2, fp);
607
608 return fpu_logn(fe);
609 }
610