Attempt to represent the S^2->S^3 header reorganisation as a series of "hg rename" operations
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)
5 #ifndef BOOST_MATH_COMPLEX_ACOS_INCLUDED
6 #define BOOST_MATH_COMPLEX_ACOS_INCLUDED
8 #ifndef BOOST_MATH_COMPLEX_DETAILS_INCLUDED
9 # include <boost/math/complex/details.hpp>
11 #ifndef BOOST_MATH_LOG1P_INCLUDED
12 # include <boost/math/special_functions/log1p.hpp>
14 #include <boost/assert.hpp>
16 #ifdef BOOST_NO_STDC_NAMESPACE
17 namespace std{ using ::sqrt; using ::fabs; using ::acos; using ::asin; using ::atan; using ::atan2; }
20 namespace boost{ namespace math{
23 std::complex<T> acos(const std::complex<T>& z)
26 // This implementation is a transcription of the pseudo-code in:
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.
34 // These static constants should really be in a maths constants library:
36 static const T one = static_cast<T>(1);
37 //static const T two = static_cast<T>(2);
38 static const T half = static_cast<T>(0.5L);
39 static const T a_crossover = static_cast<T>(1.5L);
40 static const T b_crossover = static_cast<T>(0.6417L);
41 static const T s_pi = static_cast<T>(3.141592653589793238462643383279502884197L);
42 static const T half_pi = static_cast<T>(1.57079632679489661923132169163975144L);
43 static const T log_two = static_cast<T>(0.69314718055994530941723212145817657L);
44 static const T quarter_pi = static_cast<T>(0.78539816339744830961566084581987572L);
47 // Get real and imaginary parts, discard the signs as we can
48 // figure out the sign of the result later:
50 T x = std::fabs(z.real());
51 T y = std::fabs(z.imag());
53 T real, imag; // these hold our result
56 // Handle special cases specified by the C99 standard,
57 // many of these special cases aren't really needed here,
58 // but doing it this way prevents overflow/underflow arithmetic
59 // in the main body of the logic, which may trip up some machines:
61 if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity()))
63 if(y == std::numeric_limits<T>::infinity())
66 imag = std::numeric_limits<T>::infinity();
68 else if(detail::test_is_nan(y))
70 return std::complex<T>(y, -std::numeric_limits<T>::infinity());
74 // y is not infinity or nan:
76 imag = std::numeric_limits<T>::infinity();
79 else if(detail::test_is_nan(x))
81 if(y == std::numeric_limits<T>::infinity())
82 return std::complex<T>(x, (z.imag() < 0) ? std::numeric_limits<T>::infinity() : -std::numeric_limits<T>::infinity());
83 return std::complex<T>(x, x);
85 else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity()))
88 imag = std::numeric_limits<T>::infinity();
90 else if(detail::test_is_nan(y))
92 return std::complex<T>((x == 0) ? half_pi : y, y);
97 // What follows is the regular Hull et al code,
98 // begin with the special case for real numbers:
100 if((y == 0) && (x <= one))
101 return std::complex<T>((x == 0) ? half_pi : std::acos(z.real()));
103 // Figure out if our input is within the "safe area" identified by Hull et al.
104 // This would be more efficient with portable floating point exception handling;
105 // fortunately the quantities M and u identified by Hull et al (figure 3),
106 // match with the max and min methods of numeric_limits<T>.
108 T safe_max = detail::safe_max(static_cast<T>(8));
109 T safe_min = detail::safe_min(static_cast<T>(4));
114 if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min))
117 T r = std::sqrt(xp1*xp1 + yy);
118 T s = std::sqrt(xm1*xm1 + yy);
119 T a = half * (r + s);
131 real = std::atan(std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1)))/x);
135 real = std::atan((y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1))))/x);
144 am1 = half * (yy/(r + xp1) + yy/(s - xm1));
148 am1 = half * (yy/(r + xp1) + (s + xm1));
150 imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one)));
154 imag = std::log(a + std::sqrt(a*a - one));
160 // This is the Hull et al exception handling code from Fig 6 of their paper:
162 if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1)))
167 imag = y / std::sqrt(xp1*(one-x));
172 if(((std::numeric_limits<T>::max)() / xp1) > xm1)
174 // xp1 * xm1 won't overflow:
175 imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1));
179 imag = log_two + std::log(x);
183 else if(y <= safe_min)
185 // There is an assumption in Hull et al's analysis that
186 // if we get here then x == 1. This is true for all "good"
189 // E^2 > 8*sqrt(u); with:
191 // E = std::numeric_limits<T>::epsilon()
192 // u = (std::numeric_limits<T>::min)()
194 // Hull et al provide alternative code for "bad" machines
195 // but we have no way to test that here, so for now just assert
196 // on the assumption:
198 BOOST_ASSERT(x == 1);
202 else if(std::numeric_limits<T>::epsilon() * y - one >= x)
205 imag = log_two + std::log(y);
209 real = std::atan(y/x);
211 imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy);
216 T a = std::sqrt(one + y*y);
217 imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a));
223 // Finish off by working out the sign of the result:
230 return std::complex<T>(real, imag);
235 #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED