--- /dev/null
+/*============================================================================
+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 "m68kcpu.h" // which includes softfloat.h after defining the basic types
+#include "assert.h"
+
+//#include "softfloat-specialize"
+#include "fpu_constant.h"
+
+#define floatx80_one packFloatx80(0, 0x3fff, 0x8000000000000000U)
+#define floatx80_default_nan packFloatx80(0, 0xffff, 0xffffffffffffffffU)
+
+#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'.
+*----------------------------------------------------------------------------*/
+
+static inline bits64 extractFloatx80Frac( floatx80 a )
+{
+ return a.low;
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the exponent bits of the extended double-precision floating-point
+| value `a'.
+*----------------------------------------------------------------------------*/
+
+static inline int32 extractFloatx80Exp( floatx80 a )
+{
+ return a.high & 0x7FFF;
+
+}
+
+/*----------------------------------------------------------------------------
+| Returns the sign bit of the extended double-precision floating-point value
+| `a'.
+*----------------------------------------------------------------------------*/
+
+static inline flag extractFloatx80Sign( floatx80 a )
+{
+ return a.high>>15;
+
+}
+
+/*----------------------------------------------------------------------------
+| 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.
+*----------------------------------------------------------------------------*/
+
+static inline floatx80 propagateFloatx80NaNOneArg(floatx80 a)
+{
+ if (floatx80_is_signaling_nan(a))
+ float_raise(float_flag_invalid);
+
+ a.low |= 0xC000000000000000U;
+
+ return a;
+}
+
+/*----------------------------------------------------------------------------
+| 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.
+*----------------------------------------------------------------------------*/
+
+void normalizeFloatx80Subnormal(uint64_t aSig, int32_t *zExpPtr, uint64_t *zSigPtr)
+{
+ int shiftCount = countLeadingZeros64(aSig);
+ *zSigPtr = aSig<<shiftCount;
+ *zExpPtr = 1 - shiftCount;
+}
+
+/* reduce trigonometric function argument using 128-bit precision
+ M_PI approximation */
+static uint64_t argument_reduction_kernel(uint64_t aSig0, int Exp, uint64_t *zSig0, uint64_t *zSig1)
+{
+ uint64_t term0, term1, term2;
+ uint64_t aSig1 = 0;
+
+ shortShift128Left(aSig1, aSig0, Exp, &aSig1, &aSig0);
+ uint64_t q = estimateDiv128To64(aSig1, aSig0, FLOAT_PI_HI);
+ mul128By64To192(FLOAT_PI_HI, FLOAT_PI_LO, q, &term0, &term1, &term2);
+ sub128(aSig1, aSig0, term0, term1, zSig1, zSig0);
+ while ((int64_t)(*zSig1) < 0) {
+ --q;
+ add192(*zSig1, *zSig0, term2, 0, FLOAT_PI_HI, FLOAT_PI_LO, zSig1, zSig0, &term2);
+ }
+ *zSig1 = term2;
+ return q;
+}
+
+static int reduce_trig_arg(int expDiff, int zSign, uint64_t aSig0, uint64_t aSig1)
+{
+ uint64_t term0, term1, q = 0;
+
+ if (expDiff < 0) {
+ shift128Right(aSig0, 0, 1, &aSig0, &aSig1);
+ expDiff = 0;
+ }
+ if (expDiff > 0) {
+ q = argument_reduction_kernel(aSig0, expDiff, &aSig0, &aSig1);
+ }
+ else {
+ if (FLOAT_PI_HI <= aSig0) {
+ aSig0 -= FLOAT_PI_HI;
+ q = 1;
+ }
+ }
+
+ shift128Right(FLOAT_PI_HI, FLOAT_PI_LO, 1, &term0, &term1);
+ if (! lt128(aSig0, aSig1, term0, term1))
+ {
+ int lt = lt128(term0, term1, aSig0, aSig1);
+ int eq = eq128(aSig0, aSig1, term0, term1);
+
+ if ((eq && (q & 1)) || lt) {
+ zSign = !zSign;
+ ++q;
+ }
+ if (lt) sub128(FLOAT_PI_HI, FLOAT_PI_LO, aSig0, aSig1, &aSig0, &aSig1);
+ }
+
+ return (int)(q & 3);
+}
+
+#define SIN_ARR_SIZE 11
+#define COS_ARR_SIZE 11
+
+static float128 sin_arr[SIN_ARR_SIZE] =
+{
+ PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 1 */
+ PACK_FLOAT_128(0xbffc555555555555, 0x5555555555555555), /* 3 */
+ PACK_FLOAT_128(0x3ff8111111111111, 0x1111111111111111), /* 5 */
+ PACK_FLOAT_128(0xbff2a01a01a01a01, 0xa01a01a01a01a01a), /* 7 */
+ PACK_FLOAT_128(0x3fec71de3a556c73, 0x38faac1c88e50017), /* 9 */
+ PACK_FLOAT_128(0xbfe5ae64567f544e, 0x38fe747e4b837dc7), /* 11 */
+ PACK_FLOAT_128(0x3fde6124613a86d0, 0x97ca38331d23af68), /* 13 */
+ PACK_FLOAT_128(0xbfd6ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 15 */
+ PACK_FLOAT_128(0x3fce952c77030ad4, 0xa6b2605197771b00), /* 17 */
+ PACK_FLOAT_128(0xbfc62f49b4681415, 0x724ca1ec3b7b9675), /* 19 */
+ PACK_FLOAT_128(0x3fbd71b8ef6dcf57, 0x18bef146fcee6e45) /* 21 */
+};
+
+static float128 cos_arr[COS_ARR_SIZE] =
+{
+ PACK_FLOAT_128(0x3fff000000000000, 0x0000000000000000), /* 0 */
+ PACK_FLOAT_128(0xbffe000000000000, 0x0000000000000000), /* 2 */
+ PACK_FLOAT_128(0x3ffa555555555555, 0x5555555555555555), /* 4 */
+ PACK_FLOAT_128(0xbff56c16c16c16c1, 0x6c16c16c16c16c17), /* 6 */
+ PACK_FLOAT_128(0x3fefa01a01a01a01, 0xa01a01a01a01a01a), /* 8 */
+ PACK_FLOAT_128(0xbfe927e4fb7789f5, 0xc72ef016d3ea6679), /* 10 */
+ PACK_FLOAT_128(0x3fe21eed8eff8d89, 0x7b544da987acfe85), /* 12 */
+ PACK_FLOAT_128(0xbfda93974a8c07c9, 0xd20badf145dfa3e5), /* 14 */
+ PACK_FLOAT_128(0x3fd2ae7f3e733b81, 0xf11d8656b0ee8cb0), /* 16 */
+ PACK_FLOAT_128(0xbfca6827863b97d9, 0x77bb004886a2c2ab), /* 18 */
+ PACK_FLOAT_128(0x3fc1e542ba402022, 0x507a9cad2bf8f0bb) /* 20 */
+};
+
+extern float128 OddPoly (float128 x, float128 *arr, unsigned n);
+
+/* 0 <= x <= pi/4 */
+static inline float128 poly_sin(float128 x)
+{
+ // 3 5 7 9 11 13 15
+ // x x x x x x x
+ // sin (x) ~ x - --- + --- - --- + --- - ---- + ---- - ---- =
+ // 3! 5! 7! 9! 11! 13! 15!
+ //
+ // 2 4 6 8 10 12 14
+ // x x x x x x x
+ // = x * [ 1 - --- + --- - --- + --- - ---- + ---- - ---- ] =
+ // 3! 5! 7! 9! 11! 13! 15!
+ //
+ // 3 3
+ // -- 4k -- 4k+2
+ // p(x) = > C * x > 0 q(x) = > C * x < 0
+ // -- 2k -- 2k+1
+ // k=0 k=0
+ //
+ // 2
+ // sin(x) ~ x * [ p(x) + x * q(x) ]
+ //
+
+ return OddPoly(x, sin_arr, SIN_ARR_SIZE);
+}
+
+extern float128 EvenPoly(float128 x, float128 *arr, unsigned n);
+
+/* 0 <= x <= pi/4 */
+static inline float128 poly_cos(float128 x)
+{
+ // 2 4 6 8 10 12 14
+ // x x x x x x x
+ // cos (x) ~ 1 - --- + --- - --- + --- - ---- + ---- - ----
+ // 2! 4! 6! 8! 10! 12! 14!
+ //
+ // 3 3
+ // -- 4k -- 4k+2
+ // p(x) = > C * x > 0 q(x) = > C * x < 0
+ // -- 2k -- 2k+1
+ // k=0 k=0
+ //
+ // 2
+ // cos(x) ~ [ p(x) + x * q(x) ]
+ //
+
+ return EvenPoly(x, cos_arr, COS_ARR_SIZE);
+}
+
+static inline void sincos_invalid(floatx80 *sin_a, floatx80 *cos_a, floatx80 a)
+{
+ if (sin_a) *sin_a = a;
+ if (cos_a) *cos_a = a;
+}
+
+static inline void sincos_tiny_argument(floatx80 *sin_a, floatx80 *cos_a, floatx80 a)
+{
+ if (sin_a) *sin_a = a;
+ if (cos_a) *cos_a = floatx80_one;
+}
+
+static floatx80 sincos_approximation(int neg, float128 r, uint64_t quotient)
+{
+ if (quotient & 0x1) {
+ r = poly_cos(r);
+ neg = 0;
+ } else {
+ r = poly_sin(r);
+ }
+
+ floatx80 result = float128_to_floatx80(r);
+ if (quotient & 0x2)
+ neg = ! neg;
+
+ if (neg)
+ result = floatx80_chs(result);
+
+ return result;
+}
+
+// =================================================
+// SFFSINCOS Compute sin(x) and cos(x)
+// =================================================
+
+//
+// Uses the following identities:
+// ----------------------------------------------------------
+//
+// sin(-x) = -sin(x)
+// cos(-x) = cos(x)
+//
+// sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y)
+// cos(x+y) = sin(x)*sin(y)+cos(x)*cos(y)
+//
+// sin(x+ pi/2) = cos(x)
+// sin(x+ pi) = -sin(x)
+// sin(x+3pi/2) = -cos(x)
+// sin(x+2pi) = sin(x)
+//
+
+int sf_fsincos(floatx80 a, floatx80 *sin_a, floatx80 *cos_a)
+{
+ uint64_t aSig0, aSig1 = 0;
+ int32_t aExp, zExp, expDiff;
+ int aSign, zSign;
+ int q = 0;
+
+ aSig0 = extractFloatx80Frac(a);
+ aExp = extractFloatx80Exp(a);
+ aSign = extractFloatx80Sign(a);
+
+ /* invalid argument */
+ if (aExp == 0x7FFF) {
+ if ((uint64_t) (aSig0<<1)) {
+ sincos_invalid(sin_a, cos_a, propagateFloatx80NaNOneArg(a));
+ return 0;
+ }
+
+ float_raise(float_flag_invalid);
+ sincos_invalid(sin_a, cos_a, floatx80_default_nan);
+ return 0;
+ }
+
+ if (aExp == 0) {
+ if (aSig0 == 0) {
+ sincos_tiny_argument(sin_a, cos_a, a);
+ return 0;
+ }
+
+// float_raise(float_flag_denormal);
+
+ /* handle pseudo denormals */
+ if (! (aSig0 & 0x8000000000000000U))
+ {
+ float_raise(float_flag_inexact);
+ if (sin_a)
+ float_raise(float_flag_underflow);
+ sincos_tiny_argument(sin_a, cos_a, a);
+ return 0;
+ }
+
+ normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
+ }
+
+ zSign = aSign;
+ zExp = EXP_BIAS;
+ expDiff = aExp - zExp;
+
+ /* argument is out-of-range */
+ if (expDiff >= 63)
+ return -1;
+
+ float_raise(float_flag_inexact);
+
+ if (expDiff < -1) { // doesn't require reduction
+ if (expDiff <= -68) {
+ a = packFloatx80(aSign, aExp, aSig0);
+ sincos_tiny_argument(sin_a, cos_a, a);
+ return 0;
+ }
+ zExp = aExp;
+ }
+ else {
+ q = reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
+ }
+
+ /* **************************** */
+ /* argument reduction completed */
+ /* **************************** */
+
+ /* using float128 for approximation */
+ float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
+
+ if (aSign) q = -q;
+ if (sin_a) *sin_a = sincos_approximation(zSign, r, q);
+ if (cos_a) *cos_a = sincos_approximation(zSign, r, q+1);
+
+ return 0;
+}
+
+int floatx80_fsin(floatx80 *a)
+{
+ return sf_fsincos(*a, a, 0);
+}
+
+int floatx80_fcos(floatx80 *a)
+{
+ return sf_fsincos(*a, 0, a);
+}
+
+// =================================================
+// FPTAN Compute tan(x)
+// =================================================
+
+//
+// Uses the following identities:
+//
+// 1. ----------------------------------------------------------
+//
+// sin(-x) = -sin(x)
+// cos(-x) = cos(x)
+//
+// sin(x+y) = sin(x)*cos(y)+cos(x)*sin(y)
+// cos(x+y) = sin(x)*sin(y)+cos(x)*cos(y)
+//
+// sin(x+ pi/2) = cos(x)
+// sin(x+ pi) = -sin(x)
+// sin(x+3pi/2) = -cos(x)
+// sin(x+2pi) = sin(x)
+//
+// 2. ----------------------------------------------------------
+//
+// sin(x)
+// tan(x) = ------
+// cos(x)
+//
+
+int floatx80_ftan(floatx80 *a)
+{
+ uint64_t aSig0, aSig1 = 0;
+ int32_t aExp, zExp, expDiff;
+ int aSign, zSign;
+ int q = 0;
+
+ aSig0 = extractFloatx80Frac(*a);
+ aExp = extractFloatx80Exp(*a);
+ aSign = extractFloatx80Sign(*a);
+
+ /* invalid argument */
+ if (aExp == 0x7FFF) {
+ if ((uint64_t) (aSig0<<1))
+ {
+ *a = propagateFloatx80NaNOneArg(*a);
+ return 0;
+ }
+
+ float_raise(float_flag_invalid);
+ *a = floatx80_default_nan;
+ return 0;
+ }
+
+ if (aExp == 0) {
+ if (aSig0 == 0) return 0;
+// float_raise(float_flag_denormal);
+ /* handle pseudo denormals */
+ if (! (aSig0 & 0x8000000000000000U))
+ {
+ float_raise(float_flag_inexact | float_flag_underflow);
+ return 0;
+ }
+ normalizeFloatx80Subnormal(aSig0, &aExp, &aSig0);
+ }
+
+ zSign = aSign;
+ zExp = EXP_BIAS;
+ expDiff = aExp - zExp;
+
+ /* argument is out-of-range */
+ if (expDiff >= 63)
+ return -1;
+
+ float_raise(float_flag_inexact);
+
+ if (expDiff < -1) { // doesn't require reduction
+ if (expDiff <= -68) {
+ *a = packFloatx80(aSign, aExp, aSig0);
+ return 0;
+ }
+ zExp = aExp;
+ }
+ else {
+ q = reduce_trig_arg(expDiff, zSign, aSig0, aSig1);
+ }
+
+ /* **************************** */
+ /* argument reduction completed */
+ /* **************************** */
+
+ /* using float128 for approximation */
+ float128 r = normalizeRoundAndPackFloat128(0, zExp-0x10, aSig0, aSig1);
+
+ float128 sin_r = poly_sin(r);
+ float128 cos_r = poly_cos(r);
+
+ if (q & 0x1) {
+ r = float128_div(cos_r, sin_r);
+ zSign = ! zSign;
+ } else {
+ r = float128_div(sin_r, cos_r);
+ }
+
+ *a = float128_to_floatx80(r);
+ if (zSign)
+ *a = floatx80_chs(*a);
+
+ return 0;
+}
+
+// 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) ]
+//
+
+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);
+}
+
+// 2 4 6 8 2n
+// f(x) ~ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+// 0 1 2 3 4 n
+//
+// -- 4k -- 4k+2
+// p(x) = > C * x q(x) = > C * x
+// -- 2k -- 2k+1
+//
+// 2
+// f(x) ~ [ p(x) + x * q(x) ]
+//
+
+float128 EvenPoly(float128 x, float128 *arr, unsigned n)
+{
+ return EvalPoly(float128_mul(x, x), arr, n);
+}
+
+// 3 5 7 9 2n+1
+// f(x) ~ (C * x) + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+// 0 1 2 3 4 n
+// 2 4 6 8 2n
+// = x * [ C + (C * x) + (C * x) + (C * x) + (C * x) + ... + (C * x)
+// 0 1 2 3 4 n
+//
+// -- 4k -- 4k+2
+// p(x) = > C * x q(x) = > C * x
+// -- 2k -- 2k+1
+//
+// 2
+// f(x) ~ x * [ p(x) + x * q(x) ]
+//
+
+float128 OddPoly(float128 x, float128 *arr, unsigned n)
+{
+ return float128_mul(x, EvenPoly(x, arr, n));
+}
+
+/*----------------------------------------------------------------------------
+| Scales extended double-precision floating-point value in operand `a' by
+| value `b'. The function truncates the value in the second operand 'b' to
+| an integral value and adds that value to the exponent of the operand 'a'.
+| The operation performed according to the IEC/IEEE Standard for Binary
+| Floating-Point Arithmetic.
+*----------------------------------------------------------------------------*/
+
+extern floatx80 propagateFloatx80NaN( floatx80 a, floatx80 b );
+
+floatx80 floatx80_scale(floatx80 a, floatx80 b)
+{
+ sbits32 aExp, bExp;
+ bits64 aSig, bSig;
+
+ // handle unsupported extended double-precision floating encodings
+/* if (floatx80_is_unsupported(a) || floatx80_is_unsupported(b))
+ {
+ float_raise(float_flag_invalid);
+ return floatx80_default_nan;
+ }*/
+
+ aSig = extractFloatx80Frac(a);
+ aExp = extractFloatx80Exp(a);
+ int aSign = extractFloatx80Sign(a);
+ bSig = extractFloatx80Frac(b);
+ bExp = extractFloatx80Exp(b);
+ int bSign = extractFloatx80Sign(b);
+
+ if (aExp == 0x7FFF) {
+ if ((bits64) (aSig<<1) || ((bExp == 0x7FFF) && (bits64) (bSig<<1)))
+ {
+ return propagateFloatx80NaN(a, b);
+ }
+ if ((bExp == 0x7FFF) && bSign) {
+ float_raise(float_flag_invalid);
+ return floatx80_default_nan;
+ }
+ if (bSig && (bExp == 0)) float_raise(float_flag_denormal);
+ return a;
+ }
+ if (bExp == 0x7FFF) {
+ if ((bits64) (bSig<<1)) return propagateFloatx80NaN(a, b);
+ if ((aExp | aSig) == 0) {
+ if (! bSign) {
+ float_raise(float_flag_invalid);
+ return floatx80_default_nan;
+ }
+ return a;
+ }
+ if (aSig && (aExp == 0)) float_raise(float_flag_denormal);
+ if (bSign) return packFloatx80(aSign, 0, 0);
+ return packFloatx80(aSign, 0x7FFF, 0x8000000000000000U);
+ }
+ if (aExp == 0) {
+ if (aSig == 0) return a;
+ float_raise(float_flag_denormal);
+ normalizeFloatx80Subnormal(aSig, &aExp, &aSig);
+ }
+ if (bExp == 0) {
+ if (bSig == 0) return a;
+ float_raise(float_flag_denormal);
+ normalizeFloatx80Subnormal(bSig, &bExp, &bSig);
+ }
+
+ if (bExp > 0x400E) {
+ /* generate appropriate overflow/underflow */
+ return roundAndPackFloatx80(80, aSign,
+ bSign ? -0x3FFF : 0x7FFF, aSig, 0);
+ }
+ if (bExp < 0x3FFF) return a;
+
+ int shiftCount = 0x403E - bExp;
+ bSig >>= shiftCount;
+ sbits32 scale = bSig;
+ if (bSign) scale = -scale; /* -32768..32767 */
+ return
+ roundAndPackFloatx80(80, aSign, aExp+scale, aSig, 0);
+}
projects / pistorm / blobdiff