793 lines
23 KiB
C++
793 lines
23 KiB
C++
|
|
||
|
|
// Copyright Christopher Kormanyos 2002 - 2011.
|
||
|
|
// Copyright 2011 John Maddock. Distributed under the Boost
|
||
|
|
// Distributed under the Boost Software License, Version 1.0.
|
||
|
|
// (See accompanying file LICENSE_1_0.txt or copy at
|
||
|
|
// http://www.boost.org/LICENSE_1_0.txt)
|
||
|
|
|
||
|
|
// This work is based on an earlier work:
|
||
|
|
// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
|
||
|
|
// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
|
||
|
|
//
|
||
|
|
// This file has no include guards or namespaces - it's expanded inline inside default_ops.hpp
|
||
|
|
//
|
||
|
|
|
||
|
|
#ifdef BOOST_MSVC
|
||
|
|
#pragma warning(push)
|
||
|
|
#pragma warning(disable:6326) // comparison of two constants
|
||
|
|
#endif
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void hyp0F1(T& result, const T& b, const T& x)
|
||
|
|
{
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
|
||
|
|
// Compute the series representation of Hypergeometric0F1 taken from
|
||
|
|
// http://functions.wolfram.com/HypergeometricFunctions/Hypergeometric0F1/06/01/01/
|
||
|
|
// There are no checks on input range or parameter boundaries.
|
||
|
|
|
||
|
|
T x_pow_n_div_n_fact(x);
|
||
|
|
T pochham_b (b);
|
||
|
|
T bp (b);
|
||
|
|
|
||
|
|
eval_divide(result, x_pow_n_div_n_fact, pochham_b);
|
||
|
|
eval_add(result, ui_type(1));
|
||
|
|
|
||
|
|
si_type n;
|
||
|
|
|
||
|
|
T tol;
|
||
|
|
tol = ui_type(1);
|
||
|
|
eval_ldexp(tol, tol, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
|
||
|
|
eval_multiply(tol, result);
|
||
|
|
if(eval_get_sign(tol) < 0)
|
||
|
|
tol.negate();
|
||
|
|
T term;
|
||
|
|
|
||
|
|
const int series_limit =
|
||
|
|
boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
|
||
|
|
? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
|
||
|
|
// Series expansion of hyperg_0f1(; b; x).
|
||
|
|
for(n = 2; n < series_limit; ++n)
|
||
|
|
{
|
||
|
|
eval_multiply(x_pow_n_div_n_fact, x);
|
||
|
|
eval_divide(x_pow_n_div_n_fact, n);
|
||
|
|
eval_increment(bp);
|
||
|
|
eval_multiply(pochham_b, bp);
|
||
|
|
|
||
|
|
eval_divide(term, x_pow_n_div_n_fact, pochham_b);
|
||
|
|
eval_add(result, term);
|
||
|
|
|
||
|
|
bool neg_term = eval_get_sign(term) < 0;
|
||
|
|
if(neg_term)
|
||
|
|
term.negate();
|
||
|
|
if(term.compare(tol) <= 0)
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
|
||
|
|
if(n >= series_limit)
|
||
|
|
BOOST_THROW_EXCEPTION(std::runtime_error("H0F1 Failed to Converge"));
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void eval_sin(T& result, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The sin function is only valid for floating point types.");
|
||
|
|
if(&result == &x)
|
||
|
|
{
|
||
|
|
T temp;
|
||
|
|
eval_sin(temp, x);
|
||
|
|
result = temp;
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
typedef typename mpl::front<typename T::float_types>::type fp_type;
|
||
|
|
|
||
|
|
switch(eval_fpclassify(x))
|
||
|
|
{
|
||
|
|
case FP_INFINITE:
|
||
|
|
case FP_NAN:
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
result = ui_type(0);
|
||
|
|
return;
|
||
|
|
default: ;
|
||
|
|
}
|
||
|
|
|
||
|
|
// Local copy of the argument
|
||
|
|
T xx = x;
|
||
|
|
|
||
|
|
// Analyze and prepare the phase of the argument.
|
||
|
|
// Make a local, positive copy of the argument, xx.
|
||
|
|
// The argument xx will be reduced to 0 <= xx <= pi/2.
|
||
|
|
bool b_negate_sin = false;
|
||
|
|
|
||
|
|
if(eval_get_sign(x) < 0)
|
||
|
|
{
|
||
|
|
xx.negate();
|
||
|
|
b_negate_sin = !b_negate_sin;
|
||
|
|
}
|
||
|
|
|
||
|
|
T n_pi, t;
|
||
|
|
// Remove even multiples of pi.
|
||
|
|
if(xx.compare(get_constant_pi<T>()) > 0)
|
||
|
|
{
|
||
|
|
eval_divide(n_pi, xx, get_constant_pi<T>());
|
||
|
|
eval_trunc(n_pi, n_pi);
|
||
|
|
t = ui_type(2);
|
||
|
|
eval_fmod(t, n_pi, t);
|
||
|
|
const bool b_n_pi_is_even = eval_get_sign(t) == 0;
|
||
|
|
eval_multiply(n_pi, get_constant_pi<T>());
|
||
|
|
eval_subtract(xx, n_pi);
|
||
|
|
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
|
||
|
|
|
||
|
|
// Adjust signs if the multiple of pi is not even.
|
||
|
|
if(!b_n_pi_is_even)
|
||
|
|
{
|
||
|
|
b_negate_sin = !b_negate_sin;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// Reduce the argument to 0 <= xx <= pi/2.
|
||
|
|
eval_ldexp(t, get_constant_pi<T>(), -1);
|
||
|
|
if(xx.compare(t) > 0)
|
||
|
|
{
|
||
|
|
eval_subtract(xx, get_constant_pi<T>(), xx);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
||
|
|
}
|
||
|
|
|
||
|
|
eval_subtract(t, xx);
|
||
|
|
const bool b_zero = eval_get_sign(xx) == 0;
|
||
|
|
const bool b_pi_half = eval_get_sign(t) == 0;
|
||
|
|
|
||
|
|
// Check if the reduced argument is very close to 0 or pi/2.
|
||
|
|
const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
|
||
|
|
const bool b_near_pi_half = t.compare(fp_type(1e-1)) < 0;;
|
||
|
|
|
||
|
|
if(b_zero)
|
||
|
|
{
|
||
|
|
result = ui_type(0);
|
||
|
|
}
|
||
|
|
else if(b_pi_half)
|
||
|
|
{
|
||
|
|
result = ui_type(1);
|
||
|
|
}
|
||
|
|
else if(b_near_zero)
|
||
|
|
{
|
||
|
|
eval_multiply(t, xx, xx);
|
||
|
|
eval_divide(t, si_type(-4));
|
||
|
|
T t2;
|
||
|
|
t2 = fp_type(1.5);
|
||
|
|
hyp0F1(result, t2, t);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
||
|
|
eval_multiply(result, xx);
|
||
|
|
}
|
||
|
|
else if(b_near_pi_half)
|
||
|
|
{
|
||
|
|
eval_multiply(t, t);
|
||
|
|
eval_divide(t, si_type(-4));
|
||
|
|
T t2;
|
||
|
|
t2 = fp_type(0.5);
|
||
|
|
hyp0F1(result, t2, t);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
// Scale to a small argument for an efficient Taylor series,
|
||
|
|
// implemented as a hypergeometric function. Use a standard
|
||
|
|
// divide by three identity a certain number of times.
|
||
|
|
// Here we use division by 3^9 --> (19683 = 3^9).
|
||
|
|
|
||
|
|
static const si_type n_scale = 9;
|
||
|
|
static const si_type n_three_pow_scale = static_cast<si_type>(19683L);
|
||
|
|
|
||
|
|
eval_divide(xx, n_three_pow_scale);
|
||
|
|
|
||
|
|
// Now with small arguments, we are ready for a series expansion.
|
||
|
|
eval_multiply(t, xx, xx);
|
||
|
|
eval_divide(t, si_type(-4));
|
||
|
|
T t2;
|
||
|
|
t2 = fp_type(1.5);
|
||
|
|
hyp0F1(result, t2, t);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
||
|
|
eval_multiply(result, xx);
|
||
|
|
|
||
|
|
// Convert back using multiple angle identity.
|
||
|
|
for(boost::int32_t k = static_cast<boost::int32_t>(0); k < n_scale; k++)
|
||
|
|
{
|
||
|
|
// Rescale the cosine value using the multiple angle identity.
|
||
|
|
eval_multiply(t2, result, ui_type(3));
|
||
|
|
eval_multiply(t, result, result);
|
||
|
|
eval_multiply(t, result);
|
||
|
|
eval_multiply(t, ui_type(4));
|
||
|
|
eval_subtract(result, t2, t);
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
if(b_negate_sin)
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void eval_cos(T& result, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The cos function is only valid for floating point types.");
|
||
|
|
if(&result == &x)
|
||
|
|
{
|
||
|
|
T temp;
|
||
|
|
eval_cos(temp, x);
|
||
|
|
result = temp;
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
typedef typename mpl::front<typename T::float_types>::type fp_type;
|
||
|
|
|
||
|
|
switch(eval_fpclassify(x))
|
||
|
|
{
|
||
|
|
case FP_INFINITE:
|
||
|
|
case FP_NAN:
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
result = ui_type(1);
|
||
|
|
return;
|
||
|
|
default: ;
|
||
|
|
}
|
||
|
|
|
||
|
|
// Local copy of the argument
|
||
|
|
T xx = x;
|
||
|
|
|
||
|
|
// Analyze and prepare the phase of the argument.
|
||
|
|
// Make a local, positive copy of the argument, xx.
|
||
|
|
// The argument xx will be reduced to 0 <= xx <= pi/2.
|
||
|
|
bool b_negate_cos = false;
|
||
|
|
|
||
|
|
if(eval_get_sign(x) < 0)
|
||
|
|
{
|
||
|
|
xx.negate();
|
||
|
|
}
|
||
|
|
|
||
|
|
T n_pi, t;
|
||
|
|
// Remove even multiples of pi.
|
||
|
|
if(xx.compare(get_constant_pi<T>()) > 0)
|
||
|
|
{
|
||
|
|
eval_divide(t, xx, get_constant_pi<T>());
|
||
|
|
eval_trunc(n_pi, t);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(n_pi.str(0, std::ios_base::scientific));
|
||
|
|
eval_multiply(t, n_pi, get_constant_pi<T>());
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(t.str(0, std::ios_base::scientific));
|
||
|
|
eval_subtract(xx, t);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
||
|
|
|
||
|
|
// Adjust signs if the multiple of pi is not even.
|
||
|
|
t = ui_type(2);
|
||
|
|
eval_fmod(t, n_pi, t);
|
||
|
|
const bool b_n_pi_is_even = eval_get_sign(t) == 0;
|
||
|
|
|
||
|
|
if(!b_n_pi_is_even)
|
||
|
|
{
|
||
|
|
b_negate_cos = !b_negate_cos;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
// Reduce the argument to 0 <= xx <= pi/2.
|
||
|
|
eval_ldexp(t, get_constant_pi<T>(), -1);
|
||
|
|
int com = xx.compare(t);
|
||
|
|
if(com > 0)
|
||
|
|
{
|
||
|
|
eval_subtract(xx, get_constant_pi<T>(), xx);
|
||
|
|
b_negate_cos = !b_negate_cos;
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(xx.str(0, std::ios_base::scientific));
|
||
|
|
}
|
||
|
|
|
||
|
|
const bool b_zero = eval_get_sign(xx) == 0;
|
||
|
|
const bool b_pi_half = com == 0;
|
||
|
|
|
||
|
|
// Check if the reduced argument is very close to 0.
|
||
|
|
const bool b_near_zero = xx.compare(fp_type(1e-1)) < 0;
|
||
|
|
|
||
|
|
if(b_zero)
|
||
|
|
{
|
||
|
|
result = si_type(1);
|
||
|
|
}
|
||
|
|
else if(b_pi_half)
|
||
|
|
{
|
||
|
|
result = si_type(0);
|
||
|
|
}
|
||
|
|
else if(b_near_zero)
|
||
|
|
{
|
||
|
|
eval_multiply(t, xx, xx);
|
||
|
|
eval_divide(t, si_type(-4));
|
||
|
|
n_pi = fp_type(0.5f);
|
||
|
|
hyp0F1(result, n_pi, t);
|
||
|
|
BOOST_MATH_INSTRUMENT_CODE(result.str(0, std::ios_base::scientific));
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
eval_subtract(t, xx);
|
||
|
|
eval_sin(result, t);
|
||
|
|
}
|
||
|
|
if(b_negate_cos)
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void eval_tan(T& result, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tan function is only valid for floating point types.");
|
||
|
|
if(&result == &x)
|
||
|
|
{
|
||
|
|
T temp;
|
||
|
|
eval_tan(temp, x);
|
||
|
|
result = temp;
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
T t;
|
||
|
|
eval_sin(result, x);
|
||
|
|
eval_cos(t, x);
|
||
|
|
eval_divide(result, t);
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void hyp2F1(T& result, const T& a, const T& b, const T& c, const T& x)
|
||
|
|
{
|
||
|
|
// Compute the series representation of hyperg_2f1 taken from
|
||
|
|
// Abramowitz and Stegun 15.1.1.
|
||
|
|
// There are no checks on input range or parameter boundaries.
|
||
|
|
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
|
||
|
|
T x_pow_n_div_n_fact(x);
|
||
|
|
T pochham_a (a);
|
||
|
|
T pochham_b (b);
|
||
|
|
T pochham_c (c);
|
||
|
|
T ap (a);
|
||
|
|
T bp (b);
|
||
|
|
T cp (c);
|
||
|
|
|
||
|
|
eval_multiply(result, pochham_a, pochham_b);
|
||
|
|
eval_divide(result, pochham_c);
|
||
|
|
eval_multiply(result, x_pow_n_div_n_fact);
|
||
|
|
eval_add(result, ui_type(1));
|
||
|
|
|
||
|
|
T lim;
|
||
|
|
eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
|
||
|
|
|
||
|
|
if(eval_get_sign(lim) < 0)
|
||
|
|
lim.negate();
|
||
|
|
|
||
|
|
ui_type n;
|
||
|
|
T term;
|
||
|
|
|
||
|
|
const unsigned series_limit =
|
||
|
|
boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
|
||
|
|
? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
|
||
|
|
// Series expansion of hyperg_2f1(a, b; c; x).
|
||
|
|
for(n = 2; n < series_limit; ++n)
|
||
|
|
{
|
||
|
|
eval_multiply(x_pow_n_div_n_fact, x);
|
||
|
|
eval_divide(x_pow_n_div_n_fact, n);
|
||
|
|
|
||
|
|
eval_increment(ap);
|
||
|
|
eval_multiply(pochham_a, ap);
|
||
|
|
eval_increment(bp);
|
||
|
|
eval_multiply(pochham_b, bp);
|
||
|
|
eval_increment(cp);
|
||
|
|
eval_multiply(pochham_c, cp);
|
||
|
|
|
||
|
|
eval_multiply(term, pochham_a, pochham_b);
|
||
|
|
eval_divide(term, pochham_c);
|
||
|
|
eval_multiply(term, x_pow_n_div_n_fact);
|
||
|
|
eval_add(result, term);
|
||
|
|
|
||
|
|
if(eval_get_sign(term) < 0)
|
||
|
|
term.negate();
|
||
|
|
if(lim.compare(term) >= 0)
|
||
|
|
break;
|
||
|
|
}
|
||
|
|
if(n > series_limit)
|
||
|
|
BOOST_THROW_EXCEPTION(std::runtime_error("H2F1 failed to converge."));
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void eval_asin(T& result, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The asin function is only valid for floating point types.");
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
typedef typename mpl::front<typename T::float_types>::type fp_type;
|
||
|
|
|
||
|
|
if(&result == &x)
|
||
|
|
{
|
||
|
|
T t(x);
|
||
|
|
eval_asin(result, t);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
switch(eval_fpclassify(x))
|
||
|
|
{
|
||
|
|
case FP_NAN:
|
||
|
|
case FP_INFINITE:
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
result = ui_type(0);
|
||
|
|
return;
|
||
|
|
default: ;
|
||
|
|
}
|
||
|
|
|
||
|
|
const bool b_neg = eval_get_sign(x) < 0;
|
||
|
|
|
||
|
|
T xx(x);
|
||
|
|
if(b_neg)
|
||
|
|
xx.negate();
|
||
|
|
|
||
|
|
int c = xx.compare(ui_type(1));
|
||
|
|
if(c > 0)
|
||
|
|
{
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
else if(c == 0)
|
||
|
|
{
|
||
|
|
result = get_constant_pi<T>();
|
||
|
|
eval_ldexp(result, result, -1);
|
||
|
|
if(b_neg)
|
||
|
|
result.negate();
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
if(xx.compare(fp_type(1e-4)) < 0)
|
||
|
|
{
|
||
|
|
// http://functions.wolfram.com/ElementaryFunctions/ArcSin/26/01/01/
|
||
|
|
eval_multiply(xx, xx);
|
||
|
|
T t1, t2;
|
||
|
|
t1 = fp_type(0.5f);
|
||
|
|
t2 = fp_type(1.5f);
|
||
|
|
hyp2F1(result, t1, t1, t2, xx);
|
||
|
|
eval_multiply(result, x);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
else if(xx.compare(fp_type(1 - 1e-4f)) > 0)
|
||
|
|
{
|
||
|
|
T dx1;
|
||
|
|
T t1, t2;
|
||
|
|
eval_subtract(dx1, ui_type(1), xx);
|
||
|
|
t1 = fp_type(0.5f);
|
||
|
|
t2 = fp_type(1.5f);
|
||
|
|
eval_ldexp(dx1, dx1, -1);
|
||
|
|
hyp2F1(result, t1, t1, t2, dx1);
|
||
|
|
eval_ldexp(dx1, dx1, 2);
|
||
|
|
eval_sqrt(t1, dx1);
|
||
|
|
eval_multiply(result, t1);
|
||
|
|
eval_ldexp(t1, get_constant_pi<T>(), -1);
|
||
|
|
result.negate();
|
||
|
|
eval_add(result, t1);
|
||
|
|
if(b_neg)
|
||
|
|
result.negate();
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<long double, T>::type guess_type;
|
||
|
|
#else
|
||
|
|
typedef fp_type guess_type;
|
||
|
|
#endif
|
||
|
|
// Get initial estimate using standard math function asin.
|
||
|
|
guess_type dd;
|
||
|
|
eval_convert_to(&dd, xx);
|
||
|
|
|
||
|
|
result = (guess_type)(std::asin(dd));
|
||
|
|
|
||
|
|
// Newton-Raphson iteration, we should double our precision with each iteration,
|
||
|
|
// in practice this seems to not quite work in all cases... so terminate when we
|
||
|
|
// have at least 2/3 of the digits correct on the assumption that the correction
|
||
|
|
// we've just added will finish the job...
|
||
|
|
|
||
|
|
boost::intmax_t current_precision = eval_ilogb(result);
|
||
|
|
boost::intmax_t target_precision = current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3;
|
||
|
|
|
||
|
|
// Newton-Raphson iteration
|
||
|
|
while(current_precision > target_precision)
|
||
|
|
{
|
||
|
|
T sine, cosine;
|
||
|
|
eval_sin(sine, result);
|
||
|
|
eval_cos(cosine, result);
|
||
|
|
eval_subtract(sine, xx);
|
||
|
|
eval_divide(sine, cosine);
|
||
|
|
eval_subtract(result, sine);
|
||
|
|
current_precision = eval_ilogb(sine);
|
||
|
|
#ifdef FP_ILOGB0
|
||
|
|
if(current_precision == FP_ILOGB0)
|
||
|
|
break;
|
||
|
|
#endif
|
||
|
|
}
|
||
|
|
if(b_neg)
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
inline void eval_acos(T& result, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The acos function is only valid for floating point types.");
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
|
||
|
|
switch(eval_fpclassify(x))
|
||
|
|
{
|
||
|
|
case FP_NAN:
|
||
|
|
case FP_INFINITE:
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
result = get_constant_pi<T>();
|
||
|
|
eval_ldexp(result, result, -1); // divide by two.
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
eval_abs(result, x);
|
||
|
|
int c = result.compare(ui_type(1));
|
||
|
|
|
||
|
|
if(c > 0)
|
||
|
|
{
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
else if(c == 0)
|
||
|
|
{
|
||
|
|
if(eval_get_sign(x) < 0)
|
||
|
|
result = get_constant_pi<T>();
|
||
|
|
else
|
||
|
|
result = ui_type(0);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
eval_asin(result, x);
|
||
|
|
T t;
|
||
|
|
eval_ldexp(t, get_constant_pi<T>(), -1);
|
||
|
|
eval_subtract(result, t);
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void eval_atan(T& result, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan function is only valid for floating point types.");
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::int32_t, T>::type si_type;
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
typedef typename mpl::front<typename T::float_types>::type fp_type;
|
||
|
|
|
||
|
|
switch(eval_fpclassify(x))
|
||
|
|
{
|
||
|
|
case FP_NAN:
|
||
|
|
result = x;
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
result = ui_type(0);
|
||
|
|
return;
|
||
|
|
case FP_INFINITE:
|
||
|
|
if(eval_get_sign(x) < 0)
|
||
|
|
{
|
||
|
|
eval_ldexp(result, get_constant_pi<T>(), -1);
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
else
|
||
|
|
eval_ldexp(result, get_constant_pi<T>(), -1);
|
||
|
|
return;
|
||
|
|
default: ;
|
||
|
|
}
|
||
|
|
|
||
|
|
const bool b_neg = eval_get_sign(x) < 0;
|
||
|
|
|
||
|
|
T xx(x);
|
||
|
|
if(b_neg)
|
||
|
|
xx.negate();
|
||
|
|
|
||
|
|
if(xx.compare(fp_type(0.1)) < 0)
|
||
|
|
{
|
||
|
|
T t1, t2, t3;
|
||
|
|
t1 = ui_type(1);
|
||
|
|
t2 = fp_type(0.5f);
|
||
|
|
t3 = fp_type(1.5f);
|
||
|
|
eval_multiply(xx, xx);
|
||
|
|
xx.negate();
|
||
|
|
hyp2F1(result, t1, t2, t3, xx);
|
||
|
|
eval_multiply(result, x);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
if(xx.compare(fp_type(10)) > 0)
|
||
|
|
{
|
||
|
|
T t1, t2, t3;
|
||
|
|
t1 = fp_type(0.5f);
|
||
|
|
t2 = ui_type(1u);
|
||
|
|
t3 = fp_type(1.5f);
|
||
|
|
eval_multiply(xx, xx);
|
||
|
|
eval_divide(xx, si_type(-1), xx);
|
||
|
|
hyp2F1(result, t1, t2, t3, xx);
|
||
|
|
eval_divide(result, x);
|
||
|
|
if(!b_neg)
|
||
|
|
result.negate();
|
||
|
|
eval_ldexp(t1, get_constant_pi<T>(), -1);
|
||
|
|
eval_add(result, t1);
|
||
|
|
if(b_neg)
|
||
|
|
result.negate();
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
// Get initial estimate using standard math function atan.
|
||
|
|
fp_type d;
|
||
|
|
eval_convert_to(&d, xx);
|
||
|
|
result = fp_type(std::atan(d));
|
||
|
|
|
||
|
|
// Newton-Raphson iteration, we should double our precision with each iteration,
|
||
|
|
// in practice this seems to not quite work in all cases... so terminate when we
|
||
|
|
// have at least 2/3 of the digits correct on the assumption that the correction
|
||
|
|
// we've just added will finish the job...
|
||
|
|
|
||
|
|
boost::intmax_t current_precision = eval_ilogb(result);
|
||
|
|
boost::intmax_t target_precision = current_precision - 1 - (std::numeric_limits<number<T> >::digits * 2) / 3;
|
||
|
|
|
||
|
|
T s, c, t;
|
||
|
|
while(current_precision > target_precision)
|
||
|
|
{
|
||
|
|
eval_sin(s, result);
|
||
|
|
eval_cos(c, result);
|
||
|
|
eval_multiply(t, xx, c);
|
||
|
|
eval_subtract(t, s);
|
||
|
|
eval_multiply(s, t, c);
|
||
|
|
eval_add(result, s);
|
||
|
|
current_precision = eval_ilogb(s);
|
||
|
|
#ifdef FP_ILOGB0
|
||
|
|
if(current_precision == FP_ILOGB0)
|
||
|
|
break;
|
||
|
|
#endif
|
||
|
|
}
|
||
|
|
if(b_neg)
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
|
||
|
|
template <class T>
|
||
|
|
void eval_atan2(T& result, const T& y, const T& x)
|
||
|
|
{
|
||
|
|
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The atan2 function is only valid for floating point types.");
|
||
|
|
if(&result == &y)
|
||
|
|
{
|
||
|
|
T temp(y);
|
||
|
|
eval_atan2(result, temp, x);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
else if(&result == &x)
|
||
|
|
{
|
||
|
|
T temp(x);
|
||
|
|
eval_atan2(result, y, temp);
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<boost::uint32_t, T>::type ui_type;
|
||
|
|
|
||
|
|
switch(eval_fpclassify(y))
|
||
|
|
{
|
||
|
|
case FP_NAN:
|
||
|
|
result = y;
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
{
|
||
|
|
int c = eval_get_sign(x);
|
||
|
|
if(c < 0)
|
||
|
|
result = get_constant_pi<T>();
|
||
|
|
else if(c >= 0)
|
||
|
|
result = ui_type(0); // Note we allow atan2(0,0) to be zero, even though it's mathematically undefined
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
case FP_INFINITE:
|
||
|
|
{
|
||
|
|
if(eval_fpclassify(x) == FP_INFINITE)
|
||
|
|
{
|
||
|
|
if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
|
||
|
|
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
|
||
|
|
else
|
||
|
|
BOOST_THROW_EXCEPTION(std::domain_error("Result is undefined or complex and there is no NaN for this number type."));
|
||
|
|
}
|
||
|
|
else
|
||
|
|
{
|
||
|
|
eval_ldexp(result, get_constant_pi<T>(), -1);
|
||
|
|
if(eval_get_sign(y) < 0)
|
||
|
|
result.negate();
|
||
|
|
}
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
}
|
||
|
|
|
||
|
|
switch(eval_fpclassify(x))
|
||
|
|
{
|
||
|
|
case FP_NAN:
|
||
|
|
result = x;
|
||
|
|
return;
|
||
|
|
case FP_ZERO:
|
||
|
|
{
|
||
|
|
eval_ldexp(result, get_constant_pi<T>(), -1);
|
||
|
|
if(eval_get_sign(y) < 0)
|
||
|
|
result.negate();
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
case FP_INFINITE:
|
||
|
|
if(eval_get_sign(x) > 0)
|
||
|
|
result = ui_type(0);
|
||
|
|
else
|
||
|
|
result = get_constant_pi<T>();
|
||
|
|
if(eval_get_sign(y) < 0)
|
||
|
|
result.negate();
|
||
|
|
return;
|
||
|
|
}
|
||
|
|
|
||
|
|
T xx;
|
||
|
|
eval_divide(xx, y, x);
|
||
|
|
if(eval_get_sign(xx) < 0)
|
||
|
|
xx.negate();
|
||
|
|
|
||
|
|
eval_atan(result, xx);
|
||
|
|
|
||
|
|
// Determine quadrant (sign) based on signs of x, y
|
||
|
|
const bool y_neg = eval_get_sign(y) < 0;
|
||
|
|
const bool x_neg = eval_get_sign(x) < 0;
|
||
|
|
|
||
|
|
if(y_neg != x_neg)
|
||
|
|
result.negate();
|
||
|
|
|
||
|
|
if(x_neg)
|
||
|
|
{
|
||
|
|
if(y_neg)
|
||
|
|
eval_subtract(result, get_constant_pi<T>());
|
||
|
|
else
|
||
|
|
eval_add(result, get_constant_pi<T>());
|
||
|
|
}
|
||
|
|
}
|
||
|
|
template<class T, class A>
|
||
|
|
inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const T& x, const A& a)
|
||
|
|
{
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
|
||
|
|
typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
|
||
|
|
cast_type c;
|
||
|
|
c = a;
|
||
|
|
eval_atan2(result, x, c);
|
||
|
|
}
|
||
|
|
|
||
|
|
template<class T, class A>
|
||
|
|
inline typename enable_if<is_arithmetic<A>, void>::type eval_atan2(T& result, const A& x, const T& a)
|
||
|
|
{
|
||
|
|
typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
|
||
|
|
typedef typename mpl::if_<is_same<A, canonical_type>, T, canonical_type>::type cast_type;
|
||
|
|
cast_type c;
|
||
|
|
c = x;
|
||
|
|
eval_atan2(result, c, a);
|
||
|
|
}
|
||
|
|
|
||
|
|
#ifdef BOOST_MSVC
|
||
|
|
#pragma warning(pop)
|
||
|
|
#endif
|