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
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9
> Date: Sat, 10 Oct 2026 02:19:43 +1100
> From: Jonathan Gray <jsg@jsg.id.au>
>
> 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.
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;
>
>
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
Long-standing bug in C math library function nextafterl() on OpenBSD 5.1 to 7.9