target/m68k/softfloat.c

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