summaryrefslogtreecommitdiffstatshomepage
path: root/3rdparty/softfloat
diff options
context:
space:
mode:
Diffstat (limited to '3rdparty/softfloat')
-rw-r--r--3rdparty/softfloat/f2xm1.c270
-rw-r--r--3rdparty/softfloat/fpatan.c389
-rw-r--r--3rdparty/softfloat/fyl2x.c2
-rw-r--r--3rdparty/softfloat/mamesf.h9
-rw-r--r--3rdparty/softfloat/milieu.h6
-rw-r--r--3rdparty/softfloat/softfloat.h3
6 files changed, 663 insertions, 16 deletions
diff --git a/3rdparty/softfloat/f2xm1.c b/3rdparty/softfloat/f2xm1.c
new file mode 100644
index 00000000000..bd5db4471a7
--- /dev/null
+++ b/3rdparty/softfloat/f2xm1.c
@@ -0,0 +1,270 @@
+/*============================================================================
+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
+
+#define USE_estimateDiv128To64
+#include "mamesf.h"
+#include "softfloat.h"
+#include "fpu_constant.h"
+
+#define packFloat_128(zHi, zLo) {(zHi), (zLo)}
+#define PACK_FLOAT_128(hi,lo) packFloat_128(LIT64(hi),LIT64(lo))
+
+static const floatx80 floatx80_negone = packFloatx80(1, 0x3fff, 0x8000000000000000U);
+static const floatx80 floatx80_neghalf = packFloatx80(1, 0x3ffe, 0x8000000000000000U);
+static const float128 float128_ln2 =
+ packFloat_128(0x3ffe62e42fefa39eU, 0xf35793c7673007e6U);
+
+#define LN2_SIG_HI 0xb17217f7d1cf79abU
+#define LN2_SIG_LO 0xc000000000000000U /* 67-bit precision */
+
+#define EXP_ARR_SIZE 15
+
+static float128 exp_arr[EXP_ARR_SIZE] =
+{
+ PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
+ PACK_FLOAT_128(0x3ffe000000000000, 0x0000000000000000), /* 2 */
+ PACK_FLOAT_128(0x3ffc555555555555, 0x5555555555555555), /* 3 */
+ PACK_FLOAT_128(0x3ffa555555555555, 0x5555555555555555), /* 4 */
+ PACK_FLOAT_128(0x3ff8111111111111, 0x1111111111111111), /* 5 */
+ PACK_FLOAT_128(0x3ff56c16c16c16c1, 0x6c16c16c16c16c17), /* 6 */
+ PACK_FLOAT_128(0x3ff2a01a01a01a01, 0xa01a01a01a01a01a), /* 7 */
+ PACK_FLOAT_128(0x3fefa01a01a01a01, 0xa01a01a01a01a01a), /* 8 */
+ PACK_FLOAT_128(0x3fec71de3a556c73, 0x38faac1c88e50017), /* 9 */
+ PACK_FLOAT_128(0x3fe927e4fb7789f5, 0xc72ef016d3ea6679), /* 10 */
+ PACK_FLOAT_128(0x3fe5ae64567f544e, 0x38fe747e4b837dc7), /* 11 */
+ PACK_FLOAT_128(0x3fe21eed8eff8d89, 0x7b544da987acfe85), /* 12 */
+ PACK_FLOAT_128(0x3fde6124613a86d0, 0x97ca38331d23af68), /* 13 */
+ PACK_FLOAT_128(0x3fda93974a8c07c9, 0xd20badf145dfa3e5), /* 14 */
+ PACK_FLOAT_128(0x3fd6ae7f3e733b81, 0xf11d8656b0ee8cb0) /* 15 */
+};
+
+#define EXP_BIAS 0x3FFF
+
+/*----------------------------------------------------------------------------
+| Returns the fraction bits of the extended double-precision floating-point
+| value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE bits64 extractFloatx80Frac( floatx80 a )
+{
+ return a.low;
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the exponent bits of the extended double-precision floating-point
+| value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE int32 extractFloatx80Exp( floatx80 a )
+{
+ return a.high & 0x7FFF;
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the sign bit of the extended double-precision floating-point value
+| `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE flag extractFloatx80Sign( floatx80 a )
+{
+ return a.high>>15;
+
+}
+
+/*----------------------------------------------------------------------------
+| Normalizes the subnormal extended double-precision floating-point value
+| represented by the denormalized significand `aSig'. The normalized exponent
+| and significand are stored at the locations pointed to by `zExpPtr' and
+| `zSigPtr', respectively.
+*----------------------------------------------------------------------------*/
+
+INLINE void normalizeFloatx80Subnormal(uint64_t aSig, int32_t *zExpPtr, uint64_t *zSigPtr)
+{
+ int shiftCount = countLeadingZeros64(aSig);
+ *zSigPtr = aSig<<shiftCount;
+ *zExpPtr = 1 - shiftCount;
+}
+
+
+/*----------------------------------------------------------------------------
+| Takes extended double-precision floating-point NaN `a' and returns the
+| appropriate NaN result. If `a' is a signaling NaN, the invalid exception
+| is raised.
+*----------------------------------------------------------------------------*/
+
+INLINE floatx80 propagateFloatx80NaN(floatx80 a)
+{
+ if (floatx80_is_signaling_nan(a))
+ float_raise(float_flag_invalid);
+
+ a.low |= 0xC000000000000000U;
+
+ return a;
+}
+
+// 2 3 4 n
+// f(x) ~ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+// 0 1 2 3 4 n
+//
+// -- 2k -- 2k+1
+// p(x) = > C * x q(x) = > C * x
+// -- 2k -- 2k+1
+//
+// f(x) ~ [ p(x) + x * q(x) ]
+//
+
+static float128 EvalPoly(float128 x, float128 *arr, unsigned n)
+{
+ float128 x2 = float128_mul(x, x);
+ unsigned i;
+
+ assert(n > 1);
+
+ float128 r1 = arr[--n];
+ i = n;
+ while(i >= 2) {
+ r1 = float128_mul(r1, x2);
+ i -= 2;
+ r1 = float128_add(r1, arr[i]);
+ }
+ if (i) r1 = float128_mul(r1, x);
+
+ float128 r2 = arr[--n];
+ i = n;
+ while(i >= 2) {
+ r2 = float128_mul(r2, x2);
+ i -= 2;
+ r2 = float128_add(r2, arr[i]);
+ }
+ if (i) r2 = float128_mul(r2, x);
+
+ return float128_add(r1, r2);
+}
+
+/* required -1 < x < 1 */
+static float128 poly_exp(float128 x)
+{
+/*
+ // 2 3 4 5 6 7 8 9
+ // x x x x x x x x x
+ // e - 1 ~ x + --- + --- + --- + --- + --- + --- + --- + --- + ...
+ // 2! 3! 4! 5! 6! 7! 8! 9!
+ //
+ // 2 3 4 5 6 7 8
+ // x x x x x x x x
+ // = x [ 1 + --- + --- + --- + --- + --- + --- + --- + --- + ... ]
+ // 2! 3! 4! 5! 6! 7! 8! 9!
+ //
+ // 8 8
+ // -- 2k -- 2k+1
+ // p(x) = > C * x q(x) = > C * x
+ // -- 2k -- 2k+1
+ // k=0 k=0
+ //
+ // x
+ // e - 1 ~ x * [ p(x) + x * q(x) ]
+ //
+*/
+ float128 t = EvalPoly(x, exp_arr, EXP_ARR_SIZE);
+ return float128_mul(t, x);
+}
+
+// =================================================
+// x
+// FX2M1 Compute 2 - 1
+// =================================================
+
+//
+// Uses the following identities:
+//
+// 1. ----------------------------------------------------------
+// x x*ln(2)
+// 2 = e
+//
+// 2. ----------------------------------------------------------
+// 2 3 4 5 n
+// x x x x x x x
+// e = 1 + --- + --- + --- + --- + --- + ... + --- + ...
+// 1! 2! 3! 4! 5! n!
+//
+
+floatx80 f2xm1(floatx80 a)
+{
+ bits64 zSig0, zSig1, zSig2;
+
+ bits64 aSig = extractFloatx80Frac(a);
+ sbits32 aExp = extractFloatx80Exp(a);
+ int aSign = extractFloatx80Sign(a);
+
+ if (aExp == 0x7FFF) {
+ if ((bits64) (aSig<<1))
+ return propagateFloatx80NaN(a);
+
+ return (aSign) ? floatx80_negone : a;
+ }
+
+ if (aExp == 0) {
+ if (aSig == 0) return a;
+ float_raise(float_flag_denormal | float_flag_inexact);
+ normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
+
+ tiny_argument:
+ mul128By64To192(LN2_SIG_HI, LN2_SIG_LO, aSig, &zSig0, &zSig1, &zSig2);
+ if (0 < (sbits64) zSig0) {
+ shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
+ --aExp;
+ }
+ return
+ roundAndPackFloatx80(80, aSign, aExp, zSig0, zSig1);
+ }
+
+ float_raise(float_flag_inexact);
+
+ if (aExp < 0x3FFF)
+ {
+ if (aExp < EXP_BIAS-68)
+ goto tiny_argument;
+
+ /* ******************************** */
+ /* using float128 for approximation */
+ /* ******************************** */
+
+ float128 x = floatx80_to_float128(a);
+ x = float128_mul(x, float128_ln2);
+ x = poly_exp(x);
+ return float128_to_floatx80(x);
+ }
+ else
+ {
+ if (a.high == 0xBFFF && ! (aSig<<1))
+ return floatx80_neghalf;
+
+ return a;
+ }
+}
diff --git a/3rdparty/softfloat/fpatan.c b/3rdparty/softfloat/fpatan.c
new file mode 100644
index 00000000000..db4e1f983d2
--- /dev/null
+++ b/3rdparty/softfloat/fpatan.c
@@ -0,0 +1,389 @@
+/*============================================================================
+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]
+ * Adapted for 3rdparty/softfloat in MAME by Calvin Buckley (05/2021)
+ * ==========================================================================*/
+
+#define FLOAT128
+
+#include "mamesf.h"
+#include "softfloat.h"
+//#include "softfloat-specialize"
+#include "fpu_constant.h"
+
+/* XXX: These are common w/ fsincos/fyl2x; should be moved to common header? */
+#define packFloat_128(zHi, zLo) {(zHi), (zLo)}
+#define packFloat2x128m(zHi, zLo) {(zHi), (zLo)}
+#define PACK_FLOAT_128(hi,lo) packFloat2x128m(LIT64(hi),LIT64(lo))
+#define EXP_BIAS 0x3FFF
+
+/*----------------------------------------------------------------------------
+| Returns the fraction bits of the extended double-precision floating-point
+| value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE bits64 extractFloatx80Frac( floatx80 a )
+{
+ return a.low;
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the exponent bits of the extended double-precision floating-point
+| value `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE int32 extractFloatx80Exp( floatx80 a )
+{
+ return a.high & 0x7FFF;
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the sign bit of the extended double-precision floating-point value
+| `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE flag extractFloatx80Sign( floatx80 a )
+{
+ return a.high>>15;
+
+}
+
+/*----------------------------------------------------------------------------
+| Normalizes the subnormal extended double-precision floating-point value
+| represented by the denormalized significand `aSig'. The normalized exponent
+| and significand are stored at the locations pointed to by `zExpPtr' and
+| `zSigPtr', respectively.
+*----------------------------------------------------------------------------*/
+
+INLINE void normalizeFloatx80Subnormal(uint64_t aSig, int32_t *zExpPtr, uint64_t *zSigPtr)
+{
+ int shiftCount = countLeadingZeros64(aSig);
+ *zSigPtr = aSig<<shiftCount;
+ *zExpPtr = 1 - shiftCount;
+}
+
+/*----------------------------------------------------------------------------
+| Returns 1 if the extended double-precision floating-point value `a' is a
+| NaN; otherwise returns 0.
+*----------------------------------------------------------------------------*/
+
+INLINE int floatx80_is_nan(floatx80 a)
+{
+ return ((a.high & 0x7FFF) == 0x7FFF) && (int64_t) (a.low<<1);
+}
+
+/*----------------------------------------------------------------------------
+| Takes two extended double-precision floating-point values `a' and `b', one
+| of which is a NaN, and returns the appropriate NaN result. If either `a' or
+| `b' is a signaling NaN, the invalid exception is raised.
+*----------------------------------------------------------------------------*/
+
+INLINE floatx80 propagateFloatx80NaN( floatx80 a, floatx80 b )
+{
+ flag aIsNaN, aIsSignalingNaN, bIsNaN, bIsSignalingNaN;
+
+ aIsNaN = floatx80_is_nan( a );
+ aIsSignalingNaN = floatx80_is_signaling_nan( a );
+ bIsNaN = floatx80_is_nan( b );
+ bIsSignalingNaN = floatx80_is_signaling_nan( b );
+ a.low |= LIT64( 0xC000000000000000 );
+ b.low |= LIT64( 0xC000000000000000 );
+ if ( aIsSignalingNaN | bIsSignalingNaN ) float_raise( float_flag_invalid );
+ if ( aIsNaN ) {
+ return ( aIsSignalingNaN & bIsNaN ) ? b : a;
+ }
+ else {
+ return b;
+ }
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the exponent bits of the quadruple-precision floating-point value
+| `a'.
+*----------------------------------------------------------------------------*/
+
+INLINE int32 extractFloat128Exp( float128 a )
+{
+ return ( a.high>>48 ) & 0x7FFF;
+
+}
+
+/* end copied */
+
+#define FPATAN_ARR_SIZE 11
+
+static const float128 float128_one =
+ packFloat_128(0x3fff000000000000U, 0x0000000000000000U);
+static const float128 float128_sqrt3 =
+ packFloat_128(0x3fffbb67ae8584caU, 0xa73b25742d7078b8U);
+static const floatx80 floatx80_one = packFloatx80(0, 0x3fff, 0x8000000000000000U);
+static const floatx80 floatx80_pi =
+ packFloatx80(0, 0x4000, 0xc90fdaa22168c235U);
+
+static const float128 float128_pi2 =
+ packFloat_128(0x3fff921fb54442d1U, 0x8469898CC5170416U);
+static const float128 float128_pi4 =
+ packFloat_128(0x3ffe921fb54442d1U, 0x8469898CC5170416U);
+static const float128 float128_pi6 =
+ packFloat_128(0x3ffe0c152382d736U, 0x58465BB32E0F580FU);
+
+static float128 atan_arr[FPATAN_ARR_SIZE] =
+{
+ PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
+ PACK_FLOAT_128(0xbffd555555555555, 0x5555555555555555), /* 3 */
+ PACK_FLOAT_128(0x3ffc999999999999, 0x999999999999999a), /* 5 */
+ PACK_FLOAT_128(0xbffc249249249249, 0x2492492492492492), /* 7 */
+ PACK_FLOAT_128(0x3ffbc71c71c71c71, 0xc71c71c71c71c71c), /* 9 */
+ PACK_FLOAT_128(0xbffb745d1745d174, 0x5d1745d1745d1746), /* 11 */
+ PACK_FLOAT_128(0x3ffb3b13b13b13b1, 0x3b13b13b13b13b14), /* 13 */
+ PACK_FLOAT_128(0xbffb111111111111, 0x1111111111111111), /* 15 */
+ PACK_FLOAT_128(0x3ffae1e1e1e1e1e1, 0xe1e1e1e1e1e1e1e2), /* 17 */
+ PACK_FLOAT_128(0xbffaaf286bca1af2, 0x86bca1af286bca1b), /* 19 */
+ PACK_FLOAT_128(0x3ffa861861861861, 0x8618618618618618) /* 21 */
+};
+
+extern float128 OddPoly(float128 x, float128 *arr, unsigned n);
+
+/* |x| < 1/4 */
+static float128 poly_atan(float128 x1)
+{
+/*
+ // 3 5 7 9 11 13 15 17
+ // x x x x x x x x
+ // atan(x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- + ----
+ // 3 5 7 9 11 13 15 17
+ //
+ // 2 4 6 8 10 12 14 16
+ // x x x x x x x x
+ // = x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- + ---- ]
+ // 3 5 7 9 11 13 15 17
+ //
+ // 5 5
+ // -- 4k -- 4k+2
+ // p(x) = > C * x q(x) = > C * x
+ // -- 2k -- 2k+1
+ // k=0 k=0
+ //
+ // 2
+ // atan(x) ~ x * [ p(x) + x * q(x) ]
+ //
+*/
+ return OddPoly(x1, atan_arr, FPATAN_ARR_SIZE);
+}
+
+// =================================================
+// FPATAN Compute y * log (x)
+// 2
+// =================================================
+
+//
+// Uses the following identities:
+//
+// 1. ----------------------------------------------------------
+//
+// atan(-x) = -atan(x)
+//
+// 2. ----------------------------------------------------------
+//
+// x + y
+// atan(x) + atan(y) = atan -------, xy < 1
+// 1-xy
+//
+// x + y
+// atan(x) + atan(y) = atan ------- + PI, x > 0, xy > 1
+// 1-xy
+//
+// x + y
+// atan(x) + atan(y) = atan ------- - PI, x < 0, xy > 1
+// 1-xy
+//
+// 3. ----------------------------------------------------------
+//
+// atan(x) = atan(INF) + atan(- 1/x)
+//
+// x-1
+// atan(x) = PI/4 + atan( ----- )
+// x+1
+//
+// x * sqrt(3) - 1
+// atan(x) = PI/6 + atan( ----------------- )
+// x + sqrt(3)
+//
+// 4. ----------------------------------------------------------
+// 3 5 7 9 2n+1
+// x x x x n x
+// atan(x) = x - --- + --- - --- + --- - ... + (-1) ------ + ...
+// 3 5 7 9 2n+1
+//
+
+floatx80 floatx80_fpatan(floatx80 a, floatx80 b)
+{
+ uint64_t aSig = extractFloatx80Frac(a);
+ int32_t aExp = extractFloatx80Exp(a);
+ int aSign = extractFloatx80Sign(a);
+ uint64_t bSig = extractFloatx80Frac(b);
+ int32_t bExp = extractFloatx80Exp(b);
+ int bSign = extractFloatx80Sign(b);
+
+ int zSign = aSign ^ bSign;
+
+ if (bExp == 0x7FFF)
+ {
+ if ((uint64_t) (bSig<<1))
+ return propagateFloatx80NaN(a, b);
+
+ if (aExp == 0x7FFF) {
+ if ((uint64_t) (aSig<<1))
+ return propagateFloatx80NaN(a, b);
+
+ if (aSign) { /* return 3PI/4 */
+ return roundAndPackFloatx80(80, bSign,
+ FLOATX80_3PI4_EXP, FLOAT_3PI4_HI, FLOAT_3PI4_LO);
+ }
+ else { /* return PI/4 */
+ return roundAndPackFloatx80(80, bSign,
+ FLOATX80_PI4_EXP, FLOAT_PI_HI, FLOAT_PI_LO);
+ }
+ }
+
+ if (aSig && (aExp == 0))
+ float_raise(float_flag_denormal);
+
+ /* return PI/2 */
+ return roundAndPackFloatx80(80, bSign, FLOATX80_PI2_EXP, FLOAT_PI_HI, FLOAT_PI_LO);
+ }
+ if (aExp == 0x7FFF)
+ {
+ if ((uint64_t) (aSig<<1))
+ return propagateFloatx80NaN(a, b);
+
+ if (bSig && (bExp == 0))
+ float_raise(float_flag_denormal);
+
+return_PI_or_ZERO:
+
+ if (aSign) { /* return PI */
+ return roundAndPackFloatx80(80, bSign, FLOATX80_PI_EXP, FLOAT_PI_HI, FLOAT_PI_LO);
+ } else { /* return 0 */
+ return packFloatx80(bSign, 0, 0);
+ }
+ }
+ if (bExp == 0)
+ {
+ if (bSig == 0) {
+ if (aSig && (aExp == 0)) float_raise(float_flag_denormal);
+ goto return_PI_or_ZERO;
+ }
+
+ float_raise(float_flag_denormal);
+ normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
+ }
+ if (aExp == 0)
+ {
+ if (aSig == 0) /* return PI/2 */
+ return roundAndPackFloatx80(80, bSign, FLOATX80_PI2_EXP, FLOAT_PI_HI, FLOAT_PI_LO);
+
+ float_raise(float_flag_denormal);
+ normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
+ }
+
+ float_raise(float_flag_inexact);
+
+ /* |a| = |b| ==> return PI/4 */
+ if (aSig == bSig && aExp == bExp)
+ return roundAndPackFloatx80(80, bSign, FLOATX80_PI4_EXP, FLOAT_PI_HI, FLOAT_PI_LO);
+
+ /* ******************************** */
+ /* using float128 for approximation */
+ /* ******************************** */
+
+ float128 a128 = normalizeRoundAndPackFloat128(0, aExp-0x10, aSig, 0);
+ float128 b128 = normalizeRoundAndPackFloat128(0, bExp-0x10, bSig, 0);
+ float128 x;
+ int swap = 0, add_pi6 = 0, add_pi4 = 0;
+
+ if (aExp > bExp || (aExp == bExp && aSig > bSig))
+ {
+ x = float128_div(b128, a128);
+ }
+ else {
+ x = float128_div(a128, b128);
+ swap = 1;
+ }
+
+ int32_t xExp = extractFloat128Exp(x);
+
+ if (xExp <= EXP_BIAS-40)
+ goto approximation_completed;
+
+ if (x.high >= 0x3ffe800000000000U) // 3/4 < x < 1
+ {
+ /*
+ arctan(x) = arctan((x-1)/(x+1)) + pi/4
+ */
+ float128 t1 = float128_sub(x, float128_one);
+ float128 t2 = float128_add(x, float128_one);
+ x = float128_div(t1, t2);
+ add_pi4 = 1;
+ }
+ else
+ {
+ /* argument correction */
+ if (xExp >= 0x3FFD) // 1/4 < x < 3/4
+ {
+ /*
+ arctan(x) = arctan((x*sqrt(3)-1)/(x+sqrt(3))) + pi/6
+ */
+ float128 t1 = float128_mul(x, float128_sqrt3);
+ float128 t2 = float128_add(x, float128_sqrt3);
+ x = float128_sub(t1, float128_one);
+ x = float128_div(x, t2);
+ add_pi6 = 1;
+ }
+ }
+
+ x = poly_atan(x);
+ if (add_pi6) x = float128_add(x, float128_pi6);
+ if (add_pi4) x = float128_add(x, float128_pi4);
+
+approximation_completed:
+ if (swap) x = float128_sub(float128_pi2, x);
+ floatx80 result = float128_to_floatx80(x);
+ if (zSign) result = floatx80_chs(result);
+ int rSign = extractFloatx80Sign(result);
+ if (!bSign && rSign)
+ return floatx80_add(result, floatx80_pi);
+ if (bSign && !rSign)
+ return floatx80_sub(result, floatx80_pi);
+ return result;
+}
+
+// The former function maps to x87 FPATAN, but we can simulate 68881 FATAN with
+// it by simply hardcoding one here.
+floatx80 floatx80_fatan(floatx80 a)
+{
+ return floatx80_fpatan(a, floatx80_one);
+}
diff --git a/3rdparty/softfloat/fyl2x.c b/3rdparty/softfloat/fyl2x.c
index 44c3e7610bc..ed75e3e569c 100644
--- a/3rdparty/softfloat/fyl2x.c
+++ b/3rdparty/softfloat/fyl2x.c
@@ -251,7 +251,7 @@ static float128 poly_l2p1(float128 x)
// 1-u 3 5 7 2n+1
//
-static floatx80 fyl2x(floatx80 a, floatx80 b)
+floatx80 fyl2x(floatx80 a, floatx80 b)
{
uint64_t aSig = extractFloatx80Frac(a);
int32_t aExp = extractFloatx80Exp(a);
diff --git a/3rdparty/softfloat/mamesf.h b/3rdparty/softfloat/mamesf.h
index 437f54834e3..08b03d55d75 100644
--- a/3rdparty/softfloat/mamesf.h
+++ b/3rdparty/softfloat/mamesf.h
@@ -1,13 +1,4 @@
/*----------------------------------------------------------------------------
-| One of the macros `BIGENDIAN' or `LITTLEENDIAN' must be defined.
-*----------------------------------------------------------------------------*/
-#ifdef LSB_FIRST
-#define LITTLEENDIAN
-#else
-#define BIGENDIAN
-#endif
-
-/*----------------------------------------------------------------------------
| The macro `BITS64' can be defined to indicate that 64-bit integer types are
| supported by the compiler.
*----------------------------------------------------------------------------*/
diff --git a/3rdparty/softfloat/milieu.h b/3rdparty/softfloat/milieu.h
index 10687b755d0..523710d27f1 100644
--- a/3rdparty/softfloat/milieu.h
+++ b/3rdparty/softfloat/milieu.h
@@ -34,9 +34,3 @@ these four paragraphs for those parts of this code that are retained.
| Include common integer types and flags.
*----------------------------------------------------------------------------*/
#include "mamesf.h"
-
-/*----------------------------------------------------------------------------
-| Symbolic Boolean literals.
-*----------------------------------------------------------------------------*/
-#define FALSE 0
-#define TRUE 1
diff --git a/3rdparty/softfloat/softfloat.h b/3rdparty/softfloat/softfloat.h
index 6f4389b10bc..a97408c3222 100644
--- a/3rdparty/softfloat/softfloat.h
+++ b/3rdparty/softfloat/softfloat.h
@@ -244,6 +244,9 @@ int floatx80_fsin(floatx80 &a);
int floatx80_fcos(floatx80 &a);
int floatx80_ftan(floatx80 &a);
+floatx80 floatx80_fpatan(floatx80 a, floatx80 b);
+floatx80 floatx80_fatan(floatx80 a);
+
floatx80 floatx80_flognp1(floatx80 a);
floatx80 floatx80_flogn(floatx80 a);
floatx80 floatx80_flog2(floatx80 a);