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  /* sin(x) approximation valid on the interval [-pi/4,pi/4]. */
 34  static inline double sin_piby4(double x, double xx)
 35  {
 36    /* Taylor series for sin(x) is x - x^3/3! + x^5/5! - x^7/7! ...
 37                            = x * (1 - x^2/3! + x^4/5! - x^6/7! ...
 38                            = x * f(w)
 39       where w = x*x and f(w) = (1 - w/3! + w^2/5! - w^3/7! ...
 40       We use a minimax approximation of (f(w) - 1) / w
 41       because this produces an expansion in even powers of x.
 42       If xx (the tail of x) is non-zero, we add a correction
 43       term g(x,xx) = (1-x*x/2)*xx to the result, where g(x,xx)
 44       is an approximation to cos(x)*sin(xx) valid because
 45       xx is tiny relative to x.
 46    */
 47    static const double
 48      c1 = -0.166666666666666646259241729,
 49      c2 = 0.833333333333095043065222816e-2,
 50      c3 = -0.19841269836761125688538679e-3,
 51      c4 = 0.275573161037288022676895908448e-5,
 52      c5 = -0.25051132068021699772257377197e-7,
 53      c6 = 0.159181443044859136852668200e-9;
 54    double x2, x3, r;
 55    x2 = x * x;
 56    x3 = x2 * x;
 57    r = (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6))));
 58    if (xx == 0.0)
 59      return x + x3 * (c1 + x2 * r);
 60    else
 61      return x - ((x2 * (0.5 * xx - x3 * r) - xx) - x3 * c1);
 62  }
 63  
 64  /* cos(x) approximation valid on the interval [-pi/4,pi/4]. */
 65  static inline double cos_piby4(double x, double xx)
 66  {
 67    /* Taylor series for cos(x) is 1 - x^2/2! + x^4/4! - x^6/6! ...
 68                            = f(w)
 69       where w = x*x and f(w) = (1 - w/2! + w^2/4! - w^3/6! ...
 70       We use a minimax approximation of (f(w) - 1 + w/2) / (w*w)
 71       because this produces an expansion in even powers of x.
 72       If xx (the tail of x) is non-zero, we subtract a correction
 73       term g(x,xx) = x*xx to the result, where g(x,xx)
 74       is an approximation to sin(x)*sin(xx) valid because
 75       xx is tiny relative to x.
 76    */
 77    double r, x2, t;
 78    static const double
 79      c1 = 0.41666666666666665390037e-1,
 80      c2 = -0.13888888888887398280412e-2,
 81      c3 = 0.248015872987670414957399e-4,
 82      c4 = -0.275573172723441909470836e-6,
 83      c5 = 0.208761463822329611076335e-8,
 84      c6 = -0.113826398067944859590880e-10;
 85  
 86    x2 = x * x;
 87    r = 0.5 * x2;
 88    t = 1.0 - r;
 89    return t + ((((1.0 - t) - r) - x * xx) + x2 * x2 *
 90                (c1 + x2 * (c2 + x2 * (c3 + x2 * (c4 + x2 * (c5 + x2 * c6))))));
 91  }
 92  
 93  
 94  double FN_PROTOTYPE_REF(sinpi)(double x)
 95  {
 96      double r, dx, xsgn;
 97      long long ux;
 98      const double pi = 3.1415926535897932384626433832795;
 99      GET_BITS_DP64(x, ux);
100  		ux = ux & ~SIGNBIT_DP64;
101  
102      // Handle +-Inf and NaN
103      if (ux >= PINFBITPATT_DP64)
104      {
105  		PUT_BITS_DP64(QNANBITPATT_DP64,x);
106          return x;
107  	}
108  
109      // When |x| > 2^52, the result is always an integer
110      if (ux > 0x432fffffffffffffL)
111          return x * 0.0;
112  
113      dx = x > 0.0 ? x : -x;
114      GET_BITS_DP64(dx, ux);
115      // Handle small arguments
116      if (dx <= 0.25) {
117          if (ux < 0x3f20000000000000) {
118              if (ux < 0x3e40000000000000)
119                  return x * pi;
120  
121              dx = x * pi;
122              return dx - dx*dx*dx*(1.0/6.0);
123          }
124  
125          return sin_piby4(x*pi, 0.0);
126      }
127  
128      // Remaining cases
129      ux = (long)dx;
130      r = dx - ux;
131      xsgn = (x > 0.0 ? 1.0 : -1.0) * (ux & 0x1 ? -1.0 : 1.0);
132  
133      if (r == 0.0)
134          return xsgn * 0.0;
135  
136      if (r <= 0.25)
137          return xsgn * sin_piby4(r*pi, 0.0);
138  
139      if (r < 0.5)
140          return xsgn * cos_piby4((0.5 - r)*pi, 0.0);
141  
142      if (r == 0.5)
143          return xsgn;
144  
145      if (r <= 0.75)
146          return xsgn * cos_piby4((r - 0.5)*pi, 0.0);
147  
148      // r < 1
149      return xsgn * sin_piby4((1.0 - r)*pi, 0.0);
150  }
151