From 2337abe0a213ae3cec06f052f74dc898e84e2317 Mon Sep 17 00:00:00 2001 From: Toni Wilen Date: Sat, 10 Apr 2021 21:48:47 +0300 Subject: [PATCH] Softfloat update (from Previous) --- fpp.cpp | 5 +- fpp_softfloat.cpp | 5 + include/fpp.h | 2 + softfloat/softfloat.cpp | 224 +++++++++++++++--------------- softfloat/softfloat.h | 1 + softfloat/softfloat_fpsp.cpp | 261 +++++++++++++++++++++++++++++++++++ 6 files changed, 379 insertions(+), 119 deletions(-) diff --git a/fpp.cpp b/fpp.cpp index 2a1ce6a61..1606a2093 100644 --- a/fpp.cpp +++ b/fpp.cpp @@ -127,6 +127,7 @@ FPP_AB fpp_cosh; FPP_ABP fpp_neg; FPP_AB fpp_acos; FPP_AB fpp_cos; +FPP_ABC fpp_sincos; FPP_AB fpp_getexp; FPP_AB fpp_getman; FPP_ABP fpp_div; @@ -3109,9 +3110,7 @@ static bool fp_arithmetic(fpdata *src, fpdata *dst, int extra) case 0x35: /* FSINCOS */ case 0x36: /* FSINCOS */ case 0x37: /* FSINCOS */ - fpp_cos(dst, src); - regs.fp[extra & 7] = *dst; - fpp_sin(dst, src); + fpp_sincos(dst, src, ®s.fp[extra & 7]); break; case 0x38: /* FCMP */ case 0x39: diff --git a/fpp_softfloat.cpp b/fpp_softfloat.cpp index aa059bfbc..f1253bb6d 100644 --- a/fpp_softfloat.cpp +++ b/fpp_softfloat.cpp @@ -511,6 +511,10 @@ static void fp_cos(fpdata *a, fpdata *b) { a->fpx = floatx80_cos(b->fpx, &fs); } +static void fp_sincos(fpdata * a, fpdata * b, fpdata * c) +{ + a->fpx = floatx80_sincos(b->fpx, &c->fpx, &fs); +} /* Functions for converting between float formats */ static const fptype twoto32 = 4294967296.0; @@ -814,6 +818,7 @@ void fp_init_softfloat(int fpu_model) fpp_neg = fp_neg; fpp_acos = fp_acos; fpp_cos = fp_cos; + fpp_sincos = fp_sincos; fpp_getexp = fp_getexp; fpp_getman = fp_getman; fpp_div = fp_div; diff --git a/include/fpp.h b/include/fpp.h index 2fbc49ead..f82ee4879 100644 --- a/include/fpp.h +++ b/include/fpp.h @@ -51,6 +51,7 @@ typedef void (*FPP_ABQS)(fpdata*, fpdata*, uae_u64*, uae_u8*); typedef void (*FPP_AB)(fpdata*, fpdata*); typedef void (*FPP_ABP)(fpdata*, fpdata*, int); typedef void (*FPP_A)(fpdata*); +typedef void (*FPP_ABC)(fpdata*, fpdata*, fpdata*); typedef bool (*FPP_IS)(fpdata*); typedef void (*FPP_SET_MODE)(uae_u32); @@ -151,6 +152,7 @@ extern FPP_AB fpp_cosh; extern FPP_ABP fpp_neg; extern FPP_AB fpp_acos; extern FPP_AB fpp_cos; +extern FPP_ABC fpp_sincos; extern FPP_AB fpp_getexp; extern FPP_AB fpp_getman; extern FPP_ABP fpp_div; diff --git a/softfloat/softfloat.cpp b/softfloat/softfloat.cpp index da9503379..04b730550 100644 --- a/softfloat/softfloat.cpp +++ b/softfloat/softfloat.cpp @@ -2879,87 +2879,84 @@ floatx80 floatx80_rem(floatx80 a, floatx80 b, float_status *status) } #else // 09-01-2017: Modified version for Previous -floatx80 floatx80_rem( floatx80 a, floatx80 b, uint64_t *q, flag *s, float_status *status ) +floatx80 floatx80_rem(floatx80 a, floatx80 b, uint64_t* q, flag* s, float_status* status) { flag aSign, bSign, zSign; int32_t aExp, bExp, expDiff; uint64_t aSig0, aSig1, bSig; uint64_t qTemp, term0, term1, alternateASig0, alternateASig1; - - aSig0 = extractFloatx80Frac( a ); - aExp = extractFloatx80Exp( a ); - aSign = extractFloatx80Sign( a ); - bSig = extractFloatx80Frac( b ); - bExp = extractFloatx80Exp( b ); - bSign = extractFloatx80Sign( b ); - if ( aExp == 0x7FFF ) { - if ( (uint64_t) ( aSig0<<1 ) - || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) { - return propagateFloatx80NaN( a, b, status ); + aSig0 = extractFloatx80Frac(a); + aExp = extractFloatx80Exp(a); + aSign = extractFloatx80Sign(a); + bSig = extractFloatx80Frac(b); + bExp = extractFloatx80Exp(b); + bSign = extractFloatx80Sign(b); + + *s = 0; + *q = 0; + + if (aExp == 0x7FFF) { + if ((uint64_t)(aSig0 << 1) + || ((bExp == 0x7FFF) && (uint64_t)(bSig << 1))) { + return propagateFloatx80NaN(a, b, status); } goto invalid; } - if ( bExp == 0x7FFF ) { - if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status ); + if (bExp == 0x7FFF) { + if ((uint64_t)(bSig << 1)) return propagateFloatx80NaN(a, b, status); *s = (aSign != bSign); - *q = 0; - return a; + return normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, aSign, aExp, aSig0, 0, status);; } - if ( bExp == 0 ) { - if ( bSig == 0 ) { + if (bExp == 0) { + if (bSig == 0) { invalid: - float_raise( float_flag_invalid, status ); - return floatx80_default_nan(status); + float_raise(float_flag_invalid, status); + return floatx80_default_nan(status); } - normalizeFloatx80Subnormal( bSig, &bExp, &bSig ); + normalizeFloatx80Subnormal(bSig, &bExp, &bSig); } - if ( aExp == 0 ) { -#ifdef SOFTFLOAT_68K - if ( aSig0 == 0 ) { + if (aExp == 0) { + if (aSig0 == 0) { *s = (aSign != bSign); - *q = 0; return a; } -#else - if ( (uint64_t) ( aSig0<<1 ) == 0 ) return a; -#endif - normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 ); + normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0); } - bSig |= LIT64( 0x8000000000000000 ); + bSig |= LIT64(0x8000000000000000); zSign = aSign; expDiff = aExp - bExp; *s = (aSign != bSign); aSig1 = 0; - if ( expDiff < 0 ) { - if ( expDiff < -1 ) return a; - shift128Right( aSig0, 0, 1, &aSig0, &aSig1 ); + if (expDiff < 0) { + if (expDiff < -1) return a; + shift128Right(aSig0, 0, 1, &aSig0, &aSig1); expDiff = 0; } - qTemp = ( bSig <= aSig0 ); - if ( qTemp ) aSig0 -= bSig; - *q = ( expDiff > 63 ) ? 0 : ( qTemp< 63) ? 0 : (qTemp << expDiff); expDiff -= 64; - while ( 0 < expDiff ) { - qTemp = estimateDiv128To64( aSig0, aSig1, bSig ); - qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0; - mul64To128( bSig, qTemp, &term0, &term1 ); - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 ); - *q = ( expDiff > 63 ) ? 0 : ( qTemp< 63) ? 0 : (qTemp << expDiff); expDiff -= 62; } expDiff += 64; - if ( 0 < expDiff ) { - qTemp = estimateDiv128To64( aSig0, aSig1, bSig ); - qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0; + if (0 < expDiff) { + qTemp = estimateDiv128To64(aSig0, aSig1, bSig); + qTemp = (2 < qTemp) ? qTemp - 2 : 0; qTemp >>= 64 - expDiff; - mul64To128( bSig, qTemp<<( 64 - expDiff ), &term0, &term1 ); - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 ); - while ( le128( term0, term1, aSig0, aSig1 ) ) { + mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1); + sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); + shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1); + while (le128(term0, term1, aSig0, aSig1)) { ++qTemp; - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); + sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); } *q += qTemp; } @@ -2967,118 +2964,113 @@ floatx80 floatx80_rem( floatx80 a, floatx80 b, uint64_t *q, flag *s, float_statu term1 = 0; term0 = bSig; } - sub128( term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1 ); - if ( lt128( alternateASig0, alternateASig1, aSig0, aSig1 ) - || ( eq128( alternateASig0, alternateASig1, aSig0, aSig1 ) - && ( qTemp & 1 ) ) + sub128(term0, term1, aSig0, aSig1, &alternateASig0, &alternateASig1); + if (lt128(alternateASig0, alternateASig1, aSig0, aSig1) + || (eq128(alternateASig0, alternateASig1, aSig0, aSig1) + && (qTemp & 1)) ) { aSig0 = alternateASig0; aSig1 = alternateASig1; - zSign = ! zSign; - ++*q; + zSign = !zSign; + ++* q; } return - normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, - zSign, bExp + expDiff, aSig0, aSig1, status ); - + normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, + zSign, bExp + expDiff, aSig0, aSig1, status); + } #endif // End of modification - #ifdef SOFTFLOAT_68K // 08-01-2017: Added for Previous /*---------------------------------------------------------------------------- | Returns the modulo remainder of the extended double-precision floating-point | value `a' with respect to the corresponding value `b'. *----------------------------------------------------------------------------*/ -floatx80 floatx80_mod( floatx80 a, floatx80 b, uint64_t *q, flag *s, float_status *status ) +floatx80 floatx80_mod(floatx80 a, floatx80 b, uint64_t* q, flag* s, float_status* status) { flag aSign, bSign, zSign; int32_t aExp, bExp, expDiff; uint64_t aSig0, aSig1, bSig; uint64_t qTemp, term0, term1; - - aSig0 = extractFloatx80Frac( a ); - aExp = extractFloatx80Exp( a ); - aSign = extractFloatx80Sign( a ); - bSig = extractFloatx80Frac( b ); - bExp = extractFloatx80Exp( b ); - bSign = extractFloatx80Sign( b ); - if ( aExp == 0x7FFF ) { - if ( (uint64_t) ( aSig0<<1 ) - || ( ( bExp == 0x7FFF ) && (uint64_t) ( bSig<<1 ) ) ) { - return propagateFloatx80NaN( a, b, status ); + aSig0 = extractFloatx80Frac(a); + aExp = extractFloatx80Exp(a); + aSign = extractFloatx80Sign(a); + bSig = extractFloatx80Frac(b); + bExp = extractFloatx80Exp(b); + bSign = extractFloatx80Sign(b); + + *s = 0; + *q = 0; + + if (aExp == 0x7FFF) { + if ((uint64_t)(aSig0 << 1) + || ((bExp == 0x7FFF) && (uint64_t)(bSig << 1))) { + return propagateFloatx80NaN(a, b, status); } goto invalid; } - if ( bExp == 0x7FFF ) { - if ( (uint64_t) ( bSig<<1 ) ) return propagateFloatx80NaN( a, b, status ); + if (bExp == 0x7FFF) { + if ((uint64_t)(bSig << 1)) return propagateFloatx80NaN(a, b, status); *s = (aSign != bSign); - *q = 0; - return a; + return normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, aSign, aExp, aSig0, 0, status); } - if ( bExp == 0 ) { - if ( bSig == 0 ) { + if (bExp == 0) { + if (bSig == 0) { invalid: - float_raise( float_flag_invalid, status ); - return floatx80_default_nan(status); - } - normalizeFloatx80Subnormal( bSig, &bExp, &bSig ); + float_raise(float_flag_invalid, status); + return floatx80_default_nan(status); + } + normalizeFloatx80Subnormal(bSig, &bExp, &bSig); } - if ( aExp == 0 ) { -#ifdef SOFTFLOAT_68K - if ( aSig0 == 0 ) { + if (aExp == 0) { + if (aSig0 == 0) { *s = (aSign != bSign); - *q = 0; return a; } -#else - if ( (uint64_t) ( aSig0<<1 ) == 0 ) return a; -#endif - normalizeFloatx80Subnormal( aSig0, &aExp, &aSig0 ); + normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0); } - bSig |= LIT64( 0x8000000000000000 ); + bSig |= LIT64(0x8000000000000000); zSign = aSign; expDiff = aExp - bExp; *s = (aSign != bSign); aSig1 = 0; - if ( expDiff < 0 ) return a; - qTemp = ( bSig <= aSig0 ); - if ( qTemp ) aSig0 -= bSig; - *q = ( expDiff > 63 ) ? 0 : ( qTemp<floatx80_rounding_precision, aSign, aExp, aSig0, 0, status); + qTemp = (bSig <= aSig0); + if (qTemp) aSig0 -= bSig; + *q = (expDiff > 63) ? 0 : (qTemp << expDiff); expDiff -= 64; - while ( 0 < expDiff ) { - qTemp = estimateDiv128To64( aSig0, aSig1, bSig ); - qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0; - mul64To128( bSig, qTemp, &term0, &term1 ); - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - shortShift128Left( aSig0, aSig1, 62, &aSig0, &aSig1 ); - *q = ( expDiff > 63 ) ? 0 : ( qTemp< 63) ? 0 : (qTemp << expDiff); + expDiff -= 62; } expDiff += 64; - if ( 0 < expDiff ) { - qTemp = estimateDiv128To64( aSig0, aSig1, bSig ); - qTemp = ( 2 < qTemp ) ? qTemp - 2 : 0; + if (0 < expDiff) { + qTemp = estimateDiv128To64(aSig0, aSig1, bSig); + qTemp = (2 < qTemp) ? qTemp - 2 : 0; qTemp >>= 64 - expDiff; - mul64To128( bSig, qTemp<<( 64 - expDiff ), &term0, &term1 ); - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); - shortShift128Left( 0, bSig, 64 - expDiff, &term0, &term1 ); - while ( le128( term0, term1, aSig0, aSig1 ) ) { + mul64To128(bSig, qTemp << (64 - expDiff), &term0, &term1); + sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); + shortShift128Left(0, bSig, 64 - expDiff, &term0, &term1); + while (le128(term0, term1, aSig0, aSig1)) { ++qTemp; - sub128( aSig0, aSig1, term0, term1, &aSig0, &aSig1 ); + sub128(aSig0, aSig1, term0, term1, &aSig0, &aSig1); } *q += qTemp; } return normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, - zSign, bExp + expDiff, aSig0, aSig1, status ); - + zSign, bExp + expDiff, aSig0, aSig1, status); + } #endif // end of addition for Previous - /*---------------------------------------------------------------------------- | Returns the square root of the extended double-precision floating-point | value `a'. The operation is performed according to the IEC/IEEE Standard @@ -3244,7 +3236,7 @@ floatx80 floatx80_scale(floatx80 a, floatx80 b, float_status *status) } if ( aExp == 0 ) { if ( aSig == 0 ) return packFloatx80( aSign, 0, 0); - if ( bExp < 0x3FFF ) return a; + if (bExp < 0x3FFF) return normalizeRoundAndPackFloatx80(status->floatx80_rounding_precision, aSign, aExp, aSig, 0, status); normalizeFloatx80Subnormal( aSig, &aExp, &aSig ); } diff --git a/softfloat/softfloat.h b/softfloat/softfloat.h index 727ad8885..5ec0f8125 100644 --- a/softfloat/softfloat.h +++ b/softfloat/softfloat.h @@ -459,6 +459,7 @@ floatx80 floatx80_tan(floatx80 a, float_status *status); floatx80 floatx80_tanh(floatx80 a, float_status *status); floatx80 floatx80_tentox(floatx80 a, float_status *status); floatx80 floatx80_twotox(floatx80 a, float_status *status); +floatx80 floatx80_sincos(floatx80 a, floatx80 *c, float_status *status); #endif // functions originally internal to softfloat.c diff --git a/softfloat/softfloat_fpsp.cpp b/softfloat/softfloat_fpsp.cpp index d28111852..d83d3ce3b 100644 --- a/softfloat/softfloat_fpsp.cpp +++ b/softfloat/softfloat_fpsp.cpp @@ -1574,6 +1574,267 @@ floatx80 floatx80_sin(floatx80 a, float_status *status) } } +/*---------------------------------------------------------------------------- + | Sine and cosine + *----------------------------------------------------------------------------*/ + +floatx80 floatx80_sincos(floatx80 a, floatx80 *c, float_status * status) + { + flag aSign, xSign, rSign, sSign; + int32_t aExp, xExp, rExp, sExp; + uint64_t aSig, rSig, sSig; + + int32_t compact, l, n, i, j1, j2; + floatx80 fp0, fp1, fp2, fp3, fp4, fp5, r, s, invtwopi, twopi1, twopi2; + float32 posneg1, twoto63; + flag endflag; + + aSig = extractFloatx80Frac(a); + aExp = extractFloatx80Exp(a); + aSign = extractFloatx80Sign(a); + + if (aExp == 0x7FFF) { + if ((uint64_t)(aSig << 1)) { + *c = propagateFloatx80NaNOneArg(a, status); + return *c; + + } + float_raise(float_flag_invalid, status); + *c = floatx80_default_nan(status); + return *c; + + } + + if (aExp == 0 && aSig == 0) { + *c = packFloatx80(0, one_exp, one_sig); + return packFloatx80(aSign, 0, 0); + + } + + SET_PREC; + + compact = floatx80_make_compact(aExp, aSig); + + fp0 = a; + + if (compact < 0x3FD78000 || compact > 0x4004BC7E) { // 2^(-40) > |X| > 15 PI + if (compact > 0x3FFF8000) { // |X| >= 15 PI + // REDUCEX + fp1 = packFloatx80(0, 0, 0); + if (compact == 0x7FFEFFFF) { + twopi1 = packFloatx80(aSign ^ 1, 0x7FFE, LIT64(0xC90FDAA200000000)); + twopi2 = packFloatx80(aSign ^ 1, 0x7FDC, LIT64(0x85A308D300000000)); + fp0 = floatx80_add(fp0, twopi1, status); + fp1 = fp0; + fp0 = floatx80_add(fp0, twopi2, status); + fp1 = floatx80_sub(fp1, fp0, status); + fp1 = floatx80_add(fp1, twopi2, status); + + } + loop: + xSign = extractFloatx80Sign(fp0); + xExp = extractFloatx80Exp(fp0); + xExp -= 0x3FFF; + if (xExp <= 28) { + l = 0; + endflag = 1; + } else { + l = xExp - 27; + endflag = 0; + } + invtwopi = packFloatx80(0, 0x3FFE - l, LIT64(0xA2F9836E4E44152A)); // INVTWOPI + twopi1 = packFloatx80(0, 0x3FFF + l, LIT64(0xC90FDAA200000000)); + twopi2 = packFloatx80(0, 0x3FDD + l, LIT64(0x85A308D300000000)); + + twoto63 = 0x5F000000; + twoto63 |= xSign ? 0x80000000 : 0x00000000; // SIGN(INARG)*2^63 IN SGL + + fp2 = floatx80_mul(fp0, invtwopi, status); + fp2 = floatx80_add(fp2, float32_to_floatx80(twoto63, status), status); // THE FRACTIONAL PART OF FP2 IS ROUNDED + fp2 = floatx80_sub(fp2, float32_to_floatx80(twoto63, status), status); // FP2 is N + fp4 = floatx80_mul(twopi1, fp2, status); // W = N*P1 + fp5 = floatx80_mul(twopi2, fp2, status); // w = N*P2 + fp3 = floatx80_add(fp4, fp5, status); // FP3 is P + fp4 = floatx80_sub(fp4, fp3, status); // W-P + fp0 = floatx80_sub(fp0, fp3, status); // FP0 is A := R - P + fp4 = floatx80_add(fp4, fp5, status); // FP4 is p = (W-P)+w + fp3 = fp0; // FP3 is A + fp1 = floatx80_sub(fp1, fp4, status); // FP1 is a := r - p + fp0 = floatx80_add(fp0, fp1, status); // FP0 is R := A+a + + if (endflag > 0) { + n = floatx80_to_int32(fp2, status); + goto sccont; + } + fp3 = floatx80_sub(fp3, fp0, status); // A-R + fp1 = floatx80_add(fp1, fp3, status); // FP1 is r := (A-R)+a + goto loop; + } else { + // SCSM + fp0 = float32_to_floatx80(0x3F800000, status); // 1 + + RESET_PREC; + + // COSTINY + *c = floatx80_sub(fp0, float32_to_floatx80(0x00800000, status), status); + // SINTINY + a = floatx80_move(a, status); + + float_raise(float_flag_inexact, status); + + return a; + + } + } else { + fp1 = floatx80_mul(fp0, float64_to_floatx80(LIT64(0x3FE45F306DC9C883), status), status); // X*2/PI + + n = floatx80_to_int32(fp1, status); + i = 32 + n; + + fp0 = floatx80_sub(fp0, pi_tbl[i], status); // X-Y1 + fp0 = floatx80_sub(fp0, float32_to_floatx80(pi_tbl2[i], status), status); // FP0 IS R = (X-Y1)-Y2 + + sccont: + n &= 3; // k = N mod 4 + if (n & 1) { + // NODD + j1 = n >> 1; // j1 = (k-1)/2 + j2 = j1 ^ (n & 1); // j2 = j1 EOR (k mod 2) + + rSign = extractFloatx80Sign(fp0); // R + rExp = extractFloatx80Exp(fp0); + rSig = extractFloatx80Frac(fp0); + rSign ^= j2; + + fp0 = floatx80_mul(fp0, fp0, status); // FP0 IS S = R*R + fp1 = float64_to_floatx80(LIT64(0xBD6AAA77CCC994F5), status); // A7 + fp2 = float64_to_floatx80(LIT64(0x3D2AC4D0D6011EE3), status); // B8 + fp1 = floatx80_mul(fp1, fp0, status); // FP1 IS SA7 + fp2 = floatx80_mul(fp2, fp0, status); // FP2 IS SB8 + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3DE612097AAE8DA1), status), status); // A6+SA7 + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBDA9396F9F45AC19), status), status); // B7+SB8 + fp1 = floatx80_mul(fp1, fp0, status); // S(A6+SA7) + fp2 = floatx80_mul(fp2, fp0, status); // S(B7+SB8) + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBE5AE6452A118AE4), status), status); // A5+S(A6+SA7) + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3E21EED90612C972), status), status); // B6+S(B7+SB8) + fp1 = floatx80_mul(fp1, fp0, status); // S(A5+S(A6+SA7)) + fp2 = floatx80_mul(fp2, fp0, status); // S(B6+S(B7+SB8)) + + sSign = extractFloatx80Sign(fp0); // S + sExp = extractFloatx80Exp(fp0); + sSig = extractFloatx80Frac(fp0); + sSign ^= j1; + posneg1 = 0x3F800000; + posneg1 |= j1 ? 0x80000000 : 0; + + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3EC71DE3A5341531), status), status); // A4+S(A5+S(A6+SA7)) + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBE927E4FB79D9FCF), status), status); // B5+S(B6+S(B7+SB8)) + fp1 = floatx80_mul(fp1, fp0, status); // S(A4+...) + fp2 = floatx80_mul(fp2, fp0, status); // S(B5+...) + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBF2A01A01A018B59), status), status); // A3+S(A4+...) + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EFA01A01A01D423), status), status); // B4+S(B5+...) + fp1 = floatx80_mul(fp1, fp0, status); // S(A3+...) + fp2 = floatx80_mul(fp2, fp0, status); // S(B4+...) + fp3 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); + fp4 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); + fp1 = floatx80_add(fp1, fp3, status); // A2+S(A3+...) + fp2 = floatx80_add(fp2, fp4, status); // B3+S(B4+...) + fp1 = floatx80_mul(fp1, fp0, status); // S(A2+...) + fp2 = floatx80_mul(fp2, fp0, status); // S(B3+...) + fp3 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); + fp4 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); + fp1 = floatx80_add(fp1, fp3, status); // A1+S(A2+...) + fp2 = floatx80_add(fp2, fp4, status); // B2+S(B3+...) + fp1 = floatx80_mul(fp1, fp0, status); // S(A1+...) + fp0 = floatx80_mul(fp0, fp2, status); // S(B2+...) + + r = packFloatx80(rSign, rExp, rSig); + fp1 = floatx80_mul(fp1, r, status); // R'S(A1+...) + fp0 = floatx80_add(fp0, float32_to_floatx80(0xBF000000, status), status); // B1+S(B2...) + + s = packFloatx80(sSign, sExp, sSig); + fp0 = floatx80_mul(fp0, s, status); // S'(B1+S(B2+...)) + + RESET_PREC; + + *c = floatx80_add(fp1, r, status); // COS(X) + + a = floatx80_add(fp0, float32_to_floatx80(posneg1, status), status); // SIN(X) + + float_raise(float_flag_inexact, status); + + return a; + } else { + // NEVEN + j1 = n >> 1; // j1 = k/2 + + rSign = extractFloatx80Sign(fp0); // R + rExp = extractFloatx80Exp(fp0); + rSig = extractFloatx80Frac(fp0); + rSign ^= j1; + + fp0 = floatx80_mul(fp0, fp0, status); // FP0 IS S = R*R + fp1 = float64_to_floatx80(LIT64(0x3D2AC4D0D6011EE3), status); // B8 + fp2 = float64_to_floatx80(LIT64(0xBD6AAA77CCC994F5), status); // A7 + fp1 = floatx80_mul(fp1, fp0, status); // SB8 + fp2 = floatx80_mul(fp2, fp0, status); // SA7 + + sSign = extractFloatx80Sign(fp0); // S + sExp = extractFloatx80Exp(fp0); + sSig = extractFloatx80Frac(fp0); + sSign ^= j1; + posneg1 = 0x3F800000; + posneg1 |= j1 ? 0x80000000 : 0; + + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBDA9396F9F45AC19), status), status); // B7+SB8 + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3DE612097AAE8DA1), status), status); // A6+SA7 + fp1 = floatx80_mul(fp1, fp0, status); // S(B7+SB8) + fp2 = floatx80_mul(fp2, fp0, status); // S(A6+SA7) + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3E21EED90612C972), status), status); // B6+S(B7+SB8) + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBE5AE6452A118AE4), status), status); // A5+S(A6+SA7) + fp1 = floatx80_mul(fp1, fp0, status); // S(B6+S(B7+SB8)) + fp2 = floatx80_mul(fp2, fp0, status); // S(A5+S(A6+SA7)) + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0xBE927E4FB79D9FCF), status), status); // B5+S(B6+S(B7+SB8)) + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0x3EC71DE3A5341531), status), status); // A4+S(A5+S(A6+SA7)) + fp1 = floatx80_mul(fp1, fp0, status); // S(B5+...) + fp2 = floatx80_mul(fp2, fp0, status); // S(A4+...) + fp1 = floatx80_add(fp1, float64_to_floatx80(LIT64(0x3EFA01A01A01D423), status), status); // B4+S(B5+...) + fp2 = floatx80_add(fp2, float64_to_floatx80(LIT64(0xBF2A01A01A018B59), status), status); // A3+S(A4+...) + fp1 = floatx80_mul(fp1, fp0, status); // S(B4+...) + fp2 = floatx80_mul(fp2, fp0, status); // S(A3+...) + fp3 = packFloatx80(1, 0x3FF5, LIT64(0xB60B60B60B61D438)); + fp4 = packFloatx80(0, 0x3FF8, LIT64(0x88888888888859AF)); + fp1 = floatx80_add(fp1, fp3, status); // B3+S(B4+...) + fp2 = floatx80_add(fp2, fp4, status); // A2+S(A3+...) + fp1 = floatx80_mul(fp1, fp0, status); // S(B3+...) + fp2 = floatx80_mul(fp2, fp0, status); // S(A2+...) + fp3 = packFloatx80(0, 0x3FFA, LIT64(0xAAAAAAAAAAAAAB5E)); + fp4 = packFloatx80(1, 0x3FFC, LIT64(0xAAAAAAAAAAAAAA99)); + fp1 = floatx80_add(fp1, fp3, status); // B2+S(B3+...) + fp2 = floatx80_add(fp2, fp4, status); // A1+S(A2+...) + fp1 = floatx80_mul(fp1, fp0, status); // S(B2+...) + fp0 = floatx80_mul(fp0, fp2, status); // S(A1+...) + fp1 = floatx80_add(fp1, float32_to_floatx80(0xBF000000, status), status); // B1+S(B2...) + + r = packFloatx80(rSign, rExp, rSig); + fp0 = floatx80_mul(fp0, r, status); // R'S(A1+...) + + s = packFloatx80(sSign, sExp, sSig); + fp1 = floatx80_mul(fp1, s, status); // S'(B1+S(B2+...)) + + RESET_PREC; + + *c = floatx80_add(fp1, float32_to_floatx80(posneg1, status), status); // COS(X) + + a = floatx80_add(fp0, r, status); // SIN(X) + + float_raise(float_flag_inexact, status); + + return a; + } + } +} + /*---------------------------------------------------------------------------- | Hyperbolic sine *----------------------------------------------------------------------------*/