1 /* from: FreeBSD: head/lib/msun/src/s_tanhl.c XXX */
2
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14 /*
15 * See s_tanh.c for complete comments.
16 *
17 * Converted to long double by Bruce D. Evans.
18 */
19
20 #include <float.h>
21 #ifdef __i386__
22 #include <ieeefp.h>
23 #endif
24
25 #include "math.h"
26 #include "math_private.h"
27 #include "fpmath.h"
28 #include "k_expl.h"
29
30 #if LDBL_MAX_EXP != 0x4000
31 /* We also require the usual expsign encoding. */
32 #error "Unsupported long double format"
33 #endif
34
35 #define BIAS (LDBL_MAX_EXP - 1)
36
37 static const volatile double tiny = 1.0e-300;
38 static const double one = 1.0;
39 #if LDBL_MANT_DIG == 64
40 /*
41 * Domain [-0.25, 0.25], range ~[-1.6304e-22, 1.6304e-22]:
42 * |tanh(x)/x - t(x)| < 2**-72.3
43 */
44 static const union IEEEl2bits
45 T3u = LD80C(0xaaaaaaaaaaaaaa9f, -2, -3.33333333333333333017e-1L);
46 #define T3 T3u.e
47 static const double
48 T5 = 1.3333333333333314e-1, /* 0x1111111111110a.0p-55 */
49 T7 = -5.3968253968210485e-2, /* -0x1ba1ba1ba1a1a1.0p-57 */
50 T9 = 2.1869488531393817e-2, /* 0x1664f488172022.0p-58 */
51 T11 = -8.8632352345964591e-3, /* -0x1226e34bc138d5.0p-59 */
52 T13 = 3.5921169709993771e-3, /* 0x1d6d371d3e400f.0p-61 */
53 T15 = -1.4555786415756001e-3, /* -0x17d923aa63814d.0p-62 */
54 T17 = 5.8645267876296793e-4, /* 0x13378589b85aa7.0p-63 */
55 T19 = -2.1121033571392224e-4; /* -0x1baf0af80c4090.0p-65 */
56 #elif LDBL_MANT_DIG == 113
57 /*
58 * Domain [-0.25, 0.25], range ~[-2.4211e-37, 2.4211e-37]:
59 * |tanh(x)/x - t(x)| < 2**121.6
60 */
61 static const long double
62 T3 = -3.33333333333333333333333333333332980e-1L, /* -0x1555555555555555555555555554e.0p-114L */
63 T5 = 1.33333333333333333333333333332707260e-1L, /* 0x1111111111111111111111110ab7b.0p-115L */
64 T7 = -5.39682539682539682539682535723482314e-2L, /* -0x1ba1ba1ba1ba1ba1ba1ba17b5fc98.0p-117L */
65 T9 = 2.18694885361552028218693591149061717e-2L, /* 0x1664f4882c10f9f32d6b1a12a25e5.0p-118L */
66 T11 = -8.86323552990219656883762347736381851e-3L, /* -0x1226e355e6c23c8f5a5a0f386cb4d.0p-119L */
67 T13 = 3.59212803657248101358314398220822722e-3L, /* 0x1d6d3d0e157ddfb403ad3637442c6.0p-121L */
68 T15 = -1.45583438705131796512568010348874662e-3L; /* -0x17da36452b75e150c44cc34253b34.0p-122L */
69 static const double
70 T17 = 5.9002744094556621e-4, /* 0x1355824803668e.0p-63 */
71 T19 = -2.3912911424260516e-4, /* -0x1f57d7734c8dde.0p-65 */
72 T21 = 9.6915379535512898e-5, /* 0x1967e18ad6a6ca.0p-66 */
73 T23 = -3.9278322983156353e-5, /* -0x1497d8e6b75729.0p-67 */
74 T25 = 1.5918887220143869e-5, /* 0x10b1319998cafa.0p-68 */
75 T27 = -6.4514295231630956e-6, /* -0x1b0f2b71b218eb.0p-70 */
76 T29 = 2.6120754043964365e-6, /* 0x15e963a3cf3a39.0p-71 */
77 T31 = -1.0407567231003314e-6, /* -0x1176041e656869.0p-72 */
78 T33 = 3.4744117554063574e-7; /* 0x1750fe732cab9c.0p-74 */
79 #endif /* LDBL_MANT_DIG == 64 */
80
81 static inline long double
divl(long double a,long double b,long double c,long double d,long double e,long double f)82 divl(long double a, long double b, long double c, long double d,
83 long double e, long double f)
84 {
85 long double inv, r;
86 float fr, fw;
87
88 _2sumF(a, c);
89 b = b + c;
90 _2sumF(d, f);
91 e = e + f;
92
93 inv = 1 / (d + e);
94
95 r = (a + b) * inv;
96 fr = r;
97 r = fr;
98
99 fw = d + e;
100 e = d - fw + e;
101 d = fw;
102
103 r = r + (a - d * r + b - e * r) * inv;
104
105 return r;
106 }
107
108 long double
tanhl(long double x)109 tanhl(long double x)
110 {
111 long double hi,lo,s,x2,x4,z;
112 #if LDBL_MANT_DIG == 113
113 double dx2;
114 #endif
115 int16_t jx,ix;
116
117 GET_LDBL_EXPSIGN(jx,x);
118 ix = jx&0x7fff;
119
120 /* x is INF or NaN */
121 if(ix>=0x7fff) {
122 if (jx>=0) return one/x+one; /* tanh(+-inf)=+-1 */
123 else return one/x-one; /* tanh(NaN) = NaN */
124 }
125
126 ENTERI();
127
128 /* |x| < 40 */
129 if (ix < 0x4004 || fabsl(x) < 40) { /* |x|<40 */
130 if (__predict_false(ix<BIAS-(LDBL_MANT_DIG+1)/2)) { /* |x|<TINY */
131 /* tanh(+-0) = +0; tanh(tiny) = tiny(-+) with inexact: */
132 return (x == 0 ? x : (0x1p200 * x - x) * 0x1p-200);
133 }
134 if (ix<0x3ffd) { /* |x|<0.25 */
135 x2 = x*x;
136 #if LDBL_MANT_DIG == 64
137 x4 = x2*x2;
138 RETURNI(((T19*x2 + T17)*x4 + (T15*x2 + T13))*(x2*x*x2*x4*x4) +
139 ((T11*x2 + T9)*x4 + (T7*x2 + T5))*(x2*x*x2) +
140 T3*(x2*x) + x);
141 #elif LDBL_MANT_DIG == 113
142 dx2 = x2;
143 #if 0
144 RETURNI(((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
145 T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
146 T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
147 (x2*x*x2) +
148 T3*(x2*x) + x);
149 #else
150 long double q = ((((((((((((((T33*dx2 + T31)*dx2 + T29)*dx2 + T27)*dx2 +
151 T25)*x2 + T23)*x2 + T21)*x2 + T19)*x2 + T17)*x2 +
152 T15)*x2 + T13)*x2 + T11)*x2 + T9)*x2 + T7)*x2 + T5)*
153 (x2*x*x2);
154 RETURNI(q + T3*(x2*x) + x);
155 #endif
156 #endif
157 }
158 k_hexpl(2*fabsl(x), &hi, &lo);
159 if (ix<0x4001 && fabsl(x) < 1.5) /* |x|<1.5 */
160 z = divl(hi, lo, -0.5, hi, lo, 0.5);
161 else
162 z = one - one/(lo+0.5+hi);
163 /* |x| >= 40, return +-1 */
164 } else {
165 z = one - tiny; /* raise inexact flag */
166 }
167 s = 1;
168 if (jx<0) s = -1;
169 RETURNI(s*z);
170 }
171