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