From: Mark Kettenis Subject: Re: Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9 To: Jonathan Gray Cc: beebe@math.utah.edu, tech@openbsd.org Date: Fri, 09 Oct 2026 20:10:26 +0200 > Date: Sat, 10 Oct 2026 02:19:43 +1100 > From: Jonathan Gray > > 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. That matches what FreeBSD and NetBSD do in their implementation. Linux doesn't do this, but has some differences further up. > 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; > >