[29] | 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 | // |
---|
| 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); |
---|
| 45 | |
---|
| 46 | // |
---|
| 47 | // Get real and imaginary parts, discard the signs as we can |
---|
| 48 | // figure out the sign of the result later: |
---|
| 49 | // |
---|
| 50 | T x = std::fabs(z.real()); |
---|
| 51 | T y = std::fabs(z.imag()); |
---|
| 52 | |
---|
| 53 | T real, imag; // these hold our result |
---|
| 54 | |
---|
| 55 | // |
---|
| 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: |
---|
| 60 | // |
---|
| 61 | if(std::numeric_limits<T>::has_infinity && (x == std::numeric_limits<T>::infinity())) |
---|
| 62 | { |
---|
| 63 | if(y == std::numeric_limits<T>::infinity()) |
---|
| 64 | { |
---|
| 65 | real = quarter_pi; |
---|
| 66 | imag = std::numeric_limits<T>::infinity(); |
---|
| 67 | } |
---|
| 68 | else if(detail::test_is_nan(y)) |
---|
| 69 | { |
---|
| 70 | return std::complex<T>(y, -std::numeric_limits<T>::infinity()); |
---|
| 71 | } |
---|
| 72 | else |
---|
| 73 | { |
---|
| 74 | // y is not infinity or nan: |
---|
| 75 | real = 0; |
---|
| 76 | imag = std::numeric_limits<T>::infinity(); |
---|
| 77 | } |
---|
| 78 | } |
---|
| 79 | else if(detail::test_is_nan(x)) |
---|
| 80 | { |
---|
| 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); |
---|
| 84 | } |
---|
| 85 | else if(std::numeric_limits<T>::has_infinity && (y == std::numeric_limits<T>::infinity())) |
---|
| 86 | { |
---|
| 87 | real = half_pi; |
---|
| 88 | imag = std::numeric_limits<T>::infinity(); |
---|
| 89 | } |
---|
| 90 | else if(detail::test_is_nan(y)) |
---|
| 91 | { |
---|
| 92 | return std::complex<T>((x == 0) ? half_pi : y, y); |
---|
| 93 | } |
---|
| 94 | else |
---|
| 95 | { |
---|
| 96 | // |
---|
| 97 | // What follows is the regular Hull et al code, |
---|
| 98 | // begin with the special case for real numbers: |
---|
| 99 | // |
---|
| 100 | if((y == 0) && (x <= one)) |
---|
| 101 | return std::complex<T>((x == 0) ? half_pi : std::acos(z.real())); |
---|
| 102 | // |
---|
| 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>. |
---|
| 107 | // |
---|
| 108 | T safe_max = detail::safe_max(static_cast<T>(8)); |
---|
| 109 | T safe_min = detail::safe_min(static_cast<T>(4)); |
---|
| 110 | |
---|
| 111 | T xp1 = one + x; |
---|
| 112 | T xm1 = x - one; |
---|
| 113 | |
---|
| 114 | if((x < safe_max) && (x > safe_min) && (y < safe_max) && (y > safe_min)) |
---|
| 115 | { |
---|
| 116 | T yy = y * y; |
---|
| 117 | T r = std::sqrt(xp1*xp1 + yy); |
---|
| 118 | T s = std::sqrt(xm1*xm1 + yy); |
---|
| 119 | T a = half * (r + s); |
---|
| 120 | T b = x / a; |
---|
| 121 | |
---|
| 122 | if(b <= b_crossover) |
---|
| 123 | { |
---|
| 124 | real = std::acos(b); |
---|
| 125 | } |
---|
| 126 | else |
---|
| 127 | { |
---|
| 128 | T apx = a + x; |
---|
| 129 | if(x <= one) |
---|
| 130 | { |
---|
| 131 | real = std::atan(std::sqrt(half * apx * (yy /(r + xp1) + (s-xm1)))/x); |
---|
| 132 | } |
---|
| 133 | else |
---|
| 134 | { |
---|
| 135 | real = std::atan((y * std::sqrt(half * (apx/(r + xp1) + apx/(s+xm1))))/x); |
---|
| 136 | } |
---|
| 137 | } |
---|
| 138 | |
---|
| 139 | if(a <= a_crossover) |
---|
| 140 | { |
---|
| 141 | T am1; |
---|
| 142 | if(x < one) |
---|
| 143 | { |
---|
| 144 | am1 = half * (yy/(r + xp1) + yy/(s - xm1)); |
---|
| 145 | } |
---|
| 146 | else |
---|
| 147 | { |
---|
| 148 | am1 = half * (yy/(r + xp1) + (s + xm1)); |
---|
| 149 | } |
---|
| 150 | imag = boost::math::log1p(am1 + std::sqrt(am1 * (a + one))); |
---|
| 151 | } |
---|
| 152 | else |
---|
| 153 | { |
---|
| 154 | imag = std::log(a + std::sqrt(a*a - one)); |
---|
| 155 | } |
---|
| 156 | } |
---|
| 157 | else |
---|
| 158 | { |
---|
| 159 | // |
---|
| 160 | // This is the Hull et al exception handling code from Fig 6 of their paper: |
---|
| 161 | // |
---|
| 162 | if(y <= (std::numeric_limits<T>::epsilon() * std::fabs(xm1))) |
---|
| 163 | { |
---|
| 164 | if(x < one) |
---|
| 165 | { |
---|
| 166 | real = std::acos(x); |
---|
| 167 | imag = y / std::sqrt(xp1*(one-x)); |
---|
| 168 | } |
---|
| 169 | else |
---|
| 170 | { |
---|
| 171 | real = 0; |
---|
| 172 | if(((std::numeric_limits<T>::max)() / xp1) > xm1) |
---|
| 173 | { |
---|
| 174 | // xp1 * xm1 won't overflow: |
---|
| 175 | imag = boost::math::log1p(xm1 + std::sqrt(xp1*xm1)); |
---|
| 176 | } |
---|
| 177 | else |
---|
| 178 | { |
---|
| 179 | imag = log_two + std::log(x); |
---|
| 180 | } |
---|
| 181 | } |
---|
| 182 | } |
---|
| 183 | else if(y <= safe_min) |
---|
| 184 | { |
---|
| 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" |
---|
| 187 | // machines where : |
---|
| 188 | // |
---|
| 189 | // E^2 > 8*sqrt(u); with: |
---|
| 190 | // |
---|
| 191 | // E = std::numeric_limits<T>::epsilon() |
---|
| 192 | // u = (std::numeric_limits<T>::min)() |
---|
| 193 | // |
---|
| 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: |
---|
| 197 | // |
---|
| 198 | BOOST_ASSERT(x == 1); |
---|
| 199 | real = std::sqrt(y); |
---|
| 200 | imag = std::sqrt(y); |
---|
| 201 | } |
---|
| 202 | else if(std::numeric_limits<T>::epsilon() * y - one >= x) |
---|
| 203 | { |
---|
| 204 | real = half_pi; |
---|
| 205 | imag = log_two + std::log(y); |
---|
| 206 | } |
---|
| 207 | else if(x > one) |
---|
| 208 | { |
---|
| 209 | real = std::atan(y/x); |
---|
| 210 | T xoy = x/y; |
---|
| 211 | imag = log_two + std::log(y) + half * boost::math::log1p(xoy*xoy); |
---|
| 212 | } |
---|
| 213 | else |
---|
| 214 | { |
---|
| 215 | real = half_pi; |
---|
| 216 | T a = std::sqrt(one + y*y); |
---|
| 217 | imag = half * boost::math::log1p(static_cast<T>(2)*y*(y+a)); |
---|
| 218 | } |
---|
| 219 | } |
---|
| 220 | } |
---|
| 221 | |
---|
| 222 | // |
---|
| 223 | // Finish off by working out the sign of the result: |
---|
| 224 | // |
---|
| 225 | if(z.real() < 0) |
---|
| 226 | real = s_pi - real; |
---|
| 227 | if(z.imag() > 0) |
---|
| 228 | imag = -imag; |
---|
| 229 | |
---|
| 230 | return std::complex<T>(real, imag); |
---|
| 231 | } |
---|
| 232 | |
---|
| 233 | } } // namespaces |
---|
| 234 | |
---|
| 235 | #endif // BOOST_MATH_COMPLEX_ACOS_INCLUDED |
---|