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

[Paul Zimmermann <Paul.Zimmermann@inria.fr>]

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