summaryrefslogtreecommitdiff
path: root/boost/math/special_functions
diff options
context:
space:
mode:
Diffstat (limited to 'boost/math/special_functions')
-rw-r--r--boost/math/special_functions/airy.hpp4
-rw-r--r--boost/math/special_functions/bernoulli.hpp2
-rw-r--r--boost/math/special_functions/bessel.hpp12
-rw-r--r--boost/math/special_functions/beta.hpp57
-rw-r--r--boost/math/special_functions/binomial.hpp6
-rw-r--r--boost/math/special_functions/cos_pi.hpp23
-rw-r--r--boost/math/special_functions/detail/erf_inv.hpp4
-rw-r--r--boost/math/special_functions/detail/fp_traits.hpp3
-rw-r--r--boost/math/special_functions/detail/polygamma.hpp538
-rw-r--r--boost/math/special_functions/digamma.hpp176
-rw-r--r--boost/math/special_functions/ellint_1.hpp14
-rw-r--r--boost/math/special_functions/ellint_2.hpp40
-rw-r--r--boost/math/special_functions/ellint_3.hpp406
-rw-r--r--boost/math/special_functions/ellint_d.hpp180
-rw-r--r--boost/math/special_functions/ellint_rc.hpp63
-rw-r--r--boost/math/special_functions/ellint_rd.hpp217
-rw-r--r--boost/math/special_functions/ellint_rf.hpp182
-rw-r--r--boost/math/special_functions/ellint_rg.hpp136
-rw-r--r--boost/math/special_functions/ellint_rj.hpp356
-rw-r--r--boost/math/special_functions/factorials.hpp2
-rw-r--r--boost/math/special_functions/fpclassify.hpp6
-rw-r--r--boost/math/special_functions/gamma.hpp42
-rw-r--r--boost/math/special_functions/heuman_lambda.hpp87
-rw-r--r--boost/math/special_functions/jacobi_zeta.hpp74
-rw-r--r--boost/math/special_functions/legendre.hpp4
-rw-r--r--boost/math/special_functions/math_fwd.hpp78
-rw-r--r--boost/math/special_functions/polygamma.hpp83
-rw-r--r--boost/math/special_functions/sign.hpp4
-rw-r--r--boost/math/special_functions/sin_pi.hpp11
-rw-r--r--boost/math/special_functions/trigamma.hpp469
-rw-r--r--boost/math/special_functions/zeta.hpp84
31 files changed, 2820 insertions, 543 deletions
diff --git a/boost/math/special_functions/airy.hpp b/boost/math/special_functions/airy.hpp
index 82167dc..e84a705 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 2c2ccd5..5f2b7e1 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 438b763..eab723d8 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 8486b82..98d8f7f 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 8fa2e26..9a24fc1 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 18d06c0..669a2c8 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 77aa72f..35072d5 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 63ebf11..09dc516 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 0000000..4ef503b
--- /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 785cd75..718eaf9 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 da16bc6..62a0bf3 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 72caf3e..f4f65cc 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 ac7e68c..9ab0f83 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 0000000..5bd065d
--- /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 5f6d5c6..846c752 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 61014d3..03b73b1 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 ac57257..a8a7b4b 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 0000000..bb5b7c3
--- /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 8a242f0..ac39bed 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 de24642..e36a098 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 40f6e14..8e75fae 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 b6b4c57..1db1e4c 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 0000000..6389443
--- /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 0000000..a3fa547
--- /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 79e9756..1a2ef5d 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 e952dcd..96f6072 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 0000000..6b7815d
--- /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 f5c562d..3324c90 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 16aed51..ae6b3e7 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 0000000..6fccb36
--- /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 b176f20..1ba282f 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{}
};