softfloat: Handle m68k extended precision denormals properly

Motorola treats denormals with explicit integer bit set as
having unbiased exponent 0, unlike Intel which treats it as
having unbiased exponent 1 (like all other IEEE formats).

Add a flag on FloatFmt to differentiate the behaviour.

Reported-by: Keith Packard <keithp@keithp.com>
Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
---
 fpu/softfloat.c                |  9 +++++-
 tests/tcg/m68k/denormal.c      | 53 ++++++++++++++++++++++++++++++++++
 fpu/softfloat-parts.c.inc      |  5 ++--
 tests/tcg/m68k/Makefile.target |  2 +-
 4 files changed, 65 insertions(+), 4 deletions(-)
 create mode 100644 tests/tcg/m68k/denormal.c

> Motorola treats denormals with explicit integer bit set as
> having unbiased exponent 0, unlike Intel which treats it as
> having unbiased exponent 1 (like all other IEEE formats).

Thanks for having a look at this. Your patch fixes a couple of cases,
but there are further adventures that await if you're interested.

           x:  0x1p0                      0x3fff 0x80000000 0x00000000
           y:  0x1p-16383                 0x0000 0x80000000 0x00000000
   build_mul:  0x1p-16382                 0x0000 0x80000000 0x00000000
 runtime_mul:  0x1p-16383                 0x0001 0x80000000 0x00000000

I think the enclosed additional patch fixes this. I've still got 75 fmal
failures on this target, but the obvious 'multiply is broken' problem
appears fixed.

From b722c92f8329f56f5243496eca3779f1156aff4f Mon Sep 17 00:00:00 2001
From: Keith Packard <keithp@keithp.com>
Date: Sun, 20 Aug 2023 18:20:13 -0700
Subject: [PATCH] softfloat: Handle m68k LDBL_MIN_EXP normal values

Unlike Intel 80-bit floats, Motorola allows for normal values with a
zero exponent. Handle that by not setting exponent to 1 when the value
is normal for this format.

Signed-off-by: Keith Packard <keithp@keithp.com>
---
 fpu/softfloat-parts.c.inc | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc
index d0c43c28fb..cea854cdf1 100644
--- a/fpu/softfloat-parts.c.inc
+++ b/fpu/softfloat-parts.c.inc
@@ -288,7 +288,7 @@ static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s,
             p->frac_lo &= ~round_mask;
         }
 
-        exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) != 0;
+        exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) != 0 && !fmt->m68k_denormal;
         frac_shr(p, frac_shift);
 
         if (is_tiny && (flags & float_flag_inexact)) {

On 8/20/23 18:24, Keith Packard wrote:
> 
>> Motorola treats denormals with explicit integer bit set as
>> having unbiased exponent 0, unlike Intel which treats it as
>> having unbiased exponent 1 (like all other IEEE formats).
> 
> Thanks for having a look at this. Your patch fixes a couple of cases,
> but there are further adventures that await if you're interested.
> 
>             x:  0x1p0                      0x3fff 0x80000000 0x00000000
>             y:  0x1p-16383                 0x0000 0x80000000 0x00000000
>     build_mul:  0x1p-16382                 0x0000 0x80000000 0x00000000
>   runtime_mul:  0x1p-16383                 0x0001 0x80000000 0x00000000
> 
> I think the enclosed additional patch fixes this. I've still got 75 fmal
> failures on this target, but the obvious 'multiply is broken' problem
> appears fixed.
> 
>  From b722c92f8329f56f5243496eca3779f1156aff4f Mon Sep 17 00:00:00 2001
> From: Keith Packard <keithp@keithp.com>
> Date: Sun, 20 Aug 2023 18:20:13 -0700
> Subject: [PATCH] softfloat: Handle m68k LDBL_MIN_EXP normal values
> 
> Unlike Intel 80-bit floats, Motorola allows for normal values with a
> zero exponent. Handle that by not setting exponent to 1 when the value
> is normal for this format.
> 
> Signed-off-by: Keith Packard <keithp@keithp.com>
> ---
>   fpu/softfloat-parts.c.inc | 2 +-
>   1 file changed, 1 insertion(+), 1 deletion(-)
> 
> diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc
> index d0c43c28fb..cea854cdf1 100644
> --- a/fpu/softfloat-parts.c.inc
> +++ b/fpu/softfloat-parts.c.inc
> @@ -288,7 +288,7 @@ static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s,
>               p->frac_lo &= ~round_mask;
>           }
>   
> -        exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) != 0;
> +        exp = (p->frac_hi & DECOMPOSED_IMPLICIT_BIT) != 0 && !fmt->m68k_denormal;

That does look like a correct change.  I'll fold it in.
Please let us know if you encounter anything else.


r~

> That does look like a correct change.  I'll fold it in.
> Please let us know if you encounter anything else.

Thanks so much. With these fixes, all of my long double math library
tests in picolibc are passing now (once I fixed a bunch of additional
m68k-denorm related math library bugs). That includes the fmal tests I
mentioned, which check all four rounding modes across a fairly extensive
set of inputs.

On 21/8/23 02:32, Richard Henderson wrote:
> Motorola treats denormals with explicit integer bit set as
> having unbiased exponent 0, unlike Intel which treats it as
> having unbiased exponent 1 (like all other IEEE formats).
> 
> Add a flag on FloatFmt to differentiate the behaviour.
> 
> Reported-by: Keith Packard <keithp@keithp.com>
> Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
> ---
>   fpu/softfloat.c                |  9 +++++-
>   tests/tcg/m68k/denormal.c      | 53 ++++++++++++++++++++++++++++++++++
>   fpu/softfloat-parts.c.inc      |  5 ++--
>   tests/tcg/m68k/Makefile.target |  2 +-
>   4 files changed, 65 insertions(+), 4 deletions(-)
>   create mode 100644 tests/tcg/m68k/denormal.c

Reviewed-by: Philippe Mathieu-Daudé <philmd@linaro.org>

On Sun, 20 Aug 2023, Richard Henderson wrote:
> Motorola treats denormals with explicit integer bit set as
> having unbiased exponent 0, unlike Intel which treats it as
> having unbiased exponent 1 (like all other IEEE formats).
>
> Add a flag on FloatFmt to differentiate the behaviour.
>
> Reported-by: Keith Packard <keithp@keithp.com>
> Signed-off-by: Richard Henderson <richard.henderson@linaro.org>
> ---
> fpu/softfloat.c                |  9 +++++-
> tests/tcg/m68k/denormal.c      | 53 ++++++++++++++++++++++++++++++++++
> fpu/softfloat-parts.c.inc      |  5 ++--
> tests/tcg/m68k/Makefile.target |  2 +-
> 4 files changed, 65 insertions(+), 4 deletions(-)
> create mode 100644 tests/tcg/m68k/denormal.c
>
> diff --git a/fpu/softfloat.c b/fpu/softfloat.c
> index 0cc130ae9b..3adfa8cee0 100644
> --- a/fpu/softfloat.c
> +++ b/fpu/softfloat.c
> @@ -517,6 +517,7 @@ typedef struct {
>  *   round_mask: bits below lsb which must be rounded
>  * The following optional modifiers are available:
>  *   arm_althp: handle ARM Alternative Half Precision
> + *   m68k_denormal: explicit integer bit for extended precision may be 1

Is there a better name for this such as unbiased_exp or something telling 
it's an integer and maybe also somehow avoiding the ! when using it? This 
naming seems a bit unspecific and confusing.

Regards,
BALATON Zoltan

>  */
> typedef struct {
>     int exp_size;
> @@ -526,6 +527,7 @@ typedef struct {
>     int frac_size;
>     int frac_shift;
>     bool arm_althp;
> +    bool m68k_denormal;
>     uint64_t round_mask;
> } FloatFmt;
>
> @@ -576,7 +578,12 @@ static const FloatFmt float128_params = {
> static const FloatFmt floatx80_params[3] = {
>     [floatx80_precision_s] = { FLOATX80_PARAMS(23) },
>     [floatx80_precision_d] = { FLOATX80_PARAMS(52) },
> -    [floatx80_precision_x] = { FLOATX80_PARAMS(64) },
> +    [floatx80_precision_x] = {
> +        FLOATX80_PARAMS(64),
> +#ifdef TARGET_M68K
> +        .m68k_denormal = true,
> +#endif
> +    },
> };
>
> /* Unpack a float to parts, but do not canonicalize.  */
> diff --git a/tests/tcg/m68k/denormal.c b/tests/tcg/m68k/denormal.c
> new file mode 100644
> index 0000000000..599dafa663
> --- /dev/null
> +++ b/tests/tcg/m68k/denormal.c
> @@ -0,0 +1,53 @@
> +/*
> + * Test m68k extended double denormals.
> + */
> +
> +#include <stdio.h>
> +#include <stdint.h>
> +
> +#define X0      0x1p+16383l
> +#define Y0      0x1p-16446l
> +#define X1      0x1.1p-8223l
> +#define Y1      0x1.1p-8224l
> +
> +static volatile long double test[2][3] = {
> +    { X0, Y0, X0 * Y0 },
> +    { X1, Y1, X1 * Y1 },
> +};
> +
> +static void dump_ld(const char *label, long double ld)
> +{
> +    union {
> +        long double     d;
> +        struct {
> +            uint32_t    exp:16;
> +            uint32_t    space:16;
> +            uint32_t    h;
> +            uint32_t    l;
> +        };
> +    } u;
> +
> +    u.d = ld;
> +    printf("%12s: % -27La 0x%04x 0x%08x 0x%08x\n", label, u.d, u.exp, u.h, u.l);
> +}
> +
> +int main(void)
> +{
> +    int i, err = 0;
> +
> +    for (i = 0; i < 2; ++i) {
> +        long double x = test[i][0];
> +        long double y = test[i][1];
> +        long double build_mul = test[i][2];
> +        long double runtime_mul = x * y;
> +
> +        if (runtime_mul != build_mul) {
> +            dump_ld("x", x);
> +            dump_ld("y", y);
> +            dump_ld("build_mul", build_mul);
> +            dump_ld("runtime_mul", runtime_mul);
> +            err = 1;
> +        }
> +    }
> +    return err;
> +}
> diff --git a/fpu/softfloat-parts.c.inc b/fpu/softfloat-parts.c.inc
> index 527e15e6ab..d0c43c28fb 100644
> --- a/fpu/softfloat-parts.c.inc
> +++ b/fpu/softfloat-parts.c.inc
> @@ -118,7 +118,8 @@ static void partsN(canonicalize)(FloatPartsN *p, float_status *status,
>         } else {
>             int shift = frac_normalize(p);
>             p->cls = float_class_normal;
> -            p->exp = fmt->frac_shift - fmt->exp_bias - shift + 1;
> +            p->exp = fmt->frac_shift - fmt->exp_bias
> +                   - shift + !fmt->m68k_denormal;
>         }
>     } else if (likely(p->exp < fmt->exp_max) || fmt->arm_althp) {
>         p->cls = float_class_normal;
> @@ -256,7 +257,7 @@ static void partsN(uncanon_normal)(FloatPartsN *p, float_status *s,
>             is_tiny = !frac_addi(&discard, p, inc);
>         }
>
> -        frac_shrjam(p, 1 - exp);
> +        frac_shrjam(p, !fmt->m68k_denormal - exp);
>
>         if (p->frac_lo & round_mask) {
>             /* Need to recompute round-to-even/round-to-odd. */
> diff --git a/tests/tcg/m68k/Makefile.target b/tests/tcg/m68k/Makefile.target
> index 1163c7ef03..6ff214e60a 100644
> --- a/tests/tcg/m68k/Makefile.target
> +++ b/tests/tcg/m68k/Makefile.target
> @@ -4,7 +4,7 @@
> #
>
> VPATH += $(SRC_PATH)/tests/tcg/m68k
> -TESTS += trap
> +TESTS += trap denormal
>
> # On m68k Linux supports 4k and 8k pages (but 8k is currently broken)
> EXTRA_RUNS+=run-test-mmap-4096 # run-test-mmap-8192
>