Download raw body.
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9
In testing some numerical software on OpenBSD 7.9, and many other
platforms in our large phyical/virtual test farm, I found that the
OpenBSD 80-bit nextafterl() function gives incorrect answers in the
subnormal region.
I verified that the bug exists in VMs running OpenBSD 5.1, 6.8, 6.9,
7.0, 7.1, 7.2, 7.3, 7.6, and 7.9.
Correct output is obtained on multiple OS versions on Linux (CentOS
and Ubuntu), DragonFlyBSD, FreeBSD, NetBSD, Solaris, ... on amd64
CPUs.
Below is a test program, and its output on OpenBSD 7.9. It asks for
the next 80-bit floating-point number below the minimum normal number,
which should be the largest subnormal number, +0x1.ffff_ffff_ffff_fffcp-16383
~= 3.362e-493.
% cat bug-nextafterl.c
#include <math.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>
int main(void)
{
long double x, y;
union {
long double f80;
uint32_t u32[4];
uint64_t u64[2];
} u;
(void)printf("sizeof(long double) = %lu\n", sizeof(long double));
(void)printf("\n");
x = +0x1.0p-16382L; /* smallest binary80 normal number */
y = nextafterl(x, 0.0L);
(void)printf("x = %.19Lg = %La\n", x, x);
(void)printf("y = nextafterl(x, 0.0L)\n");
(void)printf("y = %.19Lg = %La\n", y, y);
(void)printf("\n");
(void)printf("The ratio x/y should be > 1.0: %.18Lg = %La\n", x/y, x/y);
(void)printf("The ratio y/x should be < 1.0: %.18Lg = %La\n", y/x, y/x);
(void)printf("\n");
(void)printf("Floating-point value and memory layout for 1.0:\n");
memset(&u, '\0', sizeof(u));
u.f80 = 1.0L;
(void)printf("u.f80 = %.19Lg = %La\n", u.f80, u.f80);
(void)printf("u.u64 = 0x%016llx_%016llx\n",
(long long unsigned int)u.u64[1],
(long long unsigned int)u.u64[0]);
(void)printf("u.u32 = 0x%08x_%08x_%08x_%08x\n",
u.u32[3], u.u32[2], u.u32[1], u.u32[0]);
(void)printf("\n");
(void)printf("Floating-point value and memory layout for x:\n");
memset(&u, '\0', sizeof(u));
u.f80 = x;
(void)printf("u.f80 = %.19Lg = %La\n", u.f80, u.f80);
(void)printf("u.u64 = 0x%016llx_%016llx\n",
(long long unsigned int)u.u64[1],
(long long unsigned int)u.u64[0]);
(void)printf("u.u32 = 0x%08x_%08x_%08x_%08x\n",
u.u32[3], u.u32[2], u.u32[1], u.u32[0]);
(void)printf("\n");
(void)printf("Floating-point value and memory layout for y:\n");
memset(&u, '\0', sizeof(u));
u.f80 = y;
(void)printf("u.f80 = %.19Lg = %La\n", u.f80, u.f80);
(void)printf("u.u64 = 0x%016llx_%016llx\n",
(long long unsigned int)u.u64[1],
(long long unsigned int)u.u64[0]);
(void)printf("u.u32 = 0x%08x_%08x_%08x_%08x\n",
u.u32[3], u.u32[2], u.u32[1], u.u32[0]);
return (EXIT_SUCCESS);
}
Here is the OpenBSD output:
% cc bug-nextafterl.c -lm && ./a.out
sizeof(long double) = 16
x = 3.362103143112093506e-4932 = 0x8p-16385
y = nextafterl(x, 0.0L)
y = 6.724206286224187012e-4932 = 0xf.fffffffffffffffp-16385
The ratio x/y should be > 1.0: 0.5 = 0x8.000000000000001p-4
The ratio y/x should be < 1.0: 2 = 0xf.fffffffffffffffp-3
Floating-point value and memory layout for 1.0:
u.f80 = 1 = 0x8p-3
u.u64 = 0x0000000000003fff_8000000000000000
u.u32 = 0x00000000_00003fff_80000000_00000000
Floating-point value and memory layout for x:
u.f80 = 3.362103143112093506e-4932 = 0x8p-16385
u.u64 = 0x0000000000000001_8000000000000000
u.u32 = 0x00000000_00000001_80000000_00000000
Floating-point value and memory layout for y:
u.f80 = 6.724206286224187012e-4932 = 0xf.fffffffffffffffp-16385
u.u64 = 0x0000000000000000_ffffffffffffffff
u.u32 = 0x00000000_00000000_ffffffff_ffffffff
Notice that the reported ratios are completely wrong, and the number
below the minimum normal number is reported as TWICE that value,
instead of somewhat less than it.
For comparison, here is the correct output, computed on a Linux amd64
system:
sizeof(long double) = 16
x = 3.362103143112093506e-4932 = 0x8p-16385
y = nextafterl(x, 0.0L)
y = 3.362103143112093506e-4932 = 0x7.fffffffffffffffp-16385
The ratio x/y should be > 1.0: 1 = 0x8.000000000000001p-3
The ratio y/x should be < 1.0: 1 = 0xf.ffffffffffffffep-4
Floating-point value and memory layout for 1.0:
u.f80 = 1 = 0x8p-3
u.u64 = 0x0000000000003fff_8000000000000000
u.u32 = 0x00000000_00003fff_80000000_00000000
Floating-point value and memory layout for x:
u.f80 = 3.362103143112093506e-4932 = 0x8p-16385
u.u64 = 0x0000000000000001_8000000000000000
u.u32 = 0x00000000_00000001_80000000_00000000
Floating-point value and memory layout for y:
u.f80 = 3.362103143112093506e-4932 = 0x7.fffffffffffffffp-16385
u.u64 = 0x0000000000000000_7fffffffffffffff
u.u32 = 0x00000000_00000000_7fffffff_ffffffff
When OpenBSD developers repair the bug on amd64, the new nextafter*()
functions should be checked on other platforms, including those with
32-bit, 64-bit, 80-bit, and 128-bit arithmetic. Here is output from a
little-endian Debian 13 RISC-V64 machine, which has 128-bit long
double in software:
sizeof(long double) = 16
x = 3.362103143112093506e-4932 = 0x1p-16382
y = nextafterl(x, 0.0L)
y = 3.362103143112093506e-4932 = 0x0.ffffffffffffffffffffffffffffp-16382
The ratio x/y should be > 1.0: 1 = 0x1.0000000000000000000000000001p+0
The ratio y/x should be < 1.0: 1 = 0x1.fffffffffffffffffffffffffffep-1
Floating-point value and memory layout for 1.0:
u.f80 = 1 = 0x1p+0
u.u64 = 0x3fff000000000000_0000000000000000
u.u32 = 0x3fff0000_00000000_00000000_00000000
Floating-point value and memory layout for x:
u.f80 = 3.362103143112093506e-4932 = 0x1p-16382
u.u64 = 0x0001000000000000_0000000000000000
u.u32 = 0x00010000_00000000_00000000_00000000
Floating-point value and memory layout for y:
u.f80 = 3.362103143112093506e-4932 = 0x0.ffffffffffffffffffffffffffffp-16382
u.u64 = 0x0000ffffffffffff_ffffffffffffffff
u.u32 = 0x0000ffff_ffffffff_ffffffff_ffffffff
-------------------------------------------------------------------------------
- Nelson H. F. Beebe Tel: +1 801 581 5254 -
- University of Utah -
- 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: https://www.math.utah.edu/~beebe -
-------------------------------------------------------------------------------
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9