NetBSD Problem Report #58434

From www@netbsd.org  Wed Jul 17 12:38:55 2024
Return-Path: <www@netbsd.org>
Received: from mail.netbsd.org (mail.netbsd.org [199.233.217.200])
	(using TLSv1.3 with cipher TLS_AES_256_GCM_SHA384 (256/256 bits)
	 key-exchange X25519 server-signature RSA-PSS (2048 bits) server-digest SHA256
	 client-signature RSA-PSS (2048 bits) client-digest SHA256)
	(Client CN "mail.NetBSD.org", Issuer "mail.NetBSD.org CA" (not verified))
	by mollari.NetBSD.org (Postfix) with ESMTPS id AF40F1A9239
	for <gnats-bugs@gnats.NetBSD.org>; Wed, 17 Jul 2024 12:38:55 +0000 (UTC)
Message-Id: <20240717123854.453161A923A@mollari.NetBSD.org>
Date: Wed, 17 Jul 2024 12:38:54 +0000 (UTC)
From: campbell+netbsd@mumble.net
Reply-To: campbell+netbsd@mumble.net
To: gnats-bugs@NetBSD.org
Subject: single-float functions sometimes return surprisingly much precision
X-Send-Pr-Version: www-1.0

>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
>Last-Modified:  Wed Jul 17 14:55:01 +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!

>Audit-Trail:
From: Taylor R Campbell <riastradh@NetBSD.org>
To: gnats-bugs@NetBSD.org
Cc: port-i386-maintainer@NetBSD.org, gnats-admin@NetBSD.org, netbsd-bugs@NetBSD.org
Subject: Re: port-i386/58434: single-float functions sometimes return surprisingly much precision
Date: Wed, 17 Jul 2024 13:22:35 +0000

 Annotated listing of the relevant disassembly fragment for
 printf("log1pf(1)=%a=%a\n", log1pf(1), log1pf(1)), where y denotes the
 extended precision result of log(1 + x):

 	...
 	call	log1pf		# evaluate second argument first
 				# fp stack: y
 	fstps	0x28(%esp)	# pop fp stack and write (float)y to
 				# memory in single-precision
 				# fp stack: (empty)
 	...
 	call	log1pf		# evaluate first argument second
 				# fp stack: y
 	flds	0x28(%esp)	# push (float)y from memory on fp stack
 				# fp stack: (float)y y
 	fstpl	0xc(%esp)	# pop fp stack and write (double)(float)y
 				# to second printf format argument
 				# fp stack: y
 	fstpl	0x4(%esp)	# pop fp stack and write (double)y
 				# to first printf format argument
 				# fp stack: (empty)
 	movl	...,(%esp)	# set printf format string
 	call	printf

 But different compilers might do it differently!  The gcc12 shipped in
 Debian 12 saves the intermediate result with fstpl in double-precision
 rather than with fstps in single-precision.  So the output there is:

 log1pf(1)=0x1.62e42fefa39efp-1=0x1.62e42fefa39efp-1
 log1p(1)=0x1.62e42fefa39efp-1=0x1.62e42fefa39efp-1

 And of course nothing constrains the compiler's order of evaluation or
 order of explicitly rounding intermediate results, so one could also
 see

 log1pf(1)=0x1.62e43p-1=0x1.62e42fefa39efp-1
 log1p(1)=0x1.62e42fefa39efp-1=0x1.62e42fefa39efp-1

 but I haven't (yet) seen this in practice.

From: "Taylor R Campbell" <riastradh@netbsd.org>
To: gnats-bugs@gnats.NetBSD.org
Cc: 
Subject: PR/58434 CVS commit: src/tests/lib/libm
Date: Wed, 17 Jul 2024 14:52:14 +0000

 Module Name:	src
 Committed By:	riastradh
 Date:		Wed Jul 17 14:52:13 UTC 2024

 Modified Files:
 	src/tests/lib/libm: t_log.c

 Log Message:
 tests/lib/libm/t_log.c: Record xfail for PR port-i386/58434.

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


 To generate a diff of this commit:
 cvs rdiff -u -r1.18 -r1.19 src/tests/lib/libm/t_log.c

 Please note that diffs are not public domain; they are subject to the
 copyright notices on the relevant files.

Home	PR Database Search