From: Jonathan Gray Subject: Re: Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9 To: "Nelson H. F. Beebe" Cc: tech@openbsd.org Date: Sat, 10 Oct 2026 02:19:43 +1100 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 > #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 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;