asin.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_special.h" 31 32 double FN_PROTOTYPE_REF(asin)(double x) 33 { 34 /* Computes arcsin(x). 35 The argument is first reduced by noting that arcsin(x) 36 is invalid for abs(x) > 1 and arcsin(-x) = -arcsin(x). 37 For denormal and small arguments arcsin(x) = x to machine 38 accuracy. Remaining argument ranges are handled as follows. 39 For abs(x) <= 0.5 use 40 arcsin(x) = x + x^3*R(x^2) 41 where R(x^2) is a rational minimax approximation to 42 (arcsin(x) - x)/x^3. 43 For abs(x) > 0.5 exploit the identity: 44 arcsin(x) = pi/2 - 2*arcsin(sqrt(1-x)/2) 45 together with the above rational approximation, and 46 reconstruct the terms carefully. 47 */ 48 49 /* Some constants and split constants. */ 50 51 static const double 52 piby2_tail = 6.1232339957367660e-17, /* 0x3c91a62633145c07 */ 53 hpiby2_head = 7.8539816339744831e-01, /* 0x3fe921fb54442d18 */ 54 piby2 = 1.5707963267948965e+00; /* 0x3ff921fb54442d18 */ 55 double u, v, y, s=0.0, r; 56 int xexp, xnan, transform=0; 57 58 unsigned long long ux, aux, xneg; 59 GET_BITS_DP64(x, ux); 60 aux = ux & ~SIGNBIT_DP64; 61 xneg = (ux & SIGNBIT_DP64); 62 xnan = (aux > PINFBITPATT_DP64); 63 xexp = (int)((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64) - EXPBIAS_DP64; 64 65 /* Special cases */ 66 67 if (xnan) 68 { 69 #ifdef WINDOWS 70 return __amd_handle_error("asin", __amd_asin, ux|0x0008000000000000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0, 1); 71 #else 72 //return x + x; /* With invalid if it's a signalling NaN */ 73 if (ux & QNAN_MASK_64) 74 return __amd_handle_error("asin", __amd_asin, ux|0x0008000000000000, _DOMAIN, AMD_F_NONE, EDOM, x, 0.0, 1); 75 else 76 return __amd_handle_error("asin", __amd_asin, ux|0x0008000000000000, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0, 1); 77 #endif 78 } 79 else if (xexp < -28) 80 { /* y small enough that arcsin(x) = x */ 81 #ifdef WINDOWS 82 return x; //val_with_flags(x, AMD_F_INEXACT); 83 #else 84 if ((ux == SIGNBIT_DP64) || (ux == 0x0)) 85 return x; 86 else 87 return __amd_handle_error("asin", __amd_asin, ux,_UNDERFLOW, AMD_F_UNDERFLOW | AMD_F_INEXACT, ERANGE , x, 0.0, 1); 88 #endif 89 } 90 else if (xexp >= 0) 91 { /* abs(x) >= 1.0 */ 92 if (x == 1.0) 93 return piby2; //val_with_flags(piby2, AMD_F_INEXACT); 94 else if (x == -1.0) 95 return -piby2;//val_with_flags(-piby2, AMD_F_INEXACT); 96 else 97 #ifdef WINDOWS 98 return __amd_handle_error("asin", __amd_asin, INDEFBITPATT_DP64, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0, 1); 99 #else 100 //return retval_errno_edom(x); 101 return __amd_handle_error("asin", __amd_asin, INDEFBITPATT_DP64, _DOMAIN, AMD_F_INVALID, EDOM, x, 0.0, 1); 102 #endif 103 } 104 105 if (xneg) y = -x; 106 else y = x; 107 108 transform = (xexp >= -1); /* abs(x) >= 0.5 */ 109 110 if (transform) 111 { /* Transform y into the range [0,0.5) */ 112 r = 0.5*(1.0 - y); 113 #ifdef WINDOWS 114 /* VC++ intrinsic call */ 115 _mm_store_sd(&s, _mm_sqrt_sd(_mm_setzero_pd(), _mm_load_sd(&r))); 116 #else 117 /* Hammer sqrt instruction */ 118 asm volatile ("sqrtsd %1, %0" : "=x" (s) : "x" (r)); 119 #endif 120 y = s; 121 } 122 else 123 r = y*y; 124 125 /* Use a rational approximation for [0.0, 0.5] */ 126 127 u = r*(0.227485835556935010735943483075 + 128 (-0.445017216867635649900123110649 + 129 (0.275558175256937652532686256258 + 130 (-0.0549989809235685841612020091328 + 131 (0.00109242697235074662306043804220 + 132 0.0000482901920344786991880522822991*r)*r)*r)*r)*r)/ 133 (1.36491501334161032038194214209 + 134 (-3.28431505720958658909889444194 + 135 (2.76568859157270989520376345954 + 136 (-0.943639137032492685763471240072 + 137 0.105869422087204370341222318533*r)*r)*r)*r); 138 139 if (transform) 140 { /* Reconstruct asin carefully in transformed region */ 141 { 142 double c, s1, p, q; 143 unsigned long long us; 144 GET_BITS_DP64(s, us); 145 PUT_BITS_DP64(0xffffffff00000000 & us, s1); 146 c = (r-s1*s1)/(s+s1); 147 p = 2.0*s*u - (piby2_tail-2.0*c); 148 q = hpiby2_head - 2.0*s1; 149 v = hpiby2_head - (p-q); 150 } 151 } 152 else 153 { 154 #ifdef WINDOWS 155 /* Use a temporary variable to prevent VC++ rearranging 156 y + y*u 157 into 158 y * (1 + u) 159 and getting an incorrectly rounded result */ 160 double tmp; 161 tmp = y * u; 162 v = y + tmp; 163 #else 164 v = y + y*u; 165 #endif 166 } 167 168 if (xneg) return -v; 169 else return v; 170 } 171 172