ppc: softfloat float implementation issues

[Paul Clarke]

Per bug #1841491, Richard Henderson (rth) said:
> The float test failure is part of a larger problem for target/powerpc
> in which all float routines are implemented incorrectly. They are all
> implemented as double operations with rounding to float as a second
> step. Which not only produces incorrect exceptions, as in this case,
> but incorrect numerical results from the double rounding.
>
> This should probably be split to a separate bug...

[Paul Clarke (7-pc)]

-- ppc64le native:
$ gcc -c -O2 ffma.c
$ gcc -O2 test-ffma.c ffma.o -lm -o test-ffma
$ ./test-ffma $(./test-ffma)
ffma(0x1p-149, 0x1p-149, 0x1p-149)
0x0

0xa000000
FE_INEXACT FE_UNDERFLOW
0x1p-149

-- qemu-system-ppc64:
$ ./test-ffma $(./test-ffma)
ffma(0x1p-149, 0x1p-149, 0x1p-149)
0x0

0x2000000
FE_INEXACT
0x1p-149

[Alex Bennée (ajbennee)]

I'm confused by this testcase as it's not a fused multiply-add but as you say two combined operations.

[Richard Henderson (rth)]

It should be a fused multiply add; you may need to use -ffast-math or
something to get the compiler to generate the proper instruction.

However, one can see from target/ppc/translate/fp-impl.inc.c:

/* fmadd - fmadds */
GEN_FLOAT_ACB(madd, 0x1D, 1, PPC_FLOAT);

through to _GEN_FLOAT_ACB:

    gen_helper_f##op(t3, cpu_env, t0, t1, t2); \
    if (isfloat) { \
        gen_helper_frsp(t3, cpu_env, t3); \
    } \

That right there is a double-precision fma followed by a round
to single precision. This pattern is replicated for all single
precision operations, and is of course wrong.

I believe that correct results may be obtained by having
single-precision helpers that first convert the double-precision
input into a single-precision input using helper_tosingle(),
perform the required operation, then convert the result back to
double-precision using helper_todouble().

The manual says:

# For single-precision arithmetic instructions, all input values
# must be representable in single format; if they are not, the
# result placed into the target FPR, and the setting of
# status bits in the FPSCR and in the Condition Register
# (if Rc=1), are undefined.

The tosingle/todouble conversions are exact and bit-preserving.
They are used by load-single and store-single that convert a
single-precision in-memory value to the double-precision register
value. Therefore the input given to float32_add using this
conversion would be exactly the same as if we had given the
value unmollested from a memory input.

I don't know what real ppc hw does -- whether it takes all of
the double-precision input bits and rounds to 23-bits, like the
old 80387 hardware does, or truncates the input as I propose.
But for architectural results we don't have to care, because
of the UNDEFINED escape clause.

[Alex Bennée (ajbennee)]

Testing on current master shows the behavior is correct. I guess rth's patch fixed this case.

[Alex Bennée (ajbennee)]

It looks like the test case isn't properly exercising the code that is likely to be wrong. It sounds like we need a proper comprehensive testcase for fused operations (along the line of the ARM fcvt test case). This can probably be a multiarch testcase which we can build for all the various targets.

[Alex Bennée (ajbennee)]

This is a generic floating point multiply and accumulate test for
single precision floating point values. I've split of the common float
functions into a helper library so additional tests can use the same
common code.

Signed-off-by: Alex Bennée <email address hidden>
---
 tests/tcg/multiarch/Makefile.target | 7 +-
 tests/tcg/multiarch/float_helpers.c | 208 ++++++++++++++++++++++++++++
 tests/tcg/multiarch/float_helpers.h | 26 ++++
 tests/tcg/multiarch/float_madds.c | 78 +++++++++++
 4 files changed, 318 insertions(+), 1 deletion(-)
 create mode 100644 tests/tcg/multiarch/float_helpers.c
 create mode 100644 tests/tcg/multiarch/float_helpers.h
 create mode 100644 tests/tcg/multiarch/float_madds.c

diff --git a/tests/tcg/multiarch/Makefile.target b/tests/tcg/multiarch/Makefile.target
index 657a04f802d..0446b75c456 100644
--- a/tests/tcg/multiarch/Makefile.target
+++ b/tests/tcg/multiarch/Makefile.target
@@ -10,12 +10,17 @@ MULTIARCH_SRC=$(SRC_PATH)/tests/tcg/multiarch
 # Set search path for all sources
 VPATH += $(MULTIARCH_SRC)
 MULTIARCH_SRCS =$(notdir $(wildcard $(MULTIARCH_SRC)/*.c))
-MULTIARCH_TESTS =$(MULTIARCH_SRCS:.c=)
+MULTIARCH_TESTS =$(filter-out float_helpers, $(MULTIARCH_SRCS:.c=))

 #
 # The following are any additional rules needed to build things
 #

+
+float_madds: LDFLAGS+=-lm
+float_madds: float_madds.c float_helpers.c
+ $(CC) $(CFLAGS) $(EXTRA_CFLAGS) -O2 $< $(MULTIARCH_SRC)/float_helpers.c -o $@ $(LDFLAGS)
+
 testthread: LDFLAGS+=-lpthread

 # We define the runner for test-mmap after the individual
diff --git a/tests/tcg/multiarch/float_helpers.c b/tests/tcg/multiarch/float_helpers.c
new file mode 100644
index 00000000000..481d8d33317
--- /dev/null
+++ b/tests/tcg/multiarch/float_helpers.c
@@ -0,0 +1,208 @@
+/*
+ * Common Float Helpers
+ *
+ * This contains a series of useful utility routines and a set of
+ * floating point constants useful for exercising the edge cases in
+ * floating point tests.
+ *
+ * Copyright (c) 2019 Linaro
+ *
+ * SPDX-License-Identifier: GPL-3.0-or-later
+ */
+
+/* we want additional float type definitions */
+#define __STDC_WANT_IEC_60559_BFP_EXT__
+#define __STDC_WANT_IEC_60559_TYPES_EXT__
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <inttypes.h>
+#include <math.h>
+#include <float.h>
+#include <fenv.h>
+
+#include "float_helpers.h"
+
+#define ARRAY_SIZE(x) (sizeof(x) / sizeof((x)[0]))
+
+/*
+ * Half Precision Numbers
+ *
+ * Not yet well standardised so we return a plain uint16_t for now.
+ */
+
+/* no handy defines for these numbers */
+static uint16_t f16_numbers[] = {
+ 0xffff, /* -NaN / AHP -Max */
+ 0xfcff, /* -NaN / AHP */
+ 0xfc01, /* -NaN / AHP */
+ 0xfc00, /* -Inf */
+ 0xfbff, /* -Max */
+ 0xc000, /* -2 */
+ 0xbc00, /* -1 */
+ 0x8001, /* -MIN subnormal */
+ 0x8000, /* -0 */
+ 0x0000, /* +0 */
+ 0x0001, /* MIN subnormal */
+ 0x3c00, /* 1 */
+ 0x7bff, /* Max */
+ 0x7c00, /* Inf */
+ 0x7c01, /* NaN / AHP */
+ 0x7cff, /* NaN / AHP */
+ 0x7fff, /* NaN / AHP +Max*/
+};
+
+const int num_f16 = ARRAY_SIZE(f16_numbers);
+
+uint16_t get_f16(int i) {
+ return f16_numbers[i % num_f16];
+}
+
+/* only display as hex */
+c...

[Alex Bennée (ajbennee)]

Dave ran the above testcase on:

> processor : 0
> cpu : POWER8E (raw), altivec supported
> clock : 3325.000000MHz
> revision : 2.1 (pvr 004b 0201)

And there are no diffs with what you currently get from master. So I think the invalid flag fault is fixed by a previous commit and the potential precision faults don't get picked up by the test case. I guess we could be a bit smarted about testing the limits.

[Paul Clarke (7-pc)]

On 9/13/19 8:49 AM, Alex Bennée wrote:
> +static float f32_numbers[] = {
> + -SNANF,
> + -NAN,
> + -INFINITY,
> + -FLT_MAX,
> + -1.111E+31,
> + -1.111E+30,
> + -1.08700982e-12,
> + -1.78051176e-20,
> + -FLT_MIN,
> + 0.0,
> + FLT_MIN,
> + 2.98023224e-08,
> + 5.96046E-8, /* min positive FP16 subnormal */
> + 6.09756E-5, /* max subnormal FP16 */
> + 6.10352E-5, /* min positive normal FP16 */
> + 1.0,
> + 1.0009765625, /* smallest float after 1.0 FP16 */
> + 2.0,
> + M_E, M_PI,
> + 65503.0,
> + 65504.0, /* max FP16 */
> + 65505.0,
> + 131007.0,
> + 131008.0, /* max AFP */
> + 131009.0,
> + 1.111E+30,
> + FLT_MAX,
> + INFINITY,
> + NAN,
> + SNANF
> +};

I've noticed that Glibc prefers to use hex representation for float values, to ensure an accurate representation. If you care to do so, here they are:
static float f32_numbers[] = {
    -SNANF,
    -NAN,
    -INFINITY,
    -FLT_MAX,
    -0x1.1874b2p+103,
    -0x1.c0bab6p+99,
    -0x1.31f75p-40,
    -0x1.505444p-66,
    -FLT_MIN,
    0.0,
    FLT_MIN,
    0x1p-25,
    0x1.ffffe6p-25, /* min positive FP16 subnormal */
    0x1.ff801ap-15, /* max subnormal FP16 */
    0x1.00000cp-14, /* min positive normal FP16 */
    1.0,
    0x1.004p+0, /* smallest float after 1.0 FP16 */
    2.0,
    M_E, M_PI,
    0x1.ffbep+15,
    0x1.ffcp+15, /* max FP16 */
    0x1.ffc2p+15,
    0x1.ffbfp+16,
    0x1.ffcp+16, /* max AFP */
    0x1.ffc1p+16,
    0x1.c0bab6p+99,
    FLT_MAX,
    INFINITY,
    NAN,
    SNANF
};

PC

[Alex Bennée (ajbennee)]

Richard Henderson <email address hidden> writes:

> On 9/13/19 9:49 AM, Alex Bennée wrote:
>> + /* must be built with -O2 to generate fused op */
>> + r = a * b + c;
>
> I would prefer to use fmaf() or __builtin_fmaf() here.
>
> Although you'll need to link with -lm just in case the
> target doesn't support an instruction for fmaf and so
> the builtin expands to the standard library call.

I can do that - we have other tests that link to libm.

I was expecting to see more breakage but the ppc64 tests all passed (or
at least against the power8 David ran it on). What am I missing to hit
the cases you know are broken?

I've also experimented with reducing the number of iterations because if
we want to have golden references we probably don't want to dump several
hundred kilobytes of "golden" references into the source tree.

> I also like Paul's suggestion to use hex float constants.

Hmm I guess - look a bit weird but I guess that's lack of familiarity.
Is is still normalised? I guess the frac shows up (without the implicit
bit).

>
>
> r~

--
Alex Bennée

[Launchpad Janitor (janitor)]

[Expired for QEMU because there has been no activity for 60 days.]