tanf.c
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 extern void __amd_remainder_piby2d2f(unsigned long long ux, double *r, int *region); 34 35 /* tan(x) approximation valid on the interval [-pi/4,pi/4]. 36 If recip is true return -1/tan(x) instead. */ 37 static inline double tanf_piby4(double x, int recip) 38 { 39 double r, t; 40 41 /* Core Remez [1,2] approximation to tan(x) on the 42 interval [0,pi/4]. */ 43 r = x*x; 44 t = x + x*r* 45 (0.385296071263995406715129e0 - 46 0.172032480471481694693109e-1 * r) / 47 (0.115588821434688393452299e+1 + 48 (-0.51396505478854532132342e0 + 49 0.1844239256901656082986661e-1 * r) * r); 50 51 if (recip) 52 return -1.0 / t; 53 else 54 return t; 55 } 56 57 #ifdef WINDOWS 58 #pragma function(tanf) 59 #endif 60 61 float FN_PROTOTYPE_REF(tanf)(float x) 62 { 63 double r, dx; 64 int region, xneg; 65 unsigned int fux; 66 unsigned long long ux, ax; 67 GET_BITS_SP32(x, fux); 68 69 if ((fux & EXPBITS_SP32) == EXPBITS_SP32) 70 { 71 /* x is either NaN or infinity */ 72 if (fux & MANTBITS_SP32) 73 { 74 /* x is NaN */ 75 #ifdef WINDOWS 76 return __amd_handle_errorf("tanf", __amd_tan, fux | QNAN_MASK_32, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F, 1); 77 #else 78 if (fux & QNAN_MASK_32) 79 return __amd_handle_errorf("tanf", __amd_tan, fux | QNAN_MASK_32, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0F, 1); 80 else 81 return __amd_handle_errorf("tanf", __amd_tan, fux | QNAN_MASK_32, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0F, 1); 82 #endif 83 } 84 else 85 { 86 /* x is infinity. Return a NaN */ 87 return __amd_handle_errorf("tanf", __amd_tan, INDEFBITPATT_SP32, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0F, 1); 88 } 89 } 90 91 dx = (double)x; 92 GET_BITS_DP64(dx, ux); 93 ax = (ux & ~SIGNBIT_DP64); 94 if (ax <= 0x3fe921fb54442d18LL) /* abs(x) <= pi/4 */ 95 { 96 if (ax < 0x3f80000000000000LL) /* abs(x) < 2.0^(-7) */ 97 { 98 if (ax < 0x3f20000000000000LL) /* abs(x) < 2.0^(-13) */ 99 { 100 if (ax == 0x0000000000000000LL) 101 return x; 102 else 103 #ifdef WINDOWS 104 return x; //valf_with_flags(x, AMD_F_INEXACT); 105 #else 106 return __amd_handle_errorf("tanf", __amd_tan, fux, _UNDERFLOW, AMD_F_UNDERFLOW|AMD_F_INEXACT, ERANGE, x, 0.0F, 1); 107 #endif 108 } 109 else 110 return (float)(dx + dx*dx*dx*0.333333333333333333); 111 } 112 else 113 return (float)tanf_piby4(x, 0); 114 } 115 116 117 xneg = (int)(ux >> 63); 118 119 if (xneg) 120 dx = -dx; 121 122 if (dx < 5.0e5) 123 { 124 /* For these size arguments we can just carefully subtract the 125 appropriate multiple of pi/2, using extra precision where 126 dx is close to an exact multiple of pi/2 */ 127 static const double 128 twobypi = 6.36619772367581382433e-01, /* 0x3fe45f306dc9c883 */ 129 piby2_1 = 1.57079632673412561417e+00, /* 0x3ff921fb54400000 */ 130 piby2_1tail = 6.07710050650619224932e-11, /* 0x3dd0b4611a626331 */ 131 piby2_2 = 6.07710050630396597660e-11, /* 0x3dd0b4611a600000 */ 132 piby2_2tail = 2.02226624879595063154e-21, /* 0x3ba3198a2e037073 */ 133 piby2_3 = 2.02226624871116645580e-21, /* 0x3ba3198a2e000000 */ 134 piby2_3tail = 8.47842766036889956997e-32; /* 0x397b839a252049c1 */ 135 double t, rhead, rtail; 136 int npi2; 137 unsigned long long uy, xexp, expdiff; 138 xexp = ax >> EXPSHIFTBITS_DP64; 139 /* How many pi/2 is dx a multiple of? */ 140 if (ax <= 0x400f6a7a2955385eLL) /* 5pi/4 */ 141 { 142 if (ax <= 0x4002d97c7f3321d2LL) /* 3pi/4 */ 143 npi2 = 1; 144 else 145 npi2 = 2; 146 } 147 else if (ax <= 0x401c463abeccb2bbLL) /* 9pi/4 */ 148 { 149 if (ax <= 0x4015fdbbe9bba775LL) /* 7pi/4 */ 150 npi2 = 3; 151 else 152 npi2 = 4; 153 } 154 else 155 npi2 = (int)(dx * twobypi + 0.5); 156 /* Subtract the multiple from dx to get an extra-precision remainder */ 157 rhead = dx - npi2 * piby2_1; 158 rtail = npi2 * piby2_1tail; 159 GET_BITS_DP64(rhead, uy); 160 expdiff = xexp - ((uy & EXPBITS_DP64) >> EXPSHIFTBITS_DP64); 161 if (expdiff > 15) 162 { 163 /* The remainder is pretty small compared with dx, which 164 implies that dx is a near multiple of pi/2 165 (dx matches the multiple to at least 15 bits) */ 166 t = rhead; 167 rtail = npi2 * piby2_2; 168 rhead = t - rtail; 169 rtail = npi2 * piby2_2tail - ((t - rhead) - rtail); 170 if (expdiff > 48) 171 { 172 /* dx matches a pi/2 multiple to at least 48 bits */ 173 t = rhead; 174 rtail = npi2 * piby2_3; 175 rhead = t - rtail; 176 rtail = npi2 * piby2_3tail - ((t - rhead) - rtail); 177 } 178 } 179 r = rhead - rtail; 180 region = npi2 & 3; 181 } 182 else 183 { 184 /* Reduce x into range [-pi/4,pi/4] */ 185 __amd_remainder_piby2d2f(ax, &r, ®ion); 186 } 187 188 if (xneg) 189 return (float)-tanf_piby4(r, region & 1); 190 else 191 return (float)tanf_piby4(r, region & 1); 192 } 193 194