[RFC PATCH 0/8] fpu: experimental conversion of float128_addsub

[RFC PATCH 0/8] fpu: experimental conversion of float128_addsub

[Alex Bennée]

Hi Richard,

This is the current state of my experiment to convert a 128 bit float
function (in this case addsub). The actual conversion was fairly
simple (basically a copy & paste with some tweaks for using the
unint128 inline functions). However I ran into a number of stumbles
with the int128.h support including casting of values like ~(Uint128)0
and messing around to handle things like missing __builtin support for
clz. I suspect having some of the #defines expand into uint128_*
functions plays some part in the 4x growth in code compared to the old
version. However the drop in performance is a lot less than that.

In terms of total code churn we replace each deleted line in
softfloat.c with 2 lines of new code although I suspect if we pressed
on we could reduce the diffstat deficit. Debugging the actual failures
was relatively painless with rr and the new code - perhaps because I
just find it easier to follow.

I've included your early patches as that happened to be the state of
my tree when I branched off. If we want to go forward with a more
complete conversion I guess we would need:

  - a more complete int128.h conversion (including fallback for non CONFIG_INT128)
  - seeing if some of the resulting code bloat can be reduced
  - seeing what scope there is for commonality of special case handling

I'm not a fan of having so much duplication but at least I personally
find the code is more readable.

Alex Bennée (3):
  int128.h: add bunch of uint128 utility functions (INCOMPLETE)
  tests/fp: add quad support to the benchmark utility
  softfloat: implement addsub_floats128 using Uint128 and new style code

Richard Henderson (5):
  softfloat: Use mulu64 for mul64To128
  softfloat: Use int128.h for some operations
  softfloat: Tidy a * b + inf return
  softfloat: Add float_cmask and constants
  softfloat: Inline pick_nan_muladd into its caller

 include/fpu/softfloat-macros.h |  80 ++--
 include/qemu/int128.h          | 122 ++++++
 fpu/softfloat.c                | 697 ++++++++++++++++++++-------------
 tests/fp/fp-bench.c            |  88 ++++-
 fpu/softfloat-specialize.c.inc |  39 ++
 5 files changed, 711 insertions(+), 315 deletions(-)

-- 
2.20.1

[RFC PATCH 1/8] softfloat: Use mulu64 for mul64To128

[Alex Bennée]

From: Richard Henderson <richard.henderson@linaro.org>

Via host-utils.h, we use a host widening multiply for
64-bit hosts, and a common subroutine for 32-bit hosts.

Reviewed-by: David Hildenbrand <david@redhat.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Message-Id: <20200925152047.709901-2-richard.henderson@linaro.org>
---
 include/fpu/softfloat-macros.h | 24 ++++--------------------
 1 file changed, 4 insertions(+), 20 deletions(-)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index a35ec2893a..57845f8af0 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -83,6 +83,7 @@ this code that are retained.
 #define FPU_SOFTFLOAT_MACROS_H
 
 #include "fpu/softfloat-types.h"
+#include "qemu/host-utils.h"
 
 /*----------------------------------------------------------------------------
 | Shifts `a' right by the number of bits given in `count'.  If any nonzero
@@ -515,27 +516,10 @@ static inline void
 | `z0Ptr' and `z1Ptr'.
 *----------------------------------------------------------------------------*/
 
-static inline void mul64To128( uint64_t a, uint64_t b, uint64_t *z0Ptr, uint64_t *z1Ptr )
+static inline void
+mul64To128(uint64_t a, uint64_t b, uint64_t *z0Ptr, uint64_t *z1Ptr)
 {
-    uint32_t aHigh, aLow, bHigh, bLow;
-    uint64_t z0, zMiddleA, zMiddleB, z1;
-
-    aLow = a;
-    aHigh = a>>32;
-    bLow = b;
-    bHigh = b>>32;
-    z1 = ( (uint64_t) aLow ) * bLow;
-    zMiddleA = ( (uint64_t) aLow ) * bHigh;
-    zMiddleB = ( (uint64_t) aHigh ) * bLow;
-    z0 = ( (uint64_t) aHigh ) * bHigh;
-    zMiddleA += zMiddleB;
-    z0 += ( ( (uint64_t) ( zMiddleA < zMiddleB ) )<<32 ) + ( zMiddleA>>32 );
-    zMiddleA <<= 32;
-    z1 += zMiddleA;
-    z0 += ( z1 < zMiddleA );
-    *z1Ptr = z1;
-    *z0Ptr = z0;
-
+    mulu64(z1Ptr, z0Ptr, a, b);
 }
 
 /*----------------------------------------------------------------------------
-- 
2.20.1

[RFC PATCH 2/8] softfloat: Use int128.h for some operations

[Alex Bennée]

From: Richard Henderson <richard.henderson@linaro.org>

Use our Int128, which wraps the compiler's __int128_t,
instead of open-coding left shifts and arithmetic.
We'd need to extend Int128 to have unsigned operations
to replace more than these three.

Reviewed-by: David Hildenbrand <david@redhat.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Message-Id: <20200925152047.709901-3-richard.henderson@linaro.org>
---
 include/fpu/softfloat-macros.h | 39 +++++++++++++++++-----------------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index 57845f8af0..95d88d05b8 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -84,6 +84,7 @@ this code that are retained.
 
 #include "fpu/softfloat-types.h"
 #include "qemu/host-utils.h"
+#include "qemu/int128.h"
 
 /*----------------------------------------------------------------------------
 | Shifts `a' right by the number of bits given in `count'.  If any nonzero
@@ -352,13 +353,11 @@ static inline void shortShift128Left(uint64_t a0, uint64_t a1, int count,
 static inline void shift128Left(uint64_t a0, uint64_t a1, int count,
                                 uint64_t *z0Ptr, uint64_t *z1Ptr)
 {
-    if (count < 64) {
-        *z1Ptr = a1 << count;
-        *z0Ptr = count == 0 ? a0 : (a0 << count) | (a1 >> (-count & 63));
-    } else {
-        *z1Ptr = 0;
-        *z0Ptr = a1 << (count - 64);
-    }
+    Int128 a = int128_make128(a1, a0);
+    Int128 z = int128_lshift(a, count);
+
+    *z0Ptr = int128_gethi(z);
+    *z1Ptr = int128_getlo(z);
 }
 
 /*----------------------------------------------------------------------------
@@ -405,15 +404,15 @@ static inline void
 *----------------------------------------------------------------------------*/
 
 static inline void
- add128(
-     uint64_t a0, uint64_t a1, uint64_t b0, uint64_t b1, uint64_t *z0Ptr, uint64_t *z1Ptr )
+add128(uint64_t a0, uint64_t a1, uint64_t b0, uint64_t b1,
+       uint64_t *z0Ptr, uint64_t *z1Ptr)
 {
-    uint64_t z1;
-
-    z1 = a1 + b1;
-    *z1Ptr = z1;
-    *z0Ptr = a0 + b0 + ( z1 < a1 );
+    Int128 a = int128_make128(a1, a0);
+    Int128 b = int128_make128(b1, b0);
+    Int128 z = int128_add(a, b);
 
+    *z0Ptr = int128_gethi(z);
+    *z1Ptr = int128_getlo(z);
 }
 
 /*----------------------------------------------------------------------------
@@ -463,13 +462,15 @@ static inline void
 *----------------------------------------------------------------------------*/
 
 static inline void
- sub128(
-     uint64_t a0, uint64_t a1, uint64_t b0, uint64_t b1, uint64_t *z0Ptr, uint64_t *z1Ptr )
+sub128(uint64_t a0, uint64_t a1, uint64_t b0, uint64_t b1,
+       uint64_t *z0Ptr, uint64_t *z1Ptr)
 {
+    Int128 a = int128_make128(a1, a0);
+    Int128 b = int128_make128(b1, b0);
+    Int128 z = int128_sub(a, b);
 
-    *z1Ptr = a1 - b1;
-    *z0Ptr = a0 - b0 - ( a1 < b1 );
-
+    *z0Ptr = int128_gethi(z);
+    *z1Ptr = int128_getlo(z);
 }
 
 /*----------------------------------------------------------------------------
-- 
2.20.1

[RFC PATCH 3/8] softfloat: Tidy a * b + inf return

[Alex Bennée]

From: Richard Henderson <richard.henderson@linaro.org>

No reason to set values in 'a', when we already
have float_class_inf in 'c', and can flip that sign.

Reviewed-by: David Hildenbrand <david@redhat.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Message-Id: <20200925152047.709901-4-richard.henderson@linaro.org>
---
 fpu/softfloat.c | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index 67cfa0fd82..9db55d2b11 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -1380,9 +1380,8 @@ static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c,
             s->float_exception_flags |= float_flag_invalid;
             return parts_default_nan(s);
         } else {
-            a.cls = float_class_inf;
-            a.sign = c.sign ^ sign_flip;
-            return a;
+            c.sign ^= sign_flip;
+            return c;
         }
     }
 
-- 
2.20.1

[RFC PATCH 4/8] softfloat: Add float_cmask and constants

[Alex Bennée]

From: Richard Henderson <richard.henderson@linaro.org>

Testing more than one class at a time is better done with masks.
This reduces the static branch count.

Reviewed-by: David Hildenbrand <david@redhat.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Message-Id: <20200925152047.709901-5-richard.henderson@linaro.org>
---
 fpu/softfloat.c | 31 ++++++++++++++++++++++++-------
 1 file changed, 24 insertions(+), 7 deletions(-)

diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index 9db55d2b11..3e625c47cd 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -469,6 +469,20 @@ typedef enum __attribute__ ((__packed__)) {
     float_class_snan,
 } FloatClass;
 
+#define float_cmask(bit)  (1u << (bit))
+
+enum {
+    float_cmask_zero    = float_cmask(float_class_zero),
+    float_cmask_normal  = float_cmask(float_class_normal),
+    float_cmask_inf     = float_cmask(float_class_inf),
+    float_cmask_qnan    = float_cmask(float_class_qnan),
+    float_cmask_snan    = float_cmask(float_class_snan),
+
+    float_cmask_infzero = float_cmask_zero | float_cmask_inf,
+    float_cmask_anynan  = float_cmask_qnan | float_cmask_snan,
+};
+
+
 /* Simple helpers for checking if, or what kind of, NaN we have */
 static inline __attribute__((unused)) bool is_nan(FloatClass c)
 {
@@ -1335,24 +1349,27 @@ bfloat16 QEMU_FLATTEN bfloat16_mul(bfloat16 a, bfloat16 b, float_status *status)
 static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c,
                                 int flags, float_status *s)
 {
-    bool inf_zero = ((1 << a.cls) | (1 << b.cls)) ==
-                    ((1 << float_class_inf) | (1 << float_class_zero));
-    bool p_sign;
+    bool inf_zero, p_sign;
     bool sign_flip = flags & float_muladd_negate_result;
     FloatClass p_class;
     uint64_t hi, lo;
     int p_exp;
+    int ab_mask, abc_mask;
+
+    ab_mask = float_cmask(a.cls) | float_cmask(b.cls);
+    abc_mask = float_cmask(c.cls) | ab_mask;
+    inf_zero = ab_mask == float_cmask_infzero;
 
     /* It is implementation-defined whether the cases of (0,inf,qnan)
      * and (inf,0,qnan) raise InvalidOperation or not (and what QNaN
      * they return if they do), so we have to hand this information
      * off to the target-specific pick-a-NaN routine.
      */
-    if (is_nan(a.cls) || is_nan(b.cls) || is_nan(c.cls)) {
+    if (unlikely(abc_mask & float_cmask_anynan)) {
         return pick_nan_muladd(a, b, c, inf_zero, s);
     }
 
-    if (inf_zero) {
+    if (unlikely(inf_zero)) {
         s->float_exception_flags |= float_flag_invalid;
         return parts_default_nan(s);
     }
@@ -1367,9 +1384,9 @@ static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c,
         p_sign ^= 1;
     }
 
-    if (a.cls == float_class_inf || b.cls == float_class_inf) {
+    if (ab_mask & float_cmask_inf) {
         p_class = float_class_inf;
-    } else if (a.cls == float_class_zero || b.cls == float_class_zero) {
+    } else if (ab_mask & float_cmask_zero) {
         p_class = float_class_zero;
     } else {
         p_class = float_class_normal;
-- 
2.20.1

[RFC PATCH 5/8] softfloat: Inline pick_nan_muladd into its caller

[Alex Bennée]

From: Richard Henderson <richard.henderson@linaro.org>

Because of FloatParts, there will only ever be one caller.
Inlining allows us to re-use abc_mask for the snan test.

Reviewed-by: David Hildenbrand <david@redhat.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
Message-Id: <20200925152047.709901-6-richard.henderson@linaro.org>
---
 fpu/softfloat.c | 75 +++++++++++++++++++++++--------------------------
 1 file changed, 35 insertions(+), 40 deletions(-)

diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index 3e625c47cd..e038434a07 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -929,45 +929,6 @@ static FloatParts pick_nan(FloatParts a, FloatParts b, float_status *s)
     return a;
 }
 
-static FloatParts pick_nan_muladd(FloatParts a, FloatParts b, FloatParts c,
-                                  bool inf_zero, float_status *s)
-{
-    int which;
-
-    if (is_snan(a.cls) || is_snan(b.cls) || is_snan(c.cls)) {
-        s->float_exception_flags |= float_flag_invalid;
-    }
-
-    which = pickNaNMulAdd(a.cls, b.cls, c.cls, inf_zero, s);
-
-    if (s->default_nan_mode) {
-        /* Note that this check is after pickNaNMulAdd so that function
-         * has an opportunity to set the Invalid flag.
-         */
-        which = 3;
-    }
-
-    switch (which) {
-    case 0:
-        break;
-    case 1:
-        a = b;
-        break;
-    case 2:
-        a = c;
-        break;
-    case 3:
-        return parts_default_nan(s);
-    default:
-        g_assert_not_reached();
-    }
-
-    if (is_snan(a.cls)) {
-        return parts_silence_nan(a, s);
-    }
-    return a;
-}
-
 /*
  * Returns the result of adding or subtracting the values of the
  * floating-point values `a' and `b'. The operation is performed
@@ -1366,7 +1327,41 @@ static FloatParts muladd_floats(FloatParts a, FloatParts b, FloatParts c,
      * off to the target-specific pick-a-NaN routine.
      */
     if (unlikely(abc_mask & float_cmask_anynan)) {
-        return pick_nan_muladd(a, b, c, inf_zero, s);
+        int which;
+
+        if (unlikely(abc_mask & float_cmask_snan)) {
+            float_raise(float_flag_invalid, s);
+        }
+
+        which = pickNaNMulAdd(a.cls, b.cls, c.cls, inf_zero, s);
+
+        if (s->default_nan_mode) {
+            /*
+             * Note that this check is after pickNaNMulAdd so that function
+             * has an opportunity to set the Invalid flag for inf_zero.
+             */
+            which = 3;
+        }
+
+        switch (which) {
+        case 0:
+            break;
+        case 1:
+            a = b;
+            break;
+        case 2:
+            a = c;
+            break;
+        case 3:
+            return parts_default_nan(s);
+        default:
+            g_assert_not_reached();
+        }
+
+        if (is_snan(a.cls)) {
+            return parts_silence_nan(a, s);
+        }
+        return a;
     }
 
     if (unlikely(inf_zero)) {
-- 
2.20.1

[RFC PATCH 6/8] int128.h: add bunch of uint128 utility functions (INCOMPLETE)

[Alex Bennée]

These will be useful for softfloat. I've included the extract/desposit
functions with the main Int128 header and not with cutils as we need
alternate versions for systems that don't have compiler support for
Uint128. Even with compiler support some stuff we need to
hand-hack (like clz128).

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
---
 include/qemu/int128.h | 122 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 122 insertions(+)

diff --git a/include/qemu/int128.h b/include/qemu/int128.h
index 76ea405922..38c8b1ab29 100644
--- a/include/qemu/int128.h
+++ b/include/qemu/int128.h
@@ -3,8 +3,10 @@
 
 #ifdef CONFIG_INT128
 #include "qemu/bswap.h"
+#include "qemu/host-utils.h"
 
 typedef __int128_t Int128;
+typedef __uint128_t Uint128;
 
 static inline Int128 int128_make64(uint64_t a)
 {
@@ -16,6 +18,11 @@ static inline Int128 int128_make128(uint64_t lo, uint64_t hi)
     return (__uint128_t)hi << 64 | lo;
 }
 
+static inline Uint128 uint128_make128(uint64_t lo, uint64_t hi)
+{
+    return (__uint128_t)hi << 64 | lo;
+}
+
 static inline uint64_t int128_get64(Int128 a)
 {
     uint64_t r = a;
@@ -28,16 +35,31 @@ static inline uint64_t int128_getlo(Int128 a)
     return a;
 }
 
+static inline uint64_t uint128_getlo(Uint128 a)
+{
+    return a;
+}
+
 static inline int64_t int128_gethi(Int128 a)
 {
     return a >> 64;
 }
 
+static inline uint64_t uint128_gethi(Uint128 a)
+{
+    return a >> 64;
+}
+
 static inline Int128 int128_zero(void)
 {
     return 0;
 }
 
+static inline Uint128 uint128_zero(void)
+{
+    return 0;
+}
+
 static inline Int128 int128_one(void)
 {
     return 1;
@@ -58,21 +80,51 @@ static inline Int128 int128_and(Int128 a, Int128 b)
     return a & b;
 }
 
+static inline Uint128 uint128_and(Uint128 a, Uint128 b)
+{
+    return a & b;
+}
+
+static inline Int128 int128_or(Int128 a, Int128 b)
+{
+    return a | b;
+}
+
+static inline Uint128 uint128_or(Uint128 a, Uint128 b)
+{
+    return a | b;
+}
+
 static inline Int128 int128_rshift(Int128 a, int n)
 {
     return a >> n;
 }
 
+static inline Uint128 uint128_rshift(Uint128 a, int n)
+{
+    return a >> n;
+}
+
 static inline Int128 int128_lshift(Int128 a, int n)
 {
     return a << n;
 }
 
+static inline Uint128 uint128_lshift(Uint128 a, int n)
+{
+    return a << n;
+}
+
 static inline Int128 int128_add(Int128 a, Int128 b)
 {
     return a + b;
 }
 
+static inline Uint128 uint128_add(Uint128 a, Uint128 b)
+{
+    return a + b;
+}
+
 static inline Int128 int128_neg(Int128 a)
 {
     return -a;
@@ -83,6 +135,11 @@ static inline Int128 int128_sub(Int128 a, Int128 b)
     return a - b;
 }
 
+static inline Uint128 uint128_sub(Uint128 a, Uint128 b)
+{
+    return a - b;
+}
+
 static inline bool int128_nonneg(Int128 a)
 {
     return a >= 0;
@@ -93,6 +150,11 @@ static inline bool int128_eq(Int128 a, Int128 b)
     return a == b;
 }
 
+static inline bool uint128_eq(Uint128 a, Uint128 b)
+{
+    return a == b;
+}
+
 static inline bool int128_ne(Int128 a, Int128 b)
 {
     return a != b;
@@ -148,6 +210,66 @@ static inline Int128 bswap128(Int128 a)
     return int128_make128(bswap64(int128_gethi(a)), bswap64(int128_getlo(a)));
 }
 
+/**
+ * extract128:
+ * @value: the value to extract the bit field from
+ * @start: the lowest bit in the bit field (numbered from 0)
+ * @length: the length of the bit field
+ *
+ * Extract from the 128 bit input @value the bit field specified by the
+ * @start and @length parameters, and return it. The bit field must
+ * lie entirely within the 128 bit word. It is valid to request that
+ * all 128 bits are returned (ie @length 128 and @start 0).
+ *
+ * Returns: the value of the bit field extracted from the input value.
+ */
+static inline Uint128 extract128(Uint128 value, int start, int length)
+{
+    assert(start >= 0 && length > 0 && length <= 128 - start);
+    Uint128 mask = ~(Uint128)0 >> (128 - length);
+    Uint128 shifted = value >> start;
+    return shifted & mask;
+}
+
+/**
+ * deposit128:
+ * @value: initial value to insert bit field into
+ * @start: the lowest bit in the bit field (numbered from 0)
+ * @length: the length of the bit field
+ * @fieldval: the value to insert into the bit field
+ *
+ * Deposit @fieldval into the 128 bit @value at the bit field specified
+ * by the @start and @length parameters, and return the modified
+ * @value. Bits of @value outside the bit field are not modified.
+ * Bits of @fieldval above the least significant @length bits are
+ * ignored. The bit field must lie entirely within the 128 bit word.
+ * It is valid to request that all 128 bits are modified (ie @length
+ * 128 and @start 0).
+ *
+ * Returns: the modified @value.
+ */
+static inline Uint128 deposit128(Uint128 value, int start, int length,
+                                 Uint128 fieldval)
+{
+    assert(start >= 0 && length > 0 && length <= 128 - start);
+    Uint128 mask = (~(Uint128)0 >> (128 - length)) << start;
+    return (value & ~mask) | ((fieldval << start) & mask);
+}
+
+static inline int clz128(Uint128 val)
+{
+    if (val) {
+        uint64_t hi = uint128_gethi(val);
+        if (hi) {
+            return clz64(hi);
+        } else {
+            return 64 + clz64(uint128_getlo(val));
+        }
+    } else {
+        return 128;
+    }
+}
+
 #else /* !CONFIG_INT128 */
 
 typedef struct Int128 Int128;
-- 
2.20.1

[RFC PATCH 7/8] tests/fp: add quad support to the benchmark utility

[Alex Bennée]

Currently this only support softfloat calculations because working out
if the hardware supports 128 bit floats needs configure magic. The 3
op muladd operation is currently unimplemented so commented out for
now.

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
---
 tests/fp/fp-bench.c | 88 ++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 83 insertions(+), 5 deletions(-)

diff --git a/tests/fp/fp-bench.c b/tests/fp/fp-bench.c
index 4ba5e1d2d4..d319993280 100644
--- a/tests/fp/fp-bench.c
+++ b/tests/fp/fp-bench.c
@@ -14,6 +14,7 @@
 #include <math.h>
 #include <fenv.h>
 #include "qemu/timer.h"
+#include "qemu/int128.h"
 #include "fpu/softfloat.h"
 
 /* amortize the computation of random inputs */
@@ -50,8 +51,10 @@ static const char * const op_names[] = {
 enum precision {
     PREC_SINGLE,
     PREC_DOUBLE,
+    PREC_QUAD,
     PREC_FLOAT32,
     PREC_FLOAT64,
+    PREC_FLOAT128,
     PREC_MAX_NR,
 };
 
@@ -89,6 +92,7 @@ union fp {
     double d;
     float32 f32;
     float64 f64;
+    float128 f128;
     uint64_t u64;
 };
 
@@ -113,6 +117,10 @@ struct op_desc {
 static uint64_t random_ops[MAX_OPERANDS] = {
     SEED_A, SEED_B, SEED_C,
 };
+
+static float128 random_quad_ops[MAX_OPERANDS] = {
+    {SEED_A, SEED_B}, {SEED_B, SEED_C}, {SEED_C, SEED_A},
+};
 static float_status soft_status;
 static enum precision precision;
 static enum op operation;
@@ -141,25 +149,45 @@ static void update_random_ops(int n_ops, enum precision prec)
     int i;
 
     for (i = 0; i < n_ops; i++) {
-        uint64_t r = random_ops[i];
 
         switch (prec) {
         case PREC_SINGLE:
         case PREC_FLOAT32:
+        {
+            uint64_t r = random_ops[i];
             do {
                 r = xorshift64star(r);
             } while (!float32_is_normal(r));
+            random_ops[i] = r;
             break;
+        }
         case PREC_DOUBLE:
         case PREC_FLOAT64:
+        {
+            uint64_t r = random_ops[i];
             do {
                 r = xorshift64star(r);
             } while (!float64_is_normal(r));
+            random_ops[i] = r;
+            break;
+        }
+        case PREC_QUAD:
+        case PREC_FLOAT128:
+        {
+            float128 r = random_quad_ops[i];
+            uint64_t hi = r.high;
+            uint64_t lo = r.low;
+            do {
+                hi = xorshift64star(hi);
+                lo = xorshift64star(lo);
+                r = make_float128(hi, lo);
+            } while (!float128_is_normal(r));
+            random_quad_ops[i] = r;
             break;
+        }
         default:
             g_assert_not_reached();
         }
-        random_ops[i] = r;
     }
 }
 
@@ -184,6 +212,13 @@ static void fill_random(union fp *ops, int n_ops, enum precision prec,
                 ops[i].f64 = float64_chs(ops[i].f64);
             }
             break;
+        case PREC_QUAD:
+        case PREC_FLOAT128:
+            ops[i].f128 = random_quad_ops[i];
+            if (no_neg && float128_is_neg(ops[i].f128)) {
+                ops[i].f128 = float128_chs(ops[i].f128);
+            }
+            break;
         default:
             g_assert_not_reached();
         }
@@ -345,6 +380,41 @@ static void bench(enum precision prec, enum op op, int n_ops, bool no_neg)
                 }
             }
             break;
+        case PREC_FLOAT128:
+            fill_random(ops, n_ops, prec, no_neg);
+            t0 = get_clock();
+            for (i = 0; i < OPS_PER_ITER; i++) {
+                float128 a = ops[0].f128;
+                float128 b = ops[1].f128;
+                /* float128 c = ops[2].f128; */
+
+                switch (op) {
+                case OP_ADD:
+                    res.f128 = float128_add(a, b, &soft_status);
+                    break;
+                case OP_SUB:
+                    res.f128 = float128_sub(a, b, &soft_status);
+                    break;
+                case OP_MUL:
+                    res.f128 = float128_mul(a, b, &soft_status);
+                    break;
+                case OP_DIV:
+                    res.f128 = float128_div(a, b, &soft_status);
+                    break;
+                /* case OP_FMA: */
+                /*     res.f128 = float128_muladd(a, b, c, 0, &soft_status); */
+                /*     break; */
+                case OP_SQRT:
+                    res.f128 = float128_sqrt(a, &soft_status);
+                    break;
+                case OP_CMP:
+                    res.u64 = float128_compare_quiet(a, b, &soft_status);
+                    break;
+                default:
+                    g_assert_not_reached();
+                }
+            }
+            break;
         default:
             g_assert_not_reached();
         }
@@ -369,7 +439,8 @@ static void bench(enum precision prec, enum op op, int n_ops, bool no_neg)
     GEN_BENCH(bench_ ## opname ## _float, float, PREC_SINGLE, op, n_ops) \
     GEN_BENCH(bench_ ## opname ## _double, double, PREC_DOUBLE, op, n_ops) \
     GEN_BENCH(bench_ ## opname ## _float32, float32, PREC_FLOAT32, op, n_ops) \
-    GEN_BENCH(bench_ ## opname ## _float64, float64, PREC_FLOAT64, op, n_ops)
+    GEN_BENCH(bench_ ## opname ## _float64, float64, PREC_FLOAT64, op, n_ops) \
+    GEN_BENCH(bench_ ## opname ## _float128, float128, PREC_FLOAT128, op, n_ops)
 
 GEN_BENCH_ALL_TYPES(add, OP_ADD, 2)
 GEN_BENCH_ALL_TYPES(sub, OP_SUB, 2)
@@ -383,7 +454,8 @@ GEN_BENCH_ALL_TYPES(cmp, OP_CMP, 2)
     GEN_BENCH_NO_NEG(bench_ ## name ## _float, float, PREC_SINGLE, op, n) \
     GEN_BENCH_NO_NEG(bench_ ## name ## _double, double, PREC_DOUBLE, op, n) \
     GEN_BENCH_NO_NEG(bench_ ## name ## _float32, float32, PREC_FLOAT32, op, n) \
-    GEN_BENCH_NO_NEG(bench_ ## name ## _float64, float64, PREC_FLOAT64, op, n)
+    GEN_BENCH_NO_NEG(bench_ ## name ## _float64, float64, PREC_FLOAT64, op, n) \
+    GEN_BENCH_NO_NEG(bench_ ## name ## _float128, float128, PREC_FLOAT128, op, n)
 
 GEN_BENCH_ALL_TYPES_NO_NEG(sqrt, OP_SQRT, 1)
 #undef GEN_BENCH_ALL_TYPES_NO_NEG
@@ -397,6 +469,7 @@ GEN_BENCH_ALL_TYPES_NO_NEG(sqrt, OP_SQRT, 1)
         [PREC_DOUBLE]    = bench_ ## opname ## _double,         \
         [PREC_FLOAT32]   = bench_ ## opname ## _float32,        \
         [PREC_FLOAT64]   = bench_ ## opname ## _float64,        \
+        [PREC_FLOAT128]   = bench_ ## opname ## _float128,      \
     }
 
 static const bench_func_t bench_funcs[OP_MAX_NR][PREC_MAX_NR] = {
@@ -445,7 +518,7 @@ static void usage_complete(int argc, char *argv[])
     fprintf(stderr, " -h = show this help message.\n");
     fprintf(stderr, " -o = floating point operation (%s). Default: %s\n",
             op_list, op_names[0]);
-    fprintf(stderr, " -p = floating point precision (single, double). "
+    fprintf(stderr, " -p = floating point precision (single, double, quad[soft only]). "
             "Default: single\n");
     fprintf(stderr, " -r = rounding mode (even, zero, down, up, tieaway). "
             "Default: even\n");
@@ -565,6 +638,8 @@ static void parse_args(int argc, char *argv[])
                 precision = PREC_SINGLE;
             } else if (!strcmp(optarg, "double")) {
                 precision = PREC_DOUBLE;
+            } else if (!strcmp(optarg, "quad")) {
+                precision = PREC_QUAD;
             } else {
                 fprintf(stderr, "Unsupported precision '%s'\n", optarg);
                 exit(EXIT_FAILURE);
@@ -608,6 +683,9 @@ static void parse_args(int argc, char *argv[])
         case PREC_DOUBLE:
             precision = PREC_FLOAT64;
             break;
+        case PREC_QUAD:
+            precision = PREC_FLOAT128;
+            break;
         default:
             g_assert_not_reached();
         }
-- 
2.20.1

[RFC PATCH 8/8] softfloat: implement addsub_floats128 using Uint128 and new style code

[Alex Bennée]

This inevitably leads to duplication of almost identical functions for
the 128 bit and the 64/32/16 bit code. We will try and address that
later.

The expanded code is longer than the previous softfloat code and we
see a drop in performance from ~85 MFlops to ~72 MFlops.

Signed-off-by: Alex Bennée <alex.bennee@linaro.org>
---
 include/fpu/softfloat-macros.h |  17 +
 fpu/softfloat.c                | 586 ++++++++++++++++++++-------------
 fpu/softfloat-specialize.c.inc |  39 +++
 3 files changed, 421 insertions(+), 221 deletions(-)

diff --git a/include/fpu/softfloat-macros.h b/include/fpu/softfloat-macros.h
index 95d88d05b8..8e6158886d 100644
--- a/include/fpu/softfloat-macros.h
+++ b/include/fpu/softfloat-macros.h
@@ -259,6 +259,23 @@ static inline void
 
 }
 
+/* as above with Uint128 */
+static inline Uint128 shift128RightJammingUint128(Int128 in, int count)
+{
+    Uint128 out;
+
+    if (count == 0) {
+        out = in;
+    } else if ( count < 128 ) {
+        Uint128 jam_mask = uint128_sub(uint128_lshift(1, count), 1);
+        int shifted_away = uint128_and(in, jam_mask) != 0;
+        out = uint128_or(uint128_rshift(in, count), shifted_away);
+    } else {
+        out = ( in != 0 );
+    }
+    return out;
+}
+
 /*----------------------------------------------------------------------------
 | Shifts the 192-bit value formed by concatenating `a0', `a1', and `a2' right
 | by 64 _plus_ the number of bits given in `count'.  The shifted result is
diff --git a/fpu/softfloat.c b/fpu/softfloat.c
index e038434a07..694616bb3d 100644
--- a/fpu/softfloat.c
+++ b/fpu/softfloat.c
@@ -517,9 +517,20 @@ typedef struct {
     bool sign;
 } FloatParts;
 
+/* Similar for float128.  */
+typedef struct {
+    Uint128 frac;
+    int32_t exp;
+    FloatClass cls;
+    bool sign;
+} FloatParts128;
+
 #define DECOMPOSED_BINARY_POINT    (64 - 2)
 #define DECOMPOSED_IMPLICIT_BIT    (1ull << DECOMPOSED_BINARY_POINT)
 #define DECOMPOSED_OVERFLOW_BIT    (DECOMPOSED_IMPLICIT_BIT << 1)
+#define DECOMPOSED128_BINARY_POINT (128 - 2)
+#define DECOMPOSED128_IMPLICIT_BIT (uint128_lshift(1, DECOMPOSED128_BINARY_POINT))
+#define DECOMPOSED128_OVERFLOW_BIT (uint128_lshift(DECOMPOSED128_IMPLICIT_BIT, 1))
 
 /* Structure holding all of the relevant parameters for a format.
  *   exp_size: the size of the exponent field
@@ -559,6 +570,13 @@ typedef struct {
     .round_mask     = (1ull << (DECOMPOSED_BINARY_POINT - F)) - 1,   \
     .roundeven_mask = (2ull << (DECOMPOSED_BINARY_POINT - F)) - 1
 
+#define FLOAT128_PARAMS(E, F)                                           \
+    .exp_size       = E,                                             \
+    .exp_bias       = ((1 << E) - 1) >> 1,                           \
+    .exp_max        = (1 << E) - 1,                                  \
+    .frac_size      = F,                                             \
+    .frac_shift     = DECOMPOSED128_BINARY_POINT - F,
+
 static const FloatFmt float16_params = {
     FLOAT_PARAMS(5, 10)
 };
@@ -580,6 +598,16 @@ static const FloatFmt float64_params = {
     FLOAT_PARAMS(11, 52)
 };
 
+static const FloatFmt float128_params = {
+    FLOAT128_PARAMS(15, 112)
+    /* the remaining fields are not used for 128 bit ops */
+    .frac_lsb = 0,
+    .frac_lsbm1 = 0,
+    .round_mask = 0,
+    .roundeven_mask = 0,
+    .arm_althp = false
+};
+
 /* Unpack a float to parts, but do not canonicalize.  */
 static inline FloatParts unpack_raw(FloatFmt fmt, uint64_t raw)
 {
@@ -593,6 +621,18 @@ static inline FloatParts unpack_raw(FloatFmt fmt, uint64_t raw)
     };
 }
 
+static inline FloatParts128 unpack128_raw(FloatFmt fmt, Uint128 raw)
+{
+    const int sign_pos = fmt.frac_size + fmt.exp_size;
+
+    return (FloatParts128) {
+        .cls = float_class_unclassified,
+        .sign = extract128(raw, sign_pos, 1),
+        .exp = extract128(raw, fmt.frac_size, fmt.exp_size),
+        .frac = extract128(raw, 0, fmt.frac_size),
+    };
+}
+
 static inline FloatParts float16_unpack_raw(float16 f)
 {
     return unpack_raw(float16_params, f);
@@ -613,6 +653,11 @@ static inline FloatParts float64_unpack_raw(float64 f)
     return unpack_raw(float64_params, f);
 }
 
+static inline FloatParts128 float128_unpack_raw(float128 f)
+{
+    return unpack128_raw(float128_params, uint128_make128(f.low, f.high));
+}
+
 /* Pack a float from parts, but do not canonicalize.  */
 static inline uint64_t pack_raw(FloatFmt fmt, FloatParts p)
 {
@@ -621,6 +666,13 @@ static inline uint64_t pack_raw(FloatFmt fmt, FloatParts p)
     return deposit64(ret, sign_pos, 1, p.sign);
 }
 
+static inline Uint128 pack128_raw(FloatFmt fmt, FloatParts128 p)
+{
+    const int sign_pos = fmt.frac_size + fmt.exp_size;
+    Uint128 ret = deposit128(p.frac, fmt.frac_size, fmt.exp_size, p.exp);
+    return deposit128(ret, sign_pos, 1, p.sign);
+}
+
 static inline float16 float16_pack_raw(FloatParts p)
 {
     return make_float16(pack_raw(float16_params, p));
@@ -641,6 +693,12 @@ static inline float64 float64_pack_raw(FloatParts p)
     return make_float64(pack_raw(float64_params, p));
 }
 
+static inline float128 float128_pack_raw(FloatParts128 p)
+{
+    Uint128 out = pack128_raw(float128_params, p);
+    return make_float128(uint128_gethi(out), uint128_getlo(out));
+}
+
 /*----------------------------------------------------------------------------
 | Functions and definitions to determine:  (1) whether tininess for underflow
 | is detected before or after rounding by default, (2) what (if anything)
@@ -684,6 +742,41 @@ static FloatParts sf_canonicalize(FloatParts part, const FloatFmt *parm,
     return part;
 }
 
+/* Almost exactly the same as sf_canonicalize except 128 bit */
+static FloatParts128 sf128_canonicalize(FloatParts128 part, const FloatFmt *parm,
+                                        float_status *status)
+{
+    bool frac_is_zero = uint128_eq(part.frac, uint128_zero());
+    if (part.exp == parm->exp_max) {
+        if (frac_is_zero) {
+            part.cls = float_class_inf;
+        } else {
+            part.frac = uint128_lshift(part.frac, parm->frac_shift);
+            part.cls = (parts128_is_snan_frac(part.frac, status)
+                        ? float_class_snan : float_class_qnan);
+        }
+    } else if (part.exp == 0) {
+        if (likely(frac_is_zero)) {
+            part.cls = float_class_zero;
+        } else if (status->flush_inputs_to_zero) {
+            float_raise(float_flag_input_denormal, status);
+            part.cls = float_class_zero;
+            part.frac = uint128_zero();
+        } else {
+            int shift = clz128(part.frac) - 1;
+            part.cls = float_class_normal;
+            part.exp = parm->frac_shift - parm->exp_bias - shift + 1;
+            part.frac = uint128_lshift(part.frac, shift);
+        }
+    } else {
+        Uint128 imp_bit = uint128_lshift(1, DECOMPOSED128_BINARY_POINT);
+        part.cls = float_class_normal;
+        part.exp -= parm->exp_bias;
+        part.frac = uint128_add(imp_bit, uint128_lshift(part.frac, parm->frac_shift));
+    }
+    return part;
+}
+
 /* Round and uncanonicalize a floating-point number by parts. There
  * are FRAC_SHIFT bits that may require rounding at the bottom of the
  * fraction; these bits will be removed. The exponent will be biased
@@ -836,6 +929,144 @@ static FloatParts round_canonical(FloatParts p, float_status *s,
     return p;
 }
 
+/* As above but wider */
+static FloatParts128 round128_canonical(FloatParts128 p, float_status *s,
+                                        const FloatFmt *parm)
+{
+    /* Do these by hand rather than widening the FloatFmt structure */
+    const Uint128 frac_lsb = uint128_lshift(1, DECOMPOSED128_BINARY_POINT - parm->frac_size);
+    const Uint128 frac_lsbm1 = uint128_lshift(1, (DECOMPOSED128_BINARY_POINT - parm->frac_size) - 1);
+    const Uint128 round_mask = uint128_sub(frac_lsb, 1);
+    const Uint128 roundeven_mask = uint128_sub(uint128_lshift(2, DECOMPOSED128_BINARY_POINT - parm->frac_size), 1);
+    const int exp_max = parm->exp_max;
+    const int frac_shift = parm->frac_shift;
+    Uint128 frac, inc;
+    int exp, flags = 0;
+    bool overflow_norm;
+
+    frac = p.frac;
+    exp = p.exp;
+
+    switch (p.cls) {
+    case float_class_normal:
+        switch (s->float_rounding_mode) {
+        case float_round_nearest_even:
+            overflow_norm = false;
+            inc = ((frac & roundeven_mask) != frac_lsbm1 ? frac_lsbm1 : 0);
+            break;
+        case float_round_ties_away:
+            overflow_norm = false;
+            inc = frac_lsbm1;
+            break;
+        case float_round_to_zero:
+            overflow_norm = true;
+            inc = 0;
+            break;
+        case float_round_up:
+            inc = p.sign ? 0 : round_mask;
+            overflow_norm = p.sign;
+            break;
+        case float_round_down:
+            inc = p.sign ? round_mask : 0;
+            overflow_norm = !p.sign;
+            break;
+        case float_round_to_odd:
+            overflow_norm = true;
+            inc = frac & frac_lsb ? 0 : round_mask;
+            break;
+        default:
+            g_assert_not_reached();
+        }
+
+        exp += parm->exp_bias;
+        if (likely(exp > 0)) {
+            if (frac & round_mask) {
+                flags |= float_flag_inexact;
+                frac += inc;
+                if (frac & DECOMPOSED128_OVERFLOW_BIT) {
+                    frac >>= 1;
+                    exp++;
+                }
+            }
+            frac >>= frac_shift;
+
+            if (unlikely(exp >= exp_max)) {
+                flags |= float_flag_overflow | float_flag_inexact;
+                if (overflow_norm) {
+                    exp = exp_max - 1;
+                    frac = -1;
+                } else {
+                    p.cls = float_class_inf;
+                    goto do_inf;
+                }
+            }
+        } else if (s->flush_to_zero) {
+            flags |= float_flag_output_denormal;
+            p.cls = float_class_zero;
+            goto do_zero;
+        } else {
+            bool is_tiny = s->tininess_before_rounding
+                        || (exp < 0)
+                        || !(uint128_and(uint128_add(frac, inc), DECOMPOSED128_OVERFLOW_BIT));
+
+            frac = shift128RightJammingUint128(frac, 1 - exp);
+            if (frac & round_mask) {
+                /* Need to recompute round-to-even.  */
+                switch (s->float_rounding_mode) {
+                case float_round_nearest_even:
+                    inc = ((uint128_and(frac, roundeven_mask)) != frac_lsbm1
+                           ? frac_lsbm1 : 0);
+                    break;
+                case float_round_to_odd:
+                    inc = uint128_and(frac, frac_lsb ? 0 : round_mask);
+                    break;
+                default:
+                    break;
+                }
+                flags |= float_flag_inexact;
+                uint128_add(frac, inc);
+            }
+
+            exp = (uint128_and(frac, DECOMPOSED128_IMPLICIT_BIT) ? 1 : 0);
+            frac = uint128_rshift(frac, parm->frac_shift);
+
+            if (is_tiny && (flags & float_flag_inexact)) {
+                flags |= float_flag_underflow;
+            }
+            if (exp == 0 && frac == 0) {
+                p.cls = float_class_zero;
+            }
+        }
+        break;
+
+    case float_class_zero:
+    do_zero:
+        exp = 0;
+        frac = 0;
+        break;
+
+    case float_class_inf:
+    do_inf:
+        exp = exp_max;
+        frac = 0;
+        break;
+
+    case float_class_qnan:
+    case float_class_snan:
+        exp = exp_max;
+        frac = uint128_rshift(frac, parm->frac_shift);
+        break;
+
+    default:
+        g_assert_not_reached();
+    }
+
+    float_raise(flags, s);
+    p.exp = exp;
+    p.frac = frac;
+    return p;
+}
+
 /* Explicit FloatFmt version */
 static FloatParts float16a_unpack_canonical(float16 f, float_status *s,
                                             const FloatFmt *params)
@@ -889,6 +1120,16 @@ static float64 float64_round_pack_canonical(FloatParts p, float_status *s)
     return float64_pack_raw(round_canonical(p, s, &float64_params));
 }
 
+static FloatParts128 float128_unpack_canonical(float128 f, float_status *s)
+{
+    return sf128_canonicalize(float128_unpack_raw(f), &float128_params, s);
+}
+
+static float128 float128_round_pack_canonical(FloatParts128 p, float_status *s)
+{
+    return float128_pack_raw(round128_canonical(p, s, &float128_params));
+}
+
 static FloatParts return_nan(FloatParts a, float_status *s)
 {
     switch (a.cls) {
@@ -1018,6 +1259,110 @@ static FloatParts addsub_floats(FloatParts a, FloatParts b, bool subtract,
     g_assert_not_reached();
 }
 
+/* As above but with Uint128 fractions */
+static FloatParts128 pick_nan128(FloatParts128 a, FloatParts128 b, float_status *s)
+{
+    if (is_snan(a.cls) || is_snan(b.cls)) {
+        s->float_exception_flags |= float_flag_invalid;
+    }
+
+    if (s->default_nan_mode) {
+        return parts128_default_nan(s);
+    } else {
+        if (pickNaN(a.cls, b.cls,
+                    a.frac > b.frac ||
+                    (a.frac == b.frac && a.sign < b.sign), s)) {
+            a = b;
+        }
+        if (is_snan(a.cls)) {
+            return parts128_silence_nan(a, s);
+        }
+    }
+    return a;
+}
+
+static FloatParts128 addsub_floats128(FloatParts128 a, FloatParts128 b, bool subtract,
+                                   float_status *s)
+{
+    bool a_sign = a.sign;
+    bool b_sign = b.sign ^ subtract;
+
+    if (a_sign != b_sign) {
+        /* Subtraction */
+
+        if (a.cls == float_class_normal && b.cls == float_class_normal) {
+            if (a.exp > b.exp || (a.exp == b.exp && a.frac >= b.frac)) {
+                b.frac = shift128RightJammingUint128(b.frac, a.exp - b.exp);
+                a.frac = uint128_sub(a.frac, b.frac);
+            } else {
+                a.frac = shift128RightJammingUint128(a.frac, b.exp - a.exp);
+                a.frac = uint128_sub(b.frac, a.frac);
+                a.exp = b.exp;
+                a_sign ^= 1;
+            }
+
+            if (a.frac == 0) {
+                a.cls = float_class_zero;
+                a.sign = s->float_rounding_mode == float_round_down;
+            } else {
+                int shift = clz128(a.frac) - 1;
+                a.frac = uint128_lshift(a.frac, shift);
+                a.exp = a.exp - shift;
+                a.sign = a_sign;
+            }
+            return a;
+        }
+        if (is_nan(a.cls) || is_nan(b.cls)) {
+            return pick_nan128(a, b, s);
+        }
+        if (a.cls == float_class_inf) {
+            if (b.cls == float_class_inf) {
+                float_raise(float_flag_invalid, s);
+                return parts128_default_nan(s);
+            }
+            return a;
+        }
+        if (a.cls == float_class_zero && b.cls == float_class_zero) {
+            a.sign = s->float_rounding_mode == float_round_down;
+            return a;
+        }
+        if (a.cls == float_class_zero || b.cls == float_class_inf) {
+            b.sign = a_sign ^ 1;
+            return b;
+        }
+        if (b.cls == float_class_zero) {
+            return a;
+        }
+    } else {
+        /* Addition */
+        if (a.cls == float_class_normal && b.cls == float_class_normal) {
+            if (a.exp > b.exp) {
+                b.frac = shift128RightJammingUint128(b.frac, a.exp - b.exp);
+            } else if (a.exp < b.exp) {
+                a.frac = shift128RightJammingUint128(a.frac, b.exp - a.exp);
+                a.exp = b.exp;
+            }
+            a.frac = uint128_add(a.frac, b.frac);
+            if (a.frac & DECOMPOSED128_OVERFLOW_BIT) {
+                a.frac = shift128RightJammingUint128(a.frac, 1);
+                a.exp += 1;
+            }
+            return a;
+        }
+        if (is_nan(a.cls) || is_nan(b.cls)) {
+            return pick_nan128(a, b, s);
+        }
+        if (a.cls == float_class_inf || b.cls == float_class_zero) {
+            return a;
+        }
+        if (b.cls == float_class_inf || a.cls == float_class_zero) {
+            b.sign = b_sign;
+            return b;
+        }
+    }
+    g_assert_not_reached();
+}
+
 /*
  * Returns the result of adding or subtracting the floating-point
  * values `a' and `b'. The operation is performed according to the
@@ -1157,6 +1502,26 @@ float64_sub(float64 a, float64 b, float_status *s)
     return float64_addsub(a, b, s, hard_f64_sub, soft_f64_sub);
 }
 
+float128 QEMU_FLATTEN
+float128_add(float128 a, float128 b, float_status *status)
+{
+    FloatParts128 pa = float128_unpack_canonical(a, status);
+    FloatParts128 pb = float128_unpack_canonical(b, status);
+    FloatParts128 pr = addsub_floats128(pa, pb, false, status);
+
+    return float128_round_pack_canonical(pr, status);
+}
+
+float128 QEMU_FLATTEN
+float128_sub(float128 a, float128 b, float_status *status)
+{
+    FloatParts128 pa = float128_unpack_canonical(a, status);
+    FloatParts128 pb = float128_unpack_canonical(b, status);
+    FloatParts128 pr = addsub_floats128(pa, pb, true, status);
+
+    return float128_round_pack_canonical(pr, status);
+}
+
 /*
  * Returns the result of adding or subtracting the bfloat16
  * values `a' and `b'.
@@ -6921,227 +7286,6 @@ float128 float128_round_to_int(float128 a, float_status *status)
 
 }
 
-/*----------------------------------------------------------------------------
-| Returns the result of adding the absolute values of the quadruple-precision
-| floating-point values `a' and `b'.  If `zSign' is 1, the sum is negated
-| before being returned.  `zSign' is ignored if the result is a NaN.
-| The addition is performed according to the IEC/IEEE Standard for Binary
-| Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-static float128 addFloat128Sigs(float128 a, float128 b, bool zSign,
-                                float_status *status)
-{
-    int32_t aExp, bExp, zExp;
-    uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1, zSig2;
-    int32_t expDiff;
-
-    aSig1 = extractFloat128Frac1( a );
-    aSig0 = extractFloat128Frac0( a );
-    aExp = extractFloat128Exp( a );
-    bSig1 = extractFloat128Frac1( b );
-    bSig0 = extractFloat128Frac0( b );
-    bExp = extractFloat128Exp( b );
-    expDiff = aExp - bExp;
-    if ( 0 < expDiff ) {
-        if ( aExp == 0x7FFF ) {
-            if (aSig0 | aSig1) {
-                return propagateFloat128NaN(a, b, status);
-            }
-            return a;
-        }
-        if ( bExp == 0 ) {
-            --expDiff;
-        }
-        else {
-            bSig0 |= UINT64_C(0x0001000000000000);
-        }
-        shift128ExtraRightJamming(
-            bSig0, bSig1, 0, expDiff, &bSig0, &bSig1, &zSig2 );
-        zExp = aExp;
-    }
-    else if ( expDiff < 0 ) {
-        if ( bExp == 0x7FFF ) {
-            if (bSig0 | bSig1) {
-                return propagateFloat128NaN(a, b, status);
-            }
-            return packFloat128( zSign, 0x7FFF, 0, 0 );
-        }
-        if ( aExp == 0 ) {
-            ++expDiff;
-        }
-        else {
-            aSig0 |= UINT64_C(0x0001000000000000);
-        }
-        shift128ExtraRightJamming(
-            aSig0, aSig1, 0, - expDiff, &aSig0, &aSig1, &zSig2 );
-        zExp = bExp;
-    }
-    else {
-        if ( aExp == 0x7FFF ) {
-            if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
-                return propagateFloat128NaN(a, b, status);
-            }
-            return a;
-        }
-        add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
-        if ( aExp == 0 ) {
-            if (status->flush_to_zero) {
-                if (zSig0 | zSig1) {
-                    float_raise(float_flag_output_denormal, status);
-                }
-                return packFloat128(zSign, 0, 0, 0);
-            }
-            return packFloat128( zSign, 0, zSig0, zSig1 );
-        }
-        zSig2 = 0;
-        zSig0 |= UINT64_C(0x0002000000000000);
-        zExp = aExp;
-        goto shiftRight1;
-    }
-    aSig0 |= UINT64_C(0x0001000000000000);
-    add128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
-    --zExp;
-    if ( zSig0 < UINT64_C(0x0002000000000000) ) goto roundAndPack;
-    ++zExp;
- shiftRight1:
-    shift128ExtraRightJamming(
-        zSig0, zSig1, zSig2, 1, &zSig0, &zSig1, &zSig2 );
- roundAndPack:
-    return roundAndPackFloat128(zSign, zExp, zSig0, zSig1, zSig2, status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of subtracting the absolute values of the quadruple-
-| precision floating-point values `a' and `b'.  If `zSign' is 1, the
-| difference is negated before being returned.  `zSign' is ignored if the
-| result is a NaN.  The subtraction is performed according to the IEC/IEEE
-| Standard for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-static float128 subFloat128Sigs(float128 a, float128 b, bool zSign,
-                                float_status *status)
-{
-    int32_t aExp, bExp, zExp;
-    uint64_t aSig0, aSig1, bSig0, bSig1, zSig0, zSig1;
-    int32_t expDiff;
-
-    aSig1 = extractFloat128Frac1( a );
-    aSig0 = extractFloat128Frac0( a );
-    aExp = extractFloat128Exp( a );
-    bSig1 = extractFloat128Frac1( b );
-    bSig0 = extractFloat128Frac0( b );
-    bExp = extractFloat128Exp( b );
-    expDiff = aExp - bExp;
-    shortShift128Left( aSig0, aSig1, 14, &aSig0, &aSig1 );
-    shortShift128Left( bSig0, bSig1, 14, &bSig0, &bSig1 );
-    if ( 0 < expDiff ) goto aExpBigger;
-    if ( expDiff < 0 ) goto bExpBigger;
-    if ( aExp == 0x7FFF ) {
-        if ( aSig0 | aSig1 | bSig0 | bSig1 ) {
-            return propagateFloat128NaN(a, b, status);
-        }
-        float_raise(float_flag_invalid, status);
-        return float128_default_nan(status);
-    }
-    if ( aExp == 0 ) {
-        aExp = 1;
-        bExp = 1;
-    }
-    if ( bSig0 < aSig0 ) goto aBigger;
-    if ( aSig0 < bSig0 ) goto bBigger;
-    if ( bSig1 < aSig1 ) goto aBigger;
-    if ( aSig1 < bSig1 ) goto bBigger;
-    return packFloat128(status->float_rounding_mode == float_round_down,
-                        0, 0, 0);
- bExpBigger:
-    if ( bExp == 0x7FFF ) {
-        if (bSig0 | bSig1) {
-            return propagateFloat128NaN(a, b, status);
-        }
-        return packFloat128( zSign ^ 1, 0x7FFF, 0, 0 );
-    }
-    if ( aExp == 0 ) {
-        ++expDiff;
-    }
-    else {
-        aSig0 |= UINT64_C(0x4000000000000000);
-    }
-    shift128RightJamming( aSig0, aSig1, - expDiff, &aSig0, &aSig1 );
-    bSig0 |= UINT64_C(0x4000000000000000);
- bBigger:
-    sub128( bSig0, bSig1, aSig0, aSig1, &zSig0, &zSig1 );
-    zExp = bExp;
-    zSign ^= 1;
-    goto normalizeRoundAndPack;
- aExpBigger:
-    if ( aExp == 0x7FFF ) {
-        if (aSig0 | aSig1) {
-            return propagateFloat128NaN(a, b, status);
-        }
-        return a;
-    }
-    if ( bExp == 0 ) {
-        --expDiff;
-    }
-    else {
-        bSig0 |= UINT64_C(0x4000000000000000);
-    }
-    shift128RightJamming( bSig0, bSig1, expDiff, &bSig0, &bSig1 );
-    aSig0 |= UINT64_C(0x4000000000000000);
- aBigger:
-    sub128( aSig0, aSig1, bSig0, bSig1, &zSig0, &zSig1 );
-    zExp = aExp;
- normalizeRoundAndPack:
-    --zExp;
-    return normalizeRoundAndPackFloat128(zSign, zExp - 14, zSig0, zSig1,
-                                         status);
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of adding the quadruple-precision floating-point values
-| `a' and `b'.  The operation is performed according to the IEC/IEEE Standard
-| for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_add(float128 a, float128 b, float_status *status)
-{
-    bool aSign, bSign;
-
-    aSign = extractFloat128Sign( a );
-    bSign = extractFloat128Sign( b );
-    if ( aSign == bSign ) {
-        return addFloat128Sigs(a, b, aSign, status);
-    }
-    else {
-        return subFloat128Sigs(a, b, aSign, status);
-    }
-
-}
-
-/*----------------------------------------------------------------------------
-| Returns the result of subtracting the quadruple-precision floating-point
-| values `a' and `b'.  The operation is performed according to the IEC/IEEE
-| Standard for Binary Floating-Point Arithmetic.
-*----------------------------------------------------------------------------*/
-
-float128 float128_sub(float128 a, float128 b, float_status *status)
-{
-    bool aSign, bSign;
-
-    aSign = extractFloat128Sign( a );
-    bSign = extractFloat128Sign( b );
-    if ( aSign == bSign ) {
-        return subFloat128Sigs(a, b, aSign, status);
-    }
-    else {
-        return addFloat128Sigs(a, b, aSign, status);
-    }
-
-}
-
 /*----------------------------------------------------------------------------
 | Returns the result of multiplying the quadruple-precision floating-point
 | values `a' and `b'.  The operation is performed according to the IEC/IEEE
diff --git a/fpu/softfloat-specialize.c.inc b/fpu/softfloat-specialize.c.inc
index c2f87addb2..347896bf57 100644
--- a/fpu/softfloat-specialize.c.inc
+++ b/fpu/softfloat-specialize.c.inc
@@ -125,6 +125,16 @@ static bool parts_is_snan_frac(uint64_t frac, float_status *status)
     }
 }
 
+static bool parts128_is_snan_frac(Uint128 frac, float_status *status)
+{
+    if (no_signaling_nans(status)) {
+        return false;
+    } else {
+        bool msb = extract128(frac, DECOMPOSED128_BINARY_POINT - 1, 1);
+        return msb == snan_bit_is_one(status);
+    }
+}
+
 /*----------------------------------------------------------------------------
 | The pattern for a default generated deconstructed floating-point NaN.
 *----------------------------------------------------------------------------*/
@@ -169,6 +179,23 @@ static FloatParts parts_default_nan(float_status *status)
     };
 }
 
+static FloatParts128 parts128_default_nan(float_status *status)
+{
+    bool sign = 0;
+    Uint128 frac;
+
+    /* !snan_bit_is_one, set sign and msb */
+    frac = uint128_lshift(1, DECOMPOSED128_BINARY_POINT - 1);
+    sign = 1;
+
+    return (FloatParts128) {
+        .cls = float_class_qnan,
+        .sign = sign,
+        .exp = INT_MAX,
+        .frac = frac
+    };
+}
+
 /*----------------------------------------------------------------------------
 | Returns a quiet NaN from a signalling NaN for the deconstructed
 | floating-point parts.
@@ -191,6 +218,18 @@ static FloatParts parts_silence_nan(FloatParts a, float_status *status)
     return a;
 }
 
+static FloatParts128 parts128_silence_nan(FloatParts128 a, float_status *status)
+{
+    g_assert(!no_signaling_nans(status));
+    if (snan_bit_is_one(status)) {
+        return parts128_default_nan(status);
+    } else {
+        a.frac = uint128_or(a.frac, int128_lshift(1, DECOMPOSED128_BINARY_POINT - 1));
+    }
+    a.cls = float_class_qnan;
+    return a;
+}
+
 /*----------------------------------------------------------------------------
 | The pattern for a default generated extended double-precision NaN.
 *----------------------------------------------------------------------------*/
-- 
2.20.1

Re: [RFC PATCH 8/8] softfloat: implement addsub_floats128 using Uint128 and new style code

[Richard Henderson]

On 10/20/20 9:37 AM, Alex Bennée wrote:
> +static inline FloatParts128 unpack128_raw(FloatFmt fmt, Uint128 raw)
> +{
> +    const int sign_pos = fmt.frac_size + fmt.exp_size;
> +
> +    return (FloatParts128) {
> +        .cls = float_class_unclassified,
> +        .sign = extract128(raw, sign_pos, 1),
> +        .exp = extract128(raw, fmt.frac_size, fmt.exp_size),
> +        .frac = extract128(raw, 0, fmt.frac_size),
> +    };
> +}

This use of extract128 for sign and exp will not work for 32-bit.  You can't
just automatically truncate from __uint128_t to int in that case.

I don't think we should necessarily create this function, but rather leave it at

> +static inline FloatParts128 float128_unpack_raw(float128 f)
> +{
> +    return unpack128_raw(float128_params, uint128_make128(f.low, f.high));
> +}

... this one, and construct the FloatParts128 directly from the float128
components.  E.g.

    int f_size = float128_params.frac_size;
    int e_size = float128_params.exp_size;
    return (FloatParts128) {
       .sign = extract64(f.high, f_size + e_size - 64, 1);
       .exp = extract64(f.high, f_size - 64, e_size);
       .frac = extract128(int128_make128(f.low, f.high),
                          0, f_size);
    };

I don't want to over-generalize this just yet.


> +static inline Uint128 pack128_raw(FloatFmt fmt, FloatParts128 p)
> +{
> +    const int sign_pos = fmt.frac_size + fmt.exp_size;
> +    Uint128 ret = deposit128(p.frac, fmt.frac_size, fmt.exp_size, p.exp);
> +    return deposit128(ret, sign_pos, 1, p.sign);
> +}

Likewise, omit this and only have

> +static inline float128 float128_pack_raw(FloatParts128 p)
> +{
> +    Uint128 out = pack128_raw(float128_params, p);
> +    return make_float128(uint128_gethi(out), uint128_getlo(out));
> +}

this.


> +/* Almost exactly the same as sf_canonicalize except 128 bit */
> +static FloatParts128 sf128_canonicalize(FloatParts128 part, const FloatFmt *parm,
> +                                        float_status *status)

I think we may have reached the point of diminishing returns on structure
returns.  This is a 196-bit struct, and will not be passed in registers
anymore.  It might be better to do

static void sf128_canonicalize(FloatParts128 *part,
                               const FloatFmt *parm,
                               float_status *status)

and modify the FloatParts128 in place.

> +    bool frac_is_zero = uint128_eq(part.frac, uint128_zero());

With Int128, we'd use !int128_nz().

> +/* As above but wider */
> +static FloatParts128 round128_canonical(FloatParts128 p, float_status *s,
> +                                        const FloatFmt *parm)
> +{
> +    /* Do these by hand rather than widening the FloatFmt structure */
> +    const Uint128 frac_lsb = uint128_lshift(1, DECOMPOSED128_BINARY_POINT - parm->frac_size);

You can't pass constant 1 on 32-bit.
Maybe add a int128_makepow2(exp) function to make this easier?

> +        case float_round_nearest_even:
> +            overflow_norm = false;
> +            inc = ((frac & roundeven_mask) != frac_lsbm1 ? frac_lsbm1 : 0);

Can't use & or != on 32-bit.

> +            inc = frac & frac_lsb ? 0 : round_mask;
...
> +            if (frac & round_mask) {
...
> +                frac += inc;
> +                if (frac & DECOMPOSED128_OVERFLOW_BIT) {
> +                    frac >>= 1;
...
> +            frac >>= frac_shift;
...
> +                    frac = -1;
...
> +            if (frac & round_mask) {
> +                    inc = ((uint128_and(frac, roundeven_mask)) != frac_lsbm1
...
> +            if (exp == 0 && frac == 0) {
...
> +        frac = 0;
...
> +        frac = 0;

and more.  There are lots more later.

This is going to get ugly fast.  We need another solution.

> +static bool parts128_is_snan_frac(Uint128 frac, float_status *status)
> +{
> +    if (no_signaling_nans(status)) {
> +        return false;
> +    } else {
> +        bool msb = extract128(frac, DECOMPOSED128_BINARY_POINT - 1, 1);

Doesn't work for 32-bit.  Again, extract128 by itself is not the right
interface.  Do we in fact want to share code with the normal parts_is_snan_frac
by just passing in the high-part?


r~

Re: [RFC PATCH 0/8] fpu: experimental conversion of float128_addsub

[no-reply]

Patchew URL: https://patchew.org/QEMU/20201020163738.27700-1-alex.bennee@linaro.org/



Hi,

This series seems to have some coding style problems. See output below for
more information:

Type: series
Message-id: 20201020163738.27700-1-alex.bennee@linaro.org
Subject: [RFC PATCH  0/8] fpu: experimental conversion of float128_addsub

=== TEST SCRIPT BEGIN ===
#!/bin/bash
git rev-parse base > /dev/null || exit 0
git config --local diff.renamelimit 0
git config --local diff.renames True
git config --local diff.algorithm histogram
./scripts/checkpatch.pl --mailback base..
=== TEST SCRIPT END ===

Updating 3c8cf5a9c21ff8782164d1def7f44bd888713384
From https://github.com/patchew-project/qemu
 * [new tag]         patchew/20201020163738.27700-1-alex.bennee@linaro.org -> patchew/20201020163738.27700-1-alex.bennee@linaro.org
Switched to a new branch 'test'
165b2d1 softfloat: implement addsub_floats128 using Uint128 and new style code
cfb628a tests/fp: add quad support to the benchmark utility
593cf15 int128.h: add bunch of uint128 utility functions (INCOMPLETE)
83e37cd softfloat: Inline pick_nan_muladd into its caller
45f794e softfloat: Add float_cmask and constants
27be22d softfloat: Tidy a * b + inf return
e3c7e01 softfloat: Use int128.h for some operations
321df84 softfloat: Use mulu64 for mul64To128

=== OUTPUT BEGIN ===
1/8 Checking commit 321df844302f (softfloat: Use mulu64 for mul64To128)
2/8 Checking commit e3c7e01abc34 (softfloat: Use int128.h for some operations)
3/8 Checking commit 27be22d0ba4c (softfloat: Tidy a * b + inf return)
4/8 Checking commit 45f794e02b18 (softfloat: Add float_cmask and constants)
5/8 Checking commit 83e37cdf6e25 (softfloat: Inline pick_nan_muladd into its caller)
6/8 Checking commit 593cf1552545 (int128.h: add bunch of uint128 utility functions (INCOMPLETE))
7/8 Checking commit cfb628ae8542 (tests/fp: add quad support to the benchmark utility)
WARNING: line over 80 characters
#177: FILE: tests/fp/fp-bench.c:458:
+    GEN_BENCH_NO_NEG(bench_ ## name ## _float128, float128, PREC_FLOAT128, op, n)

WARNING: line over 80 characters
#194: FILE: tests/fp/fp-bench.c:521:
+    fprintf(stderr, " -p = floating point precision (single, double, quad[soft only]). "

total: 0 errors, 2 warnings, 185 lines checked

Patch 7/8 has style problems, please review.  If any of these errors
are false positives report them to the maintainer, see
CHECKPATCH in MAINTAINERS.
8/8 Checking commit 165b2d1b1292 (softfloat: implement addsub_floats128 using Uint128 and new style code)
WARNING: line over 80 characters
#101: FILE: fpu/softfloat.c:532:
+#define DECOMPOSED128_IMPLICIT_BIT (uint128_lshift(1, DECOMPOSED128_BINARY_POINT))

WARNING: line over 80 characters
#102: FILE: fpu/softfloat.c:533:
+#define DECOMPOSED128_OVERFLOW_BIT (uint128_lshift(DECOMPOSED128_IMPLICIT_BIT, 1))

WARNING: line over 80 characters
#200: FILE: fpu/softfloat.c:746:
+static FloatParts128 sf128_canonicalize(FloatParts128 part, const FloatFmt *parm,

WARNING: line over 80 characters
#229: FILE: fpu/softfloat.c:775:
+        part.frac = uint128_add(imp_bit, uint128_lshift(part.frac, parm->frac_shift));

ERROR: line over 90 characters
#246: FILE: fpu/softfloat.c:937:
+    const Uint128 frac_lsb = uint128_lshift(1, DECOMPOSED128_BINARY_POINT - parm->frac_size);

ERROR: line over 90 characters
#247: FILE: fpu/softfloat.c:938:
+    const Uint128 frac_lsbm1 = uint128_lshift(1, (DECOMPOSED128_BINARY_POINT - parm->frac_size) - 1);

ERROR: line over 90 characters
#249: FILE: fpu/softfloat.c:940:
+    const Uint128 roundeven_mask = uint128_sub(uint128_lshift(2, DECOMPOSED128_BINARY_POINT - parm->frac_size), 1);

ERROR: line over 90 characters
#319: FILE: fpu/softfloat.c:1010:
+                        || !(uint128_and(uint128_add(frac, inc), DECOMPOSED128_OVERFLOW_BIT));

WARNING: line over 80 characters
#404: FILE: fpu/softfloat.c:1263:
+static FloatParts128 pick_nan128(FloatParts128 a, FloatParts128 b, float_status *s)

WARNING: line over 80 characters
#425: FILE: fpu/softfloat.c:1284:
+static FloatParts128 addsub_floats128(FloatParts128 a, FloatParts128 b, bool subtract,

ERROR: space prohibited after that open parenthesis '('
#780: FILE: include/fpu/softfloat-macros.h:269:
+    } else if ( count < 128 ) {

ERROR: space prohibited before that close parenthesis ')'
#780: FILE: include/fpu/softfloat-macros.h:269:
+    } else if ( count < 128 ) {

ERROR: space prohibited after that open parenthesis '('
#785: FILE: include/fpu/softfloat-macros.h:274:
+        out = ( in != 0 );

ERROR: space prohibited before that close parenthesis ')'
#785: FILE: include/fpu/softfloat-macros.h:274:
+        out = ( in != 0 );

total: 8 errors, 6 warnings, 747 lines checked

Patch 8/8 has style problems, please review.  If any of these errors
are false positives report them to the maintainer, see
CHECKPATCH in MAINTAINERS.

=== OUTPUT END ===

Test command exited with code: 1


The full log is available at
http://patchew.org/logs/20201020163738.27700-1-alex.bennee@linaro.org/testing.checkpatch/?type=message.
---
Email generated automatically by Patchew [https://patchew.org/].
Please send your feedback to patchew-devel@redhat.com