diff --git a/include/boost/math/complex/acos.hpp b/include/boost/math/complex/acos.hpp index 3bced116d..48fe65693 100644 --- a/include/boost/math/complex/acos.hpp +++ b/include/boost/math/complex/acos.hpp @@ -106,8 +106,8 @@ std::complex acos(const std::complex& z) // fortunately the quantities M and u identified by Hull et al (figure 3), // match with the max() and min() methods of numeric_limits. // - T safe_max = std::sqrt((std::numeric_limits::max)()) / T(8); - T safe_min = static_cast(4) * std::sqrt((std::numeric_limits::min)()); + T safe_max = detail::safe_max(static_cast(8)); + T safe_min = detail::safe_min(static_cast(4)); T xp1 = one + x; T xm1 = x - one; diff --git a/include/boost/math/complex/asin.hpp b/include/boost/math/complex/asin.hpp index 553c8a913..2eb00f4dd 100644 --- a/include/boost/math/complex/asin.hpp +++ b/include/boost/math/complex/asin.hpp @@ -116,8 +116,8 @@ inline std::complex asin(const std::complex& z) // fortunately the quantities M and u identified by Hull et al (figure 3), // match with the max() and min() methods of numeric_limits. // - T safe_max = std::sqrt((std::numeric_limits::max)()) / T(8); - T safe_min = static_cast(4) * std::sqrt((std::numeric_limits::min)()); + T safe_max = detail::safe_max(static_cast(8)); + T safe_min = detail::safe_min(static_cast(4)); T xp1 = one + x; T xm1 = x - one; diff --git a/include/boost/math/complex/atanh.hpp b/include/boost/math/complex/atanh.hpp index 6c0a7f254..38b892330 100644 --- a/include/boost/math/complex/atanh.hpp +++ b/include/boost/math/complex/atanh.hpp @@ -50,8 +50,8 @@ std::complex atanh(const std::complex& z) T real, imag; // our results - T safe_upper = std::sqrt((std::numeric_limits::max)()) / two; - T safe_lower = std::sqrt((std::numeric_limits::min)()); + T safe_upper = detail::safe_max(two); + T safe_lower = detail::safe_min(static_cast(2)); // // Begin by handling the special cases specified in C99: @@ -131,13 +131,13 @@ std::complex atanh(const std::complex& z) // without either overflow or underflow in the squared terms. // T alpha = 0; - if(x > safe_upper) + if(x >= safe_upper) { if((x == std::numeric_limits::infinity()) || (y == std::numeric_limits::infinity())) { alpha = 0; } - else if(y > safe_upper) + else if(y >= safe_upper) { // Big x and y: divide alpha through by x*y: alpha = (two/y) / (x/y + y/x); @@ -153,7 +153,7 @@ std::complex atanh(const std::complex& z) alpha = two/x; } } - else if(y > safe_upper) + else if(y >= safe_upper) { if(x > one) { @@ -226,7 +226,10 @@ std::complex atanh(const std::complex& z) // // y^2 is negligible: // - imag = std::atan2(two*y, 1 - x*x); + if((y == zero) && (x == one)) + imag = 0; + else + imag = std::atan2(two*y, 1 - x*x); } imag /= two; if(z.imag() < 0) diff --git a/include/boost/math/complex/details.hpp b/include/boost/math/complex/details.hpp index e9b3d2459..15dd3cc6c 100644 --- a/include/boost/math/complex/details.hpp +++ b/include/boost/math/complex/details.hpp @@ -10,13 +10,17 @@ // inverse trig complex functions, it also contains all the includes // that we need to implement all these functions. // +#include #include #include #include // isnan where available #include -namespace boost{ namespace math{ namespace detail{ +#ifdef BOOST_NO_STDC_NAMESPACE +namespace std{ using ::sqrt; } +#endif +namespace boost{ namespace math{ namespace detail{ template inline bool test_is_nan(T t) @@ -48,6 +52,29 @@ inline std::complex mult_minus_i(const std::complex& t) return std::complex(t.imag(), mult_minus_one(t.real())); } +template +inline T safe_max(T t) +{ + return std::sqrt((std::numeric_limits::max)()) / t; +} +inline long double safe_max(long double t) +{ + // long double sqrt often returns infinity due to + // insufficient internal precision: + return std::sqrt((std::numeric_limits::max)()) / t; +} +template +inline T safe_min(T t) +{ + return std::sqrt((std::numeric_limits::min)()) * t; +} +inline long double safe_min(long double t) +{ + // long double sqrt often returns zero due to + // insufficient internal precision: + return std::sqrt((std::numeric_limits::min)()) * t; +} + } } } // namespaces #endif // BOOST_MATH_COMPLEX_DETAILS_INCLUDED