kc3-lang/gnulib/lib/tanl.c

• Show log

  Commit

• Hash : 99bd6129
  Author :
  Date : 2012-02-29T12:30:07
  

  math code: Add comments.
  
  * lib/acosl.c: Add comment about related glibc source files.
  * lib/asinl.c: Likewise.
  * lib/atanl.c: Likewise.
  * lib/expl.c: Likewise.
  * lib/logl.c: Likewise.
  * lib/sincosl.c: Likewise.
  * lib/sinl.c: Likewise.
  * lib/tanl.c: Likewise.
  * lib/trigl.c: Likewise.
  * lib/cosl.c: Likewise. Fix comments.
  

Download

• lib/tanl.c

• /* s_tanl.c -- long double version of s_tan.c.
   * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
   */
  
  /* @(#)s_tan.c 5.1 93/09/24 */
  /*
   * ====================================================
   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice
   * is preserved.
   * ====================================================
   */
  
  #include <config.h>
  
  /* Specification.  */
  #include <math.h>
  
  #if HAVE_SAME_LONG_DOUBLE_AS_DOUBLE
  
  long double
  tanl (long double x)
  {
    return tan (x);
  }
  
  #else
  
  /* Code based on glibc/sysdeps/ieee754/ldbl-128/s_tanl.c
     and           glibc/sysdeps/ieee754/ldbl-128/k_tanl.c.  */
  
  /* tanl(x)
   * Return tangent function of x.
   *
   * kernel function:
   *      __kernel_tanl           ... tangent function on [-pi/4,pi/4]
   *      __ieee754_rem_pio2l     ... argument reduction routine
   *
   * Method.
   *      Let S,C and T denote the sin, cos and tan respectively on
   *      [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
   *      in [-pi/4 , +pi/4], and let n = k mod 4.
   *      We have
   *
   *          n        sin(x)      cos(x)        tan(x)
   *     ----------------------------------------------------------
   *          0          S           C             T
   *          1          C          -S            -1/T
   *          2         -S          -C             T
   *          3         -C           S            -1/T
   *     ----------------------------------------------------------
   *
   * Special cases:
   *      Let trig be any of sin, cos, or tan.
   *      trig(+-INF)  is NaN, with signals;
   *      trig(NaN)    is that NaN;
   *
   * Accuracy:
   *      TRIG(x) returns trig(x) nearly rounded
   */
  
  # include "trigl.h"
  
  /*
   * ====================================================
   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice
   * is preserved.
   * ====================================================
   */
  
  /*
    Long double expansions contributed by
    Stephen L. Moshier <moshier@na-net.ornl.gov>
  */
  
  /* __kernel_tanl( x, y, k )
   * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854
   * Input x is assumed to be bounded by ~pi/4 in magnitude.
   * Input y is the tail of x.
   * Input k indicates whether tan (if k=1) or
   * -1/tan (if k= -1) is returned.
   *
   * Algorithm
   *      1. Since tan(-x) = -tan(x), we need only to consider positive x.
   *      2. if x < 2^-57, return x with inexact if x!=0.
   *      3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2)
   *          on [0,0.67433].
   *
   *         Note: tan(x+y) = tan(x) + tan'(x)*y
   *                        ~ tan(x) + (1+x*x)*y
   *         Therefore, for better accuracy in computing tan(x+y), let
   *              r = x^3 * R(x^2)
   *         then
   *              tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y))
   *
   *      4. For x in [0.67433,pi/4],  let y = pi/4 - x, then
   *              tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y))
   *                     = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y)))
   */
  
  
  static const long double
    pio4hi = 7.8539816339744830961566084581987569936977E-1L,
    pio4lo = 2.1679525325309452561992610065108379921906E-35L,
  
    /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2)
       0 <= x <= 0.6743316650390625
       Peak relative error 8.0e-36  */
   TH =  3.333333333333333333333333333333333333333E-1L,
   T0 = -1.813014711743583437742363284336855889393E7L,
   T1 =  1.320767960008972224312740075083259247618E6L,
   T2 = -2.626775478255838182468651821863299023956E4L,
   T3 =  1.764573356488504935415411383687150199315E2L,
   T4 = -3.333267763822178690794678978979803526092E-1L,
  
   U0 = -1.359761033807687578306772463253710042010E8L,
   U1 =  6.494370630656893175666729313065113194784E7L,
   U2 = -4.180787672237927475505536849168729386782E6L,
   U3 =  8.031643765106170040139966622980914621521E4L,
   U4 = -5.323131271912475695157127875560667378597E2L;
    /* 1.000000000000000000000000000000000000000E0 */
  
  
  static long double
  kernel_tanl (long double x, long double y, int iy)
  {
    long double z, r, v, w, s, u, u1;
    int invert = 0, sign;
  
    sign = 1;
    if (x < 0)
      {
        x = -x;
        y = -y;
        sign = -1;
      }
  
    if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */
      {
        if ((int) x == 0)
          {                       /* generate inexact */
            if (iy == -1 && x == 0.0)
              return 1.0L / fabs (x);
            else
              return (iy == 1) ? x : -1.0L / x;
          }
      }
    if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */
      {
        invert = 1;
  
        z = pio4hi - x;
        w = pio4lo - y;
        x = z + w;
        y = 0.0;
      }
    z = x * x;
    r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4)));
    v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z))));
    r = r / v;
  
    s = z * x;
    r = y + z * (s * r + y);
    r += TH * s;
    w = x + r;
    if (invert)
      {
        v = (long double) iy;
        w = (v - 2.0 * (x - (w * w / (w + v) - r)));
        if (sign < 0)
          w = -w;
        return w;
      }
    if (iy == 1)
      return w;
    else
      {                           /* if allow error up to 2 ulp,
                                     simply return -1.0/(x+r) here */
        /*  compute -1.0/(x+r) accurately */
        u1 = (double) w;
        v = r - (u1 - x);
        z = -1.0 / w;
        u = (double) z;
        s = 1.0 + u * u1;
        return u + z * (s + u * v);
      }
  }
  
  long double
  tanl (long double x)
  {
    long double y[2], z = 0.0L;
    int n;
  
    /* tanl(NaN) is NaN */
    if (isnanl (x))
      return x;
  
    /* |x| ~< pi/4 */
    if (x >= -0.7853981633974483096156608458198757210492 &&
        x <= 0.7853981633974483096156608458198757210492)
      return kernel_tanl (x, z, 1);
  
    /* tanl(Inf) is NaN, tanl(0) is 0 */
    else if (x + x == x)
      return x - x;               /* NaN */
  
    /* argument reduction needed */
    else
      {
        n = ieee754_rem_pio2l (x, y);
        /* 1 -- n even, -1 -- n odd */
        return kernel_tanl (y[0], y[1], 1 - ((n & 1) << 1));
      }
  }
  
  #endif
  
  #if 0
  int
  main (void)
  {
    printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492));
    printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492));
    printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 *3));
    printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492 *31));
    printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 / 2));
    printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 3/2));
    printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 5/2));
  }
  #endif
  

Download