lib/58236: nexttoward(3) is broken on subnormals

[campbell+netbsd%mumble.net@localhost]

>Number:         58236
>Category:       lib
>Synopsis:       nexttoward(3) is broken on subnormals
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    lib-bug-people
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Wed May 08 16:15:01 +0000 2024
>Originator:     Taylor R Campbell
>Release:        current, 10, 9, ...
>Organization:
The nexttowarbsd(3) Foundation
>Environment:
>Description:
nexttoward(k*n, INFINITY) should return (k + 1)*n, where n is the smallest positive subnormal and k is any integer below 2^p.  But instead it returns (k - 1)*n.

Similarly, nexttoward(k*n, -INFINITY) should return (k - 1)*n, but it returns (k + 1)*n.
>How-To-Repeat:
$ cat test.c
#include <math.h>
#include <stdio.h>

int
main(void)
{
	volatile double n = __DBL_DENORM_MIN__;
	volatile double d = INFINITY;
	printf("123*n goes up to %.0f*n\n", nextafter(123*n, d)/n);
	printf("123*n goes down to %.0f*n\n", nextafter(123*n, -d)/n);
	fflush(stdout);
	return ferror(stdout);
}
$ make test LDLIBS=-lm
cc -O2   -o test test.c -lm
$ ./test
123*n goes up to 124*n
123*n goes down to 122*n

>Fix:
Yes, please!