![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to softfloat_fpsp.c CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [HATARI the Atari ST Emulator] / hatari / src / cpu / softfloat|
Return to softfloat_fpsp.c CVS log | Up to [HATARI the Atari ST Emulator] / hatari / src / cpu / softfloat |
1.1 ! root 1: ! 2: /*============================================================================ ! 3: ! 4: This C source file is an extension to the SoftFloat IEC/IEEE Floating-point ! 5: Arithmetic Package, Release 2a. ! 6: ! 7: Written by Andreas Grabher for Previous, NeXT Computer Emulator. ! 8: ! 9: =============================================================================*/ ! 10: ! 11: #include <stdint.h> ! 12: #include <stdlib.h> ! 13: ! 14: #include "softfloat.h" ! 15: #include "softfloat-specialize.h" ! 16: #include "softfloat_fpsp_tables.h" ! 17: ! 18: ! 19: /*---------------------------------------------------------------------------- ! 20: | Algorithms for transcendental functions supported by MC68881 and MC68882 ! 21: | mathematical coprocessors. The functions are derived from FPSP library. ! 22: *----------------------------------------------------------------------------*/ ! 23: ! 24: #define pi_sig LIT64(0xc90fdaa22168c235) ! 25: #define pi_sig0 LIT64(0xc90fdaa22168c234) ! 26: #define pi_sig1 LIT64(0xc4c6628b80dc1cd1) ! 27: ! 28: #define pi_exp 0x4000 ! 29: #define piby2_exp 0x3FFF ! 30: #define piby4_exp 0x3FFE ! 31: ! 32: #define one_exp 0x3FFF ! 33: #define one_sig LIT64(0x8000000000000000) ! 34: ! 35: #define SET_PREC \ ! 36: int8_t user_rnd_mode, user_rnd_prec; \ ! 37: user_rnd_mode = status->float_rounding_mode; \ ! 38: user_rnd_prec = status->floatx80_rounding_precision; \ ! 39: status->float_rounding_mode = float_round_nearest_even; \ ! 40: status->floatx80_rounding_precision = 80 ! 41: ! 42: #define RESET_PREC \ ! 43: status->float_rounding_mode = user_rnd_mode; \ ! 44: status->floatx80_rounding_precision = user_rnd_prec ! 45: ! 46: /*---------------------------------------------------------------------------- ! 47: | Function for compactifying extended double-precision floating point values. ! 48: *----------------------------------------------------------------------------*/ ! 49: ! 50: static int32_t floatx80_make_compact(int32_t aExp, uint64_t aSig) ! 51: { ! 52: return (aExp<<16)|(aSig>>48); ! 53: } ! 54: ! 55: ! 56: /*---------------------------------------------------------------------------- ! 57: | Arc cosine ! 58: *----------------------------------------------------------------------------*/ ! 59: ! 60: floatx80 floatx80_acos(floatx80 a, float_status *status) ! 61: { ! 62: flag aSign; ! 63: int32_t aExp; ! 64: uint64_t aSig; ! 65: ! 66: ! 67: int32_t compact; ! 68: floatx80 fp0, fp1, one; ! 69: ! 70: aSig = extractFloatx80Frac(a); ! 71: aExp = extractFloatx80Exp(a); ! 72: aSign = extractFloatx80Sign(a); ! 73: ! 74: if (aExp == 0x7FFF && (uint64_t) (aSig<<1)) { ! 75: return propagateFloatx80NaNOneArg(a, status); ! 76: } ! 77: if (aExp == 0 && aSig == 0) { ! 78: float_raise(float_flag_inexact, status); ! 79: return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, piby2_exp, pi_sig, 0, status); ! 80: } ! 81: ! 82: compact = floatx80_make_compact(aExp, aSig); ! 83: ! 84: if (compact >= 0x3FFF8000) { // |X| >= 1 ! 85: if (aExp == one_exp && aSig == one_sig) { // |X| == 1 ! 86: if (aSign) { // X == -1 ! 87: a = packFloatx80(0, pi_exp, pi_sig); ! 88: float_raise(float_flag_inexact, status); ! 89: return floatx80_move(a, status); ! 90: } else { // X == +1 ! 91: return packFloatx80(0, 0, 0); ! 92: } ! 93: } else { // |X| > 1 ! 94: float_raise(float_flag_invalid, status); ! 95: a.low = floatx80_default_nan_low; ! 96: a.high = floatx80_default_nan_high; ! 97: return a; ! 98: } ! 99: } // |X| < 1 ! 100: ! 101: SET_PREC; ! 102: ! 103: one = packFloatx80(0, one_exp, one_sig); ! 104: fp0 = a; ! 105: ! 106: fp1 = floatx80_add(one, fp0, status); // 1 + X ! 107: fp0 = floatx80_sub(one, fp0, status); // 1 - X ! 108: fp0 = floatx80_div(fp0, fp1, status); // (1-X)/(1+X) ! 109: fp0 = floatx80_sqrt(fp0, status); // SQRT((1-X)/(1+X)) ! 110: fp0 = floatx80_atan(fp0, status); // ATAN(SQRT((1-X)/(1+X))) ! 111: ! 112: RESET_PREC; ! 113: ! 114: a = floatx80_add(fp0, fp0, status); // 2 * ATAN(SQRT((1-X)/(1+X))) ! 115: ! 116: float_raise(float_flag_inexact, status); ! 117: ! 118: return a; ! 119: } ! 120: ! 121: /*---------------------------------------------------------------------------- ! 122: | Arc sine ! 123: *----------------------------------------------------------------------------*/ ! 124: ! 125: floatx80 floatx80_asin(floatx80 a, float_status *status) ! 126: { ! 127: flag aSign; ! 128: int32_t aExp; ! 129: uint64_t aSig; ! 130: ! 131: int32_t compact; ! 132: floatx80 fp0, fp1, fp2, one; ! 133: ! 134: aSig = extractFloatx80Frac(a); ! 135: aExp = extractFloatx80Exp(a); ! 136: aSign = extractFloatx80Sign(a); ! 137: ! 138: if (aExp == 0x7FFF && (uint64_t) (aSig<<1)) { ! 139: return propagateFloatx80NaNOneArg(a, status); ! 140: } ! 141: ! 142: if (aExp == 0 && aSig == 0) { ! 143: return packFloatx80(aSign, 0, 0); ! 144: } ! 145: ! 146: compact = floatx80_make_compact(aExp, aSig); ! 147: ! 148: if (compact >= 0x3FFF8000) { // |X| >= 1 ! 149: if (aExp == one_exp && aSig == one_sig) { // |X| == 1 ! 150: float_raise(float_flag_inexact, status); ! 151: a = packFloatx80(aSign, piby2_exp, pi_sig); ! 152: return floatx80_move(a, status); ! 153: } else { // |X| > 1 ! 154: float_raise(float_flag_invalid, status); ! 155: a.low = floatx80_default_nan_low; ! 156: a.high = floatx80_default_nan_high; ! 157: return a; ! 158: } ! 159: ! 160: } // |X| < 1 ! 161: ! 162: SET_PREC; ! 163: ! 164: one = packFloatx80(0, one_exp, one_sig); ! 165: fp0 = a; ! 166: ! 167: fp1 = floatx80_sub(one, fp0, status); // 1 - X ! 168: fp2 = floatx80_add(one, fp0, status); // 1 + X ! 169: fp1 = floatx80_mul(fp2, fp1, status); // (1+X)*(1-X) ! 170: fp1 = floatx80_sqrt(fp1, status); // SQRT((1+X)*(1-X)) ! 171: fp0 = floatx80_div(fp0, fp1, status); // X/SQRT((1+X)*(1-X)) ! 172: ! 173: RESET_PREC; ! 174: ! 175: a = floatx80_atan(fp0, status); // ATAN(X/SQRT((1+X)*(1-X))) ! 176: ! 177: float_raise(float_flag_inexact, status); ! 178: ! 179: return a; ! 180: } ! 181: ! 182: /*---------------------------------------------------------------------------- ! 183: | Arc tangent ! 184: *----------------------------------------------------------------------------*/ ! 185: ! 186: floatx80 floatx80_atan(floatx80 a, float_status *status) ! 187: { ! 188: flag aSign; ! 189: int32_t aExp; ! 190: uint64_t aSig; ! 191: ! 192: int32_t compact, tbl_index; ! 193: floatx80 fp0, fp1, fp2, fp3, xsave; ! 194: ! 195: aSig = extractFloatx80Frac(a); ! 196: aExp = extractFloatx80Exp(a); ! 197: aSign = extractFloatx80Sign(a); ! 198: ! 199: if (aExp == 0x7FFF) { ! 200: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 201: a = packFloatx80(aSign, piby2_exp, pi_sig); ! 202: float_raise(float_flag_inexact, status); ! 203: return floatx80_move(a, status); ! 204: } ! 205: ! 206: if (aExp == 0 && aSig == 0) { ! 207: return packFloatx80(aSign, 0, 0); ! 208: } ! 209: ! 210: compact = floatx80_make_compact(aExp, aSig); ! 211: ! 212: SET_PREC; ! 213: ! 214: if (compact < 0x3FFB8000 || compact > 0x4002FFFF) { // |X| >= 16 or |X| < 1/16 ! 215: if (compact > 0x3FFF8000) { // |X| >= 16 ! 216: if (compact > 0x40638000) { // |X| > 2^(100) ! 217: fp0 = packFloatx80(aSign, piby2_exp, pi_sig); ! 218: fp1 = packFloatx80(aSign, 0x0001, one_sig); ! 219: ! 220: RESET_PREC; ! 221: ! 222: a = floatx80_sub(fp0, fp1, status); ! 223: ! 224: float_raise(float_flag_inexact, status); ! 225: ! 226: return a; ! 227: } else { ! 228: fp0 = a; ! 229: fp1 = packFloatx80(1, one_exp, one_sig); // -1 ! 230: fp1 = floatx80_div(fp1, fp0, status); // X' = -1/X ! 231: xsave = fp1; ! 232: fp0 = floatx80_mul(fp1, fp1, status); // Y = X'*X' ! 233: fp1 = floatx80_mul(fp0, fp0, status); // Z = Y*Y ! 234: fp3 = float64_to_floatx80(LIT64(0xBFB70BF398539E6A), status); // C5 ! 235: fp2 = float64_to_floatx80(LIT64(0x3FBC7187962D1D7D), status); // C4 ! 236: fp3 = floatx80_mul(fp3, fp1, status); // Z*C5 ! 237: fp2 = floatx80_mul(fp2, fp1, status); // Z*C4 ! 238: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBFC24924827107B8), status), status); // C3+Z*C5 ! 239: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FC999999996263E), status), status); // C2+Z*C4 ! 240: fp1 = floatx80_mul(fp1, fp3, status); // Z*(C3+Z*C5) ! 241: fp2 = floatx80_mul(fp2, fp0, status); // Y*(C2+Z*C4) ! 242: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBFD5555555555536), status), status); // C1+Z*(C3+Z*C5) ! 243: fp0 = floatx80_mul(fp0, xsave, status); // X'*Y ! 244: fp1 = floatx80_add(fp1, fp2, status); // [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] ! 245: fp0 = floatx80_mul(fp0, fp1, status); // X'*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ?? ! 246: fp0 = floatx80_add(fp0, xsave, status); ! 247: fp1 = packFloatx80(aSign, piby2_exp, pi_sig); ! 248: ! 249: RESET_PREC; ! 250: ! 251: a = floatx80_add(fp0, fp1, status); ! 252: ! 253: float_raise(float_flag_inexact, status); ! 254: ! 255: return a; ! 256: } ! 257: } else { // |X| < 1/16 ! 258: if (compact < 0x3FD78000) { // |X| < 2^(-40) ! 259: RESET_PREC; ! 260: ! 261: a = floatx80_move(a, status); ! 262: ! 263: float_raise(float_flag_inexact, status); ! 264: ! 265: return a; ! 266: } else { ! 267: fp0 = a; ! 268: xsave = a; ! 269: fp0 = floatx80_mul(fp0, fp0, status); // Y = X*X ! 270: fp1 = floatx80_mul(fp0, fp0, status); // Z = Y*Y ! 271: fp2 = float64_to_floatx80(LIT64(0x3FB344447F876989), status); // B6 ! 272: fp3 = float64_to_floatx80(LIT64(0xBFB744EE7FAF45DB), status); // B5 ! 273: fp2 = floatx80_mul(fp2, fp1, status); // Z*B6 ! 274: fp3 = floatx80_mul(fp3, fp1, status); // Z*B5 ! 275: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FBC71C646940220), status), status); // B4+Z*B6 ! 276: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBFC24924921872F9), status), status); // B3+Z*B5 ! 277: fp2 = floatx80_mul(fp2, fp1, status); // Z*(B4+Z*B6) ! 278: fp1 = floatx80_mul(fp1, fp3, status); // Z*(B3+Z*B5) ! 279: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FC9999999998FA9), status), status); // B2+Z*(B4+Z*B6) ! 280: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBFD5555555555555), status), status); // B1+Z*(B3+Z*B5) ! 281: fp2 = floatx80_mul(fp2, fp0, status); // Y*(B2+Z*(B4+Z*B6)) ! 282: fp0 = floatx80_mul(fp0, xsave, status); // X*Y ! 283: fp1 = floatx80_add(fp1, fp2, status); // [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] ! 284: fp0 = floatx80_mul(fp0, fp1, status); // X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) ! 285: ! 286: RESET_PREC; ! 287: ! 288: a = floatx80_add(fp0, xsave, status); ! 289: ! 290: float_raise(float_flag_inexact, status); ! 291: ! 292: return a; ! 293: } ! 294: } ! 295: } else { ! 296: aSig &= LIT64(0xF800000000000000); ! 297: aSig |= LIT64(0x0400000000000000); ! 298: xsave = packFloatx80(aSign, aExp, aSig); // F ! 299: fp0 = a; ! 300: fp1 = a; // X ! 301: fp2 = packFloatx80(0, one_exp, one_sig); // 1 ! 302: fp1 = floatx80_mul(fp1, xsave, status); // X*F ! 303: fp0 = floatx80_sub(fp0, xsave, status); // X-F ! 304: fp1 = floatx80_add(fp1, fp2, status); // 1 + X*F ! 305: fp0 = floatx80_div(fp0, fp1, status); // U = (X-F)/(1+X*F) ! 306: ! 307: tbl_index = compact; ! 308: ! 309: tbl_index &= 0x7FFF0000; ! 310: tbl_index -= 0x3FFB0000; ! 311: tbl_index >>= 1; ! 312: tbl_index += compact&0x00007800; ! 313: tbl_index >>= 11; ! 314: ! 315: fp3 = atan_tbl[tbl_index]; ! 316: ! 317: fp3.high |= aSign ? 0x8000 : 0; // ATAN(F) ! 318: ! 319: fp1 = floatx80_mul(fp0, fp0, status); // V = U*U ! 320: fp2 = float64_to_floatx80(LIT64(0xBFF6687E314987D8), status); // A3 ! 321: fp2 = floatx80_add(fp2, fp1, status); // A3+V ! 322: fp2 = floatx80_mul(fp2, fp1, status); // V*(A3+V) ! 323: fp1 = floatx80_mul(fp1, fp0, status); // U*V ! 324: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x4002AC6934A26DB3), status), status); // A2+V*(A3+V) ! 325: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0xBFC2476F4E1DA28E), status), status); // A1+U*V ! 326: fp1 = floatx80_mul(fp1, fp2, status); // A1*U*V*(A2+V*(A3+V)) ! 327: fp0 = floatx80_add(fp0, fp1, status); // ATAN(U) ! 328: ! 329: RESET_PREC; ! 330: ! 331: a = floatx80_add(fp0, fp3, status); // ATAN(X) ! 332: ! 333: float_raise(float_flag_inexact, status); ! 334: ! 335: return a; ! 336: } ! 337: } ! 338: ! 339: /*---------------------------------------------------------------------------- ! 340: | Hyperbolic arc tangent ! 341: *----------------------------------------------------------------------------*/ ! 342: ! 343: floatx80 floatx80_atanh(floatx80 a, float_status *status) ! 344: { ! 345: flag aSign; ! 346: int32_t aExp; ! 347: uint64_t aSig; ! 348: ! 349: int32_t compact; ! 350: floatx80 fp0, fp1, fp2, one; ! 351: ! 352: aSig = extractFloatx80Frac(a); ! 353: aExp = extractFloatx80Exp(a); ! 354: aSign = extractFloatx80Sign(a); ! 355: ! 356: if (aExp == 0x7FFF && (uint64_t) (aSig<<1)) { ! 357: return propagateFloatx80NaNOneArg(a, status); ! 358: } ! 359: ! 360: if (aExp == 0 && aSig == 0) { ! 361: return packFloatx80(aSign, 0, 0); ! 362: } ! 363: ! 364: compact = floatx80_make_compact(aExp, aSig); ! 365: ! 366: if (compact >= 0x3FFF8000) { // |X| >= 1 ! 367: if (aExp == one_exp && aSig == one_sig) { // |X| == 1 ! 368: float_raise(float_flag_divbyzero, status); ! 369: return packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low); ! 370: } else { // |X| > 1 ! 371: float_raise(float_flag_invalid, status); ! 372: a.low = floatx80_default_nan_low; ! 373: a.high = floatx80_default_nan_high; ! 374: return a; ! 375: } ! 376: } // |X| < 1 ! 377: ! 378: SET_PREC; ! 379: ! 380: one = packFloatx80(0, one_exp, one_sig); ! 381: fp2 = packFloatx80(aSign, 0x3FFE, one_sig); // SIGN(X) * (1/2) ! 382: fp0 = packFloatx80(0, aExp, aSig); // Y = |X| ! 383: fp1 = packFloatx80(1, aExp, aSig); // -Y ! 384: fp0 = floatx80_add(fp0, fp0, status); // 2Y ! 385: fp1 = floatx80_add(fp1, one, status); // 1-Y ! 386: fp0 = floatx80_div(fp0, fp1, status); // Z = 2Y/(1-Y) ! 387: fp0 = floatx80_lognp1(fp0, status); // LOG1P(Z) ! 388: ! 389: RESET_PREC; ! 390: ! 391: a = floatx80_mul(fp0, fp2, status); // ATANH(X) = SIGN(X) * (1/2) * LOG1P(Z) ! 392: ! 393: float_raise(float_flag_inexact, status); ! 394: ! 395: return a; ! 396: } ! 397: ! 398: /*---------------------------------------------------------------------------- ! 399: | Cosine ! 400: *----------------------------------------------------------------------------*/ ! 401: ! 402: floatx80 floatx80_cos(floatx80 a, float_status *status) ! 403: { ! 404: flag aSign, xSign; ! 405: int32_t aExp, xExp; ! 406: uint64_t aSig, xSig; ! 407: ! 408: int32_t compact, l, n, j; ! 409: floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2; ! 410: float32 posneg1, twoto63; ! 411: flag adjn, endflag; ! 412: ! 413: aSig = extractFloatx80Frac(a); ! 414: aExp = extractFloatx80Exp(a); ! 415: aSign = extractFloatx80Sign(a); ! 416: ! 417: if (aExp == 0x7FFF) { ! 418: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 419: float_raise(float_flag_invalid, status); ! 420: a.low = floatx80_default_nan_low; ! 421: a.high = floatx80_default_nan_high; ! 422: return a; ! 423: } ! 424: ! 425: if (aExp == 0 && aSig == 0) { ! 426: return packFloatx80(0, one_exp, one_sig); ! 427: } ! 428: ! 429: adjn = 1; ! 430: ! 431: SET_PREC; ! 432: ! 433: compact = floatx80_make_compact(aExp, aSig); ! 434: ! 435: fp0 = a; ! 436: ! 437: if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI ! 438: if (compact > 0x3FFF8000) { // |X| >= 15 PI ! 439: // REDUCEX ! 440: fp1 = packFloatx80(0, 0, 0); ! 441: if (compact == 0x7FFEFFFF) { ! 442: twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000)); ! 443: twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000)); ! 444: fp0 = floatx80_add(fp0, twopi1, status); ! 445: fp1 = fp0; ! 446: fp0 = floatx80_add(fp0, twopi2, status); ! 447: fp1 = floatx80_sub(fp1, fp0, status); ! 448: fp1 = floatx80_add(fp1, twopi2, status); ! 449: } ! 450: loop: ! 451: xSign = extractFloatx80Sign(fp0); ! 452: xExp = extractFloatx80Exp(fp0); ! 453: xExp -= 0x3FFF; ! 454: if (xExp <= 28) { ! 455: l = 0; ! 456: endflag = 1; ! 457: } else { ! 458: l = xExp - 27; ! 459: endflag = 0; ! 460: } ! 461: invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI ! 462: twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); ! 463: twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); ! 464: ! 465: twoto63 = 0x5F000000; ! 466: twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL ! 467: ! 468: fp2 = floatx80_mul(fp0, invtwopi, status); ! 469: fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), status); // THE FRACTIONAL PART OF FP2 IS ROUNDED ! 470: fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), status); // FP2 is N ! 471: fp4 = floatx80_mul(twopi1, fp2, status); // W = N*P1 ! 472: fp5 = floatx80_mul(twopi2, fp2, status); // w = N*P2 ! 473: fp3 = floatx80_add(fp4, fp5, status); // FP3 is P ! 474: fp4 = floatx80_sub(fp4, fp3, status); // W-P ! 475: fp0 = floatx80_sub(fp0, fp3, status); // FP0 is A := R - P ! 476: fp4 = floatx80_add(fp4, fp5, status); // FP4 is p = (W-P)+w ! 477: fp3 = fp0; // FP3 is A ! 478: fp1 = floatx80_sub(fp1, fp4, status); // FP1 is a := r - p ! 479: fp0 = floatx80_add(fp0, fp1, status); // FP0 is R := A+a ! 480: ! 481: if (endflag > 0) { ! 482: n = floatx80_to_int32(fp2, status); ! 483: goto sincont; ! 484: } ! 485: fp3 = floatx80_sub(fp3, fp0, status); // A-R ! 486: fp1 = floatx80_add(fp1, fp3, status); // FP1 is r := (A-R)+a ! 487: goto loop; ! 488: } else { ! 489: // SINSM ! 490: fp0 = float32_to_floatx80(0x3F800000, status); // 1 ! 491: ! 492: RESET_PREC; ! 493: ! 494: if (adjn) { ! 495: // COSTINY ! 496: a = floatx80_sub(fp0, float32_to_floatx80(0x00800000, status), status); ! 497: } else { ! 498: // SINTINY ! 499: a = floatx80_move(a, status); ! 500: } ! 501: float_raise(float_flag_inexact, status); ! 502: ! 503: return a; ! 504: } ! 505: } else { ! 506: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), status), status); // X*2/PI ! 507: ! 508: n = floatx80_to_int32(fp1, status); ! 509: j = 32 + n; ! 510: ! 511: fp0 = floatx80_sub(fp0, pi_tbl[j], status); // X-Y1 ! 512: fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), status); // FP0 IS R = (X-Y1)-Y2 ! 513: ! 514: sincont: ! 515: if ((n + adjn) & 1) { ! 516: // COSPOLY ! 517: fp0 = floatx80_mul(fp0, fp0, status); // FP0 IS S ! 518: fp1 = floatx80_mul(fp0, fp0, status); // FP1 IS T ! 519: fp2 = float64_to_floatx80(LIT64(0x3D2AC4D0D6011EE3), status); // B8 ! 520: fp3 = float64_to_floatx80(LIT64(0xBDA9396F9F45AC19), status); // B7 ! 521: ! 522: xSign = extractFloatx80Sign(fp0); // X IS S ! 523: xExp = extractFloatx80Exp(fp0); ! 524: xSig = extractFloatx80Frac(fp0); ! 525: ! 526: if (((n + adjn) >> 1) & 1) { ! 527: xSign ^= 1; ! 528: posneg1 = 0xBF800000; // -1 ! 529: } else { ! 530: xSign ^= 0; ! 531: posneg1 = 0x3F800000; // 1 ! 532: } // X IS NOW R'= SGN*R ! 533: ! 534: fp2 = floatx80_mul(fp2, fp1, status); // TB8 ! 535: fp3 = floatx80_mul(fp3, fp1, status); // TB7 ! 536: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3E21EED90612C972), status), status); // B6+TB8 ! 537: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE927E4FB79D9FCF), status), status); // B5+TB7 ! 538: fp2 = floatx80_mul(fp2, fp1, status); // T(B6+TB8) ! 539: fp3 = floatx80_mul(fp3, fp1, status); // T(B5+TB7) ! 540: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A01A01D423), status), status); // B4+T(B6+TB8) ! 541: fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); ! 542: fp3 = floatx80_add(fp3, fp4, status); // B3+T(B5+TB7) ! 543: fp2 = floatx80_mul(fp2, fp1, status); // T(B4+T(B6+TB8)) ! 544: fp1 = floatx80_mul(fp1, fp3, status); // T(B3+T(B5+TB7)) ! 545: fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); ! 546: fp2 = floatx80_add(fp2, fp4, status); // B2+T(B4+T(B6+TB8)) ! 547: fp1 = floatx80_add(fp1, float32_to_floatx80(0xBF000000, status), status); // B1+T(B3+T(B5+TB7)) ! 548: fp0 = floatx80_mul(fp0, fp2, status); // S(B2+T(B4+T(B6+TB8))) ! 549: fp0 = floatx80_add(fp0, fp1, status); // [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] ! 550: ! 551: x = packFloatx80(xSign, xExp, xSig); ! 552: fp0 = floatx80_mul(fp0, x, status); ! 553: ! 554: RESET_PREC; ! 555: ! 556: a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); ! 557: ! 558: float_raise(float_flag_inexact, status); ! 559: ! 560: return a; ! 561: } else { ! 562: // SINPOLY ! 563: xSign = extractFloatx80Sign(fp0); // X IS R ! 564: xExp = extractFloatx80Exp(fp0); ! 565: xSig = extractFloatx80Frac(fp0); ! 566: ! 567: xSign ^= ((n + adjn) >> 1) & 1; // X IS NOW R'= SGN*R ! 568: ! 569: fp0 = floatx80_mul(fp0, fp0, status); // FP0 IS S ! 570: fp1 = floatx80_mul(fp0, fp0, status); // FP1 IS T ! 571: fp3 = float64_to_floatx80(LIT64(0xBD6AAA77CCC994F5), status); // A7 ! 572: fp2 = float64_to_floatx80(LIT64(0x3DE612097AAE8DA1), status); // A6 ! 573: fp3 = floatx80_mul(fp3, fp1, status); // T*A7 ! 574: fp2 = floatx80_mul(fp2, fp1, status); // T*A6 ! 575: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE5AE6452A118AE4), status), status); // A5+T*A7 ! 576: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EC71DE3A5341531), status), status); // A4+T*A6 ! 577: fp3 = floatx80_mul(fp3, fp1, status); // T(A5+TA7) ! 578: fp2 = floatx80_mul(fp2, fp1, status); // T(A4+TA6) ! 579: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF2A01A01A018B59), status), status); // A3+T(A5+TA7) ! 580: fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); ! 581: fp2 = floatx80_add(fp2, fp4, status); // A2+T(A4+TA6) ! 582: fp1 = floatx80_mul(fp1, fp3, status); // T(A3+T(A5+TA7)) ! 583: fp2 = floatx80_mul(fp2, fp0, status); // S(A2+T(A4+TA6)) ! 584: fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); ! 585: fp1 = floatx80_add(fp1, fp4, status); // A1+T(A3+T(A5+TA7)) ! 586: fp1 = floatx80_add(fp1, fp2, status); // [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] ! 587: ! 588: x = packFloatx80(xSign, xExp, xSig); ! 589: fp0 = floatx80_mul(fp0, x, status); // R'*S ! 590: fp0 = floatx80_mul(fp0, fp1, status); // SIN(R')-R' ! 591: ! 592: RESET_PREC; ! 593: ! 594: a = floatx80_add(fp0, x, status); ! 595: ! 596: float_raise(float_flag_inexact, status); ! 597: ! 598: return a; ! 599: } ! 600: } ! 601: } ! 602: ! 603: /*---------------------------------------------------------------------------- ! 604: | Hyperbolic cosine ! 605: *----------------------------------------------------------------------------*/ ! 606: ! 607: floatx80 floatx80_cosh(floatx80 a, float_status *status) ! 608: { ! 609: flag aSign; ! 610: int32_t aExp; ! 611: uint64_t aSig; ! 612: ! 613: int32_t compact; ! 614: floatx80 fp0, fp1; ! 615: ! 616: aSig = extractFloatx80Frac(a); ! 617: aExp = extractFloatx80Exp(a); ! 618: aSign = extractFloatx80Sign(a); ! 619: ! 620: if (aExp == 0x7FFF) { ! 621: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 622: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 623: } ! 624: ! 625: if (aExp == 0 && aSig == 0) { ! 626: return packFloatx80(0, one_exp, one_sig); ! 627: } ! 628: ! 629: SET_PREC; ! 630: ! 631: compact = floatx80_make_compact(aExp, aSig); ! 632: ! 633: if (compact > 0x400CB167) { ! 634: if (compact > 0x400CB2B3) { ! 635: RESET_PREC; ! 636: return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, 0x8000, one_sig, 0, status); ! 637: } else { ! 638: fp0 = packFloatx80(0, aExp, aSig); ! 639: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x40C62D38D3D64634), status), status); ! 640: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x3D6F90AEB1E75CC7), status), status); ! 641: fp0 = floatx80_etox(fp0, status); ! 642: fp1 = packFloatx80(0, 0x7FFB, one_sig); ! 643: ! 644: RESET_PREC; ! 645: ! 646: a = floatx80_mul(fp0, fp1, status); ! 647: ! 648: float_raise(float_flag_inexact, status); ! 649: ! 650: return a; ! 651: } ! 652: } ! 653: ! 654: fp0 = packFloatx80(0, aExp, aSig); // |X| ! 655: fp0 = floatx80_etox(fp0, status); // EXP(|X|) ! 656: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x3F000000, status), status); // (1/2)*EXP(|X|) ! 657: fp1 = float32_to_floatx80(0x3E800000, status); // 1/4 ! 658: fp1 = floatx80_div(fp1, fp0, status); // 1/(2*EXP(|X|)) ! 659: ! 660: RESET_PREC; ! 661: ! 662: a = floatx80_add(fp0, fp1, status); ! 663: ! 664: float_raise(float_flag_inexact, status); ! 665: ! 666: return a; ! 667: } ! 668: ! 669: /*---------------------------------------------------------------------------- ! 670: | e to x ! 671: *----------------------------------------------------------------------------*/ ! 672: ! 673: floatx80 floatx80_etox(floatx80 a, float_status *status) ! 674: { ! 675: flag aSign; ! 676: int32_t aExp; ! 677: uint64_t aSig; ! 678: ! 679: int32_t compact, n, j, k, m, m1; ! 680: floatx80 fp0, fp1, fp2, fp3, l2, scale, adjscale; ! 681: flag adjflag; ! 682: ! 683: aSig = extractFloatx80Frac(a); ! 684: aExp = extractFloatx80Exp(a); ! 685: aSign = extractFloatx80Sign(a); ! 686: ! 687: if (aExp == 0x7FFF) { ! 688: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 689: if (aSign) return packFloatx80(0, 0, 0); ! 690: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 691: } ! 692: ! 693: if (aExp == 0 && aSig == 0) { ! 694: return packFloatx80(0, one_exp, one_sig); ! 695: } ! 696: ! 697: SET_PREC; ! 698: ! 699: adjflag = 0; ! 700: ! 701: if (aExp >= 0x3FBE) { // |X| >= 2^(-65) ! 702: compact = floatx80_make_compact(aExp, aSig); ! 703: ! 704: if (compact < 0x400CB167) { // |X| < 16380 log2 ! 705: fp0 = a; ! 706: fp1 = a; ! 707: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B, status), status); // 64/log2 * X ! 708: adjflag = 0; ! 709: n = floatx80_to_int32(fp0, status); // int(64/log2*X) ! 710: fp0 = int32_to_floatx80(n); ! 711: ! 712: j = n & 0x3F; // J = N mod 64 ! 713: m = n / 64; // NOTE: this is really arithmetic right shift by 6 ! 714: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity ! 715: m--; ! 716: } ! 717: m += 0x3FFF; // biased exponent of 2^(M) ! 718: ! 719: expcont1: ! 720: fp2 = fp0; // N ! 721: fp0 = floatx80_mul(fp0, float32_to_floatx80(0xBC317218, status), status); // N * L1, L1 = lead(-log2/64) ! 722: l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6)); ! 723: fp2 = floatx80_mul(fp2, l2, status); // N * L2, L1+L2 = -log2/64 ! 724: fp0 = floatx80_add(fp0, fp1, status); // X + N*L1 ! 725: fp0 = floatx80_add(fp0, fp2, status); // R ! 726: ! 727: fp1 = floatx80_mul(fp0, fp0, status); // S = R*R ! 728: fp2 = float32_to_floatx80(0x3AB60B70, status); // A5 ! 729: fp2 = floatx80_mul(fp2, fp1, status); // fp2 is S*A5 ! 730: fp3 = floatx80_mul(float32_to_floatx80(0x3C088895, status), fp1, status); // fp3 is S*A4 ! 731: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555554431), status), status); // fp2 is A3+S*A5 ! 732: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FC5555555554018), status), status); // fp3 is A2+S*A4 ! 733: fp2 = floatx80_mul(fp2, fp1, status); // fp2 is S*(A3+S*A5) ! 734: fp3 = floatx80_mul(fp3, fp1, status); // fp3 is S*(A2+S*A4) ! 735: fp2 = floatx80_add(fp2, float32_to_floatx80(0x3F000000, status), status); // fp2 is A1+S*(A3+S*A5) ! 736: fp3 = floatx80_mul(fp3, fp0, status); // fp3 IS R*S*(A2+S*A4) ! 737: fp2 = floatx80_mul(fp2, fp1, status); // fp2 IS S*(A1+S*(A3+S*A5)) ! 738: fp0 = floatx80_add(fp0, fp3, status); // fp0 IS R+R*S*(A2+S*A4) ! 739: fp0 = floatx80_add(fp0, fp2, status); // fp0 IS EXP(R) - 1 ! 740: ! 741: fp1 = exp_tbl[j]; ! 742: fp0 = floatx80_mul(fp0, fp1, status); // 2^(J/64)*(Exp(R)-1) ! 743: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status), status); // accurate 2^(J/64) ! 744: fp0 = floatx80_add(fp0, fp1, status); // 2^(J/64) + 2^(J/64)*(Exp(R)-1) ! 745: ! 746: scale = packFloatx80(0, m, one_sig); ! 747: if (adjflag) { ! 748: adjscale = packFloatx80(0, m1, one_sig); ! 749: fp0 = floatx80_mul(fp0, adjscale, status); ! 750: } ! 751: ! 752: RESET_PREC; ! 753: ! 754: a = floatx80_mul(fp0, scale, status); ! 755: ! 756: float_raise(float_flag_inexact, status); ! 757: ! 758: return a; ! 759: } else { // |X| >= 16380 log2 ! 760: if (compact > 0x400CB27C) { // |X| >= 16480 log2 ! 761: RESET_PREC; ! 762: if (aSign) { ! 763: a = roundAndPackFloatx80(status->floatx80_rounding_precision, 0, -0x1000, aSig, 0, status); ! 764: } else { ! 765: a = roundAndPackFloatx80(status->floatx80_rounding_precision, 0, 0x8000, aSig, 0, status); ! 766: } ! 767: float_raise(float_flag_inexact, status); ! 768: ! 769: return a; ! 770: } else { ! 771: fp0 = a; ! 772: fp1 = a; ! 773: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B, status), status); // 64/log2 * X ! 774: adjflag = 1; ! 775: n = floatx80_to_int32(fp0, status); // int(64/log2*X) ! 776: fp0 = int32_to_floatx80(n); ! 777: ! 778: j = n & 0x3F; // J = N mod 64 ! 779: k = n / 64; // NOTE: this is really arithmetic right shift by 6 ! 780: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity ! 781: k--; ! 782: } ! 783: m1 = k / 2; // NOTE: this is really arithmetic right shift by 1 ! 784: if (k < 0 && (k & 1)) { // arithmetic right shift is division and round towards minus infinity ! 785: m1--; ! 786: } ! 787: m = k - m1; ! 788: m1 += 0x3FFF; // biased exponent of 2^(M1) ! 789: m += 0x3FFF; // biased exponent of 2^(M) ! 790: ! 791: goto expcont1; ! 792: } ! 793: } ! 794: } else { // |X| < 2^(-65) ! 795: RESET_PREC; ! 796: ! 797: a = floatx80_add(a, float32_to_floatx80(0x3F800000, status), status); // 1 + X ! 798: ! 799: float_raise(float_flag_inexact, status); ! 800: ! 801: return a; ! 802: } ! 803: } ! 804: ! 805: /*---------------------------------------------------------------------------- ! 806: | e to x minus 1 ! 807: *----------------------------------------------------------------------------*/ ! 808: ! 809: floatx80 floatx80_etoxm1(floatx80 a, float_status *status) ! 810: { ! 811: flag aSign; ! 812: int32_t aExp; ! 813: uint64_t aSig; ! 814: ! 815: int32_t compact, n, j, m, m1; ! 816: floatx80 fp0, fp1, fp2, fp3, l2, sc, onebysc; ! 817: ! 818: aSig = extractFloatx80Frac(a); ! 819: aExp = extractFloatx80Exp(a); ! 820: aSign = extractFloatx80Sign(a); ! 821: ! 822: if (aExp == 0x7FFF) { ! 823: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 824: if (aSign) return packFloatx80(aSign, one_exp, one_sig); ! 825: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 826: } ! 827: ! 828: if (aExp == 0 && aSig == 0) { ! 829: return packFloatx80(aSign, 0, 0); ! 830: } ! 831: ! 832: SET_PREC; ! 833: ! 834: if (aExp >= 0x3FFD) { // |X| >= 1/4 ! 835: compact = floatx80_make_compact(aExp, aSig); ! 836: ! 837: if (compact <= 0x4004C215) { // |X| <= 70 log2 ! 838: fp0 = a; ! 839: fp1 = a; ! 840: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x42B8AA3B, status), status); // 64/log2 * X ! 841: n = floatx80_to_int32(fp0, status); // int(64/log2*X) ! 842: fp0 = int32_to_floatx80(n); ! 843: ! 844: j = n & 0x3F; // J = N mod 64 ! 845: m = n / 64; // NOTE: this is really arithmetic right shift by 6 ! 846: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity ! 847: m--; ! 848: } ! 849: m1 = -m; ! 850: //m += 0x3FFF; // biased exponent of 2^(M) ! 851: //m1 += 0x3FFF; // biased exponent of -2^(-M) ! 852: ! 853: fp2 = fp0; // N ! 854: fp0 = floatx80_mul(fp0, float32_to_floatx80(0xBC317218, status), status); // N * L1, L1 = lead(-log2/64) ! 855: l2 = packFloatx80(0, 0x3FDC, LIT64(0x82E308654361C4C6)); ! 856: fp2 = floatx80_mul(fp2, l2, status); // N * L2, L1+L2 = -log2/64 ! 857: fp0 = floatx80_add(fp0, fp1, status); // X + N*L1 ! 858: fp0 = floatx80_add(fp0, fp2, status); // R ! 859: ! 860: fp1 = floatx80_mul(fp0, fp0, status); // S = R*R ! 861: fp2 = float32_to_floatx80(0x3950097B, status); // A6 ! 862: fp2 = floatx80_mul(fp2, fp1, status); // fp2 is S*A6 ! 863: fp3 = floatx80_mul(float32_to_floatx80(0x3AB60B6A, status), fp1, status); // fp3 is S*A5 ! 864: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F81111111174385), status), status); // fp2 IS A4+S*A6 ! 865: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FA5555555554F5A), status), status); // fp3 is A3+S*A5 ! 866: fp2 = floatx80_mul(fp2, fp1, status); // fp2 IS S*(A4+S*A6) ! 867: fp3 = floatx80_mul(fp3, fp1, status); // fp3 IS S*(A3+S*A5) ! 868: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FC5555555555555), status), status); // fp2 IS A2+S*(A4+S*A6) ! 869: fp3 = floatx80_add(fp3, float32_to_floatx80(0x3F000000, status), status); // fp3 IS A1+S*(A3+S*A5) ! 870: fp2 = floatx80_mul(fp2, fp1, status); // fp2 IS S*(A2+S*(A4+S*A6)) ! 871: fp1 = floatx80_mul(fp1, fp3, status); // fp1 IS S*(A1+S*(A3+S*A5)) ! 872: fp2 = floatx80_mul(fp2, fp0, status); // fp2 IS R*S*(A2+S*(A4+S*A6)) ! 873: fp0 = floatx80_add(fp0, fp1, status); // fp0 IS R+S*(A1+S*(A3+S*A5)) ! 874: fp0 = floatx80_add(fp0, fp2, status); // fp0 IS EXP(R) - 1 ! 875: ! 876: fp0 = floatx80_mul(fp0, exp_tbl[j], status); // 2^(J/64)*(Exp(R)-1) ! 877: ! 878: if (m >= 64) { ! 879: fp1 = float32_to_floatx80(exp_tbl2[j], status); ! 880: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M) ! 881: fp1 = floatx80_add(fp1, onebysc, status); ! 882: fp0 = floatx80_add(fp0, fp1, status); ! 883: fp0 = floatx80_add(fp0, exp_tbl[j], status); ! 884: } else if (m < -3) { ! 885: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status), status); ! 886: fp0 = floatx80_add(fp0, exp_tbl[j], status); ! 887: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M) ! 888: fp0 = floatx80_add(fp0, onebysc, status); ! 889: } else { // -3 <= m <= 63 ! 890: fp1 = exp_tbl[j]; ! 891: fp0 = floatx80_add(fp0, float32_to_floatx80(exp_tbl2[j], status), status); ! 892: onebysc = packFloatx80(1, m1 + 0x3FFF, one_sig); // -2^(-M) ! 893: fp1 = floatx80_add(fp1, onebysc, status); ! 894: fp0 = floatx80_add(fp0, fp1, status); ! 895: } ! 896: ! 897: sc = packFloatx80(0, m + 0x3FFF, one_sig); ! 898: ! 899: RESET_PREC; ! 900: ! 901: a = floatx80_mul(fp0, sc, status); ! 902: ! 903: float_raise(float_flag_inexact, status); ! 904: ! 905: return a; ! 906: } else { // |X| > 70 log2 ! 907: if (aSign) { ! 908: fp0 = float32_to_floatx80(0xBF800000, status); // -1 ! 909: ! 910: RESET_PREC; ! 911: ! 912: a = floatx80_add(fp0, float32_to_floatx80(0x00800000, status), status); // -1 + 2^(-126) ! 913: ! 914: float_raise(float_flag_inexact, status); ! 915: ! 916: return a; ! 917: } else { ! 918: RESET_PREC; ! 919: ! 920: return floatx80_etox(a, status); ! 921: } ! 922: } ! 923: } else { // |X| < 1/4 ! 924: if (aExp >= 0x3FBE) { ! 925: fp0 = a; ! 926: fp0 = floatx80_mul(fp0, fp0, status); // S = X*X ! 927: fp1 = float32_to_floatx80(0x2F30CAA8, status); // B12 ! 928: fp1 = floatx80_mul(fp1, fp0, status); // S * B12 ! 929: fp2 = float32_to_floatx80(0x310F8290, status); // B11 ! 930: fp1 = floatx80_add(fp1, float32_to_floatx80(0x32D73220, status), status); // B10 ! 931: fp2 = floatx80_mul(fp2, fp0, status); ! 932: fp1 = floatx80_mul(fp1, fp0, status); ! 933: fp2 = floatx80_add(fp2, float32_to_floatx80(0x3493F281, status), status); // B9 ! 934: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3EC71DE3A5774682), status), status); // B8 ! 935: fp2 = floatx80_mul(fp2, fp0, status); ! 936: fp1 = floatx80_mul(fp1, fp0, status); ! 937: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A019D7CB68), status), status); // B7 ! 938: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3F2A01A01A019DF3), status), status); // B6 ! 939: fp2 = floatx80_mul(fp2, fp0, status); ! 940: fp1 = floatx80_mul(fp1, fp0, status); ! 941: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F56C16C16C170E2), status), status); // B5 ! 942: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3F81111111111111), status), status); // B4 ! 943: fp2 = floatx80_mul(fp2, fp0, status); ! 944: fp1 = floatx80_mul(fp1, fp0, status); ! 945: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555555555), status), status); // B3 ! 946: fp3 = packFloatx80(0, 0x3FFC, LIT64(0xAAAAAAAAAAAAAAAB)); ! 947: fp1 = floatx80_add(fp1, fp3, status); // B2 ! 948: fp2 = floatx80_mul(fp2, fp0, status); ! 949: fp1 = floatx80_mul(fp1, fp0, status); ! 950: ! 951: fp2 = floatx80_mul(fp2, fp0, status); ! 952: fp1 = floatx80_mul(fp1, a, status); ! 953: ! 954: fp0 = floatx80_mul(fp0, float32_to_floatx80(0x3F000000, status), status); // S*B1 ! 955: fp1 = floatx80_add(fp1, fp2, status); // Q ! 956: fp0 = floatx80_add(fp0, fp1, status); // S*B1+Q ! 957: ! 958: RESET_PREC; ! 959: ! 960: a = floatx80_add(fp0, a, status); ! 961: ! 962: float_raise(float_flag_inexact, status); ! 963: ! 964: return a; ! 965: } else { // |X| < 2^(-65) ! 966: sc = packFloatx80(1, 1, one_sig); ! 967: fp0 = a; ! 968: ! 969: if (aExp < 0x0033) { // |X| < 2^(-16382) ! 970: fp0 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x48B0000000000000), status), status); ! 971: fp0 = floatx80_add(fp0, sc, status); ! 972: ! 973: RESET_PREC; ! 974: ! 975: a = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3730000000000000), status), status); ! 976: } else { ! 977: RESET_PREC; ! 978: ! 979: a = floatx80_add(fp0, sc, status); ! 980: } ! 981: ! 982: float_raise(float_flag_inexact, status); ! 983: ! 984: return a; ! 985: } ! 986: } ! 987: } ! 988: ! 989: /*---------------------------------------------------------------------------- ! 990: | Log base 10 ! 991: *----------------------------------------------------------------------------*/ ! 992: ! 993: floatx80 floatx80_log10(floatx80 a, float_status *status) ! 994: { ! 995: flag aSign; ! 996: int32_t aExp; ! 997: uint64_t aSig; ! 998: ! 999: floatx80 fp0, fp1; ! 1000: ! 1001: aSig = extractFloatx80Frac(a); ! 1002: aExp = extractFloatx80Exp(a); ! 1003: aSign = extractFloatx80Sign(a); ! 1004: ! 1005: if (aExp == 0x7FFF) { ! 1006: if ((uint64_t) (aSig<<1)) propagateFloatx80NaNOneArg(a, status); ! 1007: if (aSign == 0) ! 1008: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 1009: } ! 1010: ! 1011: if (aExp == 0 && aSig == 0) { ! 1012: float_raise(float_flag_divbyzero, status); ! 1013: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low); ! 1014: } ! 1015: ! 1016: if (aSign) { ! 1017: float_raise(float_flag_invalid, status); ! 1018: a.low = floatx80_default_nan_low; ! 1019: a.high = floatx80_default_nan_high; ! 1020: return a; ! 1021: } ! 1022: ! 1023: SET_PREC; ! 1024: ! 1025: fp0 = floatx80_logn(a, status); ! 1026: fp1 = packFloatx80(0, 0x3FFD, LIT64(0xDE5BD8A937287195)); // INV_L10 ! 1027: ! 1028: RESET_PREC; ! 1029: ! 1030: a = floatx80_mul(fp0, fp1, status); // LOGN(X)*INV_L10 ! 1031: ! 1032: float_raise(float_flag_inexact, status); ! 1033: ! 1034: return a; ! 1035: } ! 1036: ! 1037: /*---------------------------------------------------------------------------- ! 1038: | Log base 2 ! 1039: *----------------------------------------------------------------------------*/ ! 1040: ! 1041: floatx80 floatx80_log2(floatx80 a, float_status *status) ! 1042: { ! 1043: flag aSign; ! 1044: int32_t aExp; ! 1045: uint64_t aSig; ! 1046: ! 1047: floatx80 fp0, fp1; ! 1048: ! 1049: aSig = extractFloatx80Frac(a); ! 1050: aExp = extractFloatx80Exp(a); ! 1051: aSign = extractFloatx80Sign(a); ! 1052: ! 1053: if (aExp == 0x7FFF) { ! 1054: if ((uint64_t) (aSig<<1)) propagateFloatx80NaNOneArg(a, status); ! 1055: if (aSign == 0) ! 1056: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 1057: } ! 1058: ! 1059: if (aExp == 0) { ! 1060: if (aSig == 0) { ! 1061: float_raise(float_flag_divbyzero, status); ! 1062: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low); ! 1063: } ! 1064: normalizeFloatx80Subnormal(aSig, &aExp, &aSig); ! 1065: } ! 1066: ! 1067: if (aSign) { ! 1068: float_raise(float_flag_invalid, status); ! 1069: a.low = floatx80_default_nan_low; ! 1070: a.high = floatx80_default_nan_high; ! 1071: return a; ! 1072: } ! 1073: ! 1074: SET_PREC; ! 1075: ! 1076: if (aSig == one_sig) { // X is 2^k ! 1077: RESET_PREC; ! 1078: ! 1079: a = int32_to_floatx80(aExp-0x3FFF); ! 1080: } else { ! 1081: fp0 = floatx80_logn(a, status); ! 1082: fp1 = packFloatx80(0, 0x3FFF, LIT64(0xB8AA3B295C17F0BC)); // INV_L2 ! 1083: ! 1084: RESET_PREC; ! 1085: ! 1086: a = floatx80_mul(fp0, fp1, status); // LOGN(X)*INV_L2 ! 1087: } ! 1088: ! 1089: float_raise(float_flag_inexact, status); ! 1090: ! 1091: return a; ! 1092: } ! 1093: ! 1094: /*---------------------------------------------------------------------------- ! 1095: | Log base e ! 1096: *----------------------------------------------------------------------------*/ ! 1097: ! 1098: floatx80 floatx80_logn(floatx80 a, float_status *status) ! 1099: { ! 1100: flag aSign; ! 1101: int32_t aExp; ! 1102: uint64_t aSig, fSig; ! 1103: ! 1104: int32_t compact, j, k, adjk; ! 1105: floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu; ! 1106: ! 1107: aSig = extractFloatx80Frac(a); ! 1108: aExp = extractFloatx80Exp(a); ! 1109: aSign = extractFloatx80Sign(a); ! 1110: ! 1111: if (aExp == 0x7FFF) { ! 1112: if ((uint64_t) (aSig<<1)) propagateFloatx80NaNOneArg(a, status); ! 1113: if (aSign == 0) ! 1114: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 1115: } ! 1116: ! 1117: adjk = 0; ! 1118: ! 1119: if (aExp == 0) { ! 1120: if (aSig == 0) { // zero ! 1121: float_raise(float_flag_divbyzero, status); ! 1122: return packFloatx80(1, 0x7FFF, floatx80_default_infinity_low); ! 1123: } ! 1124: #if 1 ! 1125: if ((aSig & one_sig) == 0) { // denormal ! 1126: normalizeFloatx80Subnormal(aSig, &aExp, &aSig); ! 1127: adjk = -100; ! 1128: aExp += 100; ! 1129: a = packFloatx80(aSign, aExp, aSig); ! 1130: } ! 1131: #else ! 1132: normalizeFloatx80Subnormal(aSig, &aExp, &aSig); ! 1133: #endif ! 1134: } ! 1135: ! 1136: if (aSign) { ! 1137: float_raise(float_flag_invalid, status); ! 1138: a.low = floatx80_default_nan_low; ! 1139: a.high = floatx80_default_nan_high; ! 1140: return a; ! 1141: } ! 1142: ! 1143: SET_PREC; ! 1144: ! 1145: compact = floatx80_make_compact(aExp, aSig); ! 1146: ! 1147: if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { // |X| < 15/16 or |X| > 17/16 ! 1148: k = aExp - 0x3FFF; ! 1149: k += adjk; ! 1150: fp1 = int32_to_floatx80(k); ! 1151: ! 1152: fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); ! 1153: j = (fSig >> 56) & 0x7E; // DISPLACEMENT FOR 1/F ! 1154: ! 1155: f = packFloatx80(0, 0x3FFF, fSig); // F ! 1156: fp0 = packFloatx80(0, 0x3FFF, aSig); // Y ! 1157: ! 1158: fp0 = floatx80_sub(fp0, f, status); // Y-F ! 1159: ! 1160: // LP1CONT1 ! 1161: fp0 = floatx80_mul(fp0, log_tbl[j], status); // FP0 IS U = (Y-F)/F ! 1162: logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); ! 1163: klog2 = floatx80_mul(fp1, logof2, status); // FP1 IS K*LOG2 ! 1164: fp2 = floatx80_mul(fp0, fp0, status); // FP2 IS V=U*U ! 1165: ! 1166: fp3 = fp2; ! 1167: fp1 = fp2; ! 1168: ! 1169: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3FC2499AB5E4040B), status), status); // V*A6 ! 1170: fp2 = floatx80_mul(fp2, float64_to_floatx80(LIT64(0xBFC555B5848CB7DB), status), status); // V*A5 ! 1171: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FC99999987D8730), status), status); // A4+V*A6 ! 1172: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFCFFFFFFF6F7E97), status), status); // A3+V*A5 ! 1173: fp1 = floatx80_mul(fp1, fp3, status); // V*(A4+V*A6) ! 1174: fp2 = floatx80_mul(fp2, fp3, status); // V*(A3+V*A5) ! 1175: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FD55555555555A4), status), status); // A2+V*(A4+V*A6) ! 1176: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFE0000000000008), status), status); // A1+V*(A3+V*A5) ! 1177: fp1 = floatx80_mul(fp1, fp3, status); // V*(A2+V*(A4+V*A6)) ! 1178: fp2 = floatx80_mul(fp2, fp3, status); // V*(A1+V*(A3+V*A5)) ! 1179: fp1 = floatx80_mul(fp1, fp0, status); // U*V*(A2+V*(A4+V*A6)) ! 1180: fp0 = floatx80_add(fp0, fp2, status); // U+V*(A1+V*(A3+V*A5)) ! 1181: ! 1182: fp1 = floatx80_add(fp1, log_tbl[j+1], status); // LOG(F)+U*V*(A2+V*(A4+V*A6)) ! 1183: fp0 = floatx80_add(fp0, fp1, status); // FP0 IS LOG(F) + LOG(1+U) ! 1184: ! 1185: RESET_PREC; ! 1186: ! 1187: a = floatx80_add(fp0, klog2, status); ! 1188: ! 1189: float_raise(float_flag_inexact, status); ! 1190: ! 1191: return a; ! 1192: } else { // |X-1| >= 1/16 ! 1193: fp0 = a; ! 1194: fp1 = a; ! 1195: fp1 = floatx80_sub(fp1, float32_to_floatx80(0x3F800000, status), status); // FP1 IS X-1 ! 1196: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // FP0 IS X+1 ! 1197: fp1 = floatx80_add(fp1, fp1, status); // FP1 IS 2(X-1) ! 1198: ! 1199: // LP1CONT2 ! 1200: fp1 = floatx80_div(fp1, fp0, status); // U ! 1201: saveu = fp1; ! 1202: fp0 = floatx80_mul(fp1, fp1, status); // FP0 IS V = U*U ! 1203: fp1 = floatx80_mul(fp0, fp0, status); // FP1 IS W = V*V ! 1204: ! 1205: fp3 = float64_to_floatx80(LIT64(0x3F175496ADD7DAD6), status); // B5 ! 1206: fp2 = float64_to_floatx80(LIT64(0x3F3C71C2FE80C7E0), status); // B4 ! 1207: fp3 = floatx80_mul(fp3, fp1, status); // W*B5 ! 1208: fp2 = floatx80_mul(fp2, fp1, status); // W*B4 ! 1209: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3F624924928BCCFF), status), status); // B3+W*B5 ! 1210: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F899999999995EC), status), status); // B2+W*B4 ! 1211: fp1 = floatx80_mul(fp1, fp3, status); // W*(B3+W*B5) ! 1212: fp2 = floatx80_mul(fp2, fp0, status); // V*(B2+W*B4) ! 1213: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FB5555555555555), status), status); // B1+W*(B3+W*B5) ! 1214: ! 1215: fp0 = floatx80_mul(fp0, saveu, status); // FP0 IS U*V ! 1216: fp1 = floatx80_add(fp1, fp2, status); // B1+W*(B3+W*B5) + V*(B2+W*B4) ! 1217: fp0 = floatx80_mul(fp0, fp1, status); // U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] ) ! 1218: ! 1219: RESET_PREC; ! 1220: ! 1221: a = floatx80_add(fp0, saveu, status); ! 1222: ! 1223: float_raise(float_flag_inexact, status); ! 1224: ! 1225: return a; ! 1226: } ! 1227: } ! 1228: ! 1229: /*---------------------------------------------------------------------------- ! 1230: | Log base e of x plus 1 ! 1231: *----------------------------------------------------------------------------*/ ! 1232: ! 1233: floatx80 floatx80_lognp1(floatx80 a, float_status *status) ! 1234: { ! 1235: flag aSign; ! 1236: int32_t aExp; ! 1237: uint64_t aSig, fSig; ! 1238: ! 1239: int32_t compact, j, k; ! 1240: floatx80 fp0, fp1, fp2, fp3, f, logof2, klog2, saveu; ! 1241: ! 1242: aSig = extractFloatx80Frac(a); ! 1243: aExp = extractFloatx80Exp(a); ! 1244: aSign = extractFloatx80Sign(a); ! 1245: ! 1246: if (aExp == 0x7FFF) { ! 1247: if ((uint64_t) (aSig<<1)) propagateFloatx80NaNOneArg(a, status); ! 1248: if (aSign) { ! 1249: float_raise(float_flag_invalid, status); ! 1250: a.low = floatx80_default_nan_low; ! 1251: a.high = floatx80_default_nan_high; ! 1252: return a; ! 1253: } ! 1254: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 1255: } ! 1256: ! 1257: if (aExp == 0 && aSig == 0) { ! 1258: return packFloatx80(aSign, 0, 0); ! 1259: } ! 1260: ! 1261: if (aSign && aExp >= one_exp) { ! 1262: if (aExp == one_exp && aSig == one_sig) { ! 1263: float_raise(float_flag_divbyzero, status); ! 1264: packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low); ! 1265: } ! 1266: float_raise(float_flag_invalid, status); ! 1267: a.low = floatx80_default_nan_low; ! 1268: a.high = floatx80_default_nan_high; ! 1269: return a; ! 1270: } ! 1271: ! 1272: if (aExp < 0x3f99 || (aExp == 0x3f99 && aSig == one_sig)) { // <= min threshold ! 1273: float_raise(float_flag_inexact, status); ! 1274: return floatx80_move(a, status); ! 1275: } ! 1276: ! 1277: SET_PREC; ! 1278: ! 1279: compact = floatx80_make_compact(aExp, aSig); ! 1280: ! 1281: fp0 = a; // Z ! 1282: fp1 = a; ! 1283: ! 1284: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // X = (1+Z) ! 1285: ! 1286: aExp = extractFloatx80Exp(fp0); ! 1287: aSig = extractFloatx80Frac(fp0); ! 1288: ! 1289: compact = floatx80_make_compact(aExp, aSig); ! 1290: ! 1291: if (compact < 0x3FFE8000 || compact > 0x3FFFC000) { // |X| < 1/2 or |X| > 3/2 ! 1292: k = aExp - 0x3FFF; ! 1293: fp1 = int32_to_floatx80(k); ! 1294: ! 1295: fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); ! 1296: j = (fSig >> 56) & 0x7E; // DISPLACEMENT FOR 1/F ! 1297: ! 1298: f = packFloatx80(0, 0x3FFF, fSig); // F ! 1299: fp0 = packFloatx80(0, 0x3FFF, aSig); // Y ! 1300: ! 1301: fp0 = floatx80_sub(fp0, f, status); // Y-F ! 1302: ! 1303: lp1cont1: ! 1304: // LP1CONT1 ! 1305: fp0 = floatx80_mul(fp0, log_tbl[j], status); // FP0 IS U = (Y-F)/F ! 1306: logof2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); ! 1307: klog2 = floatx80_mul(fp1, logof2, status); // FP1 IS K*LOG2 ! 1308: fp2 = floatx80_mul(fp0, fp0, status); // FP2 IS V=U*U ! 1309: ! 1310: fp3 = fp2; ! 1311: fp1 = fp2; ! 1312: ! 1313: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3FC2499AB5E4040B), status), status); // V*A6 ! 1314: fp2 = floatx80_mul(fp2, float64_to_floatx80(LIT64(0xBFC555B5848CB7DB), status), status); // V*A5 ! 1315: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FC99999987D8730), status), status); // A4+V*A6 ! 1316: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFCFFFFFFF6F7E97), status), status); // A3+V*A5 ! 1317: fp1 = floatx80_mul(fp1, fp3, status); // V*(A4+V*A6) ! 1318: fp2 = floatx80_mul(fp2, fp3, status); // V*(A3+V*A5) ! 1319: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FD55555555555A4), status), status); // A2+V*(A4+V*A6) ! 1320: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBFE0000000000008), status), status); // A1+V*(A3+V*A5) ! 1321: fp1 = floatx80_mul(fp1, fp3, status); // V*(A2+V*(A4+V*A6)) ! 1322: fp2 = floatx80_mul(fp2, fp3, status); // V*(A1+V*(A3+V*A5)) ! 1323: fp1 = floatx80_mul(fp1, fp0, status); // U*V*(A2+V*(A4+V*A6)) ! 1324: fp0 = floatx80_add(fp0, fp2, status); // U+V*(A1+V*(A3+V*A5)) ! 1325: ! 1326: fp1 = floatx80_add(fp1, log_tbl[j+1], status); // LOG(F)+U*V*(A2+V*(A4+V*A6)) ! 1327: fp0 = floatx80_add(fp0, fp1, status); // FP0 IS LOG(F) + LOG(1+U) ! 1328: ! 1329: RESET_PREC; ! 1330: ! 1331: a = floatx80_add(fp0, klog2, status); ! 1332: ! 1333: float_raise(float_flag_inexact, status); ! 1334: ! 1335: return a; ! 1336: } else if (compact < 0x3FFEF07D || compact > 0x3FFF8841) { // |X| < 1/16 or |X| > -1/16 ! 1337: // LP1CARE ! 1338: fSig = (aSig & LIT64(0xFE00000000000000)) | LIT64(0x0100000000000000); ! 1339: f = packFloatx80(0, 0x3FFF, fSig); // F ! 1340: j = (fSig >> 56) & 0x7E; // DISPLACEMENT FOR 1/F ! 1341: ! 1342: if (compact >= 0x3FFF8000) { // 1+Z >= 1 ! 1343: // KISZERO ! 1344: fp0 = floatx80_sub(float32_to_floatx80(0x3F800000, status), f, status); // 1-F ! 1345: fp0 = floatx80_add(fp0, fp1, status); // FP0 IS Y-F = (1-F)+Z ! 1346: fp1 = packFloatx80(0, 0, 0); // K = 0 ! 1347: } else { ! 1348: // KISNEG ! 1349: fp0 = floatx80_sub(float32_to_floatx80(0x40000000, status), f, status); // 2-F ! 1350: fp1 = floatx80_add(fp1, fp1, status); // 2Z ! 1351: fp0 = floatx80_add(fp0, fp1, status); // FP0 IS Y-F = (2-F)+2Z ! 1352: fp1 = packFloatx80(1, one_exp, one_sig); // K = -1 ! 1353: } ! 1354: goto lp1cont1; ! 1355: } else { ! 1356: // LP1ONE16 ! 1357: fp1 = floatx80_add(fp1, fp1, status); // FP1 IS 2Z ! 1358: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // FP0 IS 1+X ! 1359: ! 1360: // LP1CONT2 ! 1361: fp1 = floatx80_div(fp1, fp0, status); // U ! 1362: saveu = fp1; ! 1363: fp0 = floatx80_mul(fp1, fp1, status); // FP0 IS V = U*U ! 1364: fp1 = floatx80_mul(fp0, fp0, status); // FP1 IS W = V*V ! 1365: ! 1366: fp3 = float64_to_floatx80(LIT64(0x3F175496ADD7DAD6), status); // B5 ! 1367: fp2 = float64_to_floatx80(LIT64(0x3F3C71C2FE80C7E0), status); // B4 ! 1368: fp3 = floatx80_mul(fp3, fp1, status); // W*B5 ! 1369: fp2 = floatx80_mul(fp2, fp1, status); // W*B4 ! 1370: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3F624924928BCCFF), status), status); // B3+W*B5 ! 1371: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3F899999999995EC), status), status); // B2+W*B4 ! 1372: fp1 = floatx80_mul(fp1, fp3, status); // W*(B3+W*B5) ! 1373: fp2 = floatx80_mul(fp2, fp0, status); // V*(B2+W*B4) ! 1374: fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3FB5555555555555), status), status); // B1+W*(B3+W*B5) ! 1375: ! 1376: fp0 = floatx80_mul(fp0, saveu, status); // FP0 IS U*V ! 1377: fp1 = floatx80_add(fp1, fp2, status); // B1+W*(B3+W*B5) + V*(B2+W*B4) ! 1378: fp0 = floatx80_mul(fp0, fp1, status); // U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] ) ! 1379: ! 1380: RESET_PREC; ! 1381: ! 1382: a = floatx80_add(fp0, saveu, status); ! 1383: ! 1384: float_raise(float_flag_inexact, status); ! 1385: ! 1386: return a; ! 1387: } ! 1388: } ! 1389: ! 1390: /*---------------------------------------------------------------------------- ! 1391: | Sine ! 1392: *----------------------------------------------------------------------------*/ ! 1393: ! 1394: floatx80 floatx80_sin(floatx80 a, float_status *status) ! 1395: { ! 1396: flag aSign, xSign; ! 1397: int32_t aExp, xExp; ! 1398: uint64_t aSig, xSig; ! 1399: ! 1400: int32_t compact, l, n, j; ! 1401: floatx80 fp0, fp1, fp2, fp3, fp4, fp5, x, invtwopi, twopi1, twopi2; ! 1402: float32 posneg1, twoto63; ! 1403: flag adjn, endflag; ! 1404: ! 1405: aSig = extractFloatx80Frac(a); ! 1406: aExp = extractFloatx80Exp(a); ! 1407: aSign = extractFloatx80Sign(a); ! 1408: ! 1409: if (aExp == 0x7FFF) { ! 1410: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 1411: float_raise(float_flag_invalid, status); ! 1412: a.low = floatx80_default_nan_low; ! 1413: a.high = floatx80_default_nan_high; ! 1414: return a; ! 1415: } ! 1416: ! 1417: if (aExp == 0 && aSig == 0) { ! 1418: return packFloatx80(aSign, 0, 0); ! 1419: } ! 1420: ! 1421: adjn = 0; ! 1422: ! 1423: SET_PREC; ! 1424: ! 1425: compact = floatx80_make_compact(aExp, aSig); ! 1426: ! 1427: fp0 = a; ! 1428: ! 1429: if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI ! 1430: if (compact > 0x3FFF8000) { // |X| >= 15 PI ! 1431: // REDUCEX ! 1432: fp1 = packFloatx80(0, 0, 0); ! 1433: if (compact == 0x7FFEFFFF) { ! 1434: twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000)); ! 1435: twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000)); ! 1436: fp0 = floatx80_add(fp0, twopi1, status); ! 1437: fp1 = fp0; ! 1438: fp0 = floatx80_add(fp0, twopi2, status); ! 1439: fp1 = floatx80_sub(fp1, fp0, status); ! 1440: fp1 = floatx80_add(fp1, twopi2, status); ! 1441: } ! 1442: loop: ! 1443: xSign = extractFloatx80Sign(fp0); ! 1444: xExp = extractFloatx80Exp(fp0); ! 1445: xExp -= 0x3FFF; ! 1446: if (xExp <= 28) { ! 1447: l = 0; ! 1448: endflag = 1; ! 1449: } else { ! 1450: l = xExp - 27; ! 1451: endflag = 0; ! 1452: } ! 1453: invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI ! 1454: twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); ! 1455: twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); ! 1456: ! 1457: twoto63 = 0x5F000000; ! 1458: twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL ! 1459: ! 1460: fp2 = floatx80_mul(fp0, invtwopi, status); ! 1461: fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), status); // THE FRACTIONAL PART OF FP2 IS ROUNDED ! 1462: fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), status); // FP2 is N ! 1463: fp4 = floatx80_mul(twopi1, fp2, status); // W = N*P1 ! 1464: fp5 = floatx80_mul(twopi2, fp2, status); // w = N*P2 ! 1465: fp3 = floatx80_add(fp4, fp5, status); // FP3 is P ! 1466: fp4 = floatx80_sub(fp4, fp3, status); // W-P ! 1467: fp0 = floatx80_sub(fp0, fp3, status); // FP0 is A := R - P ! 1468: fp4 = floatx80_add(fp4, fp5, status); // FP4 is p = (W-P)+w ! 1469: fp3 = fp0; // FP3 is A ! 1470: fp1 = floatx80_sub(fp1, fp4, status); // FP1 is a := r - p ! 1471: fp0 = floatx80_add(fp0, fp1, status); // FP0 is R := A+a ! 1472: ! 1473: if (endflag > 0) { ! 1474: n = floatx80_to_int32(fp2, status); ! 1475: goto sincont; ! 1476: } ! 1477: fp3 = floatx80_sub(fp3, fp0, status); // A-R ! 1478: fp1 = floatx80_add(fp1, fp3, status); // FP1 is r := (A-R)+a ! 1479: goto loop; ! 1480: } else { ! 1481: // SINSM ! 1482: fp0 = float32_to_floatx80(0x3F800000, status); // 1 ! 1483: ! 1484: RESET_PREC; ! 1485: ! 1486: if (adjn) { ! 1487: // COSTINY ! 1488: a = floatx80_sub(fp0, float32_to_floatx80(0x00800000, status), status); ! 1489: } else { ! 1490: // SINTINY ! 1491: a = floatx80_move(a, status); ! 1492: } ! 1493: float_raise(float_flag_inexact, status); ! 1494: ! 1495: return a; ! 1496: } ! 1497: } else { ! 1498: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), status), status); // X*2/PI ! 1499: ! 1500: n = floatx80_to_int32(fp1, status); ! 1501: j = 32 + n; ! 1502: ! 1503: fp0 = floatx80_sub(fp0, pi_tbl[j], status); // X-Y1 ! 1504: fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[j], status), status); // FP0 IS R = (X-Y1)-Y2 ! 1505: ! 1506: sincont: ! 1507: if ((n + adjn) & 1) { ! 1508: // COSPOLY ! 1509: fp0 = floatx80_mul(fp0, fp0, status); // FP0 IS S ! 1510: fp1 = floatx80_mul(fp0, fp0, status); // FP1 IS T ! 1511: fp2 = float64_to_floatx80(LIT64(0x3D2AC4D0D6011EE3), status); // B8 ! 1512: fp3 = float64_to_floatx80(LIT64(0xBDA9396F9F45AC19), status); // B7 ! 1513: ! 1514: xSign = extractFloatx80Sign(fp0); // X IS S ! 1515: xExp = extractFloatx80Exp(fp0); ! 1516: xSig = extractFloatx80Frac(fp0); ! 1517: ! 1518: if (((n + adjn) >> 1) & 1) { ! 1519: xSign ^= 1; ! 1520: posneg1 = 0xBF800000; // -1 ! 1521: } else { ! 1522: xSign ^= 0; ! 1523: posneg1 = 0x3F800000; // 1 ! 1524: } // X IS NOW R'= SGN*R ! 1525: ! 1526: fp2 = floatx80_mul(fp2, fp1, status); // TB8 ! 1527: fp3 = floatx80_mul(fp3, fp1, status); // TB7 ! 1528: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3E21EED90612C972), status), status); // B6+TB8 ! 1529: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE927E4FB79D9FCF), status), status); // B5+TB7 ! 1530: fp2 = floatx80_mul(fp2, fp1, status); // T(B6+TB8) ! 1531: fp3 = floatx80_mul(fp3, fp1, status); // T(B5+TB7) ! 1532: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A01A01D423), status), status); // B4+T(B6+TB8) ! 1533: fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); ! 1534: fp3 = floatx80_add(fp3, fp4, status); // B3+T(B5+TB7) ! 1535: fp2 = floatx80_mul(fp2, fp1, status); // T(B4+T(B6+TB8)) ! 1536: fp1 = floatx80_mul(fp1, fp3, status); // T(B3+T(B5+TB7)) ! 1537: fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); ! 1538: fp2 = floatx80_add(fp2, fp4, status); // B2+T(B4+T(B6+TB8)) ! 1539: fp1 = floatx80_add(fp1, float32_to_floatx80(0xBF000000, status), status); // B1+T(B3+T(B5+TB7)) ! 1540: fp0 = floatx80_mul(fp0, fp2, status); // S(B2+T(B4+T(B6+TB8))) ! 1541: fp0 = floatx80_add(fp0, fp1, status); // [B1+T(B3+T(B5+TB7))]+[S(B2+T(B4+T(B6+TB8)))] ! 1542: ! 1543: x = packFloatx80(xSign, xExp, xSig); ! 1544: fp0 = floatx80_mul(fp0, x, status); ! 1545: ! 1546: RESET_PREC; ! 1547: ! 1548: a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); ! 1549: ! 1550: float_raise(float_flag_inexact, status); ! 1551: ! 1552: return a; ! 1553: } else { ! 1554: // SINPOLY ! 1555: xSign = extractFloatx80Sign(fp0); // X IS R ! 1556: xExp = extractFloatx80Exp(fp0); ! 1557: xSig = extractFloatx80Frac(fp0); ! 1558: ! 1559: xSign ^= ((n + adjn) >> 1) & 1; // X IS NOW R'= SGN*R ! 1560: ! 1561: fp0 = floatx80_mul(fp0, fp0, status); // FP0 IS S ! 1562: fp1 = floatx80_mul(fp0, fp0, status); // FP1 IS T ! 1563: fp3 = float64_to_floatx80(LIT64(0xBD6AAA77CCC994F5), status); // A7 ! 1564: fp2 = float64_to_floatx80(LIT64(0x3DE612097AAE8DA1), status); // A6 ! 1565: fp3 = floatx80_mul(fp3, fp1, status); // T*A7 ! 1566: fp2 = floatx80_mul(fp2, fp1, status); // T*A6 ! 1567: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBE5AE6452A118AE4), status), status); // A5+T*A7 ! 1568: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EC71DE3A5341531), status), status); // A4+T*A6 ! 1569: fp3 = floatx80_mul(fp3, fp1, status); // T(A5+TA7) ! 1570: fp2 = floatx80_mul(fp2, fp1, status); // T(A4+TA6) ! 1571: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF2A01A01A018B59), status), status); // A3+T(A5+TA7) ! 1572: fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); ! 1573: fp2 = floatx80_add(fp2, fp4, status); // A2+T(A4+TA6) ! 1574: fp1 = floatx80_mul(fp1, fp3, status); // T(A3+T(A5+TA7)) ! 1575: fp2 = floatx80_mul(fp2, fp0, status); // S(A2+T(A4+TA6)) ! 1576: fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); ! 1577: fp1 = floatx80_add(fp1, fp4, status); // A1+T(A3+T(A5+TA7)) ! 1578: fp1 = floatx80_add(fp1, fp2, status); // [A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))] ! 1579: ! 1580: x = packFloatx80(xSign, xExp, xSig); ! 1581: fp0 = floatx80_mul(fp0, x, status); // R'*S ! 1582: fp0 = floatx80_mul(fp0, fp1, status); // SIN(R')-R' ! 1583: ! 1584: RESET_PREC; ! 1585: ! 1586: a = floatx80_add(fp0, x, status); ! 1587: ! 1588: float_raise(float_flag_inexact, status); ! 1589: ! 1590: return a; ! 1591: } ! 1592: } ! 1593: } ! 1594: ! 1595: /*---------------------------------------------------------------------------- ! 1596: | Hyperbolic sine ! 1597: *----------------------------------------------------------------------------*/ ! 1598: ! 1599: floatx80 floatx80_sinh(floatx80 a, float_status *status) ! 1600: { ! 1601: flag aSign; ! 1602: int32_t aExp; ! 1603: uint64_t aSig; ! 1604: ! 1605: int32_t compact; ! 1606: floatx80 fp0, fp1, fp2; ! 1607: float32 fact; ! 1608: ! 1609: aSig = extractFloatx80Frac(a); ! 1610: aExp = extractFloatx80Exp(a); ! 1611: aSign = extractFloatx80Sign(a); ! 1612: ! 1613: if (aExp == 0x7FFF) { ! 1614: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 1615: return packFloatx80(aSign, 0x7FFF, floatx80_default_infinity_low); ! 1616: } ! 1617: ! 1618: if (aExp == 0 && aSig == 0) { ! 1619: return packFloatx80(aSign, 0, 0); ! 1620: } ! 1621: ! 1622: SET_PREC; ! 1623: ! 1624: compact = floatx80_make_compact(aExp, aSig); ! 1625: ! 1626: if (compact > 0x400CB167) { ! 1627: // SINHBIG ! 1628: if (compact > 0x400CB2B3) { ! 1629: RESET_PREC; ! 1630: ! 1631: return roundAndPackFloatx80(status->floatx80_rounding_precision, aSign, 0x8000, aSig, 0, status); ! 1632: } else { ! 1633: fp0 = floatx80_abs(a, status); // Y = |X| ! 1634: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x40C62D38D3D64634), status), status); // (|X|-16381LOG2_LEAD) ! 1635: fp0 = floatx80_sub(fp0, float64_to_floatx80(LIT64(0x3D6F90AEB1E75CC7), status), status); // |X| - 16381 LOG2, ACCURATE ! 1636: fp0 = floatx80_etox(fp0, status); ! 1637: fp2 = packFloatx80(aSign, 0x7FFB, one_sig); ! 1638: ! 1639: RESET_PREC; ! 1640: ! 1641: a = floatx80_mul(fp0, fp2, status); ! 1642: ! 1643: float_raise(float_flag_inexact, status); ! 1644: ! 1645: return a; ! 1646: } ! 1647: } else { // |X| < 16380 LOG2 ! 1648: fp0 = floatx80_abs(a, status); // Y = |X| ! 1649: fp0 = floatx80_etoxm1(fp0, status); // FP0 IS Z = EXPM1(Y) ! 1650: fp1 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // 1+Z ! 1651: fp2 = fp0; ! 1652: fp0 = floatx80_div(fp0, fp1, status); // Z/(1+Z) ! 1653: fp0 = floatx80_add(fp0, fp2, status); ! 1654: ! 1655: fact = 0x3F000000; ! 1656: fact |= aSign ? 0x80000000 : 0x00000000; ! 1657: ! 1658: RESET_PREC; ! 1659: ! 1660: a = floatx80_mul(fp0, float32_to_floatx80(fact, status), status); ! 1661: ! 1662: float_raise(float_flag_inexact, status); ! 1663: ! 1664: return a; ! 1665: } ! 1666: } ! 1667: ! 1668: /*---------------------------------------------------------------------------- ! 1669: | Tangent ! 1670: *----------------------------------------------------------------------------*/ ! 1671: ! 1672: floatx80 floatx80_tan(floatx80 a, float_status *status) ! 1673: { ! 1674: flag aSign, xSign; ! 1675: int32_t aExp, xExp; ! 1676: uint64_t aSig, xSig; ! 1677: ! 1678: int32_t compact, l, n, j; ! 1679: floatx80 fp0, fp1, fp2, fp3, fp4, fp5, invtwopi, twopi1, twopi2; ! 1680: float32 twoto63; ! 1681: flag endflag; ! 1682: ! 1683: aSig = extractFloatx80Frac(a); ! 1684: aExp = extractFloatx80Exp(a); ! 1685: aSign = extractFloatx80Sign(a); ! 1686: ! 1687: if (aExp == 0x7FFF) { ! 1688: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 1689: float_raise(float_flag_invalid, status); ! 1690: a.low = floatx80_default_nan_low; ! 1691: a.high = floatx80_default_nan_high; ! 1692: return a; ! 1693: } ! 1694: ! 1695: if (aExp == 0 && aSig == 0) { ! 1696: return packFloatx80(aSign, 0, 0); ! 1697: } ! 1698: ! 1699: SET_PREC; ! 1700: ! 1701: compact = floatx80_make_compact(aExp, aSig); ! 1702: ! 1703: fp0 = a; ! 1704: ! 1705: if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI ! 1706: if (compact > 0x3FFF8000) { // |X| >= 15 PI ! 1707: // REDUCEX ! 1708: fp1 = packFloatx80(0, 0, 0); ! 1709: if (compact == 0x7FFEFFFF) { ! 1710: twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000)); ! 1711: twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000)); ! 1712: fp0 = floatx80_add(fp0, twopi1, status); ! 1713: fp1 = fp0; ! 1714: fp0 = floatx80_add(fp0, twopi2, status); ! 1715: fp1 = floatx80_sub(fp1, fp0, status); ! 1716: fp1 = floatx80_add(fp1, twopi2, status); ! 1717: } ! 1718: loop: ! 1719: xSign = extractFloatx80Sign(fp0); ! 1720: xExp = extractFloatx80Exp(fp0); ! 1721: xExp -= 0x3FFF; ! 1722: if (xExp <= 28) { ! 1723: l = 0; ! 1724: endflag = 1; ! 1725: } else { ! 1726: l = xExp - 27; ! 1727: endflag = 0; ! 1728: } ! 1729: invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI ! 1730: twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); ! 1731: twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); ! 1732: ! 1733: twoto63 = 0x5F000000; ! 1734: twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL ! 1735: ! 1736: fp2 = floatx80_mul(fp0, invtwopi, status); ! 1737: fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), status); // THE FRACTIONAL PART OF FP2 IS ROUNDED ! 1738: fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), status); // FP2 is N ! 1739: fp4 = floatx80_mul(twopi1, fp2, status); // W = N*P1 ! 1740: fp5 = floatx80_mul(twopi2, fp2, status); // w = N*P2 ! 1741: fp3 = floatx80_add(fp4, fp5, status); // FP3 is P ! 1742: fp4 = floatx80_sub(fp4, fp3, status); // W-P ! 1743: fp0 = floatx80_sub(fp0, fp3, status); // FP0 is A := R - P ! 1744: fp4 = floatx80_add(fp4, fp5, status); // FP4 is p = (W-P)+w ! 1745: fp3 = fp0; // FP3 is A ! 1746: fp1 = floatx80_sub(fp1, fp4, status); // FP1 is a := r - p ! 1747: fp0 = floatx80_add(fp0, fp1, status); // FP0 is R := A+a ! 1748: ! 1749: if (endflag > 0) { ! 1750: n = floatx80_to_int32(fp2, status); ! 1751: goto tancont; ! 1752: } ! 1753: fp3 = floatx80_sub(fp3, fp0, status); // A-R ! 1754: fp1 = floatx80_add(fp1, fp3, status); // FP1 is r := (A-R)+a ! 1755: goto loop; ! 1756: } else { ! 1757: RESET_PREC; ! 1758: ! 1759: a = floatx80_move(a, status); ! 1760: ! 1761: float_raise(float_flag_inexact, status); ! 1762: ! 1763: return a; ! 1764: } ! 1765: } else { ! 1766: fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), status), 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), status); // FP0 IS R = (X-Y1)-Y2 ! 1773: ! 1774: tancont: ! 1775: if (n & 1) { ! 1776: // NODD ! 1777: fp1 = fp0; // R ! 1778: fp0 = floatx80_mul(fp0, fp0, status); // S = R*R ! 1779: fp3 = float64_to_floatx80(LIT64(0x3EA0B759F50F8688), status); // Q4 ! 1780: fp2 = float64_to_floatx80(LIT64(0xBEF2BAA5A8924F04), status); // P3 ! 1781: fp3 = floatx80_mul(fp3, fp0, status); // SQ4 ! 1782: fp2 = floatx80_mul(fp2, fp0, status); // SP3 ! 1783: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF346F59B39BA65F), status), status); // Q3+SQ4 ! 1784: fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00)); ! 1785: fp2 = floatx80_add(fp2, fp4, status); // P2+SP3 ! 1786: fp3 = floatx80_mul(fp3, fp0, status); // S(Q3+SQ4) ! 1787: fp2 = floatx80_mul(fp2, fp0, status); // S(P2+SP3) ! 1788: fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1)); ! 1789: fp3 = floatx80_add(fp3, fp4, status); // Q2+S(Q3+SQ4) ! 1790: fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA)); ! 1791: fp2 = floatx80_add(fp2, fp4, status); // P1+S(P2+SP3) ! 1792: fp3 = floatx80_mul(fp3, fp0, status); // S(Q2+S(Q3+SQ4)) ! 1793: fp2 = floatx80_mul(fp2, fp0, status); // S(P1+S(P2+SP3)) ! 1794: fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE)); ! 1795: fp3 = floatx80_add(fp3, fp4, status); // Q1+S(Q2+S(Q3+SQ4)) ! 1796: fp2 = floatx80_mul(fp2, fp1, status); // RS(P1+S(P2+SP3)) ! 1797: fp0 = floatx80_mul(fp0, fp3, status); // S(Q1+S(Q2+S(Q3+SQ4))) ! 1798: fp1 = floatx80_add(fp1, fp2, status); // R+RS(P1+S(P2+SP3)) ! 1799: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // 1+S(Q1+S(Q2+S(Q3+SQ4))) ! 1800: ! 1801: xSign = extractFloatx80Sign(fp1); ! 1802: xExp = extractFloatx80Exp(fp1); ! 1803: xSig = extractFloatx80Frac(fp1); ! 1804: xSign ^= 1; ! 1805: fp1 = packFloatx80(xSign, xExp, xSig); ! 1806: ! 1807: RESET_PREC; ! 1808: ! 1809: a = floatx80_div(fp0, fp1, status); ! 1810: ! 1811: float_raise(float_flag_inexact, status); ! 1812: ! 1813: return a; ! 1814: } else { ! 1815: fp1 = floatx80_mul(fp0, fp0, status); // S = R*R ! 1816: fp3 = float64_to_floatx80(LIT64(0x3EA0B759F50F8688), status); // Q4 ! 1817: fp2 = float64_to_floatx80(LIT64(0xBEF2BAA5A8924F04), status); // P3 ! 1818: fp3 = floatx80_mul(fp3, fp1, status); // SQ4 ! 1819: fp2 = floatx80_mul(fp2, fp1, status); // SP3 ! 1820: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0xBF346F59B39BA65F), status), status); // Q3+SQ4 ! 1821: fp4 = packFloatx80(0, 0x3FF6, LIT64(0xE073D3FC199C4A00)); ! 1822: fp2 = floatx80_add(fp2, fp4, status); // P2+SP3 ! 1823: fp3 = floatx80_mul(fp3, fp1, status); // S(Q3+SQ4) ! 1824: fp2 = floatx80_mul(fp2, fp1, status); // S(P2+SP3) ! 1825: fp4 = packFloatx80(0, 0x3FF9, LIT64(0xD23CD68415D95FA1)); ! 1826: fp3 = floatx80_add(fp3, fp4, status); // Q2+S(Q3+SQ4) ! 1827: fp4 = packFloatx80(1, 0x3FFC, LIT64(0x8895A6C5FB423BCA)); ! 1828: fp2 = floatx80_add(fp2, fp4, status); // P1+S(P2+SP3) ! 1829: fp3 = floatx80_mul(fp3, fp1, status); // S(Q2+S(Q3+SQ4)) ! 1830: fp2 = floatx80_mul(fp2, fp1, status); // S(P1+S(P2+SP3)) ! 1831: fp4 = packFloatx80(1, 0x3FFD, LIT64(0xEEF57E0DA84BC8CE)); ! 1832: fp3 = floatx80_add(fp3, fp4, status); // Q1+S(Q2+S(Q3+SQ4)) ! 1833: fp2 = floatx80_mul(fp2, fp0, status); // RS(P1+S(P2+SP3)) ! 1834: fp1 = floatx80_mul(fp1, fp3, status); // S(Q1+S(Q2+S(Q3+SQ4))) ! 1835: fp0 = floatx80_add(fp0, fp2, status); // R+RS(P1+S(P2+SP3)) ! 1836: fp1 = floatx80_add(fp1, float32_to_floatx80(0x3F800000, status), status); // 1+S(Q1+S(Q2+S(Q3+SQ4))) ! 1837: ! 1838: RESET_PREC; ! 1839: ! 1840: a = floatx80_div(fp0, fp1, status); ! 1841: ! 1842: float_raise(float_flag_inexact, status); ! 1843: ! 1844: return a; ! 1845: } ! 1846: } ! 1847: } ! 1848: ! 1849: /*---------------------------------------------------------------------------- ! 1850: | Hyperbolic tangent ! 1851: *----------------------------------------------------------------------------*/ ! 1852: ! 1853: floatx80 floatx80_tanh(floatx80 a, float_status *status) ! 1854: { ! 1855: flag aSign, vSign; ! 1856: int32_t aExp, vExp; ! 1857: uint64_t aSig, vSig; ! 1858: ! 1859: int32_t compact; ! 1860: floatx80 fp0, fp1; ! 1861: float32 sign; ! 1862: ! 1863: aSig = extractFloatx80Frac(a); ! 1864: aExp = extractFloatx80Exp(a); ! 1865: aSign = extractFloatx80Sign(a); ! 1866: ! 1867: if (aExp == 0x7FFF) { ! 1868: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 1869: return packFloatx80(aSign, one_exp, one_sig); ! 1870: } ! 1871: ! 1872: if (aExp == 0 && aSig == 0) { ! 1873: return packFloatx80(aSign, 0, 0); ! 1874: } ! 1875: ! 1876: SET_PREC; ! 1877: ! 1878: compact = floatx80_make_compact(aExp, aSig); ! 1879: ! 1880: if (compact < 0x3FD78000 || compact > 0x3FFFDDCE) { ! 1881: // TANHBORS ! 1882: if (compact < 0x3FFF8000) { ! 1883: // TANHSM ! 1884: RESET_PREC; ! 1885: ! 1886: a = floatx80_move(a, status); ! 1887: ! 1888: float_raise(float_flag_inexact, status); ! 1889: ! 1890: return a; ! 1891: } else { ! 1892: if (compact > 0x40048AA1) { ! 1893: // TANHHUGE ! 1894: sign = 0x3F800000; ! 1895: sign |= aSign ? 0x80000000 : 0x00000000; ! 1896: fp0 = float32_to_floatx80(sign, status); ! 1897: sign &= 0x80000000; ! 1898: sign ^= 0x80800000; // -SIGN(X)*EPS ! 1899: ! 1900: RESET_PREC; ! 1901: ! 1902: a = floatx80_add(fp0, float32_to_floatx80(sign, status), status); ! 1903: ! 1904: float_raise(float_flag_inexact, status); ! 1905: ! 1906: return a; ! 1907: } else { ! 1908: fp0 = packFloatx80(0, aExp+1, aSig); // Y = 2|X| ! 1909: fp0 = floatx80_etox(fp0, status); // FP0 IS EXP(Y) ! 1910: fp0 = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // EXP(Y)+1 ! 1911: sign = aSign ? 0x80000000 : 0x00000000; ! 1912: fp1 = floatx80_div(float32_to_floatx80(sign^0xC0000000, status), fp0, status); // -SIGN(X)*2 / [EXP(Y)+1] ! 1913: fp0 = float32_to_floatx80(sign | 0x3F800000, status); // SIGN ! 1914: ! 1915: RESET_PREC; ! 1916: ! 1917: a = floatx80_add(fp1, fp0, status); ! 1918: ! 1919: float_raise(float_flag_inexact, status); ! 1920: ! 1921: return a; ! 1922: } ! 1923: } ! 1924: } else { // 2**(-40) < |X| < (5/2)LOG2 ! 1925: fp0 = packFloatx80(0, aExp+1, aSig); // Y = 2|X| ! 1926: fp0 = floatx80_etoxm1(fp0, status); // FP0 IS Z = EXPM1(Y) ! 1927: fp1 = floatx80_add(fp0, float32_to_floatx80(0x40000000, status), status); // Z+2 ! 1928: ! 1929: vSign = extractFloatx80Sign(fp1); ! 1930: vExp = extractFloatx80Exp(fp1); ! 1931: vSig = extractFloatx80Frac(fp1); ! 1932: ! 1933: fp1 = packFloatx80(vSign ^ aSign, vExp, vSig); ! 1934: ! 1935: RESET_PREC; ! 1936: ! 1937: a = floatx80_div(fp0, fp1, status); ! 1938: ! 1939: float_raise(float_flag_inexact, status); ! 1940: ! 1941: return a; ! 1942: } ! 1943: } ! 1944: ! 1945: /*---------------------------------------------------------------------------- ! 1946: | 10 to x ! 1947: *----------------------------------------------------------------------------*/ ! 1948: ! 1949: floatx80 floatx80_tentox(floatx80 a, float_status *status) ! 1950: { ! 1951: flag aSign; ! 1952: int32_t aExp; ! 1953: uint64_t aSig; ! 1954: ! 1955: int32_t compact, n, j, l, m, m1; ! 1956: floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2; ! 1957: ! 1958: aSig = extractFloatx80Frac(a); ! 1959: aExp = extractFloatx80Exp(a); ! 1960: aSign = extractFloatx80Sign(a); ! 1961: ! 1962: if (aExp == 0x7FFF) { ! 1963: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 1964: if (aSign) return packFloatx80(0, 0, 0); ! 1965: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 1966: } ! 1967: ! 1968: if (aExp == 0 && aSig == 0) { ! 1969: return packFloatx80(0, one_exp, one_sig); ! 1970: } ! 1971: ! 1972: SET_PREC; ! 1973: ! 1974: fp0 = a; ! 1975: ! 1976: compact = floatx80_make_compact(aExp, aSig); ! 1977: ! 1978: if (compact < 0x3FB98000 || compact > 0x400B9B07) { // |X| > 16480 LOG2/LOG10 or |X| < 2^(-70) ! 1979: if (compact > 0x3FFF8000) { // |X| > 16480 ! 1980: RESET_PREC; ! 1981: ! 1982: if (aSign) { ! 1983: return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, -0x1000, aSig, 0, status); ! 1984: } else { ! 1985: return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, 0x8000, aSig, 0, status); ! 1986: } ! 1987: } else { // |X| < 2^(-70) ! 1988: RESET_PREC; ! 1989: ! 1990: a = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // 1 + X ! 1991: ! 1992: float_raise(float_flag_inexact, status); ! 1993: ! 1994: return a; ! 1995: } ! 1996: } else { // 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10 ! 1997: fp1 = fp0; // X ! 1998: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x406A934F0979A371), status), status); // X*64*LOG10/LOG2 ! 1999: n = floatx80_to_int32(fp1, status); // N=INT(X*64*LOG10/LOG2) ! 2000: fp1 = int32_to_floatx80(n); ! 2001: ! 2002: j = n & 0x3F; ! 2003: l = n / 64; // NOTE: this is really arithmetic right shift by 6 ! 2004: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity ! 2005: l--; ! 2006: } ! 2007: m = l / 2; // NOTE: this is really arithmetic right shift by 1 ! 2008: if (l < 0 && (l & 1)) { // arithmetic right shift is division and round towards minus infinity ! 2009: m--; ! 2010: } ! 2011: m1 = l - m; ! 2012: m1 += 0x3FFF; // ADJFACT IS 2^(M') ! 2013: ! 2014: adjfact = packFloatx80(0, m1, one_sig); ! 2015: fact1 = exp2_tbl[j]; ! 2016: fact1.high += m; ! 2017: fact2.high = exp2_tbl2[j]>>16; ! 2018: fact2.high += m; ! 2019: fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF); ! 2020: fact2.low <<= 48; ! 2021: ! 2022: fp2 = fp1; // N ! 2023: fp1 = floatx80_mul(fp1, float64_to_floatx80(LIT64(0x3F734413509F8000), status), status); // N*(LOG2/64LOG10)_LEAD ! 2024: fp3 = packFloatx80(1, 0x3FCD, LIT64(0xC0219DC1DA994FD2)); ! 2025: fp2 = floatx80_mul(fp2, fp3, status); // N*(LOG2/64LOG10)_TRAIL ! 2026: fp0 = floatx80_sub(fp0, fp1, status); // X - N L_LEAD ! 2027: fp0 = floatx80_sub(fp0, fp2, status); // X - N L_TRAIL ! 2028: fp2 = packFloatx80(0, 0x4000, LIT64(0x935D8DDDAAA8AC17)); // LOG10 ! 2029: fp0 = floatx80_mul(fp0, fp2, status); // R ! 2030: ! 2031: // EXPR ! 2032: fp1 = floatx80_mul(fp0, fp0, status); // S = R*R ! 2033: fp2 = float64_to_floatx80(LIT64(0x3F56C16D6F7BD0B2), status); // A5 ! 2034: fp3 = float64_to_floatx80(LIT64(0x3F811112302C712C), status); // A4 ! 2035: fp2 = floatx80_mul(fp2, fp1, status); // S*A5 ! 2036: fp3 = floatx80_mul(fp3, fp1, status); // S*A4 ! 2037: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555554CC1), status), status); // A3+S*A5 ! 2038: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FC5555555554A54), status), status); // A2+S*A4 ! 2039: fp2 = floatx80_mul(fp2, fp1, status); // S*(A3+S*A5) ! 2040: fp3 = floatx80_mul(fp3, fp1, status); // S*(A2+S*A4) ! 2041: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FE0000000000000), status), status); // A1+S*(A3+S*A5) ! 2042: fp3 = floatx80_mul(fp3, fp0, status); // R*S*(A2+S*A4) ! 2043: ! 2044: fp2 = floatx80_mul(fp2, fp1, status); // S*(A1+S*(A3+S*A5)) ! 2045: fp0 = floatx80_add(fp0, fp3, status); // R+R*S*(A2+S*A4) ! 2046: fp0 = floatx80_add(fp0, fp2, status); // EXP(R) - 1 ! 2047: ! 2048: fp0 = floatx80_mul(fp0, fact1, status); ! 2049: fp0 = floatx80_add(fp0, fact2, status); ! 2050: fp0 = floatx80_add(fp0, fact1, status); ! 2051: ! 2052: RESET_PREC; ! 2053: ! 2054: a = floatx80_mul(fp0, adjfact, status); ! 2055: ! 2056: float_raise(float_flag_inexact, status); ! 2057: ! 2058: return a; ! 2059: } ! 2060: } ! 2061: ! 2062: /*---------------------------------------------------------------------------- ! 2063: | 2 to x ! 2064: *----------------------------------------------------------------------------*/ ! 2065: ! 2066: floatx80 floatx80_twotox(floatx80 a, float_status *status) ! 2067: { ! 2068: flag aSign; ! 2069: int32_t aExp; ! 2070: uint64_t aSig; ! 2071: ! 2072: int32_t compact, n, j, l, m, m1; ! 2073: floatx80 fp0, fp1, fp2, fp3, adjfact, fact1, fact2; ! 2074: ! 2075: aSig = extractFloatx80Frac(a); ! 2076: aExp = extractFloatx80Exp(a); ! 2077: aSign = extractFloatx80Sign(a); ! 2078: ! 2079: if (aExp == 0x7FFF) { ! 2080: if ((uint64_t) (aSig<<1)) return propagateFloatx80NaNOneArg(a, status); ! 2081: if (aSign) return packFloatx80(0, 0, 0); ! 2082: return packFloatx80(0, 0x7FFF, floatx80_default_infinity_low); ! 2083: } ! 2084: ! 2085: if (aExp == 0 && aSig == 0) { ! 2086: return packFloatx80(0, one_exp, one_sig); ! 2087: } ! 2088: ! 2089: SET_PREC; ! 2090: ! 2091: fp0 = a; ! 2092: ! 2093: compact = floatx80_make_compact(aExp, aSig); ! 2094: ! 2095: if (compact < 0x3FB98000 || compact > 0x400D80C0) { // |X| > 16480 or |X| < 2^(-70) ! 2096: if (compact > 0x3FFF8000) { // |X| > 16480 ! 2097: RESET_PREC;; ! 2098: ! 2099: if (aSign) { ! 2100: return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, -0x1000, aSig, 0, status); ! 2101: } else { ! 2102: return roundAndPackFloatx80(status->floatx80_rounding_precision, 0, 0x8000, aSig, 0, status); ! 2103: } ! 2104: } else { // |X| < 2^(-70) ! 2105: RESET_PREC;; ! 2106: ! 2107: a = floatx80_add(fp0, float32_to_floatx80(0x3F800000, status), status); // 1 + X ! 2108: ! 2109: float_raise(float_flag_inexact, status); ! 2110: ! 2111: return a; ! 2112: } ! 2113: } else { // 2^(-70) <= |X| <= 16480 ! 2114: fp1 = fp0; // X ! 2115: fp1 = floatx80_mul(fp1, float32_to_floatx80(0x42800000, status), status); // X * 64 ! 2116: n = floatx80_to_int32(fp1, status); ! 2117: fp1 = int32_to_floatx80(n); ! 2118: j = n & 0x3F; ! 2119: l = n / 64; // NOTE: this is really arithmetic right shift by 6 ! 2120: if (n < 0 && j) { // arithmetic right shift is division and round towards minus infinity ! 2121: l--; ! 2122: } ! 2123: m = l / 2; // NOTE: this is really arithmetic right shift by 1 ! 2124: if (l < 0 && (l & 1)) { // arithmetic right shift is division and round towards minus infinity ! 2125: m--; ! 2126: } ! 2127: m1 = l - m; ! 2128: m1 += 0x3FFF; // ADJFACT IS 2^(M') ! 2129: ! 2130: adjfact = packFloatx80(0, m1, one_sig); ! 2131: fact1 = exp2_tbl[j]; ! 2132: fact1.high += m; ! 2133: fact2.high = exp2_tbl2[j]>>16; ! 2134: fact2.high += m; ! 2135: fact2.low = (uint64_t)(exp2_tbl2[j] & 0xFFFF); ! 2136: fact2.low <<= 48; ! 2137: ! 2138: fp1 = floatx80_mul(fp1, float32_to_floatx80(0x3C800000, status), status); // (1/64)*N ! 2139: fp0 = floatx80_sub(fp0, fp1, status); // X - (1/64)*INT(64 X) ! 2140: fp2 = packFloatx80(0, 0x3FFE, LIT64(0xB17217F7D1CF79AC)); // LOG2 ! 2141: fp0 = floatx80_mul(fp0, fp2, status); // R ! 2142: ! 2143: // EXPR ! 2144: fp1 = floatx80_mul(fp0, fp0, status); // S = R*R ! 2145: fp2 = float64_to_floatx80(LIT64(0x3F56C16D6F7BD0B2), status); // A5 ! 2146: fp3 = float64_to_floatx80(LIT64(0x3F811112302C712C), status); // A4 ! 2147: fp2 = floatx80_mul(fp2, fp1, status); // S*A5 ! 2148: fp3 = floatx80_mul(fp3, fp1, status); // S*A4 ! 2149: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FA5555555554CC1), status), status); // A3+S*A5 ! 2150: fp3 = floatx80_add(fp3, float64_to_floatx80(LIT64(0x3FC5555555554A54), status), status); // A2+S*A4 ! 2151: fp2 = floatx80_mul(fp2, fp1, status); // S*(A3+S*A5) ! 2152: fp3 = floatx80_mul(fp3, fp1, status); // S*(A2+S*A4) ! 2153: fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3FE0000000000000), status), status); // A1+S*(A3+S*A5) ! 2154: fp3 = floatx80_mul(fp3, fp0, status); // R*S*(A2+S*A4) ! 2155: ! 2156: fp2 = floatx80_mul(fp2, fp1, status); // S*(A1+S*(A3+S*A5)) ! 2157: fp0 = floatx80_add(fp0, fp3, status); // R+R*S*(A2+S*A4) ! 2158: fp0 = floatx80_add(fp0, fp2, status); // EXP(R) - 1 ! 2159: ! 2160: fp0 = floatx80_mul(fp0, fact1, status); ! 2161: fp0 = floatx80_add(fp0, fact2, status); ! 2162: fp0 = floatx80_add(fp0, fact1, status); ! 2163: ! 2164: RESET_PREC; ! 2165: ! 2166: a = floatx80_mul(fp0, adjfact, status); ! 2167: ! 2168: float_raise(float_flag_inexact, status); ! 2169: ! 2170: return a; ! 2171: } ! 2172: }
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.