lib/tanl.c

   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.
   3  */
   4
   5 /* @(#)s_tan.c 5.1 93/09/24 */
   6 /*
   7  * ====================================================
   8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   9  *
  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
  13  * is preserved.
  14  * ====================================================
  15  */
  16
  17 #include <config.h>
  18
  19 /* Specification.  */
  20 #include <math.h>
  21
  22 #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
  23
  24 long double
  25 tanl (long double x)
  26 {
  27   return tan (x);
  28 }
  29
  30 #else
  31
  32 /* Code based on glibc/sysdeps/ieee754/ldbl-128/s_tanl.c
  33    and           glibc/sysdeps/ieee754/ldbl-128/k_tanl.c.  */
  34
  35 /* tanl(x)
  36  * Return tangent function of x.
  37  *
  38  * kernel function:
  39  *      __kernel_tanl           ... tangent function on [-pi/4,pi/4]
  40  *      __ieee754_rem_pio2l     ... argument reduction routine
  41  *
  42  * Method.
  43  *      Let S,C and T denote the sin, cos and tan respectively on
  44  *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
  45  *      in [-pi/4 , +pi/4], and let n = k mod 4.
  46  *      We have
  47  *
  48  *          n        sin(x)      cos(x)        tan(x)
  49  *     ----------------------------------------------------------
  50  *          0          S           C             T
  51  *          1          C          -S            -1/T
  52  *          2         -S          -C             T
  53  *          3         -C           S            -1/T
  54  *     ----------------------------------------------------------
  55  *
  56  * Special cases:
  57  *      Let trig be any of sin, cos, or tan.
  58  *      trig(+-INF)  is NaN, with signals;
  59  *      trig(NaN)    is that NaN;
  60  *
  61  * Accuracy:
  62  *      TRIG(x) returns trig(x) nearly rounded
  63  */
  64
  65 # include "trigl.h"
  66
  67 /*
  68  * ====================================================
  69  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  70  *
  71  * Developed at SunPro, a Sun Microsystems, Inc. business.
  72  * Permission to use, copy, modify, and distribute this
  73  * software is freely granted, provided that this notice
  74  * is preserved.
  75  * ====================================================
  76  */
  77
  78 /*
  79   Long double expansions contributed by
  80   Stephen L. Moshier <moshier@na-net.ornl.gov>
  81 */
  82
  83 /* __kernel_tanl( x, y, k )
  84  * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
  85  * Input x is assumed to be bounded by ~pi/4 in magnitude.
  86  * Input y is the tail of x.
  87  * Input k indicates whether tan (if k=1) or
  88  * -1/tan (if k= -1) is returned.
  89  *
  90  * Algorithm
  91  *      1. Since tan(-x) = -tan(x), we need only to consider positive x.
  92  *      2. if x < 2^-57, return x with inexact if x!=0.
  93  *      3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
  94  *          on [0,0.67433].
  95  *
  96  *         Note: tan(x+y) = tan(x) + tan'(x)*y
  97  *                        ~ tan(x) + (1+x*x)*y
  98  *         Therefore, for better accuracy in computing tan(x+y), let
  99  *              r = x^3 * R(x^2)
 100  *         then
 101  *              tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
 102  *
 103  *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
 104  *              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
 105  *                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
 106  */
 107
 108
 109 static const long double
 110   pio4hi = 7.8539816339744830961566084581987569936977E-1L,
 111   pio4lo = 2.1679525325309452561992610065108379921906E-35L,
 112
 113   /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
 114      0 <= x <= 0.6743316650390625
 115      Peak relative error 8.0e-36  */
 116  TH =  3.333333333333333333333333333333333333333E-1L,
 117  T0 = -1.813014711743583437742363284336855889393E7L,
 118  T1 =  1.320767960008972224312740075083259247618E6L,
 119  T2 = -2.626775478255838182468651821863299023956E4L,
 120  T3 =  1.764573356488504935415411383687150199315E2L,
 121  T4 = -3.333267763822178690794678978979803526092E-1L,
 122
 123  U0 = -1.359761033807687578306772463253710042010E8L,
 124  U1 =  6.494370630656893175666729313065113194784E7L,
 125  U2 = -4.180787672237927475505536849168729386782E6L,
 126  U3 =  8.031643765106170040139966622980914621521E4L,
 127  U4 = -5.323131271912475695157127875560667378597E2L;
 128   /* 1.000000000000000000000000000000000000000E0 */
 129
 130
 131 static long double
 132 kernel_tanl (long double x, long double y, int iy)
 133 {
 134   long double z, r, v, w, s, u, u1;
 135   int invert = 0, sign;
 136
 137   sign = 1;
 138   if (x < 0)
 139     {
 140       x = -x;
 141       y = -y;
 142       sign = -1;
 143     }
 144
 145   if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */
 146     {
 147       if ((int) x == 0)
 148         {                       /* generate inexact */
 149           if (iy == -1 && x == 0.0)
 150             return 1.0L / fabs (x);
 151           else
 152             return (iy == 1) ? x : -1.0L / x;
 153         }
 154     }
 155   if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */
 156     {
 157       invert = 1;
 158
 159       z = pio4hi - x;
 160       w = pio4lo - y;
 161       x = z + w;
 162       y = 0.0;
 163     }
 164   z = x * x;
 165   r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
 166   v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
 167   r = r / v;
 168
 169   s = z * x;
 170   r = y + z * (s * r + y);
 171   r += TH * s;
 172   w = x + r;
 173   if (invert)
 174     {
 175       v = (long double) iy;
 176       w = (v - 2.0 * (x - (w * w / (w + v) - r)));
 177       if (sign < 0)
 178         w = -w;
 179       return w;
 180     }
 181   if (iy == 1)
 182     return w;
 183   else
 184     {                           /* if allow error up to 2 ulp,
 185                                    simply return -1.0/(x+r) here */
 186       /*  compute -1.0/(x+r) accurately */
 187       u1 = (double) w;
 188       v = r - (u1 - x);
 189       z = -1.0 / w;
 190       u = (double) z;
 191       s = 1.0 + u * u1;
 192       return u + z * (s + u * v);
 193     }
 194 }
 195
 196 long double
 197 tanl (long double x)
 198 {
 199   long double y[2], z = 0.0L;
 200   int n;
 201
 202   /* tanl(NaN) is NaN */
 203   if (isnanl (x))
 204     return x;
 205
 206   /* |x| ~< pi/4 */
 207   if (x >= -0.7853981633974483096156608458198757210492 &&
 208       x <= 0.7853981633974483096156608458198757210492)
 209     return kernel_tanl (x, z, 1);
 210
 211   /* tanl(Inf) is NaN, tanl(0) is 0 */
 212   else if (x + x == x)
 213     return x - x;               /* NaN */
 214
 215   /* argument reduction needed */
 216   else
 217     {
 218       n = ieee754_rem_pio2l (x, y);
 219       /* 1 -- n even, -1 -- n odd */
 220       return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
 221     }
 222 }
 223
 224 #endif
 225
 226 #if 0
 227 int
 228 main (void)
 229 {
 230   printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492));
 231   printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492));
 232   printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 *3));
 233   printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492 *31));
 234   printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 / 2));
 235   printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 3/2));
 236   printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 5/2));
 237 }
 238 #endif