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 | <https://www.gnu.org/licenses/>. */ |
33 | |
34 | /* __kernel_tanl( 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 <float.h> |
60 | #include <math.h> |
61 | #include <math_private.h> |
62 | #include <math-underflow.h> |
63 | #include <libc-diag.h> |
64 | |
65 | static const long double |
66 | one = 1.0L, |
67 | pio4hi = 7.8539816339744830961566084581987569936977E-1L, |
68 | pio4lo = 2.1679525325309452561992610065108379921906E-35L, |
69 | |
70 | /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2) |
71 | 0 <= x <= 0.6743316650390625 |
72 | Peak relative error 8.0e-36 */ |
73 | TH = 3.333333333333333333333333333333333333333E-1L, |
74 | T0 = -1.813014711743583437742363284336855889393E7L, |
75 | T1 = 1.320767960008972224312740075083259247618E6L, |
76 | T2 = -2.626775478255838182468651821863299023956E4L, |
77 | T3 = 1.764573356488504935415411383687150199315E2L, |
78 | T4 = -3.333267763822178690794678978979803526092E-1L, |
79 | |
80 | U0 = -1.359761033807687578306772463253710042010E8L, |
81 | U1 = 6.494370630656893175666729313065113194784E7L, |
82 | U2 = -4.180787672237927475505536849168729386782E6L, |
83 | U3 = 8.031643765106170040139966622980914621521E4L, |
84 | U4 = -5.323131271912475695157127875560667378597E2L; |
85 | /* 1.000000000000000000000000000000000000000E0 */ |
86 | |
87 | |
88 | long double |
89 | __kernel_tanl (long double x, long double y, int iy) |
90 | { |
91 | long double z, r, v, w, s; |
92 | int32_t ix, sign, hx, lx; |
93 | double xhi; |
94 | |
95 | xhi = ldbl_high (x); |
96 | EXTRACT_WORDS (hx, lx, xhi); |
97 | ix = hx & 0x7fffffff; |
98 | if (ix < 0x3c600000) /* x < 2**-57 */ |
99 | { |
100 | if ((int) x == 0) /* generate inexact */ |
101 | { |
102 | if ((ix | lx | (iy + 1)) == 0) |
103 | return one / fabs (x: x); |
104 | else if (iy == 1) |
105 | { |
106 | math_check_force_underflow (x); |
107 | return x; |
108 | } |
109 | else |
110 | return -one / x; |
111 | } |
112 | } |
113 | if (ix >= 0x3fe59420) /* |x| >= 0.6743316650390625 */ |
114 | { |
115 | if ((hx & 0x80000000) != 0) |
116 | { |
117 | x = -x; |
118 | y = -y; |
119 | sign = -1; |
120 | } |
121 | else |
122 | sign = 1; |
123 | z = pio4hi - x; |
124 | w = pio4lo - y; |
125 | x = z + w; |
126 | y = 0.0; |
127 | } |
128 | z = x * x; |
129 | r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4))); |
130 | v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z)))); |
131 | r = r / v; |
132 | |
133 | s = z * x; |
134 | r = y + z * (s * r + y); |
135 | r += TH * s; |
136 | w = x + r; |
137 | if (ix >= 0x3fe59420) |
138 | { |
139 | v = (long double) iy; |
140 | w = (v - 2.0 * (x - (w * w / (w + v) - r))); |
141 | /* SIGN is set for arguments that reach this code, but not |
142 | otherwise, resulting in warnings that it may be used |
143 | uninitialized although in the cases where it is used it has |
144 | always been set. */ |
145 | DIAG_PUSH_NEEDS_COMMENT; |
146 | DIAG_IGNORE_NEEDS_COMMENT (5, "-Wmaybe-uninitialized" ); |
147 | if (sign < 0) |
148 | w = -w; |
149 | DIAG_POP_NEEDS_COMMENT; |
150 | return w; |
151 | } |
152 | if (iy == 1) |
153 | return w; |
154 | else |
155 | { /* if allow error up to 2 ulp, |
156 | simply return -1.0/(x+r) here */ |
157 | /* compute -1.0/(x+r) accurately */ |
158 | long double u1, z1; |
159 | |
160 | u1 = ldbl_high (w); |
161 | v = r - (u1 - x); /* u1+v = r+x */ |
162 | z = -1.0 / w; |
163 | z1 = ldbl_high (z); |
164 | s = 1.0 + z1 * u1; |
165 | return z1 + z * (s + z1 * v); |
166 | } |
167 | } |
168 | |