diff options
Diffstat (limited to 'boost/math')
44 files changed, 3434 insertions, 563 deletions
diff --git a/boost/math/concepts/real_concept.hpp b/boost/math/concepts/real_concept.hpp index 2267271a00..2b105e6c85 100644 --- a/boost/math/concepts/real_concept.hpp +++ b/boost/math/concepts/real_concept.hpp @@ -328,7 +328,7 @@ namespace tools { template <> -inline concepts::real_concept make_big_value<concepts::real_concept>(boost::floatmax_t val, const char* , mpl::false_ const&, mpl::false_ const&) +inline concepts::real_concept make_big_value<concepts::real_concept>(boost::math::tools::largest_float val, const char* , mpl::false_ const&, mpl::false_ const&) { return val; // Can't use lexical_cast here, sometimes it fails.... } diff --git a/boost/math/concepts/std_real_concept.hpp b/boost/math/concepts/std_real_concept.hpp index b4f75bcadb..c565ce3773 100644 --- a/boost/math/concepts/std_real_concept.hpp +++ b/boost/math/concepts/std_real_concept.hpp @@ -342,7 +342,7 @@ namespace tools { template <> -inline concepts::std_real_concept make_big_value<concepts::std_real_concept>(boost::floatmax_t val, const char*, mpl::false_ const&, mpl::false_ const&) +inline concepts::std_real_concept make_big_value<concepts::std_real_concept>(boost::math::tools::largest_float val, const char*, mpl::false_ const&, mpl::false_ const&) { return val; // Can't use lexical_cast here, sometimes it fails.... } diff --git a/boost/math/constants/constants.hpp b/boost/math/constants/constants.hpp index e9381adeb6..59873c60d8 100644 --- a/boost/math/constants/constants.hpp +++ b/boost/math/constants/constants.hpp @@ -269,9 +269,9 @@ namespace boost{ namespace math \ \ /* Now the namespace specific versions: */ \ - } namespace float_constants{ static const float name = BOOST_JOIN(x, F); }\ - namespace double_constants{ static const double name = x; } \ - namespace long_double_constants{ static const long double name = BOOST_JOIN(x, L); }\ + } namespace float_constants{ BOOST_STATIC_CONSTEXPR float name = BOOST_JOIN(x, F); }\ + namespace double_constants{ BOOST_STATIC_CONSTEXPR double name = x; } \ + namespace long_double_constants{ BOOST_STATIC_CONSTEXPR long double name = BOOST_JOIN(x, L); }\ namespace constants{ BOOST_DEFINE_MATH_CONSTANT(half, 5.000000000000000000000000000000000000e-01, "5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01") diff --git a/boost/math/cstdfloat/cstdfloat_types.hpp b/boost/math/cstdfloat/cstdfloat_types.hpp index 3ffcce21db..9e0d00bf99 100644 --- a/boost/math/cstdfloat/cstdfloat_types.hpp +++ b/boost/math/cstdfloat/cstdfloat_types.hpp @@ -356,6 +356,15 @@ #undef BOOST_CSTDFLOAT_FLOAT_32_MAX #endif +#if (defined(__SGI_STL_PORT) || defined(_STLPORT_VERSION)) && defined(__SUNPRO_CC) +#undef BOOST_CSTDFLOAT_HAS_FLOAT80_NATIVE_TYPE +#define BOOST_CSTDFLOAT_HAS_FLOAT80_NATIVE_TYPE 0 +#undef BOOST_CSTDFLOAT_HAS_FLOAT128_NATIVE_TYPE +#define BOOST_CSTDFLOAT_HAS_FLOAT128_NATIVE_TYPE 0 +#undef BOOST_CSTDFLOAT_MAXIMUM_AVAILABLE_WIDTH +#define BOOST_CSTDFLOAT_MAXIMUM_AVAILABLE_WIDTH 64 +#endif + #if(BOOST_CSTDFLOAT_HAS_FLOAT64_NATIVE_TYPE == 1) typedef BOOST_CSTDFLOAT_FLOAT64_NATIVE_TYPE float64_t; typedef boost::float64_t float_fast64_t; diff --git a/boost/math/distributions.hpp b/boost/math/distributions.hpp index cbcc0ace88..300f2312ce 100644 --- a/boost/math/distributions.hpp +++ b/boost/math/distributions.hpp @@ -12,6 +12,7 @@ #ifndef BOOST_MATH_DISTRIBUTIONS_HPP #define BOOST_MATH_DISTRIBUTIONS_HPP +#include <boost/math/distributions/arcsine.hpp> #include <boost/math/distributions/bernoulli.hpp> #include <boost/math/distributions/beta.hpp> #include <boost/math/distributions/binomial.hpp> diff --git a/boost/math/distributions/arcsine.hpp b/boost/math/distributions/arcsine.hpp new file mode 100644 index 0000000000..8aad5b2d0b --- /dev/null +++ b/boost/math/distributions/arcsine.hpp @@ -0,0 +1,535 @@ +// boost/math/distributions/arcsine.hpp + +// Copyright John Maddock 2014. +// Copyright Paul A. Bristow 2014. + +// Use, modification and distribution are subject to 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) + +// http://en.wikipedia.org/wiki/arcsine_distribution + +// The arcsine Distribution is a continuous probability distribution. +// http://en.wikipedia.org/wiki/Arcsine_distribution +// http://www.wolframalpha.com/input/?i=ArcSinDistribution + +// Standard arcsine distribution is a special case of beta distribution with both a & b = one half, +// and 0 <= x <= 1. + +// It is generalized to include any bounded support a <= x <= b from 0 <= x <= 1 +// by Wolfram and Wikipedia, +// but using location and scale parameters by +// Virtual Laboratories in Probability and Statistics http://www.math.uah.edu/stat/index.html +// http://www.math.uah.edu/stat/special/Arcsine.html +// The end-point version is simpler and more obvious, so we implement that. +// TODO Perhaps provide location and scale functions? + + +#ifndef BOOST_MATH_DIST_ARCSINE_HPP +#define BOOST_MATH_DIST_ARCSINE_HPP + +#include <boost/math/distributions/fwd.hpp> +#include <boost/math/distributions/complement.hpp> // complements. +#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks. +#include <boost/math/constants/constants.hpp> + +#include <boost/math/special_functions/fpclassify.hpp> // isnan. + +#if defined (BOOST_MSVC) +# pragma warning(push) +# pragma warning(disable: 4702) // Unreachable code, +// in domain_error_imp in error_handling. +#endif + +#include <utility> +#include <exception> // For std::domain_error. + +namespace boost +{ + namespace math + { + namespace arcsine_detail + { + // Common error checking routines for arcsine distribution functions: + // Duplicating for x_min and x_max provides specific error messages. + template <class RealType, class Policy> + inline bool check_x_min(const char* function, const RealType& x, RealType* result, const Policy& pol) + { + if (!(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "x_min argument is %1%, but must be finite !", x, pol); + return false; + } + return true; + } // bool check_x_min + + template <class RealType, class Policy> + inline bool check_x_max(const char* function, const RealType& x, RealType* result, const Policy& pol) + { + if (!(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "x_max argument is %1%, but must be finite !", x, pol); + return false; + } + return true; + } // bool check_x_max + + + template <class RealType, class Policy> + inline bool check_x_minmax(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) + { // Check x_min < x_max + if (x_min >= x_max) + { + std::string msg = "x_max argument is %1%, but must be > x_min = " + lexical_cast<std::string>(x_min) + "!"; + *result = policies::raise_domain_error<RealType>( + function, + msg.c_str(), x_max, pol); + // "x_max argument is %1%, but must be > x_min !", x_max, pol); + // "x_max argument is %1%, but must be > x_min %2!", x_max, x_min, pol); would be better. + // But would require replication of all helpers functions in /policies/error_handling.hpp for two values, + // as well as two value versions of raise_error, raise_domain_error and do_format ... + // so use slightly hacky lexical_cast to string instead. + return false; + } + return true; + } // bool check_x_minmax + + template <class RealType, class Policy> + inline bool check_prob(const char* function, const RealType& p, RealType* result, const Policy& pol) + { + if ((p < 0) || (p > 1) || !(boost::math::isfinite)(p)) + { + *result = policies::raise_domain_error<RealType>( + function, + "Probability argument is %1%, but must be >= 0 and <= 1 !", p, pol); + return false; + } + return true; + } // bool check_prob + + template <class RealType, class Policy> + inline bool check_x(const char* function, const RealType& x_min, const RealType& x_max, const RealType& x, RealType* result, const Policy& pol) + { // Check x finite and x_min < x < x_max. + if (!(boost::math::isfinite)(x)) + { + *result = policies::raise_domain_error<RealType>( + function, + "x argument is %1%, but must be finite !", x, pol); + return false; + } + if ((x < x_min) || (x > x_max)) + { + // std::cout << x_min << ' ' << x << x_max << std::endl; + *result = policies::raise_domain_error<RealType>( + function, + "x argument is %1%, but must be x_min < x < x_max !", x, pol); + // For example: + // Error in function boost::math::pdf(arcsine_distribution<double> const&, double) : x argument is -1.01, but must be x_min < x < x_max ! + // TODO Perhaps show values of x_min and x_max? + return false; + } + return true; + } // bool check_x + + template <class RealType, class Policy> + inline bool check_dist(const char* function, const RealType& x_min, const RealType& x_max, RealType* result, const Policy& pol) + { // Check both x_min and x_max finite, and x_min < x_max. + return check_x_min(function, x_min, result, pol) + && check_x_max(function, x_max, result, pol) + && check_x_minmax(function, x_min, x_max, result, pol); + } // bool check_dist + + template <class RealType, class Policy> + inline bool check_dist_and_x(const char* function, const RealType& x_min, const RealType& x_max, RealType x, RealType* result, const Policy& pol) + { + return check_dist(function, x_min, x_max, result, pol) + && arcsine_detail::check_x(function, x_min, x_max, x, result, pol); + } // bool check_dist_and_x + + template <class RealType, class Policy> + inline bool check_dist_and_prob(const char* function, const RealType& x_min, const RealType& x_max, RealType p, RealType* result, const Policy& pol) + { + return check_dist(function, x_min, x_max, result, pol) + && check_prob(function, p, result, pol); + } // bool check_dist_and_prob + + } // namespace arcsine_detail + + template <class RealType = double, class Policy = policies::policy<> > + class arcsine_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + arcsine_distribution(RealType x_min = 0, RealType x_max = 1) : m_x_min(x_min), m_x_max(x_max) + { // Default beta (alpha = beta = 0.5) is standard arcsine with x_min = 0, x_max = 1. + // Generalized to allow x_min and x_max to be specified. + RealType result; + arcsine_detail::check_dist( + "boost::math::arcsine_distribution<%1%>::arcsine_distribution", + m_x_min, + m_x_max, + &result, Policy()); + } // arcsine_distribution constructor. + // Accessor functions: + RealType x_min() const + { + return m_x_min; + } + RealType x_max() const + { + return m_x_max; + } + + private: + RealType m_x_min; // Two x min and x max parameters of the arcsine distribution. + RealType m_x_max; + }; // template <class RealType, class Policy> class arcsine_distribution + + // Convenient typedef to construct double version. + typedef arcsine_distribution<double> arcsine; + + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const arcsine_distribution<RealType, Policy>& dist) + { // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(dist.x_min()), static_cast<RealType>(dist.x_max())); + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const arcsine_distribution<RealType, Policy>& dist) + { // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + return std::pair<RealType, RealType>(static_cast<RealType>(dist.x_min()), static_cast<RealType>(dist.x_max())); + } + + template <class RealType, class Policy> + inline RealType mean(const arcsine_distribution<RealType, Policy>& dist) + { // Mean of arcsine distribution . + RealType result; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + if (false == arcsine_detail::check_dist( + "boost::math::mean(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + return (x_min + x_max) / 2; + } // mean + + template <class RealType, class Policy> + inline RealType variance(const arcsine_distribution<RealType, Policy>& dist) + { // Variance of standard arcsine distribution = (1-0)/8 = 0.125. + RealType result; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + if (false == arcsine_detail::check_dist( + "boost::math::variance(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + return (x_max - x_min) * (x_max - x_min) / 8; + } // variance + + template <class RealType, class Policy> + inline RealType mode(const arcsine_distribution<RealType, Policy>& /* dist */) + { //There are always [*two] values for the mode, at ['x_min] and at ['x_max], default 0 and 1, + // so instead we raise the exception domain_error. + return policies::raise_domain_error<RealType>( + "boost::math::mode(arcsine_distribution<%1%>&)", + "The arcsine distribution has two modes at x_min and x_max: " + "so the return value is %1%.", + std::numeric_limits<RealType>::quiet_NaN(), Policy()); + } // mode + + template <class RealType, class Policy> + inline RealType median(const arcsine_distribution<RealType, Policy>& dist) + { // Median of arcsine distribution (a + b) / 2 == mean. + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + RealType result; + if (false == arcsine_detail::check_dist( + "boost::math::median(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + return (x_min + x_max) / 2; + } + + template <class RealType, class Policy> + inline RealType skewness(const arcsine_distribution<RealType, Policy>& dist) + { + RealType result; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + if (false == arcsine_detail::check_dist( + "boost::math::skewness(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + return 0; + } // skewness + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const arcsine_distribution<RealType, Policy>& dist) + { + RealType result; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + if (false == arcsine_detail::check_dist( + "boost::math::kurtosis_excess(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + result = -3; + return result / 2; + } // kurtosis_excess + + template <class RealType, class Policy> + inline RealType kurtosis(const arcsine_distribution<RealType, Policy>& dist) + { + RealType result; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + if (false == arcsine_detail::check_dist( + "boost::math::kurtosis(arcsine_distribution<%1%> const&, %1% )", + x_min, + x_max, + &result, Policy()) + ) + { + return result; + } + + return 3 + kurtosis_excess(dist); + } // kurtosis + + template <class RealType, class Policy> + inline RealType pdf(const arcsine_distribution<RealType, Policy>& dist, const RealType& xx) + { // Probability Density/Mass Function arcsine. + BOOST_FPU_EXCEPTION_GUARD + BOOST_MATH_STD_USING // For ADL of std functions. + + static const char* function = "boost::math::pdf(arcsine_distribution<%1%> const&, %1%)"; + + RealType lo = dist.x_min(); + RealType hi = dist.x_max(); + RealType x = xx; + + // Argument checks: + RealType result = 0; + if (false == arcsine_detail::check_dist_and_x( + function, + lo, hi, x, + &result, Policy())) + { + return result; + } + using boost::math::constants::pi; + result = static_cast<RealType>(1) / (pi<RealType>() * sqrt((x - lo) * (hi - x))); + return result; + } // pdf + + template <class RealType, class Policy> + inline RealType cdf(const arcsine_distribution<RealType, Policy>& dist, const RealType& x) + { // Cumulative Distribution Function arcsine. + BOOST_MATH_STD_USING // For ADL of std functions. + + static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; + + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + // Argument checks: + RealType result = 0; + if (false == arcsine_detail::check_dist_and_x( + function, + x_min, x_max, x, + &result, Policy())) + { + return result; + } + // Special cases: + if (x == x_min) + { + return 0; + } + else if (x == x_max) + { + return 1; + } + using boost::math::constants::pi; + result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>(); + return result; + } // arcsine cdf + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<arcsine_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function arcsine. + BOOST_MATH_STD_USING // For ADL of std functions. + static const char* function = "boost::math::cdf(arcsine_distribution<%1%> const&, %1%)"; + + RealType x = c.param; + arcsine_distribution<RealType, Policy> const& dist = c.dist; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + + // Argument checks: + RealType result = 0; + if (false == arcsine_detail::check_dist_and_x( + function, + x_min, x_max, x, + &result, Policy())) + { + return result; + } + if (x == x_min) + { + return 0; + } + else if (x == x_max) + { + return 1; + } + using boost::math::constants::pi; + // Naive version x = 1 - x; + // result = static_cast<RealType>(2) * asin(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>(); + // is less accurate, so use acos instead of asin for complement. + result = static_cast<RealType>(2) * acos(sqrt((x - x_min) / (x_max - x_min))) / pi<RealType>(); + return result; + } // arcine ccdf + + template <class RealType, class Policy> + inline RealType quantile(const arcsine_distribution<RealType, Policy>& dist, const RealType& p) + { + // Quantile or Percent Point arcsine function or + // Inverse Cumulative probability distribution function CDF. + // Return x (0 <= x <= 1), + // for a given probability p (0 <= p <= 1). + // These functions take a probability as an argument + // and return a value such that the probability that a random variable x + // will be less than or equal to that value + // is whatever probability you supplied as an argument. + BOOST_MATH_STD_USING // For ADL of std functions. + + using boost::math::constants::half_pi; + + static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; + + RealType result = 0; // of argument checks: + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + if (false == arcsine_detail::check_dist_and_prob( + function, + x_min, x_max, p, + &result, Policy())) + { + return result; + } + // Special cases: + if (p == 0) + { + return 0; + } + if (p == 1) + { + return 1; + } + + RealType sin2hpip = sin(half_pi<RealType>() * p); + RealType sin2hpip2 = sin2hpip * sin2hpip; + result = -x_min * sin2hpip2 + x_min + x_max * sin2hpip2; + + return result; + } // quantile + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<arcsine_distribution<RealType, Policy>, RealType>& c) + { + // Complement Quantile or Percent Point arcsine function. + // Return the number of expected x for a given + // complement of the probability q. + BOOST_MATH_STD_USING // For ADL of std functions. + + using boost::math::constants::half_pi; + static const char* function = "boost::math::quantile(arcsine_distribution<%1%> const&, %1%)"; + + // Error checks: + RealType q = c.param; + const arcsine_distribution<RealType, Policy>& dist = c.dist; + RealType result = 0; + RealType x_min = dist.x_min(); + RealType x_max = dist.x_max(); + if (false == arcsine_detail::check_dist_and_prob( + function, + x_min, + x_max, + q, + &result, Policy())) + { + return result; + } + // Special cases: + if (q == 1) + { + return 0; + } + if (q == 0) + { + return 1; + } + // Naive RealType p = 1 - q; result = sin(half_pi<RealType>() * p); loses accuracy, so use a cos alternative instead. + //result = cos(half_pi<RealType>() * q); // for arcsine(0,1) + //result = result * result; + // For generalized arcsine: + RealType cos2hpip = cos(half_pi<RealType>() * q); + RealType cos2hpip2 = cos2hpip * cos2hpip; + result = -x_min * cos2hpip2 + x_min + x_max * cos2hpip2; + + return result; + } // Quantile Complement + + } // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#if defined (BOOST_MSVC) +# pragma warning(pop) +#endif + +#endif // BOOST_MATH_DIST_ARCSINE_HPP diff --git a/boost/math/distributions/extreme_value.hpp b/boost/math/distributions/extreme_value.hpp index ef9fbe817d..4c6b710cbd 100644 --- a/boost/math/distributions/extreme_value.hpp +++ b/boost/math/distributions/extreme_value.hpp @@ -108,7 +108,10 @@ inline RealType pdf(const extreme_value_distribution<RealType, Policy>& dist, co return result; if(0 == detail::check_x(function, x, &result, Policy())) return result; - result = exp((a-x)/b) * exp(-exp((a-x)/b)) / b; + RealType e = (a - x) / b; + if(e < tools::log_max_value<RealType>()) + result = exp(e) * exp(-exp(e)) / b; + // else.... result *must* be zero since exp(e) is infinite... return result; } // pdf diff --git a/boost/math/distributions/fwd.hpp b/boost/math/distributions/fwd.hpp index d2b50ad010..5b212c8ec6 100644 --- a/boost/math/distributions/fwd.hpp +++ b/boost/math/distributions/fwd.hpp @@ -1,6 +1,6 @@ // fwd.hpp Forward declarations of Boost.Math distributions. -// Copyright Paul A. Bristow 2007, 2010, 2012. +// Copyright Paul A. Bristow 2007, 2010, 2012, 2014. // Copyright John Maddock 2007. // Use, modification and distribution are subject to the @@ -11,11 +11,14 @@ #ifndef BOOST_MATH_DISTRIBUTIONS_FWD_HPP #define BOOST_MATH_DISTRIBUTIONS_FWD_HPP -// 31 distributions at Boost 1.52 +// 33 distributions at Boost 1.9.1 after adding hyperexpon and arcsine namespace boost{ namespace math{ template <class RealType, class Policy> +class arcsine_distribution; + +template <class RealType, class Policy> class bernoulli_distribution; template <class RealType, class Policy> @@ -114,6 +117,7 @@ class weibull_distribution; }} // namespaces #define BOOST_MATH_DECLARE_DISTRIBUTIONS(Type, Policy)\ + typedef boost::math::arcsine_distribution<Type, Policy> arcsine;\ typedef boost::math::bernoulli_distribution<Type, Policy> bernoulli;\ typedef boost::math::beta_distribution<Type, Policy> beta;\ typedef boost::math::binomial_distribution<Type, Policy> binomial;\ diff --git a/boost/math/policies/policy.hpp b/boost/math/policies/policy.hpp index 49068a6ed5..71309fa116 100644 --- a/boost/math/policies/policy.hpp +++ b/boost/math/policies/policy.hpp @@ -850,6 +850,11 @@ inline int digits(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) typedef mpl::bool_< std::numeric_limits<T>::is_specialized > tag_type; return detail::digits_imp<T, Policy>(tag_type()); } +template <class T, class Policy> +inline int digits_base10(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) +{ + return boost::math::policies::digits<T, Policy>() * 301 / 1000L; +} template <class Policy> inline unsigned long get_max_series_iterations() diff --git a/boost/math/special_functions.hpp b/boost/math/special_functions.hpp index aa7019a1eb..52e412ed8b 100644 --- a/boost/math/special_functions.hpp +++ b/boost/math/special_functions.hpp @@ -1,4 +1,4 @@ -// Copyright John Maddock 2006, 2007, 2012. +// Copyright John Maddock 2006, 2007, 2012, 2014. // Copyright Paul A. Bristow 2006, 2007, 2012 // Use, modification and distribution are subject to the @@ -27,10 +27,14 @@ #include <boost/math/special_functions/ellint_1.hpp> #include <boost/math/special_functions/ellint_2.hpp> #include <boost/math/special_functions/ellint_3.hpp> +#include <boost/math/special_functions/ellint_d.hpp> +#include <boost/math/special_functions/jacobi_zeta.hpp> +#include <boost/math/special_functions/heuman_lambda.hpp> #include <boost/math/special_functions/ellint_rc.hpp> #include <boost/math/special_functions/ellint_rd.hpp> #include <boost/math/special_functions/ellint_rf.hpp> #include <boost/math/special_functions/ellint_rj.hpp> +#include <boost/math/special_functions/ellint_rg.hpp> #include <boost/math/special_functions/erf.hpp> #include <boost/math/special_functions/expint.hpp> #include <boost/math/special_functions/expm1.hpp> @@ -47,6 +51,7 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/special_functions/next.hpp> #include <boost/math/special_functions/owens_t.hpp> +#include <boost/math/special_functions/polygamma.hpp> #include <boost/math/special_functions/powm1.hpp> #include <boost/math/special_functions/sign.hpp> #include <boost/math/special_functions/sin_pi.hpp> diff --git a/boost/math/special_functions/airy.hpp b/boost/math/special_functions/airy.hpp index 82167dc5f0..e84a705fe8 100644 --- a/boost/math/special_functions/airy.hpp +++ b/boost/math/special_functions/airy.hpp @@ -164,7 +164,7 @@ T airy_ai_zero_imp(int m, const Policy& pol) if(m < 0) { return policies::raise_domain_error<T>("boost::math::airy_ai_zero<%1%>(%1%, int)", - "Requested the %1%'th zero, but the rank must be 1 or more !", m, pol); + "Requested the %1%'th zero, but the rank must be 1 or more !", static_cast<T>(m), pol); } // Handle case when the zero'th zero is requested. @@ -217,7 +217,7 @@ T airy_bi_zero_imp(int m, const Policy& pol) if(m < 0) { return policies::raise_domain_error<T>("boost::math::airy_bi_zero<%1%>(%1%, int)", - "Requested the %1%'th zero, but the rank must 1 or more !", m, pol); + "Requested the %1%'th zero, but the rank must 1 or more !", static_cast<T>(m), pol); } // Handle case when the zero'th zero is requested. diff --git a/boost/math/special_functions/bernoulli.hpp b/boost/math/special_functions/bernoulli.hpp index 2c2ccd5236..5f2b7e1422 100644 --- a/boost/math/special_functions/bernoulli.hpp +++ b/boost/math/special_functions/bernoulli.hpp @@ -63,7 +63,7 @@ inline T bernoulli_b2n(const int i, const Policy &pol) if(i < 0) return policies::raise_domain_error<T>("boost::math::bernoulli_b2n<%1%>", "Index should be >= 0 but got %1%", T(i), pol); - T result; + T result = 0; // The = 0 is just to silence compiler warings :-( boost::math::detail::bernoulli_number_imp<T>(&result, static_cast<std::size_t>(i), 1u, pol, tag_type()); return result; } diff --git a/boost/math/special_functions/bessel.hpp b/boost/math/special_functions/bessel.hpp index 438b763ab7..eab723d8ba 100644 --- a/boost/math/special_functions/bessel.hpp +++ b/boost/math/special_functions/bessel.hpp @@ -194,7 +194,7 @@ T cyl_bessel_i_imp(T v, T x, const Policy& pol) } if(x == 0) { - return (v == 0) ? 1 : 0; + return (v == 0) ? static_cast<T>(1) : static_cast<T>(0); } if(v == 0.5f) { @@ -386,7 +386,7 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol) if(m < 0) { // Zeros of Jv(x) with negative rank are not defined and requesting one raises a domain error. - return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol); + return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol); } // Get the absolute value of the order. @@ -402,7 +402,7 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol) if(order_is_zero) { // The zero'th zero of J0(x) is not defined and requesting it raises a domain error. - return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of J0, but the rank must be > 0 !", m, pol); + return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of J0, but the rank must be > 0 !", static_cast<T>(m), pol); } // The zero'th zero of Jv(x) for v < 0 is not defined @@ -410,7 +410,7 @@ inline T cyl_bessel_j_zero_imp(T v, int m, const Policy& pol) if(order_is_negative && (!order_is_integer)) { // For non-integer, negative order, requesting the zero'th zero raises a domain error. - return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Jv for negative, non-integer order, but the rank must be > 0 !", m, pol); + return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Jv for negative, non-integer order, but the rank must be > 0 !", static_cast<T>(m), pol); } // The zero'th zero does exist and its value is zero. @@ -462,7 +462,7 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol) // Handle negative rank. if(m < 0) { - return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", m, pol); + return policies::raise_domain_error<T>(function, "Requested the %1%'th zero, but the rank must be positive !", static_cast<T>(m), pol); } const T half_epsilon(boost::math::tools::epsilon<T>() / 2U); @@ -488,7 +488,7 @@ inline T cyl_neumann_zero_imp(T v, int m, const Policy& pol) if((m == 0) && (!order_is_negative_half_integer)) { // For non-integer, negative order, requesting the zero'th zero raises a domain error. - return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", m, pol); + return policies::raise_domain_error<T>(function, "Requested the %1%'th zero of Yv for negative, non-half-integer order, but the rank must be > 0 !", static_cast<T>(m), pol); } // For negative half-integers, use the corresponding diff --git a/boost/math/special_functions/beta.hpp b/boost/math/special_functions/beta.hpp index 8486b824f2..98d8f7fa80 100644 --- a/boost/math/special_functions/beta.hpp +++ b/boost/math/special_functions/beta.hpp @@ -13,6 +13,7 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/tools/config.hpp> #include <boost/math/special_functions/gamma.hpp> +#include <boost/math/special_functions/binomial.hpp> #include <boost/math/special_functions/factorials.hpp> #include <boost/math/special_functions/erf.hpp> #include <boost/math/special_functions/log1p.hpp> @@ -845,15 +846,54 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const Policy& po // complement of the binomial distribution cdf and use this finite sum. // template <class T> -inline T binomial_ccdf(T n, T k, T x, T y) +T binomial_ccdf(T n, T k, T x, T y) { BOOST_MATH_STD_USING // ADL of std names + T result = pow(x, n); - T term = result; - for(unsigned i = itrunc(T(n - 1)); i > k; --i) + + if(result > tools::min_value<T>()) { - term *= ((i + 1) * y) / ((n - i) * x) ; - result += term; + T term = result; + for(unsigned i = itrunc(T(n - 1)); i > k; --i) + { + term *= ((i + 1) * y) / ((n - i) * x); + result += term; + } + } + else + { + // First term underflows so we need to start at the mode of the + // distribution and work outwards: + int start = itrunc(n * x); + if(start <= k + 1) + start = itrunc(k + 2); + result = pow(x, start) * pow(y, n - start) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(start)); + if(result == 0) + { + // OK, starting slightly above the mode didn't work, + // we'll have to sum the terms the old fashioned way: + for(unsigned i = start - 1; i > k; --i) + { + result += pow(x, (int)i) * pow(y, n - i) * boost::math::binomial_coefficient<T>(itrunc(n), itrunc(i)); + } + } + else + { + T term = result; + T start_term = result; + for(unsigned i = start - 1; i > k; --i) + { + term *= ((i + 1) * y) / ((n - i) * x); + result += term; + } + term = start_term; + for(unsigned i = start + 1; i <= n; ++i) + { + term *= (n - i + 1) * x / (i * y); + result += term; + } + } } return result; @@ -1208,9 +1248,9 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de } else if(normalised) { - // the formula here for the non-normalised case is tricky to figure + // The formula here for the non-normalised case is tricky to figure // out (for me!!), and requires two pochhammer calculations rather - // than one, so leave it for now.... + // than one, so leave it for now and only use this in the normalized case.... int n = itrunc(T(floor(b)), pol); T bbar = b - n; if(bbar <= 0) @@ -1221,8 +1261,7 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de fract = ibeta_a_step(bbar, a, y, x, n, pol, normalised, static_cast<T*>(0)); fract += ibeta_a_step(a, bbar, x, y, 20, pol, normalised, static_cast<T*>(0)); if(invert) - fract -= (normalised ? 1 : boost::math::beta(a, b, pol)); - //fract = ibeta_series(a+20, bbar, x, fract, l, normalised, p_derivative, y); + fract -= 1; // Note this line would need changing if we ever enable this branch in non-normalized case fract = beta_small_b_large_a_series(T(a+20), bbar, x, y, fract, T(1), pol, normalised); if(invert) { diff --git a/boost/math/special_functions/binomial.hpp b/boost/math/special_functions/binomial.hpp index 8fa2e26d84..9a24fc15bb 100644 --- a/boost/math/special_functions/binomial.hpp +++ b/boost/math/special_functions/binomial.hpp @@ -27,12 +27,12 @@ T binomial_coefficient(unsigned n, unsigned k, const Policy& pol) return policies::raise_domain_error<T>( function, "The binomial coefficient is undefined for k > n, but got k = %1%.", - k, pol); + static_cast<T>(k), pol); T result; if((k == 0) || (k == n)) - return 1; + return static_cast<T>(1); if((k == 1) || (k == n-1)) - return n; + return static_cast<T>(n); if(n <= max_factorial<T>::value) { diff --git a/boost/math/special_functions/cos_pi.hpp b/boost/math/special_functions/cos_pi.hpp index 18d06c00df..669a2c87ae 100644 --- a/boost/math/special_functions/cos_pi.hpp +++ b/boost/math/special_functions/cos_pi.hpp @@ -25,10 +25,10 @@ T cos_pi_imp(T x, const Policy& pol) BOOST_MATH_STD_USING // ADL of std names // cos of pi*x: bool invert = false; - if(fabs(x) < 0.5) + if(fabs(x) < 0.25) return cos(constants::pi<T>() * x); - if(x < 1) + if(x < 0) { x = -x; } @@ -44,17 +44,30 @@ T cos_pi_imp(T x, const Policy& pol) if(rem == 0.5f) return 0; - rem = cos(constants::pi<T>() * rem); + if(rem > 0.25f) + { + rem = 0.5f - rem; + rem = sin(constants::pi<T>() * rem); + } + else + rem = cos(constants::pi<T>() * rem); return invert ? T(-rem) : rem; } } // namespace detail template <class T, class Policy> -inline typename tools::promote_args<T>::type cos_pi(T x, const Policy& pol) +inline typename tools::promote_args<T>::type cos_pi(T x, const Policy&) { typedef typename tools::promote_args<T>::type result_type; - return boost::math::detail::cos_pi_imp<result_type>(x, pol); + typedef typename policies::evaluation<result_type, Policy>::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast<result_type, forwarding_policy>(boost::math::detail::cos_pi_imp<value_type>(x, forwarding_policy()), "cos_pi"); } template <class T> diff --git a/boost/math/special_functions/detail/erf_inv.hpp b/boost/math/special_functions/detail/erf_inv.hpp index 77aa72fc26..35072d5155 100644 --- a/boost/math/special_functions/detail/erf_inv.hpp +++ b/boost/math/special_functions/detail/erf_inv.hpp @@ -103,7 +103,7 @@ T erf_inv_imp(const T& p, const T& q, const Policy&, const boost::mpl::int_<64>* BOOST_MATH_BIG_CONSTANT(T, 64, 1.72114765761200282724) }; T g = sqrt(-2 * log(q)); - T xs = q - 0.25; + T xs = q - 0.25f; T r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs); result = g / (Y + r); } @@ -156,7 +156,7 @@ T erf_inv_imp(const T& p, const T& q, const Policy&, const boost::mpl::int_<64>* BOOST_MATH_BIG_CONSTANT(T, 64, 0.152264338295331783612), BOOST_MATH_BIG_CONSTANT(T, 64, 0.01105924229346489121) }; - T xs = x - 1.125; + T xs = x - 1.125f; T R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs); result = Y * x + R * x; } diff --git a/boost/math/special_functions/detail/fp_traits.hpp b/boost/math/special_functions/detail/fp_traits.hpp index 63ebf11ae0..09dc5169aa 100644 --- a/boost/math/special_functions/detail/fp_traits.hpp +++ b/boost/math/special_functions/detail/fp_traits.hpp @@ -555,7 +555,8 @@ struct select_native<long double> && !defined(__SGI_STL_PORT) && !defined(_STLPORT_VERSION)\ && !defined(__FAST_MATH__)\ && !defined(BOOST_MATH_DISABLE_STD_FPCLASSIFY)\ - && !defined(BOOST_INTEL) + && !defined(BOOST_INTEL)\ + && !defined(sun) # define BOOST_MATH_USE_STD_FPCLASSIFY #endif diff --git a/boost/math/special_functions/detail/polygamma.hpp b/boost/math/special_functions/detail/polygamma.hpp new file mode 100644 index 0000000000..4ef503bf7c --- /dev/null +++ b/boost/math/special_functions/detail/polygamma.hpp @@ -0,0 +1,538 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2013 Nikhar Agrawal +// Copyright 2013 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2013 Paul Bristow +// 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) + +#ifndef _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ + #define _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ + + #include <cmath> + #include <limits> + #include <boost/cstdint.hpp> + #include <boost/math/policies/policy.hpp> + #include <boost/math/special_functions/bernoulli.hpp> + #include <boost/math/special_functions/trunc.hpp> + #include <boost/math/special_functions/zeta.hpp> + #include <boost/math/special_functions/digamma.hpp> + #include <boost/math/special_functions/sin_pi.hpp> + #include <boost/math/special_functions/cos_pi.hpp> + #include <boost/math/special_functions/pow.hpp> + #include <boost/mpl/if.hpp> + #include <boost/mpl/int.hpp> + #include <boost/static_assert.hpp> + #include <boost/type_traits/is_convertible.hpp> + + namespace boost { namespace math { namespace detail{ + + template<class T, class Policy> + T polygamma_atinfinityplus(const int n, const T& x, const Policy& pol, const char* function) // for large values of x such as for x> 400 + { + // See http://functions.wolfram.com/GammaBetaErf/PolyGamma2/06/02/0001/ + BOOST_MATH_STD_USING + // + // sum == current value of accumulated sum. + // term == value of current term to be added to sum. + // part_term == value of current term excluding the Bernoulli number part + // + if(n + x == x) + { + // x is crazy large, just concentrate on the first part of the expression and use logs: + if(n == 1) return 1 / x; + T nlx = n * log(x); + if((nlx < tools::log_max_value<T>()) && (n < max_factorial<T>::value)) + return ((n & 1) ? 1 : -1) * boost::math::factorial<T>(n - 1) * pow(x, -n); + else + return ((n & 1) ? 1 : -1) * exp(boost::math::lgamma(T(n), pol) - n * log(x)); + } + T term, sum, part_term; + T x_squared = x * x; + // + // Start by setting part_term to: + // + // (n-1)! / x^(n+1) + // + // which is common to both the first term of the series (with k = 1) + // and to the leading part. + // We can then get to the leading term by: + // + // part_term * (n + 2 * x) / 2 + // + // and to the first term in the series + // (excluding the Bernoulli number) by: + // + // part_term n * (n + 1) / (2x) + // + // If either the factorial would overflow, + // or the power term underflows, this just gets set to 0 and then we + // know that we have to use logs for the initial terms: + // + part_term = ((n > boost::math::max_factorial<T>::value) && (T(n) * n > tools::log_max_value<T>())) + ? T(0) : static_cast<T>(boost::math::factorial<T>(n - 1, pol) * pow(x, -n - 1)); + if(part_term == 0) + { + // Either n is very large, or the power term underflows, + // set the initial values of part_term, term and sum via logs: + part_term = boost::math::lgamma(n, pol) - (n + 1) * log(x); + sum = exp(part_term + log(n + 2 * x) - boost::math::constants::ln_two<T>()); + part_term += log(T(n) * (n + 1)) - boost::math::constants::ln_two<T>() - log(x); + part_term = exp(part_term); + } + else + { + sum = part_term * (n + 2 * x) / 2; + part_term *= (T(n) * (n + 1)) / 2; + part_term /= x; + } + // + // If the leading term is 0, so is the result: + // + if(sum == 0) + return sum; + + for(unsigned k = 1;;) + { + term = part_term * boost::math::bernoulli_b2n<T>(k, pol); + sum += term; + // + // Normal termination condition: + // + if(fabs(term / sum) < tools::epsilon<T>()) + break; + // + // Increment our counter, and move part_term on to the next value: + // + ++k; + part_term *= T(n + 2 * k - 2) * (n - 1 + 2 * k); + part_term /= (2 * k - 1) * 2 * k; + part_term /= x_squared; + // + // Emergency get out termination condition: + // + if(k > policies::get_max_series_iterations<Policy>()) + { + return policies::raise_evaluation_error(function, "Series did not converge, closest value was %1%", sum, pol); + } + } + + if((n - 1) & 1) + sum = -sum; + + return sum; + } + + template<class T, class Policy> + T polygamma_attransitionplus(const int n, const T& x, const Policy& pol, const char* function) + { + // See: http://functions.wolfram.com/GammaBetaErf/PolyGamma2/16/01/01/0017/ + + // Use N = (0.4 * digits) + (4 * n) for target value for x: + BOOST_MATH_STD_USING + const int d4d = static_cast<int>(0.4F * policies::digits_base10<T, Policy>()); + const int N = d4d + (4 * n); + const int m = n; + const int iter = N - itrunc(x); + + if(iter > (int)policies::get_max_series_iterations<Policy>()) + return policies::raise_evaluation_error<T>(function, ("Exceeded maximum series evaluations evaluating at n = " + boost::lexical_cast<std::string>(n) + " and x = %1%").c_str(), x, pol); + + const int minus_m_minus_one = -m - 1; + + T z(x); + T sum0(0); + T z_plus_k_pow_minus_m_minus_one(0); + + // Forward recursion to larger x, need to check for overflow first though: + if(log(z + iter) * minus_m_minus_one > -tools::log_max_value<T>()) + { + for(int k = 1; k <= iter; ++k) + { + z_plus_k_pow_minus_m_minus_one = pow(z, minus_m_minus_one); + sum0 += z_plus_k_pow_minus_m_minus_one; + z += 1; + } + sum0 *= boost::math::factorial<T>(n); + } + else + { + for(int k = 1; k <= iter; ++k) + { + T log_term = log(z) * minus_m_minus_one + boost::math::lgamma(T(n + 1), pol); + sum0 += exp(log_term); + z += 1; + } + } + if((n - 1) & 1) + sum0 = -sum0; + + return sum0 + polygamma_atinfinityplus(n, z, pol, function); + } + + template <class T, class Policy> + T polygamma_nearzero(int n, T x, const Policy& pol, const char* function) + { + BOOST_MATH_STD_USING + // + // If we take this expansion for polygamma: http://functions.wolfram.com/06.15.06.0003.02 + // and substitute in this expression for polygamma(n, 1): http://functions.wolfram.com/06.15.03.0009.01 + // we get an alternating series for polygamma when x is small in terms of zeta functions of + // integer arguments (which are easy to evaluate, at least when the integer is even). + // + // In order to avoid spurious overflow, save the n! term for later, and rescale at the end: + // + T scale = boost::math::factorial<T>(n, pol); + // + // "factorial_part" contains everything except the zeta function + // evaluations in each term: + // + T factorial_part = 1; + // + // "prefix" is what we'll be adding the accumulated sum to, it will + // be n! / z^(n+1), but since we're scaling by n! it's just + // 1 / z^(n+1) for now: + // + T prefix = pow(x, n + 1); + if(prefix == 0) + return boost::math::policies::raise_overflow_error<T>(function, 0, pol); + prefix = 1 / prefix; + // + // First term in the series is necessarily < zeta(2) < 2, so + // ignore the sum if it will have no effect on the result anyway: + // + if(prefix > 2 / policies::get_epsilon<T, Policy>()) + return ((n & 1) ? 1 : -1) * + (tools::max_value<T>() / prefix < scale ? policies::raise_overflow_error<T>(function, 0, pol) : prefix * scale); + // + // As this is an alternating series we could accelerate it using + // "Convergence Acceleration of Alternating Series", + // Henri Cohen, Fernando Rodriguez Villegas, and Don Zagier, Experimental Mathematics, 1999. + // In practice however, it appears not to make any difference to the number of terms + // required except in some edge cases which are filtered out anyway before we get here. + // + T sum = prefix; + for(unsigned k = 0;;) + { + // Get the k'th term: + T term = factorial_part * boost::math::zeta(T(k + n + 1), pol); + sum += term; + // Termination condition: + if(fabs(term) < fabs(sum * boost::math::policies::get_epsilon<T, Policy>())) + break; + // + // Move on k and factorial_part: + // + ++k; + factorial_part *= (-x * (n + k)) / k; + // + // Last chance exit: + // + if(k > policies::get_max_series_iterations<Policy>()) + return policies::raise_evaluation_error<T>(function, "Series did not converge, best value is %1%", sum, pol); + } + // + // We need to multiply by the scale, at each stage checking for oveflow: + // + if(boost::math::tools::max_value<T>() / scale < sum) + return boost::math::policies::raise_overflow_error<T>(function, 0, pol); + sum *= scale; + return n & 1 ? sum : -sum; + } + + // + // Helper function which figures out which slot our coefficient is in + // given an angle multiplier for the cosine term of power: + // + template <class Table> + typename Table::value_type::reference dereference_table(Table& table, unsigned row, unsigned power) + { + return table[row][power / 2]; + } + + + + template <class T, class Policy> + T poly_cot_pi(int n, T x, T xc, const Policy& pol, const char* function) + { + BOOST_MATH_STD_USING + // Return n'th derivative of cot(pi*x) at x, these are simply + // tabulated for up to n = 9, beyond that it is possible to + // calculate coefficients as follows: + // + // The general form of each derivative is: + // + // pi^n * SUM{k=0, n} C[k,n] * cos^k(pi * x) * csc^(n+1)(pi * x) + // + // With constant C[0,1] = -1 and all other C[k,n] = 0; + // Then for each k < n+1: + // C[k-1, n+1] -= k * C[k, n]; + // C[k+1, n+1] += (k-n-1) * C[k, n]; + // + // Note that there are many different ways of representing this derivative thanks to + // the many trigomonetric identies available. In particular, the sum of powers of + // cosines could be replaced by a sum of cosine multiple angles, and indeed if you + // plug the derivative into Mathematica this is the form it will give. The two + // forms are related via the Chebeshev polynomials of the first kind and + // T_n(cos(x)) = cos(n x). The polynomial form has the great advantage that + // all the cosine terms are zero at half integer arguments - right where this + // function has it's minumum - thus avoiding cancellation error in this region. + // + // And finally, since every other term in the polynomials is zero, we can save + // space by only storing the non-zero terms. This greatly complexifies + // subscripting the tables in the calculation, but halves the storage space + // (and complexity for that matter). + // + T s = fabs(x) < fabs(xc) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(xc, pol); + T c = boost::math::cos_pi(x, pol); + switch(n) + { + case 1: + return -constants::pi<T, Policy>() / (s * s); + case 2: + { + return 2 * constants::pi<T, Policy>() * constants::pi<T, Policy>() * c / boost::math::pow<3>(s, pol); + } + case 3: + { + int P[] = { -2, -4 }; + return boost::math::pow<3>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<4>(s, pol); + } + case 4: + { + int P[] = { 16, 8 }; + return boost::math::pow<4>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<5>(s, pol); + } + case 5: + { + int P[] = { -16, -88, -16 }; + return boost::math::pow<5>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<6>(s, pol); + } + case 6: + { + int P[] = { 272, 416, 32 }; + return boost::math::pow<6>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<7>(s, pol); + } + case 7: + { + int P[] = { -272, -2880, -1824, -64 }; + return boost::math::pow<7>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<8>(s, pol); + } + case 8: + { + int P[] = { 7936, 24576, 7680, 128 }; + return boost::math::pow<8>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<9>(s, pol); + } + case 9: + { + int P[] = { -7936, -137216, -185856, -31616, -256 }; + return boost::math::pow<9>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<10>(s, pol); + } + case 10: + { + int P[] = { 353792, 1841152, 1304832, 128512, 512 }; + return boost::math::pow<10>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<11>(s, pol); + } + case 11: + { + int P[] = { -353792, -9061376, -21253376, -8728576, -518656, -1024}; + return boost::math::pow<11>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<12>(s, pol); + } + case 12: + { + int P[] = { 22368256, 175627264, 222398464, 56520704, 2084864, 2048 }; + return boost::math::pow<12>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<13>(s, pol); + } +#ifndef BOOST_NO_LONG_LONG + case 13: + { + long long P[] = { -22368256LL, -795300864LL, -2868264960LL, -2174832640LL, -357888000LL, -8361984LL, -4096 }; + return boost::math::pow<13>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<14>(s, pol); + } + case 14: + { + long long P[] = { 1903757312LL, 21016670208LL, 41731645440LL, 20261765120LL, 2230947840LL, 33497088LL, 8192 }; + return boost::math::pow<14>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<15>(s, pol); + } + case 15: + { + long long P[] = { -1903757312LL, -89702612992LL, -460858269696LL, -559148810240LL, -182172651520LL, -13754155008LL, -134094848LL, -16384 }; + return boost::math::pow<15>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<16>(s, pol); + } + case 16: + { + long long P[] = { 209865342976LL, 3099269660672LL, 8885192097792LL, 7048869314560LL, 1594922762240LL, 84134068224LL, 536608768LL, 32768 }; + return boost::math::pow<16>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<17>(s, pol); + } + case 17: + { + long long P[] = { -209865342976LL, -12655654469632LL, -87815735738368LL, -155964390375424LL, -84842998005760LL, -13684856848384LL, -511780323328LL, -2146926592LL, -65536 }; + return boost::math::pow<17>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<18>(s, pol); + } + case 18: + { + long long P[] = { 29088885112832LL, 553753414467584LL, 2165206642589696LL, 2550316668551168LL, 985278548541440LL, 115620218667008LL, 3100738912256LL, 8588754944LL, 131072 }; + return boost::math::pow<18>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<19>(s, pol); + } + case 19: + { + long long P[] = { -29088885112832LL, -2184860175433728LL, -19686087844429824LL, -48165109676113920LL, -39471306959486976LL, -11124607890751488LL, -965271355195392LL, -18733264797696LL, -34357248000LL, -262144 }; + return boost::math::pow<19>(constants::pi<T, Policy>(), pol) * tools::evaluate_even_polynomial(P, c) / boost::math::pow<20>(s, pol); + } + case 20: + { + long long P[] = { 4951498053124096LL, 118071834535526400LL, 603968063567560704LL, 990081991141490688LL, 584901762421358592LL, 122829335169859584LL, 7984436548730880LL, 112949304754176LL, 137433710592LL, 524288 }; + return boost::math::pow<20>(constants::pi<T, Policy>(), pol) * c * tools::evaluate_even_polynomial(P, c) / boost::math::pow<21>(s, pol); + } +#endif + } + + // + // We'll have to compute the coefficients up to n, + // complexity is O(n^2) which we don't worry about for now + // as the values are computed once and then cached. + // However, if the final evaluation would have too many + // terms just bail out right away: + // + if((unsigned)n / 2u > policies::get_max_series_iterations<Policy>()) + return policies::raise_evaluation_error<T>(function, "The value of n is so large that we're unable to compute the result in reasonable time, best guess is %1%", 0, pol); +#ifdef BOOST_HAS_THREADS + static boost::detail::lightweight_mutex m; + boost::detail::lightweight_mutex::scoped_lock l(m); +#endif + static std::vector<std::vector<T> > table(1, std::vector<T>(1, T(-1))); + + int index = n - 1; + + if(index >= (int)table.size()) + { + for(int i = (int)table.size() - 1; i < index; ++i) + { + int offset = i & 1; // 1 if the first cos power is 0, otherwise 0. + int sin_order = i + 2; // order of the sin term + int max_cos_order = sin_order - 1; // largest order of the polynomial of cos terms + int max_columns = (max_cos_order - offset) / 2; // How many entries there are in the current row. + int next_offset = offset ? 0 : 1; + int next_max_columns = (max_cos_order + 1 - next_offset) / 2; // How many entries there will be in the next row + table.push_back(std::vector<T>(next_max_columns + 1, T(0))); + + for(int column = 0; column <= max_columns; ++column) + { + int cos_order = 2 * column + offset; // order of the cosine term in entry "column" + BOOST_ASSERT(column < (int)table[i].size()); + BOOST_ASSERT((cos_order + 1) / 2 < (int)table[i + 1].size()); + table[i + 1][(cos_order + 1) / 2] += ((cos_order - sin_order) * table[i][column]) / (sin_order - 1); + if(cos_order) + table[i + 1][(cos_order - 1) / 2] += (-cos_order * table[i][column]) / (sin_order - 1); + } + } + + } + T sum = boost::math::tools::evaluate_even_polynomial(&table[index][0], c, table[index].size()); + if(index & 1) + sum *= c; // First coeffient is order 1, and really an odd polynomial. + if(sum == 0) + return sum; + // + // The remaining terms are computed using logs since the powers and factorials + // get real large real quick: + // + T power_terms = n * log(boost::math::constants::pi<T>()); + if(s == 0) + return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol); + power_terms -= log(fabs(s)) * (n + 1); + power_terms += boost::math::lgamma(T(n)); + power_terms += log(fabs(sum)); + + if(power_terms > boost::math::tools::log_max_value<T>()) + return sum * boost::math::policies::raise_overflow_error<T>(function, 0, pol); + + return exp(power_terms) * ((s < 0) && ((n + 1) & 1) ? -1 : 1) * boost::math::sign(sum); + } + + template <class T, class Policy> + struct polygamma_initializer + { + struct init + { + init() + { + // Forces initialization of our table of coefficients and mutex: + boost::math::polygamma(30, T(-2.5f), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } + }; + + template <class T, class Policy> + const typename polygamma_initializer<T, Policy>::init polygamma_initializer<T, Policy>::initializer; + + template<class T, class Policy> + inline T polygamma_imp(const int n, T x, const Policy &pol) + { + BOOST_MATH_STD_USING + static const char* function = "boost::math::polygamma<%1%>(int, %1%)"; + polygamma_initializer<T, Policy>::initializer.force_instantiate(); + if(n < 0) + return policies::raise_domain_error<T>(function, "Order must be >= 0, but got %1%", static_cast<T>(n), pol); + if(x < 0) + { + if(floor(x) == x) + { + // + // Result is infinity if x is odd, and a pole error if x is even. + // + if(lltrunc(x) & 1) + return policies::raise_overflow_error<T>(function, 0, pol); + else + return policies::raise_pole_error<T>(function, "Evaluation at negative integer %1%", x, pol); + } + T z = 1 - x; + T result = polygamma_imp(n, z, pol) + constants::pi<T, Policy>() * poly_cot_pi(n, z, x, pol, function); + return n & 1 ? T(-result) : result; + } + // + // Limit for use of small-x-series is chosen + // so that the series doesn't go too divergent + // in the first few terms. Ordinarily this + // would mean setting the limit to ~ 1 / n, + // but we can tolerate a small amount of divergence: + // + T small_x_limit = std::min(T(T(5) / n), T(0.25f)); + if(x < small_x_limit) + { + return polygamma_nearzero(n, x, pol, function); + } + else if(x > 0.4F * policies::digits_base10<T, Policy>() + 4.0f * n) + { + return polygamma_atinfinityplus(n, x, pol, function); + } + else if(x == 1) + { + return (n & 1 ? 1 : -1) * boost::math::factorial<T>(n, pol) * boost::math::zeta(T(n + 1), pol); + } + else if(x == 0.5f) + { + T result = (n & 1 ? 1 : -1) * boost::math::factorial<T>(n, pol) * boost::math::zeta(T(n + 1), pol); + if(fabs(result) >= ldexp(tools::max_value<T>(), -n - 1)) + return boost::math::sign(result) * policies::raise_overflow_error<T>(function, 0, pol); + result *= ldexp(T(1), n + 1) - 1; + return result; + } + else + { + return polygamma_attransitionplus(n, x, pol, function); + } + } + +} } } // namespace boost::math::detail + +#endif // _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ + diff --git a/boost/math/special_functions/digamma.hpp b/boost/math/special_functions/digamma.hpp index 785cd75c5e..718eaf9278 100644 --- a/boost/math/special_functions/digamma.hpp +++ b/boost/math/special_functions/digamma.hpp @@ -12,6 +12,7 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/tools/rational.hpp> +#include <boost/math/tools/series.hpp> #include <boost/math/tools/promotion.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/constants/constants.hpp> @@ -27,7 +28,8 @@ namespace detail{ // inline unsigned digamma_large_lim(const mpl::int_<0>*) { return 20; } - +inline unsigned digamma_large_lim(const mpl::int_<113>*) +{ return 20; } inline unsigned digamma_large_lim(const void*) { return 10; } // @@ -41,7 +43,7 @@ inline unsigned digamma_large_lim(const void*) // This first one gives 34-digit precision for x >= 20: // template <class T> -inline T digamma_imp_large(T x, const mpl::int_<0>*) +inline T digamma_imp_large(T x, const mpl::int_<113>*) { BOOST_MATH_STD_USING // ADL of std functions. static const T P[] = { @@ -141,12 +143,47 @@ inline T digamma_imp_large(T x, const mpl::int_<24>*) return result; } // +// Fully generic asymptotic expansion in terms of Bernoulli numbers, see: +// http://functions.wolfram.com/06.14.06.0012.01 +// +template <class T> +struct digamma_series_func +{ +private: + int k; + T xx; + T term; +public: + digamma_series_func(T x) : k(1), xx(x * x), term(1 / (x * x)) {} + T operator()() + { + T result = term * boost::math::bernoulli_b2n<T>(k) / (2 * k); + term /= xx; + ++k; + return result; + } + typedef T result_type; +}; + +template <class T, class Policy> +inline T digamma_imp_large(T x, const Policy& pol, const mpl::int_<0>*) +{ + BOOST_MATH_STD_USING + digamma_series_func<T> s(x); + T result = log(x) - 1 / (2 * x); + boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); + result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, -result); + result = -result; + policies::check_series_iterations<T>("boost::math::digamma<%1%>(%1%)", max_iter, pol); + return result; +} +// // Now follow rational approximations over the range [1,2]. // // 35-digit precision: // template <class T> -T digamma_imp_1_2(T x, const mpl::int_<0>*) +T digamma_imp_1_2(T x, const mpl::int_<113>*) { // // Now the approximation, we use the form: @@ -325,16 +362,16 @@ inline T digamma_imp_1_2(T x, const mpl::int_<24>*) static const T root = 1532632.0f / 1048576; static const T root_minor = static_cast<T>(0.3700660185912626595423257213284682051735604e-6L); static const T P[] = { - 0.25479851023250261e0, - -0.44981331915268368e0, - -0.43916936919946835e0, - -0.61041765350579073e-1 + 0.25479851023250261e0f, + -0.44981331915268368e0f, + -0.43916936919946835e0f, + -0.61041765350579073e-1f }; static const T Q[] = { 0.1e1, - 0.15890202430554952e1, - 0.65341249856146947e0, - 0.63851690523355715e-1 + 0.15890202430554952e1f, + 0.65341249856146947e0f, + 0.63851690523355715e-1f }; T g = x - root; g -= root_minor; @@ -410,6 +447,98 @@ T digamma_imp(T x, const Tag* t, const Policy& pol) return result; } +template <class T, class Policy> +T digamma_imp(T x, const mpl::int_<0>* t, const Policy& pol) +{ + // + // This handles reflection of negative arguments, and all our + // error handling, then forwards to the T-specific approximation. + // + BOOST_MATH_STD_USING // ADL of std functions. + + T result = 0; + // + // Check for negative arguments and use reflection: + // + if(x <= -1) + { + // Reflect: + x = 1 - x; + // Argument reduction for tan: + T remainder = x - floor(x); + // Shift to negative if > 0.5: + if(remainder > 0.5) + { + remainder -= 1; + } + // + // check for evaluation at a negative pole: + // + if(remainder == 0) + { + return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, (1 - x), pol); + } + result = constants::pi<T>() / tan(constants::pi<T>() * remainder); + } + if(x == 0) + return policies::raise_pole_error<T>("boost::math::digamma<%1%>(%1%)", 0, x, pol); + // + // If we're above the lower-limit for the + // asymptotic expansion then use it, the + // limit is a linear interpolation with + // limit = 10 at 50 bit precision and + // limit = 250 at 1000 bit precision. + // + T lim = 10 + (tools::digits<T>() - 50) * 240 / 950; + T two_x = ldexp(x, 1); + if(x >= lim) + { + result += digamma_imp_large(x, pol, t); + } + else if(floor(x) == x) + { + // + // Special case for integer arguments, see + // http://functions.wolfram.com/06.14.03.0001.01 + // + result = -constants::euler<T, Policy>(); + T val = 1; + while(val < x) + { + result += 1 / val; + val += 1; + } + } + else if(floor(two_x) == two_x) + { + // + // Special case for half integer arguments, see: + // http://functions.wolfram.com/06.14.03.0007.01 + // + result = -2 * constants::ln_two<T, Policy>() - constants::euler<T, Policy>(); + int n = itrunc(x); + if(n) + { + for(int k = 1; k < n; ++k) + result += 1 / T(k); + for(int k = n; k <= 2 * n - 1; ++k) + result += 2 / T(k); + } + } + else + { + // + // Rescale so we can use the asymptotic expansion: + // + while(x < lim) + { + result -= 1 / x; + x += 1; + } + result += digamma_imp_large(x, pol, t); + } + return result; +} // // Initializer: ensure all our constants are initialized prior to the first call of main: // @@ -420,9 +549,15 @@ struct digamma_initializer { init() { + typedef typename policies::precision<T, Policy>::type precision_type; + do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>()); + } + void do_init(const mpl::true_&) + { boost::math::digamma(T(1.5), Policy()); boost::math::digamma(T(500), Policy()); } + void do_init(const mpl::false_&){} void force_instantiate()const{} }; static const init initializer; @@ -439,7 +574,7 @@ const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Polic template <class T, class Policy> inline typename tools::promote_args<T>::type - digamma(T x, const Policy& pol) + digamma(T x, const Policy&) { typedef typename tools::promote_args<T>::type result_type; typedef typename policies::evaluation<result_type, Policy>::type value_type; @@ -447,7 +582,7 @@ inline typename tools::promote_args<T>::type typedef typename mpl::if_< mpl::or_< mpl::less_equal<precision_type, mpl::int_<0> >, - mpl::greater<precision_type, mpl::int_<64> > + mpl::greater<precision_type, mpl::int_<114> > >, mpl::int_<0>, typename mpl::if_< @@ -456,17 +591,28 @@ inline typename tools::promote_args<T>::type typename mpl::if_< mpl::less<precision_type, mpl::int_<54> >, mpl::int_<53>, - mpl::int_<64> + typename mpl::if_< + mpl::less<precision_type, mpl::int_<65> >, + mpl::int_<64>, + mpl::int_<113> + >::type >::type >::type >::type tag_type; + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + // Force initialization of constants: - detail::digamma_initializer<result_type, Policy>::force_instantiate(); + detail::digamma_initializer<value_type, forwarding_policy>::force_instantiate(); return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp( static_cast<value_type>(x), - static_cast<const tag_type*>(0), pol), "boost::math::digamma<%1%>(%1%)"); + static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::digamma<%1%>(%1%)"); } template <class T> diff --git a/boost/math/special_functions/ellint_1.hpp b/boost/math/special_functions/ellint_1.hpp index da16bc6f26..62a0bf3fde 100644 --- a/boost/math/special_functions/ellint_1.hpp +++ b/boost/math/special_functions/ellint_1.hpp @@ -103,10 +103,22 @@ T ellint_f_imp(T phi, T k, const Policy& pol) BOOST_MATH_INSTRUMENT_VARIABLE(rphi); } T sinp = sin(rphi); + sinp *= sinp; T cosp = cos(rphi); + cosp *= cosp; + T c = 1 / sinp; BOOST_MATH_INSTRUMENT_VARIABLE(sinp); BOOST_MATH_INSTRUMENT_VARIABLE(cosp); - result = s * sinp * ellint_rf_imp(T(cosp * cosp), T(1 - k * k * sinp * sinp), T(1), pol); + if(sinp > tools::min_value<T>()) + { + // + // Use http://dlmf.nist.gov/19.25#E5, note that + // c-1 simplifies to cot^2(rphi) which avoid cancellation: + // + result = rphi == 0 ? static_cast<T>(0) : static_cast<T>(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol)); + } + else + result = s * sin(rphi); BOOST_MATH_INSTRUMENT_VARIABLE(result); if(m != 0) { diff --git a/boost/math/special_functions/ellint_2.hpp b/boost/math/special_functions/ellint_2.hpp index 72caf3eb11..f4f65cc11d 100644 --- a/boost/math/special_functions/ellint_2.hpp +++ b/boost/math/special_functions/ellint_2.hpp @@ -21,6 +21,7 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/special_functions/ellint_rf.hpp> #include <boost/math/special_functions/ellint_rd.hpp> +#include <boost/math/special_functions/ellint_rg.hpp> #include <boost/math/constants/constants.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/tools/workaround.hpp> @@ -67,6 +68,14 @@ T ellint_e_imp(T phi, T k, const Policy& pol) // just return the second part of the duplication formula: result = 2 * phi * ellint_e_imp(k, pol) / constants::pi<T>(); } + else if(k == 0) + { + return invert ? T(-phi) : phi; + } + else if(fabs(k) == 1) + { + return invert ? T(-sin(phi)) : sin(phi); + } else { // Carlson's algorithm works only for |phi| <= pi/2, @@ -87,11 +96,30 @@ T ellint_e_imp(T phi, T k, const Policy& pol) } T sinp = sin(rphi); T cosp = cos(rphi); - T x = cosp * cosp; - T t = k * k * sinp * sinp; - T y = 1 - t; - T z = 1; - result = s * sinp * (ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3); + T c = 1 / (sinp * sinp); + T cm1 = cosp * cosp / (sinp * sinp); // c - 1 + T k2 = k * k; + if(k2 > 1) + { + return policies::raise_domain_error<T>("boost::math::ellint_2<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); + } + else if(rphi == 0) + { + result = 0; + } + else if(sinp * sinp < tools::min_value<T>()) + { + T x = cosp * cosp; + T t = k * k * sinp * sinp; + T y = 1 - t; + T z = 1; + result = s * sinp * (ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3); + } + else + { + // http://dlmf.nist.gov/19.25#E10 + result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); + } if(m != 0) result += m * ellint_e_imp(k, pol); } @@ -119,7 +147,7 @@ T ellint_e_imp(T k, const Policy& pol) T t = k * k; T y = 1 - t; T z = 1; - T value = ellint_rf_imp(x, y, z, pol) - t * ellint_rd_imp(x, y, z, pol) / 3; + T value = 2 * ellint_rg_imp(x, y, z, pol); return value; } diff --git a/boost/math/special_functions/ellint_3.hpp b/boost/math/special_functions/ellint_3.hpp index ac7e68c17f..9ab0f8317a 100644 --- a/boost/math/special_functions/ellint_3.hpp +++ b/boost/math/special_functions/ellint_3.hpp @@ -24,6 +24,7 @@ #include <boost/math/special_functions/ellint_1.hpp> #include <boost/math/special_functions/ellint_2.hpp> #include <boost/math/special_functions/log1p.hpp> +#include <boost/math/special_functions/atanh.hpp> #include <boost/math/constants/constants.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/tools/workaround.hpp> @@ -43,176 +44,231 @@ T ellint_pi_imp(T v, T k, T vc, const Policy& pol); template <typename T, typename Policy> T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol) { - // Note vc = 1-v presumably without cancellation error. - T value, x, y, z, p, t; - - BOOST_MATH_STD_USING - using namespace boost::math::tools; - using namespace boost::math::constants; - - static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)"; - - if (abs(k) > 1) - { - return policies::raise_domain_error<T>(function, - "Got k = %1%, function requires |k| <= 1", k, pol); - } - - T sphi = sin(fabs(phi)); - - if(v > 1 / (sphi * sphi)) - { - // Complex result is a domain error: - return policies::raise_domain_error<T>(function, - "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol); - } - - // Special cases first: - if(v == 0) - { - // A&S 17.7.18 & 19 - return (k == 0) ? phi : ellint_f_imp(phi, k, pol); - } - if(phi == constants::pi<T>() / 2) - { - // Have to filter this case out before the next - // special case, otherwise we might get an infinity from - // tan(phi). - // Also note that since we can't represent PI/2 exactly - // in a T, this is a bit of a guess as to the users true - // intent... - // - return ellint_pi_imp(v, k, vc, pol); - } - if(k == 0) - { - // A&S 17.7.20: - if(v < 1) - { - T vcr = sqrt(vc); - return atan(vcr * tan(phi)) / vcr; - } - else if(v == 1) - { - return tan(phi); - } - else - { - // v > 1: - T vcr = sqrt(-vc); - T arg = vcr * tan(phi); - return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr); - } - } - - if(v < 0) - { - // - // If we don't shift to 0 <= v <= 1 we get - // cancellation errors later on. Use - // A&S 17.7.15/16 to shift to v > 0: - // - T k2 = k * k; - T N = (k2 - v) / (1 - v); - T Nm1 = (1 - k2) / (1 - v); - T p2 = sqrt(-v * (k2 - v) / (1 - v)); - T delta = sqrt(1 - k2 * sphi * sphi); - T result = ellint_pi_imp(N, phi, k, Nm1, pol); - - result *= sqrt(Nm1 * (1 - k2 / N)); - result += ellint_f_imp(phi, k, pol) * k2 / p2; - result += atan((p2/2) * sin(2 * phi) / delta); - result /= sqrt((1 - v) * (1 - k2 / v)); - return result; - } -#if 0 // disabled but retained for future reference: see below. - if(v > 1) - { - // - // If v > 1 we can use the identity in A&S 17.7.7/8 - // to shift to 0 <= v <= 1. Unfortunately this - // identity appears only to function correctly when - // 0 <= phi <= pi/2, but it's when phi is outside that - // range that we really need it: That's when - // Carlson's formula fails, and what's more the periodicity - // reduction used below on phi doesn't work when v > 1. - // - // So we're stuck... the code is archived here in case - // some bright spart can figure out the fix. - // - T k2 = k * k; - T N = k2 / v; - T Nm1 = (v - k2) / v; - T p1 = sqrt((-vc) * (1 - k2 / v)); - T delta = sqrt(1 - k2 * sphi * sphi); - // - // These next two terms have a large amount of cancellation - // so it's not clear if this relation is useable even if - // the issues with phi > pi/2 can be fixed: - // - T result = -ellint_pi_imp(N, phi, k, Nm1); - result += ellint_f_imp(phi, k); - // - // This log term gives the complex result when - // n > 1/sin^2(phi) - // However that case is dealt with as an error above, - // so we should always get a real result here: - // - result += log((delta + p1 * tan(phi)) / (delta - p1 * tan(phi))) / (2 * p1); - return result; - } -#endif - - // Carlson's algorithm works only for |phi| <= pi/2, - // use the integrand's periodicity to normalize phi - // - // Xiaogang's original code used a cast to long long here - // but that fails if T has more digits than a long long, - // so rewritten to use fmod instead: - // - if(fabs(phi) > 1 / tools::epsilon<T>()) - { - if(v > 1) - return policies::raise_domain_error<T>( + // Note vc = 1-v presumably without cancellation error. + BOOST_MATH_STD_USING + + static const char* function = "boost::math::ellint_3<%1%>(%1%,%1%,%1%)"; + + if(abs(k) > 1) + { + return policies::raise_domain_error<T>(function, + "Got k = %1%, function requires |k| <= 1", k, pol); + } + + T sphi = sin(fabs(phi)); + T result = 0; + + if(v > 1 / (sphi * sphi)) + { + // Complex result is a domain error: + return policies::raise_domain_error<T>(function, + "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol); + } + + // Special cases first: + if(v == 0) + { + // A&S 17.7.18 & 19 + return (k == 0) ? phi : ellint_f_imp(phi, k, pol); + } + if(v == 1) + { + // http://functions.wolfram.com/08.06.03.0008.01 + T m = k * k; + result = sqrt(1 - m * sphi * sphi) * tan(phi) - ellint_e_imp(phi, k, pol); + result /= 1 - m; + result += ellint_f_imp(phi, k, pol); + return result; + } + if(phi == constants::half_pi<T>()) + { + // Have to filter this case out before the next + // special case, otherwise we might get an infinity from + // tan(phi). + // Also note that since we can't represent PI/2 exactly + // in a T, this is a bit of a guess as to the users true + // intent... + // + return ellint_pi_imp(v, k, vc, pol); + } + if((phi > constants::half_pi<T>()) || (phi < 0)) + { + // Carlson's algorithm works only for |phi| <= pi/2, + // use the integrand's periodicity to normalize phi + // + // Xiaogang's original code used a cast to long long here + // but that fails if T has more digits than a long long, + // so rewritten to use fmod instead: + // + // See http://functions.wolfram.com/08.06.16.0002.01 + // + if(fabs(phi) > 1 / tools::epsilon<T>()) + { + if(v > 1) + return policies::raise_domain_error<T>( function, "Got v = %1%, but this is only supported for 0 <= phi <= pi/2", v, pol); - // - // Phi is so large that phi%pi is necessarily zero (or garbage), - // just return the second part of the duplication formula: - // - value = 2 * fabs(phi) * ellint_pi_imp(v, k, vc, pol) / constants::pi<T>(); - } - else - { - T rphi = boost::math::tools::fmod_workaround(T(fabs(phi)), T(constants::half_pi<T>())); - T m = boost::math::round((fabs(phi) - rphi) / constants::half_pi<T>()); - int sign = 1; - if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5) - { - m += 1; - sign = -1; - rphi = constants::half_pi<T>() - rphi; - } - T sinp = sin(rphi); - T cosp = cos(rphi); - x = cosp * cosp; - t = sinp * sinp; - y = 1 - k * k * t; - z = 1; - if(v * t < 0.5) - p = 1 - v * t; - else - p = x + vc * t; - value = sign * sinp * (ellint_rf_imp(x, y, z, pol) + v * t * ellint_rj_imp(x, y, z, p, pol) / 3); - if((m > 0) && (vc > 0)) - value += m * ellint_pi_imp(v, k, vc, pol); - } - - if (phi < 0) - { - value = -value; // odd function - } - return value; + // + // Phi is so large that phi%pi is necessarily zero (or garbage), + // just return the second part of the duplication formula: + // + result = 2 * fabs(phi) * ellint_pi_imp(v, k, vc, pol) / constants::pi<T>(); + } + else + { + T rphi = boost::math::tools::fmod_workaround(T(fabs(phi)), T(constants::half_pi<T>())); + T m = boost::math::round((fabs(phi) - rphi) / constants::half_pi<T>()); + int sign = 1; + if((m != 0) && (k >= 1)) + { + return policies::raise_domain_error<T>(function, "Got k=1 and phi=%1% but the result is complex in that domain", phi, pol); + } + if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5) + { + m += 1; + sign = -1; + rphi = constants::half_pi<T>() - rphi; + } + result = sign * ellint_pi_imp(v, rphi, k, vc, pol); + if((m > 0) && (vc > 0)) + result += m * ellint_pi_imp(v, k, vc, pol); + } + return phi < 0 ? -result : result; + } + if(k == 0) + { + // A&S 17.7.20: + if(v < 1) + { + T vcr = sqrt(vc); + return atan(vcr * tan(phi)) / vcr; + } + else if(v == 1) + { + return tan(phi); + } + else + { + // v > 1: + T vcr = sqrt(-vc); + T arg = vcr * tan(phi); + return (boost::math::log1p(arg, pol) - boost::math::log1p(-arg, pol)) / (2 * vcr); + } + } + if(v < 0) + { + // + // If we don't shift to 0 <= v <= 1 we get + // cancellation errors later on. Use + // A&S 17.7.15/16 to shift to v > 0. + // + // Mathematica simplifies the expressions + // given in A&S as follows (with thanks to + // Rocco Romeo for figuring these out!): + // + // V = (k2 - n)/(1 - n) + // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[Sqrt[(1 - V)*(1 - k2 / V)] / Sqrt[((1 - n)*(1 - k2 / n))]]] + // Result: ((-1 + k2) n) / ((-1 + n) (-k2 + n)) + // + // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[k2 / (Sqrt[-n*(k2 - n) / (1 - n)] * Sqrt[(1 - n)*(1 - k2 / n)])]] + // Result : k2 / (k2 - n) + // + // Assuming[(k2 >= 0 && k2 <= 1) && n < 0, FullSimplify[Sqrt[1 / ((1 - n)*(1 - k2 / n))]]] + // Result : Sqrt[n / ((k2 - n) (-1 + n))] + // + T k2 = k * k; + T N = (k2 - v) / (1 - v); + T Nm1 = (1 - k2) / (1 - v); + T p2 = -v * N; + T t; + if(p2 <= tools::min_value<T>()) + p2 = sqrt(-v) * sqrt(N); + else + p2 = sqrt(p2); + T delta = sqrt(1 - k2 * sphi * sphi); + if(N > k2) + { + result = ellint_pi_imp(N, phi, k, Nm1, pol); + result *= v / (v - 1); + result *= (k2 - 1) / (v - k2); + } + + if(k != 0) + { + t = ellint_f_imp(phi, k, pol); + t *= k2 / (k2 - v); + result += t; + } + t = v / ((k2 - v) * (v - 1)); + if(t > tools::min_value<T>()) + { + result += atan((p2 / 2) * sin(2 * phi) / delta) * sqrt(t); + } + else + { + result += atan((p2 / 2) * sin(2 * phi) / delta) * sqrt(fabs(1 / (k2 - v))) * sqrt(fabs(v / (v - 1))); + } + return result; + } + if(k == 1) + { + // See http://functions.wolfram.com/08.06.03.0013.01 + result = sqrt(v) * atanh(sqrt(v) * sin(phi)) - log(1 / cos(phi) + tan(phi)); + result /= v - 1; + return result; + } +#if 0 // disabled but retained for future reference: see below. + if(v > 1) + { + // + // If v > 1 we can use the identity in A&S 17.7.7/8 + // to shift to 0 <= v <= 1. In contrast to previous + // revisions of this header, this identity does now work + // but appears not to produce better error rates in + // practice. Archived here for future reference... + // + T k2 = k * k; + T N = k2 / v; + T Nm1 = (v - k2) / v; + T p1 = sqrt((-vc) * (1 - k2 / v)); + T delta = sqrt(1 - k2 * sphi * sphi); + // + // These next two terms have a large amount of cancellation + // so it's not clear if this relation is useable even if + // the issues with phi > pi/2 can be fixed: + // + result = -ellint_pi_imp(N, phi, k, Nm1, pol); + result += ellint_f_imp(phi, k, pol); + // + // This log term gives the complex result when + // n > 1/sin^2(phi) + // However that case is dealt with as an error above, + // so we should always get a real result here: + // + result += log((delta + p1 * tan(phi)) / (delta - p1 * tan(phi))) / (2 * p1); + return result; + } +#endif + // + // Carlson's algorithm works only for |phi| <= pi/2, + // by the time we get here phi should already have been + // normalised above. + // + BOOST_ASSERT(fabs(phi) < constants::half_pi<T>()); + BOOST_ASSERT(phi >= 0); + T x, y, z, p, t; + T cosp = cos(phi); + x = cosp * cosp; + t = sphi * sphi; + y = 1 - k * k * t; + z = 1; + if(v * t < 0.5) + p = 1 - v * t; + else + p = x + vc * t; + result = sphi * (ellint_rf_imp(x, y, z, pol) + v * t * ellint_rj_imp(x, y, z, p, pol) / 3); + + return result; } // Complete elliptic integral (Legendre form) of the third kind @@ -244,16 +300,16 @@ T ellint_pi_imp(T v, T k, T vc, const Policy& pol) if(v < 0) { + // Apply A&S 17.7.17: T k2 = k * k; T N = (k2 - v) / (1 - v); T Nm1 = (1 - k2) / (1 - v); - T p2 = sqrt(-v * (k2 - v) / (1 - v)); - - T result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol); - - result *= sqrt(Nm1 * (1 - k2 / N)); - result += ellint_k_imp(k, pol) * k2 / p2; - result /= sqrt((1 - v) * (1 - k2 / v)); + T result = 0; + result = boost::math::detail::ellint_pi_imp(N, k, Nm1, pol); + // This next part is split in two to avoid spurious over/underflow: + result *= -v / (1 - v); + result *= (1 - k2) / (k2 - v); + result += ellint_k_imp(k, pol) * k2 / (k2 - v); return result; } diff --git a/boost/math/special_functions/ellint_d.hpp b/boost/math/special_functions/ellint_d.hpp new file mode 100644 index 0000000000..5bd065d6e3 --- /dev/null +++ b/boost/math/special_functions/ellint_d.hpp @@ -0,0 +1,180 @@ +// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 John Maddock +// Use, modification and distribution are subject to 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) +// +// History: +// XZ wrote the original of this file as part of the Google +// Summer of Code 2006. JM modified it to fit into the +// Boost.Math conceptual framework better, and to ensure +// that the code continues to work no matter how many digits +// type T has. + +#ifndef BOOST_MATH_ELLINT_D_HPP +#define BOOST_MATH_ELLINT_D_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/special_functions/ellint_rf.hpp> +#include <boost/math/special_functions/ellint_rd.hpp> +#include <boost/math/special_functions/ellint_rg.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/tools/workaround.hpp> +#include <boost/math/special_functions/round.hpp> + +// Elliptic integrals (complete and incomplete) of the second kind +// Carlson, Numerische Mathematik, vol 33, 1 (1979) + +namespace boost { namespace math { + +template <class T1, class T2, class Policy> +typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol); + +namespace detail{ + +template <typename T, typename Policy> +T ellint_d_imp(T k, const Policy& pol); + +// Elliptic integral (Legendre form) of the second kind +template <typename T, typename Policy> +T ellint_d_imp(T phi, T k, const Policy& pol) +{ + BOOST_MATH_STD_USING + using namespace boost::math::tools; + using namespace boost::math::constants; + + bool invert = false; + if(phi < 0) + { + phi = fabs(phi); + invert = true; + } + + T result; + + if(phi >= tools::max_value<T>()) + { + // Need to handle infinity as a special case: + result = policies::raise_overflow_error<T>("boost::math::ellint_e<%1%>(%1%,%1%)", 0, pol); + } + else if(phi > 1 / tools::epsilon<T>()) + { + // Phi is so large that phi%pi is necessarily zero (or garbage), + // just return the second part of the duplication formula: + result = 2 * phi * ellint_d_imp(k, pol) / constants::pi<T>(); + } + else + { + // Carlson's algorithm works only for |phi| <= pi/2, + // use the integrand's periodicity to normalize phi + // + T rphi = boost::math::tools::fmod_workaround(phi, T(constants::half_pi<T>())); + T m = boost::math::round((phi - rphi) / constants::half_pi<T>()); + int s = 1; + if(boost::math::tools::fmod_workaround(m, T(2)) > 0.5) + { + m += 1; + s = -1; + rphi = constants::half_pi<T>() - rphi; + } + T sinp = sin(rphi); + T cosp = cos(rphi); + T c = 1 / (sinp * sinp); + T cm1 = cosp * cosp / (sinp * sinp); // c - 1 + T k2 = k * k; + if(k2 > 1) + { + return policies::raise_domain_error<T>("boost::math::ellint_d<%1%>(%1%, %1%)", "The parameter k is out of range, got k = %1%", k, pol); + } + else if(rphi == 0) + { + result = 0; + } + else + { + // http://dlmf.nist.gov/19.25#E10 + result = s * ellint_rd_imp(cm1, T(c - k2), c, pol) / 3; + } + if(m != 0) + result += m * ellint_d_imp(k, pol); + } + return invert ? T(-result) : result; +} + +// Complete elliptic integral (Legendre form) of the second kind +template <typename T, typename Policy> +T ellint_d_imp(T k, const Policy& pol) +{ + BOOST_MATH_STD_USING + using namespace boost::math::tools; + + if (abs(k) > 1) + { + return policies::raise_domain_error<T>("boost::math::ellint_e<%1%>(%1%)", + "Got k = %1%, function requires |k| <= 1", k, pol); + } + if (abs(k) == 1) + { + return static_cast<T>(1); + } + if(fabs(k) <= tools::root_epsilon<T>()) + return constants::pi<T>() / 4; + + T x = 0; + T t = k * k; + T y = 1 - t; + T z = 1; + T value = ellint_rd_imp(x, y, z, pol) / 3; + + return value; +} + +template <typename T, typename Policy> +inline typename tools::promote_args<T>::type ellint_d(T k, const Policy& pol, const mpl::true_&) +{ + typedef typename tools::promote_args<T>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_d_imp(static_cast<value_type>(k), pol), "boost::math::ellint_d<%1%>(%1%)"); +} + +// Elliptic integral (Legendre form) of the second kind +template <class T1, class T2> +inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const mpl::false_&) +{ + return boost::math::ellint_d(k, phi, policies::policy<>()); +} + +} // detail + +// Complete elliptic integral (Legendre form) of the second kind +template <typename T> +inline typename tools::promote_args<T>::type ellint_d(T k) +{ + return ellint_d(k, policies::policy<>()); +} + +// Elliptic integral (Legendre form) of the second kind +template <class T1, class T2> +inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi) +{ + typedef typename policies::is_policy<T2>::type tag_type; + return detail::ellint_d(k, phi, tag_type()); +} + +template <class T1, class T2, class Policy> +inline typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol) +{ + typedef typename tools::promote_args<T1, T2>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<result_type, Policy>(detail::ellint_d_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::ellint_2<%1%>(%1%,%1%)"); +} + +}} // namespaces + +#endif // BOOST_MATH_ELLINT_D_HPP + diff --git a/boost/math/special_functions/ellint_rc.hpp b/boost/math/special_functions/ellint_rc.hpp index 5f6d5c64bc..846c752a14 100644 --- a/boost/math/special_functions/ellint_rc.hpp +++ b/boost/math/special_functions/ellint_rc.hpp @@ -1,4 +1,4 @@ -// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock // Use, modification and distribution are subject to 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) @@ -8,6 +8,7 @@ // Summer of Code 2006. JM modified it to fit into the // Boost.Math conceptual framework better, and to correctly // handle the y < 0 case. +// Updated 2015 to use Carlson's latest methods. // #ifndef BOOST_MATH_ELLINT_RC_HPP @@ -20,6 +21,9 @@ #include <boost/math/policies/error_handling.hpp> #include <boost/math/tools/config.hpp> #include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/special_functions/log1p.hpp> +#include <boost/math/constants/constants.hpp> +#include <iostream> // Carlson's degenerate elliptic integral // R_C(x, y) = R_F(x, y, y) = 0.5 * \int_{0}^{\infty} (t+x)^{-1/2} (t+y)^{-1} dt @@ -30,11 +34,7 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> T ellint_rc_imp(T x, T y, const Policy& pol) { - T value, S, u, lambda, tolerance, prefix; - unsigned long k; - BOOST_MATH_STD_USING - using namespace boost::math::tools; static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)"; @@ -49,11 +49,9 @@ T ellint_rc_imp(T x, T y, const Policy& pol) "Argument y must not be zero but got %1%", y, pol); } - // error scales as the 6th power of tolerance - tolerance = pow(4 * tools::epsilon<T>(), T(1) / 6); - // for y < 0, the integral is singular, return Cauchy principal value - if (y < 0) + T prefix, result; + if(y < 0) { prefix = sqrt(x / (x - y)); x = x - y; @@ -62,30 +60,31 @@ T ellint_rc_imp(T x, T y, const Policy& pol) else prefix = 1; - // duplication: - k = 1; - do + if(x == 0) + { + result = constants::half_pi<T>() / sqrt(y); + } + else if(x == y) + { + result = 1 / sqrt(x); + } + else if(y > x) { - u = (x + y + y) / 3; - S = y / u - 1; // 1 - x / u = 2 * S - - if (2 * abs(S) < tolerance) - break; - - T sx = sqrt(x); - T sy = sqrt(y); - lambda = 2 * sx * sy + y; - x = (x + lambda) / 4; - y = (y + lambda) / 4; - ++k; - }while(k < policies::get_max_series_iterations<Policy>()); - // Check to see if we gave up too soon: - policies::check_series_iterations<T>(function, k, pol); - - // Taylor series expansion to the 5th order - value = (1 + S * S * (T(3) / 10 + S * (T(1) / 7 + S * (T(3) / 8 + S * T(9) / 22)))) / sqrt(u); - - return value * prefix; + result = atan(sqrt((y - x) / x)) / sqrt(y - x); + } + else + { + if(y / x > 0.5) + { + T arg = sqrt((x - y) / x); + result = (boost::math::log1p(arg) - boost::math::log1p(-arg)) / (2 * sqrt(x - y)); + } + else + { + result = log((sqrt(x) + sqrt(x - y)) / sqrt(y)) / sqrt(x - y); + } + } + return prefix * result; } } // namespace detail diff --git a/boost/math/special_functions/ellint_rd.hpp b/boost/math/special_functions/ellint_rd.hpp index 61014d3866..03b73b159f 100644 --- a/boost/math/special_functions/ellint_rd.hpp +++ b/boost/math/special_functions/ellint_rd.hpp @@ -1,4 +1,4 @@ -// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock. // Use, modification and distribution are subject to 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) @@ -7,6 +7,7 @@ // XZ wrote the original of this file as part of the Google // Summer of Code 2006. JM modified it slightly to fit into the // Boost.Math conceptual framework better. +// Updated 2015 to use Carlson's latest methods. #ifndef BOOST_MATH_ELLINT_RD_HPP #define BOOST_MATH_ELLINT_RD_HPP @@ -16,6 +17,8 @@ #endif #include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/special_functions/ellint_rc.hpp> +#include <boost/math/special_functions/pow.hpp> #include <boost/math/tools/config.hpp> #include <boost/math/policies/error_handling.hpp> @@ -28,78 +31,146 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> T ellint_rd_imp(T x, T y, T z, const Policy& pol) { - T value, u, lambda, sigma, factor, tolerance; - T X, Y, Z, EA, EB, EC, ED, EE, S1, S2; - unsigned long k; - - BOOST_MATH_STD_USING - using namespace boost::math::tools; - - static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)"; - - if (x < 0) - { - return policies::raise_domain_error<T>(function, - "Argument x must be >= 0, but got %1%", x, pol); - } - if (y < 0) - { - return policies::raise_domain_error<T>(function, - "Argument y must be >= 0, but got %1%", y, pol); - } - if (z <= 0) - { - return policies::raise_domain_error<T>(function, - "Argument z must be > 0, but got %1%", z, pol); - } - if (x + y == 0) - { - return policies::raise_domain_error<T>(function, - "At most one argument can be zero, but got, x + y = %1%", x+y, pol); - } - - // error scales as the 6th power of tolerance - tolerance = pow(tools::epsilon<T>() / 3, T(1)/6); - - // duplication - sigma = 0; - factor = 1; - k = 1; - do - { - u = (x + y + z + z + z) / 5; - X = (u - x) / u; - Y = (u - y) / u; - Z = (u - z) / u; - if ((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance) - break; - T sx = sqrt(x); - T sy = sqrt(y); - T sz = sqrt(z); - lambda = sy * (sx + sz) + sz * sx; //sqrt(x * y) + sqrt(y * z) + sqrt(z * x); - sigma += factor / (sz * (z + lambda)); - factor /= 4; - x = (x + lambda) / 4; - y = (y + lambda) / 4; - z = (z + lambda) / 4; - ++k; - } - while(k < policies::get_max_series_iterations<Policy>()); - - // Check to see if we gave up too soon: - policies::check_series_iterations<T>(function, k, pol); - - // Taylor series expansion to the 5th order - EA = X * Y; - EB = Z * Z; - EC = EA - EB; - ED = EA - 6 * EB; - EE = ED + EC + EC; - S1 = ED * (ED * T(9) / 88 - Z * EE * T(9) / 52 - T(3) / 14); - S2 = Z * (EE / 6 + Z * (-EC * T(9) / 22 + Z * EA * T(3) / 26)); - value = 3 * sigma + factor * (1 + S1 + S2) / (u * sqrt(u)); - - return value; + BOOST_MATH_STD_USING + using std::swap; + + static const char* function = "boost::math::ellint_rd<%1%>(%1%,%1%,%1%)"; + + if(x < 0) + { + return policies::raise_domain_error<T>(function, + "Argument x must be >= 0, but got %1%", x, pol); + } + if(y < 0) + { + return policies::raise_domain_error<T>(function, + "Argument y must be >= 0, but got %1%", y, pol); + } + if(z <= 0) + { + return policies::raise_domain_error<T>(function, + "Argument z must be > 0, but got %1%", z, pol); + } + if(x + y == 0) + { + return policies::raise_domain_error<T>(function, + "At most one argument can be zero, but got, x + y = %1%", x + y, pol); + } + // + // Special cases from http://dlmf.nist.gov/19.20#iv + // + using std::swap; + if(x == z) + swap(x, y); + if(y == z) + { + if(x == y) + { + return 1 / (x * sqrt(x)); + } + else if(x == 0) + { + return 3 * constants::pi<T>() / (4 * y * sqrt(y)); + } + else + { + if((std::min)(x, y) / (std::max)(x, y) > 1.3) + return 3 * (ellint_rc_imp(x, y, pol) - sqrt(x) / y) / (2 * (y - x)); + // Otherwise fall through to avoid cancellation in the above (RC(x,y) -> 1/x^0.5 as x -> y) + } + } + if(x == y) + { + if((std::min)(x, z) / (std::max)(x, z) > 1.3) + return 3 * (ellint_rc_imp(z, x, pol) - 1 / sqrt(z)) / (z - x); + // Otherwise fall through to avoid cancellation in the above (RC(x,y) -> 1/x^0.5 as x -> y) + } + if(y == 0) + swap(x, y); + if(x == 0) + { + // + // Special handling for common case, from + // Numerical Computation of Real or Complex Elliptic Integrals, eq.47 + // + T xn = sqrt(y); + T yn = sqrt(z); + T x0 = xn; + T y0 = yn; + T sum = 0; + T sum_pow = 0.25f; + + while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) + { + T t = sqrt(xn * yn); + xn = (xn + yn) / 2; + yn = t; + sum_pow *= 2; + sum += sum_pow * boost::math::pow<2>(xn - yn); + } + T RF = constants::pi<T>() / (xn + yn); + // + // This following calculation suffers from serious cancellation when y ~ z + // unless we combine terms. We have: + // + // ( ((x0 + y0)/2)^2 - z ) / (z(y-z)) + // + // Substituting y = x0^2 and z = y0^2 and simplifying we get the following: + // + T pt = (x0 + 3 * y0) / (4 * z * (x0 + y0)); + // + // Since we've moved the demoninator from eq.47 inside the expression, we + // need to also scale "sum" by the same value: + // + pt -= sum / (z * (y - z)); + return pt * RF * 3; + } + + T xn = x; + T yn = y; + T zn = z; + T An = (x + y + 3 * z) / 5; + T A0 = An; + // This has an extra 1.2 fudge factor which is really only needed when x, y and z are close in magnitude: + T Q = pow(tools::epsilon<T>() / 4, -T(1) / 8) * (std::max)((std::max)(An - x, An - y), An - z) * 1.2f; + T lambda, rx, ry, rz; + unsigned k = 0; + T fn = 1; + T RD_sum = 0; + + for(; k < policies::get_max_series_iterations<Policy>(); ++k) + { + rx = sqrt(xn); + ry = sqrt(yn); + rz = sqrt(zn); + lambda = rx * ry + rx * rz + ry * rz; + RD_sum += fn / (rz * (zn + lambda)); + An = (An + lambda) / 4; + xn = (xn + lambda) / 4; + yn = (yn + lambda) / 4; + zn = (zn + lambda) / 4; + fn /= 4; + Q /= 4; + if(Q < An) + break; + } + + policies::check_series_iterations<T, Policy>(function, k, pol); + + T X = fn * (A0 - x) / An; + T Y = fn * (A0 - y) / An; + T Z = -(X + Y) / 3; + T E2 = X * Y - 6 * Z * Z; + T E3 = (3 * X * Y - 8 * Z * Z) * Z; + T E4 = 3 * (X * Y - Z * Z) * Z * Z; + T E5 = X * Y * Z * Z * Z; + + T result = fn * pow(An, T(-3) / 2) * + (1 - 3 * E2 / 14 + E3 / 6 + 9 * E2 * E2 / 88 - 3 * E4 / 22 - 9 * E2 * E3 / 52 + 3 * E5 / 26 - E2 * E2 * E2 / 16 + + 3 * E3 * E3 / 40 + 3 * E2 * E4 / 20 + 45 * E2 * E2 * E3 / 272 - 9 * (E3 * E4 + E2 * E5) / 68); + result += 3 * RD_sum; + + return result; } } // namespace detail diff --git a/boost/math/special_functions/ellint_rf.hpp b/boost/math/special_functions/ellint_rf.hpp index ac5725783a..a8a7b4b217 100644 --- a/boost/math/special_functions/ellint_rf.hpp +++ b/boost/math/special_functions/ellint_rf.hpp @@ -1,4 +1,4 @@ -// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock // Use, modification and distribution are subject to 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) @@ -8,6 +8,7 @@ // Summer of Code 2006. JM modified it to fit into the // Boost.Math conceptual framework better, and to handle // types longer than 80-bit reals. +// Updated 2015 to use Carlson's latest methods. // #ifndef BOOST_MATH_ELLINT_RF_HPP #define BOOST_MATH_ELLINT_RF_HPP @@ -18,8 +19,9 @@ #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/tools/config.hpp> - +#include <boost/math/constants/constants.hpp> #include <boost/math/policies/error_handling.hpp> +#include <boost/math/special_functions/ellint_rc.hpp> // Carlson's elliptic integral of the first kind // R_F(x, y, z) = 0.5 * \int_{0}^{\infty} [(t+x)(t+y)(t+z)]^{-1/2} dt @@ -27,82 +29,122 @@ namespace boost { namespace math { namespace detail{ -template <typename T, typename Policy> -T ellint_rf_imp(T x, T y, T z, const Policy& pol) -{ - T value, X, Y, Z, E2, E3, u, lambda, tolerance; - unsigned long k; - - BOOST_MATH_STD_USING - using namespace boost::math::tools; + template <typename T, typename Policy> + T ellint_rf_imp(T x, T y, T z, const Policy& pol) + { + BOOST_MATH_STD_USING + using namespace boost::math; + using std::swap; - static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; + static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; - if (x < 0 || y < 0 || z < 0) - { - return policies::raise_domain_error<T>(function, + if(x < 0 || y < 0 || z < 0) + { + return policies::raise_domain_error<T>(function, "domain error, all arguments must be non-negative, " "only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol); - } - if (x + y == 0 || y + z == 0 || z + x == 0) - { - return policies::raise_domain_error<T>(function, + } + if(x + y == 0 || y + z == 0 || z + x == 0) + { + return policies::raise_domain_error<T>(function, "domain error, at most one argument can be zero, " "only sensible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol); - } - - // Carlson scales error as the 6th power of tolerance, - // but this seems not to work for types larger than - // 80-bit reals, this heuristic seems to work OK: - if(policies::digits<T, Policy>() > 64) - { - tolerance = pow(tools::epsilon<T>(), T(1)/4.25f); - BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); - } - else - { - tolerance = pow(4*tools::epsilon<T>(), T(1)/6); - BOOST_MATH_INSTRUMENT_VARIABLE(tolerance); - } - - // duplication - k = 1; - do - { - u = (x + y + z) / 3; - X = (u - x) / u; - Y = (u - y) / u; - Z = (u - z) / u; - - // Termination condition: - if ((tools::max)(abs(X), abs(Y), abs(Z)) < tolerance) - break; - - T sx = sqrt(x); - T sy = sqrt(y); - T sz = sqrt(z); - lambda = sy * (sx + sz) + sz * sx; - x = (x + lambda) / 4; - y = (y + lambda) / 4; - z = (z + lambda) / 4; - ++k; - } - while(k < policies::get_max_series_iterations<Policy>()); - - // Check to see if we gave up too soon: - policies::check_series_iterations<T>(function, k, pol); - BOOST_MATH_INSTRUMENT_VARIABLE(k); - - // Taylor series expansion to the 5th order - E2 = X * Y - Z * Z; - E3 = X * Y * Z; - value = (1 + E2*(E2/24 - E3*T(3)/44 - T(0.1)) + E3/14) / sqrt(u); - BOOST_MATH_INSTRUMENT_VARIABLE(value); - - return value; -} + } + // + // Special cases from http://dlmf.nist.gov/19.20#i + // + if(x == y) + { + if(x == z) + { + // x, y, z equal: + return 1 / sqrt(x); + } + else + { + // 2 equal, x and y: + if(z == 0) + return constants::pi<T>() / (2 * sqrt(x)); + else + return ellint_rc_imp(z, x, pol); + } + } + if(x == z) + { + if(y == 0) + return constants::pi<T>() / (2 * sqrt(x)); + else + return ellint_rc_imp(y, x, pol); + } + if(y == z) + { + if(x == 0) + return constants::pi<T>() / (2 * sqrt(y)); + else + return ellint_rc_imp(x, y, pol); + } + if(x == 0) + swap(x, z); + else if(y == 0) + swap(y, z); + if(z == 0) + { + // + // Special case for one value zero: + // + T xn = sqrt(x); + T yn = sqrt(y); + + while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) + { + T t = sqrt(xn * yn); + xn = (xn + yn) / 2; + yn = t; + } + return constants::pi<T>() / (xn + yn); + } + + T xn = x; + T yn = y; + T zn = z; + T An = (x + y + z) / 3; + T A0 = An; + T Q = pow(3 * boost::math::tools::epsilon<T>(), T(-1) / 8) * (std::max)((std::max)(fabs(An - xn), fabs(An - yn)), fabs(An - zn)); + T fn = 1; + + + // duplication + unsigned k = 1; + for(; k < boost::math::policies::get_max_series_iterations<Policy>(); ++k) + { + T root_x = sqrt(xn); + T root_y = sqrt(yn); + T root_z = sqrt(zn); + T lambda = root_x * root_y + root_x * root_z + root_y * root_z; + An = (An + lambda) / 4; + xn = (xn + lambda) / 4; + yn = (yn + lambda) / 4; + zn = (zn + lambda) / 4; + Q /= 4; + fn *= 4; + if(Q < fabs(An)) + break; + } + // Check to see if we gave up too soon: + policies::check_series_iterations<T>(function, k, pol); + BOOST_MATH_INSTRUMENT_VARIABLE(k); + + T X = (A0 - x) / (An * fn); + T Y = (A0 - y) / (An * fn); + T Z = -X - Y; + + // Taylor series expansion to the 7th order + T E2 = X * Y - Z * Z; + T E3 = X * Y * Z; + return (1 + E3 * (T(1) / 14 + 3 * E3 / 104) + E2 * (T(-1) / 10 + E2 / 24 - (3 * E3) / 44 - 5 * E2 * E2 / 208 + E2 * E3 / 16)) / sqrt(An); + } } // namespace detail diff --git a/boost/math/special_functions/ellint_rg.hpp b/boost/math/special_functions/ellint_rg.hpp new file mode 100644 index 0000000000..bb5b7c376a --- /dev/null +++ b/boost/math/special_functions/ellint_rg.hpp @@ -0,0 +1,136 @@ +// Copyright (c) 2015 John Maddock +// Use, modification and distribution are subject to 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) +// +#ifndef BOOST_MATH_ELLINT_RG_HPP +#define BOOST_MATH_ELLINT_RG_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/tools/config.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/special_functions/ellint_rd.hpp> +#include <boost/math/special_functions/ellint_rf.hpp> +#include <boost/math/special_functions/pow.hpp> + +namespace boost { namespace math { namespace detail{ + + template <typename T, typename Policy> + T ellint_rg_imp(T x, T y, T z, const Policy& pol) + { + BOOST_MATH_STD_USING + static const char* function = "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"; + + if(x < 0 || y < 0 || z < 0) + { + return policies::raise_domain_error<T>(function, + "domain error, all arguments must be non-negative, " + "only sensible result is %1%.", + std::numeric_limits<T>::quiet_NaN(), pol); + } + // + // Function is symmetric in x, y and z, but we require + // (x - z)(y - z) >= 0 to avoid cancellation error in the result + // which implies (for example) x >= z >= y + // + using std::swap; + if(x < y) + swap(x, y); + if(x < z) + swap(x, z); + if(y > z) + swap(y, z); + + BOOST_ASSERT(x >= z); + BOOST_ASSERT(z >= y); + // + // Special cases from http://dlmf.nist.gov/19.20#ii + // + if(x == z) + { + if(y == z) + { + // x = y = z + // This also works for x = y = z = 0 presumably. + return sqrt(x); + } + else if(y == 0) + { + // x = y, z = 0 + return constants::pi<T>() * sqrt(x) / 4; + } + else + { + // x = z, y != 0 + swap(x, y); + return (x == 0) ? T(sqrt(z) / 2) : T((z * ellint_rc_imp(x, z, pol) + sqrt(x)) / 2); + } + } + else if(y == z) + { + if(x == 0) + return constants::pi<T>() * sqrt(y) / 4; + else + return (y == 0) ? T(sqrt(x) / 2) : T((y * ellint_rc_imp(x, y, pol) + sqrt(x)) / 2); + } + else if(y == 0) + { + swap(y, z); + // + // Special handling for common case, from + // Numerical Computation of Real or Complex Elliptic Integrals, eq.46 + // + T xn = sqrt(x); + T yn = sqrt(y); + T x0 = xn; + T y0 = yn; + T sum = 0; + T sum_pow = 0.25f; + + while(fabs(xn - yn) >= 2.7 * tools::root_epsilon<T>() * fabs(xn)) + { + T t = sqrt(xn * yn); + xn = (xn + yn) / 2; + yn = t; + sum_pow *= 2; + sum += sum_pow * boost::math::pow<2>(xn - yn); + } + T RF = constants::pi<T>() / (xn + yn); + return ((boost::math::pow<2>((x0 + y0) / 2) - sum) * RF) / 2; + } + return (z * ellint_rf_imp(x, y, z, pol) + - (x - z) * (y - z) * ellint_rd_imp(x, y, z, pol) / 3 + + sqrt(x * y / z)) / 2; + } + +} // namespace detail + +template <class T1, class T2, class T3, class Policy> +inline typename tools::promote_args<T1, T2, T3>::type + ellint_rg(T1 x, T2 y, T3 z, const Policy& pol) +{ + typedef typename tools::promote_args<T1, T2, T3>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<result_type, Policy>( + detail::ellint_rg_imp( + static_cast<value_type>(x), + static_cast<value_type>(y), + static_cast<value_type>(z), pol), "boost::math::ellint_rf<%1%>(%1%,%1%,%1%)"); +} + +template <class T1, class T2, class T3> +inline typename tools::promote_args<T1, T2, T3>::type + ellint_rg(T1 x, T2 y, T3 z) +{ + return ellint_rg(x, y, z, policies::policy<>()); +} + +}} // namespaces + +#endif // BOOST_MATH_ELLINT_RG_HPP + diff --git a/boost/math/special_functions/ellint_rj.hpp b/boost/math/special_functions/ellint_rj.hpp index 8a242f06a4..ac39bed678 100644 --- a/boost/math/special_functions/ellint_rj.hpp +++ b/boost/math/special_functions/ellint_rj.hpp @@ -1,4 +1,4 @@ -// Copyright (c) 2006 Xiaogang Zhang +// Copyright (c) 2006 Xiaogang Zhang, 2015 John Maddock // Use, modification and distribution are subject to 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) @@ -8,6 +8,7 @@ // Summer of Code 2006. JM modified it to fit into the // Boost.Math conceptual framework better, and to correctly // handle the p < 0 case. +// Updated 2015 to use Carlson's latest methods. // #ifndef BOOST_MATH_ELLINT_RJ_HPP @@ -22,6 +23,7 @@ #include <boost/math/policies/error_handling.hpp> #include <boost/math/special_functions/ellint_rc.hpp> #include <boost/math/special_functions/ellint_rf.hpp> +#include <boost/math/special_functions/ellint_rd.hpp> // Carlson's elliptic integral of the third kind // R_J(x, y, z, p) = 1.5 * \int_{0}^{\infty} (t+p)^{-1} [(t+x)(t+y)(t+z)]^{-1/2} dt @@ -30,124 +32,244 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> +T ellint_rc1p_imp(T y, const Policy& pol) +{ + using namespace boost::math; + // Calculate RC(1, 1 + x) + BOOST_MATH_STD_USING + + static const char* function = "boost::math::ellint_rc<%1%>(%1%,%1%)"; + + if(y == -1) + { + return policies::raise_domain_error<T>(function, + "Argument y must not be zero but got %1%", y, pol); + } + + // for 1 + y < 0, the integral is singular, return Cauchy principal value + T result; + if(y < -1) + { + result = sqrt(1 / -y) * detail::ellint_rc_imp(T(-y), T(-1 - y), pol); + } + else if(y == 0) + { + result = 1; + } + else if(y > 0) + { + result = atan(sqrt(y)) / sqrt(y); + } + else + { + if(y > -0.5) + { + T arg = sqrt(-y); + result = (boost::math::log1p(arg) - boost::math::log1p(-arg)) / (2 * sqrt(-y)); + } + else + { + result = log((1 + sqrt(-y)) / sqrt(1 + y)) / sqrt(-y); + } + } + return result; +} + +template <typename T, typename Policy> T ellint_rj_imp(T x, T y, T z, T p, const Policy& pol) { - T value, u, lambda, alpha, beta, sigma, factor, tolerance; - T X, Y, Z, P, EA, EB, EC, E2, E3, S1, S2, S3; - unsigned long k; - - BOOST_MATH_STD_USING - using namespace boost::math::tools; - - static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)"; - - if (x < 0) - { - return policies::raise_domain_error<T>(function, - "Argument x must be non-negative, but got x = %1%", x, pol); - } - if(y < 0) - { - return policies::raise_domain_error<T>(function, - "Argument y must be non-negative, but got y = %1%", y, pol); - } - if(z < 0) - { - return policies::raise_domain_error<T>(function, - "Argument z must be non-negative, but got z = %1%", z, pol); - } - if(p == 0) - { - return policies::raise_domain_error<T>(function, - "Argument p must not be zero, but got p = %1%", p, pol); - } - if (x + y == 0 || y + z == 0 || z + x == 0) - { - return policies::raise_domain_error<T>(function, - "At most one argument can be zero, " - "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol); - } - - // error scales as the 6th power of tolerance - tolerance = pow(T(1) * tools::epsilon<T>() / 3, T(1) / 6); - - // for p < 0, the integral is singular, return Cauchy principal value - if (p < 0) - { - // - // We must ensure that (z - y) * (y - x) is positive. - // Since the integral is symmetrical in x, y and z - // we can just permute the values: - // - if(x > y) - std::swap(x, y); - if(y > z) - std::swap(y, z); - if(x > y) - std::swap(x, y); - - T q = -p; - T pmy = (z - y) * (y - x) / (y + q); // p - y - - BOOST_ASSERT(pmy >= 0); - - p = pmy + y; - value = boost::math::ellint_rj(x, y, z, p, pol); - value *= pmy; - value -= 3 * boost::math::ellint_rf(x, y, z, pol); - value += 3 * sqrt((x * y * z) / (x * z + p * q)) * boost::math::ellint_rc(x * z + p * q, p * q, pol); - value /= (y + q); - return value; - } - - // duplication - sigma = 0; - factor = 1; - k = 1; - do - { - u = (x + y + z + p + p) / 5; - X = (u - x) / u; - Y = (u - y) / u; - Z = (u - z) / u; - P = (u - p) / u; - - if ((tools::max)(abs(X), abs(Y), abs(Z), abs(P)) < tolerance) - break; - - T sx = sqrt(x); - T sy = sqrt(y); - T sz = sqrt(z); - - lambda = sy * (sx + sz) + sz * sx; - alpha = p * (sx + sy + sz) + sx * sy * sz; - alpha *= alpha; - beta = p * (p + lambda) * (p + lambda); - sigma += factor * boost::math::ellint_rc(alpha, beta, pol); - factor /= 4; - x = (x + lambda) / 4; - y = (y + lambda) / 4; - z = (z + lambda) / 4; - p = (p + lambda) / 4; - ++k; - } - while(k < policies::get_max_series_iterations<Policy>()); - - // Check to see if we gave up too soon: - policies::check_series_iterations<T>(function, k, pol); - - // Taylor series expansion to the 5th order - EA = X * Y + Y * Z + Z * X; - EB = X * Y * Z; - EC = P * P; - E2 = EA - 3 * EC; - E3 = EB + 2 * P * (EA - EC); - S1 = 1 + E2 * (E2 * T(9) / 88 - E3 * T(9) / 52 - T(3) / 14); - S2 = EB * (T(1) / 6 + P * (T(-6) / 22 + P * T(3) / 26)); - S3 = P * ((EA - EC) / 3 - P * EA * T(3) / 22); - value = 3 * sigma + factor * (S1 + S2 + S3) / (u * sqrt(u)); - - return value; + BOOST_MATH_STD_USING + + static const char* function = "boost::math::ellint_rj<%1%>(%1%,%1%,%1%)"; + + if(x < 0) + { + return policies::raise_domain_error<T>(function, + "Argument x must be non-negative, but got x = %1%", x, pol); + } + if(y < 0) + { + return policies::raise_domain_error<T>(function, + "Argument y must be non-negative, but got y = %1%", y, pol); + } + if(z < 0) + { + return policies::raise_domain_error<T>(function, + "Argument z must be non-negative, but got z = %1%", z, pol); + } + if(p == 0) + { + return policies::raise_domain_error<T>(function, + "Argument p must not be zero, but got p = %1%", p, pol); + } + if(x + y == 0 || y + z == 0 || z + x == 0) + { + return policies::raise_domain_error<T>(function, + "At most one argument can be zero, " + "only possible result is %1%.", std::numeric_limits<T>::quiet_NaN(), pol); + } + + // for p < 0, the integral is singular, return Cauchy principal value + if(p < 0) + { + // + // We must ensure that x < y < z. + // Since the integral is symmetrical in x, y and z + // we can just permute the values: + // + if(x > y) + std::swap(x, y); + if(y > z) + std::swap(y, z); + if(x > y) + std::swap(x, y); + + BOOST_ASSERT(x <= y); + BOOST_ASSERT(y <= z); + + T q = -p; + p = (z * (x + y + q) - x * y) / (z + q); + + BOOST_ASSERT(p >= 0); + + T value = (p - z) * ellint_rj_imp(x, y, z, p, pol); + value -= 3 * ellint_rf_imp(x, y, z, pol); + value += 3 * sqrt((x * y * z) / (x * y + p * q)) * ellint_rc_imp(T(x * y + p * q), T(p * q), pol); + value /= (z + q); + return value; + } + + // + // Special cases from http://dlmf.nist.gov/19.20#iii + // + if(x == y) + { + if(x == z) + { + if(x == p) + { + // All values equal: + return 1 / (x * sqrt(x)); + } + else + { + // x = y = z: + return 3 * (ellint_rc_imp(x, p, pol) - 1 / sqrt(x)) / (x - p); + } + } + else + { + // x = y only, permute so y = z: + using std::swap; + swap(x, z); + if(y == p) + { + return ellint_rd_imp(x, y, y, pol); + } + else if((std::max)(y, p) / (std::min)(y, p) > 1.2) + { + return 3 * (ellint_rc_imp(x, y, pol) - ellint_rc_imp(x, p, pol)) / (p - y); + } + // Otherwise fall through to normal method, special case above will suffer too much cancellation... + } + } + if(y == z) + { + if(y == p) + { + // y = z = p: + return ellint_rd_imp(x, y, y, pol); + } + else if((std::max)(y, p) / (std::min)(y, p) > 1.2) + { + // y = z: + return 3 * (ellint_rc_imp(x, y, pol) - ellint_rc_imp(x, p, pol)) / (p - y); + } + // Otherwise fall through to normal method, special case above will suffer too much cancellation... + } + if(z == p) + { + return ellint_rd_imp(x, y, z, pol); + } + + T xn = x; + T yn = y; + T zn = z; + T pn = p; + T An = (x + y + z + 2 * p) / 5; + T A0 = An; + T delta = (p - x) * (p - y) * (p - z); + T Q = pow(tools::epsilon<T>() / 5, -T(1) / 8) * (std::max)((std::max)(fabs(An - x), fabs(An - y)), (std::max)(fabs(An - z), fabs(An - p))); + + unsigned n; + T lambda; + T Dn; + T En; + T rx, ry, rz, rp; + T fmn = 1; // 4^-n + T RC_sum = 0; + + for(n = 0; n < policies::get_max_series_iterations<Policy>(); ++n) + { + rx = sqrt(xn); + ry = sqrt(yn); + rz = sqrt(zn); + rp = sqrt(pn); + Dn = (rp + rx) * (rp + ry) * (rp + rz); + En = delta / Dn; + En /= Dn; + if((En < -0.5) && (En > -1.5)) + { + // + // Occationally En ~ -1, we then have no means of calculating + // RC(1, 1+En) without terrible cancellation error, so we + // need to get to 1+En directly. By substitution we have + // + // 1+E_0 = 1 + (p-x)*(p-y)*(p-z)/((sqrt(p) + sqrt(x))*(sqrt(p)+sqrt(y))*(sqrt(p)+sqrt(z)))^2 + // = 2*sqrt(p)*(p+sqrt(x) * (sqrt(y)+sqrt(z)) + sqrt(y)*sqrt(z)) / ((sqrt(p) + sqrt(x))*(sqrt(p) + sqrt(y)*(sqrt(p)+sqrt(z)))) + // + // And since this is just an application of the duplication formula for RJ, the same + // expression works for 1+En if we use x,y,z,p_n etc. + // This branch is taken only once or twice at the start of iteration, + // after than En reverts to it's usual very small values. + // + T b = 2 * rp * (pn + rx * (ry + rz) + ry * rz) / Dn; + RC_sum += fmn / Dn * detail::ellint_rc_imp(T(1), b, pol); + } + else + { + RC_sum += fmn / Dn * ellint_rc1p_imp(En, pol); + } + lambda = rx * ry + rx * rz + ry * rz; + + // From here on we move to n+1: + An = (An + lambda) / 4; + fmn /= 4; + + if(fmn * Q < An) + break; + + xn = (xn + lambda) / 4; + yn = (yn + lambda) / 4; + zn = (zn + lambda) / 4; + pn = (pn + lambda) / 4; + delta /= 64; + } + + T X = fmn * (A0 - x) / An; + T Y = fmn * (A0 - y) / An; + T Z = fmn * (A0 - z) / An; + T P = (-X - Y - Z) / 2; + T E2 = X * Y + X * Z + Y * Z - 3 * P * P; + T E3 = X * Y * Z + 2 * E2 * P + 4 * P * P * P; + T E4 = (2 * X * Y * Z + E2 * P + 3 * P * P * P) * P; + T E5 = X * Y * Z * P * P; + T result = fmn * pow(An, T(-3) / 2) * + (1 - 3 * E2 / 14 + E3 / 6 + 9 * E2 * E2 / 88 - 3 * E4 / 22 - 9 * E2 * E3 / 52 + 3 * E5 / 26 - E2 * E2 * E2 / 16 + + 3 * E3 * E3 / 40 + 3 * E2 * E4 / 20 + 45 * E2 * E2 * E3 / 272 - 9 * (E3 * E4 + E2 * E5) / 68); + + result += 6 * RC_sum; + return result; } } // namespace detail diff --git a/boost/math/special_functions/factorials.hpp b/boost/math/special_functions/factorials.hpp index de24642ac4..e36a098bb6 100644 --- a/boost/math/special_functions/factorials.hpp +++ b/boost/math/special_functions/factorials.hpp @@ -151,7 +151,7 @@ T rising_factorial_imp(T x, int n, const Policy& pol) if((x < 1) && (x + n < 0)) { T val = boost::math::tgamma_delta_ratio(1 - x, static_cast<T>(-n), pol); - return (n & 1) ? -val : val; + return (n & 1) ? T(-val) : val; } // // We don't optimise this for small n, because diff --git a/boost/math/special_functions/fpclassify.hpp b/boost/math/special_functions/fpclassify.hpp index 40f6e14ba5..8e75fae0f2 100644 --- a/boost/math/special_functions/fpclassify.hpp +++ b/boost/math/special_functions/fpclassify.hpp @@ -309,7 +309,7 @@ namespace detail { if(std::numeric_limits<T>::is_specialized) return isfinite_impl(x, generic_tag<true>()); #endif - (void)x; // warning supression. + (void)x; // warning suppression. return true; } @@ -453,7 +453,7 @@ namespace detail { if(std::numeric_limits<T>::is_specialized) return isinf_impl(x, generic_tag<true>()); #endif - (void)x; // warning supression. + (void)x; // warning suppression. return false; } @@ -541,7 +541,7 @@ namespace detail { if(std::numeric_limits<T>::is_specialized) return isnan_impl(x, generic_tag<true>()); #endif - (void)x; // warning supression + (void)x; // warning suppression return false; } diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp index b6b4c574a9..1db1e4c2d3 100644 --- a/boost/math/special_functions/gamma.hpp +++ b/boost/math/special_functions/gamma.hpp @@ -88,10 +88,6 @@ T sinpx(T z) { z = -z; } - else - { - sign = -sign; - } T fl = floor(z); T dist; if(is_odd(fl)) @@ -1063,8 +1059,31 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, { result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative); if(result == 0) - return policies::raise_evaluation_error<T>(function, "Obtained %1% for the incomplete gamma function, but in truth we don't really know what the answer is...", result, pol); - result = log(result) + boost::math::lgamma(a, pol); + { + if(invert) + { + // Try http://functions.wolfram.com/06.06.06.0039.01 + result = 1 + 1 / (12 * a) + 1 / (288 * a * a); + result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>()); + if(p_derivative) + *p_derivative = exp(a * log(x) - x); + } + else + { + // This is method 2 below, done in logs, we're really outside the + // range of this method, but since the result is almost certainly + // infinite, we should probably be OK: + result = a * log(x) - x; + if(p_derivative) + *p_derivative = exp(result); + T init_value = 0; + result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); + } + } + else + { + result = log(result) + boost::math::lgamma(a, pol); + } } if(result > tools::log_max_value<T>()) return policies::raise_overflow_error<T>(function, 0, pol); @@ -1560,6 +1579,7 @@ T tgamma_ratio_imp(T x, T y, const Policy& pol) template <class T, class Policy> T gamma_p_derivative_imp(T a, T x, const Policy& pol) { + BOOST_MATH_STD_USING // // Usual error checks first: // @@ -1585,8 +1605,14 @@ T gamma_p_derivative_imp(T a, T x, const Policy& pol) // overflow: return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol); } - - f1 /= x; + if(f1 == 0) + { + // Underflow in calculation, use logs instead: + f1 = a * log(x) - x - lgamma(a, pol) - log(x); + f1 = exp(f1); + } + else + f1 /= x; return f1; } diff --git a/boost/math/special_functions/heuman_lambda.hpp b/boost/math/special_functions/heuman_lambda.hpp new file mode 100644 index 0000000000..6389443ee4 --- /dev/null +++ b/boost/math/special_functions/heuman_lambda.hpp @@ -0,0 +1,87 @@ +// Copyright (c) 2015 John Maddock +// Use, modification and distribution are subject to 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) + +#ifndef BOOST_MATH_ELLINT_HL_HPP +#define BOOST_MATH_ELLINT_HL_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/special_functions/ellint_rj.hpp> +#include <boost/math/special_functions/ellint_rj.hpp> +#include <boost/math/special_functions/ellint_1.hpp> +#include <boost/math/special_functions/jacobi_zeta.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/tools/workaround.hpp> + +// Elliptic integral the Jacobi Zeta function. + +namespace boost { namespace math { + +namespace detail{ + +// Elliptic integral - Jacobi Zeta +template <typename T, typename Policy> +T heuman_lambda_imp(T phi, T k, const Policy& pol) +{ + BOOST_MATH_STD_USING + using namespace boost::math::tools; + using namespace boost::math::constants; + + const char* function = "boost::math::heuman_lambda<%1%>(%1%, %1%)"; + + if(fabs(k) > 1) + return policies::raise_domain_error<T>(function, "We require |k| <= 1 but got k = %1%", k, pol); + + T result; + T sinp = sin(phi); + T cosp = cos(phi); + T s2 = sinp * sinp; + T k2 = k * k; + T kp = 1 - k2; + T delta = sqrt(1 - (kp * s2)); + if(fabs(phi) <= constants::half_pi<T>()) + { + result = kp * sinp * cosp / (delta * constants::half_pi<T>()); + result *= ellint_rf_imp(T(0), kp, T(1), pol) + k2 * ellint_rj(T(0), kp, T(1), T(1 - k2 / (delta * delta)), pol) / (3 * delta * delta); + } + else + { + T rkp = sqrt(kp); + T ratio; + if(rkp == 1) + { + return policies::raise_domain_error<T>(function, "When 1-k^2 == 1 then phi must be < Pi/2, but got phi = %1%", phi, pol); + } + else + ratio = ellint_f_imp(phi, rkp, pol) / ellint_k_imp(rkp, pol); + result = ratio + ellint_k_imp(k, pol) * jacobi_zeta_imp(phi, rkp, pol) / constants::half_pi<T>(); + } + return result; +} + +} // detail + +template <class T1, class T2, class Policy> +inline typename tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi, const Policy& pol) +{ + typedef typename tools::promote_args<T1, T2>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<result_type, Policy>(detail::heuman_lambda_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::heuman_lambda<%1%>(%1%,%1%)"); +} + +template <class T1, class T2> +inline typename tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi) +{ + return boost::math::heuman_lambda(k, phi, policies::policy<>()); +} + +}} // namespaces + +#endif // BOOST_MATH_ELLINT_D_HPP + diff --git a/boost/math/special_functions/jacobi_zeta.hpp b/boost/math/special_functions/jacobi_zeta.hpp new file mode 100644 index 0000000000..a3fa54746e --- /dev/null +++ b/boost/math/special_functions/jacobi_zeta.hpp @@ -0,0 +1,74 @@ +// Copyright (c) 2015 John Maddock +// Use, modification and distribution are subject to 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) +// + +#ifndef BOOST_MATH_ELLINT_JZ_HPP +#define BOOST_MATH_ELLINT_JZ_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/special_functions/ellint_1.hpp> +#include <boost/math/special_functions/ellint_rj.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/tools/workaround.hpp> + +// Elliptic integral the Jacobi Zeta function. + +namespace boost { namespace math { + +namespace detail{ + +// Elliptic integral - Jacobi Zeta +template <typename T, typename Policy> +T jacobi_zeta_imp(T phi, T k, const Policy& pol) +{ + BOOST_MATH_STD_USING + using namespace boost::math::tools; + using namespace boost::math::constants; + + bool invert = false; + if(phi < 0) + { + phi = fabs(phi); + invert = true; + } + + T result; + T sinp = sin(phi); + T cosp = cos(phi); + T s2 = sinp * sinp; + T k2 = k * k; + T kp = 1 - k2; + if(k == 1) + result = sinp * (boost::math::sign)(cosp); // We get here by simplifying JacobiZeta[w, 1] in Mathematica, and the fact that 0 <= phi. + else + result = k2 * sinp * cosp * sqrt(1 - k2 * s2) * ellint_rj_imp(T(0), kp, T(1), T(1 - k2 * s2), pol) / (3 * ellint_k_imp(k, pol)); + return invert ? T(-result) : result; +} + +} // detail + +template <class T1, class T2, class Policy> +inline typename tools::promote_args<T1, T2>::type jacobi_zeta(T1 k, T2 phi, const Policy& pol) +{ + typedef typename tools::promote_args<T1, T2>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<result_type, Policy>(detail::jacobi_zeta_imp(static_cast<value_type>(phi), static_cast<value_type>(k), pol), "boost::math::jacobi_zeta<%1%>(%1%,%1%)"); +} + +template <class T1, class T2> +inline typename tools::promote_args<T1, T2>::type jacobi_zeta(T1 k, T2 phi) +{ + return boost::math::jacobi_zeta(k, phi, policies::policy<>()); +} + +}} // namespaces + +#endif // BOOST_MATH_ELLINT_D_HPP + diff --git a/boost/math/special_functions/legendre.hpp b/boost/math/special_functions/legendre.hpp index 79e9756763..1a2ef5d615 100644 --- a/boost/math/special_functions/legendre.hpp +++ b/boost/math/special_functions/legendre.hpp @@ -71,7 +71,7 @@ T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false) } // namespace detail template <class T, class Policy> -inline typename tools::promote_args<T>::type +inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type legendre_p(int l, T x, const Policy& pol) { typedef typename tools::promote_args<T>::type result_type; @@ -90,7 +90,7 @@ inline typename tools::promote_args<T>::type } template <class T, class Policy> -inline typename tools::promote_args<T>::type +inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type legendre_q(unsigned l, T x, const Policy& pol) { typedef typename tools::promote_args<T>::type result_type; diff --git a/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp index e952dcdb51..96f60726a6 100644 --- a/boost/math/special_functions/math_fwd.hpp +++ b/boost/math/special_functions/math_fwd.hpp @@ -27,6 +27,7 @@ #include <boost/math/tools/promotion.hpp> // for argument promotion. #include <boost/math/policies/policy.hpp> #include <boost/mpl/comparison.hpp> +#include <boost/utility/enable_if.hpp> #include <boost/config/no_tr1/complex.hpp> #define BOOST_NO_MACRO_EXPAND /**/ @@ -145,6 +146,12 @@ namespace boost typename tools::promote_args<RT1, RT2, RT3>::type ibeta_derivative(RT1 a, RT2 b, RT3 x, const Policy& pol); // derivative of incomplete beta + // Binomial: + template <class T, class Policy> + T binomial_coefficient(unsigned n, unsigned k, const Policy& pol); + template <class T> + T binomial_coefficient(unsigned n, unsigned k); + // erf & erfc error functions. template <class RT> // Error function. typename tools::promote_args<RT>::type erf(RT z); @@ -176,7 +183,7 @@ namespace boost legendre_p(int l, T x); template <class T, class Policy> - typename tools::promote_args<T>::type + typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type legendre_p(int l, T x, const Policy& pol); template <class T> @@ -184,7 +191,7 @@ namespace boost legendre_q(unsigned l, T x); template <class T, class Policy> - typename tools::promote_args<T>::type + typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type legendre_q(unsigned l, T x, const Policy& pol); template <class T1, class T2, class T3> @@ -298,6 +305,14 @@ namespace boost typename tools::promote_args<T1, T2, T3, T4>::type ellint_rj(T1 x, T2 y, T3 z, T4 p, const Policy& pol); + template <class T1, class T2, class T3> + typename tools::promote_args<T1, T2, T3>::type + ellint_rg(T1 x, T2 y, T3 z); + + template <class T1, class T2, class T3, class Policy> + typename tools::promote_args<T1, T2, T3>::type + ellint_rg(T1 x, T2 y, T3 z, const Policy& pol); + template <typename T> typename tools::promote_args<T>::type ellint_2(T k); @@ -316,6 +331,27 @@ namespace boost template <class T1, class T2, class Policy> typename tools::promote_args<T1, T2>::type ellint_1(T1 k, T2 phi, const Policy& pol); + template <typename T> + typename tools::promote_args<T>::type ellint_d(T k); + + template <class T1, class T2> + typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi); + + template <class T1, class T2, class Policy> + typename tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi, const Policy& pol); + + template <class T1, class T2> + typename tools::promote_args<T1, T2>::type jacobi_zeta(T1 k, T2 phi); + + template <class T1, class T2, class Policy> + typename tools::promote_args<T1, T2>::type jacobi_zeta(T1 k, T2 phi, const Policy& pol); + + template <class T1, class T2> + typename tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi); + + template <class T1, class T2, class Policy> + typename tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi, const Policy& pol); + namespace detail{ template <class T, class U, class V> @@ -463,6 +499,20 @@ namespace boost template <class T, class Policy> typename tools::promote_args<T>::type digamma(T x, const Policy&); + // trigamma: + template <class T> + typename tools::promote_args<T>::type trigamma(T x); + + template <class T, class Policy> + typename tools::promote_args<T>::type trigamma(T x, const Policy&); + + // polygamma: + template <class T> + typename tools::promote_args<T>::type polygamma(int n, T x); + + template <class T, class Policy> + typename tools::promote_args<T>::type polygamma(int n, T x, const Policy&); + // Hypotenuse function sqrt(x ^ 2 + y ^ 2). template <class T1, class T2> typename tools::promote_args<T1, T2>::type @@ -1052,6 +1102,8 @@ namespace boost inline typename boost::math::tools::promote_args<RT1, RT2, RT3>::type \ ibeta_derivative(RT1 a, RT2 b, RT3 x){ return ::boost::math::ibeta_derivative(a, b, x, Policy()); }\ \ + template <class T> T binomial_coefficient(unsigned n, unsigned k){ return ::boost::math::binomial_coefficient<T, Policy>(n, k, Policy()); }\ +\ template <class RT>\ inline typename boost::math::tools::promote_args<RT>::type erf(RT z) { return ::boost::math::erf(z, Policy()); }\ \ @@ -1128,6 +1180,10 @@ namespace boost inline typename boost::math::tools::promote_args<T1, T2, T3, T4>::type \ ellint_rj(T1 x, T2 y, T3 z, T4 p){ return boost::math::ellint_rj(x, y, z, p, Policy()); }\ \ + template <class T1, class T2, class T3>\ + inline typename boost::math::tools::promote_args<T1, T2, T3>::type \ + ellint_rg(T1 x, T2 y, T3 z){ return ::boost::math::ellint_rg(x, y, z, Policy()); }\ + \ template <typename T>\ inline typename boost::math::tools::promote_args<T>::type ellint_2(T k){ return boost::math::ellint_2(k, Policy()); }\ \ @@ -1135,6 +1191,18 @@ namespace boost inline typename boost::math::tools::promote_args<T1, T2>::type ellint_2(T1 k, T2 phi){ return boost::math::ellint_2(k, phi, Policy()); }\ \ template <typename T>\ + inline typename boost::math::tools::promote_args<T>::type ellint_d(T k){ return boost::math::ellint_d(k, Policy()); }\ +\ + template <class T1, class T2>\ + inline typename boost::math::tools::promote_args<T1, T2>::type ellint_d(T1 k, T2 phi){ return boost::math::ellint_d(k, phi, Policy()); }\ +\ + template <class T1, class T2>\ + inline typename boost::math::tools::promote_args<T1, T2>::type jacobi_zeta(T1 k, T2 phi){ return boost::math::jacobi_zeta(k, phi, Policy()); }\ +\ + template <class T1, class T2>\ + inline typename boost::math::tools::promote_args<T1, T2>::type heuman_lambda(T1 k, T2 phi){ return boost::math::heuman_lambda(k, phi, Policy()); }\ +\ + template <typename T>\ inline typename boost::math::tools::promote_args<T>::type ellint_1(T k){ return boost::math::ellint_1(k, Policy()); }\ \ template <class T1, class T2>\ @@ -1205,6 +1273,12 @@ namespace boost template <class T>\ inline typename boost::math::tools::promote_args<T>::type digamma(T x){ return boost::math::digamma(x, Policy()); }\ \ + template <class T>\ + inline typename boost::math::tools::promote_args<T>::type trigamma(T x){ return boost::math::trigamma(x, Policy()); }\ +\ + template <class T>\ + inline typename boost::math::tools::promote_args<T>::type polygamma(int n, T x){ return boost::math::polygamma(n, x, Policy()); }\ + \ template <class T1, class T2>\ inline typename boost::math::tools::promote_args<T1, T2>::type \ hypot(T1 x, T2 y){ return boost::math::hypot(x, y, Policy()); }\ diff --git a/boost/math/special_functions/polygamma.hpp b/boost/math/special_functions/polygamma.hpp new file mode 100644 index 0000000000..6b7815d5e0 --- /dev/null +++ b/boost/math/special_functions/polygamma.hpp @@ -0,0 +1,83 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2013 Nikhar Agrawal +// Copyright 2013 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2013 Paul Bristow +// 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) + +#ifndef _BOOST_POLYGAMMA_2013_07_30_HPP_ + #define _BOOST_POLYGAMMA_2013_07_30_HPP_ + +#include <boost/math/special_functions/factorials.hpp> +#include <boost/math/special_functions/detail/polygamma.hpp> +#include <boost/math/special_functions/trigamma.hpp> + +namespace boost { namespace math { + + + template<class T, class Policy> + inline typename tools::promote_args<T>::type polygamma(const int n, T x, const Policy& pol) + { + // + // Filter off special cases right at the start: + // + if(n == 0) + return boost::math::digamma(x, pol); + if(n == 1) + return boost::math::trigamma(x, pol); + // + // We've found some standard library functions to misbehave if any FPU exception flags + // are set prior to their call, this code will clear those flags, then reset them + // on exit: + // + BOOST_FPU_EXCEPTION_GUARD + // + // The type of the result - the common type of T and U after + // any integer types have been promoted to double: + // + typedef typename tools::promote_args<T>::type result_type; + // + // The type used for the calculation. This may be a wider type than + // the result in order to ensure full precision: + // + typedef typename policies::evaluation<result_type, Policy>::type value_type; + // + // The type of the policy to forward to the actual implementation. + // We disable promotion of float and double as that's [possibly] + // happened already in the line above. Also reset to the default + // any policies we don't use (reduces code bloat if we're called + // multiple times with differing policies we don't actually use). + // Also normalise the type, again to reduce code bloat in case we're + // called multiple times with functionally identical policies that happen + // to be different types. + // + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + // + // Whew. Now we can make the actual call to the implementation. + // Arguments are explicitly cast to the evaluation type, and the result + // passed through checked_narrowing_cast which handles things like overflow + // according to the policy passed: + // + return policies::checked_narrowing_cast<result_type, forwarding_policy>( + detail::polygamma_imp(n, static_cast<value_type>(x), forwarding_policy()), + "boost::math::polygamma<%1%>(int, %1%)"); + } + + template<class T> + inline typename tools::promote_args<T>::type polygamma(const int n, T x) + { + return boost::math::polygamma(n, x, policies::policy<>()); + } + +} } // namespace boost::math + +#endif // _BOOST_BERNOULLI_2013_05_30_HPP_ + diff --git a/boost/math/special_functions/sign.hpp b/boost/math/special_functions/sign.hpp index f5c562d44c..3324c90a87 100644 --- a/boost/math/special_functions/sign.hpp +++ b/boost/math/special_functions/sign.hpp @@ -56,11 +56,11 @@ namespace detail { // inline int signbit_impl(long double x, generic_tag<true> const&) { - return boost::math::signbit(static_cast<double>(x)); + return (boost::math::signbit)(static_cast<double>(x)); } inline int signbit_impl(long double x, generic_tag<false> const&) { - return boost::math::signbit(static_cast<double>(x)); + return (boost::math::signbit)(static_cast<double>(x)); } #endif diff --git a/boost/math/special_functions/sin_pi.hpp b/boost/math/special_functions/sin_pi.hpp index 16aed51d2b..ae6b3e7442 100644 --- a/boost/math/special_functions/sin_pi.hpp +++ b/boost/math/special_functions/sin_pi.hpp @@ -53,10 +53,17 @@ T sin_pi_imp(T x, const Policy& pol) } // namespace detail template <class T, class Policy> -inline typename tools::promote_args<T>::type sin_pi(T x, const Policy& pol) +inline typename tools::promote_args<T>::type sin_pi(T x, const Policy&) { typedef typename tools::promote_args<T>::type result_type; - return boost::math::detail::sin_pi_imp<result_type>(x, pol); + typedef typename policies::evaluation<result_type, Policy>::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast<result_type, forwarding_policy>(boost::math::detail::sin_pi_imp<value_type>(x, forwarding_policy()), "cos_pi"); } template <class T> diff --git a/boost/math/special_functions/trigamma.hpp b/boost/math/special_functions/trigamma.hpp new file mode 100644 index 0000000000..6fccb36a3a --- /dev/null +++ b/boost/math/special_functions/trigamma.hpp @@ -0,0 +1,469 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to 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) + +#ifndef BOOST_MATH_SF_TRIGAMMA_HPP +#define BOOST_MATH_SF_TRIGAMMA_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include <boost/math/special_functions/math_fwd.hpp> +#include <boost/math/tools/rational.hpp> +#include <boost/math/tools/series.hpp> +#include <boost/math/tools/promotion.hpp> +#include <boost/math/policies/error_handling.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/mpl/comparison.hpp> +#include <boost/math/tools/big_constant.hpp> +#include <boost/math/special_functions/polygamma.hpp> + +namespace boost{ +namespace math{ +namespace detail{ + +template<class T, class Policy> +T polygamma_imp(const int n, T x, const Policy &pol); + +template <class T, class Policy> +T trigamma_prec(T x, const mpl::int_<53>*, const Policy&) +{ + // Max error in interpolated form: 3.736e-017 + static const T offset = BOOST_MATH_BIG_CONSTANT(T, 53, 2.1093254089355469); + static const T P_1_2[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, -1.1093280605946045), + BOOST_MATH_BIG_CONSTANT(T, 53, -3.8310674472619321), + BOOST_MATH_BIG_CONSTANT(T, 53, -3.3703848401898283), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.28080574467981213), + BOOST_MATH_BIG_CONSTANT(T, 53, 1.6638069578676164), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468386819102836), + }; + static const T Q_1_2[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 53, 3.4535389668541151), + BOOST_MATH_BIG_CONSTANT(T, 53, 4.5208926987851437), + BOOST_MATH_BIG_CONSTANT(T, 53, 2.7012734178351534), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.64468798399785611), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.20314516859987728e-6), + }; + // Max error in interpolated form: 1.159e-017 + static const T P_2_4[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, -0.13803835004508849e-7), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.50000049158540261), + BOOST_MATH_BIG_CONSTANT(T, 53, 1.6077979838469348), + BOOST_MATH_BIG_CONSTANT(T, 53, 2.5645435828098254), + BOOST_MATH_BIG_CONSTANT(T, 53, 2.0534873203680393), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.74566981111565923), + }; + static const T Q_2_4[] = { + BOOST_MATH_BIG_CONSTANT(T, 53, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 53, 2.8822787662376169), + BOOST_MATH_BIG_CONSTANT(T, 53, 4.1681660554090917), + BOOST_MATH_BIG_CONSTANT(T, 53, 2.7853527819234466), + BOOST_MATH_BIG_CONSTANT(T, 53, 0.74967671848044792), + BOOST_MATH_BIG_CONSTANT(T, 53, -0.00057069112416246805), + }; + // Maximum Deviation Found: 6.896e-018 + // Expected Error Term : -6.895e-018 + // Maximum Relative Change in Control Points : 8.497e-004 + static const T P_4_inf[] = { + 0.68947581948701249e-17L, + 0.49999999999998975L, + 1.0177274392923795L, + 2.498208511343429L, + 2.1921221359427595L, + 1.5897035272532764L, + 0.40154388356961734L, + }; + static const T Q_4_inf[] = { + 1.0L, + 1.7021215452463932L, + 4.4290431747556469L, + 2.9745631894384922L, + 2.3013614809773616L, + 0.28360399799075752L, + 0.022892987908906897L, + }; + + if(x <= 2) + { + return (offset + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x); + } + else if(x <= 4) + { + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_2_4, y) / tools::evaluate_polynomial(Q_2_4, y)) / x; + } + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_4_inf, y) / tools::evaluate_polynomial(Q_4_inf, y)) / x; +} + +template <class T, class Policy> +T trigamma_prec(T x, const mpl::int_<64>*, const Policy&) +{ + // Max error in interpolated form: 1.178e-020 + static const T offset_1_2 = BOOST_MATH_BIG_CONSTANT(T, 64, 2.109325408935546875); + static const T P_1_2[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -1.10932535608960258341), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.18793841543017129052), + BOOST_MATH_BIG_CONSTANT(T, 64, -4.63865531898487734531), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.919832884430500908047), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.68074038333180423012), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.21172611429185622377), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635673503366427284), + }; + static const T Q_1_2[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.77521119359546982995), + BOOST_MATH_BIG_CONSTANT(T, 64, 5.664338024578956321), + BOOST_MATH_BIG_CONSTANT(T, 64, 4.25995134879278028361), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.62956638448940402182), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.259635512844691089868), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.629642219810618032207e-8), + }; + // Max error in interpolated form: 3.912e-020 + static const T P_2_8[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -0.387540035162952880976e-11), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000276430504), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.21926880986360957306), + BOOST_MATH_BIG_CONSTANT(T, 64, 10.2550347708483445775), + BOOST_MATH_BIG_CONSTANT(T, 64, 18.9002075150709144043), + BOOST_MATH_BIG_CONSTANT(T, 64, 21.0357215832399705625), + BOOST_MATH_BIG_CONSTANT(T, 64, 13.4346512182925923978), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.98656291026448279118), + }; + static const T Q_2_8[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 64, 6.10520430478613667724), + BOOST_MATH_BIG_CONSTANT(T, 64, 18.475001060603645512), + BOOST_MATH_BIG_CONSTANT(T, 64, 31.7087534567758405638), + BOOST_MATH_BIG_CONSTANT(T, 64, 31.908814523890465398), + BOOST_MATH_BIG_CONSTANT(T, 64, 17.4175479039227084798), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.98749106958394941276), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.000115917322224411128566), + }; + // Maximum Deviation Found: 2.635e-020 + // Expected Error Term : 2.635e-020 + // Maximum Relative Change in Control Points : 1.791e-003 + static const T P_8_inf[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, -0.263527875092466899848e-19), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.500000000000000058145), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0730121433777364138677), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.94505878379957149534), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.0517092358874932620529), + BOOST_MATH_BIG_CONSTANT(T, 64, 1.07995383547483921121), + }; + static const T Q_8_inf[] = { + BOOST_MATH_BIG_CONSTANT(T, 64, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.187309046577818095504), + BOOST_MATH_BIG_CONSTANT(T, 64, 3.95255391645238842975), + BOOST_MATH_BIG_CONSTANT(T, 64, -1.14743283327078949087), + BOOST_MATH_BIG_CONSTANT(T, 64, 2.52989799376344914499), + BOOST_MATH_BIG_CONSTANT(T, 64, -0.627414303172402506396), + BOOST_MATH_BIG_CONSTANT(T, 64, 0.141554248216425512536), + }; + + if(x <= 2) + { + return (offset_1_2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x); + } + else if(x <= 8) + { + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_2_8, y) / tools::evaluate_polynomial(Q_2_8, y)) / x; + } + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_8_inf, y) / tools::evaluate_polynomial(Q_8_inf, y)) / x; +} + +template <class T, class Policy> +T trigamma_prec(T x, const mpl::int_<113>*, const Policy&) +{ + // Max error in interpolated form: 1.916e-035 + + static const T P_1_2[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.999999999999999082554457936871832533), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.71237311120865266379041700054847734), + BOOST_MATH_BIG_CONSTANT(T, 113, -7.94125711970499027763789342500817316), + BOOST_MATH_BIG_CONSTANT(T, 113, -5.74657746697664735258222071695644535), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.404213349456398905981223965160595687), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.47877781178642876561595890095758896), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.07714151702455125992166949812126433), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.858877899162360138844032265418028567), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.20499222604410032375789018837922397), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0272103140348194747360175268778415049), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0015764849020876949848954081173520686), + }; + static const T Q_1_2[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.71237311120863419878375031457715223), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.58619118655339853449127952145877467), + BOOST_MATH_BIG_CONSTANT(T, 113, 11.0940067269829372437561421279054968), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.09075424749327792073276309969037885), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.87705890159891405185343806884451286), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.22758678701914477836330837816976782), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.249092040606385004109672077814668716), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0295750413900655597027079600025569048), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157648490200498142247694709728858139), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.161264050344059471721062360645432809e-14), + }; + + // Max error in interpolated form: 8.958e-035 + static const T P_2_4[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -2.55843734739907925764326773972215085), + BOOST_MATH_BIG_CONSTANT(T, 113, -12.2830208240542011967952466273455887), + BOOST_MATH_BIG_CONSTANT(T, 113, -23.9195022162767993526575786066414403), + BOOST_MATH_BIG_CONSTANT(T, 113, -24.9256431504823483094158828285470862), + BOOST_MATH_BIG_CONSTANT(T, 113, -14.7979122765478779075108064826412285), + BOOST_MATH_BIG_CONSTANT(T, 113, -4.46654453928610666393276765059122272), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0191439033405649675717082465687845002), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.515412052554351265708917209749037352), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.195378348786064304378247325360320038), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0334761282624174313035014426794245393), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.002373665205942206348500250056602687), + }; + static const T Q_2_4[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.80098558454419907830670928248659245), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.99220727843170133895059300223445265), + BOOST_MATH_BIG_CONSTANT(T, 113, 11.8896146167631330735386697123464976), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.96613256683809091593793565879092581), + BOOST_MATH_BIG_CONSTANT(T, 113, 4.47254136149624110878909334574485751), + BOOST_MATH_BIG_CONSTANT(T, 113, 1.48600982028196527372434773913633152), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.319570735766764237068541501137990078), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0407358345787680953107374215319322066), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.00237366520593271641375755486420859837), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.239554887903526152679337256236302116e-15), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.294749244740618656265237072002026314e-17), + }; + + static const T y_offset_2_4 = BOOST_MATH_BIG_CONSTANT(T, 113, 3.558437347412109375); + + // Max error in interpolated form: 4.319e-035 + static const T P_4_8[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.166626112697021464248967707021688845e-16), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.499999999999997739552090249208808197), + BOOST_MATH_BIG_CONSTANT(T, 113, 6.40270945019053817915772473771553187), + BOOST_MATH_BIG_CONSTANT(T, 113, 41.3833374155000608013677627389343329), + BOOST_MATH_BIG_CONSTANT(T, 113, 166.803341854562809335667241074035245), + BOOST_MATH_BIG_CONSTANT(T, 113, 453.39964786925369319960722793414521), + BOOST_MATH_BIG_CONSTANT(T, 113, 851.153712317697055375935433362983944), + BOOST_MATH_BIG_CONSTANT(T, 113, 1097.70657567285059133109286478004458), + BOOST_MATH_BIG_CONSTANT(T, 113, 938.431232478455316020076349367632922), + BOOST_MATH_BIG_CONSTANT(T, 113, 487.268001604651932322080970189930074), + BOOST_MATH_BIG_CONSTANT(T, 113, 119.953445242335730062471193124820659), + }; + static const T Q_4_8[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 113, 12.4720855670474488978638945855932398), + BOOST_MATH_BIG_CONSTANT(T, 113, 78.6093129753298570701376952709727391), + BOOST_MATH_BIG_CONSTANT(T, 113, 307.470246050318322489781182863190127), + BOOST_MATH_BIG_CONSTANT(T, 113, 805.140686101151538537565264188630079), + BOOST_MATH_BIG_CONSTANT(T, 113, 1439.12019760292146454787601409644413), + BOOST_MATH_BIG_CONSTANT(T, 113, 1735.6105285756048831268586001383127), + BOOST_MATH_BIG_CONSTANT(T, 113, 1348.32500712856328019355198611280536), + BOOST_MATH_BIG_CONSTANT(T, 113, 607.225985860570846699704222144650563), + BOOST_MATH_BIG_CONSTANT(T, 113, 119.952317857277045332558673164517227), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.000140165918355036060868680809129436084), + }; + + // Maximum Deviation Found: 2.867e-035 + // Expected Error Term : 2.866e-035 + // Maximum Relative Change in Control Points : 2.662e-004 + static const T P_8_16[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, -0.184828315274146610610872315609837439e-19), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000004122475157735807738), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.02533865247313349284875558880415875), + BOOST_MATH_BIG_CONSTANT(T, 113, 13.5995927517457371243039532492642734), + BOOST_MATH_BIG_CONSTANT(T, 113, 35.3132224283087906757037999452941588), + BOOST_MATH_BIG_CONSTANT(T, 113, 67.1639424550714159157603179911505619), + BOOST_MATH_BIG_CONSTANT(T, 113, 83.5767733658513967581959839367419891), + BOOST_MATH_BIG_CONSTANT(T, 113, 71.073491212235705900866411319363501), + BOOST_MATH_BIG_CONSTANT(T, 113, 35.8621515614725564575893663483998663), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.72152231639983491987779743154333318), + }; + static const T Q_8_16[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 113, 5.71734397161293452310624822415866372), + BOOST_MATH_BIG_CONSTANT(T, 113, 25.293404179620438179337103263274815), + BOOST_MATH_BIG_CONSTANT(T, 113, 62.2619767967468199111077640625328469), + BOOST_MATH_BIG_CONSTANT(T, 113, 113.955048909238993473389714972250235), + BOOST_MATH_BIG_CONSTANT(T, 113, 130.807138328938966981862203944329408), + BOOST_MATH_BIG_CONSTANT(T, 113, 102.423146902337654110717764213057753), + BOOST_MATH_BIG_CONSTANT(T, 113, 44.0424772805245202514468199602123565), + BOOST_MATH_BIG_CONSTANT(T, 113, 8.89898032477904072082994913461386099), + BOOST_MATH_BIG_CONSTANT(T, 113, -0.0296627336872039988632793863671456398), + }; + // Maximum Deviation Found: 1.079e-035 + // Expected Error Term : -1.079e-035 + // Maximum Relative Change in Control Points : 7.884e-003 + static const T P_16_inf[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 0.0), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.500000000000000000000000000000087317), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.345625669885456215194494735902663968), + BOOST_MATH_BIG_CONSTANT(T, 113, 9.62895499360842232127552650044647769), + BOOST_MATH_BIG_CONSTANT(T, 113, 3.5936085382439026269301003761320812), + BOOST_MATH_BIG_CONSTANT(T, 113, 49.459599118438883265036646019410669), + BOOST_MATH_BIG_CONSTANT(T, 113, 7.77519237321893917784735690560496607), + BOOST_MATH_BIG_CONSTANT(T, 113, 74.4536074488178075948642351179304121), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.75209340397069050436806159297952699), + BOOST_MATH_BIG_CONSTANT(T, 113, 23.9292359711471667884504840186561598), + }; + static const T Q_16_inf[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.0), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.357918006437579097055656138920742037), + BOOST_MATH_BIG_CONSTANT(T, 113, 19.1386039850709849435325005484512944), + BOOST_MATH_BIG_CONSTANT(T, 113, 0.874349081464143606016221431763364517), + BOOST_MATH_BIG_CONSTANT(T, 113, 98.6516097434855572678195488061432509), + BOOST_MATH_BIG_CONSTANT(T, 113, -16.1051972833382893468655223662534306), + BOOST_MATH_BIG_CONSTANT(T, 113, 154.316860216253720989145047141653727), + BOOST_MATH_BIG_CONSTANT(T, 113, -40.2026880424378986053105969312264534), + BOOST_MATH_BIG_CONSTANT(T, 113, 60.1679136674264778074736441126810223), + BOOST_MATH_BIG_CONSTANT(T, 113, -13.3414844622256422644504472438320114), + BOOST_MATH_BIG_CONSTANT(T, 113, 2.53795636200649908779512969030363442), + }; + + if(x <= 2) + { + return (2 + boost::math::tools::evaluate_polynomial(P_1_2, x) / tools::evaluate_polynomial(Q_1_2, x)) / (x * x); + } + else if(x <= 4) + { + return (y_offset_2_4 + boost::math::tools::evaluate_polynomial(P_2_4, x) / tools::evaluate_polynomial(Q_2_4, x)) / (x * x); + } + else if(x <= 8) + { + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_4_8, y) / tools::evaluate_polynomial(Q_4_8, y)) / x; + } + else if(x <= 16) + { + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_8_16, y) / tools::evaluate_polynomial(Q_8_16, y)) / x; + } + T y = 1 / x; + return (1 + tools::evaluate_polynomial(P_16_inf, y) / tools::evaluate_polynomial(Q_16_inf, y)) / x; +} + +template <class T, class Tag, class Policy> +T trigamma_imp(T x, const Tag* t, const Policy& pol) +{ + // + // This handles reflection of negative arguments, and all our + // error handling, then forwards to the T-specific approximation. + // + BOOST_MATH_STD_USING // ADL of std functions. + + T result = 0; + // + // Check for negative arguments and use reflection: + // + if(x <= 0) + { + // Reflect: + T z = 1 - x; + // Argument reduction for tan: + if(floor(x) == x) + { + return policies::raise_pole_error<T>("boost::math::trigamma<%1%>(%1%)", 0, (1-x), pol); + } + T s = fabs(x) < fabs(z) ? boost::math::sin_pi(x, pol) : boost::math::sin_pi(z, pol); + return -trigamma_imp(z, t, pol) + boost::math::pow<2>(constants::pi<T>()) / (s * s); + } + if(x < 1) + { + result = 1 / (x * x); + x += 1; + } + return result + trigamma_prec(x, t, pol); +} + +template <class T, class Policy> +T trigamma_imp(T x, const mpl::int_<0>*, const Policy& pol) +{ + return polygamma_imp(1, x, pol); +} +// +// Initializer: ensure all our constants are initialized prior to the first call of main: +// +template <class T, class Policy> +struct trigamma_initializer +{ + struct init + { + init() + { + typedef typename policies::precision<T, Policy>::type precision_type; + do_init(mpl::bool_<precision_type::value && (precision_type::value <= 113)>()); + } + void do_init(const mpl::true_&) + { + boost::math::trigamma(T(2.5), Policy()); + } + void do_init(const mpl::false_&){} + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename trigamma_initializer<T, Policy>::init trigamma_initializer<T, Policy>::initializer; + +} // namespace detail + +template <class T, class Policy> +inline typename tools::promote_args<T>::type + trigamma(T x, const Policy&) +{ + typedef typename tools::promote_args<T>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + typedef typename policies::precision<T, Policy>::type precision_type; + typedef typename mpl::if_< + mpl::or_< + mpl::less_equal<precision_type, mpl::int_<0> >, + mpl::greater<precision_type, mpl::int_<114> > + >, + mpl::int_<0>, + typename mpl::if_< + mpl::less<precision_type, mpl::int_<54> >, + mpl::int_<53>, + typename mpl::if_< + mpl::less<precision_type, mpl::int_<65> >, + mpl::int_<64>, + mpl::int_<113> + >::type + >::type + >::type tag_type; + + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + // Force initialization of constants: + detail::trigamma_initializer<value_type, forwarding_policy>::force_instantiate(); + + return policies::checked_narrowing_cast<result_type, Policy>(detail::trigamma_imp( + static_cast<value_type>(x), + static_cast<const tag_type*>(0), forwarding_policy()), "boost::math::trigamma<%1%>(%1%)"); +} + +template <class T> +inline typename tools::promote_args<T>::type + trigamma(T x) +{ + return trigamma(x, policies::policy<>()); +} + +} // namespace math +} // namespace boost +#endif + diff --git a/boost/math/special_functions/zeta.hpp b/boost/math/special_functions/zeta.hpp index b176f20176..1ba282f4d1 100644 --- a/boost/math/special_functions/zeta.hpp +++ b/boost/math/special_functions/zeta.hpp @@ -1,4 +1,4 @@ -// Copyright John Maddock 2007. +// Copyright John Maddock 2007, 2014. // Use, modification and distribution are subject to 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) @@ -16,6 +16,7 @@ #include <boost/math/tools/big_constant.hpp> #include <boost/math/policies/error_handling.hpp> #include <boost/math/special_functions/gamma.hpp> +#include <boost/math/special_functions/factorials.hpp> #include <boost/math/special_functions/sin_pi.hpp> namespace boost{ namespace math{ namespace detail{ @@ -167,6 +168,8 @@ T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<0>&) { BOOST_MATH_STD_USING T result; + if(s >= policies::digits<T, Policy>()) + return 1; result = zeta_polynomial_series(s, sc, pol); #if 0 // Old code archived for future reference: @@ -863,6 +866,34 @@ T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<113>&) return result; } +template <class T, class Policy> +T zeta_imp_odd_integer(int s, const T&, const Policy&, const mpl::true_&) +{ + static const T results[] = { + BOOST_MATH_BIG_CONSTANT(T, 113, 1.2020569031595942853997381615114500), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0369277551433699263313654864570342), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0083492773819228268397975498497968), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0020083928260822144178527692324121), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0004941886041194645587022825264699), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0001227133475784891467518365263574), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000305882363070204935517285106451), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000076371976378997622736002935630), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000019082127165539389256569577951), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000004769329867878064631167196044), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000001192199259653110730677887189), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000298035035146522801860637051), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000074507117898354294919810042), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000018626597235130490064039099), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000004656629065033784072989233), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000001164155017270051977592974), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000291038504449709968692943), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000072759598350574810145209), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000018189896503070659475848), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000004547473783042154026799), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000001136868407680227849349), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000284217097688930185546), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000071054273952108527129), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000017763568435791203275), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000004440892103143813364), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000001110223025141066134), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000277555756213612417), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000069388939045441537), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000017347234760475766), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000004336808690020650), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000001084202172494241), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000271050543122347), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000067762635780452), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000016940658945098), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000004235164736273), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000001058791184068), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000264697796017), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000066174449004), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000016543612251), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000004135903063), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000001033975766), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000258493941), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000064623485), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000016155871), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000004038968), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000001009742), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000252435), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000063109), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000015777), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000003944), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000986), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000247), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000062), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000015), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000004), BOOST_MATH_BIG_CONSTANT(T, 113, 1.0000000000000000000000000000000001), + }; + return s > 113 ? 1 : results[(s - 3) / 2]; +} + +template <class T, class Policy> +T zeta_imp_odd_integer(int s, const T& sc, const Policy& pol, const mpl::false_&) +{ + static bool is_init = false; + static T results[50] = {}; + if(!is_init) + { + is_init = true; + for(int k = 0; k < sizeof(results) / sizeof(results[0]); ++k) + { + T arg = k * 2 + 3; + T c_arg = 1 - arg; + results[k] = zeta_polynomial_series(arg, c_arg, pol); + } + } + int index = (s - 3) / 2; + return index >= sizeof(results) / sizeof(results[0]) ? zeta_polynomial_series(T(s), sc, pol): results[index]; +} + template <class T, class Policy, class Tag> T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag) { @@ -874,6 +905,45 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag) "Evaluation of zeta function at pole %1%", s, pol); T result; + // + // Trivial case: + // + if(s > policies::digits<T, Policy>()) + return 1; + // + // Start by seeing if we have a simple closed form: + // + if(floor(s) == s) + { + try + { + int v = itrunc(s); + if(v == s) + { + if(v < 0) + { + if(((-v) & 1) == 0) + return 0; + int n = (-v + 1) / 2; + if(n <= boost::math::max_bernoulli_b2n<T>::value) + return T((-v & 1) ? -1 : 1) * boost::math::unchecked_bernoulli_b2n<T>(n) / (1 - v); + } + else if((v & 1) == 0) + { + if(((v / 2) <= boost::math::max_bernoulli_b2n<T>::value) && (v <= boost::math::max_factorial<T>::value)) + return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * pow(constants::pi<T, Policy>(), v) * + boost::math::unchecked_bernoulli_b2n<T>(v / 2) / boost::math::unchecked_factorial<T>(v); + return T(((v / 2 - 1) & 1) ? -1 : 1) * ldexp(T(1), v - 1) * pow(constants::pi<T, Policy>(), v) * + boost::math::bernoulli_b2n<T>(v / 2) / boost::math::factorial<T>(v); + } + else + return zeta_imp_odd_integer(v, sc, pol, mpl::bool_<(Tag::value <= 113) && Tag::value>()); + } + } + catch(const boost::math::rounding_error&){} // Just fall through, s is too large to round + catch(const std::overflow_error&){} + } + if(fabs(s) < tools::root_epsilon<T>()) { result = -0.5f - constants::log_root_two_pi<T, Policy>() * s; @@ -922,8 +992,8 @@ struct zeta_initializer { do_init(tag()); } - static void do_init(const mpl::int_<0>&){} - static void do_init(const mpl::int_<53>&){} + static void do_init(const mpl::int_<0>&){ boost::math::zeta(static_cast<T>(5), Policy()); } + static void do_init(const mpl::int_<53>&){ boost::math::zeta(static_cast<T>(5), Policy()); } static void do_init(const mpl::int_<64>&) { boost::math::zeta(static_cast<T>(0.5), Policy()); @@ -932,6 +1002,8 @@ struct zeta_initializer boost::math::zeta(static_cast<T>(6.5), Policy()); boost::math::zeta(static_cast<T>(14.5), Policy()); boost::math::zeta(static_cast<T>(40.5), Policy()); + + boost::math::zeta(static_cast<T>(5), Policy()); } static void do_init(const mpl::int_<113>&) { @@ -941,8 +1013,10 @@ struct zeta_initializer boost::math::zeta(static_cast<T>(5.5), Policy()); boost::math::zeta(static_cast<T>(9.5), Policy()); boost::math::zeta(static_cast<T>(16.5), Policy()); - boost::math::zeta(static_cast<T>(25), Policy()); - boost::math::zeta(static_cast<T>(70), Policy()); + boost::math::zeta(static_cast<T>(25.5), Policy()); + boost::math::zeta(static_cast<T>(70.5), Policy()); + + boost::math::zeta(static_cast<T>(5), Policy()); } void force_instantiate()const{} }; diff --git a/boost/math/tools/big_constant.hpp b/boost/math/tools/big_constant.hpp index 423cc1b631..3895638891 100644 --- a/boost/math/tools/big_constant.hpp +++ b/boost/math/tools/big_constant.hpp @@ -12,31 +12,49 @@ #include <boost/lexical_cast.hpp> #endif #include <boost/type_traits/is_convertible.hpp> -#include <boost/math/cstdfloat/cstdfloat_types.hpp> namespace boost{ namespace math{ namespace tools{ template <class T> -inline BOOST_CONSTEXPR_OR_CONST T make_big_value(boost::floatmax_t v, const char*, mpl::true_ const&, mpl::false_ const&) +struct numeric_traits : public std::numeric_limits< T > {}; + +#ifdef BOOST_MATH_USE_FLOAT128 +typedef __float128 largest_float; +#define BOOST_MATH_LARGEST_FLOAT_C(x) x##Q +template <> +struct numeric_traits<__float128> +{ + static const int digits = 113; + static const int digits10 = 34; + static const int max_exponent = 16384; + static const bool is_specialized = true; +}; +#else +typedef long double largest_float; +#define BOOST_MATH_LARGEST_FLOAT_C(x) x##L +#endif + +template <class T> +inline BOOST_CONSTEXPR_OR_CONST T make_big_value(largest_float v, const char*, mpl::true_ const&, mpl::false_ const&) { return static_cast<T>(v); } template <class T> -inline BOOST_CONSTEXPR_OR_CONST T make_big_value(boost::floatmax_t v, const char*, mpl::true_ const&, mpl::true_ const&) +inline BOOST_CONSTEXPR_OR_CONST T make_big_value(largest_float v, const char*, mpl::true_ const&, mpl::true_ const&) { return static_cast<T>(v); } #ifndef BOOST_MATH_NO_LEXICAL_CAST template <class T> -inline T make_big_value(boost::floatmax_t, const char* s, mpl::false_ const&, mpl::false_ const&) +inline T make_big_value(largest_float, const char* s, mpl::false_ const&, mpl::false_ const&) { return boost::lexical_cast<T>(s); } #endif template <class T> -inline BOOST_CONSTEXPR const char* make_big_value(boost::floatmax_t, const char* s, mpl::false_ const&, mpl::true_ const&) +inline BOOST_CONSTEXPR const char* make_big_value(largest_float, const char* s, mpl::false_ const&, mpl::true_ const&) { return s; } @@ -46,20 +64,20 @@ inline BOOST_CONSTEXPR const char* make_big_value(boost::floatmax_t, const char* // #define BOOST_MATH_BIG_CONSTANT(T, D, x)\ boost::math::tools::make_big_value<T>(\ - BOOST_FLOATMAX_C(x), \ + BOOST_MATH_LARGEST_FLOAT_C(x), \ BOOST_STRINGIZE(x), \ - mpl::bool_< (is_convertible<boost::floatmax_t, T>::value) && \ - ((D <= std::numeric_limits<boost::floatmax_t>::digits) \ + mpl::bool_< (is_convertible<boost::math::tools::largest_float, T>::value) && \ + ((D <= boost::math::tools::numeric_traits<boost::math::tools::largest_float>::digits) \ || is_floating_point<T>::value \ - || (std::numeric_limits<T>::is_specialized && \ - (std::numeric_limits<T>::digits10 <= std::numeric_limits<boost::floatmax_t>::digits10))) >(), \ + || (boost::math::tools::numeric_traits<T>::is_specialized && \ + (boost::math::tools::numeric_traits<T>::digits10 <= boost::math::tools::numeric_traits<boost::math::tools::largest_float>::digits10))) >(), \ boost::is_convertible<const char*, T>()) // // For constants too huge for any conceivable long double (and which generate compiler errors if we try and declare them as such): // #define BOOST_MATH_HUGE_CONSTANT(T, D, x)\ boost::math::tools::make_big_value<T>(0.0L, BOOST_STRINGIZE(x), \ - mpl::bool_<is_floating_point<T>::value || (std::numeric_limits<T>::is_specialized && std::numeric_limits<T>::max_exponent <= std::numeric_limits<long double>::max_exponent && std::numeric_limits<T>::digits <= std::numeric_limits<long double>::digits)>(), \ + mpl::bool_<is_floating_point<T>::value || (boost::math::tools::numeric_traits<T>::is_specialized && boost::math::tools::numeric_traits<T>::max_exponent <= boost::math::tools::numeric_traits<boost::math::tools::largest_float>::max_exponent && boost::math::tools::numeric_traits<T>::digits <= boost::math::tools::numeric_traits<boost::math::tools::largest_float>::digits)>(), \ boost::is_convertible<const char*, T>()) }}} // namespaces diff --git a/boost/math/tools/config.hpp b/boost/math/tools/config.hpp index 4ec5768aaf..23cba0c9a9 100644 --- a/boost/math/tools/config.hpp +++ b/boost/math/tools/config.hpp @@ -50,6 +50,10 @@ // are disabled for now. (JM 2012). # define BOOST_MATH_NO_REAL_CONCEPT_TESTS #endif +#ifdef sun +// Any use of __float128 in program startup code causes a segfault (tested JM 2015, Solaris 11). +# define BOOST_MATH_DISABLE_FLOAT128 +#endif #if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106)) && !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) // // Darwin's rather strange "double double" is rather hard to diff --git a/boost/math/tools/precision.hpp b/boost/math/tools/precision.hpp index 49e653d6a3..ed146c458f 100644 --- a/boost/math/tools/precision.hpp +++ b/boost/math/tools/precision.hpp @@ -134,7 +134,12 @@ inline T log_max_value(const mpl::int_<0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_T BOOST_ASSERT(::std::numeric_limits<T>::is_specialized); #endif BOOST_MATH_STD_USING +#ifdef __SUNPRO_CC + static const T m = (std::numeric_limits<T>::max)(); + static const T val = log(m); +#else static const T val = log((std::numeric_limits<T>::max)()); +#endif return val; } @@ -147,7 +152,12 @@ inline T log_min_value(const mpl::int_<0>& BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_T BOOST_ASSERT(::std::numeric_limits<T>::is_specialized); #endif BOOST_MATH_STD_USING +#ifdef __SUNPRO_CC + static const T m = (std::numeric_limits<T>::min)(); + static const T val = log(m); +#else static const T val = log((std::numeric_limits<T>::min)()); +#endif return val; } |