boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-13 12:32:15 +00:00
Code Issues Projects Releases Wiki Activity

Changed calculation of upper and lower safe bounds to avoid problems with std::sqrt(long double) returning infinities.

Browse Source 
[SVN r31885]
This commit is contained in:
John Maddock
2005-12-02 19:23:28 +00:00
parent cb13210a12
commit 0cf71150fb
4 changed files with 41 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
include/boost/math/complex/acos.hpp
View File

@@ -106,8 +106,8 @@ std::complex<T> acos(const std::complex<T>& 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>.
//
T safe_max = std::sqrt((std::numeric_limits<T>::max)()) / T(8);
T safe_min = static_cast<T>(4) * std::sqrt((std::numeric_limits<T>::min)());
T safe_max = detail::safe_max(static_cast<T>(8));
T safe_min = detail::safe_min(static_cast<T>(4));
T xp1 = one + x;
T xm1 = x - one;

4
include/boost/math/complex/asin.hpp
View File

@@ -116,8 +116,8 @@ inline std::complex<T> asin(const std::complex<T>& 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>.
//
T safe_max = std::sqrt((std::numeric_limits<T>::max)()) / T(8);
T safe_min = static_cast<T>(4) * std::sqrt((std::numeric_limits<T>::min)());
T safe_max = detail::safe_max(static_cast<T>(8));
T safe_min = detail::safe_min(static_cast<T>(4));
T xp1 = one + x;
T xm1 = x - one;

15
include/boost/math/complex/atanh.hpp
View File

@@ -50,8 +50,8 @@ std::complex<T> atanh(const std::complex<T>& z)
T real, imag; // our results
T safe_upper = std::sqrt((std::numeric_limits<T>::max)()) / two;
T safe_lower = std::sqrt((std::numeric_limits<T>::min)());
T safe_upper = detail::safe_max(two);
T safe_lower = detail::safe_min(static_cast<T>(2));
//
// Begin by handling the special cases specified in C99:
@@ -131,13 +131,13 @@ std::complex<T> atanh(const std::complex<T>& 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<T>::infinity()) || (y == std::numeric_limits<T>::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<T> atanh(const std::complex<T>& z)
alpha = two/x;
}
}
else if(y > safe_upper)
else if(y >= safe_upper)
{
if(x > one)
{
@@ -226,7 +226,10 @@ std::complex<T> atanh(const std::complex<T>& 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)

29
include/boost/math/complex/details.hpp
View File

@@ -10,13 +10,17 @@
// inverse trig complex functions, it also contains all the includes
// that we need to implement all these functions.
//
#include <boost/config.hpp>
#include <boost/config/no_tr1/complex.hpp>
#include <boost/limits.hpp>
#include <math.h> // isnan where available
#include <cmath>
namespace boost{ namespace math{ namespace detail{
#ifdef BOOST_NO_STDC_NAMESPACE
namespace std{ using ::sqrt; }
#endif
namespace boost{ namespace math{ namespace detail{
template <class T>
inline bool test_is_nan(T t)
@@ -48,6 +52,29 @@ inline std::complex<T> mult_minus_i(const std::complex<T>& t)
return std::complex<T>(t.imag(), mult_minus_one(t.real()));
}
template <class T>
inline T safe_max(T t)
{
return std::sqrt((std::numeric_limits<T>::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<double>::max)()) / t;
}
template <class T>
inline T safe_min(T t)
{
return std::sqrt((std::numeric_limits<T>::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<double>::min)()) * t;
}
} } } // namespaces
#endif // BOOST_MATH_COMPLEX_DETAILS_INCLUDED