|
145
|
1 /*
|
|
|
2 * ====================================================
|
|
|
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
|
4 *
|
|
|
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
|
|
|
6 * Permission to use, copy, modify, and distribute this
|
|
|
7 * software is freely granted, provided that this notice
|
|
|
8 * is preserved.
|
|
|
9 * ====================================================
|
|
|
10 */
|
|
|
11
|
|
|
12 /*
|
|
|
13 Long double expansions are
|
|
|
14 Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
|
|
|
15 and are incorporated herein by permission of the author. The author
|
|
|
16 reserves the right to distribute this material elsewhere under different
|
|
|
17 copying permissions. These modifications are distributed here under
|
|
|
18 the following terms:
|
|
|
19
|
|
|
20 This library is free software; you can redistribute it and/or
|
|
|
21 modify it under the terms of the GNU Lesser General Public
|
|
|
22 License as published by the Free Software Foundation; either
|
|
|
23 version 2.1 of the License, or (at your option) any later version.
|
|
|
24
|
|
|
25 This library is distributed in the hope that it will be useful,
|
|
|
26 but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
27 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
|
28 Lesser General Public License for more details.
|
|
|
29
|
|
|
30 You should have received a copy of the GNU Lesser General Public
|
|
|
31 License along with this library; if not, see
|
|
|
32 <http://www.gnu.org/licenses/>. */
|
|
|
33
|
|
|
34 /* __quadmath_kernel_tanq( x, y, k )
|
|
|
35 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
|
|
|
36 * Input x is assumed to be bounded by ~pi/4 in magnitude.
|
|
|
37 * Input y is the tail of x.
|
|
|
38 * Input k indicates whether tan (if k=1) or
|
|
|
39 * -1/tan (if k= -1) is returned.
|
|
|
40 *
|
|
|
41 * Algorithm
|
|
|
42 * 1. Since tan(-x) = -tan(x), we need only to consider positive x.
|
|
|
43 * 2. if x < 2^-57, return x with inexact if x!=0.
|
|
|
44 * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
|
|
|
45 * on [0,0.67433].
|
|
|
46 *
|
|
|
47 * Note: tan(x+y) = tan(x) + tan'(x)*y
|
|
|
48 * ~ tan(x) + (1+x*x)*y
|
|
|
49 * Therefore, for better accuracy in computing tan(x+y), let
|
|
|
50 * r = x^3 * R(x^2)
|
|
|
51 * then
|
|
|
52 * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
|
|
|
53 *
|
|
|
54 * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then
|
|
|
55 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
|
|
|
56 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
|
|
|
57 */
|
|
|
58
|
|
|
59 #include "quadmath-imp.h"
|
|
|
60
|
|
|
61 static const __float128
|
|
|
62 one = 1,
|
|
|
63 pio4hi = 7.8539816339744830961566084581987569936977E-1Q,
|
|
|
64 pio4lo = 2.1679525325309452561992610065108379921906E-35Q,
|
|
|
65
|
|
|
66 /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
|
|
|
67 0 <= x <= 0.6743316650390625
|
|
|
68 Peak relative error 8.0e-36 */
|
|
|
69 TH = 3.333333333333333333333333333333333333333E-1Q,
|
|
|
70 T0 = -1.813014711743583437742363284336855889393E7Q,
|
|
|
71 T1 = 1.320767960008972224312740075083259247618E6Q,
|
|
|
72 T2 = -2.626775478255838182468651821863299023956E4Q,
|
|
|
73 T3 = 1.764573356488504935415411383687150199315E2Q,
|
|
|
74 T4 = -3.333267763822178690794678978979803526092E-1Q,
|
|
|
75
|
|
|
76 U0 = -1.359761033807687578306772463253710042010E8Q,
|
|
|
77 U1 = 6.494370630656893175666729313065113194784E7Q,
|
|
|
78 U2 = -4.180787672237927475505536849168729386782E6Q,
|
|
|
79 U3 = 8.031643765106170040139966622980914621521E4Q,
|
|
|
80 U4 = -5.323131271912475695157127875560667378597E2Q;
|
|
|
81 /* 1.000000000000000000000000000000000000000E0 */
|
|
|
82
|
|
|
83
|
|
|
84 __float128
|
|
|
85 __quadmath_kernel_tanq (__float128 x, __float128 y, int iy)
|
|
|
86 {
|
|
|
87 __float128 z, r, v, w, s;
|
|
|
88 int32_t ix, sign;
|
|
|
89 ieee854_float128 u, u1;
|
|
|
90
|
|
|
91 u.value = x;
|
|
|
92 ix = u.words32.w0 & 0x7fffffff;
|
|
|
93 if (ix < 0x3fc60000) /* x < 2**-57 */
|
|
|
94 {
|
|
|
95 if ((int) x == 0)
|
|
|
96 { /* generate inexact */
|
|
|
97 if ((ix | u.words32.w1 | u.words32.w2 | u.words32.w3
|
|
|
98 | (iy + 1)) == 0)
|
|
|
99 return one / fabsq (x);
|
|
|
100 else if (iy == 1)
|
|
|
101 {
|
|
|
102 math_check_force_underflow (x);
|
|
|
103 return x;
|
|
|
104 }
|
|
|
105 else
|
|
|
106 return -one / x;
|
|
|
107 }
|
|
|
108 }
|
|
|
109 if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */
|
|
|
110 {
|
|
|
111 if ((u.words32.w0 & 0x80000000) != 0)
|
|
|
112 {
|
|
|
113 x = -x;
|
|
|
114 y = -y;
|
|
|
115 sign = -1;
|
|
|
116 }
|
|
|
117 else
|
|
|
118 sign = 1;
|
|
|
119 z = pio4hi - x;
|
|
|
120 w = pio4lo - y;
|
|
|
121 x = z + w;
|
|
|
122 y = 0.0;
|
|
|
123 }
|
|
|
124 z = x * x;
|
|
|
125 r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
|
|
|
126 v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
|
|
|
127 r = r / v;
|
|
|
128
|
|
|
129 s = z * x;
|
|
|
130 r = y + z * (s * r + y);
|
|
|
131 r += TH * s;
|
|
|
132 w = x + r;
|
|
|
133 if (ix >= 0x3ffe5942)
|
|
|
134 {
|
|
|
135 v = (__float128) iy;
|
|
|
136 w = (v - 2.0 * (x - (w * w / (w + v) - r)));
|
|
|
137 /* SIGN is set for arguments that reach this code, but not
|
|
|
138 otherwise, resulting in warnings that it may be used
|
|
|
139 uninitialized although in the cases where it is used it has
|
|
|
140 always been set. */
|
|
|
141
|
|
|
142
|
|
|
143 if (sign < 0)
|
|
|
144 w = -w;
|
|
|
145
|
|
|
146 return w;
|
|
|
147 }
|
|
|
148 if (iy == 1)
|
|
|
149 return w;
|
|
|
150 else
|
|
|
151 { /* if allow error up to 2 ulp,
|
|
|
152 simply return -1.0/(x+r) here */
|
|
|
153 /* compute -1.0/(x+r) accurately */
|
|
|
154 u1.value = w;
|
|
|
155 u1.words32.w2 = 0;
|
|
|
156 u1.words32.w3 = 0;
|
|
|
157 v = r - (u1.value - x); /* u1+v = r+x */
|
|
|
158 z = -1.0 / w;
|
|
|
159 u.value = z;
|
|
|
160 u.words32.w2 = 0;
|
|
|
161 u.words32.w3 = 0;
|
|
|
162 s = 1.0 + u.value * u1.value;
|
|
|
163 return u.value + z * (s + u.value * v);
|
|
|
164 }
|
|
|
165 }
|
