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_macros.h>
 29  #include <libm_special.h>
 30  #include <libm/amd_funcs_internal.h>
 31  #include <libm/types.h>
 32  #include <libm/constants.h>
 33  #include <libm/typehelper.h>
 34  #include <libm/compiler.h>
 35  #include <libm/poly.h>
 36  
 37  /*
 38   * ISO-IEC-10967-2: Elementary Numerical Functions
 39   * Signature:
 40   *   float tanf(float x)
 41   *
 42   * Spec:
 43   *   tanf(n· 2π + π/4)  = 1       if n∈Z and |n· 2π + π/4|   <= big_angle_F
 44   *   tanf(n· 2π + 3π/4) = −1      if n∈Z and |n· 2π + 3π/4|  <= big angle rF
 45   *   tanf(x) = x                  if x ∈ F^(2·π) and tanf(x) != tan(x)
 46   *                                               and |x| < sqrt(epsilonF/rF)
 47   *   tanf(−x) = −tanf(x)          if x ∈ F^(2·π)
 48   */
 49  
 50  static struct {
 51      float poly_tanf[7];
 52      float halfpi, invhalfpi;
 53      float huge;
 54      float halfpi1, halfpi2, halfpi3;
 55  } tanf_data = {
 56      .huge      = 0x1.8p23,
 57      .halfpi    = 0x1.921fb6p0,
 58      .invhalfpi = 0x1.45f306p-1,
 59      .halfpi1   = -0x1.921fb6p0,
 60      .halfpi2   = 0x1.777a5cp-25,
 61      .halfpi3   = 0x1.ee59dap-50,
 62      /*
 63       * Polynomial coefficients obtained using
 64       * fpminimax algorithm from Sollya
 65       */
 66      .poly_tanf = {
 67          0x1.555566p-2,
 68          0x1.110cdp-3,
 69          0x1.baf34p-5,
 70          0x1.5bf38ep-6,
 71          0x1.663acap-7,
 72          -0x1.07c6f4p-16,
 73          0x1.21cedap-8,
 74      },
 75  };
 76  
 77  #define  ALM_TANF_HUGE      tanf_data.huge
 78  #define  ALM_TANF_HALFPI    tanf_data.halfpi
 79  #define  ALM_TANF_INVHALFPI tanf_data.invhalfpi
 80  
 81  #define  ALM_TANF_HALFPI1 tanf_data.halfpi1
 82  #define  ALM_TANF_HALFPI2 tanf_data.halfpi2
 83  #define  ALM_TANF_HALFPI3 tanf_data.halfpi3
 84  
 85  #define  ALM_TANF_SIGN_MASK32 (1U<<31)
 86  
 87  #define C1 tanf_data.poly_tanf[0]
 88  #define C2 tanf_data.poly_tanf[1]
 89  #define C3 tanf_data.poly_tanf[2]
 90  #define C4 tanf_data.poly_tanf[3]
 91  #define C5 tanf_data.poly_tanf[4]
 92  #define C6 tanf_data.poly_tanf[5]
 93  #define C7 tanf_data.poly_tanf[6]
 94  
 95  
 96  /*
 97   * Implementation Notes:
 98   *  Convert given x into the form
 99   *
100   *        |x| = N*(pi/2)+f
101   *                        where N is integer and f in [-pi/4,+pi/4]
102   *
103   *             N = round(|x|*2/pi)
104   *             f = |x| - N * pi/2
105   *
106   * tan(x) = tan(f)                when N is even
107   *        = -cot(f) = -1/tan(f) when N is odd
108   *
109   * The term tan(f) can be approximated by using a polynomial
110   */
111  
112  float
113  ALM_PROTO_FAST(tanf)(float x)
114  {
115      float      F, poly, result = 0.0;
116      uint32_t   sign, n;
117      uint32_t   ux = asuint32(x);
118  
119      if (ux == 0)
120          return 0.0f;
121  
122      sign = ux & (ALM_TANF_SIGN_MASK32);
123  
124      /* x = abs(x) */
125      x    = asfloat(ux & ~ALM_TANF_SIGN_MASK32);
126  
127      /*
128       * dn = x * (2/π)
129       * would turn to fma
130       */
131      float nn =  x * ALM_TANF_INVHALFPI + ALM_TANF_HUGE;
132  
133      /* n = (int)n */
134      n  = asuint32(nn);
135  
136      nn -= ALM_TANF_HUGE;
137  
138      /* F = xd - (n * π/2)
139       * F = n*halfpi1+x ;
140       * F = n*halfpi2+F ;
141       * F = n*halfpi3+F ;
142       */
143      F = x + nn * ALM_TANF_HALFPI1;
144      F = F + nn * ALM_TANF_HALFPI2;
145      F = F + nn * ALM_TANF_HALFPI3;
146  
147      uint32_t odd = (n & 0x1);
148  
149      // Calculate the polynomial approximation
150      //  = x+C1*x^3+C2*x^5+C3*x^7+C4*x^9+C5*x^11+C6*x^13+C7*x^15
151      poly   = POLY_EVAL_ODD_15(F, C1, C2, C3, C4, C5, C6, C7);
152  
153      result = asfloat(sign ^ asuint32(poly));
154  
155      if (odd)
156          result = -1.0/result;
157  
158      return result;
159  }