]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2005. |
2 | // Distributed under the Boost Software License, Version 1.0. (See accompanying | |
3 | // file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
4 | ||
5 | #ifndef BOOST_MATH_COMPLEX_ACOS_INCLUDED | |
6 | #define BOOST_MATH_COMPLEX_ACOS_INCLUDED | |
7 | ||
8 | #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED | |
9 | # include <boost/math/complex/details.hpp> | |
10 | #endif | |
11 | #ifndef BOOST_MATH_LOG1P_INCLUDED | |
12 | # include <boost/math/special_functions/log1p.hpp> | |
13 | #endif | |
14 | #include <boost/assert.hpp> | |
15 | ||
16 | #ifdef BOOST_NO_STDC_NAMESPACE | |
17 | namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; } | |
18 | #endif | |
19 | ||
20 | namespace boost{ namespace math{ | |
21 | ||
22 | template<class T> | |
23 | std::complex<T> acos(const std::complex<T>& z) | |
24 | { | |
25 | // | |
26 | // This implementation is a transcription of the pseudo-code in: | |
27 | // | |
28 | // "Implementing the Complex Arcsine and Arccosine Functions using Exception Handling." | |
29 | // T E Hull, Thomas F Fairgrieve and Ping Tak Peter Tang. | |
30 | // ACM Transactions on Mathematical Software, Vol 23, No 3, Sept 1997. | |
31 | // | |
32 | ||
33 | // | |
34 | // These static constants should really be in a maths constants library, | |
35 | // note that we have tweaked a_crossover as per: https://svn.boost.org/trac/boost/ticket/7290 | |
36 | // | |
37 | static const T one = static_cast<T>(1); | |
38 | //static const T two = static_cast<T>(2); | |
39 | static const T half = static_cast<T>(0.5L); | |
40 | static const T a_crossover = static_cast<T>(10); | |
41 | static const T b_crossover = static_cast<T>(0.6417L); | |
42 | static const T s_pi = boost::math::constants::pi<T>(); | |
43 | static const T half_pi = s_pi / 2; | |
44 | static const T log_two = boost::math::constants::ln_two<T>(); | |
45 | static const T quarter_pi = s_pi / 4; | |
46 | ||
47 | #ifdef BOOST_MSVC | |
48 | #pragma warning(push) | |
49 | #pragma warning(disable:4127) | |
50 | #endif | |
51 | // | |
52 | // Get real and imaginary parts, discard the signs as we can | |
53 | // figure out the sign of the result later: | |
54 | // | |
55 | T x = std::fabs(z.real()); | |
56 | T y = std::fabs(z.imag()); | |
57 | ||
58 | T real, imag; // these hold our result | |
59 | ||
60 | // | |
61 | // Handle special cases specified by the C99 standard, | |
62 | // many of these special cases aren't really needed here, | |
63 | // but doing it this way prevents overflow/underflow arithmetic | |
64 | // in the main body of the logic, which may trip up some machines: | |
65 | // | |
66 | if((boost::math::isinf)(x)) | |
67 | { | |
68 | if((boost::math::isinf)(y)) | |
69 | { | |
70 | real = quarter_pi; | |
71 | imag = std::numeric_limits<T>::infinity(); | |
72 | } | |
73 | else if((boost::math::isnan)(y)) | |
74 | { | |
75 | return std::complex<T>(y, -std::numeric_limits<T>::infinity()); | |
76 | } | |
77 | else | |
78 | { | |
79 | // y is not infinity or nan: | |
80 | real = 0; | |
81 | imag = std::numeric_limits<T>::infinity(); | |
82 | } | |
83 | } | |
84 | else if((boost::math::isnan)(x)) | |
85 | { | |
86 | if((boost::math::isinf)(y)) | |
87 | return std::complex<T>(x, ((boost::math::signbit)(z.imag())) ? std::numeric_limits<T>::infinity() : -std::numeric_limits<T>::infinity()); | |
88 | return std::complex<T>(x, x); | |
89 | } | |
90 | else if((boost::math::isinf)(y)) | |
91 | { | |
92 | real = half_pi; | |
93 | imag = std::numeric_limits<T>::infinity(); | |
94 | } | |
95 | else if((boost::math::isnan)(y)) | |
96 | { | |
97 | return std::complex<T>((x == 0) ? half_pi : y, y); | |
98 | } | |
99 | else | |
100 | { | |
101 | // | |
102 | // What follows is the regular Hull et al code, | |
103 | // begin with the special case for real numbers: | |
104 | // | |
105 | if((y == 0) && (x <= one)) | |
106 | return std::complex<T>((x == 0) ? half_pi : std::acos(z.real()), (boost::math::changesign)(z.imag())); | |
107 | // | |
108 | // Figure out if our input is within the "safe area" identified by Hull et al. | |
109 | // This would be more efficient with portable floating point exception handling; | |
110 | // fortunately the quantities M and u identified by Hull et al (figure 3), | |
111 | // match with the max and min methods of numeric_limits<T>. | |
112 | // | |
113 | T safe_max = detail::safe_max(static_cast<T>(8)); | |
114 | T safe_min = detail::safe_min(static_cast<T>(4)); | |
115 | ||
116 | T xp1 = one + x; | |
117 | T xm1 = x - one; | |
118 | ||
119 | if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min)) | |
120 | { | |
121 | T yy = y * y; | |
122 | T r = std::sqrt(xp1*xp1 + yy); | |
123 | T s = std::sqrt(xm1*xm1 + yy); | |
124 | T a = half * (r + s); | |
125 | T b = x / a; | |
126 | ||
127 | if(b <= b_crossover) | |
128 | { | |
129 | real = std::acos(b); | |
130 | } | |
131 | else | |
132 | { | |
133 | T apx = a + x; | |
134 | if(x <= one) | |
135 | { | |
136 | real = std::atan(std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1)))/x); | |
137 | } | |
138 | else | |
139 | { | |
140 | real = std::atan((y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1))))/x); | |
141 | } | |
142 | } | |
143 | ||
144 | if(a <= a_crossover) | |
145 | { | |
146 | T am1; | |
147 | if(x < one) | |
148 | { | |
149 | am1 = half * (yy/(r + xp1) + yy/(s - xm1)); | |
150 | } | |
151 | else | |
152 | { | |
153 | am1 = half * (yy/(r + xp1) + (s + xm1)); | |
154 | } | |
155 | imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one))); | |
156 | } | |
157 | else | |
158 | { | |
159 | imag = std::log(a + std::sqrt(a*a - one)); | |
160 | } | |
161 | } | |
162 | else | |
163 | { | |
164 | // | |
165 | // This is the Hull et al exception handling code from Fig 6 of their paper: | |
166 | // | |
167 | if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1))) | |
168 | { | |
169 | if(x < one) | |
170 | { | |
171 | real = std::acos(x); | |
172 | imag = y / std::sqrt(xp1*(one-x)); | |
173 | } | |
174 | else | |
175 | { | |
176 | // This deviates from Hull et al's paper as per https://svn.boost.org/trac/boost/ticket/7290 | |
177 | if(((std::numeric_limits<T>::max)() / xp1) > xm1) | |
178 | { | |
179 | // xp1 * xm1 won't overflow: | |
180 | real = y / std::sqrt(xm1*xp1); | |
181 | imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1)); | |
182 | } | |
183 | else | |
184 | { | |
185 | real = y / x; | |
186 | imag = log_two + std::log(x); | |
187 | } | |
188 | } | |
189 | } | |
190 | else if(y <= safe_min) | |
191 | { | |
192 | // There is an assumption in Hull et al's analysis that | |
193 | // if we get here then x == 1. This is true for all "good" | |
194 | // machines where : | |
195 | // | |
196 | // E^2 > 8*sqrt(u); with: | |
197 | // | |
198 | // E = std::numeric_limits<T>::epsilon() | |
199 | // u = (std::numeric_limits<T>::min)() | |
200 | // | |
201 | // Hull et al provide alternative code for "bad" machines | |
202 | // but we have no way to test that here, so for now just assert | |
203 | // on the assumption: | |
204 | // | |
205 | BOOST_ASSERT(x == 1); | |
206 | real = std::sqrt(y); | |
207 | imag = std::sqrt(y); | |
208 | } | |
209 | else if(std::numeric_limits<T>::epsilon() * y - one >= x) | |
210 | { | |
211 | real = half_pi; | |
212 | imag = log_two + std::log(y); | |
213 | } | |
214 | else if(x > one) | |
215 | { | |
216 | real = std::atan(y/x); | |
217 | T xoy = x/y; | |
218 | imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy); | |
219 | } | |
220 | else | |
221 | { | |
222 | real = half_pi; | |
223 | T a = std::sqrt(one + y*y); | |
224 | imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a)); | |
225 | } | |
226 | } | |
227 | } | |
228 | ||
229 | // | |
230 | // Finish off by working out the sign of the result: | |
231 | // | |
232 | if((boost::math::signbit)(z.real())) | |
233 | real = s_pi - real; | |
234 | if(!(boost::math::signbit)(z.imag())) | |
235 | imag = (boost::math::changesign)(imag); | |
236 | ||
237 | return std::complex<T>(real, imag); | |
238 | #ifdef BOOST_MSVC | |
239 | #pragma warning(pop) | |
240 | #endif | |
241 | } | |
242 | ||
243 | } } // namespaces | |
244 | ||
245 | #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED |