1
2 /*******************************************************************************
3 MIT License
4 -----------
5
6 Copyright (c) 2002-2019 Advanced Micro Devices, Inc.
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this Software and associated documentaon files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
24 THE SOFTWARE.
25 *******************************************************************************/
26
27 #include "libm.h"
28 #include "libm_util.h"
29
30 #define USE_VAL_WITH_FLAGS
31 #define USE_NAN_WITH_FLAGS
32 #define USE_HANDLE_ERROR
33 #include "libm_inlines.h"
34 #undef USE_NAN_WITH_FLAGS
35 #undef USE_VAL_WITH_FLAGS
36 #undef USE_HANDLE_ERROR
37
38 #include "libm_errno.h"
39
40 #ifdef _MSC_VER
41 #pragma function(acos)
42 #endif
43
FN_PROTOTYPE(acos)44 double FN_PROTOTYPE(acos)(double x)
45 {
46 /* Computes arccos(x).
47 The argument is first reduced by noting that arccos(x)
48 is invalid for abs(x) > 1. For denormal and small
49 arguments arccos(x) = pi/2 to machine accuracy.
50 Remaining argument ranges are handled as follows.
51 For abs(x) <= 0.5 use
52 arccos(x) = pi/2 - arcsin(x)
53 = pi/2 - (x + x^3*R(x^2))
54 where R(x^2) is a rational minimax approximation to
55 (arcsin(x) - x)/x^3.
56 For abs(x) > 0.5 exploit the identity:
57 arccos(x) = pi - 2*arcsin(sqrt(1-x)/2)
58 together with the above rational approximation, and
59 reconstruct the terms carefully.
60 */
61
62 /* Some constants and split constants. */
63
64 static const double
65 pi = 3.1415926535897933e+00, /* 0x400921fb54442d18 */
66 piby2 = 1.5707963267948965580e+00, /* 0x3ff921fb54442d18 */
67 piby2_head = 1.5707963267948965580e+00, /* 0x3ff921fb54442d18 */
68 piby2_tail = 6.12323399573676603587e-17; /* 0x3c91a62633145c07 */
69
70 double u, y, s=0.0, r;
71 int xexp, xnan, transform=0;
72
73 unsigned long long ux, aux, xneg;
74 GET_BITS_DP64(x, ux);
75 aux = ux & ~SIGNBIT_DP64;
76 xneg = (ux & SIGNBIT_DP64);
77 xnan = (aux > PINFBITPATT_DP64);
78 xexp = (int)((ux & EXPBITS_DP64) >> EXPSHIFTBITS_DP64) - EXPBIAS_DP64;
79
80 /* Special cases */
81
82 if (xnan)
83 {
84 return _handle_error("acos", OP_ACOS, ux|0x0008000000000000, _DOMAIN,
85 0, EDOM, x, 0.0, 1);
86 }
87 else if (xexp < -56)
88 { /* y small enough that arccos(x) = pi/2 */
89 return val_with_flags(piby2, AMD_F_INEXACT);
90 }
91 else if (xexp >= 0)
92 { /* abs(x) >= 1.0 */
93 if (x == 1.0)
94 return 0.0;
95 else if (x == -1.0)
96 return val_with_flags(pi, AMD_F_INEXACT);
97 else
98 return _handle_error("acos", OP_ACOS, INDEFBITPATT_DP64, _DOMAIN,
99 AMD_F_INVALID, EDOM, x, 0.0, 1);
100 }
101
102 if (xneg) y = -x;
103 else y = x;
104
105 transform = (xexp >= -1); /* abs(x) >= 0.5 */
106
107 if (transform)
108 { /* Transform y into the range [0,0.5) */
109 r = 0.5*(1.0 - y);
110 /* VC++ intrinsic call */
111 _mm_store_sd(&s, _mm_sqrt_sd(_mm_setzero_pd(), _mm_load_sd(&r)));
112 y = s;
113 }
114 else
115 r = y*y;
116
117 /* Use a rational approximation for [0.0, 0.5] */
118
119 u = r*(0.227485835556935010735943483075 +
120 (-0.445017216867635649900123110649 +
121 (0.275558175256937652532686256258 +
122 (-0.0549989809235685841612020091328 +
123 (0.00109242697235074662306043804220 +
124 0.0000482901920344786991880522822991*r)*r)*r)*r)*r)/
125 (1.36491501334161032038194214209 +
126 (-3.28431505720958658909889444194 +
127 (2.76568859157270989520376345954 +
128 (-0.943639137032492685763471240072 +
129 0.105869422087204370341222318533*r)*r)*r)*r);
130
131 if (transform)
132 { /* Reconstruct acos carefully in transformed region */
133 if (xneg) return pi - 2.0*(s+(y*u - piby2_tail));
134 else
135 {
136 double c, s1;
137 unsigned long long us;
138 GET_BITS_DP64(s, us);
139 PUT_BITS_DP64(0xffffffff00000000 & us, s1);
140 c = (r-s1*s1)/(s+s1);
141 return 2.0*s1 + (2.0*c+2.0*y*u);
142 }
143 }
144 else
145 return piby2_head - (x - (piby2_tail - x*u));
146 }
147