port-i386/58434: single-float functions sometimes return surprisingly much precision

[campbell+netbsd%mumble.net@localhost]

>Number:         58434
>Category:       port-i386
>Synopsis:       single-float functions sometimes return surprisingly much precision
>Confidential:   no
>Severity:       serious
>Priority:       medium
>Responsible:    port-i386-maintainer
>State:          open
>Class:          sw-bug
>Submitter-Id:   net
>Arrival-Date:   Wed Jul 17 12:40:00 +0000 2024
>Originator:     Taylor R Campbell
>Release:        current, 10, 9, ...
>Organization:
The NetBSDouble Floatation
>Environment:
>Description:
https://releng.netbsd.org/b5reports/i386/2024/2024.07.16.21.10.16/test.html#lib_libm_t_log_log1p_exact

*** Check failed: /tmp/build/2024.07.16.21.10.16-i386/src/tests/lib/libm/t_log.c:384: log1pf(1) != logf(2): [3] log1pf(0x1p+0=1)=0x1.62e42fefa39efp-1=0.69314718246459961, expected 0x1.62e43p-1=0.69314718246459961

Somehow, the single-float function log1pf is returning single precision sometimes, and double precision other times -- for the same argument!
>How-To-Repeat:
$ cat foo.c
#include <math.h>
#include <stdio.h>

int
main(void)
{

	printf("log1pf(1)=%a=%a\n", log1pf(1), log1pf(1));
	printf("log1p(1)=%a=%a\n", log1p(1), log1p(1));
	fflush(stdout);
	return ferror(stdout);
}
$ make foo LDLIBS=-lm DBG=-m32\ -g\ -O2\ -fno-builtin
cc -m32 -g -O2 -fno-builtin   -o foo foo.c -lm
$ ./foo
log1pf(1)=0x1.62e42fefa39efp-1=0x1.62e43p-1
log1p(1)=0x1.62e42fefa39efp-1=0x1.62e42fefa39efp-1

Somehow, the first log1pf(1) argument to printf was evaluated in double precision, but the second one was evaluated in single precision!  And it's not because of compiler shenanigans with constant-folding log1pf or anything, because we used -fno-builtin.

What happened is roughly that log1pf(x) computes y ~ log(1 + x) in _extended precision_ internally, and leaves the result in extended precision in %st(0), the top of the x87 floating-point stack.  It is up to the caller to round the result to single precision.

The caller does that -- rounds the result to single precision -- when it is about to call log1pf again, and has to save %st(0) to memory.  But the second time around, it rounds the result to double precision to pass it to the varargs function printf, as is the rule in the ABI.  And gcc currently chooses to evaluate the arguments right-to-left, so the rightmost argument gets rounded to single precision but the leftmost argument only gets rounded to double precision.
>Fix:
just use amd64, probably, and leave i386 -- and single-precision evaluation on the x87 unit -- in the dustbin of historical textbook counterexamples

We could do something inside log1pf to force rounding to single precision.  On the other hand, leaving this bug in place will have the side effect that programs using single precision will sometimes have less error in their answers!