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