// 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
+#pragma warning(disable : 6326) // comparison of two constants
#endif
-namespace detail{
+namespace detail {
-template<typename T, typename U>
+template <typename T, typename U>
inline void pow_imp(T& result, const T& t, const U& p, const mpl::false_&)
{
// Compute the pure power of typename T t^p.
typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type;
- if(&result == &t)
+ if (&result == &t)
{
T temp;
pow_imp(temp, t, p, mpl::false_());
}
// This will store the result.
- if(U(p % U(2)) != U(0))
+ if (U(p % U(2)) != U(0))
{
result = t;
}
// The variable x stores the binary powers of t.
T x(t);
- while(U(p2 /= 2) != U(0))
+ while (U(p2 /= 2) != U(0))
{
// Square x for each binary power.
eval_multiply(x, x);
const bool has_binary_power = (U(p2 % U(2)) != U(0));
- if(has_binary_power)
+ if (has_binary_power)
{
// Multiply the result with each binary power contained in the exponent.
eval_multiply(result, x);
}
}
-template<typename T, typename U>
+template <typename T, typename U>
inline void pow_imp(T& result, const T& t, const U& p, const mpl::true_&)
{
// Signed integer power, just take care of the sign then call the unsigned version:
- typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type;
- typedef typename make_unsigned<U>::type ui_type;
+ typedef typename boost::multiprecision::detail::canonical<U, T>::type int_type;
+ typedef typename make_unsigned<U>::type ui_type;
- if(p < 0)
+ if (p < 0)
{
T temp;
temp = static_cast<int_type>(1);
} // namespace detail
-template<typename T, typename U>
+template <typename T, typename U>
inline typename enable_if_c<is_integral<U>::value>::type eval_pow(T& result, const T& t, const U& p)
{
detail::pow_imp(result, t, p, boost::is_signed<U>());
BOOST_ASSERT(&H0F0 != &x);
long tol = boost::multiprecision::detail::digits2<number<T, et_on> >::value();
- T t;
+ T t;
T x_pow_n_div_n_fact(x);
T lim;
eval_ldexp(lim, H0F0, 1 - tol);
- if(eval_get_sign(lim) < 0)
+ if (eval_get_sign(lim) < 0)
lim.negate();
ui_type n;
- const unsigned series_limit =
- boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
- ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
+ 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_0f0(; ; x).
- for(n = 2; n < series_limit; ++n)
+ 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_add(H0F0, x_pow_n_div_n_fact);
bool neg = eval_get_sign(x_pow_n_div_n_fact) < 0;
- if(neg)
+ if (neg)
x_pow_n_div_n_fact.negate();
- if(lim.compare(x_pow_n_div_n_fact) > 0)
+ if (lim.compare(x_pow_n_div_n_fact) > 0)
break;
- if(neg)
+ if (neg)
x_pow_n_div_n_fact.negate();
}
- if(n >= series_limit)
+ if (n >= series_limit)
BOOST_THROW_EXCEPTION(std::runtime_error("H0F0 failed to converge"));
}
BOOST_ASSERT(&H1F0 != &a);
T x_pow_n_div_n_fact(x);
- T pochham_a (a);
- T ap (a);
+ T pochham_a(a);
+ T ap(a);
eval_multiply(H1F0, pochham_a, x_pow_n_div_n_fact);
eval_add(H1F0, si_type(1));
T lim;
eval_ldexp(lim, H1F0, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
- if(eval_get_sign(lim) < 0)
+ if (eval_get_sign(lim) < 0)
lim.negate();
si_type n;
- T term, part;
+ T term, part;
const si_type series_limit =
- boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
- ? 100 : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
+ boost::multiprecision::detail::digits2<number<T, et_on> >::value() < 100
+ ? 100
+ : boost::multiprecision::detail::digits2<number<T, et_on> >::value();
// Series expansion of hyperg_1f0(a; ; x).
- for(n = 2; n < series_limit; n++)
+ 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_multiply(pochham_a, ap);
eval_multiply(term, pochham_a, x_pow_n_div_n_fact);
eval_add(H1F0, term);
- if(eval_get_sign(term) < 0)
+ if (eval_get_sign(term) < 0)
term.negate();
- if(lim.compare(term) >= 0)
+ if (lim.compare(term) >= 0)
break;
}
- if(n >= series_limit)
+ if (n >= series_limit)
BOOST_THROW_EXCEPTION(std::runtime_error("H1F0 failed to converge"));
}
void eval_exp(T& result, const T& x)
{
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The exp function is only valid for floating point types.");
- if(&x == &result)
+ if (&x == &result)
{
T temp;
eval_exp(temp, x);
return;
}
typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
- typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
- typedef typename T::exponent_type exp_type;
+ typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
+ typedef typename T::exponent_type exp_type;
typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
// Handle special arguments.
- int type = eval_fpclassify(x);
+ int type = eval_fpclassify(x);
bool isneg = eval_get_sign(x) < 0;
- if(type == (int)FP_NAN)
+ if (type == (int)FP_NAN)
{
result = x;
- errno = EDOM;
+ errno = EDOM;
return;
}
- else if(type == (int)FP_INFINITE)
+ else if (type == (int)FP_INFINITE)
{
- if(isneg)
+ if (isneg)
result = ui_type(0u);
- else
+ else
result = x;
return;
}
- else if(type == (int)FP_ZERO)
+ else if (type == (int)FP_ZERO)
{
result = ui_type(1);
return;
// Get local copy of argument and force it to be positive.
T xx = x;
T exp_series;
- if(isneg)
+ if (isneg)
xx.negate();
// Check the range of the argument.
- if(xx.compare(si_type(1)) <= 0)
+ if (xx.compare(si_type(1)) <= 0)
{
//
// Use series for exp(x) - 1:
//
T lim;
- if(std::numeric_limits<number<T, et_on> >::is_specialized)
+ if (std::numeric_limits<number<T, et_on> >::is_specialized)
lim = std::numeric_limits<number<T, et_on> >::epsilon().backend();
else
{
}
unsigned k = 2;
exp_series = xx;
- result = si_type(1);
- if(isneg)
+ result = si_type(1);
+ if (isneg)
eval_subtract(result, exp_series);
else
eval_add(result, exp_series);
eval_multiply(exp_series, xx);
eval_divide(exp_series, ui_type(k));
eval_add(result, exp_series);
- while(exp_series.compare(lim) > 0)
+ while (exp_series.compare(lim) > 0)
{
++k;
eval_multiply(exp_series, xx);
eval_divide(exp_series, ui_type(k));
- if(isneg && (k&1))
+ if (isneg && (k & 1))
eval_subtract(result, exp_series);
else
eval_add(result, exp_series);
typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type ll;
eval_trunc(exp_series, x);
eval_convert_to(&ll, exp_series);
- if(x.compare(ll) == 0)
+ if (x.compare(ll) == 0)
{
detail::pow_imp(result, get_constant_e<T>(), ll, mpl::true_());
return;
}
- else if(exp_series.compare(x) == 0)
+ else if (exp_series.compare(x) == 0)
{
- // We have a value that has no fractional part, but is too large to fit
+ // We have a value that has no fractional part, but is too large to fit
// in a long long, in this situation the code below will fail, so
// we're just going to assume that this will overflow:
- if(isneg)
+ if (isneg)
result = ui_type(0);
else
result = std::numeric_limits<number<T> >::has_infinity ? std::numeric_limits<number<T> >::infinity().backend() : (std::numeric_limits<number<T> >::max)().backend();
eval_ldexp(result, result, n);
eval_multiply(exp_series, result);
- if(isneg)
+ if (isneg)
eval_divide(result, ui_type(1), exp_series);
else
result = exp_series;
// log(x) = log(2) * n + log1p(1 + y)
//
typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
- typedef typename T::exponent_type exp_type;
+ typedef typename T::exponent_type exp_type;
typedef typename boost::multiprecision::detail::canonical<exp_type, T>::type canonical_exp_type;
- typedef typename mpl::front<typename T::float_types>::type fp_type;
- int s = eval_signbit(arg);
- switch(eval_fpclassify(arg))
+ typedef typename mpl::front<typename T::float_types>::type fp_type;
+ int s = eval_signbit(arg);
+ switch (eval_fpclassify(arg))
{
case FP_NAN:
result = arg;
- errno = EDOM;
+ errno = EDOM;
return;
case FP_INFINITE:
- if(s) break;
+ if (s)
+ break;
result = arg;
return;
case FP_ZERO:
errno = ERANGE;
return;
}
- if(s)
+ if (s)
{
result = std::numeric_limits<number<T> >::quiet_NaN().backend();
- errno = EDOM;
+ errno = EDOM;
return;
}
exp_type e;
- T t;
+ T t;
eval_frexp(t, arg, &e);
bool alternate = false;
- if(t.compare(fp_type(2) / fp_type(3)) <= 0)
+ if (t.compare(fp_type(2) / fp_type(3)) <= 0)
{
alternate = true;
eval_ldexp(t, t, 1);
--e;
}
-
+
eval_multiply(result, get_constant_ln2<T>(), canonical_exp_type(e));
INSTRUMENT_BACKEND(result);
eval_subtract(t, ui_type(1)); /* -0.3 <= t <= 0.3 */
- if(!alternate)
+ if (!alternate)
t.negate(); /* 0 <= t <= 0.33333 */
T pow = t;
T lim;
T t2;
- if(alternate)
+ if (alternate)
eval_add(result, t);
else
eval_subtract(result, t);
- if(std::numeric_limits<number<T, et_on> >::is_specialized)
+ if (std::numeric_limits<number<T, et_on> >::is_specialized)
eval_multiply(lim, result, std::numeric_limits<number<T, et_on> >::epsilon().backend());
else
eval_ldexp(lim, result, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
- if(eval_get_sign(lim) < 0)
+ if (eval_get_sign(lim) < 0)
lim.negate();
INSTRUMENT_BACKEND(lim);
eval_multiply(pow, t);
eval_divide(t2, pow, k);
INSTRUMENT_BACKEND(t2);
- if(alternate && ((k & 1) != 0))
+ if (alternate && ((k & 1) != 0))
eval_add(result, t2);
else
eval_subtract(result, t2);
INSTRUMENT_BACKEND(result);
- }while(lim.compare(t2) < 0);
+ } while (lim.compare(t2) < 0);
}
template <class T>
const T& get_constant_log10()
{
- static BOOST_MP_THREAD_LOCAL T result;
+ static BOOST_MP_THREAD_LOCAL T result;
static BOOST_MP_THREAD_LOCAL long digits = 0;
#ifndef BOOST_MP_USING_THREAD_LOCAL
static BOOST_MP_THREAD_LOCAL bool b = false;
{
#endif
typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
- T ten;
+ T ten;
ten = ui_type(10u);
eval_log(result, ten);
digits = boost::multiprecision::detail::digits2<number<T> >::value();
eval_divide(result, get_constant_ln2<R>());
}
-template<typename T>
+template <typename T>
inline void eval_pow(T& result, const T& x, const T& a)
{
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The pow function is only valid for floating point types.");
typedef typename boost::multiprecision::detail::canonical<int, T>::type si_type;
- typedef typename mpl::front<typename T::float_types>::type fp_type;
+ typedef typename mpl::front<typename T::float_types>::type fp_type;
- if((&result == &x) || (&result == &a))
+ if ((&result == &x) || (&result == &a))
{
T t;
eval_pow(t, x, a);
return;
}
- if((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
+ if ((a.compare(si_type(1)) == 0) || (x.compare(si_type(1)) == 0))
{
result = x;
return;
}
- if(a.compare(si_type(0)) == 0)
+ if (a.compare(si_type(0)) == 0)
{
result = si_type(1);
return;
int type = eval_fpclassify(x);
- switch(type)
+ switch (type)
{
case FP_ZERO:
- switch(eval_fpclassify(a))
+ switch (eval_fpclassify(a))
{
case FP_ZERO:
result = si_type(1);
case FP_NORMAL:
{
// Need to check for a an odd integer as a special case:
- try
+ try
{
typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type i;
eval_convert_to(&i, a);
- if(a.compare(i) == 0)
+ if (a.compare(i) == 0)
{
- if(eval_signbit(a))
+ if (eval_signbit(a))
{
- if(i & 1)
+ if (i & 1)
{
result = std::numeric_limits<number<T> >::infinity().backend();
- if(eval_signbit(x))
+ if (eval_signbit(x))
result.negate();
errno = ERANGE;
}
else
{
result = std::numeric_limits<number<T> >::infinity().backend();
- errno = ERANGE;
+ errno = ERANGE;
}
}
- else if(i & 1)
+ else if (i & 1)
{
result = x;
}
return;
}
}
- catch(const std::exception&)
+ catch (const std::exception&)
{
// fallthrough..
}
BOOST_FALLTHROUGH;
}
default:
- if(eval_signbit(a))
+ if (eval_signbit(a))
{
result = std::numeric_limits<number<T> >::infinity().backend();
- errno = ERANGE;
+ errno = ERANGE;
}
else
result = x;
return;
case FP_NAN:
result = x;
- errno = ERANGE;
+ errno = ERANGE;
return;
- default: ;
+ default:;
}
int s = eval_get_sign(a);
- if(s == 0)
+ if (s == 0)
{
result = si_type(1);
return;
}
- if(s < 0)
+ if (s < 0)
{
T t, da;
t = a;
eval_divide(result, si_type(1), da);
return;
}
-
+
typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type an;
typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type max_an =
- std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ?
- (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::max)() :
- static_cast<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>(1) << (sizeof(typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type) * CHAR_BIT - 2);
- typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type min_an =
- std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ?
- (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::min)() :
- -min_an;
-
+ std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::max)() : static_cast<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>(1) << (sizeof(typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type) * CHAR_BIT - 2);
+ typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type min_an =
+ std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::is_specialized ? (std::numeric_limits<typename boost::multiprecision::detail::canonical<boost::intmax_t, T>::type>::min)() : -min_an;
T fa;
#ifndef BOOST_NO_EXCEPTIONS
{
#endif
eval_convert_to(&an, a);
- if(a.compare(an) == 0)
+ if (a.compare(an) == 0)
{
detail::pow_imp(result, x, an, mpl::true_());
return;
}
#ifndef BOOST_NO_EXCEPTIONS
}
- catch(const std::exception&)
+ catch (const std::exception&)
{
// conversion failed, just fall through, value is not an integer.
an = (std::numeric_limits<boost::intmax_t>::max)();
}
#endif
- if((eval_get_sign(x) < 0))
+ if ((eval_get_sign(x) < 0))
{
typename boost::multiprecision::detail::canonical<boost::uintmax_t, T>::type aun;
#ifndef BOOST_NO_EXCEPTIONS
{
#endif
eval_convert_to(&aun, a);
- if(a.compare(aun) == 0)
+ if (a.compare(aun) == 0)
{
fa = x;
fa.negate();
eval_pow(result, fa, a);
- if(aun & 1u)
+ if (aun & 1u)
result.negate();
return;
}
#ifndef BOOST_NO_EXCEPTIONS
}
- catch(const std::exception&)
+ catch (const std::exception&)
{
// conversion failed, just fall through, value is not an integer.
}
#endif
eval_floor(result, a);
// -1^INF is a special case in C99:
- if((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
+ if ((x.compare(si_type(-1)) == 0) && (eval_fpclassify(a) == FP_INFINITE))
{
result = si_type(1);
}
- else if(a.compare(result) == 0)
+ else if (a.compare(result) == 0)
{
// exponent is so large we have no fractional part:
- if(x.compare(si_type(-1)) < 0)
+ if (x.compare(si_type(-1)) < 0)
{
result = std::numeric_limits<number<T, et_on> >::infinity().backend();
}
result = si_type(0);
}
}
- else if(type == FP_INFINITE)
+ else if (type == FP_INFINITE)
{
result = std::numeric_limits<number<T, et_on> >::infinity().backend();
}
- else if(std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
+ else if (std::numeric_limits<number<T, et_on> >::has_quiet_NaN)
{
result = std::numeric_limits<number<T, et_on> >::quiet_NaN().backend();
- errno = EDOM;
+ errno = EDOM;
}
else
{
eval_subtract(da, a, an);
- if((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
+ if ((x.compare(fp_type(0.5)) >= 0) && (x.compare(fp_type(0.9)) < 0) && (an < max_an) && (an > min_an))
{
- if(a.compare(fp_type(1e-5f)) <= 0)
+ if (a.compare(fp_type(1e-5f)) <= 0)
{
// Series expansion for small a.
eval_log(t, x);
{
// Series expansion for moderately sized x. Note that for large power of a,
// the power of the integer part of a is calculated using the pown function.
- if(an)
+ if (an)
{
da.negate();
t = si_type(1);
{
// Series expansion for pow(x, a). Note that for large power of a, the power
// of the integer part of a is calculated using the pown function.
- if(an)
+ if (an)
{
eval_log(t, x);
eval_multiply(t, da);
}
}
-template<class T, class A>
+template <class T, class A>
#if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
-inline typename enable_if_c<!is_integral<A>::value, void>::type
+inline typename enable_if_c<!is_integral<A>::value, void>::type
#else
-inline typename enable_if_c<is_compatible_arithmetic_type<A, number<T> >::value && !is_integral<A>::value, void>::type
+inline typename enable_if_c<is_compatible_arithmetic_type<A, number<T> >::value && !is_integral<A>::value, void>::type
#endif
- eval_pow(T& result, const T& x, const A& a)
+eval_pow(T& result, const T& x, const A& a)
{
// Note this one is restricted to float arguments since pow.hpp already has a version for
// integer powers....
- typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
+ 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;
+ cast_type c;
c = a;
eval_pow(result, x, c);
}
-template<class T, class A>
+template <class T, class A>
#if BOOST_WORKAROUND(BOOST_MSVC, < 1800)
inline void
#else
inline typename enable_if_c<is_compatible_arithmetic_type<A, number<T> >::value, void>::type
#endif
- eval_pow(T& result, const A& x, const T& a)
+eval_pow(T& result, const A& x, const T& a)
{
- typedef typename boost::multiprecision::detail::canonical<A, T>::type canonical_type;
+ 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;
+ cast_type c;
c = x;
eval_pow(result, c, a);
}
// Check for pure-integer arguments which can be either signed or unsigned.
typename boost::multiprecision::detail::canonical<typename T::exponent_type, T>::type i;
- T temp;
- try {
+ T temp;
+ try
+ {
eval_trunc(temp, arg);
eval_convert_to(&i, temp);
- if(arg.compare(i) == 0)
+ if (arg.compare(i) == 0)
{
temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
eval_ldexp(result, temp, i);
return;
}
}
- catch(const boost::math::rounding_error&)
- { /* Fallthrough */ }
- catch(const std::runtime_error&)
- { /* Fallthrough */ }
+ catch (const boost::math::rounding_error&)
+ { /* Fallthrough */
+ }
+ catch (const std::runtime_error&)
+ { /* Fallthrough */
+ }
temp = static_cast<typename mpl::front<typename T::unsigned_types>::type>(2u);
eval_pow(result, temp, arg);
}
-namespace detail{
+namespace detail {
+
+template <class T>
+void small_sinh_series(T x, T& result)
+{
+ typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
+ bool neg = eval_get_sign(x) < 0;
+ if (neg)
+ x.negate();
+ T p(x);
+ T mult(x);
+ eval_multiply(mult, x);
+ result = x;
+ ui_type k = 1;
+
+ T lim(x);
+ eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
- template <class T>
- void small_sinh_series(T x, T& result)
+ do
{
- typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
- bool neg = eval_get_sign(x) < 0;
- if(neg)
- x.negate();
- T p(x);
- T mult(x);
- eval_multiply(mult, x);
- result = x;
- ui_type k = 1;
+ eval_multiply(p, mult);
+ eval_divide(p, ++k);
+ eval_divide(p, ++k);
+ eval_add(result, p);
+ } while (p.compare(lim) >= 0);
+ if (neg)
+ result.negate();
+}
- T lim(x);
- eval_ldexp(lim, lim, 1 - boost::multiprecision::detail::digits2<number<T, et_on> >::value());
+template <class T>
+void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
+{
+ typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
+ typedef typename mpl::front<typename T::float_types>::type fp_type;
- do
+ switch (eval_fpclassify(x))
+ {
+ case FP_NAN:
+ errno = EDOM;
+ // fallthrough...
+ case FP_INFINITE:
+ if (p_sinh)
+ *p_sinh = x;
+ if (p_cosh)
{
- eval_multiply(p, mult);
- eval_divide(p, ++k);
- eval_divide(p, ++k);
- eval_add(result, p);
- }while(p.compare(lim) >= 0);
- if(neg)
- result.negate();
+ *p_cosh = x;
+ if (eval_get_sign(x) < 0)
+ p_cosh->negate();
+ }
+ return;
+ case FP_ZERO:
+ if (p_sinh)
+ *p_sinh = x;
+ if (p_cosh)
+ *p_cosh = ui_type(1);
+ return;
+ default:;
}
- template <class T>
- void sinhcosh(const T& x, T* p_sinh, T* p_cosh)
+ bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;
+
+ if (p_cosh || !small_sinh)
{
- typedef typename boost::multiprecision::detail::canonical<unsigned, T>::type ui_type;
- typedef typename mpl::front<typename T::float_types>::type fp_type;
+ T e_px, e_mx;
+ eval_exp(e_px, x);
+ eval_divide(e_mx, ui_type(1), e_px);
+ if (eval_signbit(e_mx) != eval_signbit(e_px))
+ e_mx.negate(); // Handles lack of signed zero in some types
- switch(eval_fpclassify(x))
+ if (p_sinh)
{
- case FP_NAN:
- errno = EDOM;
- // fallthrough...
- case FP_INFINITE:
- if(p_sinh)
- *p_sinh = x;
- if(p_cosh)
+ if (small_sinh)
{
- *p_cosh = x;
- if(eval_get_sign(x) < 0)
- p_cosh->negate();
+ small_sinh_series(x, *p_sinh);
}
- return;
- case FP_ZERO:
- if(p_sinh)
- *p_sinh = x;
- if(p_cosh)
- *p_cosh = ui_type(1);
- return;
- default: ;
- }
-
- bool small_sinh = eval_get_sign(x) < 0 ? x.compare(fp_type(-0.5)) > 0 : x.compare(fp_type(0.5)) < 0;
-
- if(p_cosh || !small_sinh)
- {
- T e_px, e_mx;
- eval_exp(e_px, x);
- eval_divide(e_mx, ui_type(1), e_px);
- if(eval_signbit(e_mx) != eval_signbit(e_px))
- e_mx.negate(); // Handles lack of signed zero in some types
-
- if(p_sinh)
- {
- if(small_sinh)
- {
- small_sinh_series(x, *p_sinh);
- }
- else
- {
- eval_subtract(*p_sinh, e_px, e_mx);
- eval_ldexp(*p_sinh, *p_sinh, -1);
- }
- }
- if(p_cosh)
- {
- eval_add(*p_cosh, e_px, e_mx);
- eval_ldexp(*p_cosh, *p_cosh, -1);
+ else
+ {
+ eval_subtract(*p_sinh, e_px, e_mx);
+ eval_ldexp(*p_sinh, *p_sinh, -1);
}
}
- else
+ if (p_cosh)
{
- small_sinh_series(x, *p_sinh);
+ eval_add(*p_cosh, e_px, e_mx);
+ eval_ldexp(*p_cosh, *p_cosh, -1);
}
}
+ else
+ {
+ small_sinh_series(x, *p_sinh);
+ }
+}
} // namespace detail
inline void eval_tanh(T& result, const T& x)
{
BOOST_STATIC_ASSERT_MSG(number_category<T>::value == number_kind_floating_point, "The tanh function is only valid for floating point types.");
- T c;
- detail::sinhcosh(x, &result, &c);
- if((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE))
- {
- bool s = eval_signbit(result) != eval_signbit(c);
- result = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
- if(s)
- result.negate();
- return;
- }
- eval_divide(result, c);
+ T c;
+ detail::sinhcosh(x, &result, &c);
+ if ((eval_fpclassify(result) == FP_INFINITE) && (eval_fpclassify(c) == FP_INFINITE))
+ {
+ bool s = eval_signbit(result) != eval_signbit(c);
+ result = static_cast<typename mpl::front<typename T::unsigned_types>::type>(1u);
+ if (s)
+ result.negate();
+ return;
+ }
+ eval_divide(result, c);
}
#ifdef BOOST_MSVC
projects / platform / upstream / boost.git / blobdiff