?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
qapi/error: Check format string argument in error_*prepend()
[qemu/armbru.git] / target / m68k / softfloat.c
blobb6d0ed7acfe8dd122e052d2ea9348e9ef0f750fe
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 /*
18 * Portions of this work are licensed under the terms of the GNU GPL,
19 * version 2 or later. See the COPYING file in the top-level directory.
20 */
21
22 #include "qemu/osdep.h"
23 #include "softfloat.h"
24 #include "fpu/softfloat-macros.h"
25 #include "softfloat_fpsp_tables.h"
26
27 #define pi_exp 0x4000
28 #define piby2_exp 0x3FFF
29 #define pi_sig UINT64_C(0xc90fdaa22168c235)
30
31 static floatx80 propagateFloatx80NaNOneArg(floatx80 a, float_status *status)
32 {
33 if (floatx80_is_signaling_nan(a, status)) {
34 float_raise(float_flag_invalid, status);
35 a = floatx80_silence_nan(a, status);
36 }
37
38 if (status->default_nan_mode) {
39 return floatx80_default_nan(status);
40 }
41
42 return a;
43 }
44
45 /*
46 * Returns the mantissa of the extended double-precision floating-point
47 * value `a'.
48 */
49
50 floatx80 floatx80_getman(floatx80 a, float_status *status)
51 {
52 bool aSign;
53 int32_t aExp;
54 uint64_t aSig;
55
56 aSig = extractFloatx80Frac(a);
57 aExp = extractFloatx80Exp(a);
58 aSign = extractFloatx80Sign(a);
59
60 if (aExp == 0x7FFF) {
61 if ((uint64_t) (aSig << 1)) {
62 return propagateFloatx80NaNOneArg(a , status);
63 }
64 float_raise(float_flag_invalid , status);
65 return floatx80_default_nan(status);
66 }
67
68 if (aExp == 0) {
69 if (aSig == 0) {
70 return packFloatx80(aSign, 0, 0);
71 }
72 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
73 }
74
75 return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign,
76 0x3FFF, aSig, 0, status);
77 }
78
79 /*
80 * Returns the exponent of the extended double-precision floating-point
81 * value `a' as an extended double-precision value.
82 */
83
84 floatx80 floatx80_getexp(floatx80 a, float_status *status)
85 {
86 bool aSign;
87 int32_t aExp;
88 uint64_t aSig;
89
90 aSig = extractFloatx80Frac(a);
91 aExp = extractFloatx80Exp(a);
92 aSign = extractFloatx80Sign(a);
93
94 if (aExp == 0x7FFF) {
95 if ((uint64_t) (aSig << 1)) {
96 return propagateFloatx80NaNOneArg(a , status);
97 }
98 float_raise(float_flag_invalid , status);
99 return floatx80_default_nan(status);
100 }
101
102 if (aExp == 0) {
103 if (aSig == 0) {
104 return packFloatx80(aSign, 0, 0);
105 }
106 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
107 }
108
109 return int32_to_floatx80(aExp - 0x3FFF, status);
110 }
111
112 /*
113 * Scales extended double-precision floating-point value in operand `a' by
114 * value `b'. The function truncates the value in the second operand 'b' to
115 * an integral value and adds that value to the exponent of the operand 'a'.
116 * The operation performed according to the IEC/IEEE Standard for Binary
117 * Floating-Point Arithmetic.
118 */
119
120 floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status)
121 {
122 bool aSign, bSign;
123 int32_t aExp, bExp, shiftCount;
124 uint64_t aSig, bSig;
125
126 aSig = extractFloatx80Frac(a);
127 aExp = extractFloatx80Exp(a);
128 aSign = extractFloatx80Sign(a);
129 bSig = extractFloatx80Frac(b);
130 bExp = extractFloatx80Exp(b);
131 bSign = extractFloatx80Sign(b);
132
133 if (bExp == 0x7FFF) {
134 if ((uint64_t) (bSig << 1) ||
135 ((aExp == 0x7FFF) && (uint64_t) (aSig << 1))) {
136 return propagateFloatx80NaN(a, b, status);
137 }
138 float_raise(float_flag_invalid , status);
139 return floatx80_default_nan(status);
140 }
141 if (aExp == 0x7FFF) {
142 if ((uint64_t) (aSig << 1)) {
143 return propagateFloatx80NaN(a, b, status);
144 }
145 return packFloatx80(aSign, floatx80_infinity.high,
146 floatx80_infinity.low);
147 }
148 if (aExp == 0) {
149 if (aSig == 0) {
150 return packFloatx80(aSign, 0, 0);
151 }
152 if (bExp < 0x3FFF) {
153 return a;
154 }
155 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
156 }
157
158 if (bExp < 0x3FFF) {
159 return a;
160 }
161
162 if (0x400F < bExp) {
163 aExp = bSign ? -0x6001 : 0xE000;
164 return roundAndPackFloatx80(status->floatx80_rounding_precision,
165 aSign, aExp, aSig, 0, status);
166 }
167
168 shiftCount = 0x403E - bExp;
169 bSig >>= shiftCount;
170 aExp = bSign ? (aExp - bSig) : (aExp + bSig);
171
172 return roundAndPackFloatx80(status->floatx80_rounding_precision,
173 aSign, aExp, aSig, 0, status);
174 }
175
176 floatx80 floatx80_move(floatx80 a, float_status *status)
177 {
178 bool aSign;
179 int32_t aExp;
180 uint64_t aSig;
181
182 aSig = extractFloatx80Frac(a);
183 aExp = extractFloatx80Exp(a);
184 aSign = extractFloatx80Sign(a);
185
186 if (aExp == 0x7FFF) {
187 if ((uint64_t)(aSig << 1)) {
188 return propagateFloatx80NaNOneArg(a, status);
189 }
190 return a;
191 }
192 if (aExp == 0) {
193 if (aSig == 0) {
194 return a;
195 }
196 normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision,
197 aSign, aExp, aSig, 0, status);
198 }
199 return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign,
200 aExp, aSig, 0, status);
201 }
202
203 /*
204 * Algorithms for transcendental functions supported by MC68881 and MC68882
205 * mathematical coprocessors. The functions are derived from FPSP library.
206 */
207
208 #define one_exp 0x3FFF
209 #define one_sig UINT64_C(0x8000000000000000)
210
211 /*
212 * Function for compactifying extended double-precision floating point values.
213 */
214
215 static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig)
216 {
217 return (aExp << 16) | (aSig >> 48);
218 }
219
220 /*
221 * Log base e of x plus 1
222 */
223
224 floatx80 floatx80_lognp1(floatx80 a, float_status *status)
225 {
226 bool aSign;
227 int32_t aExp;
228 uint64_t aSig, fSig;
229
230 int8_t user_rnd_mode, user_rnd_prec;
231
232 int32_t compact, j, k;
233 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
234
235 aSig = extractFloatx80Frac(a);
236 aExp = extractFloatx80Exp(a);
237 aSign = extractFloatx80Sign(a);
238
239 if (aExp == 0x7FFF) {
240 if ((uint64_t) (aSig << 1)) {
241 propagateFloatx80NaNOneArg(a, status);
242 }
243 if (aSign) {
244 float_raise(float_flag_invalid, status);
245 return floatx80_default_nan(status);
246 }
247 return packFloatx80(0, floatx80_infinity.high, floatx80_infinity.low);
248 }
249
250 if (aExp == 0 && aSig == 0) {
251 return packFloatx80(aSign, 0, 0);
252 }
253
254 if (aSign && aExp >= one_exp) {
255 if (aExp == one_exp && aSig == one_sig) {
256 float_raise(float_flag_divbyzero, status);
257 return packFloatx80(aSign, floatx80_infinity.high,
258 floatx80_infinity.low);
259 }
260 float_raise(float_flag_invalid, status);
261 return floatx80_default_nan(status);
262 }
263
264 if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) {
265 /* <= min threshold */
266 float_raise(float_flag_inexact, status);
267 return floatx80_move(a, status);
268 }
269
270 user_rnd_mode = status->float_rounding_mode;
271 user_rnd_prec = status->floatx80_rounding_precision;
272 status->float_rounding_mode = float_round_nearest_even;
273 status->floatx80_rounding_precision = 80;
274
275 compact = floatx80_make_compact(aExp, aSig);
276
277 fp0 = a; /* Z */
278 fp1 = a;
279
280 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
281 status), status); /* X = (1+Z) */
282
283 aExp = extractFloatx80Exp(fp0);
284 aSig = extractFloatx80Frac(fp0);
285
286 compact = floatx80_make_compact(aExp, aSig);
287
288 if (compact < 0x3FFE8000 || compact > 0x3FFFC000) {
289 /* |X| < 1/2 or |X| > 3/2 */
290 k = aExp - 0x3FFF;
291 fp1 = int32_to_floatx80(k, status);
292
293 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
294 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
295
296 f = packFloatx80(0, 0x3FFF, fSig); /* F */
297 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
298
299 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
300
301 lp1cont1:
302 /* LP1CONT1 */
303 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
304 logof2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC));
305 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
306 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
307
308 fp3 = fp2;
309 fp1 = fp2;
310
311 fp1 = floatx80_mul(fp1, float64_to_floatx80(
312 make_float64(0x3FC2499AB5E4040B), status),
313 status); /* V*A6 */
314 fp2 = floatx80_mul(fp2, float64_to_floatx80(
315 make_float64(0xBFC555B5848CB7DB), status),
316 status); /* V*A5 */
317 fp1 = floatx80_add(fp1, float64_to_floatx80(
318 make_float64(0x3FC99999987D8730), status),
319 status); /* A4+V*A6 */
320 fp2 = floatx80_add(fp2, float64_to_floatx80(
321 make_float64(0xBFCFFFFFFF6F7E97), status),
322 status); /* A3+V*A5 */
323 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
324 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
325 fp1 = floatx80_add(fp1, float64_to_floatx80(
326 make_float64(0x3FD55555555555A4), status),
327 status); /* A2+V*(A4+V*A6) */
328 fp2 = floatx80_add(fp2, float64_to_floatx80(
329 make_float64(0xBFE0000000000008), status),
330 status); /* A1+V*(A3+V*A5) */
331 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
332 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
333 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
334 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
335
336 fp1 = floatx80_add(fp1, log_tbl[j + 1],
337 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
338 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
339
340 status->float_rounding_mode = user_rnd_mode;
341 status->floatx80_rounding_precision = user_rnd_prec;
342
343 a = floatx80_add(fp0, klog2, status);
344
345 float_raise(float_flag_inexact, status);
346
347 return a;
348 } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
349 /* |X| < 1/16 or |X| > -1/16 */
350 /* LP1CARE */
351 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
352 f = packFloatx80(0, 0x3FFF, fSig); /* F */
353 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
354
355 if (compact >= 0x3FFF8000) { /* 1+Z >= 1 */
356 /* KISZERO */
357 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x3F800000),
358 status), f, status); /* 1-F */
359 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (1-F)+Z */
360 fp1 = packFloatx80(0, 0, 0); /* K = 0 */
361 } else {
362 /* KISNEG */
363 fp0 = floatx80_sub(float32_to_floatx80(make_float32(0x40000000),
364 status), f, status); /* 2-F */
365 fp1 = floatx80_add(fp1, fp1, status); /* 2Z */
366 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS Y-F = (2-F)+2Z */
367 fp1 = packFloatx80(1, one_exp, one_sig); /* K = -1 */
368 }
369 goto lp1cont1;
370 } else {
371 /* LP1ONE16 */
372 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2Z */
373 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
374 status), status); /* FP0 IS 1+X */
375
376 /* LP1CONT2 */
377 fp1 = floatx80_div(fp1, fp0, status); /* U */
378 saveu = fp1;
379 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
380 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
381
382 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
383 status); /* B5 */
384 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
385 status); /* B4 */
386 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
387 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
388 fp3 = floatx80_add(fp3, float64_to_floatx80(
389 make_float64(0x3F624924928BCCFF), status),
390 status); /* B3+W*B5 */
391 fp2 = floatx80_add(fp2, float64_to_floatx80(
392 make_float64(0x3F899999999995EC), status),
393 status); /* B2+W*B4 */
394 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
395 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
396 fp1 = floatx80_add(fp1, float64_to_floatx80(
397 make_float64(0x3FB5555555555555), status),
398 status); /* B1+W*(B3+W*B5) */
399
400 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
401 fp1 = floatx80_add(fp1, fp2,
402 status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
403 fp0 = floatx80_mul(fp0, fp1,
404 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
405
406 status->float_rounding_mode = user_rnd_mode;
407 status->floatx80_rounding_precision = user_rnd_prec;
408
409 a = floatx80_add(fp0, saveu, status);
410
411 /*if (!floatx80_is_zero(a)) { */
412 float_raise(float_flag_inexact, status);
413 /*} */
414
415 return a;
416 }
417 }
418
419 /*
420 * Log base e
421 */
422
423 floatx80 floatx80_logn(floatx80 a, float_status *status)
424 {
425 bool aSign;
426 int32_t aExp;
427 uint64_t aSig, fSig;
428
429 int8_t user_rnd_mode, user_rnd_prec;
430
431 int32_t compact, j, k, adjk;
432 floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu;
433
434 aSig = extractFloatx80Frac(a);
435 aExp = extractFloatx80Exp(a);
436 aSign = extractFloatx80Sign(a);
437
438 if (aExp == 0x7FFF) {
439 if ((uint64_t) (aSig << 1)) {
440 propagateFloatx80NaNOneArg(a, status);
441 }
442 if (aSign == 0) {
443 return packFloatx80(0, floatx80_infinity.high,
444 floatx80_infinity.low);
445 }
446 }
447
448 adjk = 0;
449
450 if (aExp == 0) {
451 if (aSig == 0) { /* zero */
452 float_raise(float_flag_divbyzero, status);
453 return packFloatx80(1, floatx80_infinity.high,
454 floatx80_infinity.low);
455 }
456 if ((aSig & one_sig) == 0) { /* denormal */
457 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
458 adjk = -100;
459 aExp += 100;
460 a = packFloatx80(aSign, aExp, aSig);
461 }
462 }
463
464 if (aSign) {
465 float_raise(float_flag_invalid, status);
466 return floatx80_default_nan(status);
467 }
468
469 user_rnd_mode = status->float_rounding_mode;
470 user_rnd_prec = status->floatx80_rounding_precision;
471 status->float_rounding_mode = float_round_nearest_even;
472 status->floatx80_rounding_precision = 80;
473
474 compact = floatx80_make_compact(aExp, aSig);
475
476 if (compact < 0x3FFEF07D || compact > 0x3FFF8841) {
477 /* |X| < 15/16 or |X| > 17/16 */
478 k = aExp - 0x3FFF;
479 k += adjk;
480 fp1 = int32_to_floatx80(k, status);
481
482 fSig = (aSig & UINT64_C(0xFE00000000000000)) | UINT64_C(0x0100000000000000);
483 j = (fSig >> 56) & 0x7E; /* DISPLACEMENT FOR 1/F */
484
485 f = packFloatx80(0, 0x3FFF, fSig); /* F */
486 fp0 = packFloatx80(0, 0x3FFF, aSig); /* Y */
487
488 fp0 = floatx80_sub(fp0, f, status); /* Y-F */
489
490 /* LP1CONT1 */
491 fp0 = floatx80_mul(fp0, log_tbl[j], status); /* FP0 IS U = (Y-F)/F */
492 logof2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC));
493 klog2 = floatx80_mul(fp1, logof2, status); /* FP1 IS K*LOG2 */
494 fp2 = floatx80_mul(fp0, fp0, status); /* FP2 IS V=U*U */
495
496 fp3 = fp2;
497 fp1 = fp2;
498
499 fp1 = floatx80_mul(fp1, float64_to_floatx80(
500 make_float64(0x3FC2499AB5E4040B), status),
501 status); /* V*A6 */
502 fp2 = floatx80_mul(fp2, float64_to_floatx80(
503 make_float64(0xBFC555B5848CB7DB), status),
504 status); /* V*A5 */
505 fp1 = floatx80_add(fp1, float64_to_floatx80(
506 make_float64(0x3FC99999987D8730), status),
507 status); /* A4+V*A6 */
508 fp2 = floatx80_add(fp2, float64_to_floatx80(
509 make_float64(0xBFCFFFFFFF6F7E97), status),
510 status); /* A3+V*A5 */
511 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A4+V*A6) */
512 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A3+V*A5) */
513 fp1 = floatx80_add(fp1, float64_to_floatx80(
514 make_float64(0x3FD55555555555A4), status),
515 status); /* A2+V*(A4+V*A6) */
516 fp2 = floatx80_add(fp2, float64_to_floatx80(
517 make_float64(0xBFE0000000000008), status),
518 status); /* A1+V*(A3+V*A5) */
519 fp1 = floatx80_mul(fp1, fp3, status); /* V*(A2+V*(A4+V*A6)) */
520 fp2 = floatx80_mul(fp2, fp3, status); /* V*(A1+V*(A3+V*A5)) */
521 fp1 = floatx80_mul(fp1, fp0, status); /* U*V*(A2+V*(A4+V*A6)) */
522 fp0 = floatx80_add(fp0, fp2, status); /* U+V*(A1+V*(A3+V*A5)) */
523
524 fp1 = floatx80_add(fp1, log_tbl[j + 1],
525 status); /* LOG(F)+U*V*(A2+V*(A4+V*A6)) */
526 fp0 = floatx80_add(fp0, fp1, status); /* FP0 IS LOG(F) + LOG(1+U) */
527
528 status->float_rounding_mode = user_rnd_mode;
529 status->floatx80_rounding_precision = user_rnd_prec;
530
531 a = floatx80_add(fp0, klog2, status);
532
533 float_raise(float_flag_inexact, status);
534
535 return a;
536 } else { /* |X-1| >= 1/16 */
537 fp0 = a;
538 fp1 = a;
539 fp1 = floatx80_sub(fp1, float32_to_floatx80(make_float32(0x3F800000),
540 status), status); /* FP1 IS X-1 */
541 fp0 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
542 status), status); /* FP0 IS X+1 */
543 fp1 = floatx80_add(fp1, fp1, status); /* FP1 IS 2(X-1) */
544
545 /* LP1CONT2 */
546 fp1 = floatx80_div(fp1, fp0, status); /* U */
547 saveu = fp1;
548 fp0 = floatx80_mul(fp1, fp1, status); /* FP0 IS V = U*U */
549 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS W = V*V */
550
551 fp3 = float64_to_floatx80(make_float64(0x3F175496ADD7DAD6),
552 status); /* B5 */
553 fp2 = float64_to_floatx80(make_float64(0x3F3C71C2FE80C7E0),
554 status); /* B4 */
555 fp3 = floatx80_mul(fp3, fp1, status); /* W*B5 */
556 fp2 = floatx80_mul(fp2, fp1, status); /* W*B4 */
557 fp3 = floatx80_add(fp3, float64_to_floatx80(
558 make_float64(0x3F624924928BCCFF), status),
559 status); /* B3+W*B5 */
560 fp2 = floatx80_add(fp2, float64_to_floatx80(
561 make_float64(0x3F899999999995EC), status),
562 status); /* B2+W*B4 */
563 fp1 = floatx80_mul(fp1, fp3, status); /* W*(B3+W*B5) */
564 fp2 = floatx80_mul(fp2, fp0, status); /* V*(B2+W*B4) */
565 fp1 = floatx80_add(fp1, float64_to_floatx80(
566 make_float64(0x3FB5555555555555), status),
567 status); /* B1+W*(B3+W*B5) */
568
569 fp0 = floatx80_mul(fp0, saveu, status); /* FP0 IS U*V */
570 fp1 = floatx80_add(fp1, fp2, status); /* B1+W*(B3+W*B5) + V*(B2+W*B4) */
571 fp0 = floatx80_mul(fp0, fp1,
572 status); /* U*V*([B1+W*(B3+W*B5)] + [V*(B2+W*B4)]) */
573
574 status->float_rounding_mode = user_rnd_mode;
575 status->floatx80_rounding_precision = user_rnd_prec;
576
577 a = floatx80_add(fp0, saveu, status);
578
579 /*if (!floatx80_is_zero(a)) { */
580 float_raise(float_flag_inexact, status);
581 /*} */
582
583 return a;
584 }
585 }
586
587 /*
588 * Log base 10
589 */
590
591 floatx80 floatx80_log10(floatx80 a, float_status *status)
592 {
593 bool aSign;
594 int32_t aExp;
595 uint64_t aSig;
596
597 int8_t user_rnd_mode, user_rnd_prec;
598
599 floatx80 fp0, fp1;
600
601 aSig = extractFloatx80Frac(a);
602 aExp = extractFloatx80Exp(a);
603 aSign = extractFloatx80Sign(a);
604
605 if (aExp == 0x7FFF) {
606 if ((uint64_t) (aSig << 1)) {
607 propagateFloatx80NaNOneArg(a, status);
608 }
609 if (aSign == 0) {
610 return packFloatx80(0, floatx80_infinity.high,
611 floatx80_infinity.low);
612 }
613 }
614
615 if (aExp == 0 && aSig == 0) {
616 float_raise(float_flag_divbyzero, status);
617 return packFloatx80(1, floatx80_infinity.high,
618 floatx80_infinity.low);
619 }
620
621 if (aSign) {
622 float_raise(float_flag_invalid, status);
623 return floatx80_default_nan(status);
624 }
625
626 user_rnd_mode = status->float_rounding_mode;
627 user_rnd_prec = status->floatx80_rounding_precision;
628 status->float_rounding_mode = float_round_nearest_even;
629 status->floatx80_rounding_precision = 80;
630
631 fp0 = floatx80_logn(a, status);
632 fp1 = packFloatx80(0, 0x3FFD, UINT64_C(0xDE5BD8A937287195)); /* INV_L10 */
633
634 status->float_rounding_mode = user_rnd_mode;
635 status->floatx80_rounding_precision = user_rnd_prec;
636
637 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L10 */
638
639 float_raise(float_flag_inexact, status);
640
641 return a;
642 }
643
644 /*
645 * Log base 2
646 */
647
648 floatx80 floatx80_log2(floatx80 a, float_status *status)
649 {
650 bool aSign;
651 int32_t aExp;
652 uint64_t aSig;
653
654 int8_t user_rnd_mode, user_rnd_prec;
655
656 floatx80 fp0, fp1;
657
658 aSig = extractFloatx80Frac(a);
659 aExp = extractFloatx80Exp(a);
660 aSign = extractFloatx80Sign(a);
661
662 if (aExp == 0x7FFF) {
663 if ((uint64_t) (aSig << 1)) {
664 propagateFloatx80NaNOneArg(a, status);
665 }
666 if (aSign == 0) {
667 return packFloatx80(0, floatx80_infinity.high,
668 floatx80_infinity.low);
669 }
670 }
671
672 if (aExp == 0) {
673 if (aSig == 0) {
674 float_raise(float_flag_divbyzero, status);
675 return packFloatx80(1, floatx80_infinity.high,
676 floatx80_infinity.low);
677 }
678 normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
679 }
680
681 if (aSign) {
682 float_raise(float_flag_invalid, status);
683 return floatx80_default_nan(status);
684 }
685
686 user_rnd_mode = status->float_rounding_mode;
687 user_rnd_prec = status->floatx80_rounding_precision;
688 status->float_rounding_mode = float_round_nearest_even;
689 status->floatx80_rounding_precision = 80;
690
691 if (aSig == one_sig) { /* X is 2^k */
692 status->float_rounding_mode = user_rnd_mode;
693 status->floatx80_rounding_precision = user_rnd_prec;
694
695 a = int32_to_floatx80(aExp - 0x3FFF, status);
696 } else {
697 fp0 = floatx80_logn(a, status);
698 fp1 = packFloatx80(0, 0x3FFF, UINT64_C(0xB8AA3B295C17F0BC)); /* INV_L2 */
699
700 status->float_rounding_mode = user_rnd_mode;
701 status->floatx80_rounding_precision = user_rnd_prec;
702
703 a = floatx80_mul(fp0, fp1, status); /* LOGN(X)*INV_L2 */
704 }
705
706 float_raise(float_flag_inexact, status);
707
708 return a;
709 }
710
711 /*
712 * e to x
713 */
714
715 floatx80 floatx80_etox(floatx80 a, float_status *status)
716 {
717 bool aSign;
718 int32_t aExp;
719 uint64_t aSig;
720
721 int8_t user_rnd_mode, user_rnd_prec;
722
723 int32_t compact, n, j, k, m, m1;
724 floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale;
725 bool adjflag;
726
727 aSig = extractFloatx80Frac(a);
728 aExp = extractFloatx80Exp(a);
729 aSign = extractFloatx80Sign(a);
730
731 if (aExp == 0x7FFF) {
732 if ((uint64_t) (aSig << 1)) {
733 return propagateFloatx80NaNOneArg(a, status);
734 }
735 if (aSign) {
736 return packFloatx80(0, 0, 0);
737 }
738 return packFloatx80(0, floatx80_infinity.high,
739 floatx80_infinity.low);
740 }
741
742 if (aExp == 0 && aSig == 0) {
743 return packFloatx80(0, one_exp, one_sig);
744 }
745
746 user_rnd_mode = status->float_rounding_mode;
747 user_rnd_prec = status->floatx80_rounding_precision;
748 status->float_rounding_mode = float_round_nearest_even;
749 status->floatx80_rounding_precision = 80;
750
751 adjflag = 0;
752
753 if (aExp >= 0x3FBE) { /* |X| >= 2^(-65) */
754 compact = floatx80_make_compact(aExp, aSig);
755
756 if (compact < 0x400CB167) { /* |X| < 16380 log2 */
757 fp0 = a;
758 fp1 = a;
759 fp0 = floatx80_mul(fp0, float32_to_floatx80(
760 make_float32(0x42B8AA3B), status),
761 status); /* 64/log2 * X */
762 adjflag = 0;
763 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
764 fp0 = int32_to_floatx80(n, status);
765
766 j = n & 0x3F; /* J = N mod 64 */
767 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
768 if (n < 0 && j) {
769 /*
770 * arithmetic right shift is division and
771 * round towards minus infinity
772 */
773 m--;
774 }
775 m += 0x3FFF; /* biased exponent of 2^(M) */
776
777 expcont1:
778 fp2 = fp0; /* N */
779 fp0 = floatx80_mul(fp0, float32_to_floatx80(
780 make_float32(0xBC317218), status),
781 status); /* N * L1, L1 = lead(-log2/64) */
782 l2 = packFloatx80(0, 0x3FDC, UINT64_C(0x82E308654361C4C6));
783 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
784 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
785 fp0 = floatx80_add(fp0, fp2, status); /* R */
786
787 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
788 fp2 = float32_to_floatx80(make_float32(0x3AB60B70),
789 status); /* A5 */
790 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A5 */
791 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3C088895),
792 status), fp1,
793 status); /* fp3 is S*A4 */
794 fp2 = floatx80_add(fp2, float64_to_floatx80(make_float64(
795 0x3FA5555555554431), status),
796 status); /* fp2 is A3+S*A5 */
797 fp3 = floatx80_add(fp3, float64_to_floatx80(make_float64(
798 0x3FC5555555554018), status),
799 status); /* fp3 is A2+S*A4 */
800 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*(A3+S*A5) */
801 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 is S*(A2+S*A4) */
802 fp2 = floatx80_add(fp2, float32_to_floatx80(
803 make_float32(0x3F000000), status),
804 status); /* fp2 is A1+S*(A3+S*A5) */
805 fp3 = floatx80_mul(fp3, fp0, status); /* fp3 IS R*S*(A2+S*A4) */
806 fp2 = floatx80_mul(fp2, fp1,
807 status); /* fp2 IS S*(A1+S*(A3+S*A5)) */
808 fp0 = floatx80_add(fp0, fp3, status); /* fp0 IS R+R*S*(A2+S*A4) */
809 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
810
811 fp1 = exp_tbl[j];
812 fp0 = floatx80_mul(fp0, fp1, status); /* 2^(J/64)*(Exp(R)-1) */
813 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status),
814 status); /* accurate 2^(J/64) */
815 fp0 = floatx80_add(fp0, fp1,
816 status); /* 2^(J/64) + 2^(J/64)*(Exp(R)-1) */
817
818 scale = packFloatx80(0, m, one_sig);
819 if (adjflag) {
820 adjscale = packFloatx80(0, m1, one_sig);
821 fp0 = floatx80_mul(fp0, adjscale, status);
822 }
823
824 status->float_rounding_mode = user_rnd_mode;
825 status->floatx80_rounding_precision = user_rnd_prec;
826
827 a = floatx80_mul(fp0, scale, status);
828
829 float_raise(float_flag_inexact, status);
830
831 return a;
832 } else { /* |X| >= 16380 log2 */
833 if (compact > 0x400CB27C) { /* |X| >= 16480 log2 */
834 status->float_rounding_mode = user_rnd_mode;
835 status->floatx80_rounding_precision = user_rnd_prec;
836 if (aSign) {
837 a = roundAndPackFloatx80(
838 status->floatx80_rounding_precision,
839 0, -0x1000, aSig, 0, status);
840 } else {
841 a = roundAndPackFloatx80(
842 status->floatx80_rounding_precision,
843 0, 0x8000, aSig, 0, status);
844 }
845 float_raise(float_flag_inexact, status);
846
847 return a;
848 } else {
849 fp0 = a;
850 fp1 = a;
851 fp0 = floatx80_mul(fp0, float32_to_floatx80(
852 make_float32(0x42B8AA3B), status),
853 status); /* 64/log2 * X */
854 adjflag = 1;
855 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
856 fp0 = int32_to_floatx80(n, status);
857
858 j = n & 0x3F; /* J = N mod 64 */
859 /* NOTE: this is really arithmetic right shift by 6 */
860 k = n / 64;
861 if (n < 0 && j) {
862 /* arithmetic right shift is division and
863 * round towards minus infinity
864 */
865 k--;
866 }
867 /* NOTE: this is really arithmetic right shift by 1 */
868 m1 = k / 2;
869 if (k < 0 && (k & 1)) {
870 /* arithmetic right shift is division and
871 * round towards minus infinity
872 */
873 m1--;
874 }
875 m = k - m1;
876 m1 += 0x3FFF; /* biased exponent of 2^(M1) */
877 m += 0x3FFF; /* biased exponent of 2^(M) */
878
879 goto expcont1;
880 }
881 }
882 } else { /* |X| < 2^(-65) */
883 status->float_rounding_mode = user_rnd_mode;
884 status->floatx80_rounding_precision = user_rnd_prec;
885
886 a = floatx80_add(a, float32_to_floatx80(make_float32(0x3F800000),
887 status), status); /* 1 + X */
888
889 float_raise(float_flag_inexact, status);
890
891 return a;
892 }
893 }
894
895 /*
896 * 2 to x
897 */
898
899 floatx80 floatx80_twotox(floatx80 a, float_status *status)
900 {
901 bool aSign;
902 int32_t aExp;
903 uint64_t aSig;
904
905 int8_t user_rnd_mode, user_rnd_prec;
906
907 int32_t compact, n, j, l, m, m1;
908 floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
909
910 aSig = extractFloatx80Frac(a);
911 aExp = extractFloatx80Exp(a);
912 aSign = extractFloatx80Sign(a);
913
914 if (aExp == 0x7FFF) {
915 if ((uint64_t) (aSig << 1)) {
916 return propagateFloatx80NaNOneArg(a, status);
917 }
918 if (aSign) {
919 return packFloatx80(0, 0, 0);
920 }
921 return packFloatx80(0, floatx80_infinity.high,
922 floatx80_infinity.low);
923 }
924
925 if (aExp == 0 && aSig == 0) {
926 return packFloatx80(0, one_exp, one_sig);
927 }
928
929 user_rnd_mode = status->float_rounding_mode;
930 user_rnd_prec = status->floatx80_rounding_precision;
931 status->float_rounding_mode = float_round_nearest_even;
932 status->floatx80_rounding_precision = 80;
933
934 fp0 = a;
935
936 compact = floatx80_make_compact(aExp, aSig);
937
938 if (compact < 0x3FB98000 || compact > 0x400D80C0) {
939 /* |X| > 16480 or |X| < 2^(-70) */
940 if (compact > 0x3FFF8000) { /* |X| > 16480 */
941 status->float_rounding_mode = user_rnd_mode;
942 status->floatx80_rounding_precision = user_rnd_prec;
943
944 if (aSign) {
945 return roundAndPackFloatx80(status->floatx80_rounding_precision,
946 0, -0x1000, aSig, 0, status);
947 } else {
948 return roundAndPackFloatx80(status->floatx80_rounding_precision,
949 0, 0x8000, aSig, 0, status);
950 }
951 } else { /* |X| < 2^(-70) */
952 status->float_rounding_mode = user_rnd_mode;
953 status->floatx80_rounding_precision = user_rnd_prec;
954
955 a = floatx80_add(fp0, float32_to_floatx80(
956 make_float32(0x3F800000), status),
957 status); /* 1 + X */
958
959 float_raise(float_flag_inexact, status);
960
961 return a;
962 }
963 } else { /* 2^(-70) <= |X| <= 16480 */
964 fp1 = fp0; /* X */
965 fp1 = floatx80_mul(fp1, float32_to_floatx80(
966 make_float32(0x42800000), status),
967 status); /* X * 64 */
968 n = floatx80_to_int32(fp1, status);
969 fp1 = int32_to_floatx80(n, status);
970 j = n & 0x3F;
971 l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
972 if (n < 0 && j) {
973 /*
974 * arithmetic right shift is division and
975 * round towards minus infinity
976 */
977 l--;
978 }
979 m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */
980 if (l < 0 && (l & 1)) {
981 /*
982 * arithmetic right shift is division and
983 * round towards minus infinity
984 */
985 m--;
986 }
987 m1 = l - m;
988 m1 += 0x3FFF; /* ADJFACT IS 2^(M') */
989
990 adjfact = packFloatx80(0, m1, one_sig);
991 fact1 = exp2_tbl[j];
992 fact1.high += m;
993 fact2.high = exp2_tbl2[j] >> 16;
994 fact2.high += m;
995 fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF);
996 fact2.low <<= 48;
997
998 fp1 = floatx80_mul(fp1, float32_to_floatx80(
999 make_float32(0x3C800000), status),
1000 status); /* (1/64)*N */
1001 fp0 = floatx80_sub(fp0, fp1, status); /* X - (1/64)*INT(64 X) */
1002 fp2 = packFloatx80(0, 0x3FFE, UINT64_C(0xB17217F7D1CF79AC)); /* LOG2 */
1003 fp0 = floatx80_mul(fp0, fp2, status); /* R */
1004
1005 /* EXPR */
1006 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1007 fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2),
1008 status); /* A5 */
1009 fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C),
1010 status); /* A4 */
1011 fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */
1012 fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */
1013 fp2 = floatx80_add(fp2, float64_to_floatx80(
1014 make_float64(0x3FA5555555554CC1), status),
1015 status); /* A3+S*A5 */
1016 fp3 = floatx80_add(fp3, float64_to_floatx80(
1017 make_float64(0x3FC5555555554A54), status),
1018 status); /* A2+S*A4 */
1019 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */
1020 fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */
1021 fp2 = floatx80_add(fp2, float64_to_floatx80(
1022 make_float64(0x3FE0000000000000), status),
1023 status); /* A1+S*(A3+S*A5) */
1024 fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */
1025
1026 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */
1027 fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */
1028 fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */
1029
1030 fp0 = floatx80_mul(fp0, fact1, status);
1031 fp0 = floatx80_add(fp0, fact2, status);
1032 fp0 = floatx80_add(fp0, fact1, status);
1033
1034 status->float_rounding_mode = user_rnd_mode;
1035 status->floatx80_rounding_precision = user_rnd_prec;
1036
1037 a = floatx80_mul(fp0, adjfact, status);
1038
1039 float_raise(float_flag_inexact, status);
1040
1041 return a;
1042 }
1043 }
1044
1045 /*
1046 * 10 to x
1047 */
1048
1049 floatx80 floatx80_tentox(floatx80 a, float_status *status)
1050 {
1051 bool aSign;
1052 int32_t aExp;
1053 uint64_t aSig;
1054
1055 int8_t user_rnd_mode, user_rnd_prec;
1056
1057 int32_t compact, n, j, l, m, m1;
1058 floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2;
1059
1060 aSig = extractFloatx80Frac(a);
1061 aExp = extractFloatx80Exp(a);
1062 aSign = extractFloatx80Sign(a);
1063
1064 if (aExp == 0x7FFF) {
1065 if ((uint64_t) (aSig << 1)) {
1066 return propagateFloatx80NaNOneArg(a, status);
1067 }
1068 if (aSign) {
1069 return packFloatx80(0, 0, 0);
1070 }
1071 return packFloatx80(0, floatx80_infinity.high,
1072 floatx80_infinity.low);
1073 }
1074
1075 if (aExp == 0 && aSig == 0) {
1076 return packFloatx80(0, one_exp, one_sig);
1077 }
1078
1079 user_rnd_mode = status->float_rounding_mode;
1080 user_rnd_prec = status->floatx80_rounding_precision;
1081 status->float_rounding_mode = float_round_nearest_even;
1082 status->floatx80_rounding_precision = 80;
1083
1084 fp0 = a;
1085
1086 compact = floatx80_make_compact(aExp, aSig);
1087
1088 if (compact < 0x3FB98000 || compact > 0x400B9B07) {
1089 /* |X| > 16480 LOG2/LOG10 or |X| < 2^(-70) */
1090 if (compact > 0x3FFF8000) { /* |X| > 16480 */
1091 status->float_rounding_mode = user_rnd_mode;
1092 status->floatx80_rounding_precision = user_rnd_prec;
1093
1094 if (aSign) {
1095 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1096 0, -0x1000, aSig, 0, status);
1097 } else {
1098 return roundAndPackFloatx80(status->floatx80_rounding_precision,
1099 0, 0x8000, aSig, 0, status);
1100 }
1101 } else { /* |X| < 2^(-70) */
1102 status->float_rounding_mode = user_rnd_mode;
1103 status->floatx80_rounding_precision = user_rnd_prec;
1104
1105 a = floatx80_add(fp0, float32_to_floatx80(
1106 make_float32(0x3F800000), status),
1107 status); /* 1 + X */
1108
1109 float_raise(float_flag_inexact, status);
1110
1111 return a;
1112 }
1113 } else { /* 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10 */
1114 fp1 = fp0; /* X */
1115 fp1 = floatx80_mul(fp1, float64_to_floatx80(
1116 make_float64(0x406A934F0979A371),
1117 status), status); /* X*64*LOG10/LOG2 */
1118 n = floatx80_to_int32(fp1, status); /* N=INT(X*64*LOG10/LOG2) */
1119 fp1 = int32_to_floatx80(n, status);
1120
1121 j = n & 0x3F;
1122 l = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
1123 if (n < 0 && j) {
1124 /*
1125 * arithmetic right shift is division and
1126 * round towards minus infinity
1127 */
1128 l--;
1129 }
1130 m = l / 2; /* NOTE: this is really arithmetic right shift by 1 */
1131 if (l < 0 && (l & 1)) {
1132 /*
1133 * arithmetic right shift is division and
1134 * round towards minus infinity
1135 */
1136 m--;
1137 }
1138 m1 = l - m;
1139 m1 += 0x3FFF; /* ADJFACT IS 2^(M') */
1140
1141 adjfact = packFloatx80(0, m1, one_sig);
1142 fact1 = exp2_tbl[j];
1143 fact1.high += m;
1144 fact2.high = exp2_tbl2[j] >> 16;
1145 fact2.high += m;
1146 fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF);
1147 fact2.low <<= 48;
1148
1149 fp2 = fp1; /* N */
1150 fp1 = floatx80_mul(fp1, float64_to_floatx80(
1151 make_float64(0x3F734413509F8000), status),
1152 status); /* N*(LOG2/64LOG10)_LEAD */
1153 fp3 = packFloatx80(1, 0x3FCD, UINT64_C(0xC0219DC1DA994FD2));
1154 fp2 = floatx80_mul(fp2, fp3, status); /* N*(LOG2/64LOG10)_TRAIL */
1155 fp0 = floatx80_sub(fp0, fp1, status); /* X - N L_LEAD */
1156 fp0 = floatx80_sub(fp0, fp2, status); /* X - N L_TRAIL */
1157 fp2 = packFloatx80(0, 0x4000, UINT64_C(0x935D8DDDAAA8AC17)); /* LOG10 */
1158 fp0 = floatx80_mul(fp0, fp2, status); /* R */
1159
1160 /* EXPR */
1161 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1162 fp2 = float64_to_floatx80(make_float64(0x3F56C16D6F7BD0B2),
1163 status); /* A5 */
1164 fp3 = float64_to_floatx80(make_float64(0x3F811112302C712C),
1165 status); /* A4 */
1166 fp2 = floatx80_mul(fp2, fp1, status); /* S*A5 */
1167 fp3 = floatx80_mul(fp3, fp1, status); /* S*A4 */
1168 fp2 = floatx80_add(fp2, float64_to_floatx80(
1169 make_float64(0x3FA5555555554CC1), status),
1170 status); /* A3+S*A5 */
1171 fp3 = floatx80_add(fp3, float64_to_floatx80(
1172 make_float64(0x3FC5555555554A54), status),
1173 status); /* A2+S*A4 */
1174 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A3+S*A5) */
1175 fp3 = floatx80_mul(fp3, fp1, status); /* S*(A2+S*A4) */
1176 fp2 = floatx80_add(fp2, float64_to_floatx80(
1177 make_float64(0x3FE0000000000000), status),
1178 status); /* A1+S*(A3+S*A5) */
1179 fp3 = floatx80_mul(fp3, fp0, status); /* R*S*(A2+S*A4) */
1180
1181 fp2 = floatx80_mul(fp2, fp1, status); /* S*(A1+S*(A3+S*A5)) */
1182 fp0 = floatx80_add(fp0, fp3, status); /* R+R*S*(A2+S*A4) */
1183 fp0 = floatx80_add(fp0, fp2, status); /* EXP(R) - 1 */
1184
1185 fp0 = floatx80_mul(fp0, fact1, status);
1186 fp0 = floatx80_add(fp0, fact2, status);
1187 fp0 = floatx80_add(fp0, fact1, status);
1188
1189 status->float_rounding_mode = user_rnd_mode;
1190 status->floatx80_rounding_precision = user_rnd_prec;
1191
1192 a = floatx80_mul(fp0, adjfact, status);
1193
1194 float_raise(float_flag_inexact, status);
1195
1196 return a;
1197 }
1198 }
1199
1200 /*
1201 * Tangent
1202 */
1203
1204 floatx80 floatx80_tan(floatx80 a, float_status *status)
1205 {
1206 bool aSign, xSign;
1207 int32_t aExp, xExp;
1208 uint64_t aSig, xSig;
1209
1210 int8_t user_rnd_mode, user_rnd_prec;
1211
1212 int32_t compact, l, n, j;
1213 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2;
1214 float32 twoto63;
1215 bool endflag;
1216
1217 aSig = extractFloatx80Frac(a);
1218 aExp = extractFloatx80Exp(a);
1219 aSign = extractFloatx80Sign(a);
1220
1221 if (aExp == 0x7FFF) {
1222 if ((uint64_t) (aSig << 1)) {
1223 return propagateFloatx80NaNOneArg(a, status);
1224 }
1225 float_raise(float_flag_invalid, status);
1226 return floatx80_default_nan(status);
1227 }
1228
1229 if (aExp == 0 && aSig == 0) {
1230 return packFloatx80(aSign, 0, 0);
1231 }
1232
1233 user_rnd_mode = status->float_rounding_mode;
1234 user_rnd_prec = status->floatx80_rounding_precision;
1235 status->float_rounding_mode = float_round_nearest_even;
1236 status->floatx80_rounding_precision = 80;
1237
1238 compact = floatx80_make_compact(aExp, aSig);
1239
1240 fp0 = a;
1241
1242 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1243 /* 2^(-40) > |X| > 15 PI */
1244 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1245 /* REDUCEX */
1246 fp1 = packFloatx80(0, 0, 0);
1247 if (compact == 0x7FFEFFFF) {
1248 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1249 UINT64_C(0xC90FDAA200000000));
1250 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1251 UINT64_C(0x85A308D300000000));
1252 fp0 = floatx80_add(fp0, twopi1, status);
1253 fp1 = fp0;
1254 fp0 = floatx80_add(fp0, twopi2, status);
1255 fp1 = floatx80_sub(fp1, fp0, status);
1256 fp1 = floatx80_add(fp1, twopi2, status);
1257 }
1258 loop:
1259 xSign = extractFloatx80Sign(fp0);
1260 xExp = extractFloatx80Exp(fp0);
1261 xExp -= 0x3FFF;
1262 if (xExp <= 28) {
1263 l = 0;
1264 endflag = true;
1265 } else {
1266 l = xExp - 27;
1267 endflag = false;
1268 }
1269 invtwopi = packFloatx80(0, 0x3FFE - l,
1270 UINT64_C(0xA2F9836E4E44152A)); /* INVTWOPI */
1271 twopi1 = packFloatx80(0, 0x3FFF + l, UINT64_C(0xC90FDAA200000000));
1272 twopi2 = packFloatx80(0, 0x3FDD + l, UINT64_C(0x85A308D300000000));
1273
1274 /* SIGN(INARG)*2^63 IN SGL */
1275 twoto63 = packFloat32(xSign, 0xBE, 0);
1276
1277 fp2 = floatx80_mul(fp0, invtwopi, status);
1278 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1279 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1280 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1281 status); /* FP2 is N */
1282 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1283 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1284 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1285 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1286 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1287 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1288 fp3 = fp0; /* FP3 is A */
1289 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1290 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1291
1292 if (endflag) {
1293 n = floatx80_to_int32(fp2, status);
1294 goto tancont;
1295 }
1296 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1297 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1298 goto loop;
1299 } else {
1300 status->float_rounding_mode = user_rnd_mode;
1301 status->floatx80_rounding_precision = user_rnd_prec;
1302
1303 a = floatx80_move(a, status);
1304
1305 float_raise(float_flag_inexact, status);
1306
1307 return a;
1308 }
1309 } else {
1310 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1311 make_float64(0x3FE45F306DC9C883), status),
1312 status); /* X*2/PI */
1313
1314 n = floatx80_to_int32(fp1, status);
1315 j = 32 + n;
1316
1317 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1318 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1319 status); /* FP0 IS R = (X-Y1)-Y2 */
1320
1321 tancont:
1322 if (n & 1) {
1323 /* NODD */
1324 fp1 = fp0; /* R */
1325 fp0 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1326 fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688),
1327 status); /* Q4 */
1328 fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04),
1329 status); /* P3 */
1330 fp3 = floatx80_mul(fp3, fp0, status); /* SQ4 */
1331 fp2 = floatx80_mul(fp2, fp0, status); /* SP3 */
1332 fp3 = floatx80_add(fp3, float64_to_floatx80(
1333 make_float64(0xBF346F59B39BA65F), status),
1334 status); /* Q3+SQ4 */
1335 fp4 = packFloatx80(0, 0x3FF6, UINT64_C(0xE073D3FC199C4A00));
1336 fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */
1337 fp3 = floatx80_mul(fp3, fp0, status); /* S(Q3+SQ4) */
1338 fp2 = floatx80_mul(fp2, fp0, status); /* S(P2+SP3) */
1339 fp4 = packFloatx80(0, 0x3FF9, UINT64_C(0xD23CD68415D95FA1));
1340 fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */
1341 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0x8895A6C5FB423BCA));
1342 fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */
1343 fp3 = floatx80_mul(fp3, fp0, status); /* S(Q2+S(Q3+SQ4)) */
1344 fp2 = floatx80_mul(fp2, fp0, status); /* S(P1+S(P2+SP3)) */
1345 fp4 = packFloatx80(1, 0x3FFD, UINT64_C(0xEEF57E0DA84BC8CE));
1346 fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */
1347 fp2 = floatx80_mul(fp2, fp1, status); /* RS(P1+S(P2+SP3)) */
1348 fp0 = floatx80_mul(fp0, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */
1349 fp1 = floatx80_add(fp1, fp2, status); /* R+RS(P1+S(P2+SP3)) */
1350 fp0 = floatx80_add(fp0, float32_to_floatx80(
1351 make_float32(0x3F800000), status),
1352 status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */
1353
1354 xSign = extractFloatx80Sign(fp1);
1355 xExp = extractFloatx80Exp(fp1);
1356 xSig = extractFloatx80Frac(fp1);
1357 xSign ^= 1;
1358 fp1 = packFloatx80(xSign, xExp, xSig);
1359
1360 status->float_rounding_mode = user_rnd_mode;
1361 status->floatx80_rounding_precision = user_rnd_prec;
1362
1363 a = floatx80_div(fp0, fp1, status);
1364
1365 float_raise(float_flag_inexact, status);
1366
1367 return a;
1368 } else {
1369 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
1370 fp3 = float64_to_floatx80(make_float64(0x3EA0B759F50F8688),
1371 status); /* Q4 */
1372 fp2 = float64_to_floatx80(make_float64(0xBEF2BAA5A8924F04),
1373 status); /* P3 */
1374 fp3 = floatx80_mul(fp3, fp1, status); /* SQ4 */
1375 fp2 = floatx80_mul(fp2, fp1, status); /* SP3 */
1376 fp3 = floatx80_add(fp3, float64_to_floatx80(
1377 make_float64(0xBF346F59B39BA65F), status),
1378 status); /* Q3+SQ4 */
1379 fp4 = packFloatx80(0, 0x3FF6, UINT64_C(0xE073D3FC199C4A00));
1380 fp2 = floatx80_add(fp2, fp4, status); /* P2+SP3 */
1381 fp3 = floatx80_mul(fp3, fp1, status); /* S(Q3+SQ4) */
1382 fp2 = floatx80_mul(fp2, fp1, status); /* S(P2+SP3) */
1383 fp4 = packFloatx80(0, 0x3FF9, UINT64_C(0xD23CD68415D95FA1));
1384 fp3 = floatx80_add(fp3, fp4, status); /* Q2+S(Q3+SQ4) */
1385 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0x8895A6C5FB423BCA));
1386 fp2 = floatx80_add(fp2, fp4, status); /* P1+S(P2+SP3) */
1387 fp3 = floatx80_mul(fp3, fp1, status); /* S(Q2+S(Q3+SQ4)) */
1388 fp2 = floatx80_mul(fp2, fp1, status); /* S(P1+S(P2+SP3)) */
1389 fp4 = packFloatx80(1, 0x3FFD, UINT64_C(0xEEF57E0DA84BC8CE));
1390 fp3 = floatx80_add(fp3, fp4, status); /* Q1+S(Q2+S(Q3+SQ4)) */
1391 fp2 = floatx80_mul(fp2, fp0, status); /* RS(P1+S(P2+SP3)) */
1392 fp1 = floatx80_mul(fp1, fp3, status); /* S(Q1+S(Q2+S(Q3+SQ4))) */
1393 fp0 = floatx80_add(fp0, fp2, status); /* R+RS(P1+S(P2+SP3)) */
1394 fp1 = floatx80_add(fp1, float32_to_floatx80(
1395 make_float32(0x3F800000), status),
1396 status); /* 1+S(Q1+S(Q2+S(Q3+SQ4))) */
1397
1398 status->float_rounding_mode = user_rnd_mode;
1399 status->floatx80_rounding_precision = user_rnd_prec;
1400
1401 a = floatx80_div(fp0, fp1, status);
1402
1403 float_raise(float_flag_inexact, status);
1404
1405 return a;
1406 }
1407 }
1408 }
1409
1410 /*
1411 * Sine
1412 */
1413
1414 floatx80 floatx80_sin(floatx80 a, float_status *status)
1415 {
1416 bool aSign, xSign;
1417 int32_t aExp, xExp;
1418 uint64_t aSig, xSig;
1419
1420 int8_t user_rnd_mode, user_rnd_prec;
1421
1422 int32_t compact, l, n, j;
1423 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1424 float32 posneg1, twoto63;
1425 bool endflag;
1426
1427 aSig = extractFloatx80Frac(a);
1428 aExp = extractFloatx80Exp(a);
1429 aSign = extractFloatx80Sign(a);
1430
1431 if (aExp == 0x7FFF) {
1432 if ((uint64_t) (aSig << 1)) {
1433 return propagateFloatx80NaNOneArg(a, status);
1434 }
1435 float_raise(float_flag_invalid, status);
1436 return floatx80_default_nan(status);
1437 }
1438
1439 if (aExp == 0 && aSig == 0) {
1440 return packFloatx80(aSign, 0, 0);
1441 }
1442
1443 user_rnd_mode = status->float_rounding_mode;
1444 user_rnd_prec = status->floatx80_rounding_precision;
1445 status->float_rounding_mode = float_round_nearest_even;
1446 status->floatx80_rounding_precision = 80;
1447
1448 compact = floatx80_make_compact(aExp, aSig);
1449
1450 fp0 = a;
1451
1452 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1453 /* 2^(-40) > |X| > 15 PI */
1454 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1455 /* REDUCEX */
1456 fp1 = packFloatx80(0, 0, 0);
1457 if (compact == 0x7FFEFFFF) {
1458 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1459 UINT64_C(0xC90FDAA200000000));
1460 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1461 UINT64_C(0x85A308D300000000));
1462 fp0 = floatx80_add(fp0, twopi1, status);
1463 fp1 = fp0;
1464 fp0 = floatx80_add(fp0, twopi2, status);
1465 fp1 = floatx80_sub(fp1, fp0, status);
1466 fp1 = floatx80_add(fp1, twopi2, status);
1467 }
1468 loop:
1469 xSign = extractFloatx80Sign(fp0);
1470 xExp = extractFloatx80Exp(fp0);
1471 xExp -= 0x3FFF;
1472 if (xExp <= 28) {
1473 l = 0;
1474 endflag = true;
1475 } else {
1476 l = xExp - 27;
1477 endflag = false;
1478 }
1479 invtwopi = packFloatx80(0, 0x3FFE - l,
1480 UINT64_C(0xA2F9836E4E44152A)); /* INVTWOPI */
1481 twopi1 = packFloatx80(0, 0x3FFF + l, UINT64_C(0xC90FDAA200000000));
1482 twopi2 = packFloatx80(0, 0x3FDD + l, UINT64_C(0x85A308D300000000));
1483
1484 /* SIGN(INARG)*2^63 IN SGL */
1485 twoto63 = packFloat32(xSign, 0xBE, 0);
1486
1487 fp2 = floatx80_mul(fp0, invtwopi, status);
1488 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1489 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1490 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1491 status); /* FP2 is N */
1492 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1493 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1494 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1495 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1496 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1497 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1498 fp3 = fp0; /* FP3 is A */
1499 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1500 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1501
1502 if (endflag) {
1503 n = floatx80_to_int32(fp2, status);
1504 goto sincont;
1505 }
1506 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1507 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1508 goto loop;
1509 } else {
1510 /* SINSM */
1511 fp0 = float32_to_floatx80(make_float32(0x3F800000),
1512 status); /* 1 */
1513
1514 status->float_rounding_mode = user_rnd_mode;
1515 status->floatx80_rounding_precision = user_rnd_prec;
1516
1517 /* SINTINY */
1518 a = floatx80_move(a, status);
1519 float_raise(float_flag_inexact, status);
1520
1521 return a;
1522 }
1523 } else {
1524 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1525 make_float64(0x3FE45F306DC9C883), status),
1526 status); /* X*2/PI */
1527
1528 n = floatx80_to_int32(fp1, status);
1529 j = 32 + n;
1530
1531 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1532 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1533 status); /* FP0 IS R = (X-Y1)-Y2 */
1534
1535 sincont:
1536 if (n & 1) {
1537 /* COSPOLY */
1538 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1539 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1540 fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3),
1541 status); /* B8 */
1542 fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19),
1543 status); /* B7 */
1544
1545 xSign = extractFloatx80Sign(fp0); /* X IS S */
1546 xExp = extractFloatx80Exp(fp0);
1547 xSig = extractFloatx80Frac(fp0);
1548
1549 if ((n >> 1) & 1) {
1550 xSign ^= 1;
1551 posneg1 = make_float32(0xBF800000); /* -1 */
1552 } else {
1553 xSign ^= 0;
1554 posneg1 = make_float32(0x3F800000); /* 1 */
1555 } /* X IS NOW R'= SGN*R */
1556
1557 fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */
1558 fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */
1559 fp2 = floatx80_add(fp2, float64_to_floatx80(
1560 make_float64(0x3E21EED90612C972), status),
1561 status); /* B6+TB8 */
1562 fp3 = floatx80_add(fp3, float64_to_floatx80(
1563 make_float64(0xBE927E4FB79D9FCF), status),
1564 status); /* B5+TB7 */
1565 fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */
1566 fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */
1567 fp2 = floatx80_add(fp2, float64_to_floatx80(
1568 make_float64(0x3EFA01A01A01D423), status),
1569 status); /* B4+T(B6+TB8) */
1570 fp4 = packFloatx80(1, 0x3FF5, UINT64_C(0xB60B60B60B61D438));
1571 fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */
1572 fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */
1573 fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */
1574 fp4 = packFloatx80(0, 0x3FFA, UINT64_C(0xAAAAAAAAAAAAAB5E));
1575 fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */
1576 fp1 = floatx80_add(fp1, float32_to_floatx80(
1577 make_float32(0xBF000000), status),
1578 status); /* B1+T(B3+T(B5+TB7)) */
1579 fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */
1580 fp0 = floatx80_add(fp0, fp1, status); /* [B1+T(B3+T(B5+TB7))]+
1581 * [S(B2+T(B4+T(B6+TB8)))]
1582 */
1583
1584 x = packFloatx80(xSign, xExp, xSig);
1585 fp0 = floatx80_mul(fp0, x, status);
1586
1587 status->float_rounding_mode = user_rnd_mode;
1588 status->floatx80_rounding_precision = user_rnd_prec;
1589
1590 a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status);
1591
1592 float_raise(float_flag_inexact, status);
1593
1594 return a;
1595 } else {
1596 /* SINPOLY */
1597 xSign = extractFloatx80Sign(fp0); /* X IS R */
1598 xExp = extractFloatx80Exp(fp0);
1599 xSig = extractFloatx80Frac(fp0);
1600
1601 xSign ^= (n >> 1) & 1; /* X IS NOW R'= SGN*R */
1602
1603 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1604 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1605 fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5),
1606 status); /* A7 */
1607 fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1),
1608 status); /* A6 */
1609 fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */
1610 fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */
1611 fp3 = floatx80_add(fp3, float64_to_floatx80(
1612 make_float64(0xBE5AE6452A118AE4), status),
1613 status); /* A5+T*A7 */
1614 fp2 = floatx80_add(fp2, float64_to_floatx80(
1615 make_float64(0x3EC71DE3A5341531), status),
1616 status); /* A4+T*A6 */
1617 fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */
1618 fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */
1619 fp3 = floatx80_add(fp3, float64_to_floatx80(
1620 make_float64(0xBF2A01A01A018B59), status),
1621 status); /* A3+T(A5+TA7) */
1622 fp4 = packFloatx80(0, 0x3FF8, UINT64_C(0x88888888888859AF));
1623 fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */
1624 fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */
1625 fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */
1626 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAA99));
1627 fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */
1628 fp1 = floatx80_add(fp1, fp2,
1629 status); /* [A1+T(A3+T(A5+TA7))]+
1630 * [S(A2+T(A4+TA6))]
1631 */
1632
1633 x = packFloatx80(xSign, xExp, xSig);
1634 fp0 = floatx80_mul(fp0, x, status); /* R'*S */
1635 fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */
1636
1637 status->float_rounding_mode = user_rnd_mode;
1638 status->floatx80_rounding_precision = user_rnd_prec;
1639
1640 a = floatx80_add(fp0, x, status);
1641
1642 float_raise(float_flag_inexact, status);
1643
1644 return a;
1645 }
1646 }
1647 }
1648
1649 /*
1650 * Cosine
1651 */
1652
1653 floatx80 floatx80_cos(floatx80 a, float_status *status)
1654 {
1655 bool aSign, xSign;
1656 int32_t aExp, xExp;
1657 uint64_t aSig, xSig;
1658
1659 int8_t user_rnd_mode, user_rnd_prec;
1660
1661 int32_t compact, l, n, j;
1662 floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2;
1663 float32 posneg1, twoto63;
1664 bool endflag;
1665
1666 aSig = extractFloatx80Frac(a);
1667 aExp = extractFloatx80Exp(a);
1668 aSign = extractFloatx80Sign(a);
1669
1670 if (aExp == 0x7FFF) {
1671 if ((uint64_t) (aSig << 1)) {
1672 return propagateFloatx80NaNOneArg(a, status);
1673 }
1674 float_raise(float_flag_invalid, status);
1675 return floatx80_default_nan(status);
1676 }
1677
1678 if (aExp == 0 && aSig == 0) {
1679 return packFloatx80(0, one_exp, one_sig);
1680 }
1681
1682 user_rnd_mode = status->float_rounding_mode;
1683 user_rnd_prec = status->floatx80_rounding_precision;
1684 status->float_rounding_mode = float_round_nearest_even;
1685 status->floatx80_rounding_precision = 80;
1686
1687 compact = floatx80_make_compact(aExp, aSig);
1688
1689 fp0 = a;
1690
1691 if (compact < 0x3FD78000 || compact > 0x4004BC7E) {
1692 /* 2^(-40) > |X| > 15 PI */
1693 if (compact > 0x3FFF8000) { /* |X| >= 15 PI */
1694 /* REDUCEX */
1695 fp1 = packFloatx80(0, 0, 0);
1696 if (compact == 0x7FFEFFFF) {
1697 twopi1 = packFloatx80(aSign ^ 1, 0x7FFE,
1698 UINT64_C(0xC90FDAA200000000));
1699 twopi2 = packFloatx80(aSign ^ 1, 0x7FDC,
1700 UINT64_C(0x85A308D300000000));
1701 fp0 = floatx80_add(fp0, twopi1, status);
1702 fp1 = fp0;
1703 fp0 = floatx80_add(fp0, twopi2, status);
1704 fp1 = floatx80_sub(fp1, fp0, status);
1705 fp1 = floatx80_add(fp1, twopi2, status);
1706 }
1707 loop:
1708 xSign = extractFloatx80Sign(fp0);
1709 xExp = extractFloatx80Exp(fp0);
1710 xExp -= 0x3FFF;
1711 if (xExp <= 28) {
1712 l = 0;
1713 endflag = true;
1714 } else {
1715 l = xExp - 27;
1716 endflag = false;
1717 }
1718 invtwopi = packFloatx80(0, 0x3FFE - l,
1719 UINT64_C(0xA2F9836E4E44152A)); /* INVTWOPI */
1720 twopi1 = packFloatx80(0, 0x3FFF + l, UINT64_C(0xC90FDAA200000000));
1721 twopi2 = packFloatx80(0, 0x3FDD + l, UINT64_C(0x85A308D300000000));
1722
1723 /* SIGN(INARG)*2^63 IN SGL */
1724 twoto63 = packFloat32(xSign, 0xBE, 0);
1725
1726 fp2 = floatx80_mul(fp0, invtwopi, status);
1727 fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status),
1728 status); /* THE FRACT PART OF FP2 IS ROUNDED */
1729 fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status),
1730 status); /* FP2 is N */
1731 fp4 = floatx80_mul(twopi1, fp2, status); /* W = N*P1 */
1732 fp5 = floatx80_mul(twopi2, fp2, status); /* w = N*P2 */
1733 fp3 = floatx80_add(fp4, fp5, status); /* FP3 is P */
1734 fp4 = floatx80_sub(fp4, fp3, status); /* W-P */
1735 fp0 = floatx80_sub(fp0, fp3, status); /* FP0 is A := R - P */
1736 fp4 = floatx80_add(fp4, fp5, status); /* FP4 is p = (W-P)+w */
1737 fp3 = fp0; /* FP3 is A */
1738 fp1 = floatx80_sub(fp1, fp4, status); /* FP1 is a := r - p */
1739 fp0 = floatx80_add(fp0, fp1, status); /* FP0 is R := A+a */
1740
1741 if (endflag) {
1742 n = floatx80_to_int32(fp2, status);
1743 goto sincont;
1744 }
1745 fp3 = floatx80_sub(fp3, fp0, status); /* A-R */
1746 fp1 = floatx80_add(fp1, fp3, status); /* FP1 is r := (A-R)+a */
1747 goto loop;
1748 } else {
1749 /* SINSM */
1750 fp0 = float32_to_floatx80(make_float32(0x3F800000), status); /* 1 */
1751
1752 status->float_rounding_mode = user_rnd_mode;
1753 status->floatx80_rounding_precision = user_rnd_prec;
1754
1755 /* COSTINY */
1756 a = floatx80_sub(fp0, float32_to_floatx80(
1757 make_float32(0x00800000), status),
1758 status);
1759 float_raise(float_flag_inexact, status);
1760
1761 return a;
1762 }
1763 } else {
1764 fp1 = floatx80_mul(fp0, float64_to_floatx80(
1765 make_float64(0x3FE45F306DC9C883), status),
1766 status); /* X*2/PI */
1767
1768 n = floatx80_to_int32(fp1, status);
1769 j = 32 + n;
1770
1771 fp0 = floatx80_sub(fp0, pi_tbl[j], status); /* X-Y1 */
1772 fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status),
1773 status); /* FP0 IS R = (X-Y1)-Y2 */
1774
1775 sincont:
1776 if ((n + 1) & 1) {
1777 /* COSPOLY */
1778 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1779 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1780 fp2 = float64_to_floatx80(make_float64(0x3D2AC4D0D6011EE3),
1781 status); /* B8 */
1782 fp3 = float64_to_floatx80(make_float64(0xBDA9396F9F45AC19),
1783 status); /* B7 */
1784
1785 xSign = extractFloatx80Sign(fp0); /* X IS S */
1786 xExp = extractFloatx80Exp(fp0);
1787 xSig = extractFloatx80Frac(fp0);
1788
1789 if (((n + 1) >> 1) & 1) {
1790 xSign ^= 1;
1791 posneg1 = make_float32(0xBF800000); /* -1 */
1792 } else {
1793 xSign ^= 0;
1794 posneg1 = make_float32(0x3F800000); /* 1 */
1795 } /* X IS NOW R'= SGN*R */
1796
1797 fp2 = floatx80_mul(fp2, fp1, status); /* TB8 */
1798 fp3 = floatx80_mul(fp3, fp1, status); /* TB7 */
1799 fp2 = floatx80_add(fp2, float64_to_floatx80(
1800 make_float64(0x3E21EED90612C972), status),
1801 status); /* B6+TB8 */
1802 fp3 = floatx80_add(fp3, float64_to_floatx80(
1803 make_float64(0xBE927E4FB79D9FCF), status),
1804 status); /* B5+TB7 */
1805 fp2 = floatx80_mul(fp2, fp1, status); /* T(B6+TB8) */
1806 fp3 = floatx80_mul(fp3, fp1, status); /* T(B5+TB7) */
1807 fp2 = floatx80_add(fp2, float64_to_floatx80(
1808 make_float64(0x3EFA01A01A01D423), status),
1809 status); /* B4+T(B6+TB8) */
1810 fp4 = packFloatx80(1, 0x3FF5, UINT64_C(0xB60B60B60B61D438));
1811 fp3 = floatx80_add(fp3, fp4, status); /* B3+T(B5+TB7) */
1812 fp2 = floatx80_mul(fp2, fp1, status); /* T(B4+T(B6+TB8)) */
1813 fp1 = floatx80_mul(fp1, fp3, status); /* T(B3+T(B5+TB7)) */
1814 fp4 = packFloatx80(0, 0x3FFA, UINT64_C(0xAAAAAAAAAAAAAB5E));
1815 fp2 = floatx80_add(fp2, fp4, status); /* B2+T(B4+T(B6+TB8)) */
1816 fp1 = floatx80_add(fp1, float32_to_floatx80(
1817 make_float32(0xBF000000), status),
1818 status); /* B1+T(B3+T(B5+TB7)) */
1819 fp0 = floatx80_mul(fp0, fp2, status); /* S(B2+T(B4+T(B6+TB8))) */
1820 fp0 = floatx80_add(fp0, fp1, status);
1821 /* [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] */
1822
1823 x = packFloatx80(xSign, xExp, xSig);
1824 fp0 = floatx80_mul(fp0, x, status);
1825
1826 status->float_rounding_mode = user_rnd_mode;
1827 status->floatx80_rounding_precision = user_rnd_prec;
1828
1829 a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status);
1830
1831 float_raise(float_flag_inexact, status);
1832
1833 return a;
1834 } else {
1835 /* SINPOLY */
1836 xSign = extractFloatx80Sign(fp0); /* X IS R */
1837 xExp = extractFloatx80Exp(fp0);
1838 xSig = extractFloatx80Frac(fp0);
1839
1840 xSign ^= ((n + 1) >> 1) & 1; /* X IS NOW R'= SGN*R */
1841
1842 fp0 = floatx80_mul(fp0, fp0, status); /* FP0 IS S */
1843 fp1 = floatx80_mul(fp0, fp0, status); /* FP1 IS T */
1844 fp3 = float64_to_floatx80(make_float64(0xBD6AAA77CCC994F5),
1845 status); /* A7 */
1846 fp2 = float64_to_floatx80(make_float64(0x3DE612097AAE8DA1),
1847 status); /* A6 */
1848 fp3 = floatx80_mul(fp3, fp1, status); /* T*A7 */
1849 fp2 = floatx80_mul(fp2, fp1, status); /* T*A6 */
1850 fp3 = floatx80_add(fp3, float64_to_floatx80(
1851 make_float64(0xBE5AE6452A118AE4), status),
1852 status); /* A5+T*A7 */
1853 fp2 = floatx80_add(fp2, float64_to_floatx80(
1854 make_float64(0x3EC71DE3A5341531), status),
1855 status); /* A4+T*A6 */
1856 fp3 = floatx80_mul(fp3, fp1, status); /* T(A5+TA7) */
1857 fp2 = floatx80_mul(fp2, fp1, status); /* T(A4+TA6) */
1858 fp3 = floatx80_add(fp3, float64_to_floatx80(
1859 make_float64(0xBF2A01A01A018B59), status),
1860 status); /* A3+T(A5+TA7) */
1861 fp4 = packFloatx80(0, 0x3FF8, UINT64_C(0x88888888888859AF));
1862 fp2 = floatx80_add(fp2, fp4, status); /* A2+T(A4+TA6) */
1863 fp1 = floatx80_mul(fp1, fp3, status); /* T(A3+T(A5+TA7)) */
1864 fp2 = floatx80_mul(fp2, fp0, status); /* S(A2+T(A4+TA6)) */
1865 fp4 = packFloatx80(1, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAA99));
1866 fp1 = floatx80_add(fp1, fp4, status); /* A1+T(A3+T(A5+TA7)) */
1867 fp1 = floatx80_add(fp1, fp2, status);
1868 /* [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] */
1869
1870 x = packFloatx80(xSign, xExp, xSig);
1871 fp0 = floatx80_mul(fp0, x, status); /* R'*S */
1872 fp0 = floatx80_mul(fp0, fp1, status); /* SIN(R')-R' */
1873
1874 status->float_rounding_mode = user_rnd_mode;
1875 status->floatx80_rounding_precision = user_rnd_prec;
1876
1877 a = floatx80_add(fp0, x, status);
1878
1879 float_raise(float_flag_inexact, status);
1880
1881 return a;
1882 }
1883 }
1884 }
1885
1886 /*
1887 * Arc tangent
1888 */
1889
1890 floatx80 floatx80_atan(floatx80 a, float_status *status)
1891 {
1892 bool aSign;
1893 int32_t aExp;
1894 uint64_t aSig;
1895
1896 int8_t user_rnd_mode, user_rnd_prec;
1897
1898 int32_t compact, tbl_index;
1899 floatx80 fp0, fp1, fp2, fp3, xsave;
1900
1901 aSig = extractFloatx80Frac(a);
1902 aExp = extractFloatx80Exp(a);
1903 aSign = extractFloatx80Sign(a);
1904
1905 if (aExp == 0x7FFF) {
1906 if ((uint64_t) (aSig << 1)) {
1907 return propagateFloatx80NaNOneArg(a, status);
1908 }
1909 a = packFloatx80(aSign, piby2_exp, pi_sig);
1910 float_raise(float_flag_inexact, status);
1911 return floatx80_move(a, status);
1912 }
1913
1914 if (aExp == 0 && aSig == 0) {
1915 return packFloatx80(aSign, 0, 0);
1916 }
1917
1918 compact = floatx80_make_compact(aExp, aSig);
1919
1920 user_rnd_mode = status->float_rounding_mode;
1921 user_rnd_prec = status->floatx80_rounding_precision;
1922 status->float_rounding_mode = float_round_nearest_even;
1923 status->floatx80_rounding_precision = 80;
1924
1925 if (compact < 0x3FFB8000 || compact > 0x4002FFFF) {
1926 /* |X| >= 16 or |X| < 1/16 */
1927 if (compact > 0x3FFF8000) { /* |X| >= 16 */
1928 if (compact > 0x40638000) { /* |X| > 2^(100) */
1929 fp0 = packFloatx80(aSign, piby2_exp, pi_sig);
1930 fp1 = packFloatx80(aSign, 0x0001, one_sig);
1931
1932 status->float_rounding_mode = user_rnd_mode;
1933 status->floatx80_rounding_precision = user_rnd_prec;
1934
1935 a = floatx80_sub(fp0, fp1, status);
1936
1937 float_raise(float_flag_inexact, status);
1938
1939 return a;
1940 } else {
1941 fp0 = a;
1942 fp1 = packFloatx80(1, one_exp, one_sig); /* -1 */
1943 fp1 = floatx80_div(fp1, fp0, status); /* X' = -1/X */
1944 xsave = fp1;
1945 fp0 = floatx80_mul(fp1, fp1, status); /* Y = X'*X' */
1946 fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */
1947 fp3 = float64_to_floatx80(make_float64(0xBFB70BF398539E6A),
1948 status); /* C5 */
1949 fp2 = float64_to_floatx80(make_float64(0x3FBC7187962D1D7D),
1950 status); /* C4 */
1951 fp3 = floatx80_mul(fp3, fp1, status); /* Z*C5 */
1952 fp2 = floatx80_mul(fp2, fp1, status); /* Z*C4 */
1953 fp3 = floatx80_add(fp3, float64_to_floatx80(
1954 make_float64(0xBFC24924827107B8), status),
1955 status); /* C3+Z*C5 */
1956 fp2 = floatx80_add(fp2, float64_to_floatx80(
1957 make_float64(0x3FC999999996263E), status),
1958 status); /* C2+Z*C4 */
1959 fp1 = floatx80_mul(fp1, fp3, status); /* Z*(C3+Z*C5) */
1960 fp2 = floatx80_mul(fp2, fp0, status); /* Y*(C2+Z*C4) */
1961 fp1 = floatx80_add(fp1, float64_to_floatx80(
1962 make_float64(0xBFD5555555555536), status),
1963 status); /* C1+Z*(C3+Z*C5) */
1964 fp0 = floatx80_mul(fp0, xsave, status); /* X'*Y */
1965 /* [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] */
1966 fp1 = floatx80_add(fp1, fp2, status);
1967 /* X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ?? */
1968 fp0 = floatx80_mul(fp0, fp1, status);
1969 fp0 = floatx80_add(fp0, xsave, status);
1970 fp1 = packFloatx80(aSign, piby2_exp, pi_sig);
1971
1972 status->float_rounding_mode = user_rnd_mode;
1973 status->floatx80_rounding_precision = user_rnd_prec;
1974
1975 a = floatx80_add(fp0, fp1, status);
1976
1977 float_raise(float_flag_inexact, status);
1978
1979 return a;
1980 }
1981 } else { /* |X| < 1/16 */
1982 if (compact < 0x3FD78000) { /* |X| < 2^(-40) */
1983 status->float_rounding_mode = user_rnd_mode;
1984 status->floatx80_rounding_precision = user_rnd_prec;
1985
1986 a = floatx80_move(a, status);
1987
1988 float_raise(float_flag_inexact, status);
1989
1990 return a;
1991 } else {
1992 fp0 = a;
1993 xsave = a;
1994 fp0 = floatx80_mul(fp0, fp0, status); /* Y = X*X */
1995 fp1 = floatx80_mul(fp0, fp0, status); /* Z = Y*Y */
1996 fp2 = float64_to_floatx80(make_float64(0x3FB344447F876989),
1997 status); /* B6 */
1998 fp3 = float64_to_floatx80(make_float64(0xBFB744EE7FAF45DB),
1999 status); /* B5 */
2000 fp2 = floatx80_mul(fp2, fp1, status); /* Z*B6 */
2001 fp3 = floatx80_mul(fp3, fp1, status); /* Z*B5 */
2002 fp2 = floatx80_add(fp2, float64_to_floatx80(
2003 make_float64(0x3FBC71C646940220), status),
2004 status); /* B4+Z*B6 */
2005 fp3 = floatx80_add(fp3, float64_to_floatx80(
2006 make_float64(0xBFC24924921872F9),
2007 status), status); /* B3+Z*B5 */
2008 fp2 = floatx80_mul(fp2, fp1, status); /* Z*(B4+Z*B6) */
2009 fp1 = floatx80_mul(fp1, fp3, status); /* Z*(B3+Z*B5) */
2010 fp2 = floatx80_add(fp2, float64_to_floatx80(
2011 make_float64(0x3FC9999999998FA9), status),
2012 status); /* B2+Z*(B4+Z*B6) */
2013 fp1 = floatx80_add(fp1, float64_to_floatx80(
2014 make_float64(0xBFD5555555555555), status),
2015 status); /* B1+Z*(B3+Z*B5) */
2016 fp2 = floatx80_mul(fp2, fp0, status); /* Y*(B2+Z*(B4+Z*B6)) */
2017 fp0 = floatx80_mul(fp0, xsave, status); /* X*Y */
2018 /* [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] */
2019 fp1 = floatx80_add(fp1, fp2, status);
2020 /* X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) */
2021 fp0 = floatx80_mul(fp0, fp1, status);
2022
2023 status->float_rounding_mode = user_rnd_mode;
2024 status->floatx80_rounding_precision = user_rnd_prec;
2025
2026 a = floatx80_add(fp0, xsave, status);
2027
2028 float_raise(float_flag_inexact, status);
2029
2030 return a;
2031 }
2032 }
2033 } else {
2034 aSig &= UINT64_C(0xF800000000000000);
2035 aSig |= UINT64_C(0x0400000000000000);
2036 xsave = packFloatx80(aSign, aExp, aSig); /* F */
2037 fp0 = a;
2038 fp1 = a; /* X */
2039 fp2 = packFloatx80(0, one_exp, one_sig); /* 1 */
2040 fp1 = floatx80_mul(fp1, xsave, status); /* X*F */
2041 fp0 = floatx80_sub(fp0, xsave, status); /* X-F */
2042 fp1 = floatx80_add(fp1, fp2, status); /* 1 + X*F */
2043 fp0 = floatx80_div(fp0, fp1, status); /* U = (X-F)/(1+X*F) */
2044
2045 tbl_index = compact;
2046
2047 tbl_index &= 0x7FFF0000;
2048 tbl_index -= 0x3FFB0000;
2049 tbl_index >>= 1;
2050 tbl_index += compact & 0x00007800;
2051 tbl_index >>= 11;
2052
2053 fp3 = atan_tbl[tbl_index];
2054
2055 fp3.high |= aSign ? 0x8000 : 0; /* ATAN(F) */
2056
2057 fp1 = floatx80_mul(fp0, fp0, status); /* V = U*U */
2058 fp2 = float64_to_floatx80(make_float64(0xBFF6687E314987D8),
2059 status); /* A3 */
2060 fp2 = floatx80_add(fp2, fp1, status); /* A3+V */
2061 fp2 = floatx80_mul(fp2, fp1, status); /* V*(A3+V) */
2062 fp1 = floatx80_mul(fp1, fp0, status); /* U*V */
2063 fp2 = floatx80_add(fp2, float64_to_floatx80(
2064 make_float64(0x4002AC6934A26DB3), status),
2065 status); /* A2+V*(A3+V) */
2066 fp1 = floatx80_mul(fp1, float64_to_floatx80(
2067 make_float64(0xBFC2476F4E1DA28E), status),
2068 status); /* A1+U*V */
2069 fp1 = floatx80_mul(fp1, fp2, status); /* A1*U*V*(A2+V*(A3+V)) */
2070 fp0 = floatx80_add(fp0, fp1, status); /* ATAN(U) */
2071
2072 status->float_rounding_mode = user_rnd_mode;
2073 status->floatx80_rounding_precision = user_rnd_prec;
2074
2075 a = floatx80_add(fp0, fp3, status); /* ATAN(X) */
2076
2077 float_raise(float_flag_inexact, status);
2078
2079 return a;
2080 }
2081 }
2082
2083 /*
2084 * Arc sine
2085 */
2086
2087 floatx80 floatx80_asin(floatx80 a, float_status *status)
2088 {
2089 bool aSign;
2090 int32_t aExp;
2091 uint64_t aSig;
2092
2093 int8_t user_rnd_mode, user_rnd_prec;
2094
2095 int32_t compact;
2096 floatx80 fp0, fp1, fp2, one;
2097
2098 aSig = extractFloatx80Frac(a);
2099 aExp = extractFloatx80Exp(a);
2100 aSign = extractFloatx80Sign(a);
2101
2102 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2103 return propagateFloatx80NaNOneArg(a, status);
2104 }
2105
2106 if (aExp == 0 && aSig == 0) {
2107 return packFloatx80(aSign, 0, 0);
2108 }
2109
2110 compact = floatx80_make_compact(aExp, aSig);
2111
2112 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2113 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2114 float_raise(float_flag_inexact, status);
2115 a = packFloatx80(aSign, piby2_exp, pi_sig);
2116 return floatx80_move(a, status);
2117 } else { /* |X| > 1 */
2118 float_raise(float_flag_invalid, status);
2119 return floatx80_default_nan(status);
2120 }
2121
2122 } /* |X| < 1 */
2123
2124 user_rnd_mode = status->float_rounding_mode;
2125 user_rnd_prec = status->floatx80_rounding_precision;
2126 status->float_rounding_mode = float_round_nearest_even;
2127 status->floatx80_rounding_precision = 80;
2128
2129 one = packFloatx80(0, one_exp, one_sig);
2130 fp0 = a;
2131
2132 fp1 = floatx80_sub(one, fp0, status); /* 1 - X */
2133 fp2 = floatx80_add(one, fp0, status); /* 1 + X */
2134 fp1 = floatx80_mul(fp2, fp1, status); /* (1+X)*(1-X) */
2135 fp1 = floatx80_sqrt(fp1, status); /* SQRT((1+X)*(1-X)) */
2136 fp0 = floatx80_div(fp0, fp1, status); /* X/SQRT((1+X)*(1-X)) */
2137
2138 status->float_rounding_mode = user_rnd_mode;
2139 status->floatx80_rounding_precision = user_rnd_prec;
2140
2141 a = floatx80_atan(fp0, status); /* ATAN(X/SQRT((1+X)*(1-X))) */
2142
2143 float_raise(float_flag_inexact, status);
2144
2145 return a;
2146 }
2147
2148 /*
2149 * Arc cosine
2150 */
2151
2152 floatx80 floatx80_acos(floatx80 a, float_status *status)
2153 {
2154 bool aSign;
2155 int32_t aExp;
2156 uint64_t aSig;
2157
2158 int8_t user_rnd_mode, user_rnd_prec;
2159
2160 int32_t compact;
2161 floatx80 fp0, fp1, one;
2162
2163 aSig = extractFloatx80Frac(a);
2164 aExp = extractFloatx80Exp(a);
2165 aSign = extractFloatx80Sign(a);
2166
2167 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2168 return propagateFloatx80NaNOneArg(a, status);
2169 }
2170 if (aExp == 0 && aSig == 0) {
2171 float_raise(float_flag_inexact, status);
2172 return roundAndPackFloatx80(status->floatx80_rounding_precision, 0,
2173 piby2_exp, pi_sig, 0, status);
2174 }
2175
2176 compact = floatx80_make_compact(aExp, aSig);
2177
2178 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2179 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2180 if (aSign) { /* X == -1 */
2181 a = packFloatx80(0, pi_exp, pi_sig);
2182 float_raise(float_flag_inexact, status);
2183 return floatx80_move(a, status);
2184 } else { /* X == +1 */
2185 return packFloatx80(0, 0, 0);
2186 }
2187 } else { /* |X| > 1 */
2188 float_raise(float_flag_invalid, status);
2189 return floatx80_default_nan(status);
2190 }
2191 } /* |X| < 1 */
2192
2193 user_rnd_mode = status->float_rounding_mode;
2194 user_rnd_prec = status->floatx80_rounding_precision;
2195 status->float_rounding_mode = float_round_nearest_even;
2196 status->floatx80_rounding_precision = 80;
2197
2198 one = packFloatx80(0, one_exp, one_sig);
2199 fp0 = a;
2200
2201 fp1 = floatx80_add(one, fp0, status); /* 1 + X */
2202 fp0 = floatx80_sub(one, fp0, status); /* 1 - X */
2203 fp0 = floatx80_div(fp0, fp1, status); /* (1-X)/(1+X) */
2204 fp0 = floatx80_sqrt(fp0, status); /* SQRT((1-X)/(1+X)) */
2205 fp0 = floatx80_atan(fp0, status); /* ATAN(SQRT((1-X)/(1+X))) */
2206
2207 status->float_rounding_mode = user_rnd_mode;
2208 status->floatx80_rounding_precision = user_rnd_prec;
2209
2210 a = floatx80_add(fp0, fp0, status); /* 2 * ATAN(SQRT((1-X)/(1+X))) */
2211
2212 float_raise(float_flag_inexact, status);
2213
2214 return a;
2215 }
2216
2217 /*
2218 * Hyperbolic arc tangent
2219 */
2220
2221 floatx80 floatx80_atanh(floatx80 a, float_status *status)
2222 {
2223 bool aSign;
2224 int32_t aExp;
2225 uint64_t aSig;
2226
2227 int8_t user_rnd_mode, user_rnd_prec;
2228
2229 int32_t compact;
2230 floatx80 fp0, fp1, fp2, one;
2231
2232 aSig = extractFloatx80Frac(a);
2233 aExp = extractFloatx80Exp(a);
2234 aSign = extractFloatx80Sign(a);
2235
2236 if (aExp == 0x7FFF && (uint64_t) (aSig << 1)) {
2237 return propagateFloatx80NaNOneArg(a, status);
2238 }
2239
2240 if (aExp == 0 && aSig == 0) {
2241 return packFloatx80(aSign, 0, 0);
2242 }
2243
2244 compact = floatx80_make_compact(aExp, aSig);
2245
2246 if (compact >= 0x3FFF8000) { /* |X| >= 1 */
2247 if (aExp == one_exp && aSig == one_sig) { /* |X| == 1 */
2248 float_raise(float_flag_divbyzero, status);
2249 return packFloatx80(aSign, floatx80_infinity.high,
2250 floatx80_infinity.low);
2251 } else { /* |X| > 1 */
2252 float_raise(float_flag_invalid, status);
2253 return floatx80_default_nan(status);
2254 }
2255 } /* |X| < 1 */
2256
2257 user_rnd_mode = status->float_rounding_mode;
2258 user_rnd_prec = status->floatx80_rounding_precision;
2259 status->float_rounding_mode = float_round_nearest_even;
2260 status->floatx80_rounding_precision = 80;
2261
2262 one = packFloatx80(0, one_exp, one_sig);
2263 fp2 = packFloatx80(aSign, 0x3FFE, one_sig); /* SIGN(X) * (1/2) */
2264 fp0 = packFloatx80(0, aExp, aSig); /* Y = |X| */
2265 fp1 = packFloatx80(1, aExp, aSig); /* -Y */
2266 fp0 = floatx80_add(fp0, fp0, status); /* 2Y */
2267 fp1 = floatx80_add(fp1, one, status); /* 1-Y */
2268 fp0 = floatx80_div(fp0, fp1, status); /* Z = 2Y/(1-Y) */
2269 fp0 = floatx80_lognp1(fp0, status); /* LOG1P(Z) */
2270
2271 status->float_rounding_mode = user_rnd_mode;
2272 status->floatx80_rounding_precision = user_rnd_prec;
2273
2274 a = floatx80_mul(fp0, fp2,
2275 status); /* ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z) */
2276
2277 float_raise(float_flag_inexact, status);
2278
2279 return a;
2280 }
2281
2282 /*
2283 * e to x minus 1
2284 */
2285
2286 floatx80 floatx80_etoxm1(floatx80 a, float_status *status)
2287 {
2288 bool aSign;
2289 int32_t aExp;
2290 uint64_t aSig;
2291
2292 int8_t user_rnd_mode, user_rnd_prec;
2293
2294 int32_t compact, n, j, m, m1;
2295 floatx80 fp0, fp1, fp2, fp3, l2, sc, onebysc;
2296
2297 aSig = extractFloatx80Frac(a);
2298 aExp = extractFloatx80Exp(a);
2299 aSign = extractFloatx80Sign(a);
2300
2301 if (aExp == 0x7FFF) {
2302 if ((uint64_t) (aSig << 1)) {
2303 return propagateFloatx80NaNOneArg(a, status);
2304 }
2305 if (aSign) {
2306 return packFloatx80(aSign, one_exp, one_sig);
2307 }
2308 return packFloatx80(0, floatx80_infinity.high,
2309 floatx80_infinity.low);
2310 }
2311
2312 if (aExp == 0 && aSig == 0) {
2313 return packFloatx80(aSign, 0, 0);
2314 }
2315
2316 user_rnd_mode = status->float_rounding_mode;
2317 user_rnd_prec = status->floatx80_rounding_precision;
2318 status->float_rounding_mode = float_round_nearest_even;
2319 status->floatx80_rounding_precision = 80;
2320
2321 if (aExp >= 0x3FFD) { /* |X| >= 1/4 */
2322 compact = floatx80_make_compact(aExp, aSig);
2323
2324 if (compact <= 0x4004C215) { /* |X| <= 70 log2 */
2325 fp0 = a;
2326 fp1 = a;
2327 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2328 make_float32(0x42B8AA3B), status),
2329 status); /* 64/log2 * X */
2330 n = floatx80_to_int32(fp0, status); /* int(64/log2*X) */
2331 fp0 = int32_to_floatx80(n, status);
2332
2333 j = n & 0x3F; /* J = N mod 64 */
2334 m = n / 64; /* NOTE: this is really arithmetic right shift by 6 */
2335 if (n < 0 && j) {
2336 /*
2337 * arithmetic right shift is division and
2338 * round towards minus infinity
2339 */
2340 m--;
2341 }
2342 m1 = -m;
2343 /*m += 0x3FFF; // biased exponent of 2^(M) */
2344 /*m1 += 0x3FFF; // biased exponent of -2^(-M) */
2345
2346 fp2 = fp0; /* N */
2347 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2348 make_float32(0xBC317218), status),
2349 status); /* N * L1, L1 = lead(-log2/64) */
2350 l2 = packFloatx80(0, 0x3FDC, UINT64_C(0x82E308654361C4C6));
2351 fp2 = floatx80_mul(fp2, l2, status); /* N * L2, L1+L2 = -log2/64 */
2352 fp0 = floatx80_add(fp0, fp1, status); /* X + N*L1 */
2353 fp0 = floatx80_add(fp0, fp2, status); /* R */
2354
2355 fp1 = floatx80_mul(fp0, fp0, status); /* S = R*R */
2356 fp2 = float32_to_floatx80(make_float32(0x3950097B),
2357 status); /* A6 */
2358 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 is S*A6 */
2359 fp3 = floatx80_mul(float32_to_floatx80(make_float32(0x3AB60B6A),
2360 status), fp1, status); /* fp3 is S*A5 */
2361 fp2 = floatx80_add(fp2, float64_to_floatx80(
2362 make_float64(0x3F81111111174385), status),
2363 status); /* fp2 IS A4+S*A6 */
2364 fp3 = floatx80_add(fp3, float64_to_floatx80(
2365 make_float64(0x3FA5555555554F5A), status),
2366 status); /* fp3 is A3+S*A5 */
2367 fp2 = floatx80_mul(fp2, fp1, status); /* fp2 IS S*(A4+S*A6) */
2368 fp3 = floatx80_mul(fp3, fp1, status); /* fp3 IS S*(A3+S*A5) */
2369 fp2 = floatx80_add(fp2, float64_to_floatx80(
2370 make_float64(0x3FC5555555555555), status),
2371 status); /* fp2 IS A2+S*(A4+S*A6) */
2372 fp3 = floatx80_add(fp3, float32_to_floatx80(
2373 make_float32(0x3F000000), status),
2374 status); /* fp3 IS A1+S*(A3+S*A5) */
2375 fp2 = floatx80_mul(fp2, fp1,
2376 status); /* fp2 IS S*(A2+S*(A4+S*A6)) */
2377 fp1 = floatx80_mul(fp1, fp3,
2378 status); /* fp1 IS S*(A1+S*(A3+S*A5)) */
2379 fp2 = floatx80_mul(fp2, fp0,
2380 status); /* fp2 IS R*S*(A2+S*(A4+S*A6)) */
2381 fp0 = floatx80_add(fp0, fp1,
2382 status); /* fp0 IS R+S*(A1+S*(A3+S*A5)) */
2383 fp0 = floatx80_add(fp0, fp2, status); /* fp0 IS EXP(R) - 1 */
2384
2385 fp0 = floatx80_mul(fp0, exp_tbl[j],
2386 status); /* 2^(J/64)*(Exp(R)-1) */
2387
2388 if (m >= 64) {
2389 fp1 = float32_to_floatx80(exp_tbl2[j], status);
2390 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2391 fp1 = floatx80_add(fp1, onebysc, status);
2392 fp0 = floatx80_add(fp0, fp1, status);
2393 fp0 = floatx80_add(fp0, exp_tbl[j], status);
2394 } else if (m < -3) {
2395 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j],
2396 status), status);
2397 fp0 = floatx80_add(fp0, exp_tbl[j], status);
2398 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2399 fp0 = floatx80_add(fp0, onebysc, status);
2400 } else { /* -3 <= m <= 63 */
2401 fp1 = exp_tbl[j];
2402 fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j],
2403 status), status);
2404 onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); /* -2^(-M) */
2405 fp1 = floatx80_add(fp1, onebysc, status);
2406 fp0 = floatx80_add(fp0, fp1, status);
2407 }
2408
2409 sc = packFloatx80(0, m + 0x3FFF, one_sig);
2410
2411 status->float_rounding_mode = user_rnd_mode;
2412 status->floatx80_rounding_precision = user_rnd_prec;
2413
2414 a = floatx80_mul(fp0, sc, status);
2415
2416 float_raise(float_flag_inexact, status);
2417
2418 return a;
2419 } else { /* |X| > 70 log2 */
2420 if (aSign) {
2421 fp0 = float32_to_floatx80(make_float32(0xBF800000),
2422 status); /* -1 */
2423
2424 status->float_rounding_mode = user_rnd_mode;
2425 status->floatx80_rounding_precision = user_rnd_prec;
2426
2427 a = floatx80_add(fp0, float32_to_floatx80(
2428 make_float32(0x00800000), status),
2429 status); /* -1 + 2^(-126) */
2430
2431 float_raise(float_flag_inexact, status);
2432
2433 return a;
2434 } else {
2435 status->float_rounding_mode = user_rnd_mode;
2436 status->floatx80_rounding_precision = user_rnd_prec;
2437
2438 return floatx80_etox(a, status);
2439 }
2440 }
2441 } else { /* |X| < 1/4 */
2442 if (aExp >= 0x3FBE) {
2443 fp0 = a;
2444 fp0 = floatx80_mul(fp0, fp0, status); /* S = X*X */
2445 fp1 = float32_to_floatx80(make_float32(0x2F30CAA8),
2446 status); /* B12 */
2447 fp1 = floatx80_mul(fp1, fp0, status); /* S * B12 */
2448 fp2 = float32_to_floatx80(make_float32(0x310F8290),
2449 status); /* B11 */
2450 fp1 = floatx80_add(fp1, float32_to_floatx80(
2451 make_float32(0x32D73220), status),
2452 status); /* B10 */
2453 fp2 = floatx80_mul(fp2, fp0, status);
2454 fp1 = floatx80_mul(fp1, fp0, status);
2455 fp2 = floatx80_add(fp2, float32_to_floatx80(
2456 make_float32(0x3493F281), status),
2457 status); /* B9 */
2458 fp1 = floatx80_add(fp1, float64_to_floatx80(
2459 make_float64(0x3EC71DE3A5774682), status),
2460 status); /* B8 */
2461 fp2 = floatx80_mul(fp2, fp0, status);
2462 fp1 = floatx80_mul(fp1, fp0, status);
2463 fp2 = floatx80_add(fp2, float64_to_floatx80(
2464 make_float64(0x3EFA01A019D7CB68), status),
2465 status); /* B7 */
2466 fp1 = floatx80_add(fp1, float64_to_floatx80(
2467 make_float64(0x3F2A01A01A019DF3), status),
2468 status); /* B6 */
2469 fp2 = floatx80_mul(fp2, fp0, status);
2470 fp1 = floatx80_mul(fp1, fp0, status);
2471 fp2 = floatx80_add(fp2, float64_to_floatx80(
2472 make_float64(0x3F56C16C16C170E2), status),
2473 status); /* B5 */
2474 fp1 = floatx80_add(fp1, float64_to_floatx80(
2475 make_float64(0x3F81111111111111), status),
2476 status); /* B4 */
2477 fp2 = floatx80_mul(fp2, fp0, status);
2478 fp1 = floatx80_mul(fp1, fp0, status);
2479 fp2 = floatx80_add(fp2, float64_to_floatx80(
2480 make_float64(0x3FA5555555555555), status),
2481 status); /* B3 */
2482 fp3 = packFloatx80(0, 0x3FFC, UINT64_C(0xAAAAAAAAAAAAAAAB));
2483 fp1 = floatx80_add(fp1, fp3, status); /* B2 */
2484 fp2 = floatx80_mul(fp2, fp0, status);
2485 fp1 = floatx80_mul(fp1, fp0, status);
2486
2487 fp2 = floatx80_mul(fp2, fp0, status);
2488 fp1 = floatx80_mul(fp1, a, status);
2489
2490 fp0 = floatx80_mul(fp0, float32_to_floatx80(
2491 make_float32(0x3F000000), status),
2492 status); /* S*B1 */
2493 fp1 = floatx80_add(fp1, fp2, status); /* Q */
2494 fp0 = floatx80_add(fp0, fp1, status); /* S*B1+Q */
2495
2496 status->float_rounding_mode = user_rnd_mode;
2497 status->floatx80_rounding_precision = user_rnd_prec;
2498
2499 a = floatx80_add(fp0, a, status);
2500
2501 float_raise(float_flag_inexact, status);
2502
2503 return a;
2504 } else { /* |X| < 2^(-65) */
2505 sc = packFloatx80(1, 1, one_sig);
2506 fp0 = a;
2507
2508 if (aExp < 0x0033) { /* |X| < 2^(-16382) */
2509 fp0 = floatx80_mul(fp0, float64_to_floatx80(
2510 make_float64(0x48B0000000000000), status),
2511 status);
2512 fp0 = floatx80_add(fp0, sc, status);
2513
2514 status->float_rounding_mode = user_rnd_mode;
2515 status->floatx80_rounding_precision = user_rnd_prec;
2516
2517 a = floatx80_mul(fp0, float64_to_floatx80(
2518 make_float64(0x3730000000000000), status),
2519 status);
2520 } else {
2521 status->float_rounding_mode = user_rnd_mode;
2522 status->floatx80_rounding_precision = user_rnd_prec;
2523
2524 a = floatx80_add(fp0, sc, status);
2525 }
2526
2527 float_raise(float_flag_inexact, status);
2528
2529 return a;
2530 }
2531 }
2532 }
2533
2534 /*
2535 * Hyperbolic tangent
2536 */
2537
2538 floatx80 floatx80_tanh(floatx80 a, float_status *status)
2539 {
2540 bool aSign, vSign;
2541 int32_t aExp, vExp;
2542 uint64_t aSig, vSig;
2543
2544 int8_t user_rnd_mode, user_rnd_prec;
2545
2546 int32_t compact;
2547 floatx80 fp0, fp1;
2548 uint32_t sign;
2549
2550 aSig = extractFloatx80Frac(a);
2551 aExp = extractFloatx80Exp(a);
2552 aSign = extractFloatx80Sign(a);
2553
2554 if (aExp == 0x7FFF) {
2555 if ((uint64_t) (aSig << 1)) {
2556 return propagateFloatx80NaNOneArg(a, status);
2557 }
2558 return packFloatx80(aSign, one_exp, one_sig);
2559 }
2560
2561 if (aExp == 0 && aSig == 0) {
2562 return packFloatx80(aSign, 0, 0);
2563 }
2564
2565 user_rnd_mode = status->float_rounding_mode;
2566 user_rnd_prec = status->floatx80_rounding_precision;
2567 status->float_rounding_mode = float_round_nearest_even;
2568 status->floatx80_rounding_precision = 80;
2569
2570 compact = floatx80_make_compact(aExp, aSig);
2571
2572 if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) {
2573 /* TANHBORS */
2574 if (compact < 0x3FFF8000) {
2575 /* TANHSM */
2576 status->float_rounding_mode = user_rnd_mode;
2577 status->floatx80_rounding_precision = user_rnd_prec;
2578
2579 a = floatx80_move(a, status);
2580
2581 float_raise(float_flag_inexact, status);
2582
2583 return a;
2584 } else {
2585 if (compact > 0x40048AA1) {
2586 /* TANHHUGE */
2587 sign = 0x3F800000;
2588 sign |= aSign ? 0x80000000 : 0x00000000;
2589 fp0 = float32_to_floatx80(make_float32(sign), status);
2590 sign &= 0x80000000;
2591 sign ^= 0x80800000; /* -SIGN(X)*EPS */
2592
2593 status->float_rounding_mode = user_rnd_mode;
2594 status->floatx80_rounding_precision = user_rnd_prec;
2595
2596 a = floatx80_add(fp0, float32_to_floatx80(make_float32(sign),
2597 status), status);
2598
2599 float_raise(float_flag_inexact, status);
2600
2601 return a;
2602 } else {
2603 fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */
2604 fp0 = floatx80_etox(fp0, status); /* FP0 IS EXP(Y) */
2605 fp0 = floatx80_add(fp0, float32_to_floatx80(
2606 make_float32(0x3F800000),
2607 status), status); /* EXP(Y)+1 */
2608 sign = aSign ? 0x80000000 : 0x00000000;
2609 fp1 = floatx80_div(float32_to_floatx80(make_float32(
2610 sign ^ 0xC0000000), status), fp0,
2611 status); /* -SIGN(X)*2 / [EXP(Y)+1] */
2612 fp0 = float32_to_floatx80(make_float32(sign | 0x3F800000),
2613 status); /* SIGN */
2614
2615 status->float_rounding_mode = user_rnd_mode;
2616 status->floatx80_rounding_precision = user_rnd_prec;
2617
2618 a = floatx80_add(fp1, fp0, status);
2619
2620 float_raise(float_flag_inexact, status);
2621
2622 return a;
2623 }
2624 }
2625 } else { /* 2**(-40) < |X| < (5/2)LOG2 */
2626 fp0 = packFloatx80(0, aExp + 1, aSig); /* Y = 2|X| */
2627 fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */
2628 fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x40000000),
2629 status),
2630 status); /* Z+2 */
2631
2632 vSign = extractFloatx80Sign(fp1);
2633 vExp = extractFloatx80Exp(fp1);
2634 vSig = extractFloatx80Frac(fp1);
2635
2636 fp1 = packFloatx80(vSign ^ aSign, vExp, vSig);
2637
2638 status->float_rounding_mode = user_rnd_mode;
2639 status->floatx80_rounding_precision = user_rnd_prec;
2640
2641 a = floatx80_div(fp0, fp1, status);
2642
2643 float_raise(float_flag_inexact, status);
2644
2645 return a;
2646 }
2647 }
2648
2649 /*
2650 * Hyperbolic sine
2651 */
2652
2653 floatx80 floatx80_sinh(floatx80 a, float_status *status)
2654 {
2655 bool aSign;
2656 int32_t aExp;
2657 uint64_t aSig;
2658
2659 int8_t user_rnd_mode, user_rnd_prec;
2660
2661 int32_t compact;
2662 floatx80 fp0, fp1, fp2;
2663 float32 fact;
2664
2665 aSig = extractFloatx80Frac(a);
2666 aExp = extractFloatx80Exp(a);
2667 aSign = extractFloatx80Sign(a);
2668
2669 if (aExp == 0x7FFF) {
2670 if ((uint64_t) (aSig << 1)) {
2671 return propagateFloatx80NaNOneArg(a, status);
2672 }
2673 return packFloatx80(aSign, floatx80_infinity.high,
2674 floatx80_infinity.low);
2675 }
2676
2677 if (aExp == 0 && aSig == 0) {
2678 return packFloatx80(aSign, 0, 0);
2679 }
2680
2681 user_rnd_mode = status->float_rounding_mode;
2682 user_rnd_prec = status->floatx80_rounding_precision;
2683 status->float_rounding_mode = float_round_nearest_even;
2684 status->floatx80_rounding_precision = 80;
2685
2686 compact = floatx80_make_compact(aExp, aSig);
2687
2688 if (compact > 0x400CB167) {
2689 /* SINHBIG */
2690 if (compact > 0x400CB2B3) {
2691 status->float_rounding_mode = user_rnd_mode;
2692 status->floatx80_rounding_precision = user_rnd_prec;
2693
2694 return roundAndPackFloatx80(status->floatx80_rounding_precision,
2695 aSign, 0x8000, aSig, 0, status);
2696 } else {
2697 fp0 = floatx80_abs(a); /* Y = |X| */
2698 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2699 make_float64(0x40C62D38D3D64634), status),
2700 status); /* (|X|-16381LOG2_LEAD) */
2701 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2702 make_float64(0x3D6F90AEB1E75CC7), status),
2703 status); /* |X| - 16381 LOG2, ACCURATE */
2704 fp0 = floatx80_etox(fp0, status);
2705 fp2 = packFloatx80(aSign, 0x7FFB, one_sig);
2706
2707 status->float_rounding_mode = user_rnd_mode;
2708 status->floatx80_rounding_precision = user_rnd_prec;
2709
2710 a = floatx80_mul(fp0, fp2, status);
2711
2712 float_raise(float_flag_inexact, status);
2713
2714 return a;
2715 }
2716 } else { /* |X| < 16380 LOG2 */
2717 fp0 = floatx80_abs(a); /* Y = |X| */
2718 fp0 = floatx80_etoxm1(fp0, status); /* FP0 IS Z = EXPM1(Y) */
2719 fp1 = floatx80_add(fp0, float32_to_floatx80(make_float32(0x3F800000),
2720 status), status); /* 1+Z */
2721 fp2 = fp0;
2722 fp0 = floatx80_div(fp0, fp1, status); /* Z/(1+Z) */
2723 fp0 = floatx80_add(fp0, fp2, status);
2724
2725 fact = packFloat32(aSign, 0x7E, 0);
2726
2727 status->float_rounding_mode = user_rnd_mode;
2728 status->floatx80_rounding_precision = user_rnd_prec;
2729
2730 a = floatx80_mul(fp0, float32_to_floatx80(fact, status), status);
2731
2732 float_raise(float_flag_inexact, status);
2733
2734 return a;
2735 }
2736 }
2737
2738 /*
2739 * Hyperbolic cosine
2740 */
2741
2742 floatx80 floatx80_cosh(floatx80 a, float_status *status)
2743 {
2744 int32_t aExp;
2745 uint64_t aSig;
2746
2747 int8_t user_rnd_mode, user_rnd_prec;
2748
2749 int32_t compact;
2750 floatx80 fp0, fp1;
2751
2752 aSig = extractFloatx80Frac(a);
2753 aExp = extractFloatx80Exp(a);
2754
2755 if (aExp == 0x7FFF) {
2756 if ((uint64_t) (aSig << 1)) {
2757 return propagateFloatx80NaNOneArg(a, status);
2758 }
2759 return packFloatx80(0, floatx80_infinity.high,
2760 floatx80_infinity.low);
2761 }
2762
2763 if (aExp == 0 && aSig == 0) {
2764 return packFloatx80(0, one_exp, one_sig);
2765 }
2766
2767 user_rnd_mode = status->float_rounding_mode;
2768 user_rnd_prec = status->floatx80_rounding_precision;
2769 status->float_rounding_mode = float_round_nearest_even;
2770 status->floatx80_rounding_precision = 80;
2771
2772 compact = floatx80_make_compact(aExp, aSig);
2773
2774 if (compact > 0x400CB167) {
2775 if (compact > 0x400CB2B3) {
2776 status->float_rounding_mode = user_rnd_mode;
2777 status->floatx80_rounding_precision = user_rnd_prec;
2778 return roundAndPackFloatx80(status->floatx80_rounding_precision, 0,
2779 0x8000, one_sig, 0, status);
2780 } else {
2781 fp0 = packFloatx80(0, aExp, aSig);
2782 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2783 make_float64(0x40C62D38D3D64634), status),
2784 status);
2785 fp0 = floatx80_sub(fp0, float64_to_floatx80(
2786 make_float64(0x3D6F90AEB1E75CC7), status),
2787 status);
2788 fp0 = floatx80_etox(fp0, status);
2789 fp1 = packFloatx80(0, 0x7FFB, one_sig);
2790
2791 status->float_rounding_mode = user_rnd_mode;
2792 status->floatx80_rounding_precision = user_rnd_prec;
2793
2794 a = floatx80_mul(fp0, fp1, status);
2795
2796 float_raise(float_flag_inexact, status);
2797
2798 return a;
2799 }
2800 }
2801
2802 fp0 = packFloatx80(0, aExp, aSig); /* |X| */
2803 fp0 = floatx80_etox(fp0, status); /* EXP(|X|) */
2804 fp0 = floatx80_mul(fp0, float32_to_floatx80(make_float32(0x3F000000),
2805 status), status); /* (1/2)*EXP(|X|) */
2806 fp1 = float32_to_floatx80(make_float32(0x3E800000), status); /* 1/4 */
2807 fp1 = floatx80_div(fp1, fp0, status); /* 1/(2*EXP(|X|)) */
2808
2809 status->float_rounding_mode = user_rnd_mode;
2810 status->floatx80_rounding_precision = user_rnd_prec;
2811
2812 a = floatx80_add(fp0, fp1, status);
2813
2814 float_raise(float_flag_inexact, status);
2815
2816 return a;
2817 }
armbru's QEMU hackery