1  /*
  2   * Copyright (C) 2008-2020 Advanced Micro Devices, Inc. All rights reserved.
  3   *
  4   * Redistribution and use in source and binary forms, with or without modification,
  5   * are permitted provided that the following conditions are met:
  6   * 1. Redistributions of source code must retain the above copyright notice,
  7   *    this list of conditions and the following disclaimer.
  8   * 2. Redistributions in binary form must reproduce the above copyright notice,
  9   *    this list of conditions and the following disclaimer in the documentation
 10   *    and/or other materials provided with the distribution.
 11   * 3. Neither the name of the copyright holder nor the names of its contributors
 12   *    may be used to endorse or promote products derived from this software without
 13   *    specific prior written permission.
 14   *
 15   * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 16   * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 17   * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
 18   * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
 19   * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 20   * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
 21   * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 22   * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 23   * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 24   * POSSIBILITY OF SUCH DAMAGE.
 25   *
 26   */
 27  
 28  #include "libm_amd.h"
 29  #include "libm_util_amd.h"
 30  #include "libm_inlines_amd.h"
 31  #include "libm_special.h"
 32  
 33  
 34  // tan(x + xx) approximation valid on the interval [-pi/4,pi/4].
 35  // If recip is true return -1/tan(x + xx) instead.
 36  static inline double tan_piby4(double x, double xx, int recip)
 37  {
 38    double r, t1, t2, xl;
 39    int transform = 0;
 40    static const double
 41       piby4_lead = 7.85398163397448278999e-01, /* 0x3fe921fb54442d18 */
 42       piby4_tail = 3.06161699786838240164e-17; /* 0x3c81a62633145c06 */
 43  
 44    /* In order to maintain relative precision transform using the identity:
 45       tan(pi/4-x) = (1-tan(x))/(1+tan(x)) for arguments close to pi/4.
 46       Similarly use tan(x-pi/4) = (tan(x)-1)/(tan(x)+1) close to -pi/4. */
 47  
 48    if (x > 0.68)
 49      {
 50        transform = 1;
 51        x = piby4_lead - x;
 52        xl = piby4_tail - xx;
 53        x += xl;
 54        xx = 0.0;
 55      }
 56    else if (x < -0.68)
 57      {
 58        transform = -1;
 59        x = piby4_lead + x;
 60        xl = piby4_tail + xx;
 61        x += xl;
 62        xx = 0.0;
 63      }
 64  
 65    /* Core Remez [2,3] approximation to tan(x+xx) on the
 66       interval [0,0.68]. */
 67  
 68    r = x*x + 2.0 * x * xx;
 69    t1 = x;
 70    t2 = xx + x*r*
 71      (0.372379159759792203640806338901e0 +
 72       (-0.229345080057565662883358588111e-1 +
 73        0.224044448537022097264602535574e-3*r)*r)/
 74      (0.111713747927937668539901657944e1 +
 75       (-0.515658515729031149329237816945e0 +
 76        (0.260656620398645407524064091208e-1 -
 77         0.232371494088563558304549252913e-3*r)*r)*r);
 78  
 79    /* Reconstruct tan(x) in the transformed case. */
 80  
 81    if (transform)
 82      {
 83        double t;
 84        t = t1 + t2;
 85        if (recip)
 86           return transform*(2*t/(t-1) - 1.0);
 87        else
 88           return transform*(1.0 - 2*t/(1+t));
 89      }
 90  
 91    if (recip)
 92      {
 93        /* Compute -1.0/(t1 + t2) accurately */
 94        double trec, trec_top, z1, z2, t;
 95        unsigned long long u;
 96        t = t1 + t2;
 97        GET_BITS_DP64(t, u);
 98        u &= 0xffffffff00000000;
 99        PUT_BITS_DP64(u, z1);
100        z2 = t2 - (z1 - t1);
101        trec = -1.0 / t;
102        GET_BITS_DP64(trec, u);
103        u &= 0xffffffff00000000;
104        PUT_BITS_DP64(u, trec_top);
105        return trec_top + trec * ((1.0 + trec_top * z1) + trec_top * z2);
106  
107      }
108    else
109      return t1 + t2;
110  }
111  
112  
113  
114  double FN_PROTOTYPE_REF(tanpi)(double x)
115  {
116      double r, dx, xsgn, xodd;
117      long long ux;
118      const double pi = 3.1415926535897932384626433832795;
119  
120      xsgn = x > 0.0 ? 1.0 : -1.0;
121      GET_BITS_DP64(x, ux);
122      ux = ux & ~SIGNBIT_DP64;
123  
124      // Handle +-Inf and NaN
125      if (ux >= PINFBITPATT_DP64)
126      {
127  		PUT_BITS_DP64(INDEFBITPATT_DP64,x);
128  		return x;
129  	}
130  
131      // If x >= 2^53, it is an even integer
132      if (ux > 0x433fffffffffffffL)
133  	return xsgn * 0.0;
134  
135      // If x >= 2^52, it is an integer
136      // Given the previous test, we also know that the last bits of ux are the
137      // same as the last bits of (long)ux
138      if (ux > 0x432fffffffffffffL)
139  	return xsgn * (ux & 0x1L ? -0.0 : 0.0);
140  
141      dx = x * xsgn;
142  
143      GET_BITS_DP64(dx, ux);
144      // Take care of small arguments
145      if (dx <= 0.25) {
146          if (ux < 0x3f10000000000000) {
147  	    if (ux < 0x3e40000000000000)
148  		return x * pi;
149  	    dx = x * pi;
150  	    return dx + dx*dx*dx*(1.0 / 3.0);
151  	}
152  	return tan_piby4(x*pi, 0.0, 0);
153      }
154  
155      ux = (long)dx;
156      r = dx - (double)ux;
157      xodd = xsgn * (ux & 0x1 ? -1.0 : 1.0);
158  
159      if (r == 0.0)
160  	return xodd * 0.0;
161  
162      if (r <= 0.25)
163  	return xsgn * tan_piby4(r*pi, 0.0, 0);
164  
165      if (r < 0.5)
166  	return -xsgn * tan_piby4((0.5 - r)*pi, 0.0, 1);
167  
168      if (r == 0.5)
169      {
170  		PUT_BITS_DP64(PINFBITPATT_DP64,x);
171  		return xodd * x;
172  	}
173  
174      if (r <= 0.75)
175  	return xsgn * tan_piby4((r - 0.5)*pi, 0.0, 1);
176  
177      // r < 1
178      return -xsgn * tan_piby4((1.0 - r)*pi, 0.0, 0);
179  }
180  
181  
182  
183