xref: /openbmc/qemu/target/m68k/softfloat.c (revision 992861fb)
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 }
2818