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