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