ext: Include libfputils

This changeset includes libfputils from revision bbf0d61d75. This
library can be used to convert to and from 80-bit floats and query the
type of an 80-bit float, which is needed to support the x87 FPU.

1 files changed, 247 insertions, 0 deletions

diff --git a/ext/fputils/fp80.c b/ext/fputils/fp80.c
new file mode 100644
index 000000000..703b07ac9
--- /dev/null
+++ b/ext/fputils/fp80.c
@@ -0,0 +1,247 @@
+/*
+ * Copyright (c) 2013, Andreas Sandberg
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ *    notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above
+ *    copyright notice, this list of conditions and the following
+ *    disclaimer in the documentation and/or other materials provided
+ *    with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+ * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+ * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+ * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
+ * OF THE POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#include <fputils/fp80.h>
+#include "fpbits.h"
+
+#include <assert.h>
+#include <stdint.h>
+
+#include <stdio.h>
+
+typedef union {
+    union {
+        uint64_t bits;
+        double value;
+    };
+} fp64_t;
+
+
+const fp80_t fp80_pinf = BUILD_FP80(0, 0, FP80_EXP_SPECIAL);
+const fp80_t fp80_ninf = BUILD_FP80(1, 0, FP80_EXP_SPECIAL);
+const fp80_t fp80_qnan = BUILD_FP80(0,  FP80_FRAC_QNAN, FP80_EXP_SPECIAL);
+const fp80_t fp80_qnani = BUILD_FP80(1, FP80_FRAC_QNANI, FP80_EXP_SPECIAL);
+const fp80_t fp80_snan = BUILD_FP80(0, FP80_FRAC_SNAN, FP80_EXP_SPECIAL);
+const fp80_t fp80_nan = BUILD_FP80(0, FP80_FRAC_QNAN, FP80_EXP_SPECIAL);
+
+static const fp64_t fp64_pinf = BUILD_FP64(0, 0, FP64_EXP_SPECIAL);
+static const fp64_t fp64_ninf = BUILD_FP64(1, 0, FP64_EXP_SPECIAL);
+static const fp64_t fp64_qnan = BUILD_FP64(0,  FP64_FRAC_QNAN,
+                                           FP64_EXP_SPECIAL);
+static const fp64_t fp64_nqnan = BUILD_FP64(1,  FP64_FRAC_QNAN,
+                                            FP64_EXP_SPECIAL);
+static const fp64_t fp64_qnani = BUILD_FP64(1, FP64_FRAC_QNANI,
+                                            FP64_EXP_SPECIAL);
+static const fp64_t fp64_snan = BUILD_FP64(0, FP64_FRAC_SNAN,
+                                           FP64_EXP_SPECIAL);
+static const fp64_t fp64_nsnan = BUILD_FP64(1, FP64_FRAC_SNAN,
+                                            FP64_EXP_SPECIAL);
+
+static double
+build_fp64(int sign, uint64_t frac, int exp)
+{
+    const fp64_t f = BUILD_FP64(sign, frac, exp);
+
+    return f.value;
+}
+
+int
+fp80_sgn(fp80_t fp80)
+{
+    return (fp80.u.repr.se & FP80_SIGN_BIT) ? -1 : 1;
+}
+
+int
+fp80_isspecial(fp80_t fp80)
+{
+    const int exp = FP80_EXP(fp80);
+
+    return exp == FP80_EXP_SPECIAL;
+}
+
+int
+fp80_isinf(fp80_t fp80)
+{
+    const uint64_t frac = FP80_FRAC(fp80);
+
+    return fp80_isspecial(fp80) && frac == 0 ? fp80_sgn(fp80) : 0;
+}
+
+
+int
+fp80_isqnan(fp80_t fp80)
+{
+    const uint64_t frac = FP80_FRAC(fp80);
+
+    return fp80_isspecial(fp80) && (frac & FP80_QNAN_BIT);
+}
+
+int
+fp80_isqnani(fp80_t fp80)
+{
+    const uint64_t frac_low = fp80.u.repr.fi & (FP80_FRAC_MASK >> 1);
+
+    return fp80_isqnan(fp80) && (fp80.u.repr.se & FP80_SIGN_BIT) && !frac_low;
+}
+
+int
+fp80_issnan(fp80_t fp80)
+{
+    const uint64_t frac = FP80_FRAC(fp80);
+
+    return fp80_isspecial(fp80) && !(frac & FP80_QNAN_BIT) && frac;
+}
+
+int
+fp80_isfinite(fp80_t fp80)
+{
+    return !fp80_isnan(fp80) && !fp80_isinf(fp80);
+}
+
+int
+fp80_isnan(fp80_t fp80)
+{
+    return fp80_issnan(fp80) || fp80_isqnan(fp80) ? fp80_sgn(fp80) : 0;
+}
+
+int
+fp80_iszero(fp80_t fp80)
+{
+    return fp80.u.repr.fi == 0 && FP80_EXP(fp80) == 0 ? fp80_sgn(fp80) : 0;
+}
+
+int
+fp80_isnormal(fp80_t fp80)
+{
+    return FP80_EXP(fp80) != 0 && !fp80_isspecial(fp80) ?
+        fp80_sgn(fp80) : 0;
+}
+
+int
+fp80_issubnormal(fp80_t fp80)
+{
+    return FP80_FRAC(fp80) && FP80_EXP(fp80) == 0 ? fp80_sgn(fp80) : 0;
+}
+
+int
+fp80_classify(fp80_t fp80)
+{
+    if (fp80_issubnormal(fp80)) {
+        return FP_SUBNORMAL;
+    } else if (fp80_iszero(fp80)) {
+        return FP_ZERO;
+    } else if (fp80_isinf(fp80)) {
+        return FP_INFINITE;
+    } else if (fp80_isnan(fp80)) {
+        return FP_NAN;
+    } else {
+        assert(fp80_isfinite(fp80));
+        return FP_NORMAL;
+    }
+}
+
+double
+fp80_cvtd(fp80_t fp80)
+{
+    const int sign = fp80.u.repr.se & FP80_SIGN_BIT;
+
+    if (!fp80_isspecial(fp80)) {
+        const uint64_t frac = fp80.u.repr.fi;
+        const int unb_exp = FP80_EXP(fp80) - FP80_EXP_BIAS;
+        const int fp64_exp = unb_exp + FP64_EXP_BIAS;
+        const uint64_t fp64_frac = frac >> (FP80_FRAC_BITS - FP64_FRAC_BITS);
+
+        if (fp64_exp > 0 && fp64_exp < FP64_EXP_SPECIAL) {
+            /* These numbers fall in the range of what we can express
+             * as normals */
+            return build_fp64(sign, fp64_frac, fp64_exp);
+        } else if (fp64_exp <= 0) {
+            uint64_t fp64_denormal_frac = fp64_frac >> (-fp64_exp);
+            /* Generate a denormal or zero */
+            return build_fp64(sign, fp64_denormal_frac, 0);
+        } else {
+            /* Infinity */
+            return build_fp64(sign, 0, FP64_EXP_SPECIAL);
+        }
+    } else {
+        if (fp80_isinf(fp80)) {
+            return build_fp64(sign, 0, FP64_EXP_SPECIAL);
+        } else if (fp80_issnan(fp80)) {
+            return fp80_sgn(fp80) > 0 ? fp64_snan.value : fp64_nsnan.value;
+        } else if (fp80_isqnani(fp80)) {
+            return fp64_qnani.value;
+        } else {
+            assert(fp80_isqnan(fp80));
+            return fp80_sgn(fp80) > 0 ? fp64_qnan.value : fp64_nqnan.value;
+        }
+    }
+}
+
+fp80_t
+fp80_cvfd(double value)
+{
+    const fp64_t fp64 = { .value = value };
+    const uint64_t frac = FP64_FRAC(fp64);
+    const unsigned exp = FP64_EXP(fp64);
+    const int unb_exp = exp - FP64_EXP_BIAS;
+    const uint64_t fp80_frac = frac << (FP80_FRAC_BITS - FP64_FRAC_BITS);
+
+    if (exp != 0) {
+        // Normal, inf, nan
+        const unsigned fp80_exp = exp == FP64_EXP_SPECIAL ?
+            FP80_EXP_SPECIAL : (unb_exp + FP80_EXP_BIAS);
+        const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT,
+                                       fp80_frac, fp80_exp);
+        return fp80;
+    } else if (exp == 0 && frac == 0) {
+        // Zero
+        const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT, 0, 0);
+        return fp80;
+    } else {
+        // Denormal
+        uint64_t fp80_fi = fp80_frac;
+        int shift_amt = 0;
+        while (!(fp80_fi & FP80_INT_BIT)) {
+            fp80_fi <<= 1;
+            ++shift_amt;
+        }
+        const unsigned fp80_exp = (unb_exp - shift_amt) + FP80_EXP_BIAS;
+        const fp80_t fp80 = BUILD_FP80(fp64.bits & FP64_SIGN_BIT,
+                                       fp80_fi, fp80_exp);
+        return fp80;
+    }
+}
+
+void
+fp80_debug_dump(FILE *fout, fp80_t fp80)
+{
+    fprintf(fout, "sgn: %i, int: %i, frac: 0x%llx, exp: 0x%x (%i)\n",
+            fp80_sgn(fp80), !!(fp80.u.repr.fi & FP80_INT_BIT), FP80_FRAC(fp80),
+            FP80_EXP(fp80), FP80_EXP(fp80) - FP80_EXP_BIAS);
+}