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