[PATCH v2 23/25] math: Use coshf from CORE-MATH

[DJ Delorie]

Adhemerval Zanella <adhemerval.zanella@linaro.org> writes:
> The CORE-MATH implementation is correctly rounded (for any rounding mode),
> although it should worse performance than current one.  The current
> implementation performance comes mainly from the internal usage of
> the optimize expf implementation, and shows a maximum ULPs of 2 for
> FE_TONEAREST and 3 for other rounding modes.

To what extent are we willing to sacrifice speed for accuracy?  I would
hope that gcc's -ffast-math option might work sufficiently well to
bypass our slow-but-accurate implementation on the x86 systems that are
most affected...

LGTM
Reviewed-by: DJ Delorie <dj@redhat.com>

> diff --git a/SHARED-FILES b/SHARED-FILES
> index d32c837b46..320e0b3be9 100644
> --- a/SHARED-FILES
> +++ b/SHARED-FILES
> @@ -322,3 +322,7 @@ sysdeps/ieee754/flt-32/e_atanhf.c:
>    (src/binary32/atanh/atanhf.c in CORE-MATH)
>    - The code was adapted to use glibc code style and internal
>      functions to handle errno, overflow, and underflow.
> +sysdeps/ieee754/flt-32/e_coshf.c:
> +  (src/binary32/cosh/coshf.c in CORE-MATH)
> +  - the code was adapted to use glibc code style and internal
> +    functions to handle errno, overflow, and underflow.

Ok.

> diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps
> diff --git a/sysdeps/alpha/fpu/libm-test-ulps b/sysdeps/alpha/fpu/libm-test-ulps
> diff --git a/sysdeps/arc/fpu/libm-test-ulps b/sysdeps/arc/fpu/libm-test-ulps
> diff --git a/sysdeps/arc/nofpu/libm-test-ulps b/sysdeps/arc/nofpu/libm-test-ulps
> diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps
> diff --git a/sysdeps/csky/fpu/libm-test-ulps b/sysdeps/csky/fpu/libm-test-ulps
> diff --git a/sysdeps/csky/nofpu/libm-test-ulps b/sysdeps/csky/nofpu/libm-test-ulps
> diff --git a/sysdeps/hppa/fpu/libm-test-ulps b/sysdeps/hppa/fpu/libm-test-ulps
> diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
> diff --git a/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps b/sysdeps/i386/i686/fpu/multiarch/libm-test-ulps

Ok.

> diff --git a/sysdeps/ieee754/flt-32/e_coshf.c b/sysdeps/ieee754/flt-32/e_coshf.c
> index 052d387e42..81a8ac4b99 100644
> --- a/sysdeps/ieee754/flt-32/e_coshf.c
> +++ b/sysdeps/ieee754/flt-32/e_coshf.c
> @@ -1,63 +1,117 @@
> -/* e_coshf.c -- float version of e_cosh.c.
> - */
> +/* Correctly-rounded hyperbolic cosine function for binary32 value.
>  
> -/*
> - * ====================================================
> - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
> - *
> - * Developed at SunPro, a Sun Microsystems, Inc. business.
> - * Permission to use, copy, modify, and distribute this
> - * software is freely granted, provided that this notice
> - * is preserved.
> - * ====================================================
> - */
> +Copyright (c) 2022-2024 Alexei Sibidanov.
>  
> +The original version of this file was copied from the CORE-MATH
> +project (file src/binary32/cosh/coshf.c, revision 572ecec).
> +
> +Permission is hereby granted, free of charge, to any person obtaining a copy
> +of this software and associated documentation files (the "Software"), to deal
> +in the Software without restriction, including without limitation the rights
> +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> +copies of the Software, and to permit persons to whom the Software is
> +furnished to do so, subject to the following conditions:
> +
> +The above copyright notice and this permission notice shall be included in all
> +copies or substantial portions of the Software.
> +
> +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> +SOFTWARE.
> +*/

Ok.

> +#include <stdint.h>
>  #include <math.h>
> -#include <math-narrow-eval.h>
> -#include <math_private.h>
>  #include <libm-alias-finite.h>
> -
> -static const float huge = 1.0e30;
> -static const float one = 1.0, half=0.5;
> +#include "math_config.h"
>  

Ok.

>  float
>  __ieee754_coshf (float x)
>  {
> -	float t,w;
> -	int32_t ix;
> -
> -	GET_FLOAT_WORD(ix,x);
> -	ix &= 0x7fffffff;
> -
> -    /* |x| in [0,22] */
> -	if (ix < 0x41b00000) {
> -	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
> -		if(ix<0x3eb17218) {
> -		    if (ix<0x24000000) return one;	/* cosh(tiny) = 1 */
> -		    t = __expm1f(fabsf(x));
> -		    w = one+t;
> -		    return one+(t*t)/(w+w);
> -		}
> -
> -	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
> -		t = __ieee754_expf(fabsf(x));
> -		return half*t+half/t;

Ok.

> +  static const double c[] =
> +    {
> +      1, 0x1.62e42fef4c4e7p-6, 0x1.ebfd1b232f475p-13, 0x1.c6b19384ecd93p-20
> +    };
> +  static const double ch[] =
> +    {
> +      1,                     0x1.62e42fefa39efp-6,  0x1.ebfbdff82c58fp-13,
> +      0x1.c6b08d702e0edp-20, 0x1.3b2ab6fb92e5ep-27, 0x1.5d886e6d54203p-35,
> +      0x1.430976b8ce6efp-43
> +    };
> +  static const uint64_t tb[] =
> +    {
> +      0x3fe0000000000000, 0x3fe059b0d3158574, 0x3fe0b5586cf9890f,
> +      0x3fe11301d0125b51, 0x3fe172b83c7d517b, 0x3fe1d4873168b9aa,
> +      0x3fe2387a6e756238, 0x3fe29e9df51fdee1, 0x3fe306fe0a31b715,
> +      0x3fe371a7373aa9cb, 0x3fe3dea64c123422, 0x3fe44e086061892d,
> +      0x3fe4bfdad5362a27, 0x3fe5342b569d4f82, 0x3fe5ab07dd485429,
> +      0x3fe6247eb03a5585, 0x3fe6a09e667f3bcd, 0x3fe71f75e8ec5f74,
> +      0x3fe7a11473eb0187, 0x3fe82589994cce13, 0x3fe8ace5422aa0db,
> +      0x3fe93737b0cdc5e5, 0x3fe9c49182a3f090, 0x3fea5503b23e255d,
> +      0x3feae89f995ad3ad, 0x3feb7f76f2fb5e47, 0x3fec199bdd85529c,
> +      0x3fecb720dcef9069, 0x3fed5818dcfba487, 0x3fedfc97337b9b5f,
> +      0x3feea4afa2a490da, 0x3fef50765b6e4540
> +    };
> +  const double iln2 = 0x1.71547652b82fep+5;

Ok.

> +  double z = x;
> +  uint32_t ax = asuint (x) << 1;
> +  if (__glibc_unlikely (ax > 0x8565a9f8u))
> +    { /* |x| >~ 89.4 */
> +      if (ax >= 0xff000000u)
> +	{
> +	  if (ax << 8)
> +	    return x + x;  /* nan */
> +	  return INFINITY; /* +-inf */
>  	}

Ok.

> -
> -    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
> -	if (ix < 0x42b17180)  return half*__ieee754_expf(fabsf(x));
> -
> -    /* |x| in [log(maxdouble), overflowthresold] */
> -	if (ix<=0x42b2d4fc) {
> -	    w = __ieee754_expf(half*fabsf(x));
> -	    t = half*w;
> -	    return t*w;

Ok.

> +      return __math_oflowf (0);
> +    }
> +  if (__glibc_unlikely (ax < 0x7c000000u))
> +    { /* |x| < 0.125 */
> +      if (__glibc_unlikely (ax < 0x74000000u))
> +	{					   /* |x| < 0x1p-11 */
> +	  if (__glibc_unlikely (ax < 0x66000000u)) /* |x| < 0x1p-24 */
> +	    return fmaf (fabsf (x), 0x1p-25, 1.0f);
> +	  return (0.5f * x) * x + 1.0f;
>  	}

Ok.

> -
> -    /* x is INF or NaN */
> -	if(ix>=0x7f800000) return x*x;
> -
> -    /* |x| > overflowthresold, cosh(x) overflow */
> -	return math_narrow_eval (huge*huge);

Ok.

> +      static const double cp[] =
> +	{
> +	  0x1.fffffffffffe3p-2,  0x1.55555555723cfp-5,
> +	  0x1.6c16bee4a5986p-10, 0x1.a0483fc0328f7p-16
> +	};
> +      double z2 = z * z;
> +      double z4 = z2 * z2;
> +      return 1 + z2 * ((cp[0] + z2 * cp[1]) + z4 * (cp[2] + z2 * (cp[3])));
> +    }

Ok.  The "1" instead of "1.0f" is inconsistent, but doesn't matter...

> +  double a = iln2 * z;
> +  double ia = roundeven_finite (a);
> +  double h = a - ia;
> +  double h2 = h * h;
> +  int64_t jp = asuint64 (ia + 0x1.8p52);
> +  int64_t jm = -jp;
> +  double sp = asdouble (tb[jp & 31] + ((jp >> 5) << 52));
> +  double sm = asdouble (tb[jm & 31] + ((jm >> 5) << 52));
> +  double te = c[0] + h2 * c[2];
> +  double to = (c[1] + h2 * c[3]);
> +  double rp = sp * (te + h * to);
> +  double rm = sm * (te - h * to);
> +  double r = rp + rm;
> +  float ub = r;
> +  double lb = r - 1.45e-10 * r;

Ok.

> +  if (__glibc_unlikely (ub != lb))
> +    {
> +      const double iln2h = 0x1.7154765p+5;
> +      const double iln2l = 0x1.5c17f0bbbe88p-26;
> +      h = (iln2h * z - ia) + iln2l * z;
> +      h2 = h * h;
> +      te = ch[0] + h2 * ch[2] + (h2 * h2) * (ch[4] + h2 * ch[6]);
> +      to = ch[1] + h2 * (ch[3] + h2 * ch[5]);
> +      r = sp * (te + h * to) + sm * (te - h * to);
> +      ub = r;
> +    }
> +  return ub;
>  }
>  libm_alias_finite (__ieee754_coshf, __coshf)

Ok.

> diff --git a/sysdeps/loongarch/lp64/libm-test-ulps b/sysdeps/loongarch/lp64/libm-test-ulps
> diff --git a/sysdeps/microblaze/libm-test-ulps b/sysdeps/microblaze/libm-test-ulps
> diff --git a/sysdeps/mips/mips32/libm-test-ulps b/sysdeps/mips/mips32/libm-test-ulps
> diff --git a/sysdeps/mips/mips64/libm-test-ulps b/sysdeps/mips/mips64/libm-test-ulps
> diff --git a/sysdeps/or1k/fpu/libm-test-ulps b/sysdeps/or1k/fpu/libm-test-ulps
> diff --git a/sysdeps/or1k/nofpu/libm-test-ulps b/sysdeps/or1k/nofpu/libm-test-ulps
> diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
> diff --git a/sysdeps/powerpc/nofpu/libm-test-ulps b/sysdeps/powerpc/nofpu/libm-test-ulps
> diff --git a/sysdeps/riscv/nofpu/libm-test-ulps b/sysdeps/riscv/nofpu/libm-test-ulps
> diff --git a/sysdeps/riscv/rvd/libm-test-ulps b/sysdeps/riscv/rvd/libm-test-ulps
> diff --git a/sysdeps/s390/fpu/libm-test-ulps b/sysdeps/s390/fpu/libm-test-ulps
> diff --git a/sysdeps/sh/libm-test-ulps b/sysdeps/sh/libm-test-ulps
> diff --git a/sysdeps/sparc/fpu/libm-test-ulps b/sysdeps/sparc/fpu/libm-test-ulps
> diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps

Ok.