?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
target/m68k: implement fasin
[qemu/ar7.git] / target / m68k / softfloat.c
blob91f0435ac36dd9c09cd7f09c04e747b1b01faa18
1 /*
2 * Ported from a work by Andreas Grabher for Previous, NeXT Computer Emulator,
3 * derived from NetBSD M68040 FPSP functions,
4 * derived from release 2a of the SoftFloat IEC/IEEE Floating-point Arithmetic
5 * Package. Those parts of the code (and some later contributions) are
6 * provided under that license, as detailed below.
7 * It has subsequently been modified by contributors to the QEMU Project,
8 * so some portions are provided under:
9 * the SoftFloat-2a license
10 * the BSD license
11 * GPL-v2-or-later
12 *
13 * Any future contributions to this file will be taken to be licensed under
14 * the Softfloat-2a license unless specifically indicated otherwise.
15 */
16
17 /* Portions of this work are licensed under the terms of the GNU GPL,
18 * version 2 or later. See the COPYING file in the top-level directory.
19 */
20
21 #include "qemu/osdep.h"
22 #include "softfloat.h"
23 #include "fpu/softfloat-macros.h"
24 #include "softfloat_fpsp_tables.h"
25
26 #define piby2_exp 0x3FFF
27 #define pi_sig LIT64(0xc90fdaa22168c235)
28
29 static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status)
30 {
31 if (floatx80_is_signaling_nan(a, status)) {
32 float_raise(float_flag_invalid, status);
33 }
34
35 if (status->default_nan_mode) {
36 return floatx80_default_nan(status);
37 }
38
39 return floatx80_maybe_silence_nan(a, status);
40 }
41
42 /*----------------------------------------------------------------------------
43 | Returns the modulo remainder of the extended double-precision floating-point
44 | value `a' with respect to the corresponding value `b'.
45 *----------------------------------------------------------------------------*/
46
47 floatx80 floatx80_mod(floatx80 a, floatx80 b, float_status *status)
48 {
49 flag aSign, zSign;
50 int32_t aExp, bExp, expDiff;
51 uint64_t aSig0, aSig1, bSig;
52 uint64_t qTemp, term0, term1;
53
54 aSig0 = extractFloatx80Frac(a);
55 aExp = extractFloatx80Exp(a);
56 aSign = extractFloatx80Sign(a);
57 bSig = extractFloatx80Frac(b);
58 bExp = extractFloatx80Exp(b);
59
60 if (aExp == 0x7FFF) {
61 if ((uint64_t) (aSig0 << 1)
62 || ((bExp == 0x7FFF) && (uint64_t) (bSig << 1))) {
63 return propagateFloatx80NaN(a, b, status);
64 }
65 goto invalid;
66 }
67 if (bExp == 0x7FFF) {
68 if ((uint64_t) (bSig << 1)) {
69 return propagateFloatx80NaN(a, b, status);
70 }
71 return a;
72 }
73 if (bExp == 0) {
74 if (bSig == 0) {
75 invalid:
76 float_raise(float_flag_invalid, status);
77 return floatx80_default_nan(status);
78 }
79 normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
80 }
81 if (aExp == 0) {
82 if ((uint64_t) (aSig0 << 1) == 0) {
83 return a;
84 }
85 normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
86 }
87 bSig |= LIT64(0x8000000000000000);
88 zSign = aSign;
89 expDiff = aExp - bExp;
90 aSig1 = 0;
91 if (expDiff < 0) {
92 return a;
93 }
94 qTemp = (bSig <= aSig0);
95 if (qTemp) {
96 aSig0 -= bSig;
97 }
98 expDiff -= 64;
99 while (0 < expDiff) {
100 qTemp = estimateDiv128To64(aSig0, aSig1, bSig);
101 qTemp = (2 < qTemp) ? qTemp - 2 : 0;
102 mul64To128(bSig, qTemp, &term0, &term1);
103 sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1);
104 shortShift128Left(aSig0, aSig1, 62, &aSig0, &aSig1);
105 }
106 expDiff += 64;
107 if (0 < expDiff) {
108 qTemp = estimateDiv128To64(aSig0, aSig1, bSig);
109 qTemp = (2 < qTemp) ? qTemp - 2 : 0;
110 qTemp >>= 64 - expDiff;
111 mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1);
112 sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1);
113 shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1);
114 while (le128(term0, term1, aSig0, aSig1)) {
115 ++qTemp;
116 sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1);
117 }
118 }
119 return
120 normalizeRoundAndPackFloatx80(
121 80, zSign, bExp + expDiff, aSig0, aSig1, status);
122 }
123
124 /*----------------------------------------------------------------------------
125 | Returns the mantissa of the extended double-precision floating-point
126 | value `a'.
127 *----------------------------------------------------------------------------*/
128
129 floatx80 floatx80_getman(floatx80 a, float_status *status)
130 {
131 flag aSign;
132 int32_t aExp;
133 uint64_t aSig;
134
135 aSig = extractFloatx80Frac(a);
136 aExp = extractFloatx80Exp(a);
137 aSign = extractFloatx80Sign(a);
138
139 if (aExp == 0x7FFF) {
140 if ((uint64_t) (aSig << 1)) {
141 return propagateFloatx80NaNOneArg(a , status);
142 }
143 float_raise(float_flag_invalid , status);
144 return floatx80_default_nan(status);
145 }
146
147 if (aExp == 0) {
148 if (aSig == 0) {
149 return packFloatx80(aSign, 0, 0);
150 }
151 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
152 }
153
154 return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign,
155 0x3FFF, aSig, 0, status);
156 }
157
158 /*----------------------------------------------------------------------------
159 | Returns the exponent of the extended double-precision floating-point
160 | value `a' as an extended double-precision value.
161 *----------------------------------------------------------------------------*/
162
163 floatx80 floatx80_getexp(floatx80 a, float_status *status)
164 {
165 flag aSign;
166 int32_t aExp;
167 uint64_t aSig;
168
169 aSig = extractFloatx80Frac(a);
170 aExp = extractFloatx80Exp(a);
171 aSign = extractFloatx80Sign(a);
172
173 if (aExp == 0x7FFF) {
174 if ((uint64_t) (aSig << 1)) {
175 return propagateFloatx80NaNOneArg(a , status);
176 }
177 float_raise(float_flag_invalid , status);
178 return floatx80_default_nan(status);
179 }
180
181 if (aExp == 0) {
182 if (aSig == 0) {
183 return packFloatx80(aSign, 0, 0);
184 }
185 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
186 }
187
188 return int32_to_floatx80(aExp - 0x3FFF, status);
189 }
190
191 /*----------------------------------------------------------------------------
192 | Scales extended double-precision floating-point value in operand `a' by
193 | value `b'. The function truncates the value in the second operand 'b' to
194 | an integral value and adds that value to the exponent of the operand 'a'.
195 | The operation performed according to the IEC/IEEE Standard for Binary
196 | Floating-Point Arithmetic.
197 *----------------------------------------------------------------------------*/
198
199 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status)
200 {
201 flag aSign, bSign;
202 int32_t aExp, bExp, shiftCount;
203 uint64_t aSig, bSig;
204
205 aSig = extractFloatx80Frac(a);
206 aExp = extractFloatx80Exp(a);
207 aSign = extractFloatx80Sign(a);
208 bSig = extractFloatx80Frac(b);
209 bExp = extractFloatx80Exp(b);
210 bSign = extractFloatx80Sign(b);
211
212 if (bExp == 0x7FFF) {
213 if ((uint64_t) (bSig << 1) ||
214 ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) {
215 return propagateFloatx80NaN(a, b, status);
216 }
217 float_raise(float_flag_invalid , status);
218 return floatx80_default_nan(status);
219 }
220 if (aExp == 0x7FFF) {
221 if ((uint64_t) (aSig << 1)) {
222 return propagateFloatx80NaN(a, b, status);
223 }
224 return packFloatx80(aSign, floatx80_infinity.high,
225 floatx80_infinity.low);
226 }
227 if (aExp == 0) {
228 if (aSig == 0) {
229 return packFloatx80(aSign, 0, 0);
230 }
231 if (bExp < 0x3FFF) {
232 return a;
233 }
234 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
235 }
236
237 if (bExp < 0x3FFF) {
238 return a;
239 }
240
241 if (0x400F < bExp) {
242 aExp = bSign ? -0x6001 : 0xE000;
243 return roundAndPackFloatx80(status->floatx80_rounding_precision,
244 aSign, aExp, aSig, 0, status);
245 }
246
247 shiftCount = 0x403E - bExp;
248 bSig >>= shiftCount;
249 aExp = bSign ? (aExp - bSig) : (aExp + bSig);
250
251 return roundAndPackFloatx80(status->floatx80_rounding_precision,
252 aSign, aExp, aSig, 0, status);
253 }
254
255 floatx80 floatx80_move(floatx80 a, float_status *status)
256 {
257 flag aSign;
258 int32_t aExp;
259 uint64_t aSig;
260
261 aSig = extractFloatx80Frac(a);
262 aExp = extractFloatx80Exp(a);
263 aSign = extractFloatx80Sign(a);
264
265 if (aExp == 0x7FFF) {
266 if ((uint64_t)(aSig << 1)) {
267 return propagateFloatx80NaNOneArg(a, status);
268 }
269 return a;
270 }
271 if (aExp == 0) {
272 if (aSig == 0) {
273 return a;
274 }
275 normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision,
276 aSign, aExp, aSig, 0, status);
277 }
278 return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign,
279 aExp, aSig, 0, status);
280 }
281
282 /*----------------------------------------------------------------------------
283 | Algorithms for transcendental functions supported by MC68881 and MC68882
284 | mathematical coprocessors. The functions are derived from FPSP library.
285 *----------------------------------------------------------------------------*/
286
287 #define one_exp 0x3FFF
288 #define one_sig LIT64(0x8000000000000000)
289
290 /*----------------------------------------------------------------------------
291 | Function for compactifying extended double-precision floating point values.
292 *----------------------------------------------------------------------------*/
293
294 static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig)
295 {
296 return (aExp << 16) | (aSig >> 48);
297 }
298
299 /*----------------------------------------------------------------------------
300 | Log base e of x plus 1
301 *----------------------------------------------------------------------------*/
302
303 floatx80 floatx80_lognp1(floatx80 a, float_status *status)
304 {
305 flag aSign;
306 int32_t aExp;
307 uint64_t aSig, fSig;
308
309 int8_t user_rnd_mode, user_rnd_prec;
310
311 int32_t compact, j, k;
312 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
313
314 aSig = extractFloatx80Frac(a);
315 aExp = extractFloatx80Exp(a);
316 aSign = extractFloatx80Sign(a);
317
318 if (aExp == 0x7FFF) {
319 if ((uint64_t) (aSig << 1)) {
320 propagateFloatx80NaNOneArg(a, status);
321 }
322 if (aSign) {
323 float_raise(float_flag_invalid, status);
324 return floatx80_default_nan(status);
325 }
326 return packFloatx80(0, floatx80_infinity.high, floatx80_infinity.low);
327 }
328
329 if (aExp == 0 && aSig == 0) {
330 return packFloatx80(aSign, 0, 0);
331 }
332
333 if (aSign && aExp >= one_exp) {
334 if (aExp == one_exp && aSig == one_sig) {
335 float_raise(float_flag_divbyzero, status);
336 packFloatx80(aSign, floatx80_infinity.high, floatx80_infinity.low);
337 }
338 float_raise(float_flag_invalid, status);
339 return floatx80_default_nan(status);
340 }
341
342 if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) {
343 /* <= min threshold */
344 float_raise(float_flag_inexact, status);
345 return floatx80_move(a, status);
346 }
347
348 user_rnd_mode = status->float_rounding_mode;
349 user_rnd_prec = status->floatx80_rounding_precision;
350 status->float_rounding_mode = float_round_nearest_even;
351 status->floatx80_rounding_precision = 80;
352
353 compact = floatx80_make_compact(aExp, aSig);
354
355 fp0 = a; /* Z */
356 fp1 = a;
357
358 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
359 status), status); /* X = (1+Z) */
360
361 aExp = extractFloatx80Exp(fp0);
362 aSig = extractFloatx80Frac(fp0);
363
364 compact = floatx80_make_compact(aExp, aSig);
365
366 if (compact < 0x3FFE8000 || compact > 0x3FFFC000) {
367 /* |X| < 1/2 or |X| > 3/2 */
368 k = aExp - 0x3FFF;
369 fp1 = int32_to_floatx80(k, status);
370
371 fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000);
372 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
373
374 f = packFloatx80(0, 0x3FFF, fSig); /* F */
375 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
376
377 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
378
379 lp1cont1:
380 /* LP1CONT1 */
381 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
382 logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC));
383 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
384 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
385
386 fp3 = fp2;
387 fp1 = fp2;
388
389 fp1 = floatx80_mul(fp1, float64_to_floatx80(
390 make_float64(0x3FC2499AB5E4040B), status),
391 status); /* V*A6 */
392 fp2 = floatx80_mul(fp2, float64_to_floatx80(
393 make_float64(0xBFC555B5848CB7DB), status),
394 status); /* V*A5 */
395 fp1 = floatx80_add(fp1, float64_to_floatx80(
396 make_float64(0x3FC99999987D8730), status),
397 status); /* A4+V*A6 */
398 fp2 = floatx80_add(fp2, float64_to_floatx80(
399 make_float64(0xBFCFFFFFFF6F7E97), status),
400 status); /* A3+V*A5 */
401 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
402 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
403 fp1 = floatx80_add(fp1, float64_to_floatx80(
404 make_float64(0x3FD55555555555A4), status),
405 status); /* A2+V*(A4+V*A6) */
406 fp2 = floatx80_add(fp2, float64_to_floatx80(
407 make_float64(0xBFE0000000000008), status),
408 status); /* A1+V*(A3+V*A5) */
409 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
410 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
411 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
412 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
413
414 fp1 = floatx80_add(fp1, log_tbl[j + 1],
415 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
416 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
417
418 status->float_rounding_mode = user_rnd_mode;
419 status->floatx80_rounding_precision = user_rnd_prec;
420
421 a = floatx80_add(fp0, klog2, status);
422
423 float_raise(float_flag_inexact, status);
424
425 return a;
426 } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
427 /* |X| < 1/16 or |X| > -1/16 */
428 /* LP1CARE */
429 fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000);
430 f = packFloatx80(0, 0x3FFF, fSig); /* F */
431 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
432
433 if (compact >= 0x3FFF8000) { /* 1+Z >= 1 */
434 /* KISZERO */
435 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x3F800000),
436 status), f, status); /* 1-F */
437 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (1-F)+Z */
438 fp1 = packFloatx80(0, 0, 0); /* K = 0 */
439 } else {
440 /* KISNEG */
441 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x40000000),
442 status), f, status); /* 2-F */
443 fp1 = floatx80_add(fp1, fp1, status); /* 2Z */
444 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (2-F)+2Z */
445 fp1 = packFloatx80(1, one_exp, one_sig); /* K = -1 */
446 }
447 goto lp1cont1;
448 } else {
449 /* LP1ONE16 */
450 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2Z */
451 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
452 status), status); /* FP0 IS 1+X */
453
454 /* LP1CONT2 */
455 fp1 = floatx80_div(fp1, fp0, status); /* U */
456 saveu = fp1;
457 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
458 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
459
460 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
461 status); /* B5 */
462 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
463 status); /* B4 */
464 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
465 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
466 fp3 = floatx80_add(fp3, float64_to_floatx80(
467 make_float64(0x3F624924928BCCFF), status),
468 status); /* B3+W*B5 */
469 fp2 = floatx80_add(fp2, float64_to_floatx80(
470 make_float64(0x3F899999999995EC), status),
471 status); /* B2+W*B4 */
472 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
473 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
474 fp1 = floatx80_add(fp1, float64_to_floatx80(
475 make_float64(0x3FB5555555555555), status),
476 status); /* B1+W*(B3+W*B5) */
477
478 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
479 fp1 = floatx80_add(fp1, fp2,
480 status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
481 fp0 = floatx80_mul(fp0, fp1,
482 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
483
484 status->float_rounding_mode = user_rnd_mode;
485 status->floatx80_rounding_precision = user_rnd_prec;
486
487 a = floatx80_add(fp0, saveu, status);
488
489 /*if (!floatx80_is_zero(a)) { */
490 float_raise(float_flag_inexact, status);
491 /*} */
492
493 return a;
494 }
495 }
496
497 /*----------------------------------------------------------------------------
498 | Log base e
499 *----------------------------------------------------------------------------*/
500
501 floatx80 floatx80_logn(floatx80 a, float_status *status)
502 {
503 flag aSign;
504 int32_t aExp;
505 uint64_t aSig, fSig;
506
507 int8_t user_rnd_mode, user_rnd_prec;
508
509 int32_t compact, j, k, adjk;
510 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
511
512 aSig = extractFloatx80Frac(a);
513 aExp = extractFloatx80Exp(a);
514 aSign = extractFloatx80Sign(a);
515
516 if (aExp == 0x7FFF) {
517 if ((uint64_t) (aSig << 1)) {
518 propagateFloatx80NaNOneArg(a, status);
519 }
520 if (aSign == 0) {
521 return packFloatx80(0, floatx80_infinity.high,
522 floatx80_infinity.low);
523 }
524 }
525
526 adjk = 0;
527
528 if (aExp == 0) {
529 if (aSig == 0) { /* zero */
530 float_raise(float_flag_divbyzero, status);
531 return packFloatx80(1, floatx80_infinity.high,
532 floatx80_infinity.low);
533 }
534 if ((aSig & one_sig) == 0) { /* denormal */
535 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
536 adjk = -100;
537 aExp += 100;
538 a = packFloatx80(aSign, aExp, aSig);
539 }
540 }
541
542 if (aSign) {
543 float_raise(float_flag_invalid, status);
544 return floatx80_default_nan(status);
545 }
546
547 user_rnd_mode = status->float_rounding_mode;
548 user_rnd_prec = status->floatx80_rounding_precision;
549 status->float_rounding_mode = float_round_nearest_even;
550 status->floatx80_rounding_precision = 80;
551
552 compact = floatx80_make_compact(aExp, aSig);
553
554 if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
555 /* |X| < 15/16 or |X| > 17/16 */
556 k = aExp - 0x3FFF;
557 k += adjk;
558 fp1 = int32_to_floatx80(k, status);
559
560 fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000);
561 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
562
563 f = packFloatx80(0, 0x3FFF, fSig); /* F */
564 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
565
566 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
567
568 /* LP1CONT1 */
569 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
570 logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC));
571 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
572 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
573
574 fp3 = fp2;
575 fp1 = fp2;
576
577 fp1 = floatx80_mul(fp1, float64_to_floatx80(
578 make_float64(0x3FC2499AB5E4040B), status),
579 status); /* V*A6 */
580 fp2 = floatx80_mul(fp2, float64_to_floatx80(
581 make_float64(0xBFC555B5848CB7DB), status),
582 status); /* V*A5 */
583 fp1 = floatx80_add(fp1, float64_to_floatx80(
584 make_float64(0x3FC99999987D8730), status),
585 status); /* A4+V*A6 */
586 fp2 = floatx80_add(fp2, float64_to_floatx80(
587 make_float64(0xBFCFFFFFFF6F7E97), status),
588 status); /* A3+V*A5 */
589 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
590 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
591 fp1 = floatx80_add(fp1, float64_to_floatx80(
592 make_float64(0x3FD55555555555A4), status),
593 status); /* A2+V*(A4+V*A6) */
594 fp2 = floatx80_add(fp2, float64_to_floatx80(
595 make_float64(0xBFE0000000000008), status),
596 status); /* A1+V*(A3+V*A5) */
597 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
598 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
599 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
600 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
601
602 fp1 = floatx80_add(fp1, log_tbl[j + 1],
603 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
604 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
605
606 status->float_rounding_mode = user_rnd_mode;
607 status->floatx80_rounding_precision = user_rnd_prec;
608
609 a = floatx80_add(fp0, klog2, status);
610
611 float_raise(float_flag_inexact, status);
612
613 return a;
614 } else { /* |X-1| >= 1/16 */
615 fp0 = a;
616 fp1 = a;
617 fp1 = floatx80_sub(fp1, float32_to_floatx80(make_float32(0x3F800000),
618 status), status); /* FP1 IS X-1 */
619 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
620 status), status); /* FP0 IS X+1 */
621 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2(X-1) */
622
623 /* LP1CONT2 */
624 fp1 = floatx80_div(fp1, fp0, status); /* U */
625 saveu = fp1;
626 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
627 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
628
629 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
630 status); /* B5 */
631 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
632 status); /* B4 */
633 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
634 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
635 fp3 = floatx80_add(fp3, float64_to_floatx80(
636 make_float64(0x3F624924928BCCFF), status),
637 status); /* B3+W*B5 */
638 fp2 = floatx80_add(fp2, float64_to_floatx80(
639 make_float64(0x3F899999999995EC), status),
640 status); /* B2+W*B4 */
641 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
642 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
643 fp1 = floatx80_add(fp1, float64_to_floatx80(
644 make_float64(0x3FB5555555555555), status),
645 status); /* B1+W*(B3+W*B5) */
646
647 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
648 fp1 = floatx80_add(fp1, fp2, status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
649 fp0 = floatx80_mul(fp0, fp1,
650 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
651
652 status->float_rounding_mode = user_rnd_mode;
653 status->floatx80_rounding_precision = user_rnd_prec;
654
655 a = floatx80_add(fp0, saveu, status);
656
657 /*if (!floatx80_is_zero(a)) { */
658 float_raise(float_flag_inexact, status);
659 /*} */
660
661 return a;
662 }
663 }
664
665 /*----------------------------------------------------------------------------
666 | Log base 10
667 *----------------------------------------------------------------------------*/
668
669 floatx80 floatx80_log10(floatx80 a, float_status *status)
670 {
671 flag aSign;
672 int32_t aExp;
673 uint64_t aSig;
674
675 int8_t user_rnd_mode, user_rnd_prec;
676
677 floatx80 fp0, fp1;
678
679 aSig = extractFloatx80Frac(a);
680 aExp = extractFloatx80Exp(a);
681 aSign = extractFloatx80Sign(a);
682
683 if (aExp == 0x7FFF) {
684 if ((uint64_t) (aSig << 1)) {
685 propagateFloatx80NaNOneArg(a, status);
686 }
687 if (aSign == 0) {
688 return packFloatx80(0, floatx80_infinity.high,
689 floatx80_infinity.low);
690 }
691 }
692
693 if (aExp == 0 && aSig == 0) {
694 float_raise(float_flag_divbyzero, status);
695 return packFloatx80(1, floatx80_infinity.high,
696 floatx80_infinity.low);
697 }
698
699 if (aSign) {
700 float_raise(float_flag_invalid, status);
701 return floatx80_default_nan(status);
702 }
703
704 user_rnd_mode = status->float_rounding_mode;
705 user_rnd_prec = status->floatx80_rounding_precision;
706 status->float_rounding_mode = float_round_nearest_even;
707 status->floatx80_rounding_precision = 80;
708
709 fp0 = floatx80_logn(a, status);
710 fp1 = packFloatx80(0, 0x3FFD, LIT64(0xDE5BD8A937287195)); /* INV_L10 */
711
712 status->float_rounding_mode = user_rnd_mode;
713 status->floatx80_rounding_precision = user_rnd_prec;
714
715 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L10 */
716
717 float_raise(float_flag_inexact, status);
718
719 return a;
720 }
721
722 /*----------------------------------------------------------------------------
723 | Log base 2
724 *----------------------------------------------------------------------------*/
725
726 floatx80 floatx80_log2(floatx80 a, float_status *status)
727 {
728 flag aSign;
729 int32_t aExp;
730 uint64_t aSig;
731
732 int8_t user_rnd_mode, user_rnd_prec;
733
734 floatx80 fp0, fp1;
735
736 aSig = extractFloatx80Frac(a);
737 aExp = extractFloatx80Exp(a);
738 aSign = extractFloatx80Sign(a);
739
740 if (aExp == 0x7FFF) {
741 if ((uint64_t) (aSig << 1)) {
742 propagateFloatx80NaNOneArg(a, status);
743 }
744 if (aSign == 0) {
745 return packFloatx80(0, floatx80_infinity.high,
746 floatx80_infinity.low);
747 }
748 }
749
750 if (aExp == 0) {
751 if (aSig == 0) {
752 float_raise(float_flag_divbyzero, status);
753 return packFloatx80(1, floatx80_infinity.high,
754 floatx80_infinity.low);
755 }
756 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
757 }
758
759 if (aSign) {
760 float_raise(float_flag_invalid, status);
761 return floatx80_default_nan(status);
762 }
763
764 user_rnd_mode = status->float_rounding_mode;
765 user_rnd_prec = status->floatx80_rounding_precision;
766 status->float_rounding_mode = float_round_nearest_even;
767 status->floatx80_rounding_precision = 80;
768
769 if (aSig == one_sig) { /* X is 2^k */
770 status->float_rounding_mode = user_rnd_mode;
771 status->floatx80_rounding_precision = user_rnd_prec;
772
773 a = int32_to_floatx80(aExp - 0x3FFF, status);
774 } else {
775 fp0 = floatx80_logn(a, status);
776 fp1 = packFloatx80(0, 0x3FFF, LIT64(0xB8AA3B295C17F0BC)); /* INV_L2 */
777
778 status->float_rounding_mode = user_rnd_mode;
779 status->floatx80_rounding_precision = user_rnd_prec;
780
781 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L2 */
782 }
783
784 float_raise(float_flag_inexact, status);
785
786 return a;
787 }
788
789 /*----------------------------------------------------------------------------
790 | e to x
791 *----------------------------------------------------------------------------*/
792
793 floatx80 floatx80_etox(floatx80 a, float_status *status)
794 {
795 flag aSign;
796 int32_t aExp;
797 uint64_t aSig;
798
799 int8_t user_rnd_mode, user_rnd_prec;
800
801 int32_t compact, n, j, k, m, m1;
802 floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale;
803 flag adjflag;
804
805 aSig = extractFloatx80Frac(a);
806 aExp = extractFloatx80Exp(a);
807 aSign = extractFloatx80Sign(a);
808
809 if (aExp == 0x7FFF) {
810 if ((uint64_t) (aSig << 1)) {
811 return propagateFloatx80NaNOneArg(a, status);
812 }
813 if (aSign) {
814 return packFloatx80(0, 0, 0);
815 }
816 return packFloatx80(0, floatx80_infinity.high,
817 floatx80_infinity.low);
818 }
819
820 if (aExp == 0 && aSig == 0) {
821 return packFloatx80(0, one_exp, one_sig);
822 }
823
824 user_rnd_mode = status->float_rounding_mode;
825 user_rnd_prec = status->floatx80_rounding_precision;
826 status->float_rounding_mode = float_round_nearest_even;
827 status->floatx80_rounding_precision = 80;
828
829 adjflag = 0;
830
831 if (aExp >= 0x3FBE) { /* |X| >= 2^(-65) */
832 compact = floatx80_make_compact(aExp, aSig);
833
834 if (compact < 0x400CB167) { /* |X| < 16380 log2 */
835 fp0 = a;
836 fp1 = a;
837 fp0 = floatx80_mul(fp0, float32_to_floatx80(
838 make_float32(0x42B8AA3B), status),
839 status); /* 64/log2 * X */
840 adjflag = 0;
841 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
842 fp0 = int32_to_floatx80(n, status);
843
844 j = n & 0x3F; /* J = N mod 64 */
845 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
846 if (n < 0 && j) {
847 /* arithmetic right shift is division and
848 * round towards minus infinity
849 */
850 m--;
851 }
852 m += 0x3FFF; /* biased exponent of 2^(M) */
853
854 expcont1:
855 fp2 = fp0; /* N */
856 fp0 = floatx80_mul(fp0, float32_to_floatx80(
857 make_float32(0xBC317218), status),
858 status); /* N * L1, L1 = lead(-log2/64) */
859 l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6));
860 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
861 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
862 fp0 = floatx80_add(fp0, fp2, status); /* R */
863
864 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
865 fp2 = float32_to_floatx80(make_float32(0x3AB60B70),
866 status); /* A5 */
867 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A5 */
868 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3C088895),
869 status), fp1,
870 status); /* fp3 is S*A4 */
871 fp2 = floatx80_add(fp2, float64_to_floatx80(make_float64(
872 0x3FA5555555554431), status),
873 status); /* fp2 is A3+S*A5 */
874 fp3 = floatx80_add(fp3, float64_to_floatx80(make_float64(
875 0x3FC5555555554018), status),
876 status); /* fp3 is A2+S*A4 */
877 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*(A3+S*A5) */
878 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 is S*(A2+S*A4) */
879 fp2 = floatx80_add(fp2, float32_to_floatx80(
880 make_float32(0x3F000000), status),
881 status); /* fp2 is A1+S*(A3+S*A5) */
882 fp3 = floatx80_mul(fp3, fp0, status); /* fp3 IS R*S*(A2+S*A4) */
883 fp2 = floatx80_mul(fp2, fp1,
884 status); /* fp2 IS S*(A1+S*(A3+S*A5)) */
885 fp0 = floatx80_add(fp0, fp3, status); /* fp0 IS R+R*S*(A2+S*A4) */
886 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
887
888 fp1 = exp_tbl[j];
889 fp0 = floatx80_mul(fp0, fp1, status); /* 2^(J/64)*(Exp(R)-1) */
890 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status),
891 status); /* accurate 2^(J/64) */
892 fp0 = floatx80_add(fp0, fp1,
893 status); /* 2^(J/64) + 2^(J/64)*(Exp(R)-1) */
894
895 scale = packFloatx80(0, m, one_sig);
896 if (adjflag) {
897 adjscale = packFloatx80(0, m1, one_sig);
898 fp0 = floatx80_mul(fp0, adjscale, status);
899 }
900
901 status->float_rounding_mode = user_rnd_mode;
902 status->floatx80_rounding_precision = user_rnd_prec;
903
904 a = floatx80_mul(fp0, scale, status);
905
906 float_raise(float_flag_inexact, status);
907
908 return a;
909 } else { /* |X| >= 16380 log2 */
910 if (compact > 0x400CB27C) { /* |X| >= 16480 log2 */
911 status->float_rounding_mode = user_rnd_mode;
912 status->floatx80_rounding_precision = user_rnd_prec;
913 if (aSign) {
914 a = roundAndPackFloatx80(
915 status->floatx80_rounding_precision,
916 0, -0x1000, aSig, 0, status);
917 } else {
918 a = roundAndPackFloatx80(
919 status->floatx80_rounding_precision,
920 0, 0x8000, aSig, 0, status);
921 }
922 float_raise(float_flag_inexact, status);
923
924 return a;
925 } else {
926 fp0 = a;
927 fp1 = a;
928 fp0 = floatx80_mul(fp0, float32_to_floatx80(
929 make_float32(0x42B8AA3B), status),
930 status); /* 64/log2 * X */
931 adjflag = 1;
932 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
933 fp0 = int32_to_floatx80(n, status);
934
935 j = n & 0x3F; /* J = N mod 64 */
936 /* NOTE: this is really arithmetic right shift by 6 */
937 k = n / 64;
938 if (n < 0 && j) {
939 /* arithmetic right shift is division and
940 * round towards minus infinity
941 */
942 k--;
943 }
944 /* NOTE: this is really arithmetic right shift by 1 */
945 m1 = k / 2;
946 if (k < 0 && (k & 1)) {
947 /* arithmetic right shift is division and
948 * round towards minus infinity
949 */
950 m1--;
951 }
952 m = k - m1;
953 m1 += 0x3FFF; /* biased exponent of 2^(M1) */
954 m += 0x3FFF; /* biased exponent of 2^(M) */
955
956 goto expcont1;
957 }
958 }
959 } else { /* |X| < 2^(-65) */
960 status->float_rounding_mode = user_rnd_mode;
961 status->floatx80_rounding_precision = user_rnd_prec;
962
963 a = floatx80_add(a, float32_to_floatx80(make_float32(0x3F800000),
964 status), status); /* 1 + X */
965
966 float_raise(float_flag_inexact, status);
967
968 return a;
969 }
970 }
971
972 /*----------------------------------------------------------------------------
973 | 2 to x
974 *----------------------------------------------------------------------------*/
975
976 floatx80 floatx80_twotox(floatx80 a, float_status *status)
977 {
978 flag aSign;
979 int32_t aExp;
980 uint64_t aSig;
981
982 int8_t user_rnd_mode, user_rnd_prec;
983
984 int32_t compact, n, j, l, m, m1;
985 floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
986
987 aSig = extractFloatx80Frac(a);
988 aExp = extractFloatx80Exp(a);
989 aSign = extractFloatx80Sign(a);
990
991 if (aExp == 0x7FFF) {
992 if ((uint64_t) (aSig << 1)) {
993 return propagateFloatx80NaNOneArg(a, status);
994 }
995 if (aSign) {
996 return packFloatx80(0, 0, 0);
997 }
998 return packFloatx80(0, floatx80_infinity.high,
999 floatx80_infinity.low);
1000 }
1001
1002 if (aExp == 0 && aSig == 0) {
1003 return packFloatx80(0, one_exp, one_sig);
1004 }
1005
1006 user_rnd_mode = status->float_rounding_mode;
1007 user_rnd_prec = status->floatx80_rounding_precision;
1008 status->float_rounding_mode = float_round_nearest_even;
1009 status->floatx80_rounding_precision = 80;
1010
1011 fp0 = a;
1012
1013 compact = floatx80_make_compact(aExp, aSig);
1014
1015 if (compact < 0x3FB98000 || compact > 0x400D80C0) {
1016 /* |X| > 16480 or |X| < 2^(-70) */
1017 if (compact > 0x3FFF8000) { /* |X| > 16480 */
1018 status->float_rounding_mode = user_rnd_mode;
1019 status->floatx80_rounding_precision = user_rnd_prec;
1020
1021 if (aSign) {
1022 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1023 0, -0x1000, aSig, 0, status);
1024 } else {
1025 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1026 0, 0x8000, aSig, 0, status);
1027 }
1028 } else { /* |X| < 2^(-70) */
1029 status->float_rounding_mode = user_rnd_mode;
1030 status->floatx80_rounding_precision = user_rnd_prec;
1031
1032 a = floatx80_add(fp0, float32_to_floatx80(
1033 make_float32(0x3F800000), status),
1034 status); /* 1 + X */
1035
1036 float_raise(float_flag_inexact, status);
1037
1038 return a;
1039 }
1040 } else { /* 2^(-70) <= |X| <= 16480 */
1041 fp1 = fp0; /* X */
1042 fp1 = floatx80_mul(fp1, float32_to_floatx80(
1043 make_float32(0x42800000), status),
1044 status); /* X * 64 */
1045 n = floatx80_to_int32(fp1, status);
1046 fp1 = int32_to_floatx80(n, status);
1047 j = n & 0x3F;
1048 l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
1049 if (n < 0 && j) {
1050 /* arithmetic right shift is division and
1051 * round towards minus infinity
1052 */
1053 l--;
1054 }
1055 m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */
1056 if (l < 0 && (l & 1)) {
1057 /* arithmetic right shift is division and
1058 * round towards minus infinity
1059 */
1060 m--;
1061 }
1062 m1 = l - m;
1063 m1 += 0x3FFF; /* ADJFACT IS 2^(M') */
1064
1065 adjfact = packFloatx80(0, m1, one_sig);
1066 fact1 = exp2_tbl[j];
1067 fact1.high += m;
1068 fact2.high = exp2_tbl2[j] >> 16;
1069 fact2.high += m;
1070 fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF);
1071 fact2.low <<= 48;
1072
1073 fp1 = floatx80_mul(fp1, float32_to_floatx80(
1074 make_float32(0x3C800000), status),
1075 status); /* (1/64)*N */
1076 fp0 = floatx80_sub(fp0, fp1, status); /* X - (1/64)*INT(64 X) */
1077 fp2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); /* LOG2 */
1078 fp0 = floatx80_mul(fp0, fp2, status); /* R */
1079
1080 /* EXPR */
1081 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1082 fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2),
1083 status); /* A5 */
1084 fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C),
1085 status); /* A4 */
1086 fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */
1087 fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */
1088 fp2 = floatx80_add(fp2, float64_to_floatx80(
1089 make_float64(0x3FA5555555554CC1), status),
1090 status); /* A3+S*A5 */
1091 fp3 = floatx80_add(fp3, float64_to_floatx80(
1092 make_float64(0x3FC5555555554A54), status),
1093 status); /* A2+S*A4 */
1094 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */
1095 fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */
1096 fp2 = floatx80_add(fp2, float64_to_floatx80(
1097 make_float64(0x3FE0000000000000), status),
1098 status); /* A1+S*(A3+S*A5) */
1099 fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */
1100
1101 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */
1102 fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */
1103 fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */
1104
1105 fp0 = floatx80_mul(fp0, fact1, status);
1106 fp0 = floatx80_add(fp0, fact2, status);
1107 fp0 = floatx80_add(fp0, fact1, status);
1108
1109 status->float_rounding_mode = user_rnd_mode;
1110 status->floatx80_rounding_precision = user_rnd_prec;
1111
1112 a = floatx80_mul(fp0, adjfact, status);
1113
1114 float_raise(float_flag_inexact, status);
1115
1116 return a;
1117 }
1118 }
1119
1120 /*----------------------------------------------------------------------------
1121 | 10 to x
1122 *----------------------------------------------------------------------------*/
1123
1124 floatx80 floatx80_tentox(floatx80 a, float_status *status)
1125 {
1126 flag aSign;
1127 int32_t aExp;
1128 uint64_t aSig;
1129
1130 int8_t user_rnd_mode, user_rnd_prec;
1131
1132 int32_t compact, n, j, l, m, m1;
1133 floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
1134
1135 aSig = extractFloatx80Frac(a);
1136 aExp = extractFloatx80Exp(a);
1137 aSign = extractFloatx80Sign(a);
1138
1139 if (aExp == 0x7FFF) {
1140 if ((uint64_t) (aSig << 1)) {
1141 return propagateFloatx80NaNOneArg(a, status);
1142 }
1143 if (aSign) {
1144 return packFloatx80(0, 0, 0);
1145 }
1146 return packFloatx80(0, floatx80_infinity.high,
1147 floatx80_infinity.low);
1148 }
1149
1150 if (aExp == 0 && aSig == 0) {
1151 return packFloatx80(0, one_exp, one_sig);
1152 }
1153
1154 user_rnd_mode = status->float_rounding_mode;
1155 user_rnd_prec = status->floatx80_rounding_precision;
1156 status->float_rounding_mode = float_round_nearest_even;
1157 status->floatx80_rounding_precision = 80;
1158
1159 fp0 = a;
1160
1161 compact = floatx80_make_compact(aExp, aSig);
1162
1163 if (compact < 0x3FB98000 || compact > 0x400B9B07) {
1164 /* |X| > 16480 LOG2/LOG10 or |X| < 2^(-70) */
1165 if (compact > 0x3FFF8000) { /* |X| > 16480 */
1166 status->float_rounding_mode = user_rnd_mode;
1167 status->floatx80_rounding_precision = user_rnd_prec;
1168
1169 if (aSign) {
1170 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1171 0, -0x1000, aSig, 0, status);
1172 } else {
1173 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1174 0, 0x8000, aSig, 0, status);
1175 }
1176 } else { /* |X| < 2^(-70) */
1177 status->float_rounding_mode = user_rnd_mode;
1178 status->floatx80_rounding_precision = user_rnd_prec;
1179
1180 a = floatx80_add(fp0, float32_to_floatx80(
1181 make_float32(0x3F800000), status),
1182 status); /* 1 + X */
1183
1184 float_raise(float_flag_inexact, status);
1185
1186 return a;
1187 }
1188 } else { /* 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10 */
1189 fp1 = fp0; /* X */
1190 fp1 = floatx80_mul(fp1, float64_to_floatx80(
1191 make_float64(0x406A934F0979A371),
1192 status), status); /* X*64*LOG10/LOG2 */
1193 n = floatx80_to_int32(fp1, status); /* N=INT(X*64*LOG10/LOG2) */
1194 fp1 = int32_to_floatx80(n, status);
1195
1196 j = n & 0x3F;
1197 l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
1198 if (n < 0 && j) {
1199 /* arithmetic right shift is division and
1200 * round towards minus infinity
1201 */
1202 l--;
1203 }
1204 m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */
1205 if (l < 0 && (l & 1)) {
1206 /* arithmetic right shift is division and
1207 * round towards minus infinity
1208 */
1209 m--;
1210 }
1211 m1 = l - m;
1212 m1 += 0x3FFF; /* ADJFACT IS 2^(M') */
1213
1214 adjfact = packFloatx80(0, m1, one_sig);
1215 fact1 = exp2_tbl[j];
1216 fact1.high += m;
1217 fact2.high = exp2_tbl2[j] >> 16;
1218 fact2.high += m;
1219 fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF);
1220 fact2.low <<= 48;
1221
1222 fp2 = fp1; /* N */
1223 fp1 = floatx80_mul(fp1, float64_to_floatx80(
1224 make_float64(0x3F734413509F8000), status),
1225 status); /* N*(LOG2/64LOG10)_LEAD */
1226 fp3 = packFloatx80(1, 0x3FCD, LIT64(0xC0219DC1DA994FD2));
1227 fp2 = floatx80_mul(fp2, fp3, status); /* N*(LOG2/64LOG10)_TRAIL */
1228 fp0 = floatx80_sub(fp0, fp1, status); /* X - N L_LEAD */
1229 fp0 = floatx80_sub(fp0, fp2, status); /* X - N L_TRAIL */
1230 fp2 = packFloatx80(0, 0x4000, LIT64(0x935D8DDDAAA8AC17)); /* LOG10 */
1231 fp0 = floatx80_mul(fp0, fp2, status); /* R */
1232
1233 /* EXPR */
1234 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1235 fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2),
1236 status); /* A5 */
1237 fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C),
1238 status); /* A4 */
1239 fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */
1240 fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */
1241 fp2 = floatx80_add(fp2, float64_to_floatx80(
1242 make_float64(0x3FA5555555554CC1), status),
1243 status); /* A3+S*A5 */
1244 fp3 = floatx80_add(fp3, float64_to_floatx80(
1245 make_float64(0x3FC5555555554A54), status),
1246 status); /* A2+S*A4 */
1247 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */
1248 fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */
1249 fp2 = floatx80_add(fp2, float64_to_floatx80(
1250 make_float64(0x3FE0000000000000), status),
1251 status); /* A1+S*(A3+S*A5) */
1252 fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */
1253
1254 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */
1255 fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */
1256 fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */
1257
1258 fp0 = floatx80_mul(fp0, fact1, status);
1259 fp0 = floatx80_add(fp0, fact2, status);
1260 fp0 = floatx80_add(fp0, fact1, status);
1261
1262 status->float_rounding_mode = user_rnd_mode;
1263 status->floatx80_rounding_precision = user_rnd_prec;
1264
1265 a = floatx80_mul(fp0, adjfact, status);
1266
1267 float_raise(float_flag_inexact, status);
1268
1269 return a;
1270 }
1271 }
1272
1273 /*----------------------------------------------------------------------------
1274 | Tangent
1275 *----------------------------------------------------------------------------*/
1276
1277 floatx80 floatx80_tan(floatx80 a, float_status *status)
1278 {
1279 flag aSign, xSign;
1280 int32_t aExp, xExp;
1281 uint64_t aSig, xSig;
1282
1283 int8_t user_rnd_mode, user_rnd_prec;
1284
1285 int32_t compact, l, n, j;
1286 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2;
1287 float32 twoto63;
1288 flag endflag;
1289
1290 aSig = extractFloatx80Frac(a);
1291 aExp = extractFloatx80Exp(a);
1292 aSign = extractFloatx80Sign(a);
1293
1294 if (aExp == 0x7FFF) {
1295 if ((uint64_t) (aSig << 1)) {
1296 return propagateFloatx80NaNOneArg(a, status);
1297 }
1298 float_raise(float_flag_invalid, status);
1299 return floatx80_default_nan(status);
1300 }
1301
1302 if (aExp == 0 && aSig == 0) {
1303 return packFloatx80(aSign, 0, 0);
1304 }
1305
1306 user_rnd_mode = status->float_rounding_mode;
1307 user_rnd_prec = status->floatx80_rounding_precision;
1308 status->float_rounding_mode = float_round_nearest_even;
1309 status->floatx80_rounding_precision = 80;
1310
1311 compact = floatx80_make_compact(aExp, aSig);
1312
1313 fp0 = a;
1314
1315 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1316 /* 2^(-40) > |X| > 15 PI */
1317 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1318 /* REDUCEX */
1319 fp1 = packFloatx80(0, 0, 0);
1320 if (compact == 0x7FFEFFFF) {
1321 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1322 LIT64(0xC90FDAA200000000));
1323 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1324 LIT64(0x85A308D300000000));
1325 fp0 = floatx80_add(fp0, twopi1, status);
1326 fp1 = fp0;
1327 fp0 = floatx80_add(fp0, twopi2, status);
1328 fp1 = floatx80_sub(fp1, fp0, status);
1329 fp1 = floatx80_add(fp1, twopi2, status);
1330 }
1331 loop:
1332 xSign = extractFloatx80Sign(fp0);
1333 xExp = extractFloatx80Exp(fp0);
1334 xExp -= 0x3FFF;
1335 if (xExp <= 28) {
1336 l = 0;
1337 endflag = 1;
1338 } else {
1339 l = xExp - 27;
1340 endflag = 0;
1341 }
1342 invtwopi = packFloatx80(0, 0x3FFE - l,
1343 LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */
1344 twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000));
1345 twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000));
1346
1347 /* SIGN(INARG)*2^63 IN SGL */
1348 twoto63 = packFloat32(xSign, 0xBE, 0);
1349
1350 fp2 = floatx80_mul(fp0, invtwopi, status);
1351 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1352 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1353 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1354 status); /* FP2 is N */
1355 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1356 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1357 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1358 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1359 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1360 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1361 fp3 = fp0; /* FP3 is A */
1362 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1363 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1364
1365 if (endflag > 0) {
1366 n = floatx80_to_int32(fp2, status);
1367 goto tancont;
1368 }
1369 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1370 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1371 goto loop;
1372 } else {
1373 status->float_rounding_mode = user_rnd_mode;
1374 status->floatx80_rounding_precision = user_rnd_prec;
1375
1376 a = floatx80_move(a, status);
1377
1378 float_raise(float_flag_inexact, status);
1379
1380 return a;
1381 }
1382 } else {
1383 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1384 make_float64(0x3FE45F306DC9C883), status),
1385 status); /* X*2/PI */
1386
1387 n = floatx80_to_int32(fp1, status);
1388 j = 32 + n;
1389
1390 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1391 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1392 status); /* FP0 IS R = (X-Y1)-Y2 */
1393
1394 tancont:
1395 if (n & 1) {
1396 /* NODD */
1397 fp1 = fp0; /* R */
1398 fp0 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1399 fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688),
1400 status); /* Q4 */
1401 fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04),
1402 status); /* P3 */
1403 fp3 = floatx80_mul(fp3, fp0, status); /* SQ4 */
1404 fp2 = floatx80_mul(fp2, fp0, status); /* SP3 */
1405 fp3 = floatx80_add(fp3, float64_to_floatx80(
1406 make_float64(0xBF346F59B39BA65F), status),
1407 status); /* Q3+SQ4 */
1408 fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00));
1409 fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */
1410 fp3 = floatx80_mul(fp3, fp0, status); /* S(Q3+SQ4) */
1411 fp2 = floatx80_mul(fp2, fp0, status); /* S(P2+SP3) */
1412 fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1));
1413 fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */
1414 fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA));
1415 fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */
1416 fp3 = floatx80_mul(fp3, fp0, status); /* S(Q2+S(Q3+SQ4)) */
1417 fp2 = floatx80_mul(fp2, fp0, status); /* S(P1+S(P2+SP3)) */
1418 fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE));
1419 fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */
1420 fp2 = floatx80_mul(fp2, fp1, status); /* RS(P1+S(P2+SP3)) */
1421 fp0 = floatx80_mul(fp0, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */
1422 fp1 = floatx80_add(fp1, fp2, status); /* R+RS(P1+S(P2+SP3)) */
1423 fp0 = floatx80_add(fp0, float32_to_floatx80(
1424 make_float32(0x3F800000), status),
1425 status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */
1426
1427 xSign = extractFloatx80Sign(fp1);
1428 xExp = extractFloatx80Exp(fp1);
1429 xSig = extractFloatx80Frac(fp1);
1430 xSign ^= 1;
1431 fp1 = packFloatx80(xSign, xExp, xSig);
1432
1433 status->float_rounding_mode = user_rnd_mode;
1434 status->floatx80_rounding_precision = user_rnd_prec;
1435
1436 a = floatx80_div(fp0, fp1, status);
1437
1438 float_raise(float_flag_inexact, status);
1439
1440 return a;
1441 } else {
1442 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1443 fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688),
1444 status); /* Q4 */
1445 fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04),
1446 status); /* P3 */
1447 fp3 = floatx80_mul(fp3, fp1, status); /* SQ4 */
1448 fp2 = floatx80_mul(fp2, fp1, status); /* SP3 */
1449 fp3 = floatx80_add(fp3, float64_to_floatx80(
1450 make_float64(0xBF346F59B39BA65F), status),
1451 status); /* Q3+SQ4 */
1452 fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00));
1453 fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */
1454 fp3 = floatx80_mul(fp3, fp1, status); /* S(Q3+SQ4) */
1455 fp2 = floatx80_mul(fp2, fp1, status); /* S(P2+SP3) */
1456 fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1));
1457 fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */
1458 fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA));
1459 fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */
1460 fp3 = floatx80_mul(fp3, fp1, status); /* S(Q2+S(Q3+SQ4)) */
1461 fp2 = floatx80_mul(fp2, fp1, status); /* S(P1+S(P2+SP3)) */
1462 fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE));
1463 fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */
1464 fp2 = floatx80_mul(fp2, fp0, status); /* RS(P1+S(P2+SP3)) */
1465 fp1 = floatx80_mul(fp1, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */
1466 fp0 = floatx80_add(fp0, fp2, status); /* R+RS(P1+S(P2+SP3)) */
1467 fp1 = floatx80_add(fp1, float32_to_floatx80(
1468 make_float32(0x3F800000), status),
1469 status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */
1470
1471 status->float_rounding_mode = user_rnd_mode;
1472 status->floatx80_rounding_precision = user_rnd_prec;
1473
1474 a = floatx80_div(fp0, fp1, status);
1475
1476 float_raise(float_flag_inexact, status);
1477
1478 return a;
1479 }
1480 }
1481 }
1482
1483 /*----------------------------------------------------------------------------
1484 | Sine
1485 *----------------------------------------------------------------------------*/
1486
1487 floatx80 floatx80_sin(floatx80 a, float_status *status)
1488 {
1489 flag aSign, xSign;
1490 int32_t aExp, xExp;
1491 uint64_t aSig, xSig;
1492
1493 int8_t user_rnd_mode, user_rnd_prec;
1494
1495 int32_t compact, l, n, j;
1496 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1497 float32 posneg1, twoto63;
1498 flag adjn, endflag;
1499
1500 aSig = extractFloatx80Frac(a);
1501 aExp = extractFloatx80Exp(a);
1502 aSign = extractFloatx80Sign(a);
1503
1504 if (aExp == 0x7FFF) {
1505 if ((uint64_t) (aSig << 1)) {
1506 return propagateFloatx80NaNOneArg(a, status);
1507 }
1508 float_raise(float_flag_invalid, status);
1509 return floatx80_default_nan(status);
1510 }
1511
1512 if (aExp == 0 && aSig == 0) {
1513 return packFloatx80(aSign, 0, 0);
1514 }
1515
1516 adjn = 0;
1517
1518 user_rnd_mode = status->float_rounding_mode;
1519 user_rnd_prec = status->floatx80_rounding_precision;
1520 status->float_rounding_mode = float_round_nearest_even;
1521 status->floatx80_rounding_precision = 80;
1522
1523 compact = floatx80_make_compact(aExp, aSig);
1524
1525 fp0 = a;
1526
1527 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1528 /* 2^(-40) > |X| > 15 PI */
1529 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1530 /* REDUCEX */
1531 fp1 = packFloatx80(0, 0, 0);
1532 if (compact == 0x7FFEFFFF) {
1533 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1534 LIT64(0xC90FDAA200000000));
1535 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1536 LIT64(0x85A308D300000000));
1537 fp0 = floatx80_add(fp0, twopi1, status);
1538 fp1 = fp0;
1539 fp0 = floatx80_add(fp0, twopi2, status);
1540 fp1 = floatx80_sub(fp1, fp0, status);
1541 fp1 = floatx80_add(fp1, twopi2, status);
1542 }
1543 loop:
1544 xSign = extractFloatx80Sign(fp0);
1545 xExp = extractFloatx80Exp(fp0);
1546 xExp -= 0x3FFF;
1547 if (xExp <= 28) {
1548 l = 0;
1549 endflag = 1;
1550 } else {
1551 l = xExp - 27;
1552 endflag = 0;
1553 }
1554 invtwopi = packFloatx80(0, 0x3FFE - l,
1555 LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */
1556 twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000));
1557 twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000));
1558
1559 /* SIGN(INARG)*2^63 IN SGL */
1560 twoto63 = packFloat32(xSign, 0xBE, 0);
1561
1562 fp2 = floatx80_mul(fp0, invtwopi, status);
1563 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1564 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1565 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1566 status); /* FP2 is N */
1567 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1568 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1569 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1570 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1571 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1572 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1573 fp3 = fp0; /* FP3 is A */
1574 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1575 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1576
1577 if (endflag > 0) {
1578 n = floatx80_to_int32(fp2, status);
1579 goto sincont;
1580 }
1581 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1582 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1583 goto loop;
1584 } else {
1585 /* SINSM */
1586 fp0 = float32_to_floatx80(make_float32(0x3F800000),
1587 status); /* 1 */
1588
1589 status->float_rounding_mode = user_rnd_mode;
1590 status->floatx80_rounding_precision = user_rnd_prec;
1591
1592 if (adjn) {
1593 /* COSTINY */
1594 a = floatx80_sub(fp0, float32_to_floatx80(
1595 make_float32(0x00800000), status), status);
1596 } else {
1597 /* SINTINY */
1598 a = floatx80_move(a, status);
1599 }
1600 float_raise(float_flag_inexact, status);
1601
1602 return a;
1603 }
1604 } else {
1605 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1606 make_float64(0x3FE45F306DC9C883), status),
1607 status); /* X*2/PI */
1608
1609 n = floatx80_to_int32(fp1, status);
1610 j = 32 + n;
1611
1612 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1613 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1614 status); /* FP0 IS R = (X-Y1)-Y2 */
1615
1616 sincont:
1617 if ((n + adjn) & 1) {
1618 /* COSPOLY */
1619 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1620 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1621 fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3),
1622 status); /* B8 */
1623 fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19),
1624 status); /* B7 */
1625
1626 xSign = extractFloatx80Sign(fp0); /* X IS S */
1627 xExp = extractFloatx80Exp(fp0);
1628 xSig = extractFloatx80Frac(fp0);
1629
1630 if (((n + adjn) >> 1) & 1) {
1631 xSign ^= 1;
1632 posneg1 = make_float32(0xBF800000); /* -1 */
1633 } else {
1634 xSign ^= 0;
1635 posneg1 = make_float32(0x3F800000); /* 1 */
1636 } /* X IS NOW R'= SGN*R */
1637
1638 fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */
1639 fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */
1640 fp2 = floatx80_add(fp2, float64_to_floatx80(
1641 make_float64(0x3E21EED90612C972), status),
1642 status); /* B6+TB8 */
1643 fp3 = floatx80_add(fp3, float64_to_floatx80(
1644 make_float64(0xBE927E4FB79D9FCF), status),
1645 status); /* B5+TB7 */
1646 fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */
1647 fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */
1648 fp2 = floatx80_add(fp2, float64_to_floatx80(
1649 make_float64(0x3EFA01A01A01D423), status),
1650 status); /* B4+T(B6+TB8) */
1651 fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438));
1652 fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */
1653 fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */
1654 fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */
1655 fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E));
1656 fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */
1657 fp1 = floatx80_add(fp1, float32_to_floatx80(
1658 make_float32(0xBF000000), status),
1659 status); /* B1+T(B3+T(B5+TB7)) */
1660 fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */
1661 fp0 = floatx80_add(fp0, fp1, status); /* [B1+T(B3+T(B5+TB7))]+
1662 * [S(B2+T(B4+T(B6+TB8)))]
1663 */
1664
1665 x = packFloatx80(xSign, xExp, xSig);
1666 fp0 = floatx80_mul(fp0, x, status);
1667
1668 status->float_rounding_mode = user_rnd_mode;
1669 status->floatx80_rounding_precision = user_rnd_prec;
1670
1671 a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status);
1672
1673 float_raise(float_flag_inexact, status);
1674
1675 return a;
1676 } else {
1677 /* SINPOLY */
1678 xSign = extractFloatx80Sign(fp0); /* X IS R */
1679 xExp = extractFloatx80Exp(fp0);
1680 xSig = extractFloatx80Frac(fp0);
1681
1682 xSign ^= ((n + adjn) >> 1) & 1; /* X IS NOW R'= SGN*R */
1683
1684 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1685 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1686 fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5),
1687 status); /* A7 */
1688 fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1),
1689 status); /* A6 */
1690 fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */
1691 fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */
1692 fp3 = floatx80_add(fp3, float64_to_floatx80(
1693 make_float64(0xBE5AE6452A118AE4), status),
1694 status); /* A5+T*A7 */
1695 fp2 = floatx80_add(fp2, float64_to_floatx80(
1696 make_float64(0x3EC71DE3A5341531), status),
1697 status); /* A4+T*A6 */
1698 fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */
1699 fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */
1700 fp3 = floatx80_add(fp3, float64_to_floatx80(
1701 make_float64(0xBF2A01A01A018B59), status),
1702 status); /* A3+T(A5+TA7) */
1703 fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF));
1704 fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */
1705 fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */
1706 fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */
1707 fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99));
1708 fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */
1709 fp1 = floatx80_add(fp1, fp2,
1710 status); /* [A1+T(A3+T(A5+TA7))]+
1711 * [S(A2+T(A4+TA6))]
1712 */
1713
1714 x = packFloatx80(xSign, xExp, xSig);
1715 fp0 = floatx80_mul(fp0, x, status); /* R'*S */
1716 fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */
1717
1718 status->float_rounding_mode = user_rnd_mode;
1719 status->floatx80_rounding_precision = user_rnd_prec;
1720
1721 a = floatx80_add(fp0, x, status);
1722
1723 float_raise(float_flag_inexact, status);
1724
1725 return a;
1726 }
1727 }
1728 }
1729
1730 /*----------------------------------------------------------------------------
1731 | Cosine
1732 *----------------------------------------------------------------------------*/
1733
1734 floatx80 floatx80_cos(floatx80 a, float_status *status)
1735 {
1736 flag aSign, xSign;
1737 int32_t aExp, xExp;
1738 uint64_t aSig, xSig;
1739
1740 int8_t user_rnd_mode, user_rnd_prec;
1741
1742 int32_t compact, l, n, j;
1743 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1744 float32 posneg1, twoto63;
1745 flag adjn, endflag;
1746
1747 aSig = extractFloatx80Frac(a);
1748 aExp = extractFloatx80Exp(a);
1749 aSign = extractFloatx80Sign(a);
1750
1751 if (aExp == 0x7FFF) {
1752 if ((uint64_t) (aSig << 1)) {
1753 return propagateFloatx80NaNOneArg(a, status);
1754 }
1755 float_raise(float_flag_invalid, status);
1756 return floatx80_default_nan(status);
1757 }
1758
1759 if (aExp == 0 && aSig == 0) {
1760 return packFloatx80(0, one_exp, one_sig);
1761 }
1762
1763 adjn = 1;
1764
1765 user_rnd_mode = status->float_rounding_mode;
1766 user_rnd_prec = status->floatx80_rounding_precision;
1767 status->float_rounding_mode = float_round_nearest_even;
1768 status->floatx80_rounding_precision = 80;
1769
1770 compact = floatx80_make_compact(aExp, aSig);
1771
1772 fp0 = a;
1773
1774 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1775 /* 2^(-40) > |X| > 15 PI */
1776 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1777 /* REDUCEX */
1778 fp1 = packFloatx80(0, 0, 0);
1779 if (compact == 0x7FFEFFFF) {
1780 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1781 LIT64(0xC90FDAA200000000));
1782 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1783 LIT64(0x85A308D300000000));
1784 fp0 = floatx80_add(fp0, twopi1, status);
1785 fp1 = fp0;
1786 fp0 = floatx80_add(fp0, twopi2, status);
1787 fp1 = floatx80_sub(fp1, fp0, status);
1788 fp1 = floatx80_add(fp1, twopi2, status);
1789 }
1790 loop:
1791 xSign = extractFloatx80Sign(fp0);
1792 xExp = extractFloatx80Exp(fp0);
1793 xExp -= 0x3FFF;
1794 if (xExp <= 28) {
1795 l = 0;
1796 endflag = 1;
1797 } else {
1798 l = xExp - 27;
1799 endflag = 0;
1800 }
1801 invtwopi = packFloatx80(0, 0x3FFE - l,
1802 LIT64(0xA2F9836E4E44152A)); /* INVTWOPI */
1803 twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000));
1804 twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000));
1805
1806 /* SIGN(INARG)*2^63 IN SGL */
1807 twoto63 = packFloat32(xSign, 0xBE, 0);
1808
1809 fp2 = floatx80_mul(fp0, invtwopi, status);
1810 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1811 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1812 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1813 status); /* FP2 is N */
1814 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1815 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1816 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1817 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1818 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1819 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1820 fp3 = fp0; /* FP3 is A */
1821 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1822 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1823
1824 if (endflag > 0) {
1825 n = floatx80_to_int32(fp2, status);
1826 goto sincont;
1827 }
1828 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1829 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1830 goto loop;
1831 } else {
1832 /* SINSM */
1833 fp0 = float32_to_floatx80(make_float32(0x3F800000), status); /* 1 */
1834
1835 status->float_rounding_mode = user_rnd_mode;
1836 status->floatx80_rounding_precision = user_rnd_prec;
1837
1838 if (adjn) {
1839 /* COSTINY */
1840 a = floatx80_sub(fp0, float32_to_floatx80(
1841 make_float32(0x00800000), status),
1842 status);
1843 } else {
1844 /* SINTINY */
1845 a = floatx80_move(a, status);
1846 }
1847 float_raise(float_flag_inexact, status);
1848
1849 return a;
1850 }
1851 } else {
1852 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1853 make_float64(0x3FE45F306DC9C883), status),
1854 status); /* X*2/PI */
1855
1856 n = floatx80_to_int32(fp1, status);
1857 j = 32 + n;
1858
1859 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1860 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1861 status); /* FP0 IS R = (X-Y1)-Y2 */
1862
1863 sincont:
1864 if ((n + adjn) & 1) {
1865 /* COSPOLY */
1866 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1867 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1868 fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3),
1869 status); /* B8 */
1870 fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19),
1871 status); /* B7 */
1872
1873 xSign = extractFloatx80Sign(fp0); /* X IS S */
1874 xExp = extractFloatx80Exp(fp0);
1875 xSig = extractFloatx80Frac(fp0);
1876
1877 if (((n + adjn) >> 1) & 1) {
1878 xSign ^= 1;
1879 posneg1 = make_float32(0xBF800000); /* -1 */
1880 } else {
1881 xSign ^= 0;
1882 posneg1 = make_float32(0x3F800000); /* 1 */
1883 } /* X IS NOW R'= SGN*R */
1884
1885 fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */
1886 fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */
1887 fp2 = floatx80_add(fp2, float64_to_floatx80(
1888 make_float64(0x3E21EED90612C972), status),
1889 status); /* B6+TB8 */
1890 fp3 = floatx80_add(fp3, float64_to_floatx80(
1891 make_float64(0xBE927E4FB79D9FCF), status),
1892 status); /* B5+TB7 */
1893 fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */
1894 fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */
1895 fp2 = floatx80_add(fp2, float64_to_floatx80(
1896 make_float64(0x3EFA01A01A01D423), status),
1897 status); /* B4+T(B6+TB8) */
1898 fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438));
1899 fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */
1900 fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */
1901 fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */
1902 fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E));
1903 fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */
1904 fp1 = floatx80_add(fp1, float32_to_floatx80(
1905 make_float32(0xBF000000), status),
1906 status); /* B1+T(B3+T(B5+TB7)) */
1907 fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */
1908 fp0 = floatx80_add(fp0, fp1, status);
1909 /* [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] */
1910
1911 x = packFloatx80(xSign, xExp, xSig);
1912 fp0 = floatx80_mul(fp0, x, status);
1913
1914 status->float_rounding_mode = user_rnd_mode;
1915 status->floatx80_rounding_precision = user_rnd_prec;
1916
1917 a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status);
1918
1919 float_raise(float_flag_inexact, status);
1920
1921 return a;
1922 } else {
1923 /* SINPOLY */
1924 xSign = extractFloatx80Sign(fp0); /* X IS R */
1925 xExp = extractFloatx80Exp(fp0);
1926 xSig = extractFloatx80Frac(fp0);
1927
1928 xSign ^= ((n + adjn) >> 1) & 1; /* X IS NOW R'= SGN*R */
1929
1930 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1931 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1932 fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5),
1933 status); /* A7 */
1934 fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1),
1935 status); /* A6 */
1936 fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */
1937 fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */
1938 fp3 = floatx80_add(fp3, float64_to_floatx80(
1939 make_float64(0xBE5AE6452A118AE4), status),
1940 status); /* A5+T*A7 */
1941 fp2 = floatx80_add(fp2, float64_to_floatx80(
1942 make_float64(0x3EC71DE3A5341531), status),
1943 status); /* A4+T*A6 */
1944 fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */
1945 fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */
1946 fp3 = floatx80_add(fp3, float64_to_floatx80(
1947 make_float64(0xBF2A01A01A018B59), status),
1948 status); /* A3+T(A5+TA7) */
1949 fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF));
1950 fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */
1951 fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */
1952 fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */
1953 fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99));
1954 fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */
1955 fp1 = floatx80_add(fp1, fp2, status);
1956 /* [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] */
1957
1958 x = packFloatx80(xSign, xExp, xSig);
1959 fp0 = floatx80_mul(fp0, x, status); /* R'*S */
1960 fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */
1961
1962 status->float_rounding_mode = user_rnd_mode;
1963 status->floatx80_rounding_precision = user_rnd_prec;
1964
1965 a = floatx80_add(fp0, x, status);
1966
1967 float_raise(float_flag_inexact, status);
1968
1969 return a;
1970 }
1971 }
1972 }
1973
1974 /*----------------------------------------------------------------------------
1975 | Arc tangent
1976 *----------------------------------------------------------------------------*/
1977
1978 floatx80 floatx80_atan(floatx80 a, float_status *status)
1979 {
1980 flag aSign;
1981 int32_t aExp;
1982 uint64_t aSig;
1983
1984 int8_t user_rnd_mode, user_rnd_prec;
1985
1986 int32_t compact, tbl_index;
1987 floatx80 fp0, fp1, fp2, fp3, xsave;
1988
1989 aSig = extractFloatx80Frac(a);
1990 aExp = extractFloatx80Exp(a);
1991 aSign = extractFloatx80Sign(a);
1992
1993 if (aExp == 0x7FFF) {
1994 if ((uint64_t) (aSig << 1)) {
1995 return propagateFloatx80NaNOneArg(a, status);
1996 }
1997 a = packFloatx80(aSign, piby2_exp, pi_sig);
1998 float_raise(float_flag_inexact, status);
1999 return floatx80_move(a, status);
2000 }
2001
2002 if (aExp == 0 && aSig == 0) {
2003 return packFloatx80(aSign, 0, 0);
2004 }
2005
2006 compact = floatx80_make_compact(aExp, aSig);
2007
2008 user_rnd_mode = status->float_rounding_mode;
2009 user_rnd_prec = status->floatx80_rounding_precision;
2010 status->float_rounding_mode = float_round_nearest_even;
2011 status->floatx80_rounding_precision = 80;
2012
2013 if (compact < 0x3FFB8000 || compact > 0x4002FFFF) {
2014 /* |X| >= 16 or |X| < 1/16 */
2015 if (compact > 0x3FFF8000) { /* |X| >= 16 */
2016 if (compact > 0x40638000) { /* |X| > 2^(100) */
2017 fp0 = packFloatx80(aSign, piby2_exp, pi_sig);
2018 fp1 = packFloatx80(aSign, 0x0001, one_sig);
2019
2020 status->float_rounding_mode = user_rnd_mode;
2021 status->floatx80_rounding_precision = user_rnd_prec;
2022
2023 a = floatx80_sub(fp0, fp1, status);
2024
2025 float_raise(float_flag_inexact, status);
2026
2027 return a;
2028 } else {
2029 fp0 = a;
2030 fp1 = packFloatx80(1, one_exp, one_sig); /* -1 */
2031 fp1 = floatx80_div(fp1, fp0, status); /* X' = -1/X */
2032 xsave = fp1;
2033 fp0 = floatx80_mul(fp1, fp1, status); /* Y = X'*X' */
2034 fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */
2035 fp3 = float64_to_floatx80(make_float64(0xBFB70BF398539E6A),
2036 status); /* C5 */
2037 fp2 = float64_to_floatx80(make_float64(0x3FBC7187962D1D7D),
2038 status); /* C4 */
2039 fp3 = floatx80_mul(fp3, fp1, status); /* Z*C5 */
2040 fp2 = floatx80_mul(fp2, fp1, status); /* Z*C4 */
2041 fp3 = floatx80_add(fp3, float64_to_floatx80(
2042 make_float64(0xBFC24924827107B8), status),
2043 status); /* C3+Z*C5 */
2044 fp2 = floatx80_add(fp2, float64_to_floatx80(
2045 make_float64(0x3FC999999996263E), status),
2046 status); /* C2+Z*C4 */
2047 fp1 = floatx80_mul(fp1, fp3, status); /* Z*(C3+Z*C5) */
2048 fp2 = floatx80_mul(fp2, fp0, status); /* Y*(C2+Z*C4) */
2049 fp1 = floatx80_add(fp1, float64_to_floatx80(
2050 make_float64(0xBFD5555555555536), status),
2051 status); /* C1+Z*(C3+Z*C5) */
2052 fp0 = floatx80_mul(fp0, xsave, status); /* X'*Y */
2053 /* [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] */
2054 fp1 = floatx80_add(fp1, fp2, status);
2055 /* X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ?? */
2056 fp0 = floatx80_mul(fp0, fp1, status);
2057 fp0 = floatx80_add(fp0, xsave, status);
2058 fp1 = packFloatx80(aSign, piby2_exp, pi_sig);
2059
2060 status->float_rounding_mode = user_rnd_mode;
2061 status->floatx80_rounding_precision = user_rnd_prec;
2062
2063 a = floatx80_add(fp0, fp1, status);
2064
2065 float_raise(float_flag_inexact, status);
2066
2067 return a;
2068 }
2069 } else { /* |X| < 1/16 */
2070 if (compact < 0x3FD78000) { /* |X| < 2^(-40) */
2071 status->float_rounding_mode = user_rnd_mode;
2072 status->floatx80_rounding_precision = user_rnd_prec;
2073
2074 a = floatx80_move(a, status);
2075
2076 float_raise(float_flag_inexact, status);
2077
2078 return a;
2079 } else {
2080 fp0 = a;
2081 xsave = a;
2082 fp0 = floatx80_mul(fp0, fp0, status); /* Y = X*X */
2083 fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */
2084 fp2 = float64_to_floatx80(make_float64(0x3FB344447F876989),
2085 status); /* B6 */
2086 fp3 = float64_to_floatx80(make_float64(0xBFB744EE7FAF45DB),
2087 status); /* B5 */
2088 fp2 = floatx80_mul(fp2, fp1, status); /* Z*B6 */
2089 fp3 = floatx80_mul(fp3, fp1, status); /* Z*B5 */
2090 fp2 = floatx80_add(fp2, float64_to_floatx80(
2091 make_float64(0x3FBC71C646940220), status),
2092 status); /* B4+Z*B6 */
2093 fp3 = floatx80_add(fp3, float64_to_floatx80(
2094 make_float64(0xBFC24924921872F9),
2095 status), status); /* B3+Z*B5 */
2096 fp2 = floatx80_mul(fp2, fp1, status); /* Z*(B4+Z*B6) */
2097 fp1 = floatx80_mul(fp1, fp3, status); /* Z*(B3+Z*B5) */
2098 fp2 = floatx80_add(fp2, float64_to_floatx80(
2099 make_float64(0x3FC9999999998FA9), status),
2100 status); /* B2+Z*(B4+Z*B6) */
2101 fp1 = floatx80_add(fp1, float64_to_floatx80(
2102 make_float64(0xBFD5555555555555), status),
2103 status); /* B1+Z*(B3+Z*B5) */
2104 fp2 = floatx80_mul(fp2, fp0, status); /* Y*(B2+Z*(B4+Z*B6)) */
2105 fp0 = floatx80_mul(fp0, xsave, status); /* X*Y */
2106 /* [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] */
2107 fp1 = floatx80_add(fp1, fp2, status);
2108 /* X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) */
2109 fp0 = floatx80_mul(fp0, fp1, status);
2110
2111 status->float_rounding_mode = user_rnd_mode;
2112 status->floatx80_rounding_precision = user_rnd_prec;
2113
2114 a = floatx80_add(fp0, xsave, status);
2115
2116 float_raise(float_flag_inexact, status);
2117
2118 return a;
2119 }
2120 }
2121 } else {
2122 aSig &= LIT64(0xF800000000000000);
2123 aSig |= LIT64(0x0400000000000000);
2124 xsave = packFloatx80(aSign, aExp, aSig); /* F */
2125 fp0 = a;
2126 fp1 = a; /* X */
2127 fp2 = packFloatx80(0, one_exp, one_sig); /* 1 */
2128 fp1 = floatx80_mul(fp1, xsave, status); /* X*F */
2129 fp0 = floatx80_sub(fp0, xsave, status); /* X-F */
2130 fp1 = floatx80_add(fp1, fp2, status); /* 1 + X*F */
2131 fp0 = floatx80_div(fp0, fp1, status); /* U = (X-F)/(1+X*F) */
2132
2133 tbl_index = compact;
2134
2135 tbl_index &= 0x7FFF0000;
2136 tbl_index -= 0x3FFB0000;
2137 tbl_index >>= 1;
2138 tbl_index += compact & 0x00007800;
2139 tbl_index >>= 11;
2140
2141 fp3 = atan_tbl[tbl_index];
2142
2143 fp3.high |= aSign ? 0x8000 : 0; /* ATAN(F) */
2144
2145 fp1 = floatx80_mul(fp0, fp0, status); /* V = U*U */
2146 fp2 = float64_to_floatx80(make_float64(0xBFF6687E314987D8),
2147 status); /* A3 */
2148 fp2 = floatx80_add(fp2, fp1, status); /* A3+V */
2149 fp2 = floatx80_mul(fp2, fp1, status); /* V*(A3+V) */
2150 fp1 = floatx80_mul(fp1, fp0, status); /* U*V */
2151 fp2 = floatx80_add(fp2, float64_to_floatx80(
2152 make_float64(0x4002AC6934A26DB3), status),
2153 status); /* A2+V*(A3+V) */
2154 fp1 = floatx80_mul(fp1, float64_to_floatx80(
2155 make_float64(0xBFC2476F4E1DA28E), status),
2156 status); /* A1+U*V */
2157 fp1 = floatx80_mul(fp1, fp2, status); /* A1*U*V*(A2+V*(A3+V)) */
2158 fp0 = floatx80_add(fp0, fp1, status); /* ATAN(U) */
2159
2160 status->float_rounding_mode = user_rnd_mode;
2161 status->floatx80_rounding_precision = user_rnd_prec;
2162
2163 a = floatx80_add(fp0, fp3, status); /* ATAN(X) */
2164
2165 float_raise(float_flag_inexact, status);
2166
2167 return a;
2168 }
2169 }
2170
2171 /*----------------------------------------------------------------------------
2172 | Arc sine
2173 *----------------------------------------------------------------------------*/
2174
2175 floatx80 floatx80_asin(floatx80 a, float_status *status)
2176 {
2177 flag aSign;
2178 int32_t aExp;
2179 uint64_t aSig;
2180
2181 int8_t user_rnd_mode, user_rnd_prec;
2182
2183 int32_t compact;
2184 floatx80 fp0, fp1, fp2, one;
2185
2186 aSig = extractFloatx80Frac(a);
2187 aExp = extractFloatx80Exp(a);
2188 aSign = extractFloatx80Sign(a);
2189
2190 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2191 return propagateFloatx80NaNOneArg(a, status);
2192 }
2193
2194 if (aExp == 0 && aSig == 0) {
2195 return packFloatx80(aSign, 0, 0);
2196 }
2197
2198 compact = floatx80_make_compact(aExp, aSig);
2199
2200 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2201 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2202 float_raise(float_flag_inexact, status);
2203 a = packFloatx80(aSign, piby2_exp, pi_sig);
2204 return floatx80_move(a, status);
2205 } else { /* |X| > 1 */
2206 float_raise(float_flag_invalid, status);
2207 return floatx80_default_nan(status);
2208 }
2209
2210 } /* |X| < 1 */
2211
2212 user_rnd_mode = status->float_rounding_mode;
2213 user_rnd_prec = status->floatx80_rounding_precision;
2214 status->float_rounding_mode = float_round_nearest_even;
2215 status->floatx80_rounding_precision = 80;
2216
2217 one = packFloatx80(0, one_exp, one_sig);
2218 fp0 = a;
2219
2220 fp1 = floatx80_sub(one, fp0, status); /* 1 - X */
2221 fp2 = floatx80_add(one, fp0, status); /* 1 + X */
2222 fp1 = floatx80_mul(fp2, fp1, status); /* (1+X)*(1-X) */
2223 fp1 = floatx80_sqrt(fp1, status); /* SQRT((1+X)*(1-X)) */
2224 fp0 = floatx80_div(fp0, fp1, status); /* X/SQRT((1+X)*(1-X)) */
2225
2226 status->float_rounding_mode = user_rnd_mode;
2227 status->floatx80_rounding_precision = user_rnd_prec;
2228
2229 a = floatx80_atan(fp0, status); /* ATAN(X/SQRT((1+X)*(1-X))) */
2230
2231 float_raise(float_flag_inexact, status);
2232
2233 return a;
2234 }
AR7 emulation for QEMU