80-bit subnormals printed incorrectly on Debian 11 M68K
I run a large farm of physical and virtual machines that we use for
software testing. We have multiple versions of most of the major
operating systems, covering the major CPU families of the past 30
years, including M68K.
In testing some numerical software on Debian 11 on M68k (emulated by
QEMU 4.2.1), I discovered that 80-bit subnormals are printed
incorrectly: they are exactly HALF their correct values.
A test program is provided below, and a snippet of its identical and
correct output on x86_64 and IA-64 (Itanium) physical hardware looks
like this around the transition from tiny normal numbers to subnormal
numbers:
k = -16380 x = 0x8.0000000000000000p-16383 = 1.344841257244837403e-4931 = 0x0003_80000000_00000000
k = -16381 x = 0x8.0000000000000000p-16384 = 6.724206286224187013e-4932 = 0x0002_80000000_00000000
k = -16382 x = 0x8.0000000000000000p-16385 = 3.362103143112093506e-4932 = 0x0001_80000000_00000000
---------- begin subnormals ----------
k = -16383 x = 0x4.0000000000000000p-16385 = 1.681051571556046753e-4932 = 0x0000_40000000_00000000
k = -16384 x = 0x2.0000000000000000p-16385 = 8.405257857780233766e-4933 = 0x0000_20000000_00000000
k = -16385 x = 0x1.0000000000000000p-16385 = 4.202628928890116883e-4933 = 0x0000_10000000_00000000
Here is the output from Debian 11 on M68k (identical with both gcc-9
and gcc-10):
k = -16380 x = 0x8.0000000000000000p-16383 = 1.344841257244837403e-4931 = 0x0003_80000000_00000000
k = -16381 x = 0x8.0000000000000000p-16384 = 6.724206286224187013e-4932 = 0x0002_80000000_00000000
k = -16382 x = 0x8.0000000000000000p-16385 = 3.362103143112093506e-4932 = 0x0001_80000000_00000000
---------- begin subnormals ----------
k = -16383 x = 0x4.0000000000000000p-16386 = 8.405257857780233766e-4933 = 0x0000_40000000_00000000
k = -16384 x = 0x2.0000000000000000p-16386 = 4.202628928890116883e-4933 = 0x0000_20000000_00000000
k = -16385 x = 0x1.0000000000000000p-16386 = 2.101314464445058441e-4933 = 0x0000_10000000_00000000
In the output, k is the power of 2. The M68K normals are correct, but
the subnormals are half their correct size in both hexadecimal and
decimal. The storage values shown in hex on the right prove that the
hardware, and QEMU, are doing the right thing, so it is definitely not
a QEMU bug.
The cause MIGHT be the incorrect value in <float.h> of LDBL_MIN_EXP:
the M68K system has -16382, whereas test output from every other
system in our farm that supports an 80-bit IEEE 754 format has -16381.
For another possibly independent check, I tried to install clang on
M68K: clang is listed in the Debian 11 package repertoire, but it
fails to install because of a missing dependency.
It is unclear to me what the appropriate bug reporting address is for
this problem: the wrong output is produced by the printf() family in
libc, but the <float.h> file is provided by the compiler. Both need
fixing.
A local patch to <float.h> won't fix the problem, because the value of
the macro LDBL_MIN_EXP has already been compiled and frozen into code
in libc.
Perhaps a seasoned debian-68k list member might be kind enough to
suggest an appropriate bug-reporting address.
At present, I have no other operating system than Debian 11 on M68K.
Web searches indicate that OpenBSD 5.1 ran on that CPU, but its
package archives have been been deleted. NetBSD 9.2 has an ISO image
for M68K, but I have not yet successfully created a VM for it.
Suggestions for other O/Ses to try are welcome.
The test code below can be run like this:
$ cc bug-float80.c && ./a.out
$ cat bug-float80.c
/***********************************************************************
** Demonstrate a bug in the Debian 11 M68K display of subnormal numbers,
** which are shown in both decimal and hexadecimal at half their correct
** values. The stored memory patterns agree on M68K, Intel x86, and
** Intel IA-64, showing that the hardware (and its QEMU emulation) is
** consistent and correct.
**
** The bug MIGHT be traceable to the off-by-one error in <float.h> of
** LDBL_MIN_EXP, which is incorrectly defined as -16382, instead of the
** correct -16381 found on the Intel systems.
**
** [20-Jul-2021]
***********************************************************************/
#include <float.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#if defined(__m68k__)
typedef long double __float80;
#endif
#define PRINTF (void)printf
int
main(void)
{
__float80 x;
int k, big_endian;
union { long double v; uint16_t i[6]; uint32_t j[3]; } u;
u.v = 0x1p0;
big_endian = (u.i[4] == 0);
PRINTF("Addressing is %s-endian\n", big_endian ? "big" : "little");
PRINTF("sizeof(long double) = %d\n", (int)sizeof(long double));
PRINTF("\n");
PRINTF("LDBL_MANT_DIG = %d\n", (int)LDBL_MANT_DIG);
PRINTF("LDBL_MIN_EXP = %d\n", (int)LDBL_MIN_EXP);
PRINTF("LDBL_MIN = %0.16La = %0.19Lg\n", (long double)LDBL_MIN, (long double)LDBL_MIN);
PRINTF("\n");
x = 0x0.deadbeefcafefeedp-16381L;
u.v = x;
k = -16381;
/*
** From the M68000 Family Programmerâ??s Reference Manual, page 1-16,
** M68000PM/A Rev. 1, 1992, memory storage of an 80-bit floating-point value
** is in a 96-bit (12-byte) big-endian-addressed field that looks like this:
**
** |-- sign |||||||||||||||--always zero
** v vvvvvvvvvvvvvvv
** 94 79 63 0
** seeeeeeeeeeeeeee0000000000000000mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm
** ^^^^^^^^^^^^^^^
** |||||||||||||||-- exponent
** ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
** ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||-- significand (formerly, mantissa)
**
** aaaaaaaaaaaaaaaa0000000000000000bbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbbcccccccccccccccccccccccccccccccc <- 32-bit fullwords
** 0 1 2 <- fullword number
**
** aaaaaaaaaaaaaaaa0000000000000000bbbbbbbbbbbbbbbbccccccccccccccccddddddddddddddddeeeeeeeeeeeeeeee <- 16-bit halfword
** 0 1 2 3 4 5 <- halfword number
**
** aaaaaaaabbbbbbbb0000000000000000ccccccccddddddddeeeeeeeeffffffffgggggggghhhhhhhhiiiiiiiijjjjjjjj <- 8-bit bytes
** 0 1 2 3 4 5 6 7 8 9 10 11 <- byte number
**
** Notice that storage bits 64--80 are unused, and set to zero.
**
** On the Intel x86 and IA-64 and the MC 68K, the binary point lies
** AFTER the first (high-order) significand digit, which is always
** explicitly stored (never hidden). The Intel formats store an
** 80-bit value right-adjusted in a 16-byte little-endian-addressed field,
** with the high order six bytes set to zero.
**
** On page 1-17, the PRM says:
**
** There is a subtle difference between the definition of an
** extended-precision number with an exponent equal to zero and a
** single- or double-precision number with an exponent equal to
** zero. The zero exponent of a single- or double-precision number
** denormalizes the number's definition, and the implied integer bit is
** zero. An extended-precision number with an exponent of zero may have
** an explicit integer bit equal to one. This results in a normalized
** number, though the exponent is equal to the minimum value. For
** simplicity, the following discussion treats all three floating-point
** formats in the same manner, where an exponent value of zero
** identifies a denormalized number. However, remember the
** extended-precision format can deviate from this rule.
**
** The code below exhibits the stored data via both 32-bit and 16-bit integer
** overlays, being careful to discard bytes in positions 64..80
*/
if (big_endian)
PRINTF("k = %6d\tx = %0.16La = %#0.19Lg = 0x%04x_%08x_%08x\n", k, x, x, u.j[0] >> 16, u.j[1], u.j[2]);
else
PRINTF("k = %6d\tx = %0.16La = %#0.19Lg = 0x%04x_%08x_%08x\n", k, x, x, u.j[2], u.j[1], u.j[0]);
if (big_endian)
PRINTF("k = %6d\tx = %0.16La = %#0.19Lg = 0x%04hx_%04hx%04hx_%04hx%04hx\n", k, x, x, u.i[0], u.i[2], u.i[3], u.i[4], u.i[5]);
else
PRINTF("k = %6d\tx = %0.16La = %#0.19Lg = 0x%04hx_%04hx%04hx_%04hx%04hx\n", k, x, x, u.i[4], u.i[3], u.i[2], u.i[1], u.i[0]);
PRINTF("\n");
x = 0x1p-16376L;
for (k = -16376; k > -16390; --k)
{
u.v = x;
if (big_endian)
PRINTF("k = %6d\tx = %0.16La = %#0.19Lg = 0x%04hx_%04hx%04hx_%04hx%04hx\n", k, x, x, u.i[0], u.i[2], u.i[3], u.i[4], u.i[5]);
else
PRINTF("k = %6d\tx = %0.16La = %#0.19Lg = 0x%04hx_%04hx%04hx_%04hx%04hx\n", k, x, x, u.i[4], u.i[3], u.i[2], u.i[1], u.i[0]);
if (k == -16382)
PRINTF("\n---------- begin subnormals ----------\n\n");
x *= 0.5L;
}
return (EXIT_SUCCESS);
}
-------------------------------------------------------------------------------
- Nelson H. F. Beebe Tel: +1 801 581 5254 -
- University of Utah FAX: +1 801 581 4148 -
- Department of Mathematics, 110 LCB Internet e-mail: beebe@math.utah.edu -
- 155 S 1400 E RM 233 beebe@acm.org beebe@computer.org -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------
Reply to: