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 }