Back to index

glibc  2.9
k_tanl.c
Go to the documentation of this file.
00001 /*
00002  * ====================================================
00003  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00004  *
00005  * Developed at SunPro, a Sun Microsystems, Inc. business.
00006  * Permission to use, copy, modify, and distribute this
00007  * software is freely granted, provided that this notice
00008  * is preserved.
00009  * ====================================================
00010  */
00011 
00012 /*
00013   Long double expansions are
00014   Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
00015   and are incorporated herein by permission of the author.  The author
00016   reserves the right to distribute this material elsewhere under different
00017   copying permissions.  These modifications are distributed here under
00018   the following terms:
00019 
00020     This library is free software; you can redistribute it and/or
00021     modify it under the terms of the GNU Lesser General Public
00022     License as published by the Free Software Foundation; either
00023     version 2.1 of the License, or (at your option) any later version.
00024 
00025     This library is distributed in the hope that it will be useful,
00026     but WITHOUT ANY WARRANTY; without even the implied warranty of
00027     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00028     Lesser General Public License for more details.
00029 
00030     You should have received a copy of the GNU Lesser General Public
00031     License along with this library; if not, write to the Free Software
00032     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
00033 
00034 /* __kernel_tanl( x, y, k )
00035  * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
00036  * Input x is assumed to be bounded by ~pi/4 in magnitude.
00037  * Input y is the tail of x.
00038  * Input k indicates whether tan (if k=1) or
00039  * -1/tan (if k= -1) is returned.
00040  *
00041  * Algorithm
00042  *     1. Since tan(-x) = -tan(x), we need only to consider positive x.
00043  *     2. if x < 2^-57, return x with inexact if x!=0.
00044  *     3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
00045  *          on [0,0.67433].
00046  *
00047  *        Note: tan(x+y) = tan(x) + tan'(x)*y
00048  *                      ~ tan(x) + (1+x*x)*y
00049  *        Therefore, for better accuracy in computing tan(x+y), let
00050  *            r = x^3 * R(x^2)
00051  *        then
00052  *            tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
00053  *
00054  *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
00055  *            tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
00056  *                   = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
00057  */
00058 
00059 #include "math.h"
00060 #include "math_private.h"
00061 #ifdef __STDC__
00062 static const long double
00063 #else
00064 static long double
00065 #endif
00066   one = 1.0L,
00067   pio4hi = 7.8539816339744830961566084581987569936977E-1L,
00068   pio4lo = 2.1679525325309452561992610065108379921906E-35L,
00069 
00070   /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
00071      0 <= x <= 0.6743316650390625
00072      Peak relative error 8.0e-36  */
00073  TH =  3.333333333333333333333333333333333333333E-1L,
00074  T0 = -1.813014711743583437742363284336855889393E7L,
00075  T1 =  1.320767960008972224312740075083259247618E6L,
00076  T2 = -2.626775478255838182468651821863299023956E4L,
00077  T3 =  1.764573356488504935415411383687150199315E2L,
00078  T4 = -3.333267763822178690794678978979803526092E-1L,
00079 
00080  U0 = -1.359761033807687578306772463253710042010E8L,
00081  U1 =  6.494370630656893175666729313065113194784E7L,
00082  U2 = -4.180787672237927475505536849168729386782E6L,
00083  U3 =  8.031643765106170040139966622980914621521E4L,
00084  U4 = -5.323131271912475695157127875560667378597E2L;
00085   /* 1.000000000000000000000000000000000000000E0 */
00086 
00087 
00088 #ifdef __STDC__
00089 long double
00090 __kernel_tanl (long double x, long double y, int iy)
00091 #else
00092 long double
00093 __kernel_tanl (x, y, iy)
00094      long double x, y;
00095      int iy;
00096 #endif
00097 {
00098   long double z, r, v, w, s;
00099   int32_t ix, sign;
00100   ieee854_long_double_shape_type u, u1;
00101 
00102   u.value = x;
00103   ix = u.parts32.w0 & 0x7fffffff;
00104   if (ix < 0x3c600000)             /* x < 2**-57 */
00105     {
00106       if ((int) x == 0)
00107        {                    /* generate inexact */
00108          if ((ix | u.parts32.w1 | (u.parts32.w2 & 0x7fffffff) | u.parts32.w3
00109               | (iy + 1)) == 0)
00110            return one / fabs (x);
00111          else
00112            return (iy == 1) ? x : -one / x;
00113        }
00114     }
00115   if (ix >= 0x3fe59420) /* |x| >= 0.6743316650390625 */
00116     {
00117       if ((u.parts32.w0 & 0x80000000) != 0)
00118        {
00119          x = -x;
00120          y = -y;
00121          sign = -1;
00122        }
00123       else
00124        sign = 1;
00125       z = pio4hi - x;
00126       w = pio4lo - y;
00127       x = z + w;
00128       y = 0.0;
00129     }
00130   z = x * x;
00131   r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
00132   v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
00133   r = r / v;
00134 
00135   s = z * x;
00136   r = y + z * (s * r + y);
00137   r += TH * s;
00138   w = x + r;
00139   if (ix >= 0x3fe59420)
00140     {
00141       v = (long double) iy;
00142       w = (v - 2.0 * (x - (w * w / (w + v) - r)));
00143       if (sign < 0)
00144        w = -w;
00145       return w;
00146     }
00147   if (iy == 1)
00148     return w;
00149   else
00150     {                       /* if allow error up to 2 ulp,
00151                                simply return -1.0/(x+r) here */
00152       /*  compute -1.0/(x+r) accurately */
00153       u1.value = w;
00154       u1.parts32.w2 = 0;
00155       u1.parts32.w3 = 0;
00156       v = r - (u1.value - x);             /* u1+v = r+x */
00157       z = -1.0 / w;
00158       u.value = z;
00159       u.parts32.w2 = 0;
00160       u.parts32.w3 = 0;
00161       s = 1.0 + u.value * u1.value;
00162       return u.value + z * (s + u.value * v);
00163     }
00164 }