math: Use coshf from CORE-MATH

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.

The code was adapted to glibc style and to use the definition of
math_config.h (to handle errno, overflow, and underflow).

Benchtest on x64_64 (Ryzen 9 5900X, gcc 14.2.1), aarch64 (Neoverse-N1,
gcc 13.3.1), and powerpc (POWER10, gcc 13.2.1):

Latency                      master        patched   improvement
x86_64                      40.6995        49.0737       -20.58%
x86_64v2                    40.5841        44.3604        -9.30%
x86_64v3                    39.3879        39.7502        -0.92%
i686                       112.3380       129.8570       -15.59%
aarch64 (Neoverse)          18.6914        17.0946         8.54%
power10                     11.1343        9.3245         16.25%

reciprocal-throughput        master        patched   improvement
x86_64                      18.6471        24.1077       -29.28%
x86_64v2                    17.7501        20.2946       -14.34%
x86_64v3                    17.8262        17.1877         3.58%
i686                        64.1454        86.5645       -34.95%
aarch64 (Neoverse)          9.77226        12.2314       -25.16%
power10                      4.0200        5.3316        -32.63%

Signed-off-by: Alexei Sibidanov <sibid@uvic.ca>
Signed-off-by: Paul Zimmermann <Paul.Zimmermann@inria.fr>
Signed-off-by: Adhemerval Zanella <adhemerval.zanella@linaro.org>
Reviewed-by: DJ Delorie <dj@redhat.com>

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.

sysdeps/aarch64/libm-test-ulps

@@ -774,7 +774,6 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_advsimd":
@@ -783,7 +782,6 @@ float: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_sve":
@@ -792,12 +790,10 @@ float: 2
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: "cospi":

sysdeps/alpha/fpu/libm-test-ulps

@@ -621,22 +621,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: Real part of "cpow":

sysdeps/arc/fpu/libm-test-ulps

@@ -493,19 +493,15 @@ float: 2
 
 Function: "cosh":
 double: 3
-float: 3
 
 Function: "cosh_downward":
 double: 3
-float: 1
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: "cosh_upward":
 double: 3
-float: 2
 
 Function: Real part of "cpow":
 double: 9

sysdeps/arc/nofpu/libm-test-ulps

@@ -123,7 +123,6 @@ float: 1
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: Real part of "cpow":
 double: 2

sysdeps/arm/libm-test-ulps

@@ -84,8 +84,6 @@ float: 1
 Function: "atan":
 double: 1
 
-Function: "atan2":
-
 Function: "atan2_downward":
 double: 1
 
@@ -553,19 +551,15 @@ float: 2
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: "cosh_upward":
 double: 2
-float: 2
 
 Function: "cospi":
 double: 1
@@ -750,7 +744,7 @@ float: 2
 
 Function: Real part of "ctanh_downward":
 double: 4
-float: 2
+float: 3
 
 Function: Imaginary part of "ctanh_downward":
 double: 6

sysdeps/csky/fpu/libm-test-ulps

@@ -485,19 +485,15 @@ float: 1
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: "cosh_upward":
 double: 2
-float: 2
 
 Function: Real part of "cpow":
 double: 2

sysdeps/csky/nofpu/libm-test-ulps

@@ -483,19 +483,15 @@ float: 2
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 1
-float: 1
 
 Function: "cosh_towardzero":
 double: 1
-float: 1
 
 Function: "cosh_upward":
 double: 1
-float: 2
 
 Function: Real part of "cpow":
 double: 2

sysdeps/hppa/fpu/libm-test-ulps

@@ -563,19 +563,15 @@ float: 2
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: "cosh_upward":
 double: 2
-float: 2
 
 Function: "cospi":
 double: 1

sysdeps/i386/fpu/libm-test-ulps

@@ -850,25 +850,21 @@ ldouble: 2
 
 Function: "cosh":
 double: 1
-float: 2
 float128: 2
 ldouble: 3
 
 Function: "cosh_downward":
 double: 3
-float: 1
 float128: 3
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 float128: 3
 ldouble: 3
 
 Function: "cosh_upward":
 double: 4
-float: 2
 float128: 3
 ldouble: 3
 

sysdeps/i386/i686/fpu/multiarch/libm-test-ulps

@@ -850,25 +850,21 @@ ldouble: 2
 
 Function: "cosh":
 double: 1
-float: 2
 float128: 2
 ldouble: 3
 
 Function: "cosh_downward":
 double: 3
-float: 1
 float128: 3
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 float128: 3
 ldouble: 3
 
 Function: "cosh_upward":
 double: 4
-float: 2
 float128: 3
 ldouble: 3
 

sysdeps/ieee754/flt-32/e_atan2f.c

@@ -130,9 +130,9 @@ __ieee754_atan2f (float y, float x)
       if (yinf & xinf)
 	{
 	  if (ux >> 31)
-	    return 0x1.2d97c7f3321d2p+1 * sgn[uy >> 31];
+	    return 0x1.2d97c7f3321d2p+1 * sgn[uy >> 31]; /* +/-3pi/4 */
 	  else
-	    return 0x1.921fb54442d18p-1 * sgn[uy >> 31];
+	    return 0x1.921fb54442d18p-1 * sgn[uy >> 31]; /* +/-pi/4 */
 	}
       if (xinf)
 	{
@@ -146,7 +146,7 @@ __ieee754_atan2f (float y, float x)
     }
   if (__glibc_unlikely (ay == 0))
     {
-      if (__glibc_unlikely (!(ay | ax)))
+      if (__glibc_unlikely (!ax))
 	{
 	  uint32_t i = (uy >> 31) * 4 + (ux >> 31) * 2;
 	  if (ux >> 31)

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) 2022-2024 Alexei Sibidanov.
- * ====================================================
- * 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.
- * ====================================================
- */
 
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/cosh/coshf.c, revision 5c58ea1).
+
+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.
+*/
+
+#include <stdint.h>
 #include <math.h>
-#include <math-narrow-eval.h>
-#include <math_private.h>
 #include <libm-alias-finite.h>
-
+#include "math_config.h"
-static const float huge = 1.0e30;
-static const float one = 1.0, half=0.5;
 
 float
 __ieee754_coshf (float x)
 {
-	float t,w;
+  static const double c[] =
-	int32_t ix;
+    {
-
+      1, 0x1.62e42fef4c4e7p-6, 0x1.ebfd1b232f475p-13, 0x1.c6b19384ecd93p-20
-	GET_FLOAT_WORD(ix,x);
+    };
-	ix &= 0x7fffffff;
+  static const double ch[] =
-
+    {
-    /* |x| in [0,22] */
+      1,                     0x1.62e42fefa39efp-6,  0x1.ebfbdff82c58fp-13,
-	if (ix < 0x41b00000) {
+      0x1.c6b08d702e0edp-20, 0x1.3b2ab6fb92e5ep-27, 0x1.5d886e6d54203p-35,
-	    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
+      0x1.430976b8ce6efp-43
-		if(ix<0x3eb17218) {
+    };
-		    if (ix<0x24000000) return one;	/* cosh(tiny) = 1 */
+  static const uint64_t tb[] =
-		    t = __expm1f(fabsf(x));
+    {
-		    w = one+t;
+      0x3fe0000000000000, 0x3fe059b0d3158574, 0x3fe0b5586cf9890f,
-		    return one+(t*t)/(w+w);
+      0x3fe11301d0125b51, 0x3fe172b83c7d517b, 0x3fe1d4873168b9aa,
-		}
+      0x3fe2387a6e756238, 0x3fe29e9df51fdee1, 0x3fe306fe0a31b715,
-
+      0x3fe371a7373aa9cb, 0x3fe3dea64c123422, 0x3fe44e086061892d,
-	    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
+      0x3fe4bfdad5362a27, 0x3fe5342b569d4f82, 0x3fe5ab07dd485429,
-		t = __ieee754_expf(fabsf(x));
+      0x3fe6247eb03a5585, 0x3fe6a09e667f3bcd, 0x3fe71f75e8ec5f74,
-		return half*t+half/t;
+      0x3fe7a11473eb0187, 0x3fe82589994cce13, 0x3fe8ace5422aa0db,
+      0x3fe93737b0cdc5e5, 0x3fe9c49182a3f090, 0x3fea5503b23e255d,
+      0x3feae89f995ad3ad, 0x3feb7f76f2fb5e47, 0x3fec199bdd85529c,
+      0x3fecb720dcef9069, 0x3fed5818dcfba487, 0x3fedfc97337b9b5f,
+      0x3feea4afa2a490da, 0x3fef50765b6e4540
+    };
+  const double iln2 = 0x1.71547652b82fep+5;
+  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 */
 	}
-
+      return __math_oflowf (0);
-    /* |x| in [22, log(maxdouble)] return half*exp(|x|) */
+    }
-	if (ix < 0x42b17180)  return half*__ieee754_expf(fabsf(x));
+  if (__glibc_unlikely (ax < 0x7c000000u))
-
+    { /* |x| < 0.125 */
-    /* |x| in [log(maxdouble), overflowthresold] */
+      if (__glibc_unlikely (ax < 0x74000000u))
-	if (ix<=0x42b2d4fc) {
+	{					   /* |x| < 0x1p-11 */
-	    w = __ieee754_expf(half*fabsf(x));
+	  if (__glibc_unlikely (ax < 0x66000000u)) /* |x| < 0x1p-24 */
-	    t = half*w;
+	    return fmaf (fabsf (x), 0x1p-25, 1.0f);
-	    return t*w;
+	  return (0.5f * x) * x + 1.0f;
 	}
-
+      static const double cp[] =
-    /* x is INF or NaN */
+	{
-	if(ix>=0x7f800000) return x*x;
+	  0x1.fffffffffffe3p-2,  0x1.55555555723cfp-5,
-
+	  0x1.6c16bee4a5986p-10, 0x1.a0483fc0328f7p-16
-    /* |x| > overflowthresold, cosh(x) overflow */
+	};
-	return math_narrow_eval (huge*huge);
+      double z2 = z * z;
+      double z4 = z2 * z2;
+      return 1.0 + z2 * ((cp[0] + z2 * cp[1]) + z4 * (cp[2] + z2 * (cp[3])));
+    }
+  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;
+  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)

sysdeps/loongarch/lp64/libm-test-ulps

@@ -621,22 +621,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: Real part of "cpow":

sysdeps/microblaze/libm-test-ulps

@@ -118,7 +118,6 @@ float: 1
 
 Function: "cosh":
 double: 1
-float: 1
 
 Function: Real part of "cpow":
 double: 2

sysdeps/mips/mips32/libm-test-ulps

@@ -489,19 +489,15 @@ float: 2
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: "cosh_upward":
 double: 2
-float: 2
 
 Function: Real part of "cpow":
 double: 2

sysdeps/mips/mips64/libm-test-ulps

@@ -621,22 +621,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: Real part of "cpow":

sysdeps/or1k/fpu/libm-test-ulps

@@ -489,19 +489,15 @@ float: 1
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: "cosh_upward":
 double: 2
-float: 2
 
 Function: Real part of "cpow":
 double: 2

sysdeps/or1k/nofpu/libm-test-ulps

@@ -505,19 +505,15 @@ float: 1
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_downward":
 double: 2
-float: 1
 
 Function: "cosh_towardzero":
 double: 2
-float: 1
 
 Function: "cosh_upward":
 double: 2
-float: 2
 
 Function: "cospi":
 double: 1

sysdeps/powerpc/fpu/libm-test-ulps

@@ -854,25 +854,21 @@ ldouble: 5
 
 Function: "cosh":
 double: 2
-float: 2
 float128: 2
 ldouble: 3
 
 Function: "cosh_downward":
 double: 3
-float: 1
 float128: 3
 ldouble: 6
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 float128: 3
 ldouble: 6
 
 Function: "cosh_upward":
 double: 2
-float: 2
 float128: 3
 ldouble: 2
 

sysdeps/powerpc/nofpu/libm-test-ulps

@@ -625,22 +625,18 @@ ldouble: 5
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 3
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 6
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 6
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 2
 
 Function: Real part of "cpow":

sysdeps/riscv/nofpu/libm-test-ulps

@@ -618,22 +618,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 1
-float: 1
 ldouble: 2
 
 Function: "cosh_towardzero":
 double: 1
-float: 1
 ldouble: 2
 
 Function: "cosh_upward":
 double: 1
-float: 2
 ldouble: 3
 
 Function: Real part of "cpow":

sysdeps/riscv/rvd/libm-test-ulps

@@ -621,22 +621,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: Real part of "cpow":

sysdeps/s390/fpu/libm-test-ulps

@@ -621,22 +621,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: "cospi":

sysdeps/sh/libm-test-ulps

@@ -241,11 +241,9 @@ float: 1
 
 Function: "cosh":
 double: 2
-float: 2
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 
 Function: Real part of "cpow":
 double: 2

sysdeps/sparc/fpu/libm-test-ulps

@@ -621,22 +621,18 @@ ldouble: 2
 
 Function: "cosh":
 double: 2
-float: 2
 ldouble: 2
 
 Function: "cosh_downward":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 ldouble: 3
 
 Function: "cospi":

sysdeps/x86_64/fpu/libm-test-ulps

@@ -1026,25 +1026,21 @@ float: 1
 
 Function: "cosh":
 double: 2
-float: 2
 float128: 2
 ldouble: 3
 
 Function: "cosh_downward":
 double: 3
-float: 1
 float128: 3
 ldouble: 3
 
 Function: "cosh_towardzero":
 double: 3
-float: 1
 float128: 3
 ldouble: 3
 
 Function: "cosh_upward":
 double: 2
-float: 2
 float128: 3
 ldouble: 3
 
@@ -1342,7 +1338,7 @@ ldouble: 2
 
 Function: Real part of "ctanh_downward":
 double: 4
-float: 2
+float: 3
 float128: 5
 ldouble: 4