Funny results with long double denorms on m68k

* Funny results with long double denorms on m68k
@ 2023-08-20 20:55 Keith Packard via
  2023-08-21  8:34 ` Laurent Vivier
  0 siblings, 1 reply; 4+ messages in thread
From: Keith Packard via @ 2023-08-20 20:55 UTC (permalink / raw)
  To: qemu-devel; +Cc: Laurent Vivier

[-- Attachment #1: Type: text/plain, Size: 3264 bytes --]

I'm doing some testing of an fmal implementation and discovered some
"odd" results on m68k where the emulated 80-bit FPU is generating
results that don't match how GCC computes things. Assuming gcc is
correct, this means there are some subtle bugs in how qemu is handling
denorms for this platform.

From what I can discern from examining gcc output, the 68881 80-bit
format handles denorms differently than the 8087. It allows normal
numbers to have an exponent field of zero -- with an explicit 1
bit. Furthermore, denorms don't have their exponent biased by 1, the
natural 0 value is correct.

As a simple example of this issue, I've attached a sample program which
multiplies the smallest denorm (__LDBL_DENORM_MIN__, or 0x1p-16446l) by
the largest value (__LDBL_MAX__, or 0x1p+16383l). This should be
0x1p-63, but qemu computes this as 0x1p-62. In raw form:

           x:  0x1p+16383                 0x7ffe 0x80000000 0x00000000
           y:  0x1p-16446                 0x0000 0x00000000 0x00000001
   build_mul:  0x1p-63                    0x3fc0 0x80000000 0x00000000
 runtime_mul:  0x1p-62                    0x3fc1 0x80000000 0x00000000

Looking just at the exponents, we see that the runtime computed value
has an exponent of 0x3fc1. Subtracting the bias of 0x3fff, we get -62.

This particular fault comes from converting the denorm into canonical
form in fpu/softfloat-parts.c.inc:partsN(canonicalize):

        } else {
            int shift = frac_normalize(p);
            p->cls = float_class_normal;
            p->exp = fmt->frac_shift - fmt->exp_bias - shift + 1;
        }

the extra '1' added there is the exponent bias required for standard
denorm format values.

This is only one of a number of related faults; there are similar issues
when converting back from canonical form to 68881 form. As that involves
complicated rounding semantics, it's a lot more difficult to figure out
how to fix it. For instance:

           x:  0x1.1p-8223                0x1fe0 0x88000000 0x00000000
           y:  0x1.1p-8224                0x1fdf 0x88000000 0x00000000
   build_mul:  0x1p-16446                 0x0000 0x00000000 0x00000001
 runtime_mul:  0x0p+0                     0x0000 0x00000000 0x00000000

In this case, the multiplication results in a value just larger than 1/2
of the smallest denorm. That should round up to the smallest denorm, but
qemu generates zero instead.

---------

#include <stdio.h>
#include <stdint.h>

#define X       0x1p+16383l
#define Y       0x1p-16446l

static long double build_mul = X * Y;
static volatile long double x = X;
static volatile long double y = Y;

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)
{
    long double runtime_mul = x * y;

    dump_ld("x", x);
    dump_ld("y", y);
    dump_ld("build_mul", build_mul);
    dump_ld("runtime_mul", runtime_mul);
    return 0;
}

-- 
-keith

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

^ permalink raw reply	[flat|nested] 4+ messages in thread