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
signature.asc
Description: PGP signature
