mailing list of musl libc
 help / color / mirror / Atom feed
* [musl] issue with sinhf
@ 2021-02-05  8:01 Paul Zimmermann
  2021-02-05 20:09 ` Szabolcs Nagy
  0 siblings, 1 reply; 5+ messages in thread
From: Paul Zimmermann @ 2021-02-05  8:01 UTC (permalink / raw)
  To: musl

       Hi,

a similar issue happens with sinhf:

$ 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);
}

With GNU libc:
$ gcc -fno-builtin test_sinh_musl.c -lm
$ ./a.out
x=0x1.62e4p+6 y=0x1.ffe808p+126

With musl-1.2.2:
$ gcc -fno-builtin test_sinh_musl.c $FILES
$ ./a.out
x=0x1.62e4p+6 y=-nan

Here we get NaN whereas the result is perfectly valid.

Best regards,
Paul

^ permalink raw reply	[flat|nested] 5+ messages in thread

* Re: [musl] issue with sinhf
  2021-02-05  8:01 [musl] issue with sinhf Paul Zimmermann
@ 2021-02-05 20:09 ` Szabolcs Nagy
  2021-02-05 20:12   ` Szabolcs Nagy
  2021-02-06  6:29   ` Paul Zimmermann
  0 siblings, 2 replies; 5+ messages in thread
From: Szabolcs Nagy @ 2021-02-05 20:09 UTC (permalink / raw)
  To: Paul Zimmermann; +Cc: musl

[-- Attachment #1: Type: text/plain, Size: 401 bytes --]

* 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


[-- Attachment #2: 0001-math-fix-acoshf-for-negative-inputs.patch --]
[-- Type: text/x-diff, Size: 1641 bytes --]

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


^ permalink raw reply	[flat|nested] 5+ messages in thread

* Re: [musl] issue with sinhf
  2021-02-05 20:09 ` Szabolcs Nagy
@ 2021-02-05 20:12   ` Szabolcs Nagy
  2021-02-06  6:34     ` Paul Zimmermann
  2021-02-06  6:29   ` Paul Zimmermann
  1 sibling, 1 reply; 5+ messages in thread
From: Szabolcs Nagy @ 2021-02-05 20:12 UTC (permalink / raw)
  To: Paul Zimmermann, musl

[-- Attachment #1: Type: text/plain, Size: 704 bytes --]

* Szabolcs Nagy <nsz@port70.net> [2021-02-05 21:09:34 +0100]:

> * 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

that one was for the acoshf bug, the sinhf issue was a bug
in expm1f (overflowed too early, although sinhf could have
used a less tight threshold that would have prevented this
causing an issue there)


[-- Attachment #2: 0001-math-fix-expm1f-overflow-threshold.patch --]
[-- Type: text/x-diff, Size: 1126 bytes --]

From 09f8837310082aeee019828b9aa805cf4528b3d1 Mon Sep 17 00:00:00 2001
From: Szabolcs Nagy <nsz@port70.net>
Date: Fri, 5 Feb 2021 19:51:36 +0000
Subject: [PATCH] math: fix expm1f overflow threshold

the threshold was wrong so expm1f overflowed to inf a bit too early
and on most targets uint32_t compare is faster than float compare so
use that.

this also fixes sinhf incorrectly returning nan for some values where
the internal expm1f overflowed.
---
 src/math/expm1f.c | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/math/expm1f.c b/src/math/expm1f.c
index 297e0b44..09a41afe 100644
--- a/src/math/expm1f.c
+++ b/src/math/expm1f.c
@@ -16,7 +16,6 @@
 #include "libm.h"
 
 static const float
-o_threshold = 8.8721679688e+01, /* 0x42b17180 */
 ln2_hi      = 6.9313812256e-01, /* 0x3f317180 */
 ln2_lo      = 9.0580006145e-06, /* 0x3717f7d1 */
 invln2      = 1.4426950216e+00, /* 0x3fb8aa3b */
@@ -41,7 +40,7 @@ float expm1f(float x)
 			return x;
 		if (sign)
 			return -1;
-		if (x > o_threshold) {
+		if (hx > 0x42b17217) { /* x > log(FLT_MAX) */
 			x *= 0x1p127f;
 			return x;
 		}
-- 
2.29.2


^ permalink raw reply	[flat|nested] 5+ messages in thread

* Re: [musl] issue with sinhf
  2021-02-05 20:09 ` Szabolcs Nagy
  2021-02-05 20:12   ` Szabolcs Nagy
@ 2021-02-06  6:29   ` Paul Zimmermann
  1 sibling, 0 replies; 5+ messages in thread
From: Paul Zimmermann @ 2021-02-06  6:29 UTC (permalink / raw)
  To: Szabolcs Nagy; +Cc: musl

       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



^ permalink raw reply	[flat|nested] 5+ messages in thread

* Re: [musl] issue with sinhf
  2021-02-05 20:12   ` Szabolcs Nagy
@ 2021-02-06  6:34     ` Paul Zimmermann
  0 siblings, 0 replies; 5+ messages in thread
From: Paul Zimmermann @ 2021-02-06  6:34 UTC (permalink / raw)
  To: Szabolcs Nagy; +Cc: musl

       Dear Szabolcs,

> that one was for the acoshf bug, the sinhf issue was a bug
> in expm1f (overflowed too early, although sinhf could have
> used a less tight threshold that would have prevented this
> causing an issue there)
> 
> 
> [2:text/x-diff Show Save:0001-math-fix-expm1f-overflow-threshold.patch (1kB)]

and that one fixes the sinhf issue:

GLIBC :
$ gcc -fno-builtin test_sinh_musl.c -lm
$ ./a.out
x=0x1.62e4p+6 y=0x1.ffe808p+126

musl-1.2.2 + above patch:
$ gcc -fno-builtin test_sinh_musl.c $FILES
$ ./a.out
x=0x1.62e4p+6 y=0x1.ffe808p+126

Thank you!
Paul



^ permalink raw reply	[flat|nested] 5+ messages in thread

end of thread, other threads:[~2021-02-06  6:34 UTC | newest]

Thread overview: 5+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2021-02-05  8:01 [musl] issue with sinhf Paul Zimmermann
2021-02-05 20:09 ` Szabolcs Nagy
2021-02-05 20:12   ` Szabolcs Nagy
2021-02-06  6:34     ` Paul Zimmermann
2021-02-06  6:29   ` Paul Zimmermann

mailing list of musl libc

This inbox may be cloned and mirrored by anyone:

	git clone --mirror http://inbox.vuxu.org/musl

	# If you have public-inbox 1.1+ installed, you may
	# initialize and index your mirror using the following commands:
	public-inbox-init -V1 musl musl/ http://inbox.vuxu.org/musl \
		musl@inbox.vuxu.org
	public-inbox-index musl

Example config snippet for mirrors.
Newsgroup available over NNTP:
	nntp://inbox.vuxu.org/vuxu.archive.musl


code repositories for the project(s) associated with this inbox:

	https://git.vuxu.org/mirror/musl/

AGPL code for this site: git clone https://public-inbox.org/public-inbox.git