80-bit subnormals printed incorrectly on Debian 11 M68K

[Nelson H. F. Beebe]

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: [phone] -
- University of Utah FAX: [phone] -
- Department of Mathematics, 110 LCB Internet e-mail: ***@math.utah.edu -
- 155 S 1400 E RM 233 ***@acm.org ***@computer.org -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------

[Andreas Schwab]

Post by Nelson H. F. Beebe
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.

This is correct. In the m68881 extended float format, denormals have an
exponent bias of 0x3fff, whereas in the i387 extended float format, the
bias is 0x3fffe. That means that a normalized number in m68881 ext
format can have a biased exponent of zero.

Andreas.

[Nelson H. F. Beebe]

Andreas Schwab responds to my problem report posted at

https://lists.debian.org/debian-68k/2021/07/msg00000.html

Post by Andreas Schwab
...

Post by Nelson H. F. Beebe
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.

This is correct. In the m68881 extended float format, denormals have an
exponent bias of 0x3fff, whereas in the i387 extended float format, the
bias is 0x3fffe. That means that a normalized number in m68881 ext
format can have a biased exponent of zero.
...

I must disagree.

From the cited M68000 Family Programmer's Reference Manual, on page

1-23, in Table 1-6. Extended-Precision Real Format Summary, the biased
exponent is defined as 15 bits wide, with a bias of +16383 (0x3fff).

In IA-32 Intel Architecture Software Developer’s Manual Volume 1:
Basic Architecture, 2001, page 4-6, below Table 4-3, Floating-Point

Post by Andreas Schwab
...
The biasing constant is 127 for the single-precision format, 1023 for
the double-precision format, and 16,383 for the double
extended-precision format.
...

The Intel IA-64 Application Developer's Architecture Guide, May 1999,
(Order Number: 245188-001) on page 5-1, in Table 5-1 has

Post by Andreas Schwab
...
Total memory format width (bits) 32 64 80 128
Exponent bias +127 +1023 +16383 +16383
...

IEEE Std 754-1985, IEEE Standard for Floating-Point Arithmetic on page
3, says that the exponent bias for the single and double extended
formats is unspecified. The later 2008 and 2019 IEEE Standards do not
mention the 80-bit format.

The IEEE 754 designers wrote a series of papers in the IEEE journal
Computer in 1980, one of which is

Jerome T. Coonen
An Implementation Guide to a Proposed Standard for Floating-Point Arithmetic
Computer 13(1) 68--79 (January 1980)
https://doi.org/10.1109/MC.1980.1653344

Its Table 2 on page 70 shows a minimum exponent of the 80-bit format
as less than or equal to -16383 (at the time, CPU designs for IEEE 754
were in progress at Intel and Motorola, and possibly also National
Semiconductor, but those companies had not yet released details, or
chips). The earliest hardware implementing the IEEE 754 draft was the
Intel 8087 floating-point co-processor, which hit the market in the
summer of 1980. It is described in detail in the book

John F. Palmer and Stephen P. Morse
The 8087 Primer
Wiley (1984)
ISBN 0-471-87569-4

The Motorola 68881 floating-point co-processor reached the market in
1984. It is described in

MC68881 Floating-Point Coprocessor User's Manual, second edition
(1985)

See also

https://en.wikipedia.org/wiki/IA-64
https://en.wikipedia.org/wiki/Intel_8087
https://en.wikipedia.org/wiki/Motorola_68881
https://en.wikipedia.org/wiki/NS32000

More information about floating-point designs can be found in this
extensive bibliography of about 6900 literature references:

http://www.math.utah.edu/pub/tex/bib/fparith.bib
http://www.math.utah.edu/pub/tex/bib/fparith.bib

Thus, three independent hardware designs --- M68K, IA-32, and IA-64
--- agree that the bias in the 80-bit format is 16363 (0x3fff). That
in turn means that their in-CPU-register formats for the 80-bit size
are identical, except for the question of what bit patterns are used
for NaNs, which are not relevant for my report.

The test output in my report shows that all three architectures have
the same bit patterns for numbers on either side of the transition
from normals to subnormals. The M68K subnormals in the output should
have values of 2**k, but are in fact printed as half that. The test
loop halves the x value on each iteration (an EXACT operation in
binary floating-point), so the tabulated numbers must fall by exactly
half. They do for x86 and IA-64, but not for M68K once subnormals are
reached.

-------------------------------------------------------------------------------
- Nelson H. F. Beebe Tel: [phone] -
- University of Utah FAX: [phone] -
- Department of Mathematics, 110 LCB Internet e-mail: ***@math.utah.edu -
- 155 S 1400 E RM 233 ***@acm.org ***@computer.org -
- Salt Lake City, UT 84112-0090, USA URL: http://www.math.utah.edu/~beebe/ -
-------------------------------------------------------------------------------

[Andreas Schwab]

Post by Nelson H. F. Beebe
The Intel IA-64 Application Developer's Architecture Guide, May 1999,
(Order Number: 245188-001) on page 5-1, in Table 5-1 has

Post by Andreas Schwab
...
Total memory format width (bits) 32 64 80 128
Exponent bias +127 +1023 +16383 +16383

That are the biases for normalized numbers. What does it say about
denormalized numbers? Note that the i387 format does not allow for a
biased exponent of zero when the explicit integer bit is one, unlike the
m68881 format.

Andreas.

[Vincent Lefevre]

Post by Andreas Schwab

Post by Nelson H. F. Beebe
The Intel IA-64 Application Developer's Architecture Guide, May 1999,
(Order Number: 245188-001) on page 5-1, in Table 5-1 has

...
Total memory format width (bits) 32 64 80 128
Exponent bias +127 +1023 +16383 +16383

That are the biases for normalized numbers. What does it say about
denormalized numbers? Note that the i387 format does not allow for a
biased exponent of zero when the explicit integer bit is one, unlike the
m68881 format.

As this has been discussed in the GNU MPFR mailing-list a few days ago,
I've found that https://www.nxp.com/docs/en/reference-manual/M68000PM.pdf
(MOTOROLA M68000 FAMILY – Programmer's Reference Manual) confirms the
difference with the i387 format:

Page 1-18, Section 1.6.1 "Normalized Numbers" says:

In extended precision, the mantissa's MSB, the explicit integer bit,
can only be a one (see Figure 1-13); and the exponent can be zero.

and Section 1.6.2 "Denormalized Numbers" says:

In extended precision, the mantissa's MSB, the explicit integer bit,
can only be a zero (see Figure 1-14).

So on m86k, emin is one less the Intel x87 emin, and the printf output
is correct.

[Finn Thain]

Post by Nelson H. F. Beebe
Here is the output from Debian 11 on M68k (identical with both gcc-9
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

Here's the output of your program from a Mac IIci running Debian SID
-----
$ cat /proc/cpuinfo
CPU: 68030
MMU: 68030
FPU: 68882

I wonder if that hardware should be expected to give the same result as
68040 hardware (?) Both QEMU and Aranym emulate the latter:

CPU: 68040
MMU: 68040
FPU: 68040
Clocking: 211.0MHz
BogoMips: 140.69
Calibration: 703488 loops

Clocking: 23.1MHz
BogoMips: 5.78
Calibration: 28928 loops
$ cc bug-float80.c
$ ./a.out
Addressing is big-endian
sizeof(long double) = 12
LDBL_MANT_DIG = 64
LDBL_MIN_EXP = -16382
LDBL_MIN = 0x8.0000000000000000p-16386 = 1.681051571556046753e-4932
k = -16381 x = 0xd.eadbeefcafefeed0p-16385 = 5.848974526544159967e-4932 = 0x0001_deadbeef_cafefeed
k = -16381 x = 0xd.eadbeefcafefeed0p-16385 = 5.848974526544159967e-4932 = 0x0001_deadbeef_cafefeed
k = -16376 x = 0x8.0000000000000000p-16379 = 2.151746011591739844e-4930 = 0x0007_80000000_00000000
k = -16377 x = 0x8.0000000000000000p-16380 = 1.075873005795869922e-4930 = 0x0006_80000000_00000000
k = -16378 x = 0x8.0000000000000000p-16381 = 5.379365028979349610e-4931 = 0x0005_80000000_00000000
k = -16379 x = 0x8.0000000000000000p-16382 = 2.689682514489674805e-4931 = 0x0004_80000000_00000000
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 = 0x8.0000000000000000p-16386 = 1.681051571556046753e-4932 = 0x0000_80000000_00000000
k = -16384 x = 0x4.0000000000000000p-16386 = 8.405257857780233766e-4933 = 0x0000_40000000_00000000
k = -16385 x = 0x2.0000000000000000p-16386 = 4.202628928890116883e-4933 = 0x0000_20000000_00000000
k = -16386 x = 0x1.0000000000000000p-16386 = 2.101314464445058441e-4933 = 0x0000_10000000_00000000
k = -16387 x = 0x0.8000000000000000p-16386 = 1.050657232222529221e-4933 = 0x0000_08000000_00000000
k = -16388 x = 0x0.4000000000000000p-16386 = 5.253286161112646104e-4934 = 0x0000_04000000_00000000
k = -16389 x = 0x0.2000000000000000p-16386 = 2.626643080556323052e-4934 = 0x0000_02000000_00000000

My Aranym 68040 test agrees with the Motorola 68882 result from Stan.

My QEMU 68040 test agrees with Nelson's results.

It's suprising to see a discrepancy between the two emulators, but it
seems to be real (same rootfs image).

$ qemu-system-m68k --version
QEMU emulator version 6.0.0
Copyright (c) 2003-2021 Fabrice Bellard and the QEMU Project developers

$ aranym-mmu --version
ARAnyM 1.1.0 2020/06/16 (git:3ab938a4)
Using config file: '/home/fthain/.aranym/config'
Configuration:
SDL (compiled) : 2.0.10
SDL (linked) : 2.0.14
CPU JIT compiler : disabled
FPU JIT compiler : disabled
Addressing mode : direct
Memory check : page
Full MMU : enabled
FPU : MPFR
DSP : enabled
DSP disassembler : disabled
OpenGL support : enabled
Native features : audio bootstrap xhdi ethernet hostfs cdrom(linux) scsi jpeg vdi(opengl,sw) exec config clipboard

...

Post by Nelson H. F. Beebe
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.

NetBSD runs on m68k systems; see http://www.netbsd.org. You could also
try an earlier version of Debian (3.0 or 4.0) on m68k. And you might
want to compare the musl libc to glibc; see
https://wiki.musl-libc.org/functional-differences-from-glibc.html

Another possibility would be A/UX, if more results are needed. There is an
emulator that's intented to run A/UX called Shoebill.
https://github.com/pruten/shoebill

I don't know whether any of the many available emulators is capable of
running a NetBSD/m68k port. That might be a question for a different
mailing list.

[Brad Boyer]

Post by Finn Thain

$ cat /proc/cpuinfo
CPU: 68030
MMU: 68030
FPU: 68882

I wonder if that hardware should be expected to give the same result as
CPU: 68040
MMU: 68040
FPU: 68040

The m68k PRM does document some minor differences between the 68881/68882
and the built-in FPU in the 68040 (other than the obvious unimplemented
instructions in the 68040), but I don't think any of it would rise to
this level. They're almost entirely compatible. My first guess would be
an emulation bug. This is the sort of thing that would likely be easy to
get wrong.

My apologies for not having any of my 68040 systems available for a test
on the real hardware. I'm not even sure if any of them still work.

Brad Boyer
***@allandria.com

[Finn Thain]

Post by Brad Boyer

Post by Finn Thain

$ cat /proc/cpuinfo
CPU: 68030
MMU: 68030
FPU: 68882

I wonder if that hardware should be expected to give the same result as
CPU: 68040
MMU: 68040
FPU: 68040

The m68k PRM does document some minor differences between the 68881/68882
and the built-in FPU in the 68040 (other than the obvious unimplemented
instructions in the 68040), but I don't think any of it would rise to
this level. They're almost entirely compatible. My first guess would be
an emulation bug. This is the sort of thing that would likely be easy to
get wrong.
My apologies for not having any of my 68040 systems available for a test
on the real hardware. I'm not even sure if any of them still work.
Brad Boyer

1) 68040, Centris 650, Debian SID, gcc 9.2.1

It appears that your Motorola 68040 result agrees with your Motorola 68882
result, as Brad predicted.

2) 68040, Centris 650, NetBSD 9.1, gcc 7.5.0
3) 68030, Mac SE/30, NetBSD 9.1, gcc 7.5.0

The NetBSD test results are in agreement, but they differ from Linux. I
wonder why?

The bug-float80.c program doesn't compile in its current form in A/UX;
not only does stdint.h not exist there, but both Aople's C compiler and
an early gcc (2.7.2) repported syntax errors.

The program can probably be ported to System V Release 2 without too much
pain. You'll have to drop stdint.h. You may need to include limits.h. And
you may need to build it with "gcc -D__m68k__ bug-float80.c".

Debian/m68k 3 "woody" has gcc 2.95.4, and it fails like this:

sh-2.05a# cc -D__m68k__ bug-float80.c
bug-float80.c: In function `main':
bug-float80.c:45: hexadecimal floating constant has no exponent
bug-float80.c:45: missing white space after number `0x0.deadbee'
bug-float80.c:45: parse error before `cafefeedp'

So I think you'll want to start with a patch like this:

--- a/bug-float80.c 2021-07-22 19:08:30.000000000 +1000
+++ b/bug-float80.c 2021-07-24 23:40:27.000000000 +1000
@@ -13,9 +13,10 @@
***********************************************************************/

#include <float.h>
-#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
+#include <limits.h>

#if defined(__m68k__)
typedef long double __float80;
@@ -28,7 +29,7 @@
{
__float80 x;
int k, big_endian;
- union { long double v; uint16_t i[6]; uint32_t j[3]; } u;
+ union { long double v; unsigned short i[6]; unsigned long j[3]; } u;

u.v = 0x1p0;
big_endian = (u.i[4] == 0);
@@ -42,7 +43,7 @@
PRINTF("LDBL_MIN = %0.16La = %0.19Lg\n", (long double)LDBL_MIN, (long double)LDBL_MIN);
PRINTF("\n");

- x = 0x0.deadbeefcafefeedp-16381L;
+ memcpy(&x, "\x00\x01\x00\x00\xDE\xAD\xBE\xEF\xCA\xFE\xFE\xED", sizeof(x));
u.v = x;
k = -16381;


BTW, limits.h in A/UX 3.0.1 has this:

#if !defined(DBL_MIN)
/* Minimum normalised double */
#ifndef __STDC__
# define DBL_MIN (2.2250738585072018e-308)
#else
# define DBL_MIN (2.2250738585072014e-308)
#endif
#endif

However, the LDBL_* definitions seem to be independent of other macros.

[Antonio Vargas Gonzalez]

I could run this on my physical A1200 / 68060 compiling via VBCC for
AmigaOS if you are interested in further info.

( FWIW: I'd really expect that 68882 68040 and 68060 will agree on all
encondings, and also could be interesting to double check against WinUAE or
FS-UAE)

Post by Brad Boyer

Post by Brad Boyer

Post by Finn Thain

$ cat /proc/cpuinfo
CPU: 68030
MMU: 68030
FPU: 68882

I wonder if that hardware should be expected to give the same result

as

Post by Brad Boyer

Post by Finn Thain
CPU: 68040
MMU: 68040
FPU: 68040

The m68k PRM does document some minor differences between the

68881/68882

Post by Brad Boyer
and the built-in FPU in the 68040 (other than the obvious unimplemented
instructions in the 68040), but I don't think any of it would rise to
this level. They're almost entirely compatible. My first guess would be
an emulation bug. This is the sort of thing that would likely be easy

to

Post by Brad Boyer
get wrong.
My apologies for not having any of my 68040 systems available for a

test

Post by Brad Boyer
on the real hardware. I'm not even sure if any of them still work.
Brad Boyer

1) 68040, Centris 650, Debian SID, gcc 9.2.1

It appears that your Motorola 68040 result agrees with your Motorola 68882
result, as Brad predicted.

2) 68040, Centris 650, NetBSD 9.1, gcc 7.5.0
3) 68030, Mac SE/30, NetBSD 9.1, gcc 7.5.0

The NetBSD test results are in agreement, but they differ from Linux. I
wonder why?

The bug-float80.c program doesn't compile in its current form in A/UX;
not only does stdint.h not exist there, but both Aople's C compiler and
an early gcc (2.7.2) repported syntax errors.

The program can probably be ported to System V Release 2 without too much
pain. You'll have to drop stdint.h. You may need to include limits.h. And
you may need to build it with "gcc -D__m68k__ bug-float80.c".
sh-2.05a# cc -D__m68k__ bug-float80.c
bug-float80.c:45: hexadecimal floating constant has no exponent
bug-float80.c:45: missing white space after number `0x0.deadbee'
bug-float80.c:45: parse error before `cafefeedp'
--- a/bug-float80.c 2021-07-22 19:08:30.000000000 +1000
+++ b/bug-float80.c 2021-07-24 23:40:27.000000000 +1000
@@ -13,9 +13,10 @@
***********************************************************************/
#include <float.h>
-#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
+#include <string.h>
+#include <limits.h>
#if defined(__m68k__)
typedef long double __float80;
@@ -28,7 +29,7 @@
{
__float80 x;
int k, big_endian;
- union { long double v; uint16_t i[6]; uint32_t j[3]; } u;
+ union { long double v; unsigned short i[6]; unsigned long j[3]; } u;
u.v = 0x1p0;
big_endian = (u.i[4] == 0);
@@ -42,7 +43,7 @@
PRINTF("LDBL_MIN = %0.16La = %0.19Lg\n", (long double)LDBL_MIN,
(long double)LDBL_MIN);
PRINTF("\n");
- x = 0x0.deadbeefcafefeedp-16381L;
+ memcpy(&x, "\x00\x01\x00\x00\xDE\xAD\xBE\xEF\xCA\xFE\xFE\xED", sizeof(x));
u.v = x;
k = -16381;
#if !defined(DBL_MIN)
/* Minimum normalised double */
#ifndef __STDC__
# define DBL_MIN (2.2250738585072018e-308)
#else
# define DBL_MIN (2.2250738585072014e-308)
#endif
#endif
However, the LDBL_* definitions seem to be independent of other macros.

[Finn Thain]

Hi Nelson,

The results from your bug-float80.c program demonstrate two discrepancies:
results from NetBSD and Linux (running on Motorola processors) are
inconsistent, as are results from Aranym and QEMU (with Linux guests). Do
you know of any official bug reports about these two issues?

Regards,
Finn

[Eero Tamminen]

Hi,

Post by Finn Thain
results from NetBSD and Linux (running on Motorola processors) are
inconsistent, as are results from Aranym and QEMU (with Linux guests). Do
you know of any official bug reports about these two issues?

(Git versions) of WinUAE and Hatari emulators
could also be tried, with their software FPU
emulation options.

(In Hatari case, "--fpu-softfloat on" option.)


- Eero

[John Paul Adrian Glaubitz]

Hi Vincent!

Here, the number has been divided by 4 instead of 2. The printf output
is correct. This is a QEMU bug in the multiplication with subnormals.
https://sympa.inria.fr/sympa/arc/mpfr/2022-12/msg00036.html

I'm actually the admin of the machine QEMU machine "mitchy" that was used by the
original poster. I have actually not updated the QEMU version for some time and
I think I should actually do that :-).

I will use the sample code from the thread to see if the update made any difference.

Adrian

[Vincent Lefevre]

Post by Nelson H. F. Beebe
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
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
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

Here, the number has been divided by 4 instead of 2. The printf output
is correct. This is a QEMU bug in the multiplication with subnormals.
See the discussion in the MPFR mailing-list:
https://sympa.inria.fr/sympa/arc/mpfr/2022-12/msg00036.html