1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_tanhl.c */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
6  * Developed at SunPro, a Sun Microsystems, Inc. business.
7  * Permission to use, copy, modify, and distribute this
8  * software is freely granted, provided that this notice
9  * is preserved.
10  * ====================================================
11  */
12 /* tanhl(x)
13  * Return the Hyperbolic Tangent of x
14  *
15  * Method :
16  *                                      x    -x
17  *                                     e  - e
18  *      0. tanhl(x) is defined to be -----------
19  *                                      x    -x
20  *                                     e  + e
21  *      1. reduce x to non-negative by tanhl(-x) = -tanhl(x).
22  *      2.  0      <= x <= 2**-55 : tanhl(x) := x*(one+x)
23  *                                               -t
24  *          2**-55 <  x <=  1     : tanhl(x) := -----; t = expm1l(-2x)
25  *                                              t + 2
26  *                                                    2
27  *          1      <= x <=  23.0  : tanhl(x) := 1-  ----- ; t=expm1l(2x)
28  *                                                  t + 2
29  *          23.0   <  x <= INF    : tanhl(x) := 1.
30  *
31  * Special cases:
32  *      tanhl(NaN) is NaN;
33  *      only tanhl(0)=0 is exact for finite argument.
34  */
36 #include "libm.h"
38 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
39 long double tanhl(long double x)
40 {
41         return tanh(x);
42 }
43 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
44 static const long double one=1.0, two=2.0, tiny = 1.0e-4900L;
46 long double tanhl(long double x)
47 {
48         long double t,z;
49         int32_t se;
50         uint32_t jj0,jj1,ix;
52         /* High word of |x|. */
53         GET_LDOUBLE_WORDS(se, jj0, jj1, x);
54         ix = se & 0x7fff;
56         /* x is INF or NaN */
57         if (ix == 0x7fff) {
58                 /* for NaN it's not important which branch: tanhl(NaN) = NaN */
59                 if (se & 0x8000)
60                         return one/x-one;  /* tanhl(-inf)= -1; */
61                 return one/x+one;          /* tanhl(+inf)= +1 */
62         }
64         /* |x| < 23 */
65         if (ix < 0x4003 || (ix == 0x4003 && jj0 < 0xb8000000u)) {
66                 if ((ix|jj0|jj1) == 0) /* x == +- 0 */
67                         return x;
68                 if (ix < 0x3fc8)       /* |x| < 2**-55 */
69                         return x*(one+tiny);  /* tanh(small) = small */
70                 if (ix >= 0x3fff) {    /* |x| >= 1  */
71                         t = expm1l(two*fabsl(x));
72                         z = one - two/(t+two);
73                 } else {
74                         t = expm1l(-two*fabsl(x));
75                         z = -t/(t+two);
76                 }
77         /* |x| > 23, return +-1 */
78         } else {
79                 z = one - tiny;  /* raise inexact flag */
80         }
81         return se & 0x8000 ? -z : z;
82 }
83 #endif