-/* s_tanhf.c -- float version of s_tanh.c.
- */
+/* Correctly-rounded hyperbolic tangent 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.
-#if defined(LIBM_SCCS) && !defined(lint)
-static char rcsid[] = "$NetBSD: s_tanhf.c,v 1.4 1995/05/10 20:48:24 jtc Exp $";
-#endif
+The original version of this file was copied from the CORE-MATH
+project (file src/binary32/tanh/tanhf.c, revision b1ecd83).
-#include <float.h>
-#include <math.h>
-#include <math_private.h>
-#include <math-underflow.h>
-#include <libm-alias-float.h>
+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:
-static const float one=1.0, two=2.0, tiny = 1.0e-30;
-
-float __tanhf(float x)
-{
- float t,z;
- int32_t jx,ix;
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
- GET_FLOAT_WORD(jx,x);
- ix = jx&0x7fffffff;
+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.
+*/
- /* x is INF or NaN */
- if(ix>=0x7f800000) {
- if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */
- else return one/x-one; /* tanh(NaN) = NaN */
- }
+#include <math.h>
+#include <stdint.h>
+#include <libm-alias-float.h>
+#include "math_config.h"
- /* |x| < 22 */
- if (ix < 0x41b00000) { /* |x|<22 */
- if (ix == 0)
- return x; /* x == +-0 */
- if (ix<0x24000000) /* |x|<2**-55 */
- {
- math_check_force_underflow (x);
- return x*(one+x); /* tanh(small) = small */
- }
- if (ix>=0x3f800000) { /* |x|>=1 */
- t = __expm1f(two*fabsf(x));
- z = one - two/(t+two);
- } else {
- t = __expm1f(-two*fabsf(x));
- z= -t/(t+two);
- }
- /* |x| > 22, return +-1 */
- } else {
- z = one - tiny; /* raised inexact flag */
+float
+__tanhf (float x)
+{
+ double z = x;
+ uint32_t ux = asuint (x);
+ int e = (ux >> 23) & 0xff;
+ if (__glibc_unlikely (e == 0xff))
+ {
+ if (ux << 9)
+ return x + x; /* x = nan */
+ static const float ir[] = { 1.0f, -1.0f };
+ return ir[ux >> 31]; /* x = +-inf */
+ }
+ if (__glibc_unlikely (e < 115))
+ {
+ if (__glibc_unlikely (e < 102))
+ {
+ if (__glibc_unlikely ((ux << 1) == 0))
+ return x;
+ return fmaf (-x, fabsf (x), x);
}
- return (jx>=0)? z: -z;
+ float x2 = x * x;
+ return fmaf (x, -0x1.555556p-2f * x2, x);
+ }
+ if ((ux << 1) > (0x41102cb3u << 1))
+ return copysignf (1.0f, x) - copysignf (0x1p-25f, x);
+ double z2 = z * z;
+ double z4 = z2 * z2;
+ double z8 = z4 * z4;
+ static const double cn[] =
+ {
+ 0x1p+0, 0x1.30877b8b72d33p-3, 0x1.694aa09ae9e5ep-8,
+ 0x1.4101377abb729p-14, 0x1.e0392b1db0018p-22, 0x1.2533756e546f7p-30,
+ 0x1.d62e5abe6ae8ap-41, 0x1.b06be534182dep-54
+ };
+ static const double cd[] =
+ {
+ 0x1p+0, 0x1.ed99131b0ebeap-2, 0x1.0d27ed6c95a69p-5,
+ 0x1.7cbdaca0e9fccp-11, 0x1.b4e60b892578ep-18, 0x1.a6f707c5c71abp-26,
+ 0x1.35a8b6e2cd94cp-35, 0x1.ca8230677aa01p-47
+ };
+ double n0 = cn[0] + z2 * cn[1];
+ double n2 = cn[2] + z2 * cn[3];
+ double n4 = cn[4] + z2 * cn[5];
+ double n6 = cn[6] + z2 * cn[7];
+ n0 += z4 * n2;
+ n4 += z4 * n6;
+ n0 += z8 * n4;
+ double d0 = cd[0] + z2 * cd[1];
+ double d2 = cd[2] + z2 * cd[3];
+ double d4 = cd[4] + z2 * cd[5];
+ double d6 = cd[6] + z2 * cd[7];
+ d0 += z4 * d2;
+ d4 += z4 * d6;
+ d0 += z8 * d4;
+ double r = z * n0 / d0;
+ return r;
}
libm_alias_float (__tanh, tanh)
projects / thirdparty / glibc.git / commitdiff
