1 /* s_tanl.c -- long double version of s_tan.c.
2 * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
5 /* @(#)s_tan.c 5.1 93/09/24 */
7 * ====================================================
8 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
10 * Developed at SunPro, a Sun Microsystems, Inc. business.
11 * Permission to use, copy, modify, and distribute this
12 * software is freely granted, provided that this notice
14 * ====================================================
22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
33 * Return tangent function of x.
36 * __kernel_tanl ... tangent function on [-pi/4,pi/4]
37 * __ieee754_rem_pio2l ... argument reduction routine
40 * Let S,C and T denote the sin, cos and tan respectively on
41 * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
42 * in [-pi/4 , +pi/4], and let n = k mod 4.
45 * n sin(x) cos(x) tan(x)
46 * ----------------------------------------------------------
51 * ----------------------------------------------------------
54 * Let trig be any of sin, cos, or tan.
55 * trig(+-INF) is NaN, with signals;
56 * trig(NaN) is that NaN;
59 * TRIG(x) returns trig(x) nearly rounded
65 * ====================================================
66 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
68 * Developed at SunPro, a Sun Microsystems, Inc. business.
69 * Permission to use, copy, modify, and distribute this
70 * software is freely granted, provided that this notice
72 * ====================================================
76 Long double expansions contributed by
77 Stephen L. Moshier <moshier@na-net.ornl.gov>
80 /* __kernel_tanl( x, y, k )
81 * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
82 * Input x is assumed to be bounded by ~pi/4 in magnitude.
83 * Input y is the tail of x.
84 * Input k indicates whether tan (if k=1) or
85 * -1/tan (if k= -1) is returned.
88 * 1. Since tan(-x) = -tan(x), we need only to consider positive x.
89 * 2. if x < 2^-57, return x with inexact if x!=0.
90 * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
93 * Note: tan(x+y) = tan(x) + tan'(x)*y
94 * ~ tan(x) + (1+x*x)*y
95 * Therefore, for better accuracy in computing tan(x+y), let
98 * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
100 * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then
101 * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
102 * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
106 static const long double
107 pio4hi = 7.8539816339744830961566084581987569936977E-1L,
108 pio4lo = 2.1679525325309452561992610065108379921906E-35L,
110 /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
111 0 <= x <= 0.6743316650390625
112 Peak relative error 8.0e-36 */
113 TH = 3.333333333333333333333333333333333333333E-1L,
114 T0 = -1.813014711743583437742363284336855889393E7L,
115 T1 = 1.320767960008972224312740075083259247618E6L,
116 T2 = -2.626775478255838182468651821863299023956E4L,
117 T3 = 1.764573356488504935415411383687150199315E2L,
118 T4 = -3.333267763822178690794678978979803526092E-1L,
120 U0 = -1.359761033807687578306772463253710042010E8L,
121 U1 = 6.494370630656893175666729313065113194784E7L,
122 U2 = -4.180787672237927475505536849168729386782E6L,
123 U3 = 8.031643765106170040139966622980914621521E4L,
124 U4 = -5.323131271912475695157127875560667378597E2L;
125 /* 1.000000000000000000000000000000000000000E0 */
129 kernel_tanl (long double x, long double y, int iy)
131 long double z, r, v, w, s, u, u1;
132 int invert = 0, sign;
142 if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */
145 { /* generate inexact */
146 if (iy == -1 && x == 0.0)
147 return 1.0L / fabs (x);
149 return (iy == 1) ? x : -1.0L / x;
152 if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */
162 r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
163 v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
167 r = y + z * (s * r + y);
172 v = (long double) iy;
173 w = (v - 2.0 * (x - (w * w / (w + v) - r)));
181 { /* if allow error up to 2 ulp,
182 simply return -1.0/(x+r) here */
183 /* compute -1.0/(x+r) accurately */
189 return u + z * (s + u * v);
196 long double y[2], z = 0.0L;
199 /* tanl(NaN) is NaN */
204 if (x >= -0.7853981633974483096156608458198757210492 &&
205 x <= 0.7853981633974483096156608458198757210492)
206 return kernel_tanl (x, z, 1);
208 /* tanl(Inf) is NaN, tanl(0) is 0 */
210 return x - x; /* NaN */
212 /* argument reduction needed */
215 n = ieee754_rem_pio2l (x, y);
216 /* 1 -- n even, -1 -- n odd */
217 return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
227 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492));
228 printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492));
229 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 *3));
230 printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492 *31));
231 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 / 2));
232 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 3/2));
233 printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 5/2));