issue with sinhf

Submitted by Szabolcs Nagy on Feb. 5, 2021, 8:09 p.m.

Details

Message ID 20210205200934.GC2447252@port70.net
State New
Series "issue with sinhf"
Headers show

Commit Message

Szabolcs Nagy Feb. 5, 2021, 8:09 p.m.
* Paul Zimmermann <Paul.Zimmermann@inria.fr> [2021-02-05 09:01:09 +0100]:
> $ cat test_sinh_musl.c
> #include <stdio.h>
> #include <stdlib.h>
> #include <math.h>
> 
> int
> main ()
> {
>   float x = 0x1.62e4p+6;
>   float y = sinhf (x);
>   printf ("x=%a y=%a\n", x, y);
> }
...
> $ gcc -fno-builtin test_sinh_musl.c $FILES
> $ ./a.out
> x=0x1.62e4p+6 y=-nan

this seems to be a bug, attaching a fix

Patch hide | download patch | download mbox

From 786b3725ff2f88a42a42c1c1fc98a34ef4c391ae Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Fri, 5 Feb 2021 18:48:19 +0000
Subject: [PATCH] math: fix acoshf for negative inputs

on some negative inputs (e.g. -0x1.1e6ae8p+5) acoshf failed to return nan.
ensure that negative inputs result nan without introducing new branches.
this was tried before in

  commit 101e6012856918440b5d7474739c3fc22a8d3b85
  math: fix acoshf on negative values

but that fix was wrong. there are 3 formulas used:

  log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1)))
  log(2*x - 1/(x+sqrt(x*x-1)))
  log(x) + 0.693147180559945309417232121458176568

the first fails on large negative inputs (may compute log1p(0) or log1p(inf)),
the second one fails on some mid range or large negative inputs (may compute
log(large) or log(inf)) and the last one fails on -0 (returns -inf).
---
 src/math/acoshf.c | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/math/acoshf.c b/src/math/acoshf.c
index 8a4ec4d5..b773d48e 100644
--- a/src/math/acoshf.c
+++ b/src/math/acoshf.c
@@ -15,12 +15,12 @@  float acoshf(float x)
 	uint32_t a = u.i & 0x7fffffff;
 
 	if (a < 0x3f800000+(1<<23))
-		/* |x| < 2, invalid if x < 1 or nan */
+		/* |x| < 2, invalid if x < 1 */
 		/* up to 2ulp error in [1,1.125] */
 		return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1)));
-	if (a < 0x3f800000+(12<<23))
-		/* |x| < 0x1p12 */
+	if (u.i < 0x3f800000+(12<<23))
+		/* 2 <= x < 0x1p12 */
 		return logf(2*x - 1/(x+sqrtf(x*x-1)));
-	/* x >= 0x1p12 */
+	/* x >= 0x1p12 or x <= -2 or nan */
 	return logf(x) + 0.693147180559945309417232121458176568f;
 }
-- 
2.29.2


Comments

Paul Zimmermann Feb. 6, 2021, 6:29 a.m.
Dear Szabolcs,

> this seems to be a bug, attaching a fix
> 
> 
> [2:text/x-diff Show Save:0001-math-fix-acoshf-for-negative-inputs.patch (2kB)]

I confirm this solves the acoshf issue:

GNU libc:
$ gcc -fno-builtin test_acosh_musl.c -lm
$ ./a.out
x=-0x1.1e6ae8p+5 y=-nan

musl-1.2.2 + above patch:
$ gcc -fno-builtin test_acosh_musl.c $FILES
$ ./a.out
x=-0x1.1e6ae8p+5 y=-nan

Thank you!
Paul