tan implementation for 128-bit long double.
[kopensolaris-gnu/glibc.git] / sysdeps / ieee754 / ldbl-128 / k_tanl.c
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 contributed by
14   Stephen L. Moshier <moshier@na-net.ornl.gov>
15 */
16
17 /* __kernel_tanl( x, y, k )
18  * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
19  * Input x is assumed to be bounded by ~pi/4 in magnitude.
20  * Input y is the tail of x.
21  * Input k indicates whether tan (if k=1) or
22  * -1/tan (if k= -1) is returned.
23  *
24  * Algorithm
25  *      1. Since tan(-x) = -tan(x), we need only to consider positive x.
26  *      2. if x < 2^-57, return x with inexact if x!=0.
27  *      3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
28  *          on [0,0.67433].
29  *
30  *         Note: tan(x+y) = tan(x) + tan'(x)*y
31  *                        ~ tan(x) + (1+x*x)*y
32  *         Therefore, for better accuracy in computing tan(x+y), let
33  *              r = x^3 * R(x^2)
34  *         then
35  *              tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
36  *
37  *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
38  *              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
39  *                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
40  */
41
42 #include "math.h"
43 #include "math_private.h"
44 #ifdef __STDC__
45 static const long double
46 #else
47 static long double
48 #endif
49   one = 1.0L,
50   pio4hi = 7.8539816339744830961566084581987569936977E-1L,
51   pio4lo = 2.1679525325309452561992610065108379921906E-35L,
52
53   /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
54      0 <= x <= 0.6743316650390625
55      Peak relative error 8.0e-36  */
56  TH =  3.333333333333333333333333333333333333333E-1L,
57  T0 = -1.813014711743583437742363284336855889393E7L,
58  T1 =  1.320767960008972224312740075083259247618E6L,
59  T2 = -2.626775478255838182468651821863299023956E4L,
60  T3 =  1.764573356488504935415411383687150199315E2L,
61  T4 = -3.333267763822178690794678978979803526092E-1L,
62
63  U0 = -1.359761033807687578306772463253710042010E8L,
64  U1 =  6.494370630656893175666729313065113194784E7L,
65  U2 = -4.180787672237927475505536849168729386782E6L,
66  U3 =  8.031643765106170040139966622980914621521E4L,
67  U4 = -5.323131271912475695157127875560667378597E2L;
68   /* 1.000000000000000000000000000000000000000E0 */
69
70
71 #ifdef __STDC__
72 long double
73 __kernel_tanl (long double x, long double y, int iy)
74 #else
75 long double
76 __kernel_tanl (x, y, iy)
77      long double x, y;
78      int iy;
79 #endif
80 {
81   long double z, r, v, w, s;
82   int32_t ix, sign;
83   ieee854_long_double_shape_type u, u1;
84
85   u.value = x;
86   ix = u.parts32.w0 & 0x7fffffff;
87   if (ix < 0x3fc60000)          /* x < 2**-57 */
88     {
89       if ((int) x == 0)
90         {                       /* generate inexact */
91           if ((ix | u.parts32.w1 | u.parts32.w2 | u.parts32.w3
92                | (iy + 1)) == 0)
93             return one / fabs (x);
94           else
95             return (iy == 1) ? x : -one / x;
96         }
97     }
98   if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */
99     {
100       if ((u.parts32.w0 & 0x80000000) != 0)
101         {
102           x = -x;
103           y = -y;
104           sign = -1;
105         }
106       else
107         sign = 1;
108       z = pio4hi - x;
109       w = pio4lo - y;
110       x = z + w;
111       y = 0.0;
112     }
113   z = x * x;
114   r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
115   v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
116   r = r / v;
117
118   s = z * x;
119   r = y + z * (s * r + y);
120   r += TH * s;
121   w = x + r;
122   if (ix >= 0x3ffe5942)
123     {
124       v = (long double) iy;
125       w = (v - 2.0 * (x - (w * w / (w + v) - r)));
126       if (sign < 0)
127         w = -w;
128       return w;
129     }
130   if (iy == 1)
131     return w;
132   else
133     {                           /* if allow error up to 2 ulp,
134                                    simply return -1.0/(x+r) here */
135       /*  compute -1.0/(x+r) accurately */
136       u1.value = w;
137       u1.parts32.w2 = 0;
138       u1.parts32.w3 = 0;
139       v = r - (u1.value - x);           /* u1+v = r+x */
140       z = -1.0 / w;
141       u.value = z;
142       u.parts32.w2 = 0;
143       u.parts32.w3 = 0;
144       s = 1.0 + u.value * u1.value;
145       return u.value + z * (s + u.value * v);
146     }
147 }