1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
Wrong multiplication result of 'long double' on m68k
Description of problem:
In both x86 and m68k, 'long double' is an 80-bit format consisting of
- 1 bit sign, 15 bits exponent,
- 1 explicit 1 bit, 63 fraction bits.
According to <https://en.wikipedia.org/wiki/Extended_precision> and
<https://www.nxp.com/docs/en/reference-manual/M68000PRM.pdf> table 1-6 (page 1-23), with two differences:
- In m68k, there are 16 zero bits as filler after the sign/exponent
word, so that the total size is 96 bits.
- In x86, the minimum exponent of normalized numbers is 1;
in m68k, the minimum exponent of normalized numbers is 0.
The latter difference is reflected in the values of LDBL_MIN_EXP and
LDBL_MIN in gcc:
In x86:
```
$ echo '#include <float.h>' | gcc -E -dM - | grep __LDBL_MIN_EXP_
#define LDBL_MIN_EXP __LDBL_MIN_EXP__
#define __LDBL_MIN_EXP__ (-16381)
$ echo '#include <float.h>' | gcc -E -dM - | grep __LDBL_MIN__
#define __LDBL_MIN__ 3.36210314311209350626267781732175260e-4932L
#define LDBL_MIN __LDBL_MIN__
```
In m68k (I use Debian 12/Linux):
```
$ echo '#include <float.h>' | gcc -E -dM - | grep __LDBL_MIN_EXP_
#define LDBL_MIN_EXP __LDBL_MIN_EXP__
#define __LDBL_MIN_EXP__ (-16382)
$ echo '#include <float.h>' | gcc -E -dM - | grep __LDBL_MIN__
#define __LDBL_MIN__ 1.68105157155604675313e-4932L
#define LDBL_MIN __LDBL_MIN__
```
Steps to reproduce:
Take this program, foo.c:
```
/* Show extended-precision https://en.wikipedia.org/wiki/Extended_precision
multiplication bug in QEMU. */
#include <stdio.h>
static void
show (const long double *p)
{
#ifdef __m68k__
printf("<S,E: 0x%08X M: 0x%08X%08X>",
((const unsigned int *) p)[0],
((const unsigned int *) p)[1],
((const unsigned int *) p)[2]);
#else /* x86 */
printf("<S,E: 0x%04X M: 0x%08X%08X>",
((const unsigned short *) p)[4],
((const unsigned int *) p)[1],
((const unsigned int *) p)[0]);
#endif
printf (" = %La = %Lg", *p, *p);
}
static void
show_mult (long double a, long double b)
{
printf ("Factors: ");
show (&a);
printf ("\n and: ");
show (&b);
long double c = a * b;
printf ("\nProduct: ");
show (&c);
printf ("\n\n");
}
/* Return 2^n. */
static long double
pow2l (int n)
{
int k = n;
volatile long double x = 1;
volatile long double y = 2;
/* Invariant: 2^n == x * y^k. */
if (k < 0)
{
y = 0.5L;
k = - k;
}
while (k > 0)
{
if (k != 2 * (k / 2))
{
x = x * y;
k = k - 1;
}
if (k == 0)
break;
y = y * y;
k = k / 2;
}
/* Now k == 0, hence x == 2^n. */
return x;
}
int main ()
{
show_mult (pow2l (-16382), 0.5L);
show_mult (pow2l (-16381), 0.25L);
return 0;
}
```
Its output on x86:
```
$ ./a.out
Factors: <S,E: 0x0001 M: 0x8000000000000000> = 0x8p-16385 = 3.3621e-4932
and: <S,E: 0x3FFE M: 0x8000000000000000> = 0x8p-4 = 0.5
Product: <S,E: 0x0000 M: 0x4000000000000000> = 0x4p-16385 = 1.68105e-4932
Factors: <S,E: 0x0002 M: 0x8000000000000000> = 0x8p-16384 = 6.72421e-4932
and: <S,E: 0x3FFD M: 0x8000000000000000> = 0x8p-5 = 0.25
Product: <S,E: 0x0000 M: 0x4000000000000000> = 0x4p-16385 = 1.68105e-4932
```
Its output on m68k:
```
$ ./a.out
Factors: <S,E: 0x00010000 M: 0x8000000000000000> = 0x8p-16385 = 3.3621e-4932
and: <S,E: 0x3FFE0000 M: 0x8000000000000000> = 0x8p-4 = 0.5
Product: <S,E: 0x00000000 M: 0x4000000000000000> = 0x4p-16386 = 8.40526e-4933
Factors: <S,E: 0x00020000 M: 0x8000000000000000> = 0x8p-16384 = 6.72421e-4932
and: <S,E: 0x3FFD0000 M: 0x8000000000000000> = 0x8p-5 = 0.25
Product: <S,E: 0x00000000 M: 0x4000000000000000> = 0x4p-16386 = 8.40526e-4933
```
The product, computed by QEMU, is incorrect. It is only half as large as the
correct value. The expected output should be:
```
Factors: <S,E: 0x00010000 M: 0x8000000000000000> = 0x8p-16385 = 3.3621e-4932
and: <S,E: 0x3FFE0000 M: 0x8000000000000000> = 0x8p-4 = 0.5
Product: <S,E: 0x00000000 M: 0x8000000000000000> = 0x8p-16386 = 1.68105e-4932
Factors: <S,E: 0x00020000 M: 0x8000000000000000> = 0x8p-16384 = 6.72421e-4932
and: <S,E: 0x3FFD0000 M: 0x8000000000000000> = 0x8p-5 = 0.25
Product: <S,E: 0x00000000 M: 0x8000000000000000> = 0x8p-16386 = 1.68105e-4932
```
Additional information:
In QEMU's source code, I would guess that this multiplication is performed by the `floatx80_mul` function.
|