From: "Nelson H. F. Beebe" Subject: Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9 To: tech@openbsd.org Cc: beebe@math.utah.edu Date: Thu, 8 Oct 2026 16:36:20 -0600 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 #include #include #include #include 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 - -------------------------------------------------------------------------------