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>
acos(const std::complex<T> & z)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
246