[musl] musl 1.2.4

[musl] musl 1.2.4

[Paul Zimmermann]

       Hi,

while updating our comparison of mathematical functions with Musl 1.2.4,
I noticed the following two issues with acoshl and exp10l in double extended
precision:

zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl acosh 1000
Checking acosh with musl-1.2.4
Using seed 2804715
Using 6 threads
NEW acosh 0 -1 -0x6.e2368c0ed74e5698p+16l [-nan] [inf] inf inf
libm gives -0x4.b4d6a621e8e631f8p+0l
mpfr gives nanl

zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl exp10 1000
Checking exp10 with musl-1.2.4
Using seed 2807610
Using 6 threads
exp10 0 -1 0x2.68826a13ef3fde64p+16376l [-nan] [inf] inf inf
libm gives nanl
mpfr gives infl

These issues were already present in 1.2.2 at least, but only detected
recently by our program.

Best regards,
Paul

Re: [musl] musl 1.2.4

[Szabolcs Nagy]

* Paul Zimmermann <Paul.Zimmermann@inria.fr> [2023-08-16 16:39:25 +0200]:
> while updating our comparison of mathematical functions with Musl 1.2.4,
> I noticed the following two issues with acoshl and exp10l in double extended
> precision:
> 
> zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl acosh 1000
> Checking acosh with musl-1.2.4
> Using seed 2804715
> Using 6 threads
> NEW acosh 0 -1 -0x6.e2368c0ed74e5698p+16l [-nan] [inf] inf inf
> libm gives -0x4.b4d6a621e8e631f8p+0l
> mpfr gives nanl


yes it seems acoshl(x) does not handle x<0
(and can go wrong for x in [-2,-0x1p32])

> 
> zimmerma@coriandre:~/svn/tbd/20/src/binary80$ VERBOSE=-v ./doit.musl exp10 1000
> Checking exp10 with musl-1.2.4
> Using seed 2807610
> Using 6 threads
> exp10 0 -1 0x2.68826a13ef3fde64p+16376l [-nan] [inf] inf inf
> libm gives nanl
> mpfr gives infl

this seems to be a bug in powl: the reducl function can overflow.


thanks for the reports.

[musl] [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0

[Szabolcs Nagy]

acosh(x) is nan for x < 1, but x < 0 cases were not handled specially
and acoshl gave wrong result for some -0x1p32 < x < -2 values, e.g.:

acoshl(-0x1p20) returned -inf,
acoshl(-0x1.4p20) returned -0x1.db365758403aa9acp+0L,

fixed by checking the sign bit and handling it specially.

reported by Paul Zimmermann.
---
 src/math/acoshl.c | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/math/acoshl.c b/src/math/acoshl.c
index 8d4b43f6..943cec17 100644
--- a/src/math/acoshl.c
+++ b/src/math/acoshl.c
@@ -10,14 +10,18 @@ long double acoshl(long double x)
 long double acoshl(long double x)
 {
 	union ldshape u = {x};
-	int e = u.i.se & 0x7fff;
+	int e = u.i.se;
 
 	if (e < 0x3fff + 1)
-		/* |x| < 2, invalid if x < 1 or nan */
+		/* 0 <= x < 2, invalid if x < 1 */
 		return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
 	if (e < 0x3fff + 32)
-		/* |x| < 0x1p32 */
+		/* 2 <= x < 0x1p32 */
 		return logl(2*x - 1/(x+sqrtl(x*x-1)));
+	if (e & 0x8000)
+		/* x < 0 or x = -0, invalid */
+		return (x - x) / (x - x);
+	/* 0x1p32 <= x or nan */
 	return logl(x) + 0.693147180559945309417232121458176568L;
 }
 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
-- 
2.41.0

[musl] Re: [PATCH 1/2] math: fix ld80 acoshl(x) for x < 0

[Paul Zimmermann]

       Dear Szabolcs,

ok for me, I now get a largest error of 2.99 ulp:

NEW acosh 0 -1 0x1.1ecdb5b8f0c5d79p+0l [3] [2.99] 2.98085 2.980840623325726

Thank you for fixing that rapidly!

Paul

> Date: Fri, 18 Aug 2023 23:16:00 +0200
> From: Szabolcs Nagy <nsz@port70.net>
> Cc: musl@lists.openwall.com
> 
> acosh(x) is nan for x < 1, but x < 0 cases were not handled specially
> and acoshl gave wrong result for some -0x1p32 < x < -2 values, e.g.:
> 
> acoshl(-0x1p20) returned -inf,
> acoshl(-0x1.4p20) returned -0x1.db365758403aa9acp+0L,
> 
> fixed by checking the sign bit and handling it specially.
> 
> reported by Paul Zimmermann.
> ---
>  src/math/acoshl.c | 10 +++++++---
>  1 file changed, 7 insertions(+), 3 deletions(-)
> 
> diff --git a/src/math/acoshl.c b/src/math/acoshl.c
> index 8d4b43f6..943cec17 100644
> --- a/src/math/acoshl.c
> +++ b/src/math/acoshl.c
> @@ -10,14 +10,18 @@ long double acoshl(long double x)
>  long double acoshl(long double x)
>  {
>  	union ldshape u = {x};
> -	int e = u.i.se & 0x7fff;
> +	int e = u.i.se;
>  
>  	if (e < 0x3fff + 1)
> -		/* |x| < 2, invalid if x < 1 or nan */
> +		/* 0 <= x < 2, invalid if x < 1 */
>  		return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1)));
>  	if (e < 0x3fff + 32)
> -		/* |x| < 0x1p32 */
> +		/* 2 <= x < 0x1p32 */
>  		return logl(2*x - 1/(x+sqrtl(x*x-1)));
> +	if (e & 0x8000)
> +		/* x < 0 or x = -0, invalid */
> +		return (x - x) / (x - x);
> +	/* 0x1p32 <= x or nan */
>  	return logl(x) + 0.693147180559945309417232121458176568L;
>  }
>  #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
> -- 
> 2.41.0
> 
> 

[musl] [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)

[Szabolcs Nagy]

powl used >= LDBL_MAX as infinity check, but LDBL_MAX is finite, so
this can cause wrong results e.g. powl(LDBL_MAX, 0.5) returned inf
or powl(2, LDBL_MAX) returned inf without raising overflow.

huge y values (close to LDBL_MAX) could cause intermediate results to
overflow (computing y * log2(x) with more than long double precision)
and e.g. powl(0.5, 0x1p16380L) or powl(10, 0x1p16380L) returned nan.
this is fixed by handling huge y early since that always overflows or
underflows.

reported by Paul Zimmermann against expl10 (which uses powl).
---
 src/math/powl.c | 34 +++++++++++++++++++++-------------
 1 file changed, 21 insertions(+), 13 deletions(-)

diff --git a/src/math/powl.c b/src/math/powl.c
index 5b6da07b..6f64ea71 100644
--- a/src/math/powl.c
+++ b/src/math/powl.c
@@ -212,25 +212,33 @@ long double powl(long double x, long double y)
 	}
 	if (x == 1.0)
 		return 1.0; /* 1**y = 1, even if y is nan */
-	if (x == -1.0 && !isfinite(y))
-		return 1.0; /* -1**inf = 1 */
 	if (y == 0.0)
 		return 1.0; /* x**0 = 1, even if x is nan */
 	if (y == 1.0)
 		return x;
-	if (y >= LDBL_MAX) {
-		if (x > 1.0 || x < -1.0)
-			return INFINITY;
-		if (x != 0.0)
-			return 0.0;
-	}
-	if (y <= -LDBL_MAX) {
-		if (x > 1.0 || x < -1.0)
+	/* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0
+	   if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows
+	   if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so
+	   x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */
+	if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) {
+		/* y is not an odd int */
+		if (x == -1.0)
+			return 1.0;
+		if (y == INFINITY) {
+			if (x > 1.0 || x < -1.0)
+				return INFINITY;
 			return 0.0;
-		if (x != 0.0 || y == -INFINITY)
+		}
+		if (y == -INFINITY) {
+			if (x > 1.0 || x < -1.0)
+				return 0.0;
 			return INFINITY;
+		}
+		if ((x > 1.0 || x < -1.0) == (y > 0))
+			return huge * huge;
+		return twom10000 * twom10000;
 	}
-	if (x >= LDBL_MAX) {
+	if (x == INFINITY) {
 		if (y > 0.0)
 			return INFINITY;
 		return 0.0;
@@ -253,7 +261,7 @@ long double powl(long double x, long double y)
 			yoddint = 1;
 	}
 
-	if (x <= -LDBL_MAX) {
+	if (x == -INFINITY) {
 		if (y > 0.0) {
 			if (yoddint)
 				return -INFINITY;
-- 
2.41.0

[musl] Re: [PATCH 2/2] math: fix ld80 powl(x,huge) and powl(LDBL_MAX,small)

[Paul Zimmermann]

       Dear Szabolcs,

looks good to me, I now get a largest error of 40.1 ulps for exp10l:

NEW exp10 0 0 0xd.4135bbaf96d9716p+8l [40] [40.1] 40.0525 40.05240557189884

Paul

> Date: Fri, 18 Aug 2023 23:17:38 +0200
> From: Szabolcs Nagy <nsz@port70.net>
> Cc: musl@lists.openwall.com
> 
> powl used >= LDBL_MAX as infinity check, but LDBL_MAX is finite, so
> this can cause wrong results e.g. powl(LDBL_MAX, 0.5) returned inf
> or powl(2, LDBL_MAX) returned inf without raising overflow.
> 
> huge y values (close to LDBL_MAX) could cause intermediate results to
> overflow (computing y * log2(x) with more than long double precision)
> and e.g. powl(0.5, 0x1p16380L) or powl(10, 0x1p16380L) returned nan.
> this is fixed by handling huge y early since that always overflows or
> underflows.
> 
> reported by Paul Zimmermann against expl10 (which uses powl).
> ---
>  src/math/powl.c | 34 +++++++++++++++++++++-------------
>  1 file changed, 21 insertions(+), 13 deletions(-)
> 
> diff --git a/src/math/powl.c b/src/math/powl.c
> index 5b6da07b..6f64ea71 100644
> --- a/src/math/powl.c
> +++ b/src/math/powl.c
> @@ -212,25 +212,33 @@ long double powl(long double x, long double y)
>  	}
>  	if (x == 1.0)
>  		return 1.0; /* 1**y = 1, even if y is nan */
> -	if (x == -1.0 && !isfinite(y))
> -		return 1.0; /* -1**inf = 1 */
>  	if (y == 0.0)
>  		return 1.0; /* x**0 = 1, even if x is nan */
>  	if (y == 1.0)
>  		return x;
> -	if (y >= LDBL_MAX) {
> -		if (x > 1.0 || x < -1.0)
> -			return INFINITY;
> -		if (x != 0.0)
> -			return 0.0;
> -	}
> -	if (y <= -LDBL_MAX) {
> -		if (x > 1.0 || x < -1.0)
> +	/* if y*log2(x) < log2(LDBL_TRUE_MIN)-1 then x^y uflows to 0
> +	   if y*log2(x) > -log2(LDBL_TRUE_MIN)+1 > LDBL_MAX_EXP then x^y oflows
> +	   if |x|!=1 then |log2(x)| > |log(x)| > LDBL_EPSILON/2 so
> +	   x^y oflows/uflows if |y|*LDBL_EPSILON/2 > -log2(LDBL_TRUE_MIN)+1 */
> +	if (fabsl(y) > 2*(-LDBL_MIN_EXP+LDBL_MANT_DIG+1)/LDBL_EPSILON) {
> +		/* y is not an odd int */
> +		if (x == -1.0)
> +			return 1.0;
> +		if (y == INFINITY) {
> +			if (x > 1.0 || x < -1.0)
> +				return INFINITY;
>  			return 0.0;
> -		if (x != 0.0 || y == -INFINITY)
> +		}
> +		if (y == -INFINITY) {
> +			if (x > 1.0 || x < -1.0)
> +				return 0.0;
>  			return INFINITY;
> +		}
> +		if ((x > 1.0 || x < -1.0) == (y > 0))
> +			return huge * huge;
> +		return twom10000 * twom10000;
>  	}
> -	if (x >= LDBL_MAX) {
> +	if (x == INFINITY) {
>  		if (y > 0.0)
>  			return INFINITY;
>  		return 0.0;
> @@ -253,7 +261,7 @@ long double powl(long double x, long double y)
>  			yoddint = 1;
>  	}
>  
> -	if (x <= -LDBL_MAX) {
> +	if (x == -INFINITY) {
>  		if (y > 0.0) {
>  			if (yoddint)
>  				return -INFINITY;
> -- 
> 2.41.0
> 
>