diff options
Diffstat (limited to '3rdparty/softfloat3/bochs_ext/fyl2x.c')
-rw-r--r-- | 3rdparty/softfloat3/bochs_ext/fyl2x.c | 409 |
1 files changed, 409 insertions, 0 deletions
diff --git a/3rdparty/softfloat3/bochs_ext/fyl2x.c b/3rdparty/softfloat3/bochs_ext/fyl2x.c new file mode 100644 index 00000000000..5b9f21bbf11 --- /dev/null +++ b/3rdparty/softfloat3/bochs_ext/fyl2x.c @@ -0,0 +1,409 @@ +/*============================================================================ +This source file is an extension to the SoftFloat IEC/IEEE Floating-point +Arithmetic Package, Release 2b, written for Bochs (x86 achitecture simulator) +floating point emulation. + +THIS SOFTWARE IS DISTRIBUTED AS IS, FOR FREE. Although reasonable effort has +been made to avoid it, THIS SOFTWARE MAY CONTAIN FAULTS THAT WILL AT TIMES +RESULT IN INCORRECT BEHAVIOR. USE OF THIS SOFTWARE IS RESTRICTED TO PERSONS +AND ORGANIZATIONS WHO CAN AND WILL TAKE FULL RESPONSIBILITY FOR ALL LOSSES, +COSTS, OR OTHER PROBLEMS THEY INCUR DUE TO THE SOFTWARE, AND WHO FURTHERMORE +EFFECTIVELY INDEMNIFY JOHN HAUSER AND THE INTERNATIONAL COMPUTER SCIENCE +INSTITUTE (possibly via similar legal warning) AGAINST ALL LOSSES, COSTS, OR +OTHER PROBLEMS INCURRED BY THEIR CUSTOMERS AND CLIENTS DUE TO THE SOFTWARE. + +Derivative works are acceptable, even for commercial purposes, so long as +(1) the source code for the derivative work includes prominent notice that +the work is derivative, and (2) the source code includes prominent notice with +these four paragraphs for those parts of this code that are retained. +=============================================================================*/ + +/*============================================================================ + * Written for Bochs (x86 achitecture simulator) by + * Stanislav Shwartsman [sshwarts at sourceforge net] + * ==========================================================================*/ + +#define FLOAT128 + +#include "../build/MAME/platform.h" +#include "../source/include/internals.h" +#include "../source/include/softfloat.h" + +#include "fpu_constant.h" +#include "softfloat-extra.h" +#include "softfloat-helpers.h" +#include "softfloat-specialize.h" + +static const extFloat80_t floatx80_one = packFloatx80(0, 0x3fff, uint64_t(0x8000000000000000)); +static const extFloat80_t floatx80_ln_2 = packFloatx80(0, 0x3ffe, 0xb17217f7d1cf79acU); +static const extFloat80_t floatx80_log10_2 = packFloatx80(0, 0x3ffd, 0x9a209a84fbcff798U); + +static const float128_t float128_one = + packFloat128(uint64_t(0x3fff000000000000), uint64_t(0x0000000000000000)); +static const float128_t float128_two = + packFloat128(uint64_t(0x4000000000000000), uint64_t(0x0000000000000000)); + +static const float128_t float128_ln2inv2 = + packFloat128(uint64_t(0x400071547652b82f), uint64_t(0xe1777d0ffda0d23a)); + +#define SQRT2_HALF_SIG uint64_t(0xb504f333f9de6484) + +extern float128_t OddPoly(float128_t x, const float128_t *arr, int n); + +#define L2_ARR_SIZE 9 + +static float128_t ln_arr[L2_ARR_SIZE] = +{ + PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */ + PACK_FLOAT_128(0x3ffd555555555555, 0x5555555555555555), /* 3 */ + PACK_FLOAT_128(0x3ffc999999999999, 0x999999999999999a), /* 5 */ + PACK_FLOAT_128(0x3ffc249249249249, 0x2492492492492492), /* 7 */ + PACK_FLOAT_128(0x3ffbc71c71c71c71, 0xc71c71c71c71c71c), /* 9 */ + PACK_FLOAT_128(0x3ffb745d1745d174, 0x5d1745d1745d1746), /* 11 */ + PACK_FLOAT_128(0x3ffb3b13b13b13b1, 0x3b13b13b13b13b14), /* 13 */ + PACK_FLOAT_128(0x3ffb111111111111, 0x1111111111111111), /* 15 */ + PACK_FLOAT_128(0x3ffae1e1e1e1e1e1, 0xe1e1e1e1e1e1e1e2) /* 17 */ +}; + +static float128_t poly_ln(float128_t x1) +{ +/* + // + // 3 5 7 9 11 13 15 + // 1+u u u u u u u u + // 1/2 ln --- ~ u + --- + --- + --- + --- + ---- + ---- + ---- = + // 1-u 3 5 7 9 11 13 15 + // + // 2 4 6 8 10 12 14 + // u u u u u u u + // = u * [ 1 + --- + --- + --- + --- + ---- + ---- + ---- ] = + // 3 5 7 9 11 13 15 + // + // 3 3 + // -- 4k -- 4k+2 + // p(u) = > C * u q(u) = > C * u + // -- 2k -- 2k+1 + // k=0 k=0 + // + // 1+u 2 + // 1/2 ln --- ~ u * [ p(u) + u * q(u) ] + // 1-u + // +*/ + return OddPoly(x1, (const float128_t*) ln_arr, L2_ARR_SIZE); +} + +/* required sqrt(2)/2 < x < sqrt(2) */ +static float128_t poly_l2(float128_t x) +{ + /* using float128 for approximation */ + float128_t x_p1 = f128_add(x, float128_one); + float128_t x_m1 = f128_sub(x, float128_one); + x = f128_div(x_m1, x_p1); + x = poly_ln(x); + x = f128_mul(x, float128_ln2inv2); + return x; +} + +static float128_t poly_l2p1(float128_t x) +{ + /* using float128 for approximation */ + float128_t x_plus2 = f128_add(x, float128_two); + x = f128_div(x, x_plus2); + x = poly_ln(x); + x = f128_mul(x, float128_ln2inv2); + return x; +} + +// ================================================= +// FYL2X Compute y * log (x) +// 2 +// ================================================= + +// +// Uses the following identities: +// +// 1. ---------------------------------------------------------- +// ln(x) +// log (x) = -------, ln (x*y) = ln(x) + ln(y) +// 2 ln(2) +// +// 2. ---------------------------------------------------------- +// 1+u x-1 +// ln (x) = ln -----, when u = ----- +// 1-u x+1 +// +// 3. ---------------------------------------------------------- +// 3 5 7 2n+1 +// 1+u u u u u +// ln ----- = 2 [ u + --- + --- + --- + ... + ------ + ... ] +// 1-u 3 5 7 2n+1 +// + +extFloat80_t extFloat80_fyl2x(extFloat80_t a, extFloat80_t b) +{ + // handle unsupported extended double-precision floating encodings + if (extF80_isUnsupported(a) || extF80_isUnsupported(b)) { +invalid: + softfloat_exceptionFlags |= softfloat_flag_invalid; + return floatx80_default_nan; + } + + uint64_t aSig = extF80_fraction(a); + int32_t aExp = extF80_exp(a); + int aSign = extF80_sign(a); + uint64_t bSig = extF80_fraction(b); + int32_t bExp = extF80_exp(b); + int bSign = extF80_sign(b); + + int zSign = bSign ^ 1; + + if (aExp == 0x7FFF) { + if ((aSig<<1) || ((bExp == 0x7FFF) && (bSig<<1))) { + const uint128 nan = softfloat_propagateNaNExtF80UI(a.signExp, aSig, b.signExp, bSig); + extFloat80_t rv; + rv.signExp = nan.v64; + rv.signif = nan.v0; + return rv; + } + if (aSign) goto invalid; + else { + if (! bExp) { + if (! bSig) goto invalid; + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + } + return packFloatx80(bSign, 0x7FFF, uint64_t(0x8000000000000000)); + } + } + if (bExp == 0x7FFF) { + if (bSig << 1) { + const uint128 nan = softfloat_propagateNaNExtF80UI(a.signExp, aSig, b.signExp, bSig); + extFloat80_t rv; + rv.signExp = nan.v64; + rv.signif = nan.v0; + return rv; + } + if (aSign && (uint64_t)(aExp | aSig)) goto invalid; + if (aSig && ! aExp) + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + if (aExp < 0x3FFF) + { + return packFloatx80(zSign, 0x7FFF, uint64_t(0x8000000000000000)); + } + if (aExp == 0x3FFF && ! (aSig<<1)) goto invalid; + return packFloatx80(bSign, 0x7FFF, uint64_t(0x8000000000000000)); + } + if (! aExp) { + if (! aSig) { + if ((bExp | bSig) == 0) goto invalid; + softfloat_exceptionFlags |= softfloat_flag_invalid; // divide by zero + return packFloatx80(zSign, 0x7FFF, uint64_t(0x8000000000000000)); + } + if (aSign) goto invalid; + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + struct exp32_sig64 normExpSig = softfloat_normSubnormalExtF80Sig(aSig); + aExp = normExpSig.exp + 1; + aSig = normExpSig.sig; + } + if (aSign) goto invalid; + if (! bExp) { + if (! bSig) { + if (aExp < 0x3FFF) return packFloatx80(zSign, 0, 0); + return packFloatx80(bSign, 0, 0); + } + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + struct exp32_sig64 normExpSig = softfloat_normSubnormalExtF80Sig(bSig); + bExp = normExpSig.exp + 1; + bSig = normExpSig.sig; + } + if (aExp == 0x3FFF && ! (aSig<<1)) + return packFloatx80(bSign, 0, 0); + + softfloat_exceptionFlags |= softfloat_flag_inexact; + + int ExpDiff = aExp - 0x3FFF; + aExp = 0; + if (aSig >= SQRT2_HALF_SIG) { + ExpDiff++; + aExp--; + } + + /* ******************************** */ + /* using float128 for approximation */ + /* ******************************** */ + + float128_t b128 = softfloat_normRoundPackToF128(bSign, bExp-0x10, bSig, 0); + + uint64_t zSig0, zSig1; + shortShift128Right(aSig<<1, 0, 16, &zSig0, &zSig1); + float128_t x = packFloat128(0, aExp+0x3FFF, zSig0, zSig1); + x = poly_l2(x); + x = f128_add(x, i32_to_f128(ExpDiff)); + x = f128_mul(x, b128); + return f128_to_extF80(x); +} + +// ================================================= +// FYL2XP1 Compute y * log (x + 1) +// 2 +// ================================================= + +// +// Uses the following identities: +// +// 1. ---------------------------------------------------------- +// ln(x) +// log (x) = ------- +// 2 ln(2) +// +// 2. ---------------------------------------------------------- +// 1+u x +// ln (x+1) = ln -----, when u = ----- +// 1-u x+2 +// +// 3. ---------------------------------------------------------- +// 3 5 7 2n+1 +// 1+u u u u u +// ln ----- = 2 [ u + --- + --- + --- + ... + ------ + ... ] +// 1-u 3 5 7 2n+1 +// + +extFloat80_t extFloat80_fyl2xp1(extFloat80_t a, extFloat80_t b) +{ + int32_t aExp, bExp; + uint64_t aSig, bSig, zSig0, zSig1, zSig2; + int aSign, bSign; + + // handle unsupported extended double-precision floating encodings + if (extF80_isUnsupported(a) || extF80_isUnsupported(b)) { +invalid: + softfloat_exceptionFlags |= softfloat_flag_invalid; + return floatx80_default_nan; + } + + + aSig = extF80_fraction(a); + aExp = extF80_exp(a); + aSign = extF80_sign(a); + bSig = extF80_fraction(b); + bExp = extF80_exp(b); + bSign = extF80_sign(b); + int zSign = aSign ^ bSign; + + if (aExp == 0x7FFF) { + if ((aSig<<1) != 0 || ((bExp == 0x7FFF) && (bSig<<1) != 0)) { + const uint128 nan = softfloat_propagateNaNExtF80UI(a.signExp, aSig, b.signExp, bSig); + extFloat80_t rv; + rv.signExp = nan.v64; + rv.signif = nan.v0; + return rv; + } + if (aSign) goto invalid; + else { + if (! bExp) { + if (! bSig) goto invalid; + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + } + return packFloatx80(bSign, 0x7FFF, uint64_t(0x8000000000000000)); + } + } + if (bExp == 0x7FFF) + { + if (bSig << 1) { + const uint128 nan = softfloat_propagateNaNExtF80UI(a.signExp, aSig, b.signExp, bSig); + extFloat80_t rv; + rv.signExp = nan.v64; + rv.signif = nan.v0; + return rv; + } + if (! aExp) { + if (! aSig) goto invalid; + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + } + + return packFloatx80(zSign, 0x7FFF, uint64_t(0x8000000000000000)); + } + if (! aExp) { + if (! aSig) { + if (bSig && ! bExp) + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + + return packFloatx80(zSign, 0, 0); + } + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + struct exp32_sig64 normExpSig = softfloat_normSubnormalExtF80Sig(aSig); + aExp = normExpSig.exp + 1; + aSig = normExpSig.sig; + } + if (! bExp) { + if (! bSig) return packFloatx80(zSign, 0, 0); + softfloat_exceptionFlags |= softfloat_flag_invalid; // actually denormal + struct exp32_sig64 normExpSig = softfloat_normSubnormalExtF80Sig(bSig); + bExp = normExpSig.exp + 1; + bSig = normExpSig.sig; + } + + softfloat_exceptionFlags |= softfloat_flag_inexact; + + if (aSign && aExp >= 0x3FFF) + return a; + + if (aExp >= 0x3FFC) // big argument + { + return extFloat80_fyl2x(extF80_add(a, floatx80_one), b); + } + + // handle tiny argument + if (aExp < FLOATX80_EXP_BIAS-70) + { + // first order approximation, return (a*b)/ln(2) + int32_t zExp = aExp + FLOAT_LN2INV_EXP - 0x3FFE; + + mul128By64To192(FLOAT_LN2INV_HI, FLOAT_LN2INV_LO, aSig, &zSig0, &zSig1, &zSig2); + if (0 < (int64_t) zSig0) { + shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1); + --zExp; + } + + zExp = zExp + bExp - 0x3FFE; + mul128By64To192(zSig0, zSig1, bSig, &zSig0, &zSig1, &zSig2); + if (0 < (int64_t) zSig0) { + shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1); + --zExp; + } + + return softfloat_roundPackToExtF80(aSign ^ bSign, zExp, zSig0, zSig1, 80); + } + + /* ******************************** */ + /* using float128 for approximation */ + /* ******************************** */ + + float128_t b128 = softfloat_normRoundPackToF128(bSign, bExp-0x10, bSig, 0); + + shortShift128Right(aSig<<1, 0, 16, &zSig0, &zSig1); + float128_t x = packFloat128(aSign, aExp, zSig0, zSig1); + x = poly_l2p1(x); + x = f128_mul(x, b128); + return f128_to_extF80(x); +} + +extFloat80_t extFloat80_lognp1(extFloat80_t a) +{ + return extFloat80_fyl2x(a, floatx80_ln_2); +} + +extFloat80_t extFloat80_logn(extFloat80_t a) +{ + return extFloat80_fyl2x(a, floatx80_ln_2); +} + +extFloat80_t extFloat80_log2(extFloat80_t a) +{ + return extFloat80_fyl2x(a, floatx80_one); +} + +extFloat80_t extFloat80_log10(extFloat80_t a) +{ + return extFloat80_fyl2x(a, floatx80_log10_2); +} |