![[first]](../icons/first.png){kind=icon}
![[previous]](../icons/n_left.png){kind=icon}
![[next]](../icons/n_right.png){kind=icon}
![[last]](../icons/last.png){kind=icon}
![[top]](../icons/n_top.png){kind=icon}
![[bottom]](../icons/bottom.png){kind=icon}
![[index]](../icons/index.png){kind=icon}
![[help]](../icons/help.png){kind=icon}
*/
1 /* tan (tangent) function with 'long double' argument.
2
3 Copyright (C) 2003-2021 Free Software Foundation, Inc.
4
5 This file is free software: you can redistribute it and/or modify
6 it under the terms of the GNU Lesser General Public License as
7 published by the Free Software Foundation; either version 3 of the
8 License, or (at your option) any later version.
9
10 This file is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU Lesser General Public License for more details.
14
15 You should have received a copy of the GNU Lesser General Public License
16 along with this program. If not, see <https://www.gnu.org/licenses/>. */
17
18 /* s_tanl.c -- long double version of s_tan.c.
19 * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
20 */
21
22 /* @(#)s_tan.c 5.1 93/09/24 */
23 /*
24 * ====================================================
25 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
26 *
27 * Developed at SunPro, a Sun Microsystems, Inc. business.
28 * Permission to use, copy, modify, and distribute this
29 * software is freely granted, provided that this notice
30 * is preserved.
31 * ====================================================
32 */
33
34 #include <config.h>
35
36 /* Specification. */
37 #include <math.h>
38
39 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
40
41 long double
42 tanl (long double x)
/* ![[next]](../icons/right.png)
![[first]](../icons/n_first.png)
![[top]](../icons/top.png)
*/
43 {
44 return tan (x);
45 }
46
47 #else
48
49 /* Code based on glibc/sysdeps/ieee754/ldbl-128/s_tanl.c
50 and glibc/sysdeps/ieee754/ldbl-128/k_tanl.c. */
51
52 /* tanl(x)
53 * Return tangent function of x.
54 *
55 * kernel function:
56 * __kernel_tanl ... tangent function on [-pi/4,pi/4]
57 * __ieee754_rem_pio2l ... argument reduction routine
58 *
59 * Method.
60 * Let S,C and T denote the sin, cos and tan respectively on
61 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
62 * in [-pi/4 , +pi/4], and let n = k mod 4.
63 * We have
64 *
65 * n sin(x) cos(x) tan(x)
66 * ----------------------------------------------------------
67 * 0 S C T
68 * 1 C -S -1/T
69 * 2 -S -C T
70 * 3 -C S -1/T
71 * ----------------------------------------------------------
72 *
73 * Special cases:
74 * Let trig be any of sin, cos, or tan.
75 * trig(+-INF) is NaN, with signals;
76 * trig(NaN) is that NaN;
77 *
78 * Accuracy:
79 * TRIG(x) returns trig(x) nearly rounded
80 */
81
82 # include "trigl.h"
83
84 /*
85 * ====================================================
86 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
87 *
88 * Developed at SunPro, a Sun Microsystems, Inc. business.
89 * Permission to use, copy, modify, and distribute this
90 * software is freely granted, provided that this notice
91 * is preserved.
92 * ====================================================
93 */
94
95 /*
96 Long double expansions contributed by
97 Stephen L. Moshier <moshier@na-net.ornl.gov>
98 */
99
100 /* __kernel_tanl( x, y, k )
101 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
102 * Input x is assumed to be bounded by ~pi/4 in magnitude.
103 * Input y is the tail of x.
104 * Input k indicates whether tan (if k=1) or
105 * -1/tan (if k= -1) is returned.
106 *
107 * Algorithm
108 * 1. Since tan(-x) = -tan(x), we need only to consider positive x.
109 * 2. if x < 2^-57, return x with inexact if x!=0.
110 * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
111 * on [0,0.67433].
112 *
113 * Note: tan(x+y) = tan(x) + tan'(x)*y
114 * ~ tan(x) + (1+x*x)*y
115 * Therefore, for better accuracy in computing tan(x+y), let
116 * r = x^3 * R(x^2)
117 * then
118 * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
119 *
120 * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then
121 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
122 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
123 */
124
125
126 static const long double
127 pio4hi = 7.8539816339744830961566084581987569936977E-1L,
128 pio4lo = 2.1679525325309452561992610065108379921906E-35L,
129
130 /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
131 0 <= x <= 0.6743316650390625
132 Peak relative error 8.0e-36 */
133 TH = 3.333333333333333333333333333333333333333E-1L,
134 T0 = -1.813014711743583437742363284336855889393E7L,
135 T1 = 1.320767960008972224312740075083259247618E6L,
136 T2 = -2.626775478255838182468651821863299023956E4L,
137 T3 = 1.764573356488504935415411383687150199315E2L,
138 T4 = -3.333267763822178690794678978979803526092E-1L,
139
140 U0 = -1.359761033807687578306772463253710042010E8L,
141 U1 = 6.494370630656893175666729313065113194784E7L,
142 U2 = -4.180787672237927475505536849168729386782E6L,
143 U3 = 8.031643765106170040139966622980914621521E4L,
144 U4 = -5.323131271912475695157127875560667378597E2L;
145 /* 1.000000000000000000000000000000000000000E0 */
146
147
148 static long double
149 kernel_tanl (long double x, long double y, int iy)
/* ![[previous]](../icons/left.png)
*/
150 {
151 long double z, r, v, w, s, u, u1;
152 int invert = 0, sign;
153
154 sign = 1;
155 if (x < 0)
156 {
157 x = -x;
158 y = -y;
159 sign = -1;
160 }
161
162 if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */
163 {
164 if ((int) x == 0)
165 { /* generate inexact */
166 if (iy == -1 && x == 0.0)
167 return 1.0L / fabs (x);
168 else
169 return (iy == 1) ? x : -1.0L / x;
170 }
171 }
172 if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */
173 {
174 invert = 1;
175
176 z = pio4hi - x;
177 w = pio4lo - y;
178 x = z + w;
179 y = 0.0;
180 }
181 z = x * x;
182 r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
183 v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
184 r = r / v;
185
186 s = z * x;
187 r = y + z * (s * r + y);
188 r += TH * s;
189 w = x + r;
190 if (invert)
191 {
192 v = (long double) iy;
193 w = (v - 2.0 * (x - (w * w / (w + v) - r)));
194 if (sign < 0)
195 w = -w;
196 return w;
197 }
198 if (iy == 1)
199 return w;
200 else
201 { /* if allow error up to 2 ulp,
202 simply return -1.0/(x+r) here */
203 /* compute -1.0/(x+r) accurately */
204 u1 = (double) w;
205 v = r - (u1 - x);
206 z = -1.0 / w;
207 u = (double) z;
208 s = 1.0 + u * u1;
209 return u + z * (s + u * v);
210 }
211 }
212
213 long double
214 tanl (long double x)
/* */
215 {
216 long double y[2], z = 0.0L;
217 int n;
218
219 /* tanl(NaN) is NaN */
220 if (isnanl (x))
221 return x;
222
223 /* |x| ~< pi/4 */
224 if (x >= -0.7853981633974483096156608458198757210492 &&
225 x <= 0.7853981633974483096156608458198757210492)
226 return kernel_tanl (x, z, 1);
227
228 /* tanl(Inf) is NaN, tanl(0) is 0 */
229 else if (x + x == x)
230 return x - x; /* NaN */
231
232 /* argument reduction needed */
233 else
234 {
235 n = ieee754_rem_pio2l (x, y);
236 /* 1 -- n even, -1 -- n odd */
237 return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
238 }
239 }
240
241 #endif
242
243 #if 0
244 int
245 main (void)
/* ![[last]](../icons/n_last.png)
*/
246 {
247 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492));
248 printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492));
249 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 *3));
250 printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492 *31));
251 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 / 2));
252 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 3/2));
253 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 5/2));
254 }
255 #endif
![[bottom]](../icons/n_bottom.png){kind=icon}
*/