Download raw body.
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9
On Thu, Oct 08, 2026 at 04:36:20PM -0600, Nelson H. F. Beebe wrote:
> 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
It seems this comes down to clearing a bit.
Index: lib/libm/src/ld80/s_nextafterl.c
===================================================================
RCS file: /cvs/src/lib/libm/src/ld80/s_nextafterl.c,v
diff -u -p -r1.8 s_nextafterl.c
--- lib/libm/src/ld80/s_nextafterl.c 20 Jan 2017 23:19:07 -0000 1.8
+++ lib/libm/src/ld80/s_nextafterl.c 9 Oct 2026 15:10:16 -0000
@@ -79,6 +79,7 @@ nextafterl(long double x, long double y)
if(esy==0x7fff) return x+x; /* overflow */
if(esy==0) {
volatile long double u = x*x; /* underflow */
+ hx &= ~0x80000000;
if(u==x) {
SET_LDOUBLE_WORDS(x,esx,hx,lx);
return x;
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9