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