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