-/* origin: FreeBSD /usr/src/lib/msun/src/s_tanhf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
- */
-/*
- * ====================================================
- * 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.
- * ====================================================
- */
-
#include "libm.h"
-static const float one = 1.0, two = 2.0, tiny = 1.0e-30, huge = 1.0e30;
-
float tanhf(float x)
{
- float t,z;
- int32_t jx,ix;
+ union {float f; uint32_t i;} u = {.f = x};
+ uint32_t w;
+ int sign;
+ float t;
- GET_FLOAT_WORD(jx, x);
- ix = jx & 0x7fffffff;
+ /* x = |x| */
+ sign = u.i >> 31;
+ u.i &= 0x7fffffff;
+ x = u.f;
+ w = u.i;
- /* 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 */
- }
-
- if (ix < 0x41100000) { /* |x| < 9 */
- if (ix < 0x39800000) { /* |x| < 2**-12 */
- /* tanh(tiny) = tiny with inexact */
- if (huge+x > one)
- return x;
- }
- if (ix >= 0x3f800000) { /* |x|>=1 */
- t = expm1f(two*fabsf(x));
- z = one - two/(t+two);
+ if (w > 0x3f0c9f54) {
+ /* |x| > log(3)/2 ~= 0.5493 or nan */
+ if (w > 0x41200000) {
+ /* |x| > 10 */
+ t = 1 + 0/x;
} else {
- t = expm1f(-two*fabsf(x));
- z = -t/(t+two);
+ t = expm1f(2*x);
+ t = 1 - 2/(t+2);
}
- } else { /* |x| >= 9, return +-1 */
- z = one - tiny; /* raise inexact */
+ } else if (w > 0x3e82c578) {
+ /* |x| > log(5/3)/2 ~= 0.2554 */
+ t = expm1f(2*x);
+ t = t/(t+2);
+ } else if (w >= 0x00800000) {
+ /* |x| >= 0x1p-126 */
+ t = expm1f(-2*x);
+ t = -t/(t+2);
+ } else {
+ /* |x| is subnormal */
+ FORCE_EVAL(x*x);
+ t = x;
}
- return jx >= 0 ? z : -z;
+ return sign ? -t : t;
}