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 double dx2;
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