summaryrefslogtreecommitdiff
path: root/boost/math/special_functions/detail
diff options
context:
space:
mode:
authorChanho Park <chanho61.park@samsung.com>2014-12-11 09:55:56 (GMT)
committerChanho Park <chanho61.park@samsung.com>2014-12-11 09:55:56 (GMT)
commit08c1e93fa36a49f49325a07fe91ff92c964c2b6c (patch)
tree7a7053ceb8874b28ec4b868d4c49b500008a102e /boost/math/special_functions/detail
parentbb4dd8289b351fae6b55e303f189127a394a1edd (diff)
downloadboost-08c1e93fa36a49f49325a07fe91ff92c964c2b6c.zip
boost-08c1e93fa36a49f49325a07fe91ff92c964c2b6c.tar.gz
boost-08c1e93fa36a49f49325a07fe91ff92c964c2b6c.tar.bz2
Imported Upstream version 1.57.0upstream/1.57.0
Diffstat (limited to 'boost/math/special_functions/detail')
-rw-r--r--boost/math/special_functions/detail/airy_ai_bi_zero.hpp160
-rw-r--r--boost/math/special_functions/detail/bernoulli_details.hpp653
-rw-r--r--boost/math/special_functions/detail/bessel_derivatives_linear.hpp75
-rw-r--r--boost/math/special_functions/detail/bessel_i0.hpp5
-rw-r--r--boost/math/special_functions/detail/bessel_i1.hpp5
-rw-r--r--boost/math/special_functions/detail/bessel_ik.hpp23
-rw-r--r--boost/math/special_functions/detail/bessel_j0.hpp16
-rw-r--r--boost/math/special_functions/detail/bessel_j1.hpp17
-rw-r--r--boost/math/special_functions/detail/bessel_jn.hpp13
-rw-r--r--boost/math/special_functions/detail/bessel_jy.hpp1042
-rw-r--r--boost/math/special_functions/detail/bessel_jy_asym.hpp176
-rw-r--r--boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp141
-rw-r--r--boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp220
-rw-r--r--boost/math/special_functions/detail/bessel_jy_series.hpp8
-rw-r--r--boost/math/special_functions/detail/bessel_jy_zero.hpp617
-rw-r--r--boost/math/special_functions/detail/bessel_kn.hpp1
-rw-r--r--boost/math/special_functions/detail/bessel_y0.hpp17
-rw-r--r--boost/math/special_functions/detail/bessel_y1.hpp16
-rw-r--r--boost/math/special_functions/detail/bessel_yn.hpp3
-rw-r--r--boost/math/special_functions/detail/erf_inv.hpp51
-rw-r--r--boost/math/special_functions/detail/fp_traits.hpp11
-rw-r--r--boost/math/special_functions/detail/gamma_inva.hpp8
-rw-r--r--boost/math/special_functions/detail/ibeta_inv_ab.hpp24
-rw-r--r--boost/math/special_functions/detail/ibeta_inverse.hpp97
-rw-r--r--boost/math/special_functions/detail/igamma_inverse.hpp22
-rw-r--r--boost/math/special_functions/detail/lanczos_sse2.hpp20
-rw-r--r--boost/math/special_functions/detail/lgamma_small.hpp12
-rw-r--r--boost/math/special_functions/detail/round_fwd.hpp21
-rw-r--r--boost/math/special_functions/detail/t_distribution_inv.hpp19
-rw-r--r--boost/math/special_functions/detail/unchecked_bernoulli.hpp700
-rw-r--r--boost/math/special_functions/detail/unchecked_factorial.hpp210
31 files changed, 3636 insertions, 767 deletions
diff --git a/boost/math/special_functions/detail/airy_ai_bi_zero.hpp b/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
new file mode 100644
index 0000000..dbb7388
--- /dev/null
+++ b/boost/math/special_functions/detail/airy_ai_bi_zero.hpp
@@ -0,0 +1,160 @@
+// Copyright (c) 2013 Christopher Kormanyos
+// 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)
+//
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+//
+// This header contains implementation details for estimating the zeros
+// of the Airy functions airy_ai and airy_bi on the negative real axis.
+//
+#ifndef _AIRY_AI_BI_ZERO_2013_01_20_HPP_
+ #define _AIRY_AI_BI_ZERO_2013_01_20_HPP_
+
+ #include <boost/math/constants/constants.hpp>
+ #include <boost/math/special_functions/cbrt.hpp>
+
+ namespace boost { namespace math {
+ namespace detail
+ {
+ // Forward declarations of the needed Airy function implementations.
+ template <class T, class Policy>
+ T airy_ai_imp(T x, const Policy& pol);
+ template <class T, class Policy>
+ T airy_bi_imp(T x, const Policy& pol);
+ template <class T, class Policy>
+ T airy_ai_prime_imp(T x, const Policy& pol);
+ template <class T, class Policy>
+ T airy_bi_prime_imp(T x, const Policy& pol);
+
+ namespace airy_zero
+ {
+ template<class T>
+ T equation_as_10_4_105(const T& z)
+ {
+ const T one_over_z (T(1) / z);
+ const T one_over_z_squared(one_over_z * one_over_z);
+
+ const T z_pow_third (boost::math::cbrt(z));
+ const T z_pow_two_thirds(z_pow_third * z_pow_third);
+
+ // Implement the top line of Eq. 10.4.105.
+ const T fz(z_pow_two_thirds * ((((( + (T(162375596875.0) / 334430208UL)
+ * one_over_z_squared - ( T(108056875.0) / 6967296UL))
+ * one_over_z_squared + ( T(77125UL) / 82944UL))
+ * one_over_z_squared - ( T(5U) / 36U))
+ * one_over_z_squared + ( T(5U) / 48U))
+ * one_over_z_squared + (1)));
+
+ return fz;
+ }
+
+ namespace airy_ai_zero_detail
+ {
+ template<class T>
+ T initial_guess(const int m)
+ {
+ T guess;
+
+ switch(m)
+ {
+ case 0: { guess = T(0); break; }
+ case 1: { guess = T(-2.33810741045976703849); break; }
+ case 2: { guess = T(-4.08794944413097061664); break; }
+ case 3: { guess = T(-5.52055982809555105913); break; }
+ case 4: { guess = T(-6.78670809007175899878); break; }
+ case 5: { guess = T(-7.94413358712085312314); break; }
+ case 6: { guess = T(-9.02265085334098038016); break; }
+ case 7: { guess = T(-10.0401743415580859306); break; }
+ case 8: { guess = T(-11.0085243037332628932); break; }
+ case 9: { guess = T(-11.9360155632362625170); break; }
+ case 10:{ guess = T(-12.8287767528657572004); break; }
+ default:
+ {
+ const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 1)) / 8);
+ guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+ break;
+ }
+ }
+
+ return guess;
+ }
+
+ template<class T, class Policy>
+ class function_object_ai_and_ai_prime
+ {
+ public:
+ function_object_ai_and_ai_prime(const Policy pol) : my_pol(pol) { }
+
+ boost::math::tuple<T, T> operator()(const T& x) const
+ {
+ // Return a tuple containing both Ai(x) and Ai'(x).
+ return boost::math::make_tuple(
+ boost::math::detail::airy_ai_imp (x, my_pol),
+ boost::math::detail::airy_ai_prime_imp(x, my_pol));
+ }
+
+ private:
+ const Policy& my_pol;
+ const function_object_ai_and_ai_prime& operator=(const function_object_ai_and_ai_prime&);
+ };
+ } // namespace airy_ai_zero_detail
+
+ namespace airy_bi_zero_detail
+ {
+ template<class T>
+ T initial_guess(const int m)
+ {
+ T guess;
+
+ switch(m)
+ {
+ case 0: { guess = T(0); break; }
+ case 1: { guess = T(-1.17371322270912792492); break; }
+ case 2: { guess = T(-3.27109330283635271568); break; }
+ case 3: { guess = T(-4.83073784166201593267); break; }
+ case 4: { guess = T(-6.16985212831025125983); break; }
+ case 5: { guess = T(-7.37676207936776371360); break; }
+ case 6: { guess = T(-8.49194884650938801345); break; }
+ case 7: { guess = T(-9.53819437934623888663); break; }
+ case 8: { guess = T(-10.5299135067053579244); break; }
+ case 9: { guess = T(-11.4769535512787794379); break; }
+ case 10: { guess = T(-12.3864171385827387456); break; }
+ default:
+ {
+ const T t(((boost::math::constants::pi<T>() * 3) * ((T(m) * 4) - 3)) / 8);
+ guess = -boost::math::detail::airy_zero::equation_as_10_4_105(t);
+ break;
+ }
+ }
+
+ return guess;
+ }
+
+ template<class T, class Policy>
+ class function_object_bi_and_bi_prime
+ {
+ public:
+ function_object_bi_and_bi_prime(const Policy pol) : my_pol(pol) { }
+
+ boost::math::tuple<T, T> operator()(const T& x) const
+ {
+ // Return a tuple containing both Bi(x) and Bi'(x).
+ return boost::math::make_tuple(
+ boost::math::detail::airy_bi_imp (x, my_pol),
+ boost::math::detail::airy_bi_prime_imp(x, my_pol));
+ }
+
+ private:
+ const Policy& my_pol;
+ const function_object_bi_and_bi_prime& operator=(const function_object_bi_and_bi_prime&);
+ };
+ } // namespace airy_bi_zero_detail
+ } // namespace airy_zero
+ } // namespace detail
+ } // namespace math
+ } // namespaces boost
+
+#endif // _AIRY_AI_BI_ZERO_2013_01_20_HPP_
diff --git a/boost/math/special_functions/detail/bernoulli_details.hpp b/boost/math/special_functions/detail/bernoulli_details.hpp
new file mode 100644
index 0000000..f2d3c65
--- /dev/null
+++ b/boost/math/special_functions/detail/bernoulli_details.hpp
@@ -0,0 +1,653 @@
+///////////////////////////////////////////////////////////////////////////////
+// Copyright 2013 John Maddock
+// 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_MATH_BERNOULLI_DETAIL_HPP
+#define BOOST_MATH_BERNOULLI_DETAIL_HPP
+
+#include <boost/config.hpp>
+#include <boost/detail/lightweight_mutex.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <boost/math/tools/toms748_solve.hpp>
+
+#ifdef BOOST_HAS_THREADS
+
+#ifndef BOOST_NO_CXX11_HDR_ATOMIC
+# include <atomic>
+# define BOOST_MATH_ATOMIC_NS std
+#if ATOMIC_INT_LOCK_FREE == 2
+typedef std::atomic<int> atomic_counter_type;
+typedef int atomic_integer_type;
+#elif ATOMIC_SHORT_LOCK_FREE == 2
+typedef std::atomic<short> atomic_counter_type;
+typedef short atomic_integer_type;
+#elif ATOMIC_LONG_LOCK_FREE == 2
+typedef std::atomic<long> atomic_counter_type;
+typedef long atomic_integer_type;
+#elif ATOMIC_LLONG_LOCK_FREE == 2
+typedef std::atomic<long long> atomic_counter_type;
+typedef long long atomic_integer_type;
+#else
+# define BOOST_MATH_NO_ATOMIC_INT
+#endif
+
+#else // BOOST_NO_CXX11_HDR_ATOMIC
+//
+// We need Boost.Atomic, but on any platform that supports auto-linking we do
+// not need to link against a separate library:
+//
+#define BOOST_ATOMIC_NO_LIB
+#include <boost/atomic.hpp>
+# define BOOST_MATH_ATOMIC_NS boost
+
+namespace boost{ namespace math{ namespace detail{
+
+//
+// We need a type to use as an atomic counter:
+//
+#if BOOST_ATOMIC_INT_LOCK_FREE == 2
+typedef boost::atomic<int> atomic_counter_type;
+typedef int atomic_integer_type;
+#elif BOOST_ATOMIC_SHORT_LOCK_FREE == 2
+typedef boost::atomic<short> atomic_counter_type;
+typedef short atomic_integer_type;
+#elif BOOST_ATOMIC_LONG_LOCK_FREE == 2
+typedef boost::atomic<long> atomic_counter_type;
+typedef long atomic_integer_type;
+#elif BOOST_ATOMIC_LLONG_LOCK_FREE == 2
+typedef boost::atomic<long long> atomic_counter_type;
+typedef long long atomic_integer_type;
+#else
+# define BOOST_MATH_NO_ATOMIC_INT
+#endif
+
+}}} // namespaces
+
+#endif // BOOST_NO_CXX11_HDR_ATOMIC
+
+#endif // BOOST_HAS_THREADS
+
+namespace boost{ namespace math{ namespace detail{
+//
+// Asymptotic expansion for B2n due to
+// Luschny LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
+//
+template <class T, class Policy>
+T b2n_asymptotic(int n)
+{
+ BOOST_MATH_STD_USING
+ const T nx = static_cast<T>(n);
+ const T nx2(nx * nx);
+
+ const T approximate_log_of_bernoulli_bn =
+ ((boost::math::constants::half<T>() + nx) * log(nx))
+ + ((boost::math::constants::half<T>() - nx) * log(boost::math::constants::pi<T>()))
+ + (((T(3) / 2) - nx) * boost::math::constants::ln_two<T>())
+ + ((nx * (T(2) - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
+ return ((n / 2) & 1 ? 1 : -1) * (approximate_log_of_bernoulli_bn > tools::log_max_value<T>()
+ ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, nx, Policy())
+ : static_cast<T>(exp(approximate_log_of_bernoulli_bn)));
+}
+
+template <class T, class Policy>
+T t2n_asymptotic(int n)
+{
+ BOOST_MATH_STD_USING
+ // Just get B2n and convert to a Tangent number:
+ T t2n = fabs(b2n_asymptotic<T, Policy>(2 * n)) / (2 * n);
+ T p2 = ldexp(T(1), n);
+ if(tools::max_value<T>() / p2 < t2n)
+ return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", 0, T(n), Policy());
+ t2n *= p2;
+ p2 -= 1;
+ if(tools::max_value<T>() / p2 < t2n)
+ return policies::raise_overflow_error<T>("boost::math::tangent_t2n<%1%>(std::size_t)", 0, Policy());
+ t2n *= p2;
+ return t2n;
+}
+//
+// We need to know the approximate value of /n/ which will
+// cause bernoulli_b2n<T>(n) to return infinity - this allows
+// us to elude a great deal of runtime checking for values below
+// n, and only perform the full overflow checks when we know that we're
+// getting close to the point where our calculations will overflow.
+// We use Luschny's LogB3 formula (http://www.luschny.de/math/primes/bernincl.html)
+// to find the limit, and since we're dealing with the log of the Bernoulli numbers
+// we need only perform the calculation at double precision and not with T
+// (which may be a multiprecision type). The limit returned is within 1 of the true
+// limit for all the types tested. Note that although the code below is basically
+// the same as b2n_asymptotic above, it has been recast as a continuous real-valued
+// function as this makes the root finding go smoother/faster. It also omits the
+// sign of the Bernoulli number.
+//
+struct max_bernoulli_root_functor
+{
+ max_bernoulli_root_functor(long long t) : target(static_cast<double>(t)) {}
+ double operator()(double n)
+ {
+ BOOST_MATH_STD_USING
+
+ // Luschny LogB3(n) formula.
+
+ const double nx2(n * n);
+
+ const double approximate_log_of_bernoulli_bn
+ = ((boost::math::constants::half<double>() + n) * log(n))
+ + ((boost::math::constants::half<double>() - n) * log(boost::math::constants::pi<double>()))
+ + (((double(3) / 2) - n) * boost::math::constants::ln_two<double>())
+ + ((n * (2 - (nx2 * 7) * (1 + ((nx2 * 30) * ((nx2 * 12) - 1))))) / (((nx2 * nx2) * nx2) * 2520));
+
+ return approximate_log_of_bernoulli_bn - target;
+ }
+private:
+ double target;
+};
+
+template <class T, class Policy>
+inline std::size_t find_bernoulli_overflow_limit(const mpl::false_&)
+{
+ long long t = lltrunc(boost::math::tools::log_max_value<T>());
+ max_bernoulli_root_functor fun(t);
+ boost::math::tools::equal_floor tol;
+ boost::uintmax_t max_iter = boost::math::policies::get_max_root_iterations<Policy>();
+ return static_cast<std::size_t>(boost::math::tools::toms748_solve(fun, sqrt(double(t)), double(t), tol, max_iter).first) / 2;
+}
+
+template <class T, class Policy>
+inline std::size_t find_bernoulli_overflow_limit(const mpl::true_&)
+{
+ return max_bernoulli_index<bernoulli_imp_variant<T>::value>::value;
+}
+
+template <class T, class Policy>
+std::size_t b2n_overflow_limit()
+{
+ // This routine is called at program startup if it's called at all:
+ // that guarantees safe initialization of the static variable.
+ typedef mpl::bool_<(bernoulli_imp_variant<T>::value >= 1) && (bernoulli_imp_variant<T>::value <= 3)> tag_type;
+ static const std::size_t lim = find_bernoulli_overflow_limit<T, Policy>(tag_type());
+ return lim;
+}
+
+//
+// The tangent numbers grow larger much more rapidly than the Bernoulli numbers do....
+// so to compute the Bernoulli numbers from the tangent numbers, we need to avoid spurious
+// overflow in the calculation, we can do this by scaling all the tangent number by some scale factor:
+//
+template <class T>
+inline typename enable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
+{
+ BOOST_MATH_STD_USING
+ return ldexp(T(1), std::numeric_limits<T>::min_exponent + 5);
+}
+template <class T>
+inline typename disable_if_c<std::numeric_limits<T>::is_specialized && (std::numeric_limits<T>::radix == 2), T>::type tangent_scale_factor()
+{
+ return tools::min_value<T>() * 16;
+}
+//
+// Initializer: ensure all our constants are initialized prior to the first call of main:
+//
+template <class T, class Policy>
+struct bernoulli_initializer
+{
+ struct init
+ {
+ init()
+ {
+ //
+ // We call twice, once to initialize our static table, and once to
+ // initialize our dymanic table:
+ //
+ boost::math::bernoulli_b2n<T>(2, Policy());
+ try{
+ boost::math::bernoulli_b2n<T>(max_bernoulli_b2n<T>::value + 1, Policy());
+ } catch(const std::overflow_error&){}
+ boost::math::tangent_t2n<T>(2, Policy());
+ }
+ void force_instantiate()const{}
+ };
+ static const init initializer;
+ static void force_instantiate()
+ {
+ initializer.force_instantiate();
+ }
+};
+
+template <class T, class Policy>
+const typename bernoulli_initializer<T, Policy>::init bernoulli_initializer<T, Policy>::initializer;
+
+//
+// We need something to act as a cache for our calculated Bernoulli numbers. In order to
+// ensure both fast access and thread safety, we need a stable table which may be extended
+// in size, but which never reallocates: that way values already calculated may be accessed
+// concurrently with another thread extending the table with new values.
+//
+// Very very simple vector class that will never allocate more than once, we could use
+// boost::container::static_vector here, but that allocates on the stack, which may well
+// cause issues for the amount of memory we want in the extreme case...
+//
+template <class T>
+struct fixed_vector : private std::allocator<T>
+{
+ typedef unsigned size_type;
+ typedef T* iterator;
+ typedef const T* const_iterator;
+ fixed_vector() : m_used(0)
+ {
+ std::size_t overflow_limit = 5 + b2n_overflow_limit<T, policies::policy<> >();
+ m_capacity = static_cast<unsigned>((std::min)(overflow_limit, static_cast<std::size_t>(100000u)));
+ m_data = this->allocate(m_capacity);
+ }
+ ~fixed_vector()
+ {
+ for(unsigned i = 0; i < m_used; ++i)
+ this->destroy(&m_data[i]);
+ this->deallocate(m_data, m_capacity);
+ }
+ T& operator[](unsigned n) { BOOST_ASSERT(n < m_used); return m_data[n]; }
+ const T& operator[](unsigned n)const { BOOST_ASSERT(n < m_used); return m_data[n]; }
+ unsigned size()const { return m_used; }
+ unsigned size() { return m_used; }
+ void resize(unsigned n, const T& val)
+ {
+ if(n > m_capacity)
+ throw std::runtime_error("Exhausted storage for Bernoulli numbers.");
+ for(unsigned i = m_used; i < n; ++i)
+ new (m_data + i) T(val);
+ m_used = n;
+ }
+ void resize(unsigned n) { resize(n, T()); }
+ T* begin() { return m_data; }
+ T* end() { return m_data + m_used; }
+ T* begin()const { return m_data; }
+ T* end()const { return m_data + m_used; }
+ unsigned capacity()const { return m_capacity; }
+private:
+ T* m_data;
+ unsigned m_used, m_capacity;
+};
+
+template <class T, class Policy>
+class bernoulli_numbers_cache
+{
+public:
+ bernoulli_numbers_cache() : m_overflow_limit((std::numeric_limits<std::size_t>::max)())
+#if defined(BOOST_HAS_THREADS) && !defined(BOOST_MATH_NO_ATOMIC_INT)
+ , m_counter(0)
+#endif
+ {}
+
+ typedef fixed_vector<T> container_type;
+
+ void tangent(std::size_t m)
+ {
+ static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
+ tn.resize(static_cast<typename container_type::size_type>(m), T(0U));
+
+ BOOST_MATH_INSTRUMENT_VARIABLE(min_overflow_index);
+
+ std::size_t prev_size = m_intermediates.size();
+ m_intermediates.resize(m, T(0U));
+
+ if(prev_size == 0)
+ {
+ m_intermediates[1] = tangent_scale_factor<T>() /*T(1U)*/;
+ tn[0U] = T(0U);
+ tn[1U] = tangent_scale_factor<T>()/* T(1U)*/;
+ BOOST_MATH_INSTRUMENT_VARIABLE(tn[0]);
+ BOOST_MATH_INSTRUMENT_VARIABLE(tn[1]);
+ }
+
+ for(std::size_t i = std::max<size_t>(2, prev_size); i < m; i++)
+ {
+ bool overflow_check = false;
+ if(i >= min_overflow_index && (boost::math::tools::max_value<T>() / (i-1) < m_intermediates[1]) )
+ {
+ std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
+ break;
+ }
+ m_intermediates[1] = m_intermediates[1] * (i-1);
+ for(std::size_t j = 2; j <= i; j++)
+ {
+ overflow_check =
+ (i >= min_overflow_index) && (
+ (boost::math::tools::max_value<T>() / (i - j) < m_intermediates[j])
+ || (boost::math::tools::max_value<T>() / (i - j + 2) < m_intermediates[j-1])
+ || (boost::math::tools::max_value<T>() - m_intermediates[j] * (i - j) < m_intermediates[j-1] * (i - j + 2))
+ || ((boost::math::isinf)(m_intermediates[j]))
+ );
+
+ if(overflow_check)
+ {
+ std::fill(tn.begin() + i, tn.end(), boost::math::tools::max_value<T>());
+ break;
+ }
+ m_intermediates[j] = m_intermediates[j] * (i - j) + m_intermediates[j-1] * (i - j + 2);
+ }
+ if(overflow_check)
+ break; // already filled the tn...
+ tn[static_cast<typename container_type::size_type>(i)] = m_intermediates[i];
+ BOOST_MATH_INSTRUMENT_VARIABLE(i);
+ BOOST_MATH_INSTRUMENT_VARIABLE(tn[static_cast<typename container_type::size_type>(i)]);
+ }
+ }
+
+ void tangent_numbers_series(const std::size_t m)
+ {
+ BOOST_MATH_STD_USING
+ static const std::size_t min_overflow_index = b2n_overflow_limit<T, Policy>() - 1;
+
+ typename container_type::size_type old_size = bn.size();
+
+ tangent(m);
+ bn.resize(static_cast<typename container_type::size_type>(m));
+
+ if(!old_size)
+ {
+ bn[0] = 1;
+ old_size = 1;
+ }
+
+ T power_two(ldexp(T(1), static_cast<int>(2 * old_size)));
+
+ for(std::size_t i = old_size; i < m; i++)
+ {
+ T b(static_cast<T>(i * 2));
+ //
+ // Not only do we need to take care to avoid spurious over/under flow in
+ // the calculation, but we also need to avoid overflow altogether in case
+ // we're calculating with a type where "bad things" happen in that case:
+ //
+ b = b / (power_two * tangent_scale_factor<T>());
+ b /= (power_two - 1);
+ bool overflow_check = (i >= min_overflow_index) && (tools::max_value<T>() / tn[static_cast<typename container_type::size_type>(i)] < b);
+ if(overflow_check)
+ {
+ m_overflow_limit = i;
+ while(i < m)
+ {
+ b = std::numeric_limits<T>::has_infinity ? std::numeric_limits<T>::infinity() : tools::max_value<T>();
+ bn[static_cast<typename container_type::size_type>(i)] = ((i % 2U) ? b : T(-b));
+ ++i;
+ }
+ break;
+ }
+ else
+ {
+ b *= tn[static_cast<typename container_type::size_type>(i)];
+ }
+
+ power_two = ldexp(power_two, 2);
+
+ const bool b_neg = i % 2 == 0;
+
+ bn[static_cast<typename container_type::size_type>(i)] = ((!b_neg) ? b : T(-b));
+ }
+ }
+
+ template <class OutputIterator>
+ OutputIterator copy_bernoulli_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
+ {
+ //
+ // There are basically 3 thread safety options:
+ //
+ // 1) There are no threads (BOOST_HAS_THREADS is not defined).
+ // 2) There are threads, but we do not have a true atomic integer type,
+ // in this case we just use a mutex to guard against race conditions.
+ // 3) There are threads, and we have an atomic integer: in this case we can
+ // use the double-checked locking pattern to avoid thread synchronisation
+ // when accessing values already in the cache.
+ //
+ // First off handle the common case for overflow and/or asymptotic expansion:
+ //
+ if(start + n > bn.capacity())
+ {
+ if(start < bn.capacity())
+ {
+ out = copy_bernoulli_numbers(out, start, bn.capacity() - start, pol);
+ n -= bn.capacity() - start;
+ start = static_cast<std::size_t>(bn.capacity());
+ }
+ if(start < b2n_overflow_limit<T, Policy>() + 2u)
+ {
+ for(; n; ++start, --n)
+ {
+ *out = b2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start * 2U));
+ ++out;
+ }
+ }
+ for(; n; ++start, --n)
+ {
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
+ ++out;
+ }
+ return out;
+ }
+ #if !defined(BOOST_HAS_THREADS)
+ //
+ // Single threaded code, very simple:
+ //
+ if(start + n >= bn.size())
+ {
+ std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ tangent_numbers_series(new_size);
+ }
+
+ for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
+ {
+ *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
+ ++out;
+ }
+ #elif defined(BOOST_MATH_NO_ATOMIC_INT)
+ //
+ // We need to grab a mutex every time we get here, for both readers and writers:
+ //
+ boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
+ if(start + n >= bn.size())
+ {
+ std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ tangent_numbers_series(new_size);
+ }
+
+ for(std::size_t i = (std::max)(max_bernoulli_b2n<T>::value + 1, start); i < start + n; ++i)
+ {
+ *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[i];
+ ++out;
+ }
+
+ #else
+ //
+ // Double-checked locking pattern, lets us access cached already cached values
+ // without locking:
+ //
+ // Get the counter and see if we need to calculate more constants:
+ //
+ if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
+ {
+ boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
+
+ if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
+ {
+ if(start + n >= bn.size())
+ {
+ std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ tangent_numbers_series(new_size);
+ }
+ m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
+ }
+ }
+
+ for(std::size_t i = (std::max)(static_cast<std::size_t>(max_bernoulli_b2n<T>::value + 1), start); i < start + n; ++i)
+ {
+ *out = (i >= m_overflow_limit) ? policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol) : bn[static_cast<typename container_type::size_type>(i)];
+ ++out;
+ }
+
+ #endif
+ return out;
+ }
+
+ template <class OutputIterator>
+ OutputIterator copy_tangent_numbers(OutputIterator out, std::size_t start, std::size_t n, const Policy& pol)
+ {
+ //
+ // There are basically 3 thread safety options:
+ //
+ // 1) There are no threads (BOOST_HAS_THREADS is not defined).
+ // 2) There are threads, but we do not have a true atomic integer type,
+ // in this case we just use a mutex to guard against race conditions.
+ // 3) There are threads, and we have an atomic integer: in this case we can
+ // use the double-checked locking pattern to avoid thread synchronisation
+ // when accessing values already in the cache.
+ //
+ //
+ // First off handle the common case for overflow and/or asymptotic expansion:
+ //
+ if(start + n > bn.capacity())
+ {
+ if(start < bn.capacity())
+ {
+ out = copy_tangent_numbers(out, start, bn.capacity() - start, pol);
+ n -= bn.capacity() - start;
+ start = static_cast<std::size_t>(bn.capacity());
+ }
+ if(start < b2n_overflow_limit<T, Policy>() + 2u)
+ {
+ for(; n; ++start, --n)
+ {
+ *out = t2n_asymptotic<T, Policy>(static_cast<typename container_type::size_type>(start));
+ ++out;
+ }
+ }
+ for(; n; ++start, --n)
+ {
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(start), pol);
+ ++out;
+ }
+ return out;
+ }
+ #if !defined(BOOST_HAS_THREADS)
+ //
+ // Single threaded code, very simple:
+ //
+ if(start + n >= bn.size())
+ {
+ std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ tangent_numbers_series(new_size);
+ }
+
+ for(std::size_t i = start; i < start + n; ++i)
+ {
+ if(i >= m_overflow_limit)
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
+ else
+ {
+ if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
+ else
+ *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
+ }
+ ++out;
+ }
+ #elif defined(BOOST_MATH_NO_ATOMIC_INT)
+ //
+ // We need to grab a mutex every time we get here, for both readers and writers:
+ //
+ boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
+ if(start + n >= bn.size())
+ {
+ std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ tangent_numbers_series(new_size);
+ }
+
+ for(std::size_t i = start; i < start + n; ++i)
+ {
+ if(i >= m_overflow_limit)
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
+ else
+ {
+ if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
+ else
+ *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
+ }
+ ++out;
+ }
+
+ #else
+ //
+ // Double-checked locking pattern, lets us access cached already cached values
+ // without locking:
+ //
+ // Get the counter and see if we need to calculate more constants:
+ //
+ if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
+ {
+ boost::detail::lightweight_mutex::scoped_lock l(m_mutex);
+
+ if(static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n)
+ {
+ if(start + n >= bn.size())
+ {
+ std::size_t new_size = (std::min)((std::max)((std::max)(start + n, std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity()));
+ tangent_numbers_series(new_size);
+ }
+ m_counter.store(static_cast<atomic_integer_type>(bn.size()), BOOST_MATH_ATOMIC_NS::memory_order_release);
+ }
+ }
+
+ for(std::size_t i = start; i < start + n; ++i)
+ {
+ if(i >= m_overflow_limit)
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
+ else
+ {
+ if(tools::max_value<T>() * tangent_scale_factor<T>() < tn[static_cast<typename container_type::size_type>(i)])
+ *out = policies::raise_overflow_error<T>("boost::math::bernoulli_b2n<%1%>(std::size_t)", 0, T(i), pol);
+ else
+ *out = tn[static_cast<typename container_type::size_type>(i)] / tangent_scale_factor<T>();
+ }
+ ++out;
+ }
+
+ #endif
+ return out;
+ }
+
+private:
+ //
+ // The caches for Bernoulli and tangent numbers, once allocated,
+ // these must NEVER EVER reallocate as it breaks our thread
+ // safety guarentees:
+ //
+ fixed_vector<T> bn, tn;
+ std::vector<T> m_intermediates;
+ // The value at which we know overflow has already occured for the Bn:
+ std::size_t m_overflow_limit;
+#if !defined(BOOST_HAS_THREADS)
+#elif defined(BOOST_MATH_NO_ATOMIC_INT)
+ boost::detail::lightweight_mutex m_mutex;
+#else
+ boost::detail::lightweight_mutex m_mutex;
+ atomic_counter_type m_counter;
+#endif
+};
+
+template <class T, class Policy>
+inline bernoulli_numbers_cache<T, Policy>& get_bernoulli_numbers_cache()
+{
+ //
+ // Force this function to be called at program startup so all the static variables
+ // get initailzed then (thread safety).
+ //
+ bernoulli_initializer<T, Policy>::force_instantiate();
+ static bernoulli_numbers_cache<T, Policy> data;
+ return data;
+}
+
+}}}
+
+#endif // BOOST_MATH_BERNOULLI_DETAIL_HPP
diff --git a/boost/math/special_functions/detail/bessel_derivatives_linear.hpp b/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
new file mode 100644
index 0000000..2ee86a0
--- /dev/null
+++ b/boost/math/special_functions/detail/bessel_derivatives_linear.hpp
@@ -0,0 +1,75 @@
+// Copyright (c) 2013 Anton Bikineev
+// 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)
+
+//
+// This is a partial header, do not include on it's own!!!
+//
+// Linear combination for bessel derivatives are defined here
+#ifndef BOOST_MATH_SF_DETAIL_BESSEL_DERIVATIVES_LINEAR_HPP
+#define BOOST_MATH_SF_DETAIL_BESSEL_DERIVATIVES_LINEAR_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class T, class Tag, class Policy>
+inline T bessel_j_derivative_linear(T v, T x, Tag tag, Policy pol)
+{
+ return (boost::math::detail::cyl_bessel_j_imp<T>(v-1, x, tag, pol) - boost::math::detail::cyl_bessel_j_imp<T>(v+1, x, tag, pol)) / 2;
+}
+
+template <class T, class Policy>
+inline T bessel_j_derivative_linear(T v, T x, const bessel_int_tag& tag, Policy pol)
+{
+ return (boost::math::detail::cyl_bessel_j_imp<T>(itrunc(v-1), x, tag, pol) - boost::math::detail::cyl_bessel_j_imp<T>(itrunc(v+1), x, tag, pol)) / 2;
+}
+
+template <class T, class Policy>
+inline T sph_bessel_j_derivative_linear(unsigned v, T x, Policy pol)
+{
+ return (v / x) * boost::math::detail::sph_bessel_j_imp<T>(v, x, pol) - boost::math::detail::sph_bessel_j_imp<T>(v+1, x, pol);
+}
+
+template <class T, class Policy>
+inline T bessel_i_derivative_linear(T v, T x, Policy pol)
+{
+ return (boost::math::detail::cyl_bessel_i_imp<T>(v-1, x, pol) + boost::math::detail::cyl_bessel_i_imp<T>(v+1, x, pol)) / 2;
+}
+
+template <class T, class Tag, class Policy>
+inline T bessel_k_derivative_linear(T v, T x, Tag tag, Policy pol)
+{
+ return (boost::math::detail::cyl_bessel_k_imp<T>(v-1, x, tag, pol) + boost::math::detail::cyl_bessel_k_imp<T>(v+1, x, tag, pol)) / -2;
+}
+
+template <class T, class Policy>
+inline T bessel_k_derivative_linear(T v, T x, const bessel_int_tag& tag, Policy pol)
+{
+ return (boost::math::detail::cyl_bessel_k_imp<T>(itrunc(v-1), x, tag, pol) + boost::math::detail::cyl_bessel_k_imp<T>(itrunc(v+1), x, tag, pol)) / -2;
+}
+
+template <class T, class Tag, class Policy>
+inline T bessel_y_derivative_linear(T v, T x, Tag tag, Policy pol)
+{
+ return (boost::math::detail::cyl_neumann_imp<T>(v-1, x, tag, pol) - boost::math::detail::cyl_neumann_imp<T>(v+1, x, tag, pol)) / 2;
+}
+
+template <class T, class Policy>
+inline T bessel_y_derivative_linear(T v, T x, const bessel_int_tag& tag, Policy pol)
+{
+ return (boost::math::detail::cyl_neumann_imp<T>(itrunc(v-1), x, tag, pol) - boost::math::detail::cyl_neumann_imp<T>(itrunc(v+1), x, tag, pol)) / 2;
+}
+
+template <class T, class Policy>
+inline T sph_neumann_derivative_linear(unsigned v, T x, Policy pol)
+{
+ return (v / x) * boost::math::detail::sph_neumann_imp<T>(v, x, pol) - boost::math::detail::sph_neumann_imp<T>(v+1, x, pol);
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_SF_DETAIL_BESSEL_DERIVATIVES_LINEAR_HPP
diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp
index 7dc65d1..676eb71 100644
--- a/boost/math/special_functions/detail/bessel_i0.hpp
+++ b/boost/math/special_functions/detail/bessel_i0.hpp
@@ -102,10 +102,7 @@ T bessel_i0(T x)
BOOST_MATH_STD_USING
using namespace boost::math::tools;
- if (x < 0)
- {
- x = -x; // even function
- }
+ BOOST_ASSERT(x >= 0); // negative x is handled before we get here
if (x == 0)
{
return static_cast<T>(1);
diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp
index 47f1b79..b85bc67 100644
--- a/boost/math/special_functions/detail/bessel_i1.hpp
+++ b/boost/math/special_functions/detail/bessel_i1.hpp
@@ -103,6 +103,7 @@ T bessel_i1(T x)
BOOST_MATH_STD_USING
using namespace boost::math::tools;
+ BOOST_ASSERT(x >= 0); // negative x is handled before we get here
w = abs(x);
if (x == 0)
{
@@ -123,10 +124,6 @@ T bessel_i1(T x)
value = factor * r;
}
- if (x < 0)
- {
- value *= -value; // odd function
- }
return value;
}
diff --git a/boost/math/special_functions/detail/bessel_ik.hpp b/boost/math/special_functions/detail/bessel_ik.hpp
index a589673..10118d9 100644
--- a/boost/math/special_functions/detail/bessel_ik.hpp
+++ b/boost/math/special_functions/detail/bessel_ik.hpp
@@ -234,6 +234,7 @@ int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)
BOOST_MATH_INSTRUMENT_VARIABLE(b);
BOOST_MATH_INSTRUMENT_VARIABLE(D);
BOOST_MATH_INSTRUMENT_VARIABLE(f);
+
for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) // starting from 2
{
// continued fraction f = z1 / z0
@@ -250,10 +251,27 @@ int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)
C *= -a / k;
Q += C * q;
S += Q * delta;
+ //
+ // Under some circumstances q can grow very small and C very
+ // large, leading to under/overflow. This is particularly an
+ // issue for types which have many digits precision but a narrow
+ // exponent range. A typical example being a "double double" type.
+ // To avoid this situation we can normalise q (and related prev/current)
+ // and C. All other variables remain unchanged in value. A typical
+ // test case occurs when x is close to 2, for example cyl_bessel_k(9.125, 2.125).
+ //
+ if(q < tools::epsilon<T>())
+ {
+ C *= q;
+ prev /= q;
+ current /= q;
+ q = 1;
+ }
// S converges slower than f
BOOST_MATH_INSTRUMENT_VARIABLE(Q * delta);
BOOST_MATH_INSTRUMENT_VARIABLE(abs(S) * tolerance);
+ BOOST_MATH_INSTRUMENT_VARIABLE(S);
if (abs(Q * delta) < abs(S) * tolerance)
{
break;
@@ -261,7 +279,10 @@ int CF2_ik(T v, T x, T* Kv, T* Kv1, const Policy& pol)
}
policies::check_series_iterations<T>("boost::math::bessel_ik<%1%>(%1%,%1%) in CF2_ik", k, pol);
- *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
+ if(x >= tools::log_max_value<T>())
+ *Kv = exp(0.5f * log(pi<T>() / (2 * x)) - x - log(S));
+ else
+ *Kv = sqrt(pi<T>() / (2 * x)) * exp(-x) / S;
*Kv1 = *Kv * (0.5f + v + x + (v * v - 0.25f) * f) / x;
BOOST_MATH_INSTRUMENT_VARIABLE(*Kv);
BOOST_MATH_INSTRUMENT_VARIABLE(*Kv1);
diff --git a/boost/math/special_functions/detail/bessel_j0.hpp b/boost/math/special_functions/detail/bessel_j0.hpp
index a07052d..ebcab17 100644
--- a/boost/math/special_functions/detail/bessel_j0.hpp
+++ b/boost/math/special_functions/detail/bessel_j0.hpp
@@ -165,13 +165,23 @@ T bessel_j0(T x)
{
T y = 8 / x;
T y2 = y * y;
- T z = x - 0.25f * pi<T>();
BOOST_ASSERT(sizeof(PC) == sizeof(QC));
BOOST_ASSERT(sizeof(PS) == sizeof(QS));
rc = evaluate_rational(PC, QC, y2);
rs = evaluate_rational(PS, QS, y2);
- factor = sqrt(2 / (x * pi<T>()));
- value = factor * (rc * cos(z) - y * rs * sin(z));
+ factor = constants::one_div_root_pi<T>() / sqrt(x);
+ //
+ // What follows is really just:
+ //
+ // T z = x - pi/4;
+ // value = factor * (rc * cos(z) - y * rs * sin(z));
+ //
+ // But using the addition formulae for sin and cos, plus
+ // the special values for sin/cos of pi/4.
+ //
+ T sx = sin(x);
+ T cx = cos(x);
+ value = factor * (rc * (cx + sx) - y * rs * (sx - cx));
}
return value;
diff --git a/boost/math/special_functions/detail/bessel_j1.hpp b/boost/math/special_functions/detail/bessel_j1.hpp
index 09d862c..91ecd28 100644
--- a/boost/math/special_functions/detail/bessel_j1.hpp
+++ b/boost/math/special_functions/detail/bessel_j1.hpp
@@ -166,13 +166,24 @@ T bessel_j1(T x)
{
T y = 8 / w;
T y2 = y * y;
- T z = w - 0.75f * pi<T>();
BOOST_ASSERT(sizeof(PC) == sizeof(QC));
BOOST_ASSERT(sizeof(PS) == sizeof(QS));
rc = evaluate_rational(PC, QC, y2);
rs = evaluate_rational(PS, QS, y2);
- factor = sqrt(2 / (w * pi<T>()));
- value = factor * (rc * cos(z) - y * rs * sin(z));
+ factor = 1 / (sqrt(w) * constants::root_pi<T>());
+ //
+ // What follows is really just:
+ //
+ // T z = w - 0.75f * pi<T>();
+ // value = factor * (rc * cos(z) - y * rs * sin(z));
+ //
+ // but using the sin/cos addition rules plus constants
+ // for the values of sin/cos of 3PI/4 which then cancel
+ // out with corresponding terms in "factor".
+ //
+ T sx = sin(x);
+ T cx = cos(x);
+ value = factor * (rc * (sx - cx) + y * rs * (sx + cx));
}
if (x < 0)
diff --git a/boost/math/special_functions/detail/bessel_jn.hpp b/boost/math/special_functions/detail/bessel_jn.hpp
index 2bf8d78..3f15f9c 100644
--- a/boost/math/special_functions/detail/bessel_jn.hpp
+++ b/boost/math/special_functions/detail/bessel_jn.hpp
@@ -42,6 +42,11 @@ T bessel_jn(int n, T x, const Policy& pol)
{
factor = 1;
}
+ if(x < 0)
+ {
+ factor *= (n & 0x1) ? -1 : 1; // J_{n}(-z) = (-1)^n J_n(z)
+ x = -x;
+ }
//
// Special cases:
//
@@ -59,8 +64,7 @@ T bessel_jn(int n, T x, const Policy& pol)
return static_cast<T>(0);
}
- typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
- if(fabs(x) > asymptotic_bessel_j_limit<T>(n, tag_type()))
+ if(asymptotic_bessel_large_x_limit(T(n), x))
return factor * asymptotic_bessel_j_large_x_2<T>(n, x);
BOOST_ASSERT(n > 1);
@@ -69,6 +73,7 @@ T bessel_jn(int n, T x, const Policy& pol)
{
prev = bessel_j0(x);
current = bessel_j1(x);
+ policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
for (int k = 1; k < n; k++)
{
T fact = 2 * k / x;
@@ -86,7 +91,7 @@ T bessel_jn(int n, T x, const Policy& pol)
current = value;
}
}
- else if(x < 1)
+ else if((x < 1) || (n > x * x / 4) || (x < 5))
{
return factor * bessel_j_small_z_series(T(n), x, pol);
}
@@ -97,6 +102,8 @@ T bessel_jn(int n, T x, const Policy& pol)
boost::math::detail::CF1_jy(static_cast<T>(n), x, &fn, &s, pol);
prev = fn;
current = 1;
+ // Check recursion won't go on too far:
+ policies::check_series_iterations<T>("boost::math::bessel_j_n<%1%>(%1%,%1%)", n, pol);
for (int k = n; k > 0; k--)
{
T fact = 2 * k / x;
diff --git a/boost/math/special_functions/detail/bessel_jy.hpp b/boost/math/special_functions/detail/bessel_jy.hpp
index d60dda2..b67d989 100644
--- a/boost/math/special_functions/detail/bessel_jy.hpp
+++ b/boost/math/special_functions/detail/bessel_jy.hpp
@@ -28,440 +28,464 @@
namespace boost { namespace math {
-namespace detail {
-
-//
-// Simultaneous calculation of A&S 9.2.9 and 9.2.10
-// for use in A&S 9.2.5 and 9.2.6.
-// This series is quick to evaluate, but divergent unless
-// x is very large, in fact it's pretty hard to figure out
-// with any degree of precision when this series actually
-// *will* converge!! Consequently, we may just have to
-// try it and see...
-//
-template <class T, class Policy>
-bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )
-{
- BOOST_MATH_STD_USING
- T tolerance = 2 * policies::get_epsilon<T, Policy>();
- *p = 1;
- *q = 0;
- T k = 1;
- T z8 = 8 * x;
- T sq = 1;
- T mu = 4 * v * v;
- T term = 1;
- bool ok = true;
- do
- {
- term *= (mu - sq * sq) / (k * z8);
- *q += term;
- k += 1;
- sq += 2;
- T mult = (sq * sq - mu) / (k * z8);
- ok = fabs(mult) < 0.5f;
- term *= mult;
- *p += term;
- k += 1;
- sq += 2;
- }
- while((fabs(term) > tolerance * *p) && ok);
- return ok;
-}
-
-// Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
-// Temme, Journal of Computational Physics, vol 21, 343 (1976)
-template <typename T, typename Policy>
-int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
-{
- T g, h, p, q, f, coef, sum, sum1, tolerance;
- T a, d, e, sigma;
- unsigned long k;
-
- BOOST_MATH_STD_USING
- using namespace boost::math::tools;
- using namespace boost::math::constants;
-
- BOOST_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine
-
- T gp = boost::math::tgamma1pm1(v, pol);
- T gm = boost::math::tgamma1pm1(-v, pol);
- T spv = boost::math::sin_pi(v, pol);
- T spv2 = boost::math::sin_pi(v/2, pol);
- T xp = pow(x/2, v);
-
- a = log(x / 2);
- sigma = -a * v;
- d = abs(sigma) < tools::epsilon<T>() ?
- T(1) : sinh(sigma) / sigma;
- e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
- : T(2 * spv2 * spv2 / v);
-
- T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
- T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
- T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
- f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
-
- p = vspv / (xp * (1 + gm));
- q = vspv * xp / (1 + gp);
-
- g = f + e * q;
- h = p;
- coef = 1;
- sum = coef * g;
- sum1 = coef * h;
-
- T v2 = v * v;
- T coef_mult = -x * x / 4;
-
- // series summation
- tolerance = policies::get_epsilon<T, Policy>();
- for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
- {
- f = (k * f + p + q) / (k*k - v2);
- p /= k - v;
- q /= k + v;
- g = f + e * q;
- h = p - k * g;
- coef *= coef_mult / k;
- sum += coef * g;
- sum1 += coef * h;
- if (abs(coef * g) < abs(sum) * tolerance)
- {
- break;
- }
- }
- policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
- *Y = -sum;
- *Y1 = -2 * sum1 / x;
-
- return 0;
-}
-
-// Evaluate continued fraction fv = J_(v+1) / J_v, see
-// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
-template <typename T, typename Policy>
-int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
-{
- T C, D, f, a, b, delta, tiny, tolerance;
- unsigned long k;
- int s = 1;
-
- BOOST_MATH_STD_USING
-
- // |x| <= |v|, CF1_jy converges rapidly
- // |x| > |v|, CF1_jy needs O(|x|) iterations to converge
-
- // modified Lentz's method, see
- // Lentz, Applied Optics, vol 15, 668 (1976)
- tolerance = 2 * policies::get_epsilon<T, Policy>();;
- tiny = sqrt(tools::min_value<T>());
- C = f = tiny; // b0 = 0, replace with tiny
- D = 0;
- for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
- {
- a = -1;
- b = 2 * (v + k) / x;
- C = b + a / C;
- D = b + a * D;
- if (C == 0) { C = tiny; }
- if (D == 0) { D = tiny; }
- D = 1 / D;
- delta = C * D;
- f *= delta;
- if (D < 0) { s = -s; }
- if (abs(delta - 1) < tolerance)
- { break; }
- }
- policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
- *fv = -f;
- *sign = s; // sign of denominator
-
- return 0;
-}
-//
-// This algorithm was originally written by Xiaogang Zhang
-// using std::complex to perform the complex arithmetic.
-// However, that turns out to 10x or more slower than using
-// all real-valued arithmetic, so it's been rewritten using
-// real values only.
-//
-template <typename T, typename Policy>
-int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
-{
- BOOST_MATH_STD_USING
-
- T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
- T tiny;
- unsigned long k;
-
- // |x| >= |v|, CF2_jy converges rapidly
- // |x| -> 0, CF2_jy fails to converge
- BOOST_ASSERT(fabs(x) > 1);
-
- // modified Lentz's method, complex numbers involved, see
- // Lentz, Applied Optics, vol 15, 668 (1976)
- T tolerance = 2 * policies::get_epsilon<T, Policy>();
- tiny = sqrt(tools::min_value<T>());
- Cr = fr = -0.5f / x;
- Ci = fi = 1;
- //Dr = Di = 0;
- T v2 = v * v;
- a = (0.25f - v2) / x; // Note complex this one time only!
- br = 2 * x;
- bi = 2;
- temp = Cr * Cr + 1;
- Ci = bi + a * Cr / temp;
- Cr = br + a / temp;
- Dr = br;
- Di = bi;
- if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
- if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
- temp = Dr * Dr + Di * Di;
- Dr = Dr / temp;
- Di = -Di / temp;
- delta_r = Cr * Dr - Ci * Di;
- delta_i = Ci * Dr + Cr * Di;
- temp = fr;
- fr = temp * delta_r - fi * delta_i;
- fi = temp * delta_i + fi * delta_r;
- for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
- {
- a = k - 0.5f;
- a *= a;
- a -= v2;
- bi += 2;
- temp = Cr * Cr + Ci * Ci;
- Cr = br + a * Cr / temp;
- Ci = bi - a * Ci / temp;
- Dr = br + a * Dr;
- Di = bi + a * Di;
- if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
- if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
- temp = Dr * Dr + Di * Di;
- Dr = Dr / temp;
- Di = -Di / temp;
- delta_r = Cr * Dr - Ci * Di;
- delta_i = Ci * Dr + Cr * Di;
- temp = fr;
- fr = temp * delta_r - fi * delta_i;
- fi = temp * delta_i + fi * delta_r;
- if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
- break;
- }
- policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
- *p = fr;
- *q = fi;
-
- return 0;
-}
-
-enum
-{
- need_j = 1, need_y = 2
-};
-
-// Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
-// Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
-template <typename T, typename Policy>
-int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
-{
- BOOST_ASSERT(x >= 0);
-
- T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
- T W, p, q, gamma, current, prev, next;
- bool reflect = false;
- unsigned n, k;
- int s;
- int org_kind = kind;
- T cp = 0;
- T sp = 0;
-
- static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
-
- BOOST_MATH_STD_USING
- using namespace boost::math::tools;
- using namespace boost::math::constants;
-
- if (v < 0)
- {
- reflect = true;
- v = -v; // v is non-negative from here
- kind = need_j|need_y; // need both for reflection formula
- }
- n = iround(v, pol);
- u = v - n; // -1/2 <= u < 1/2
-
- if(reflect)
- {
- T z = (u + n % 2);
- cp = boost::math::cos_pi(z, pol);
- sp = boost::math::sin_pi(z, pol);
- }
-
- if (x == 0)
- {
- *J = *Y = policies::raise_overflow_error<T>(
- function, 0, pol);
- return 1;
- }
-
- // x is positive until reflection
- W = T(2) / (x * pi<T>()); // Wronskian
- T Yv_scale = 1;
- if((x > 8) && (x < 1000) && hankel_PQ(v, x, &p, &q, pol))
- {
- //
- // Hankel approximation: note that this method works best when x
- // is large, but in that case we end up calculating sines and cosines
- // of large values, with horrendous resulting accuracy. It is fast though
- // when it works....
- //
- T chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
- T sc = sin(chi);
- T cc = cos(chi);
- chi = sqrt(2 / (boost::math::constants::pi<T>() * x));
- Yv = chi * (p * sc + q * cc);
- Jv = chi * (p * cc - q * sc);
- }
- else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
- {
- // Evaluate using series representations.
- // This is particularly important for x << v as in this
- // area temme_jy may be slow to converge, if it converges at all.
- // Requires x is not an integer.
- if(kind&need_j)
- Jv = bessel_j_small_z_series(v, x, pol);
- else
- Jv = std::numeric_limits<T>::quiet_NaN();
- if((org_kind&need_y && (!reflect || (cp != 0)))
- || (org_kind & need_j && (reflect && (sp != 0))))
- {
- // Only calculate if we need it, and if the reflection formula will actually use it:
- Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN();
- }
- else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
- {
- // Truncated series evaluation for small x and v an integer,
- // much quicker in this area than temme_jy below.
- if(kind&need_j)
- Jv = bessel_j_small_z_series(v, x, pol);
- else
- Jv = std::numeric_limits<T>::quiet_NaN();
- if((org_kind&need_y && (!reflect || (cp != 0)))
- || (org_kind & need_j && (reflect && (sp != 0))))
- {
- // Only calculate if we need it, and if the reflection formula will actually use it:
- Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN();
- }
- else if (x <= 2) // x in (0, 2]
- {
- if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series
- {
- // domain error:
- *J = *Y = Yu;
- return 1;
- }
- prev = Yu;
- current = Yu1;
- T scale = 1;
- for (k = 1; k <= n; k++) // forward recurrence for Y
- {
- T fact = 2 * (u + k) / x;
- if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+ namespace detail {
+
+ //
+ // Simultaneous calculation of A&S 9.2.9 and 9.2.10
+ // for use in A&S 9.2.5 and 9.2.6.
+ // This series is quick to evaluate, but divergent unless
+ // x is very large, in fact it's pretty hard to figure out
+ // with any degree of precision when this series actually
+ // *will* converge!! Consequently, we may just have to
+ // try it and see...
+ //
+ template <class T, class Policy>
+ bool hankel_PQ(T v, T x, T* p, T* q, const Policy& )
+ {
+ BOOST_MATH_STD_USING
+ T tolerance = 2 * policies::get_epsilon<T, Policy>();
+ *p = 1;
+ *q = 0;
+ T k = 1;
+ T z8 = 8 * x;
+ T sq = 1;
+ T mu = 4 * v * v;
+ T term = 1;
+ bool ok = true;
+ do
+ {
+ term *= (mu - sq * sq) / (k * z8);
+ *q += term;
+ k += 1;
+ sq += 2;
+ T mult = (sq * sq - mu) / (k * z8);
+ ok = fabs(mult) < 0.5f;
+ term *= mult;
+ *p += term;
+ k += 1;
+ sq += 2;
+ }
+ while((fabs(term) > tolerance * *p) && ok);
+ return ok;
+ }
+
+ // Calculate Y(v, x) and Y(v+1, x) by Temme's method, see
+ // Temme, Journal of Computational Physics, vol 21, 343 (1976)
+ template <typename T, typename Policy>
+ int temme_jy(T v, T x, T* Y, T* Y1, const Policy& pol)
+ {
+ T g, h, p, q, f, coef, sum, sum1, tolerance;
+ T a, d, e, sigma;
+ unsigned long k;
+
+ BOOST_MATH_STD_USING
+ using namespace boost::math::tools;
+ using namespace boost::math::constants;
+
+ BOOST_ASSERT(fabs(v) <= 0.5f); // precondition for using this routine
+
+ T gp = boost::math::tgamma1pm1(v, pol);
+ T gm = boost::math::tgamma1pm1(-v, pol);
+ T spv = boost::math::sin_pi(v, pol);
+ T spv2 = boost::math::sin_pi(v/2, pol);
+ T xp = pow(x/2, v);
+
+ a = log(x / 2);
+ sigma = -a * v;
+ d = abs(sigma) < tools::epsilon<T>() ?
+ T(1) : sinh(sigma) / sigma;
+ e = abs(v) < tools::epsilon<T>() ? T(v*pi<T>()*pi<T>() / 2)
+ : T(2 * spv2 * spv2 / v);
+
+ T g1 = (v == 0) ? T(-euler<T>()) : T((gp - gm) / ((1 + gp) * (1 + gm) * 2 * v));
+ T g2 = (2 + gp + gm) / ((1 + gp) * (1 + gm) * 2);
+ T vspv = (fabs(v) < tools::epsilon<T>()) ? T(1/constants::pi<T>()) : T(v / spv);
+ f = (g1 * cosh(sigma) - g2 * a * d) * 2 * vspv;
+
+ p = vspv / (xp * (1 + gm));
+ q = vspv * xp / (1 + gp);
+
+ g = f + e * q;
+ h = p;
+ coef = 1;
+ sum = coef * g;
+ sum1 = coef * h;
+
+ T v2 = v * v;
+ T coef_mult = -x * x / 4;
+
+ // series summation
+ tolerance = policies::get_epsilon<T, Policy>();
+ for (k = 1; k < policies::get_max_series_iterations<Policy>(); k++)
+ {
+ f = (k * f + p + q) / (k*k - v2);
+ p /= k - v;
+ q /= k + v;
+ g = f + e * q;
+ h = p - k * g;
+ coef *= coef_mult / k;
+ sum += coef * g;
+ sum1 += coef * h;
+ if (abs(coef * g) < abs(sum) * tolerance)
+ {
+ break;
+ }
+ }
+ policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in temme_jy", k, pol);
+ *Y = -sum;
+ *Y1 = -2 * sum1 / x;
+
+ return 0;
+ }
+
+ // Evaluate continued fraction fv = J_(v+1) / J_v, see
+ // Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
+ template <typename T, typename Policy>
+ int CF1_jy(T v, T x, T* fv, int* sign, const Policy& pol)
+ {
+ T C, D, f, a, b, delta, tiny, tolerance;
+ unsigned long k;
+ int s = 1;
+
+ BOOST_MATH_STD_USING
+
+ // |x| <= |v|, CF1_jy converges rapidly
+ // |x| > |v|, CF1_jy needs O(|x|) iterations to converge
+
+ // modified Lentz's method, see
+ // Lentz, Applied Optics, vol 15, 668 (1976)
+ tolerance = 2 * policies::get_epsilon<T, Policy>();;
+ tiny = sqrt(tools::min_value<T>());
+ C = f = tiny; // b0 = 0, replace with tiny
+ D = 0;
+ for (k = 1; k < policies::get_max_series_iterations<Policy>() * 100; k++)
+ {
+ a = -1;
+ b = 2 * (v + k) / x;
+ C = b + a / C;
+ D = b + a * D;
+ if (C == 0) { C = tiny; }
+ if (D == 0) { D = tiny; }
+ D = 1 / D;
+ delta = C * D;
+ f *= delta;
+ if (D < 0) { s = -s; }
+ if (abs(delta - 1) < tolerance)
+ { break; }
+ }
+ policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF1_jy", k / 100, pol);
+ *fv = -f;
+ *sign = s; // sign of denominator
+
+ return 0;
+ }
+ //
+ // This algorithm was originally written by Xiaogang Zhang
+ // using std::complex to perform the complex arithmetic.
+ // However, that turns out to 10x or more slower than using
+ // all real-valued arithmetic, so it's been rewritten using
+ // real values only.
+ //
+ template <typename T, typename Policy>
+ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol)
+ {
+ BOOST_MATH_STD_USING
+
+ T Cr, Ci, Dr, Di, fr, fi, a, br, bi, delta_r, delta_i, temp;
+ T tiny;
+ unsigned long k;
+
+ // |x| >= |v|, CF2_jy converges rapidly
+ // |x| -> 0, CF2_jy fails to converge
+ BOOST_ASSERT(fabs(x) > 1);
+
+ // modified Lentz's method, complex numbers involved, see
+ // Lentz, Applied Optics, vol 15, 668 (1976)
+ T tolerance = 2 * policies::get_epsilon<T, Policy>();
+ tiny = sqrt(tools::min_value<T>());
+ Cr = fr = -0.5f / x;
+ Ci = fi = 1;
+ //Dr = Di = 0;
+ T v2 = v * v;
+ a = (0.25f - v2) / x; // Note complex this one time only!
+ br = 2 * x;
+ bi = 2;
+ temp = Cr * Cr + 1;
+ Ci = bi + a * Cr / temp;
+ Cr = br + a / temp;
+ Dr = br;
+ Di = bi;
+ if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
+ if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
+ temp = Dr * Dr + Di * Di;
+ Dr = Dr / temp;
+ Di = -Di / temp;
+ delta_r = Cr * Dr - Ci * Di;
+ delta_i = Ci * Dr + Cr * Di;
+ temp = fr;
+ fr = temp * delta_r - fi * delta_i;
+ fi = temp * delta_i + fi * delta_r;
+ for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++)
+ {
+ a = k - 0.5f;
+ a *= a;
+ a -= v2;
+ bi += 2;
+ temp = Cr * Cr + Ci * Ci;
+ Cr = br + a * Cr / temp;
+ Ci = bi - a * Ci / temp;
+ Dr = br + a * Dr;
+ Di = bi + a * Di;
+ if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; }
+ if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; }
+ temp = Dr * Dr + Di * Di;
+ Dr = Dr / temp;
+ Di = -Di / temp;
+ delta_r = Cr * Dr - Ci * Di;
+ delta_i = Ci * Dr + Cr * Di;
+ temp = fr;
+ fr = temp * delta_r - fi * delta_i;
+ fi = temp * delta_i + fi * delta_r;
+ if (fabs(delta_r - 1) + fabs(delta_i) < tolerance)
+ break;
+ }
+ policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol);
+ *p = fr;
+ *q = fi;
+
+ return 0;
+ }
+
+ static const int need_j = 1;
+ static const int need_y = 2;
+
+ // Compute J(v, x) and Y(v, x) simultaneously by Steed's method, see
+ // Barnett et al, Computer Physics Communications, vol 8, 377 (1974)
+ template <typename T, typename Policy>
+ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
+ {
+ BOOST_ASSERT(x >= 0);
+
+ T u, Jv, Ju, Yv, Yv1, Yu, Yu1(0), fv, fu;
+ T W, p, q, gamma, current, prev, next;
+ bool reflect = false;
+ unsigned n, k;
+ int s;
+ int org_kind = kind;
+ T cp = 0;
+ T sp = 0;
+
+ static const char* function = "boost::math::bessel_jy<%1%>(%1%,%1%)";
+
+ BOOST_MATH_STD_USING
+ using namespace boost::math::tools;
+ using namespace boost::math::constants;
+
+ if (v < 0)
+ {
+ reflect = true;
+ v = -v; // v is non-negative from here
+ }
+ if (v > static_cast<T>((std::numeric_limits<int>::max)()))
+ {
+ *J = *Y = policies::raise_evaluation_error<T>(function, "Order of Bessel function is too large to evaluate: got %1%", v, pol);
+ return 1;
+ }
+ n = iround(v, pol);
+ u = v - n; // -1/2 <= u < 1/2
+
+ if(reflect)
+ {
+ T z = (u + n % 2);
+ cp = boost::math::cos_pi(z, pol);
+ sp = boost::math::sin_pi(z, pol);
+ if(u != 0)
+ kind = need_j|need_y; // need both for reflection formula
+ }
+
+ if(x == 0)
+ {
+ if(v == 0)
+ *J = 1;
+ else if((u == 0) || !reflect)
+ *J = 0;
+ else if(kind & need_j)
+ *J = policies::raise_domain_error<T>(function, "Value of Bessel J_v(x) is complex-infinity at %1%", x, pol); // complex infinity
+ else
+ *J = std::numeric_limits<T>::quiet_NaN(); // any value will do, not using J.
+
+ if((kind & need_y) == 0)
+ *Y = std::numeric_limits<T>::quiet_NaN(); // any value will do, not using Y.
+ else if(v == 0)
+ *Y = -policies::raise_overflow_error<T>(function, 0, pol);
+ else
+ *Y = policies::raise_domain_error<T>(function, "Value of Bessel Y_v(x) is complex-infinity at %1%", x, pol); // complex infinity
+ return 1;
+ }
+
+ // x is positive until reflection
+ W = T(2) / (x * pi<T>()); // Wronskian
+ T Yv_scale = 1;
+ if(((kind & need_y) == 0) && ((x < 1) || (v > x * x / 4) || (x < 5)))
+ {
+ //
+ // This series will actually converge rapidly for all small
+ // x - say up to x < 20 - but the first few terms are large
+ // and divergent which leads to large errors :-(
+ //
+ Jv = bessel_j_small_z_series(v, x, pol);
+ Yv = std::numeric_limits<T>::quiet_NaN();
+ }
+ else if((x < 1) && (u != 0) && (log(policies::get_epsilon<T, Policy>() / 2) > v * log((x/2) * (x/2) / v)))
+ {
+ // Evaluate using series representations.
+ // This is particularly important for x << v as in this
+ // area temme_jy may be slow to converge, if it converges at all.
+ // Requires x is not an integer.
+ if(kind&need_j)
+ Jv = bessel_j_small_z_series(v, x, pol);
+ else
+ Jv = std::numeric_limits<T>::quiet_NaN();
+ if((org_kind&need_y && (!reflect || (cp != 0)))
+ || (org_kind & need_j && (reflect && (sp != 0))))
{
- scale /= current;
- prev /= current;
- current = 1;
+ // Only calculate if we need it, and if the reflection formula will actually use it:
+ Yv = bessel_y_small_z_series(v, x, &Yv_scale, pol);
}
- next = fact * current - prev;
- prev = current;
- current = next;
- }
- Yv = prev;
- Yv1 = current;
- if(kind&need_j)
+ else
+ Yv = std::numeric_limits<T>::quiet_NaN();
+ }
+ else if((u == 0) && (x < policies::get_epsilon<T, Policy>()))
{
- CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy
- Jv = scale * W / (Yv * fv - Yv1); // Wronskian relation
+ // Truncated series evaluation for small x and v an integer,
+ // much quicker in this area than temme_jy below.
+ if(kind&need_j)
+ Jv = bessel_j_small_z_series(v, x, pol);
+ else
+ Jv = std::numeric_limits<T>::quiet_NaN();
+ if((org_kind&need_y && (!reflect || (cp != 0)))
+ || (org_kind & need_j && (reflect && (sp != 0))))
+ {
+ // Only calculate if we need it, and if the reflection formula will actually use it:
+ Yv = bessel_yn_small_z(n, x, &Yv_scale, pol);
+ }
+ else
+ Yv = std::numeric_limits<T>::quiet_NaN();
}
- else
- Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
- Yv_scale = scale;
- }
- else // x in (2, \infty)
- {
- // Get Y(u, x):
- // define tag type that will dispatch to right limits:
- typedef typename bessel_asymptotic_tag<T, Policy>::type tag_type;
-
- T lim, ratio;
- switch(kind)
- {
- case need_j:
- lim = asymptotic_bessel_j_limit<T>(v, tag_type());
- break;
- case need_y:
- lim = asymptotic_bessel_y_limit<T>(tag_type());
- break;
- default:
- lim = (std::max)(
- asymptotic_bessel_j_limit<T>(v, tag_type()),
- asymptotic_bessel_y_limit<T>(tag_type()));
- break;
- }
- if(x > lim)
- {
- if(kind&need_y)
- {
- Yu = asymptotic_bessel_y_large_x_2(u, x);
- Yu1 = asymptotic_bessel_y_large_x_2(T(u + 1), x);
- }
- else
- Yu = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
- if(kind&need_j)
- {
- Jv = asymptotic_bessel_j_large_x_2(v, x);
- }
- else
- Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
- }
- else
- {
- CF1_jy(v, x, &fv, &s, pol);
- // tiny initial value to prevent overflow
- T init = sqrt(tools::min_value<T>());
- prev = fv * s * init;
- current = s * init;
- if(v < max_factorial<T>::value)
- {
- for (k = n; k > 0; k--) // backward recurrence for J
- {
+ else if(asymptotic_bessel_large_x_limit(v, x))
+ {
+ if(kind&need_y)
+ {
+ Yv = asymptotic_bessel_y_large_x_2(v, x);
+ }
+ else
+ Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+ if(kind&need_j)
+ {
+ Jv = asymptotic_bessel_j_large_x_2(v, x);
+ }
+ else
+ Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+ }
+ else if((x > 8) && hankel_PQ(v, x, &p, &q, pol))
+ {
+ //
+ // Hankel approximation: note that this method works best when x
+ // is large, but in that case we end up calculating sines and cosines
+ // of large values, with horrendous resulting accuracy. It is fast though
+ // when it works....
+ //
+ // Normally we calculate sin/cos(chi) where:
+ //
+ // chi = x - fmod(T(v / 2 + 0.25f), T(2)) * boost::math::constants::pi<T>();
+ //
+ // But this introduces large errors, so use sin/cos addition formulae to
+ // improve accuracy:
+ //
+ T mod_v = fmod(T(v / 2 + 0.25f), T(2));
+ T sx = sin(x);
+ T cx = cos(x);
+ T sv = sin_pi(mod_v);
+ T cv = cos_pi(mod_v);
+
+ T sc = sx * cv - sv * cx; // == sin(chi);
+ T cc = cx * cv + sx * sv; // == cos(chi);
+ T chi = boost::math::constants::root_two<T>() / (boost::math::constants::root_pi<T>() * sqrt(x)); //sqrt(2 / (boost::math::constants::pi<T>() * x));
+ Yv = chi * (p * sc + q * cc);
+ Jv = chi * (p * cc - q * sc);
+ }
+ else if (x <= 2) // x in (0, 2]
+ {
+ if(temme_jy(u, x, &Yu, &Yu1, pol)) // Temme series
+ {
+ // domain error:
+ *J = *Y = Yu;
+ return 1;
+ }
+ prev = Yu;
+ current = Yu1;
+ T scale = 1;
+ policies::check_series_iterations<T>(function, n, pol);
+ for (k = 1; k <= n; k++) // forward recurrence for Y
+ {
+ T fact = 2 * (u + k) / x;
+ if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+ {
+ scale /= current;
+ prev /= current;
+ current = 1;
+ }
+ next = fact * current - prev;
+ prev = current;
+ current = next;
+ }
+ Yv = prev;
+ Yv1 = current;
+ if(kind&need_j)
+ {
+ CF1_jy(v, x, &fv, &s, pol); // continued fraction CF1_jy
+ Jv = scale * W / (Yv * fv - Yv1); // Wronskian relation
+ }
+ else
+ Jv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+ Yv_scale = scale;
+ }
+ else // x in (2, \infty)
+ {
+ // Get Y(u, x):
+
+ T ratio;
+ CF1_jy(v, x, &fv, &s, pol);
+ // tiny initial value to prevent overflow
+ T init = sqrt(tools::min_value<T>());
+ BOOST_MATH_INSTRUMENT_VARIABLE(init);
+ prev = fv * s * init;
+ current = s * init;
+ if(v < max_factorial<T>::value)
+ {
+ policies::check_series_iterations<T>(function, n, pol);
+ for (k = n; k > 0; k--) // backward recurrence for J
+ {
next = 2 * (u + k) * current / x - prev;
prev = current;
current = next;
- }
- ratio = (s * init) / current; // scaling ratio
- // can also call CF1_jy() to get fu, not much difference in precision
- fu = prev / current;
- }
- else
- {
- //
- // When v is large we may get overflow in this calculation
- // leading to NaN's and other nasty surprises:
- //
- bool over = false;
- for (k = n; k > 0; k--) // backward recurrence for J
- {
+ }
+ ratio = (s * init) / current; // scaling ratio
+ // can also call CF1_jy() to get fu, not much difference in precision
+ fu = prev / current;
+ }
+ else
+ {
+ //
+ // When v is large we may get overflow in this calculation
+ // leading to NaN's and other nasty surprises:
+ //
+ policies::check_series_iterations<T>(function, n, pol);
+ bool over = false;
+ for (k = n; k > 0; k--) // backward recurrence for J
+ {
T t = 2 * (u + k) / x;
- if(tools::max_value<T>() / t < current)
+ if((t > 1) && (tools::max_value<T>() / t < current))
{
over = true;
break;
@@ -469,87 +493,95 @@ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol)
next = t * current - prev;
prev = current;
current = next;
- }
- if(!over)
- {
- ratio = (s * init) / current; // scaling ratio
- // can also call CF1_jy() to get fu, not much difference in precision
- fu = prev / current;
- }
- else
- {
- ratio = 0;
- fu = 1;
- }
- }
- CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy
- T t = u / x - fu; // t = J'/J
- gamma = (p - t) / q;
- //
- // We can't allow gamma to cancel out to zero competely as it messes up
- // the subsequent logic. So pretend that one bit didn't cancel out
- // and set to a suitably small value. The only test case we've been able to
- // find for this, is when v = 8.5 and x = 4*PI.
- //
- if(gamma == 0)
- {
- gamma = u * tools::epsilon<T>() / x;
- }
- Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
-
- Jv = Ju * ratio; // normalization
-
- Yu = gamma * Ju;
- Yu1 = Yu * (u/x - p - q/gamma);
- }
- if(kind&need_y)
- {
- // compute Y:
- prev = Yu;
- current = Yu1;
- for (k = 1; k <= n; k++) // forward recurrence for Y
- {
- T fact = 2 * (u + k) / x;
- if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+ }
+ if(!over)
{
- prev /= current;
- Yv_scale /= current;
- current = 1;
+ ratio = (s * init) / current; // scaling ratio
+ // can also call CF1_jy() to get fu, not much difference in precision
+ fu = prev / current;
}
- next = fact * current - prev;
- prev = current;
- current = next;
- }
- Yv = prev;
- }
- else
- Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
- }
-
- if (reflect)
- {
- if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
- *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
- else
- *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale)); // reflection formula
- if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
- *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
- else
- *Y = sp * Jv + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
- }
- else
- {
- *J = Jv;
- if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
- *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
- else
- *Y = Yv / Yv_scale;
- }
-
- return 0;
-}
-
-} // namespace detail
+ else
+ {
+ ratio = 0;
+ fu = 1;
+ }
+ }
+ CF2_jy(u, x, &p, &q, pol); // continued fraction CF2_jy
+ T t = u / x - fu; // t = J'/J
+ gamma = (p - t) / q;
+ //
+ // We can't allow gamma to cancel out to zero competely as it messes up
+ // the subsequent logic. So pretend that one bit didn't cancel out
+ // and set to a suitably small value. The only test case we've been able to
+ // find for this, is when v = 8.5 and x = 4*PI.
+ //
+ if(gamma == 0)
+ {
+ gamma = u * tools::epsilon<T>() / x;
+ }
+ BOOST_MATH_INSTRUMENT_VARIABLE(current);
+ BOOST_MATH_INSTRUMENT_VARIABLE(W);
+ BOOST_MATH_INSTRUMENT_VARIABLE(q);
+ BOOST_MATH_INSTRUMENT_VARIABLE(gamma);
+ BOOST_MATH_INSTRUMENT_VARIABLE(p);
+ BOOST_MATH_INSTRUMENT_VARIABLE(t);
+ Ju = sign(current) * sqrt(W / (q + gamma * (p - t)));
+ BOOST_MATH_INSTRUMENT_VARIABLE(Ju);
+
+ Jv = Ju * ratio; // normalization
+
+ Yu = gamma * Ju;
+ Yu1 = Yu * (u/x - p - q/gamma);
+
+ if(kind&need_y)
+ {
+ // compute Y:
+ prev = Yu;
+ current = Yu1;
+ policies::check_series_iterations<T>(function, n, pol);
+ for (k = 1; k <= n; k++) // forward recurrence for Y
+ {
+ T fact = 2 * (u + k) / x;
+ if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+ {
+ prev /= current;
+ Yv_scale /= current;
+ current = 1;
+ }
+ next = fact * current - prev;
+ prev = current;
+ current = next;
+ }
+ Yv = prev;
+ }
+ else
+ Yv = std::numeric_limits<T>::quiet_NaN(); // any value will do, we're not using it.
+ }
+
+ if (reflect)
+ {
+ if((sp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(sp * Yv)))
+ *J = org_kind & need_j ? T(-sign(sp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+ else
+ *J = cp * Jv - (sp == 0 ? T(0) : T((sp * Yv) / Yv_scale)); // reflection formula
+ if((cp != 0) && (tools::max_value<T>() * fabs(Yv_scale) < fabs(cp * Yv)))
+ *Y = org_kind & need_y ? T(-sign(cp) * sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+ else
+ *Y = (sp != 0 ? sp * Jv : T(0)) + (cp == 0 ? T(0) : T((cp * Yv) / Yv_scale));
+ }
+ else
+ {
+ *J = Jv;
+ if(tools::max_value<T>() * fabs(Yv_scale) < fabs(Yv))
+ *Y = org_kind & need_y ? T(sign(Yv) * sign(Yv_scale) * policies::raise_overflow_error<T>(function, 0, pol)) : T(0);
+ else
+ *Y = Yv / Yv_scale;
+ }
+
+ return 0;
+ }
+
+ } // namespace detail
}} // namespaces
diff --git a/boost/math/special_functions/detail/bessel_jy_asym.hpp b/boost/math/special_functions/detail/bessel_jy_asym.hpp
index 0021f8c..81f6238 100644
--- a/boost/math/special_functions/detail/bessel_jy_asym.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_asym.hpp
@@ -21,61 +21,6 @@
namespace boost{ namespace math{ namespace detail{
template <class T>
-inline T asymptotic_bessel_j_large_x_P(T v, T x)
-{
- // A&S 9.2.9
- T s = 1;
- T mu = 4 * v * v;
- T ez2 = 8 * x;
- ez2 *= ez2;
- s -= (mu-1) * (mu-9) / (2 * ez2);
- s += (mu-1) * (mu-9) * (mu-25) * (mu - 49) / (24 * ez2 * ez2);
- return s;
-}
-
-template <class T>
-inline T asymptotic_bessel_j_large_x_Q(T v, T x)
-{
- // A&S 9.2.10
- T s = 0;
- T mu = 4 * v * v;
- T ez = 8*x;
- s += (mu-1) / ez;
- s -= (mu-1) * (mu-9) * (mu-25) / (6 * ez*ez*ez);
- return s;
-}
-
-template <class T>
-inline T asymptotic_bessel_j_large_x(T v, T x)
-{
- //
- // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
- //
- // Also A&S 9.2.5
- //
- BOOST_MATH_STD_USING // ADL of std names
- T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
- return sqrt(2 / (constants::pi<T>() * x))
- * (asymptotic_bessel_j_large_x_P(v, x) * cos(chi)
- - asymptotic_bessel_j_large_x_Q(v, x) * sin(chi));
-}
-
-template <class T>
-inline T asymptotic_bessel_y_large_x(T v, T x)
-{
- //
- // See http://functions.wolfram.com/BesselAiryStruveFunctions/BesselJ/06/02/02/0001/
- //
- // Also A&S 9.2.5
- //
- BOOST_MATH_STD_USING // ADL of std names
- T chi = fabs(x) - constants::pi<T>() * (2 * v + 1) / 4;
- return sqrt(2 / (constants::pi<T>() * x))
- * (asymptotic_bessel_j_large_x_P(v, x) * sin(chi)
- - asymptotic_bessel_j_large_x_Q(v, x) * cos(chi));
-}
-
-template <class T>
inline T asymptotic_bessel_amplitude(T v, T x)
{
// Calculate the amplitude of J(v, x) and Y(v, x) for large
@@ -99,13 +44,14 @@ T asymptotic_bessel_phase_mx(T v, T x)
//
// Calculate the phase of J(v, x) and Y(v, x) for large x.
// See A&S 9.2.29.
- // Note that the result returned is the phase less x.
+ // Note that the result returned is the phase less (x - PI(v/2 + 1/4))
+ // which we'll factor in later when we calculate the sines/cosines of the result:
//
T mu = 4 * v * v;
T denom = 4 * x;
T denom_mult = denom * denom;
- T s = -constants::pi<T>() * (v / 2 + 0.25f);
+ T s = 0;
s += (mu - 1) / (2 * denom);
denom *= denom_mult;
s += (mu - 1) * (mu - 25) / (6 * denom);
@@ -127,10 +73,16 @@ inline T asymptotic_bessel_y_large_x_2(T v, T x)
BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
BOOST_MATH_INSTRUMENT_VARIABLE(phase);
//
- // Calculate the sine of the phase, using:
- // sin(x+p) = sin(x)cos(p) + cos(x)sin(p)
+ // Calculate the sine of the phase, using
+ // sine/cosine addition rules to factor in
+ // the x - PI(v/2 + 1/4) term not added to the
+ // phase when we calculated it.
//
- T sin_phase = sin(phase) * cos(x) + cos(phase) * sin(x);
+ T cx = cos(x);
+ T sx = sin(x);
+ T ci = cos_pi(v / 2 + 0.25f);
+ T si = sin_pi(v / 2 + 0.25f);
+ T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
BOOST_MATH_INSTRUMENT_CODE(sin(phase));
BOOST_MATH_INSTRUMENT_CODE(cos(x));
BOOST_MATH_INSTRUMENT_CODE(cos(phase));
@@ -149,101 +101,39 @@ inline T asymptotic_bessel_j_large_x_2(T v, T x)
BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
BOOST_MATH_INSTRUMENT_VARIABLE(phase);
//
- // Calculate the sine of the phase, using:
- // cos(x+p) = cos(x)cos(p) - sin(x)sin(p)
+ // Calculate the sine of the phase, using
+ // sine/cosine addition rules to factor in
+ // the x - PI(v/2 + 1/4) term not added to the
+ // phase when we calculated it.
//
BOOST_MATH_INSTRUMENT_CODE(cos(phase));
BOOST_MATH_INSTRUMENT_CODE(cos(x));
BOOST_MATH_INSTRUMENT_CODE(sin(phase));
BOOST_MATH_INSTRUMENT_CODE(sin(x));
- T sin_phase = cos(phase) * cos(x) - sin(phase) * sin(x);
+ T cx = cos(x);
+ T sx = sin(x);
+ T ci = cos_pi(v / 2 + 0.25f);
+ T si = sin_pi(v / 2 + 0.25f);
+ T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
return sin_phase * ampl;
}
-//
-// Various limits for the J and Y asymptotics
-// (the asympotic expansions are safe to use if
-// x is less than the limit given).
-// We assume that if we don't use these expansions then the
-// error will likely be >100eps, so the limits given are chosen
-// to lead to < 100eps truncation error.
-//
template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<0>&)
+inline bool asymptotic_bessel_large_x_limit(const T& v, const T& x)
{
- // default case:
BOOST_MATH_STD_USING
- return 2.25 / pow(100 * tools::epsilon<T>() / T(0.001f), T(0.2f));
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<53>&)
-{
- // double case:
- return 304 /*780*/;
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<64>&)
-{
- // 80-bit extended-double case:
- return 1552 /*3500*/;
-}
-template <class T>
-inline T asymptotic_bessel_y_limit(const mpl::int_<113>&)
-{
- // 128-bit long double case:
- return 1245243 /*3128000*/;
-}
-
-template <class T, class Policy>
-struct bessel_asymptotic_tag
-{
- typedef typename policies::precision<T, Policy>::type precision_type;
- typedef typename mpl::if_<
- mpl::or_<
- mpl::equal_to<precision_type, mpl::int_<0> >,
- mpl::greater<precision_type, mpl::int_<113> > >,
- mpl::int_<0>,
- typename mpl::if_<
- mpl::greater<precision_type, mpl::int_<64> >,
- mpl::int_<113>,
- typename mpl::if_<
- mpl::greater<precision_type, mpl::int_<53> >,
- mpl::int_<64>,
- mpl::int_<53>
- >::type
- >::type
- >::type type;
-};
-
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<0>&)
-{
- // default case:
- BOOST_MATH_STD_USING
- T v2 = (std::max)(T(3), T(v * v));
- return v2 / pow(100 * tools::epsilon<T>() / T(2e-5f), T(0.17f));
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<53>&)
-{
- // double case:
- T v2 = (std::max)(T(3), T(v * v));
- return v2 * 33 /*73*/;
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<64>&)
-{
- // 80-bit extended-double case:
- T v2 = (std::max)(T(3), T(v * v));
- return v2 * 121 /*266*/;
-}
-template <class T>
-inline T asymptotic_bessel_j_limit(const T& v, const mpl::int_<113>&)
-{
- // 128-bit long double case:
- T v2 = (std::max)(T(3), T(v * v));
- return v2 * 39154 /*85700*/;
+ //
+ // Determines if x is large enough compared to v to take the asymptotic
+ // forms above. From A&S 9.2.28 we require:
+ // v < x * eps^1/8
+ // and from A&S 9.2.29 we require:
+ // v^12/10 < 1.5 * x * eps^1/10
+ // using the former seems to work OK in practice with broadly similar
+ // error rates either side of the divide for v < 10000.
+ // At double precision eps^1/8 ~= 0.01.
+ //
+ return (std::max)(T(fabs(v)), T(1)) < x * sqrt(tools::forth_root_epsilon<T>());
}
template <class T, class Policy>
diff --git a/boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp b/boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp
new file mode 100644
index 0000000..bdbfb9d
--- /dev/null
+++ b/boost/math/special_functions/detail/bessel_jy_derivatives_asym.hpp
@@ -0,0 +1,141 @@
+// Copyright (c) 2013 Anton Bikineev
+// 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)
+
+//
+// This is a partial header, do not include on it's own!!!
+//
+// Contains asymptotic expansions for derivatives of Bessel J(v,x) and Y(v,x)
+// functions, as x -> INF.
+#ifndef BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
+#define BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class T>
+inline T asymptotic_bessel_derivative_amplitude(T v, T x)
+{
+ // Calculate the amplitude for J'(v,x) and I'(v,x)
+ // for large x: see A&S 9.2.30.
+ BOOST_MATH_STD_USING
+ T s = 1;
+ const T mu = 4 * v * v;
+ T txq = 2 * x;
+ txq *= txq;
+
+ s -= (mu - 3) / (2 * txq);
+ s -= ((mu - 1) * (mu - 45)) / (txq * txq * 8);
+
+ return sqrt(s * 2 / (boost::math::constants::pi<T>() * x));
+}
+
+template <class T>
+inline T asymptotic_bessel_derivative_phase_mx(T v, T x)
+{
+ // Calculate the phase of J'(v, x) and Y'(v, x) for large x.
+ // See A&S 9.2.31.
+ // Note that the result returned is the phase less (x - PI(v/2 - 1/4))
+ // which we'll factor in later when we calculate the sines/cosines of the result:
+ const T mu = 4 * v * v;
+ const T mu2 = mu * mu;
+ const T mu3 = mu2 * mu;
+ T denom = 4 * x;
+ T denom_mult = denom * denom;
+
+ T s = 0;
+ s += (mu + 3) / (2 * denom);
+ denom *= denom_mult;
+ s += (mu2 + (46 * mu) - 63) / (6 * denom);
+ denom *= denom_mult;
+ s += (mu3 + (185 * mu2) - (2053 * mu) + 1899) / (5 * denom);
+ return s;
+}
+
+template <class T>
+inline T asymptotic_bessel_y_derivative_large_x_2(T v, T x)
+{
+ // See A&S 9.2.20.
+ BOOST_MATH_STD_USING
+ // Get the phase and amplitude:
+ const T ampl = asymptotic_bessel_derivative_amplitude(v, x);
+ const T phase = asymptotic_bessel_derivative_phase_mx(v, x);
+ BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
+ BOOST_MATH_INSTRUMENT_VARIABLE(phase);
+ //
+ // Calculate the sine of the phase, using
+ // sine/cosine addition rules to factor in
+ // the x - PI(v/2 - 1/4) term not added to the
+ // phase when we calculated it.
+ //
+ const T cx = cos(x);
+ const T sx = sin(x);
+ const T vd2shifted = (v / 2) - 0.25f;
+ const T ci = cos_pi(vd2shifted);
+ const T si = sin_pi(vd2shifted);
+ const T sin_phase = sin(phase) * (cx * ci + sx * si) + cos(phase) * (sx * ci - cx * si);
+ BOOST_MATH_INSTRUMENT_CODE(sin(phase));
+ BOOST_MATH_INSTRUMENT_CODE(cos(x));
+ BOOST_MATH_INSTRUMENT_CODE(cos(phase));
+ BOOST_MATH_INSTRUMENT_CODE(sin(x));
+ return sin_phase * ampl;
+}
+
+template <class T>
+inline T asymptotic_bessel_j_derivative_large_x_2(T v, T x)
+{
+ // See A&S 9.2.20.
+ BOOST_MATH_STD_USING
+ // Get the phase and amplitude:
+ const T ampl = asymptotic_bessel_derivative_amplitude(v, x);
+ const T phase = asymptotic_bessel_derivative_phase_mx(v, x);
+ BOOST_MATH_INSTRUMENT_VARIABLE(ampl);
+ BOOST_MATH_INSTRUMENT_VARIABLE(phase);
+ //
+ // Calculate the sine of the phase, using
+ // sine/cosine addition rules to factor in
+ // the x - PI(v/2 - 1/4) term not added to the
+ // phase when we calculated it.
+ //
+ BOOST_MATH_INSTRUMENT_CODE(cos(phase));
+ BOOST_MATH_INSTRUMENT_CODE(cos(x));
+ BOOST_MATH_INSTRUMENT_CODE(sin(phase));
+ BOOST_MATH_INSTRUMENT_CODE(sin(x));
+ const T cx = cos(x);
+ const T sx = sin(x);
+ const T vd2shifted = (v / 2) - 0.25f;
+ const T ci = cos_pi(vd2shifted);
+ const T si = sin_pi(vd2shifted);
+ const T sin_phase = cos(phase) * (cx * ci + sx * si) - sin(phase) * (sx * ci - cx * si);
+ BOOST_MATH_INSTRUMENT_VARIABLE(sin_phase);
+ return sin_phase * ampl;
+}
+
+template <class T>
+inline bool asymptotic_bessel_derivative_large_x_limit(const T& v, const T& x)
+{
+ BOOST_MATH_STD_USING
+ //
+ // This function is the copy of math::asymptotic_bessel_large_x_limit
+ // It means that we use the same rules for determining how x is large
+ // compared to v.
+ //
+ // Determines if x is large enough compared to v to take the asymptotic
+ // forms above. From A&S 9.2.28 we require:
+ // v < x * eps^1/8
+ // and from A&S 9.2.29 we require:
+ // v^12/10 < 1.5 * x * eps^1/10
+ // using the former seems to work OK in practice with broadly similar
+ // error rates either side of the divide for v < 10000.
+ // At double precision eps^1/8 ~= 0.01.
+ //
+ return (std::max)(T(fabs(v)), T(1)) < x * sqrt(boost::math::tools::forth_root_epsilon<T>());
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_SF_DETAIL_BESSEL_JY_DERIVATIVES_ASYM_HPP
diff --git a/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp b/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
new file mode 100644
index 0000000..0dc68fc
--- /dev/null
+++ b/boost/math/special_functions/detail/bessel_jy_derivatives_series.hpp
@@ -0,0 +1,220 @@
+// Copyright (c) 2013 Anton Bikineev
+// 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_BESSEL_JY_DERIVATIVES_SERIES_HPP
+#define BOOST_MATH_BESSEL_JY_DERIVATIVES_SERIES_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class T, class Policy>
+struct bessel_j_derivative_small_z_series_term
+{
+ typedef T result_type;
+
+ bessel_j_derivative_small_z_series_term(T v_, T x)
+ : N(0), v(v_), term(1), mult(x / 2)
+ {
+ mult *= -mult;
+ // iterate if v == 0; otherwise result of
+ // first term is 0 and tools::sum_series stops
+ if (v == 0)
+ iterate();
+ }
+ T operator()()
+ {
+ T r = term * (v + 2 * N);
+ iterate();
+ return r;
+ }
+private:
+ void iterate()
+ {
+ ++N;
+ term *= mult / (N * (N + v));
+ }
+ unsigned N;
+ T v;
+ T term;
+ T mult;
+};
+//
+// Series evaluation for BesselJ'(v, z) as z -> 0.
+// It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselJ/06/01/04/01/01/0003/
+// Converges rapidly for all z << v.
+//
+template <class T, class Policy>
+inline T bessel_j_derivative_small_z_series(T v, T x, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+ T prefix;
+ if (v < boost::math::max_factorial<T>::value)
+ {
+ prefix = pow(x / 2, v - 1) / 2 / boost::math::tgamma(v + 1, pol);
+ }
+ else
+ {
+ prefix = (v - 1) * log(x / 2) - constants::ln_two<T>() - boost::math::lgamma(v + 1, pol);
+ prefix = exp(prefix);
+ }
+ if (0 == prefix)
+ return prefix;
+
+ bessel_j_derivative_small_z_series_term<T, Policy> s(v, x);
+ boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ T zero = 0;
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+ boost::math::policies::check_series_iterations<T>("boost::math::bessel_j_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
+ return prefix * result;
+}
+
+template <class T, class Policy>
+struct bessel_y_derivative_small_z_series_term_a
+{
+ typedef T result_type;
+
+ bessel_y_derivative_small_z_series_term_a(T v_, T x)
+ : N(0), v(v_)
+ {
+ mult = x / 2;
+ mult *= -mult;
+ term = 1;
+ }
+ T operator()()
+ {
+ T r = term * (-v + 2 * N);
+ ++N;
+ term *= mult / (N * (N - v));
+ return r;
+ }
+private:
+ unsigned N;
+ T v;
+ T mult;
+ T term;
+};
+
+template <class T, class Policy>
+struct bessel_y_derivative_small_z_series_term_b
+{
+ typedef T result_type;
+
+ bessel_y_derivative_small_z_series_term_b(T v_, T x)
+ : N(0), v(v_)
+ {
+ mult = x / 2;
+ mult *= -mult;
+ term = 1;
+ }
+ T operator()()
+ {
+ T r = term * (v + 2 * N);
+ ++N;
+ term *= mult / (N * (N + v));
+ return r;
+ }
+private:
+ unsigned N;
+ T v;
+ T mult;
+ T term;
+};
+//
+// Series form for BesselY' as z -> 0,
+// It's derivative of http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/01/0003/
+// This series is only useful when the second term is small compared to the first
+// otherwise we get catestrophic cancellation errors.
+//
+// Approximating tgamma(v) by v^v, and assuming |tgamma(-z)| < eps we end up requiring:
+// eps/2 * v^v(x/2)^-v > (x/2)^v or log(eps/2) > v log((x/2)^2/v)
+//
+template <class T, class Policy>
+inline T bessel_y_derivative_small_z_series(T v, T x, const Policy& pol)
+{
+ BOOST_MATH_STD_USING
+ static const char* function = "bessel_y_derivative_small_z_series<%1%>(%1%,%1%)";
+ T prefix;
+ T gam;
+ T p = log(x / 2);
+ T scale = 1;
+ bool need_logs = (v >= boost::math::max_factorial<T>::value) || (boost::math::tools::log_max_value<T>() / v < fabs(p));
+ if (!need_logs)
+ {
+ gam = boost::math::tgamma(v, pol);
+ p = pow(x / 2, v + 1) * 2;
+ if (boost::math::tools::max_value<T>() * p < gam)
+ {
+ scale /= gam;
+ gam = 1;
+ if (boost::math::tools::max_value<T>() * p < gam)
+ {
+ return -boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+ }
+ }
+ prefix = -gam / (boost::math::constants::pi<T>() * p);
+ }
+ else
+ {
+ gam = boost::math::lgamma(v, pol);
+ p = (v + 1) * p + constants::ln_two<T>();
+ prefix = gam - log(boost::math::constants::pi<T>()) - p;
+ if (boost::math::tools::log_max_value<T>() < prefix)
+ {
+ prefix -= log(boost::math::tools::max_value<T>() / 4);
+ scale /= (boost::math::tools::max_value<T>() / 4);
+ if (boost::math::tools::log_max_value<T>() < prefix)
+ {
+ return -boost::math::policies::raise_overflow_error<T>(function, 0, pol);
+ }
+ }
+ prefix = -exp(prefix);
+ }
+ bessel_y_derivative_small_z_series_term_a<T, Policy> s(v, x);
+ boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ T zero = 0;
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+ T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+ boost::math::policies::check_series_iterations<T>("boost::math::bessel_y_derivative_small_z_series<%1%>(%1%,%1%)", max_iter, pol);
+ result *= prefix;
+
+ p = pow(x / 2, v - 1) / 2;
+ if (!need_logs)
+ {
+ prefix = boost::math::tgamma(-v, pol) * boost::math::cos_pi(v) * p / boost::math::constants::pi<T>();
+ }
+ else
+ {
+ int sgn;
+ prefix = boost::math::lgamma(-v, &sgn, pol) + (v - 1) * log(x / 2) - constants::ln_two<T>();
+ prefix = exp(prefix) * sgn / boost::math::constants::pi<T>();
+ }
+ bessel_y_derivative_small_z_series_term_b<T, Policy> s2(v, x);
+ max_iter = boost::math::policies::get_max_series_iterations<Policy>();
+#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter, zero);
+#else
+ T b = boost::math::tools::sum_series(s2, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
+#endif
+ result += scale * prefix * b;
+ return result;
+}
+
+// Calculating of BesselY'(v,x) with small x (x < epsilon) and integer x using derivatives
+// of formulas in http://functions.wolfram.com/Bessel-TypeFunctions/BesselY/06/01/04/01/02/
+// seems to lose precision. Instead using linear combination of regular Bessel is preferred.
+
+}}} // namespaces
+
+#endif // BOOST_MATH_BESSEL_JY_DERIVATVIES_SERIES_HPP
diff --git a/boost/math/special_functions/detail/bessel_jy_series.hpp b/boost/math/special_functions/detail/bessel_jy_series.hpp
index b926366..d50bef8 100644
--- a/boost/math/special_functions/detail/bessel_jy_series.hpp
+++ b/boost/math/special_functions/detail/bessel_jy_series.hpp
@@ -194,9 +194,9 @@ inline T bessel_y_small_z_series(T v, T x, T* pscale, const Policy& pol)
}
else
{
- int s;
- prefix = boost::math::lgamma(-v, &s, pol) + p;
- prefix = exp(prefix) * s / constants::pi<T>();
+ int sgn;
+ prefix = boost::math::lgamma(-v, &sgn, pol) + p;
+ prefix = exp(prefix) * sgn / constants::pi<T>();
}
bessel_y_small_z_series_term_b<T, Policy> s2(v, x);
max_iter = policies::get_max_series_iterations<Policy>();
@@ -235,7 +235,7 @@ T bessel_yn_small_z(int n, T z, T* scale, const Policy& pol)
{
return (z * z) / (4 * constants::pi<T>()) * log(z / 2)
- (4 / (constants::pi<T>() * z * z))
- - ((z * z) / (8 * constants::pi<T>())) * (3/2 - 2 * constants::euler<T>());
+ - ((z * z) / (8 * constants::pi<T>())) * (T(3)/2 - 2 * constants::euler<T>());
}
else
{
diff --git a/boost/math/special_functions/detail/bessel_jy_zero.hpp b/boost/math/special_functions/detail/bessel_jy_zero.hpp
new file mode 100644
index 0000000..ecd8696
--- /dev/null
+++ b/boost/math/special_functions/detail/bessel_jy_zero.hpp
@@ -0,0 +1,617 @@
+// Copyright (c) 2013 Christopher Kormanyos
+// 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)
+//
+// This work is based on an earlier work:
+// "Algorithm 910: A Portable C++ Multiple-Precision System for Special-Function Calculations",
+// in ACM TOMS, {VOL 37, ISSUE 4, (February 2011)} (C) ACM, 2011. http://doi.acm.org/10.1145/1916461.1916469
+//
+// This header contains implementation details for estimating the zeros
+// of cylindrical Bessel and Neumann functions on the positive real axis.
+// Support is included for both positive as well as negative order.
+// Various methods are used to estimate the roots. These include
+// empirical curve fitting and McMahon's asymptotic approximation
+// for small order, uniform asymptotic expansion for large order,
+// and iteration and root interlacing for negative order.
+//
+#ifndef _BESSEL_JY_ZERO_2013_01_18_HPP_
+ #define _BESSEL_JY_ZERO_2013_01_18_HPP_
+
+ #include <algorithm>
+ #include <boost/math/constants/constants.hpp>
+ #include <boost/math/special_functions/math_fwd.hpp>
+ #include <boost/math/special_functions/cbrt.hpp>
+ #include <boost/math/special_functions/detail/airy_ai_bi_zero.hpp>
+
+ namespace boost { namespace math {
+ namespace detail
+ {
+ namespace bessel_zero
+ {
+ template<class T>
+ T equation_nist_10_21_19(const T& v, const T& a)
+ {
+ // Get the initial estimate of the m'th root of Jv or Yv.
+ // This subroutine is used for the order m with m > 1.
+ // The order m has been used to create the input parameter a.
+
+ // This is Eq. 10.21.19 in the NIST Handbook.
+ const T mu = (v * v) * 4U;
+ const T mu_minus_one = mu - T(1);
+ const T eight_a_inv = T(1) / (a * 8U);
+ const T eight_a_inv_squared = eight_a_inv * eight_a_inv;
+
+ const T term3 = ((mu_minus_one * 4U) * ((mu * 7U) - T(31U) )) / 3U;
+ const T term5 = ((mu_minus_one * 32U) * ((((mu * 83U) - T(982U) ) * mu) + T(3779U) )) / 15U;
+ const T term7 = ((mu_minus_one * 64U) * ((((((mu * 6949U) - T(153855UL)) * mu) + T(1585743UL)) * mu) - T(6277237UL))) / 105U;
+
+ return a + (((( - term7
+ * eight_a_inv_squared - term5)
+ * eight_a_inv_squared - term3)
+ * eight_a_inv_squared - mu_minus_one)
+ * eight_a_inv);
+ }
+
+ template<typename T>
+ class equation_as_9_3_39_and_its_derivative
+ {
+ public:
+ equation_as_9_3_39_and_its_derivative(const T& zt) : zeta(zt) { }
+
+ boost::math::tuple<T, T> operator()(const T& z) const
+ {
+ BOOST_MATH_STD_USING // ADL of std names, needed for acos, sqrt.
+
+ // Return the function of zeta that is implicitly defined
+ // in A&S Eq. 9.3.39 as a function of z. The function is
+ // returned along with its derivative with respect to z.
+
+ const T zsq_minus_one_sqrt = sqrt((z * z) - T(1));
+
+ const T the_function(
+ zsq_minus_one_sqrt
+ - ( acos(T(1) / z) + ((T(2) / 3U) * (zeta * sqrt(zeta)))));
+
+ const T its_derivative(zsq_minus_one_sqrt / z);
+
+ return boost::math::tuple<T, T>(the_function, its_derivative);
+ }
+
+ private:
+ const equation_as_9_3_39_and_its_derivative& operator=(const equation_as_9_3_39_and_its_derivative&);
+ const T zeta;
+ };
+
+ template<class T>
+ static T equation_as_9_5_26(const T& v, const T& ai_bi_root)
+ {
+ BOOST_MATH_STD_USING // ADL of std names, needed for log, sqrt.
+
+ // Obtain the estimate of the m'th zero of Jv or Yv.
+ // The order m has been used to create the input parameter ai_bi_root.
+ // Here, v is larger than about 2.2. The estimate is computed
+ // from Abramowitz and Stegun Eqs. 9.5.22 and 9.5.26, page 371.
+ //
+ // The inversion of z as a function of zeta is mentioned in the text
+ // following A&S Eq. 9.5.26. Here, we accomplish the inversion by
+ // performing a Taylor expansion of Eq. 9.3.39 for large z to order 2
+ // and solving the resulting quadratic equation, thereby taking
+ // the positive root of the quadratic.
+ // In other words: (2/3)(-zeta)^(3/2) approx = z + 1/(2z) - pi/2.
+ // This leads to: z^2 - [(2/3)(-zeta)^(3/2) + pi/2]z + 1/2 = 0.
+ //
+ // With this initial estimate, Newton-Raphson iteration is used
+ // to refine the value of the estimate of the root of z
+ // as a function of zeta.
+
+ const T v_pow_third(boost::math::cbrt(v));
+ const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
+
+ // Obtain zeta using the order v combined with the m'th root of
+ // an airy function, as shown in A&S Eq. 9.5.22.
+ const T zeta = v_pow_minus_two_thirds * (-ai_bi_root);
+
+ const T zeta_sqrt = sqrt(zeta);
+
+ // Set up a quadratic equation based on the Taylor series
+ // expansion mentioned above.
+ const T b = -((((zeta * zeta_sqrt) * 2U) / 3U) + boost::math::constants::half_pi<T>());
+
+ // Solve the quadratic equation, taking the positive root.
+ const T z_estimate = (-b + sqrt((b * b) - T(2))) / 2U;
+
+ // Establish the range, the digits, and the iteration limit
+ // for the upcoming root-finding.
+ const T range_zmin = (std::max<T>)(z_estimate - T(1), T(1));
+ const T range_zmax = z_estimate + T(1);
+
+ const int my_digits10 = static_cast<int>(static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
+
+ // Select the maximum allowed iterations based on the number
+ // of decimal digits in the numeric type T, being at least 12.
+ const boost::uintmax_t iterations_allowed = static_cast<boost::uintmax_t>((std::max)(12, my_digits10 * 2));
+
+ boost::uintmax_t iterations_used = iterations_allowed;
+
+ // Calculate the root of z as a function of zeta.
+ const T z = boost::math::tools::newton_raphson_iterate(
+ boost::math::detail::bessel_zero::equation_as_9_3_39_and_its_derivative<T>(zeta),
+ z_estimate,
+ range_zmin,
+ range_zmax,
+ (std::min)(boost::math::tools::digits<T>(), boost::math::tools::digits<float>()),
+ iterations_used);
+
+ static_cast<void>(iterations_used);
+
+ // Continue with the implementation of A&S Eq. 9.3.39.
+ const T zsq_minus_one = (z * z) - T(1);
+ const T zsq_minus_one_sqrt = sqrt(zsq_minus_one);
+
+ // This is A&S Eq. 9.3.42.
+ const T b0_term_5_24 = T(5) / ((zsq_minus_one * zsq_minus_one_sqrt) * 24U);
+ const T b0_term_1_8 = T(1) / ( zsq_minus_one_sqrt * 8U);
+ const T b0_term_5_48 = T(5) / ((zeta * zeta) * 48U);
+
+ const T b0 = -b0_term_5_48 + ((b0_term_5_24 + b0_term_1_8) / zeta_sqrt);
+
+ // This is the second line of A&S Eq. 9.5.26 for f_k with k = 1.
+ const T f1 = ((z * zeta_sqrt) * b0) / zsq_minus_one_sqrt;
+
+ // This is A&S Eq. 9.5.22 expanded to k = 1 (i.e., one term in the series).
+ return (v * z) + (f1 / v);
+ }
+
+ namespace cyl_bessel_j_zero_detail
+ {
+ template<class T>
+ T equation_nist_10_21_40_a(const T& v)
+ {
+ const T v_pow_third(boost::math::cbrt(v));
+ const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
+
+ return v * ((((( + T(0.043)
+ * v_pow_minus_two_thirds - T(0.0908))
+ * v_pow_minus_two_thirds - T(0.00397))
+ * v_pow_minus_two_thirds + T(1.033150))
+ * v_pow_minus_two_thirds + T(1.8557571))
+ * v_pow_minus_two_thirds + T(1));
+ }
+
+ template<class T, class Policy>
+ class function_object_jv
+ {
+ public:
+ function_object_jv(const T& v,
+ const Policy& pol) : my_v(v),
+ my_pol(pol) { }
+
+ T operator()(const T& x) const
+ {
+ return boost::math::cyl_bessel_j(my_v, x, my_pol);
+ }
+
+ private:
+ const T my_v;
+ const Policy& my_pol;
+ const function_object_jv& operator=(const function_object_jv&);
+ };
+
+ template<class T, class Policy>
+ class function_object_jv_and_jv_prime
+ {
+ public:
+ function_object_jv_and_jv_prime(const T& v,
+ const bool order_is_zero,
+ const Policy& pol) : my_v(v),
+ my_order_is_zero(order_is_zero),
+ my_pol(pol) { }
+
+ boost::math::tuple<T, T> operator()(const T& x) const
+ {
+ // Obtain Jv(x) and Jv'(x).
+ // Chris's original code called the Bessel function implementation layer direct,
+ // but that circumvented optimizations for integer-orders. Call the documented
+ // top level functions instead, and let them sort out which implementation to use.
+ T j_v;
+ T j_v_prime;
+
+ if(my_order_is_zero)
+ {
+ j_v = boost::math::cyl_bessel_j(0, x, my_pol);
+ j_v_prime = -boost::math::cyl_bessel_j(1, x, my_pol);
+ }
+ else
+ {
+ j_v = boost::math::cyl_bessel_j( my_v, x, my_pol);
+ const T j_v_m1 (boost::math::cyl_bessel_j(T(my_v - 1), x, my_pol));
+ j_v_prime = j_v_m1 - ((my_v * j_v) / x);
+ }
+
+ // Return a tuple containing both Jv(x) and Jv'(x).
+ return boost::math::make_tuple(j_v, j_v_prime);
+ }
+
+ private:
+ const T my_v;
+ const bool my_order_is_zero;
+ const Policy& my_pol;
+ const function_object_jv_and_jv_prime& operator=(const function_object_jv_and_jv_prime&);
+ };
+
+ template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
+
+ template<class T, class Policy>
+ T initial_guess(const T& v, const int m, const Policy& pol)
+ {
+ BOOST_MATH_STD_USING // ADL of std names, needed for floor.
+
+ // Compute an estimate of the m'th root of cyl_bessel_j.
+
+ T guess;
+
+ // There is special handling for negative order.
+ if(v < 0)
+ {
+ if((m == 1) && (v > -0.5F))
+ {
+ // For small, negative v, use the results of empirical curve fitting.
+ // Mathematica(R) session for the coefficients:
+ // Table[{n, BesselJZero[n, 1]}, {n, -(1/2), 0, 1/10}]
+ // N[%, 20]
+ // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
+ guess = ((((( - T(0.2321156900729)
+ * v - T(0.1493247777488))
+ * v - T(0.15205419167239))
+ * v + T(0.07814930561249))
+ * v - T(0.17757573537688))
+ * v + T(1.542805677045663))
+ * v + T(2.40482555769577277);
+
+ return guess;
+ }
+
+ // Create the positive order and extract its positive floor integer part.
+ const T vv(-v);
+ const T vv_floor(floor(vv));
+
+ // The to-be-found root is bracketed by the roots of the
+ // Bessel function whose reflected, positive integer order
+ // is less than, but nearest to vv.
+
+ T root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m, pol);
+ T root_lo;
+
+ if(m == 1)
+ {
+ // The estimate of the first root for negative order is found using
+ // an adaptive range-searching algorithm.
+ root_lo = T(root_hi - 0.1F);
+
+ const bool hi_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_hi, pol) < 0);
+
+ while((root_lo > boost::math::tools::epsilon<T>()))
+ {
+ const bool lo_end_of_bracket_is_negative = (boost::math::cyl_bessel_j(v, root_lo, pol) < 0);
+
+ if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
+ {
+ break;
+ }
+
+ root_hi = root_lo;
+
+ // Decrease the lower end of the bracket using an adaptive algorithm.
+ if(root_lo > 0.5F)
+ {
+ root_lo -= 0.5F;
+ }
+ else
+ {
+ root_lo *= 0.75F;
+ }
+ }
+ }
+ else
+ {
+ root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(vv_floor, m - 1, pol);
+ }
+
+ // Perform several steps of bisection iteration to refine the guess.
+ boost::uintmax_t number_of_iterations(12U);
+
+ // Do the bisection iteration.
+ const boost::math::tuple<T, T> guess_pair =
+ boost::math::tools::bisect(
+ boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::function_object_jv<T, Policy>(v, pol),
+ root_lo,
+ root_hi,
+ boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::my_bisection_unreachable_tolerance<T>,
+ number_of_iterations);
+
+ return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
+ }
+
+ if(m == 1U)
+ {
+ // Get the initial estimate of the first root.
+
+ if(v < 2.2F)
+ {
+ // For small v, use the results of empirical curve fitting.
+ // Mathematica(R) session for the coefficients:
+ // Table[{n, BesselJZero[n, 1]}, {n, 0, 22/10, 1/10}]
+ // N[%, 20]
+ // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
+ guess = ((((( - T(0.0008342379046010)
+ * v + T(0.007590035637410))
+ * v - T(0.030640914772013))
+ * v + T(0.078232088020106))
+ * v - T(0.169668712590620))
+ * v + T(1.542187960073750))
+ * v + T(2.4048359915254634);
+ }
+ else
+ {
+ // For larger v, use the first line of Eqs. 10.21.40 in the NIST Handbook.
+ guess = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::equation_nist_10_21_40_a(v);
+ }
+ }
+ else
+ {
+ if(v < 2.2F)
+ {
+ // Use Eq. 10.21.19 in the NIST Handbook.
+ const T a(((v + T(m * 2U)) - T(0.5)) * boost::math::constants::half_pi<T>());
+
+ guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
+ }
+ else
+ {
+ // Get an estimate of the m'th root of airy_ai.
+ const T airy_ai_root(boost::math::detail::airy_zero::airy_ai_zero_detail::initial_guess<T>(m));
+
+ // Use Eq. 9.5.26 in the A&S Handbook.
+ guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_ai_root);
+ }
+ }
+
+ return guess;
+ }
+ } // namespace cyl_bessel_j_zero_detail
+
+ namespace cyl_neumann_zero_detail
+ {
+ template<class T>
+ T equation_nist_10_21_40_b(const T& v)
+ {
+ const T v_pow_third(boost::math::cbrt(v));
+ const T v_pow_minus_two_thirds(T(1) / (v_pow_third * v_pow_third));
+
+ return v * ((((( - T(0.001)
+ * v_pow_minus_two_thirds - T(0.0060))
+ * v_pow_minus_two_thirds + T(0.01198))
+ * v_pow_minus_two_thirds + T(0.260351))
+ * v_pow_minus_two_thirds + T(0.9315768))
+ * v_pow_minus_two_thirds + T(1));
+ }
+
+ template<class T, class Policy>
+ class function_object_yv
+ {
+ public:
+ function_object_yv(const T& v,
+ const Policy& pol) : my_v(v),
+ my_pol(pol) { }
+
+ T operator()(const T& x) const
+ {
+ return boost::math::cyl_neumann(my_v, x, my_pol);
+ }
+
+ private:
+ const T my_v;
+ const Policy& my_pol;
+ const function_object_yv& operator=(const function_object_yv&);
+ };
+
+ template<class T, class Policy>
+ class function_object_yv_and_yv_prime
+ {
+ public:
+ function_object_yv_and_yv_prime(const T& v,
+ const Policy& pol) : my_v(v),
+ my_pol(pol) { }
+
+ boost::math::tuple<T, T> operator()(const T& x) const
+ {
+ const T half_epsilon(boost::math::tools::epsilon<T>() / 2U);
+
+ const bool order_is_zero = ((my_v > -half_epsilon) && (my_v < +half_epsilon));
+
+ // Obtain Yv(x) and Yv'(x).
+ // Chris's original code called the Bessel function implementation layer direct,
+ // but that circumvented optimizations for integer-orders. Call the documented
+ // top level functions instead, and let them sort out which implementation to use.
+ T y_v;
+ T y_v_prime;
+
+ if(order_is_zero)
+ {
+ y_v = boost::math::cyl_neumann(0, x, my_pol);
+ y_v_prime = -boost::math::cyl_neumann(1, x, my_pol);
+ }
+ else
+ {
+ y_v = boost::math::cyl_neumann( my_v, x, my_pol);
+ const T y_v_m1 (boost::math::cyl_neumann(T(my_v - 1), x, my_pol));
+ y_v_prime = y_v_m1 - ((my_v * y_v) / x);
+ }
+
+ // Return a tuple containing both Yv(x) and Yv'(x).
+ return boost::math::make_tuple(y_v, y_v_prime);
+ }
+
+ private:
+ const T my_v;
+ const Policy& my_pol;
+ const function_object_yv_and_yv_prime& operator=(const function_object_yv_and_yv_prime&);
+ };
+
+ template<class T> bool my_bisection_unreachable_tolerance(const T&, const T&) { return false; }
+
+ template<class T, class Policy>
+ T initial_guess(const T& v, const int m, const Policy& pol)
+ {
+ BOOST_MATH_STD_USING // ADL of std names, needed for floor.
+
+ // Compute an estimate of the m'th root of cyl_neumann.
+
+ T guess;
+
+ // There is special handling for negative order.
+ if(v < 0)
+ {
+ // Create the positive order and extract its positive floor and ceiling integer parts.
+ const T vv(-v);
+ const T vv_floor(floor(vv));
+
+ // The to-be-found root is bracketed by the roots of the
+ // Bessel function whose reflected, positive integer order
+ // is less than, but nearest to vv.
+
+ // The special case of negative, half-integer order uses
+ // the relation between Yv and spherical Bessel functions
+ // in order to obtain the bracket for the root.
+ // In these special cases, cyl_neumann(-n/2, x) = sph_bessel_j(+n/2, x)
+ // for v = -n/2.
+
+ T root_hi;
+ T root_lo;
+
+ if(m == 1)
+ {
+ // The estimate of the first root for negative order is found using
+ // an adaptive range-searching algorithm.
+ // Take special precautions for the discontinuity at negative,
+ // half-integer orders and use different brackets above and below these.
+ if(T(vv - vv_floor) < 0.5F)
+ {
+ root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
+ }
+ else
+ {
+ root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
+ }
+
+ root_lo = T(root_hi - 0.1F);
+
+ const bool hi_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_hi, pol) < 0);
+
+ while((root_lo > boost::math::tools::epsilon<T>()))
+ {
+ const bool lo_end_of_bracket_is_negative = (boost::math::cyl_neumann(v, root_lo, pol) < 0);
+
+ if(hi_end_of_bracket_is_negative != lo_end_of_bracket_is_negative)
+ {
+ break;
+ }
+
+ root_hi = root_lo;
+
+ // Decrease the lower end of the bracket using an adaptive algorithm.
+ if(root_lo > 0.5F)
+ {
+ root_lo -= 0.5F;
+ }
+ else
+ {
+ root_lo *= 0.75F;
+ }
+ }
+ }
+ else
+ {
+ if(T(vv - vv_floor) < 0.5F)
+ {
+ root_lo = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m - 1, pol);
+ root_hi = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::initial_guess(vv_floor, m, pol);
+ root_lo += 0.01F;
+ root_hi += 0.01F;
+ }
+ else
+ {
+ root_lo = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m - 1, pol);
+ root_hi = boost::math::detail::bessel_zero::cyl_bessel_j_zero_detail::initial_guess(T(vv_floor + 0.5F), m, pol);
+ root_lo += 0.01F;
+ root_hi += 0.01F;
+ }
+ }
+
+ // Perform several steps of bisection iteration to refine the guess.
+ boost::uintmax_t number_of_iterations(12U);
+
+ // Do the bisection iteration.
+ const boost::math::tuple<T, T> guess_pair =
+ boost::math::tools::bisect(
+ boost::math::detail::bessel_zero::cyl_neumann_zero_detail::function_object_yv<T, Policy>(v, pol),
+ root_lo,
+ root_hi,
+ boost::math::detail::bessel_zero::cyl_neumann_zero_detail::my_bisection_unreachable_tolerance<T>,
+ number_of_iterations);
+
+ return (boost::math::get<0>(guess_pair) + boost::math::get<1>(guess_pair)) / 2U;
+ }
+
+ if(m == 1U)
+ {
+ // Get the initial estimate of the first root.
+
+ if(v < 2.2F)
+ {
+ // For small v, use the results of empirical curve fitting.
+ // Mathematica(R) session for the coefficients:
+ // Table[{n, BesselYZero[n, 1]}, {n, 0, 22/10, 1/10}]
+ // N[%, 20]
+ // Fit[%, {n^0, n^1, n^2, n^3, n^4, n^5, n^6}, n]
+ guess = ((((( - T(0.0025095909235652)
+ * v + T(0.021291887049053))
+ * v - T(0.076487785486526))
+ * v + T(0.159110268115362))
+ * v - T(0.241681668765196))
+ * v + T(1.4437846310885244))
+ * v + T(0.89362115190200490);
+ }
+ else
+ {
+ // For larger v, use the second line of Eqs. 10.21.40 in the NIST Handbook.
+ guess = boost::math::detail::bessel_zero::cyl_neumann_zero_detail::equation_nist_10_21_40_b(v);
+ }
+ }
+ else
+ {
+ if(v < 2.2F)
+ {
+ // Use Eq. 10.21.19 in the NIST Handbook.
+ const T a(((v + T(m * 2U)) - T(1.5)) * boost::math::constants::half_pi<T>());
+
+ guess = boost::math::detail::bessel_zero::equation_nist_10_21_19(v, a);
+ }
+ else
+ {
+ // Get an estimate of the m'th root of airy_bi.
+ const T airy_bi_root(boost::math::detail::airy_zero::airy_bi_zero_detail::initial_guess<T>(m));
+
+ // Use Eq. 9.5.26 in the A&S Handbook.
+ guess = boost::math::detail::bessel_zero::equation_as_9_5_26(v, airy_bi_root);
+ }
+ }
+
+ return guess;
+ }
+ } // namespace cyl_neumann_zero_detail
+ } // namespace bessel_zero
+ } } } // namespace boost::math::detail
+
+#endif // _BESSEL_JY_ZERO_2013_01_18_HPP_
diff --git a/boost/math/special_functions/detail/bessel_kn.hpp b/boost/math/special_functions/detail/bessel_kn.hpp
index 5f01460..e3a5023 100644
--- a/boost/math/special_functions/detail/bessel_kn.hpp
+++ b/boost/math/special_functions/detail/bessel_kn.hpp
@@ -22,6 +22,7 @@ namespace boost { namespace math { namespace detail{
template <typename T, typename Policy>
T bessel_kn(int n, T x, const Policy& pol)
{
+ BOOST_MATH_STD_USING
T value, current, prev;
using namespace boost::math::tools;
diff --git a/boost/math/special_functions/detail/bessel_y0.hpp b/boost/math/special_functions/detail/bessel_y0.hpp
index 289bda5..533ab7c 100644
--- a/boost/math/special_functions/detail/bessel_y0.hpp
+++ b/boost/math/special_functions/detail/bessel_y0.hpp
@@ -197,11 +197,22 @@ T bessel_y0(T x, const Policy& pol)
{
T y = 8 / x;
T y2 = y * y;
- T z = x - 0.25f * pi<T>();
rc = evaluate_rational(PC, QC, y2);
rs = evaluate_rational(PS, QS, y2);
- factor = sqrt(2 / (x * pi<T>()));
- value = factor * (rc * sin(z) + y * rs * cos(z));
+ factor = constants::one_div_root_pi<T>() / sqrt(x);
+ //
+ // The following code is really just:
+ //
+ // T z = x - 0.25f * pi<T>();
+ // value = factor * (rc * sin(z) + y * rs * cos(z));
+ //
+ // But using the sin/cos addition formulae and constant values for
+ // sin/cos of PI/4 which then cancel part of the "factor" term as they're all
+ // 1 / sqrt(2):
+ //
+ T sx = sin(x);
+ T cx = cos(x);
+ value = factor * (rc * (sx - cx) + y * rs * (cx + sx));
}
return value;
diff --git a/boost/math/special_functions/detail/bessel_y1.hpp b/boost/math/special_functions/detail/bessel_y1.hpp
index caf09ff..8396f8f 100644
--- a/boost/math/special_functions/detail/bessel_y1.hpp
+++ b/boost/math/special_functions/detail/bessel_y1.hpp
@@ -170,11 +170,21 @@ T bessel_y1(T x, const Policy& pol)
{
T y = 8 / x;
T y2 = y * y;
- T z = x - 0.75f * pi<T>();
rc = evaluate_rational(PC, QC, y2);
rs = evaluate_rational(PS, QS, y2);
- factor = sqrt(2 / (x * pi<T>()));
- value = factor * (rc * sin(z) + y * rs * cos(z));
+ factor = 1 / (sqrt(x) * root_pi<T>());
+ //
+ // This code is really just:
+ //
+ // T z = x - 0.75f * pi<T>();
+ // value = factor * (rc * sin(z) + y * rs * cos(z));
+ //
+ // But using the sin/cos addition rules, plus constants for sin/cos of 3PI/4
+ // which then cancel out with corresponding terms in "factor".
+ //
+ T sx = sin(x);
+ T cx = cos(x);
+ value = factor * (y * rs * (sx - cx) - rc * (sx + cx));
}
return value;
diff --git a/boost/math/special_functions/detail/bessel_yn.hpp b/boost/math/special_functions/detail/bessel_yn.hpp
index b4f9855..0509062 100644
--- a/boost/math/special_functions/detail/bessel_yn.hpp
+++ b/boost/math/special_functions/detail/bessel_yn.hpp
@@ -75,10 +75,11 @@ T bessel_yn(int n, T x, const Policy& pol)
current = bessel_y1(x, pol);
int k = 1;
BOOST_ASSERT(k < n);
+ policies::check_series_iterations<T>("boost::math::bessel_y_n<%1%>(%1%,%1%)", n, pol);
do
{
T fact = 2 * k / x;
- if((tools::max_value<T>() - fabs(prev)) / fact < fabs(current))
+ if((fact > 1) && ((tools::max_value<T>() - fabs(prev)) / fact < fabs(current)))
{
prev /= current;
factor /= current;
diff --git a/boost/math/special_functions/detail/erf_inv.hpp b/boost/math/special_functions/detail/erf_inv.hpp
index d51db9d..77aa72f 100644
--- a/boost/math/special_functions/detail/erf_inv.hpp
+++ b/boost/math/special_functions/detail/erf_inv.hpp
@@ -50,7 +50,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.00538772965071242932965)
};
static const T Q[] = {
- 1,
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, -0.970005043303290640362),
BOOST_MATH_BIG_CONSTANT(T, 64, -1.56574558234175846809),
BOOST_MATH_BIG_CONSTANT(T, 64, 1.56221558398423026363),
@@ -92,7 +92,7 @@ T erf_inv_imp(const T& p, const T& q, const Policy&, const boost::mpl::int_<64>*
BOOST_MATH_BIG_CONSTANT(T, 64, -3.67192254707729348546)
};
static const T Q[] = {
- BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, 6.24264124854247537712),
BOOST_MATH_BIG_CONSTANT(T, 64, 3.9713437953343869095),
BOOST_MATH_BIG_CONSTANT(T, 64, -28.6608180499800029974),
@@ -147,7 +147,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.681149956853776992068e-9)
};
static const T Q[] = {
- 1,
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, 3.46625407242567245975),
BOOST_MATH_BIG_CONSTANT(T, 64, 5.38168345707006855425),
BOOST_MATH_BIG_CONSTANT(T, 64, 4.77846592945843778382),
@@ -176,7 +176,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.266339227425782031962e-11)
};
static const T Q[] = {
- BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, 1.3653349817554063097),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.762059164553623404043),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.220091105764131249824),
@@ -204,7 +204,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.99055709973310326855e-16)
};
static const T Q[] = {
- BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.591429344886417493481),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.138151865749083321638),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.0160746087093676504695),
@@ -231,7 +231,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.116765012397184275695e-17)
};
static const T Q[] = {
- BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.207123112214422517181),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.0169410838120975906478),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.000690538265622684595676),
@@ -258,7 +258,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.348890393399948882918e-21)
};
static const T Q[] = {
- BOOST_MATH_BIG_CONSTANT(T, 64, 1),
+ BOOST_MATH_BIG_CONSTANT(T, 64, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.0845746234001899436914),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.00282092984726264681981),
BOOST_MATH_BIG_CONSTANT(T, 64, 0.468292921940894236786e-4),
@@ -282,7 +282,7 @@ struct erf_roots
BOOST_MATH_STD_USING
T derivative = sign * (2 / sqrt(constants::pi<T>())) * exp(-(guess * guess));
T derivative2 = -2 * guess * derivative;
- return boost::math::make_tuple(((sign > 0) ? boost::math::erf(guess, Policy()) : boost::math::erfc(guess, Policy())) - target, derivative, derivative2);
+ return boost::math::make_tuple(((sign > 0) ? static_cast<T>(boost::math::erf(guess, Policy()) - target) : static_cast<T>(boost::math::erfc(guess, Policy())) - target), derivative, derivative2);
}
erf_roots(T z, int s) : target(z), sign(s) {}
private:
@@ -331,18 +331,34 @@ struct erf_inv_initializer
{
do_init();
}
+ static bool is_value_non_zero(T);
static void do_init()
{
boost::math::erf_inv(static_cast<T>(0.25), Policy());
boost::math::erf_inv(static_cast<T>(0.55), Policy());
boost::math::erf_inv(static_cast<T>(0.95), Policy());
boost::math::erfc_inv(static_cast<T>(1e-15), Policy());
- if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)) != 0)
+ // These following initializations must not be called if
+ // type T can not hold the relevant values without
+ // underflow to zero. We check this at runtime because
+ // some tools such as valgrind silently change the precision
+ // of T at runtime, and numeric_limits basically lies!
+ if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130))))
boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-130)), Policy());
- if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)) != 0)
+
+ // Some compilers choke on constants that would underflow, even in code that isn't instantiated
+ // so try and filter these cases out in the preprocessor:
+#if LDBL_MAX_10_EXP >= 800
+ if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800))))
boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-800)), Policy());
- if(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)) != 0)
+ if(is_value_non_zero(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900))))
boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy());
+#else
+ if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800))))
+ boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-800)), Policy());
+ if(is_value_non_zero(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900))))
+ boost::math::erfc_inv(static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 64, 1e-900)), Policy());
+#endif
}
void force_instantiate()const{}
};
@@ -356,6 +372,15 @@ struct erf_inv_initializer
template <class T, class Policy>
const typename erf_inv_initializer<T, Policy>::init erf_inv_initializer<T, Policy>::initializer;
+template <class T, class Policy>
+bool erf_inv_initializer<T, Policy>::init::is_value_non_zero(T v)
+{
+ // This needs to be non-inline to detect whether v is non zero at runtime
+ // rather than at compile time, only relevant when running under valgrind
+ // which changes long double's to double's on the fly.
+ return v != 0;
+}
+
} // namespace detail
template <class T, class Policy>
@@ -368,7 +393,7 @@ typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol)
//
static const char* function = "boost::math::erfc_inv<%1%>(%1%, %1%)";
if((z < 0) || (z > 2))
- policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);
+ return policies::raise_domain_error<result_type>(function, "Argument outside range [0,2] in inverse erfc function (got p=%1%).", z, pol);
if(z == 0)
return policies::raise_overflow_error<result_type>(function, 0, pol);
if(z == 2)
@@ -432,7 +457,7 @@ typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol)
//
static const char* function = "boost::math::erf_inv<%1%>(%1%, %1%)";
if((z < -1) || (z > 1))
- policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);
+ return policies::raise_domain_error<result_type>(function, "Argument outside range [-1, 1] in inverse erf function (got p=%1%).", z, pol);
if(z == 1)
return policies::raise_overflow_error<result_type>(function, 0, pol);
if(z == -1)
diff --git a/boost/math/special_functions/detail/fp_traits.hpp b/boost/math/special_functions/detail/fp_traits.hpp
index 50c034d..63ebf11 100644
--- a/boost/math/special_functions/detail/fp_traits.hpp
+++ b/boost/math/special_functions/detail/fp_traits.hpp
@@ -351,6 +351,13 @@ struct fp_traits_non_native<long double, extended_double_precision>
// the Intel extended double precision format (80 bits) and
// the IEEE extended double precision format with 15 exponent bits (128 bits).
+#elif defined(__GNUC__) && (LDBL_MANT_DIG == 106)
+
+//
+// Define nothing here and fall though to generic_tag:
+// We have GCC's "double double" in effect, and any attempt
+// to handle it via bit-fiddling is pretty much doomed to fail...
+//
// long double (>64 bits), PowerPC ---------------------------------------------
@@ -546,7 +553,9 @@ struct select_native<long double>
&& !defined(__DECCXX)\
&& !defined(__osf__) \
&& !defined(__SGI_STL_PORT) && !defined(_STLPORT_VERSION)\
- && !defined(BOOST_MATH_DISABLE_STD_FPCLASSIFY)
+ && !defined(__FAST_MATH__)\
+ && !defined(BOOST_MATH_DISABLE_STD_FPCLASSIFY)\
+ && !defined(BOOST_INTEL)
# define BOOST_MATH_USE_STD_FPCLASSIFY
#endif
diff --git a/boost/math/special_functions/detail/gamma_inva.hpp b/boost/math/special_functions/detail/gamma_inva.hpp
index 549bc3d..7c32d29 100644
--- a/boost/math/special_functions/detail/gamma_inva.hpp
+++ b/boost/math/special_functions/detail/gamma_inva.hpp
@@ -75,7 +75,7 @@ T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
//
if(p == 0)
{
- return tools::max_value<T>();
+ return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy());
}
if(q == 0)
{
@@ -144,7 +144,7 @@ T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
//
std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);
if(max_iter >= policies::get_max_root_iterations<Policy>())
- policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
+ return policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
return (r.first + r.second) / 2;
}
@@ -165,7 +165,7 @@ inline typename tools::promote_args<T1, T2>::type
if(p == 0)
{
- return tools::max_value<result_type>();
+ policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy());
}
if(p == 1)
{
@@ -195,7 +195,7 @@ inline typename tools::promote_args<T1, T2>::type
if(q == 1)
{
- return tools::max_value<result_type>();
+ policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", 0, Policy());
}
if(q == 0)
{
diff --git a/boost/math/special_functions/detail/ibeta_inv_ab.hpp b/boost/math/special_functions/detail/ibeta_inv_ab.hpp
index 8318a28..f5735a8 100644
--- a/boost/math/special_functions/detail/ibeta_inv_ab.hpp
+++ b/boost/math/special_functions/detail/ibeta_inv_ab.hpp
@@ -153,7 +153,7 @@ T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab,
boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
if(max_iter >= policies::get_max_root_iterations<Policy>())
- policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
+ return policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
return (r.first + r.second) / 2;
}
@@ -172,9 +172,10 @@ typename tools::promote_args<RT1, RT2, RT3>::type
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ static const char* function = "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)";
if(p == 0)
{
- return tools::max_value<result_type>();
+ return policies::raise_overflow_error<result_type>(function, 0, Policy());
}
if(p == 1)
{
@@ -188,7 +189,7 @@ typename tools::promote_args<RT1, RT2, RT3>::type
static_cast<value_type>(p),
static_cast<value_type>(1 - static_cast<value_type>(p)),
false, pol),
- "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)");
+ function);
}
template <class RT1, class RT2, class RT3, class Policy>
@@ -204,9 +205,10 @@ typename tools::promote_args<RT1, RT2, RT3>::type
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ static const char* function = "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)";
if(q == 1)
{
- return tools::max_value<result_type>();
+ return policies::raise_overflow_error<result_type>(function, 0, Policy());
}
if(q == 0)
{
@@ -220,7 +222,7 @@ typename tools::promote_args<RT1, RT2, RT3>::type
static_cast<value_type>(1 - static_cast<value_type>(q)),
static_cast<value_type>(q),
false, pol),
- "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)");
+ function);
}
template <class RT1, class RT2, class RT3, class Policy>
@@ -236,13 +238,14 @@ typename tools::promote_args<RT1, RT2, RT3>::type
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
+ static const char* function = "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)";
if(p == 0)
{
return tools::min_value<result_type>();
}
if(p == 1)
{
- return tools::max_value<result_type>();
+ return policies::raise_overflow_error<result_type>(function, 0, Policy());
}
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
@@ -252,13 +255,14 @@ typename tools::promote_args<RT1, RT2, RT3>::type
static_cast<value_type>(p),
static_cast<value_type>(1 - static_cast<value_type>(p)),
true, pol),
- "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)");
+ function);
}
template <class RT1, class RT2, class RT3, class Policy>
typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
{
+ static const char* function = "boost::math::ibeta_invb<%1%>(%1%, %1%, %1%)";
typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
typedef typename policies::normalise<
@@ -274,7 +278,7 @@ typename tools::promote_args<RT1, RT2, RT3>::type
}
if(q == 0)
{
- return tools::max_value<result_type>();
+ return policies::raise_overflow_error<result_type>(function, 0, Policy());
}
return policies::checked_narrowing_cast<result_type, forwarding_policy>(
@@ -282,9 +286,9 @@ typename tools::promote_args<RT1, RT2, RT3>::type
static_cast<value_type>(a),
static_cast<value_type>(x),
static_cast<value_type>(1 - static_cast<value_type>(q)),
- static_cast<value_type>(q),
+ static_cast<value_type>(q),
true, pol),
- "boost::math::ibetac_invb<%1%>(%1%,%1%,%1%)");
+ function);
}
template <class RT1, class RT2, class RT3>
diff --git a/boost/math/special_functions/detail/ibeta_inverse.hpp b/boost/math/special_functions/detail/ibeta_inverse.hpp
index ccfa919..a9fe8cd 100644
--- a/boost/math/special_functions/detail/ibeta_inverse.hpp
+++ b/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -35,12 +35,12 @@ struct temme_root_finder
if(y == 0)
{
T big = tools::max_value<T>() / 4;
- return boost::math::make_tuple(-big, -big);
+ return boost::math::make_tuple(static_cast<T>(-big), static_cast<T>(-big));
}
if(x == 0)
{
T big = tools::max_value<T>() / 4;
- return boost::math::make_tuple(-big, big);
+ return boost::math::make_tuple(static_cast<T>(-big), big);
}
T f = log(x) + a * log(y) + t;
T f1 = (1 / x) - (a / (y));
@@ -455,6 +455,11 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
BOOST_MATH_STD_USING // For ADL of math functions.
//
+ // The flag invert is set to true if we swap a for b and p for q,
+ // in which case the result has to be subtracted from 1:
+ //
+ bool invert = false;
+ //
// Handle trivial cases first:
//
if(q == 0)
@@ -467,17 +472,19 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if(py) *py = 1;
return 0;
}
- else if((a == 1) && (b == 1))
+ else if(a == 1)
{
- if(py) *py = 1 - p;
- return p;
+ if(b == 1)
+ {
+ if(py) *py = 1 - p;
+ return p;
+ }
+ // Change things around so we can handle as b == 1 special case below:
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = true;
}
//
- // The flag invert is set to true if we swap a for b and p for q,
- // in which case the result has to be subtracted from 1:
- //
- bool invert = false;
- //
// Depending upon which approximation method we use, we may end up
// calculating either x or y initially (where y = 1-x):
//
@@ -495,21 +502,61 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// Student's T with b = 0.5 gets handled as a special case, swap
// around if the arguments are in the "wrong" order:
//
- if((a == 0.5f) && (b >= 0.5f))
+ if(a == 0.5f)
{
- std::swap(a, b);
- std::swap(p, q);
- invert = !invert;
+ if(b == 0.5f)
+ {
+ x = sin(p * constants::half_pi<T>());
+ x *= x;
+ if(py)
+ {
+ *py = sin(q * constants::half_pi<T>());
+ *py *= *py;
+ }
+ return x;
+ }
+ else if(b > 0.5f)
+ {
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = !invert;
+ }
}
//
// Select calculation method for the initial estimate:
//
- if((b == 0.5f) && (a >= 0.5f))
+ if((b == 0.5f) && (a >= 0.5f) && (p != 1))
{
//
// We have a Student's T distribution:
x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol);
}
+ else if(b == 1)
+ {
+ if(p < q)
+ {
+ if(a > 1)
+ {
+ x = pow(p, 1 / a);
+ y = -boost::math::expm1(log(p) / a, pol);
+ }
+ else
+ {
+ x = pow(p, 1 / a);
+ y = 1 - x;
+ }
+ }
+ else
+ {
+ x = exp(boost::math::log1p(-q, pol) / a);
+ y = -boost::math::expm1(boost::math::log1p(-q, pol) / a, pol);
+ }
+ if(invert)
+ std::swap(x, y);
+ if(py)
+ *py = y;
+ return x;
+ }
else if(a + b > 5)
{
//
@@ -866,14 +913,16 @@ template <class T1, class T2, class T3>
inline typename tools::promote_args<T1, T2, T3>::type
ibeta_inv(T1 a, T2 b, T3 p)
{
- return ibeta_inv(a, b, p, static_cast<T1*>(0), policies::policy<>());
+ typedef typename tools::promote_args<T1, T2, T3>::type result_type;
+ return ibeta_inv(a, b, p, static_cast<result_type*>(0), policies::policy<>());
}
template <class T1, class T2, class T3, class Policy>
inline typename tools::promote_args<T1, T2, T3>::type
ibeta_inv(T1 a, T2 b, T3 p, const Policy& pol)
{
- return ibeta_inv(a, b, p, static_cast<T1*>(0), pol);
+ typedef typename tools::promote_args<T1, T2, T3>::type result_type;
+ return ibeta_inv(a, b, p, static_cast<result_type*>(0), pol);
}
template <class T1, class T2, class T3, class T4, class Policy>
@@ -892,11 +941,11 @@ inline typename tools::promote_args<T1, T2, T3, T4>::type
policies::assert_undefined<> >::type forwarding_policy;
if(a <= 0)
- policies::raise_domain_error<result_type>(function, "The argument a to the incomplete beta function inverse must be greater than zero (got a=%1%).", a, pol);
+ return policies::raise_domain_error<result_type>(function, "The argument a to the incomplete beta function inverse must be greater than zero (got a=%1%).", a, pol);
if(b <= 0)
- policies::raise_domain_error<result_type>(function, "The argument b to the incomplete beta function inverse must be greater than zero (got b=%1%).", b, pol);
+ return policies::raise_domain_error<result_type>(function, "The argument b to the incomplete beta function inverse must be greater than zero (got b=%1%).", b, pol);
if((q < 0) || (q > 1))
- policies::raise_domain_error<result_type>(function, "Argument q outside the range [0,1] in the incomplete beta function inverse (got q=%1%).", q, pol);
+ return policies::raise_domain_error<result_type>(function, "Argument q outside the range [0,1] in the incomplete beta function inverse (got q=%1%).", q, pol);
value_type rx, ry;
@@ -922,16 +971,16 @@ template <class RT1, class RT2, class RT3>
inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_inv(RT1 a, RT2 b, RT3 q)
{
- typedef typename remove_cv<RT1>::type dummy;
- return ibetac_inv(a, b, q, static_cast<dummy*>(0), policies::policy<>());
+ typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
+ return ibetac_inv(a, b, q, static_cast<result_type*>(0), policies::policy<>());
}
template <class RT1, class RT2, class RT3, class Policy>
inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_inv(RT1 a, RT2 b, RT3 q, const Policy& pol)
{
- typedef typename remove_cv<RT1>::type dummy;
- return ibetac_inv(a, b, q, static_cast<dummy*>(0), pol);
+ typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
+ return ibetac_inv(a, b, q, static_cast<result_type*>(0), pol);
}
} // namespace math
diff --git a/boost/math/special_functions/detail/igamma_inverse.hpp b/boost/math/special_functions/detail/igamma_inverse.hpp
index 53875ff..fd0189c 100644
--- a/boost/math/special_functions/detail/igamma_inverse.hpp
+++ b/boost/math/special_functions/detail/igamma_inverse.hpp
@@ -281,11 +281,11 @@ T find_inverse_gamma(T a, T p, T q, const Policy& pol, bool* p_has_10_digits)
// DiDonato and Morris Eq 35:
T v = log(p) + boost::math::lgamma(ap1, pol);
z = exp((v + w) / a);
- s = boost::math::log1p(z / ap1 * (1 + z / ap2));
+ s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
z = exp((v + z - s) / a);
- s = boost::math::log1p(z / ap1 * (1 + z / ap2));
+ s = boost::math::log1p(z / ap1 * (1 + z / ap2), pol);
z = exp((v + z - s) / a);
- s = boost::math::log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))));
+ s = boost::math::log1p(z / ap1 * (1 + z / ap2 * (1 + z / (a + 3))), pol);
z = exp((v + z - s) / a);
BOOST_MATH_INSTRUMENT_VARIABLE(z);
}
@@ -341,7 +341,7 @@ struct gamma_p_inverse_func
// flag is set, then Q(x) - q and it's derivatives.
//
typedef typename policies::evaluation<T, Policy>::type value_type;
- typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
+ // typedef typename lanczos::lanczos<T, Policy>::type evaluation_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
@@ -378,7 +378,7 @@ struct gamma_p_inverse_func
f2 = -f2;
}
- return boost::math::make_tuple(f - p, f1, f2);
+ return boost::math::make_tuple(static_cast<T>(f - p), f1, f2);
}
private:
T a, p;
@@ -396,11 +396,11 @@ T gamma_p_inv_imp(T a, T p, const Policy& pol)
BOOST_MATH_INSTRUMENT_VARIABLE(p);
if(a <= 0)
- policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
+ return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
if((p < 0) || (p > 1))
- policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
+ return policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got p=%1%).", p, pol);
if(p == 1)
- return tools::max_value<T>();
+ return policies::raise_overflow_error<T>(function, 0, Policy());
if(p == 0)
return 0;
bool has_10_digits;
@@ -456,11 +456,11 @@ T gamma_q_inv_imp(T a, T q, const Policy& pol)
static const char* function = "boost::math::gamma_q_inv<%1%>(%1%, %1%)";
if(a <= 0)
- policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
+ return policies::raise_domain_error<T>(function, "Argument a in the incomplete gamma function inverse must be >= 0 (got a=%1%).", a, pol);
if((q < 0) || (q > 1))
- policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
+ return policies::raise_domain_error<T>(function, "Probabilty must be in the range [0,1] in the incomplete gamma function inverse (got q=%1%).", q, pol);
if(q == 0)
- return tools::max_value<T>();
+ return policies::raise_overflow_error<T>(function, 0, Policy());
if(q == 1)
return 0;
bool has_10_digits;
diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp
index f8846bf..edef3a0 100644
--- a/boost/math/special_functions/detail/lanczos_sse2.hpp
+++ b/boost/math/special_functions/detail/lanczos_sse2.hpp
@@ -51,11 +51,11 @@ inline double lanczos13m53::lanczos_sum<double>(const double& x)
static_cast<double>(23531376880.41075968857200767445163675473L),
static_cast<double>(0u)
};
- register __m128d vx = _mm_load1_pd(&x);
- register __m128d sum_even = _mm_load_pd(coeff);
- register __m128d sum_odd = _mm_load_pd(coeff+2);
- register __m128d nc_odd, nc_even;
- register __m128d vx2 = _mm_mul_pd(vx, vx);
+ __m128d vx = _mm_load1_pd(&x);
+ __m128d sum_even = _mm_load_pd(coeff);
+ __m128d sum_odd = _mm_load_pd(coeff+2);
+ __m128d nc_odd, nc_even;
+ __m128d vx2 = _mm_mul_pd(vx, vx);
sum_even = _mm_mul_pd(sum_even, vx2);
nc_even = _mm_load_pd(coeff + 4);
@@ -136,11 +136,11 @@ inline double lanczos13m53::lanczos_sum_expG_scaled<double>(const double& x)
static_cast<double>(56906521.91347156388090791033559122686859L),
static_cast<double>(0u)
};
- register __m128d vx = _mm_load1_pd(&x);
- register __m128d sum_even = _mm_load_pd(coeff);
- register __m128d sum_odd = _mm_load_pd(coeff+2);
- register __m128d nc_odd, nc_even;
- register __m128d vx2 = _mm_mul_pd(vx, vx);
+ __m128d vx = _mm_load1_pd(&x);
+ __m128d sum_even = _mm_load_pd(coeff);
+ __m128d sum_odd = _mm_load_pd(coeff+2);
+ __m128d nc_odd, nc_even;
+ __m128d vx2 = _mm_mul_pd(vx, vx);
sum_even = _mm_mul_pd(sum_even, vx2);
nc_even = _mm_load_pd(coeff + 4);
diff --git a/boost/math/special_functions/detail/lgamma_small.hpp b/boost/math/special_functions/detail/lgamma_small.hpp
index ec28ed2..e65f8b7 100644
--- a/boost/math/special_functions/detail/lgamma_small.hpp
+++ b/boost/math/special_functions/detail/lgamma_small.hpp
@@ -87,7 +87,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l *
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.324588649825948492091e-4))
};
static const T Q[] = {
- static_cast<T>(0.1e1),
+ static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.196202987197795200688e1)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.148019669424231326694e1)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.541391432071720958364e0)),
@@ -198,7 +198,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<64>&, const Policy& /* l *
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.431171342679297331241e-3))
};
static const T Q[] = {
- static_cast<T>(0.1e1),
+ static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.1e1)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.150169356054485044494e1)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 0.846973248876495016101e0)),
static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -0.220095151814995745555e0)),
@@ -278,7 +278,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l
BOOST_MATH_BIG_CONSTANT(T, 113, -0.70529798686542184668416911331718963364e-8)
};
static const T Q[] = {
- BOOST_MATH_BIG_CONSTANT(T, 113, 1),
+ BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 113, 2.5877485070422317542808137697939233685),
BOOST_MATH_BIG_CONSTANT(T, 113, 2.8797959228352591788629602533153837126),
BOOST_MATH_BIG_CONSTANT(T, 113, 1.8030885955284082026405495275461180977),
@@ -357,7 +357,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l
BOOST_MATH_BIG_CONSTANT(T, 113, 0.13680157145361387405588201461036338274e-8)
};
static const T Q[] = {
- 1,
+ BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 113, 4.9106336261005990534095838574132225599),
BOOST_MATH_BIG_CONSTANT(T, 113, 10.258804800866438510889341082793078432),
BOOST_MATH_BIG_CONSTANT(T, 113, 11.88588976846826108836629960537466889),
@@ -408,7 +408,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l
BOOST_MATH_BIG_CONSTANT(T, 113, 0.8207548771933585614380644961342925976e-6)
};
static const T Q[] = {
- 1,
+ BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 113, -2.9629552288944259229543137757200262073),
BOOST_MATH_BIG_CONSTANT(T, 113, 3.7118380799042118987185957298964772755),
BOOST_MATH_BIG_CONSTANT(T, 113, -2.5569815272165399297600586376727357187),
@@ -449,7 +449,7 @@ T lgamma_small_imp(T z, T zm1, T zm2, const mpl::int_<113>&, const Policy& /* l
BOOST_MATH_BIG_CONSTANT(T, 113, 0.13240510580220763969511741896361984162e-6)
};
static const T Q[] = {
- 1,
+ BOOST_MATH_BIG_CONSTANT(T, 113, 1.0),
BOOST_MATH_BIG_CONSTANT(T, 113, -2.4240003754444040525462170802796471996),
BOOST_MATH_BIG_CONSTANT(T, 113, 2.4868383476933178722203278602342786002),
BOOST_MATH_BIG_CONSTANT(T, 113, -1.4047068395206343375520721509193698547),
diff --git a/boost/math/special_functions/detail/round_fwd.hpp b/boost/math/special_functions/detail/round_fwd.hpp
index 952259a..8c45a7d 100644
--- a/boost/math/special_functions/detail/round_fwd.hpp
+++ b/boost/math/special_functions/detail/round_fwd.hpp
@@ -9,6 +9,7 @@
#define BOOST_MATH_SPECIAL_ROUND_FWD_HPP
#include <boost/config.hpp>
+#include <boost/math/tools/promotion.hpp>
#ifdef _MSC_VER
#pragma once
@@ -20,9 +21,9 @@ namespace boost
{
template <class T, class Policy>
- T trunc(const T& v, const Policy& pol);
+ typename tools::promote_args<T>::type trunc(const T& v, const Policy& pol);
template <class T>
- T trunc(const T& v);
+ typename tools::promote_args<T>::type trunc(const T& v);
template <class T, class Policy>
int itrunc(const T& v, const Policy& pol);
template <class T>
@@ -38,9 +39,9 @@ namespace boost
boost::long_long_type lltrunc(const T& v);
#endif
template <class T, class Policy>
- T round(const T& v, const Policy& pol);
+ typename tools::promote_args<T>::type round(const T& v, const Policy& pol);
template <class T>
- T round(const T& v);
+ typename tools::promote_args<T>::type round(const T& v);
template <class T, class Policy>
int iround(const T& v, const Policy& pol);
template <class T>
@@ -76,5 +77,17 @@ namespace boost
}
}
+
+#undef BOOST_MATH_STD_USING
+#define BOOST_MATH_STD_USING BOOST_MATH_STD_USING_CORE\
+ using boost::math::round;\
+ using boost::math::iround;\
+ using boost::math::lround;\
+ using boost::math::trunc;\
+ using boost::math::itrunc;\
+ using boost::math::ltrunc;\
+ using boost::math::modf;
+
+
#endif // BOOST_MATH_SPECIAL_ROUND_FWD_HPP
diff --git a/boost/math/special_functions/detail/t_distribution_inv.hpp b/boost/math/special_functions/detail/t_distribution_inv.hpp
index 4e0d2d1..72f6f0c 100644
--- a/boost/math/special_functions/detail/t_distribution_inv.hpp
+++ b/boost/math/special_functions/detail/t_distribution_inv.hpp
@@ -372,7 +372,13 @@ T inverse_students_t(T df, T u, T v, const Policy& pol, bool* pexact = 0)
else
{
calculate_real:
- if(df < 3)
+ if(df > 0x10000000)
+ {
+ result = -boost::math::erfc_inv(2 * u, pol) * constants::root_two<T>();
+ if((pexact) && (df >= 1e20))
+ *pexact = true;
+ }
+ else if(df < 3)
{
//
// Use a roughly linear scheme to choose between Shaw's
@@ -395,7 +401,7 @@ calculate_real:
// where we use Shaw's tail series.
// The crossover point is roughly exponential in -df:
//
- T crossover = ldexp(1.0f, iround(T(df / -0.654f), pol));
+ T crossover = ldexp(1.0f, iround(T(df / -0.654f), typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type()));
if(u > crossover)
{
result = boost::math::detail::inverse_students_t_hill(df, u, pol);
@@ -410,15 +416,14 @@ calculate_real:
}
template <class T, class Policy>
-inline T find_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const Policy& pol)
+inline T find_ibeta_inv_from_t_dist(T a, T p, T /*q*/, T* py, const Policy& pol)
{
- T u = (p > q) ? T(0.5f - q) / T(2) : T(p / 2);
- T v = 1 - u; // u < 0.5 so no cancellation error
+ T u = p / 2;
+ T v = 1 - u;
T df = a * 2;
T t = boost::math::detail::inverse_students_t(df, u, v, pol);
- T x = df / (df + t * t);
*py = t * t / (df + t * t);
- return x;
+ return df / (df + t * t);
}
template <class T, class Policy>
diff --git a/boost/math/special_functions/detail/unchecked_bernoulli.hpp b/boost/math/special_functions/detail/unchecked_bernoulli.hpp
new file mode 100644
index 0000000..03c3766
--- /dev/null
+++ b/boost/math/special_functions/detail/unchecked_bernoulli.hpp
@@ -0,0 +1,700 @@
+
+///////////////////////////////////////////////////////////////////////////////
+// Copyright 2013 Nikhar Agrawal
+// Copyright 2013 Christopher Kormanyos
+// Copyright 2013 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_MATH_UNCHECKED_BERNOULLI_HPP
+#define BOOST_MATH_UNCHECKED_BERNOULLI_HPP
+
+#include <limits>
+#include <cmath>
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/constants/constants.hpp>
+#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/mpl/int.hpp>
+#include <boost/type_traits/is_convertible.hpp>
+
+namespace boost { namespace math {
+
+namespace detail {
+
+template <unsigned N>
+struct max_bernoulli_index
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 17);
+};
+
+template <>
+struct max_bernoulli_index<1>
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 32);
+};
+
+template <>
+struct max_bernoulli_index<2>
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 129);
+};
+
+template <>
+struct max_bernoulli_index<3>
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 1156);
+};
+
+template <>
+struct max_bernoulli_index<4>
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 11);
+};
+
+template <class T>
+struct bernoulli_imp_variant
+{
+ static const unsigned value =
+ (std::numeric_limits<T>::max_exponent == 128)
+ && (std::numeric_limits<T>::radix == 2)
+ && (std::numeric_limits<T>::digits <= std::numeric_limits<float>::digits)
+ && (boost::is_convertible<float, T>::value) ? 1 :
+ (
+ (std::numeric_limits<T>::max_exponent == 1024)
+ && (std::numeric_limits<T>::radix == 2)
+ && (std::numeric_limits<T>::digits <= std::numeric_limits<double>::digits)
+ && (boost::is_convertible<double, T>::value) ? 2 :
+ (
+ (std::numeric_limits<T>::max_exponent == 16384)
+ && (std::numeric_limits<T>::radix == 2)
+ && (std::numeric_limits<T>::digits <= std::numeric_limits<long double>::digits)
+ && (boost::is_convertible<long double, T>::value) ? 3 : (!is_convertible<boost::int64_t, T>::value ? 4 : 0)
+ )
+ );
+};
+
+} // namespace detail
+
+template <class T>
+struct max_bernoulli_b2n : public detail::max_bernoulli_index<detail::bernoulli_imp_variant<T>::value>{};
+
+namespace detail{
+
+template <class T>
+inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<0>& )
+{
+ static const boost::array<boost::int64_t, 1 + max_bernoulli_b2n<T>::value> numerators =
+ {{
+ boost::int64_t( +1LL),
+ boost::int64_t( +1LL),
+ boost::int64_t( -1LL),
+ boost::int64_t( +1LL),
+ boost::int64_t( -1LL),
+ boost::int64_t( +5LL),
+ boost::int64_t( -691LL),
+ boost::int64_t( +7LL),
+ boost::int64_t( -3617LL),
+ boost::int64_t( +43867LL),
+ boost::int64_t( -174611LL),
+ boost::int64_t( +854513LL),
+ boost::int64_t( -236364091LL),
+ boost::int64_t( +8553103LL),
+ boost::int64_t( -23749461029LL),
+ boost::int64_t(+8615841276005LL),
+ boost::int64_t(-7709321041217LL),
+ boost::int64_t(+2577687858367LL)
+ }};
+
+ static const boost::array<boost::int64_t, 1 + max_bernoulli_b2n<T>::value> denominators =
+ {{
+ boost::int64_t( 1LL),
+ boost::int64_t( 6LL),
+ boost::int64_t( 30LL),
+ boost::int64_t( 42LL),
+ boost::int64_t( 30LL),
+ boost::int64_t( 66LL),
+ boost::int64_t( 2730LL),
+ boost::int64_t( 6LL),
+ boost::int64_t( 510LL),
+ boost::int64_t( 798LL),
+ boost::int64_t( 330LL),
+ boost::int64_t( 138LL),
+ boost::int64_t( 2730LL),
+ boost::int64_t( 6LL),
+ boost::int64_t( 870LL),
+ boost::int64_t( 14322LL),
+ boost::int64_t( 510LL),
+ boost::int64_t( 6LL)
+ }};
+ return T(numerators[n]) / denominators[n];
+}
+
+template <class T>
+inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<1>& )
+{
+ static const boost::array<float, 1 + max_bernoulli_b2n<T>::value> bernoulli_data =
+ {{
+ +1.00000000000000000000000000000000000000000F,
+ +0.166666666666666666666666666666666666666667F,
+ -0.0333333333333333333333333333333333333333333F,
+ +0.0238095238095238095238095238095238095238095F,
+ -0.0333333333333333333333333333333333333333333F,
+ +0.0757575757575757575757575757575757575757576F,
+ -0.253113553113553113553113553113553113553114F,
+ +1.16666666666666666666666666666666666666667F,
+ -7.09215686274509803921568627450980392156863F,
+ +54.9711779448621553884711779448621553884712F,
+ -529.124242424242424242424242424242424242424F,
+ +6192.12318840579710144927536231884057971014F,
+ -86580.2531135531135531135531135531135531136F,
+ +1.42551716666666666666666666666666666666667e6F,
+ -2.72982310678160919540229885057471264367816e7F,
+ +6.01580873900642368384303868174835916771401e8F,
+ -1.51163157670921568627450980392156862745098e10F,
+ +4.29614643061166666666666666666666666666667e11F,
+ -1.37116552050883327721590879485616327721591e13F,
+ +4.88332318973593166666666666666666666666667e14F,
+ -1.92965793419400681486326681448632668144863e16F,
+ +8.41693047573682615000553709856035437430786e17F,
+ -4.03380718540594554130768115942028985507246e19F,
+ +2.11507486380819916056014539007092198581560e21F,
+ -1.20866265222965259346027311937082525317819e23F,
+ +7.50086674607696436685572007575757575757576e24F,
+ -5.03877810148106891413789303052201257861635e26F,
+ +3.65287764848181233351104308429711779448622e28F,
+ -2.84987693024508822262691464329106781609195e30F,
+ +2.38654274996836276446459819192192149717514e32F,
+ -2.13999492572253336658107447651910973926742e34F,
+ +2.05009757234780975699217330956723102516667e36F,
+ -2.09380059113463784090951852900279701847092e38F,
+ }};
+
+ return bernoulli_data[n];
+}
+
+
+template <class T>
+inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<2>& )
+{
+ static const boost::array<double, 1 + max_bernoulli_b2n<T>::value> bernoulli_data =
+ {{
+ +1.00000000000000000000000000000000000000000,
+ +0.166666666666666666666666666666666666666667,
+ -0.0333333333333333333333333333333333333333333,
+ +0.0238095238095238095238095238095238095238095,
+ -0.0333333333333333333333333333333333333333333,
+ +0.0757575757575757575757575757575757575757576,
+ -0.253113553113553113553113553113553113553114,
+ +1.16666666666666666666666666666666666666667,
+ -7.09215686274509803921568627450980392156863,
+ +54.9711779448621553884711779448621553884712,
+ -529.124242424242424242424242424242424242424,
+ +6192.12318840579710144927536231884057971014,
+ -86580.2531135531135531135531135531135531136,
+ +1.42551716666666666666666666666666666666667e6,
+ -2.72982310678160919540229885057471264367816e7,
+ +6.01580873900642368384303868174835916771401e8,
+ -1.51163157670921568627450980392156862745098e10,
+ +4.29614643061166666666666666666666666666667e11,
+ -1.37116552050883327721590879485616327721591e13,
+ +4.88332318973593166666666666666666666666667e14,
+ -1.92965793419400681486326681448632668144863e16,
+ +8.41693047573682615000553709856035437430786e17,
+ -4.03380718540594554130768115942028985507246e19,
+ +2.11507486380819916056014539007092198581560e21,
+ -1.20866265222965259346027311937082525317819e23,
+ +7.50086674607696436685572007575757575757576e24,
+ -5.03877810148106891413789303052201257861635e26,
+ +3.65287764848181233351104308429711779448622e28,
+ -2.84987693024508822262691464329106781609195e30,
+ +2.38654274996836276446459819192192149717514e32,
+ -2.13999492572253336658107447651910973926742e34,
+ +2.05009757234780975699217330956723102516667e36,
+ -2.09380059113463784090951852900279701847092e38,
+ +2.27526964884635155596492603527692645814700e40,
+ -2.62577102862395760473030497361582020814490e42,
+ +3.21250821027180325182047923042649852435219e44,
+ -4.15982781667947109139170744952623589366896e46,
+ +5.69206954820352800238834562191210586444805e48,
+ -8.21836294197845756922906534686173330145509e50,
+ +1.25029043271669930167323398297028955241772e53,
+ -2.00155832332483702749253291988132987687242e55,
+ +3.36749829153643742333966769033387530162196e57,
+ -5.94709705031354477186604968440515408405791e59,
+ +1.10119103236279775595641307904376916046305e62,
+ -2.13552595452535011886583850190410656789733e64,
+ +4.33288969866411924196166130593792062184514e66,
+ -9.18855282416693282262005552155018971389604e68,
+ +2.03468967763290744934550279902200200659751e71,
+ -4.70038339580357310785752555350060606545967e73,
+ +1.13180434454842492706751862577339342678904e76,
+ -2.83822495706937069592641563364817647382847e78,
+ +7.40642489796788506297508271409209841768797e80,
+ -2.00964548027566044834656196727153631868673e83,
+ +5.66571700508059414457193460305193569614195e85,
+ -1.65845111541362169158237133743199123014950e88,
+ +5.03688599504923774192894219151801548124424e90,
+ -1.58614682376581863693634015729664387827410e93,
+ +5.17567436175456269840732406825071225612408e95,
+ -1.74889218402171173396900258776181591451415e98,
+ +6.11605199949521852558245252642641677807677e100,
+ -2.21227769127078349422883234567129324455732e103,
+ +8.27227767987709698542210624599845957312047e105,
+ -3.19589251114157095835916343691808148735263e108,
+ +1.27500822233877929823100243029266798669572e111,
+ -5.25009230867741338994028246245651754469199e113,
+ +2.23018178942416252098692981988387281437383e116,
+ -9.76845219309552044386335133989802393011669e118,
+ +4.40983619784529542722726228748131691918758e121,
+ -2.05085708864640888397293377275830154864566e124,
+ +9.82144332797912771075729696020975210414919e126,
+ -4.84126007982088805087891967099634127611305e129,
+ +2.45530888014809826097834674040886903996737e132,
+ -1.28069268040847475487825132786017857218118e135,
+ +6.86761671046685811921018885984644004360924e137,
+ -3.78464685819691046949789954163795568144895e140,
+ +2.14261012506652915508713231351482720966602e143,
+ -1.24567271371836950070196429616376072194583e146,
+ +7.43457875510001525436796683940520613117807e148,
+ -4.55357953046417048940633332233212748767721e151,
+ +2.86121128168588683453638472510172325229190e154,
+ -1.84377235520338697276882026536287854875414e157,
+ +1.21811545362210466995013165065995213558174e160,
+ -8.24821871853141215484818457296893447301419e162,
+ +5.72258779378329433296516498142978615918685e165,
+ -4.06685305250591047267679693831158655602196e168,
+ +2.95960920646420500628752695815851870426379e171,
+ -2.20495225651894575090311752273445984836379e174,
+ +1.68125970728895998058311525151360665754464e177,
+ -1.31167362135569576486452806355817153004431e180,
+ +1.04678940094780380821832853929823089643829e183,
+ -8.54328935788337077185982546299082774593270e185,
+ +7.12878213224865423522884066771438224721245e188,
+ -6.08029314555358993000847118686477458461988e191,
+ +5.29967764248499239300942910043247266228490e194,
+ -4.71942591687458626443646229013379911103761e197,
+ +4.29284137914029810894168296541074669045521e200,
+ -3.98767449682322074434477655542938795106651e203,
+ +3.78197804193588827138944181161393327898220e206,
+ -3.66142336836811912436858082151197348755196e209,
+ +3.61760902723728623488554609298914089477541e212,
+ -3.64707726451913543621383088655499449048682e215,
+ +3.75087554364544090983452410104814189306842e218,
+ -3.93458672964390282694891288533713429355657e221,
+ +4.20882111481900820046571171111494898242731e224,
+ -4.59022962206179186559802940573325591059371e227,
+ +5.10317257726295759279198185106496768539760e230,
+ -5.78227623036569554015377271242917142512200e233,
+ +6.67624821678358810322637794412809363451080e236,
+ -7.85353076444504163225916259639312444428230e239,
+ +9.41068940670587255245443288258762485293948e242,
+ -1.14849338734651839938498599206805592548354e246,
+ +1.42729587428487856771416320087122499897180e249,
+ -1.80595595869093090142285728117654560926719e252,
+ +2.32615353076608052161297985184708876161736e255,
+ -3.04957517154995947681942819261542593785327e258,
+ +4.06858060764339734424012124124937318633684e261,
+ -5.52310313219743616252320044093186392324280e264,
+ +7.62772793964343924869949690204961215533859e267,
+ -1.07155711196978863132793524001065396932667e271,
+ +1.53102008959691884453440916153355334355847e274,
+ -2.22448916821798346676602348865048510824835e277,
+ +3.28626791906901391668189736436895275365183e280,
+ -4.93559289559603449020711938191575963496999e283,
+ +7.53495712008325067212266049779283956727824e286,
+ -1.16914851545841777278088924731655041783900e290,
+ +1.84352614678389394126646201597702232396492e293,
+ -2.95368261729680829728014917350525183485207e296,
+ +4.80793212775015697668878704043264072227967e299,
+ -7.95021250458852528538243631671158693036798e302,
+ +1.33527841873546338750122832017820518292039e306
+ }};
+
+ return bernoulli_data[n];
+}
+
+template <class T>
+inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<3>& )
+{
+ static const boost::array<long double, 1 + max_bernoulli_b2n<T>::value> bernoulli_data =
+ {{
+ +1.00000000000000000000000000000000000000000L,
+ +0.166666666666666666666666666666666666666667L,
+ -0.0333333333333333333333333333333333333333333L,
+ +0.0238095238095238095238095238095238095238095L,
+ -0.0333333333333333333333333333333333333333333L,
+ +0.0757575757575757575757575757575757575757576L,
+ -0.253113553113553113553113553113553113553114L,
+ +1.16666666666666666666666666666666666666667L,
+ -7.09215686274509803921568627450980392156863L,
+ +54.9711779448621553884711779448621553884712L,
+ -529.124242424242424242424242424242424242424L,
+ +6192.12318840579710144927536231884057971014L,
+ -86580.2531135531135531135531135531135531136L,
+ +1.42551716666666666666666666666666666666667E6L,
+ -2.72982310678160919540229885057471264367816E7L,
+ +6.01580873900642368384303868174835916771401E8L,
+ -1.51163157670921568627450980392156862745098E10L,
+ +4.29614643061166666666666666666666666666667E11L,
+ -1.37116552050883327721590879485616327721591E13L,
+ +4.88332318973593166666666666666666666666667E14L,
+ -1.92965793419400681486326681448632668144863E16L,
+ +8.41693047573682615000553709856035437430786E17L,
+ -4.03380718540594554130768115942028985507246E19L,
+ +2.11507486380819916056014539007092198581560E21L,
+ -1.20866265222965259346027311937082525317819E23L,
+ +7.50086674607696436685572007575757575757576E24L,
+ -5.03877810148106891413789303052201257861635E26L,
+ +3.65287764848181233351104308429711779448622E28L,
+ -2.84987693024508822262691464329106781609195E30L,
+ +2.38654274996836276446459819192192149717514E32L,
+ -2.13999492572253336658107447651910973926742E34L,
+ +2.05009757234780975699217330956723102516667E36L,
+ -2.09380059113463784090951852900279701847092E38L,
+ +2.27526964884635155596492603527692645814700E40L,
+ -2.62577102862395760473030497361582020814490E42L,
+ +3.21250821027180325182047923042649852435219E44L,
+ -4.15982781667947109139170744952623589366896E46L,
+ +5.69206954820352800238834562191210586444805E48L,
+ -8.21836294197845756922906534686173330145509E50L,
+ +1.25029043271669930167323398297028955241772E53L,
+ -2.00155832332483702749253291988132987687242E55L,
+ +3.36749829153643742333966769033387530162196E57L,
+ -5.94709705031354477186604968440515408405791E59L,
+ +1.10119103236279775595641307904376916046305E62L,
+ -2.13552595452535011886583850190410656789733E64L,
+ +4.33288969866411924196166130593792062184514E66L,
+ -9.18855282416693282262005552155018971389604E68L,
+ +2.03468967763290744934550279902200200659751E71L,
+ -4.70038339580357310785752555350060606545967E73L,
+ +1.13180434454842492706751862577339342678904E76L,
+ -2.83822495706937069592641563364817647382847E78L,
+ +7.40642489796788506297508271409209841768797E80L,
+ -2.00964548027566044834656196727153631868673E83L,
+ +5.66571700508059414457193460305193569614195E85L,
+ -1.65845111541362169158237133743199123014950E88L,
+ +5.03688599504923774192894219151801548124424E90L,
+ -1.58614682376581863693634015729664387827410E93L,
+ +5.17567436175456269840732406825071225612408E95L,
+ -1.74889218402171173396900258776181591451415E98L,
+ +6.11605199949521852558245252642641677807677E100L,
+ -2.21227769127078349422883234567129324455732E103L,
+ +8.27227767987709698542210624599845957312047E105L,
+ -3.19589251114157095835916343691808148735263E108L,
+ +1.27500822233877929823100243029266798669572E111L,
+ -5.25009230867741338994028246245651754469199E113L,
+ +2.23018178942416252098692981988387281437383E116L,
+ -9.76845219309552044386335133989802393011669E118L,
+ +4.40983619784529542722726228748131691918758E121L,
+ -2.05085708864640888397293377275830154864566E124L,
+ +9.82144332797912771075729696020975210414919E126L,
+ -4.84126007982088805087891967099634127611305E129L,
+ +2.45530888014809826097834674040886903996737E132L,
+ -1.28069268040847475487825132786017857218118E135L,
+ +6.86761671046685811921018885984644004360924E137L,
+ -3.78464685819691046949789954163795568144895E140L,
+ +2.14261012506652915508713231351482720966602E143L,
+ -1.24567271371836950070196429616376072194583E146L,
+ +7.43457875510001525436796683940520613117807E148L,
+ -4.55357953046417048940633332233212748767721E151L,
+ +2.86121128168588683453638472510172325229190E154L,
+ -1.84377235520338697276882026536287854875414E157L,
+ +1.21811545362210466995013165065995213558174E160L,
+ -8.24821871853141215484818457296893447301419E162L,
+ +5.72258779378329433296516498142978615918685E165L,
+ -4.06685305250591047267679693831158655602196E168L,
+ +2.95960920646420500628752695815851870426379E171L,
+ -2.20495225651894575090311752273445984836379E174L,
+ +1.68125970728895998058311525151360665754464E177L,
+ -1.31167362135569576486452806355817153004431E180L,
+ +1.04678940094780380821832853929823089643829E183L,
+ -8.54328935788337077185982546299082774593270E185L,
+ +7.12878213224865423522884066771438224721245E188L,
+ -6.08029314555358993000847118686477458461988E191L,
+ +5.29967764248499239300942910043247266228490E194L,
+ -4.71942591687458626443646229013379911103761E197L,
+ +4.29284137914029810894168296541074669045521E200L,
+ -3.98767449682322074434477655542938795106651E203L,
+ +3.78197804193588827138944181161393327898220E206L,
+ -3.66142336836811912436858082151197348755196E209L,
+ +3.61760902723728623488554609298914089477541E212L,
+ -3.64707726451913543621383088655499449048682E215L,
+ +3.75087554364544090983452410104814189306842E218L,
+ -3.93458672964390282694891288533713429355657E221L,
+ +4.20882111481900820046571171111494898242731E224L,
+ -4.59022962206179186559802940573325591059371E227L,
+ +5.10317257726295759279198185106496768539760E230L,
+ -5.78227623036569554015377271242917142512200E233L,
+ +6.67624821678358810322637794412809363451080E236L,
+ -7.85353076444504163225916259639312444428230E239L,
+ +9.41068940670587255245443288258762485293948E242L,
+ -1.14849338734651839938498599206805592548354E246L,
+ +1.42729587428487856771416320087122499897180E249L,
+ -1.80595595869093090142285728117654560926719E252L,
+ +2.32615353076608052161297985184708876161736E255L,
+ -3.04957517154995947681942819261542593785327E258L,
+ +4.06858060764339734424012124124937318633684E261L,
+ -5.52310313219743616252320044093186392324280E264L,
+ +7.62772793964343924869949690204961215533859E267L,
+ -1.07155711196978863132793524001065396932667E271L,
+ +1.53102008959691884453440916153355334355847E274L,
+ -2.22448916821798346676602348865048510824835E277L,
+ +3.28626791906901391668189736436895275365183E280L,
+ -4.93559289559603449020711938191575963496999E283L,
+ +7.53495712008325067212266049779283956727824E286L,
+ -1.16914851545841777278088924731655041783900E290L,
+ +1.84352614678389394126646201597702232396492E293L,
+ -2.95368261729680829728014917350525183485207E296L,
+ +4.80793212775015697668878704043264072227967E299L,
+ -7.95021250458852528538243631671158693036798E302L,
+ +1.33527841873546338750122832017820518292039E306L,
+#if LDBL_MAX_EXP == 16384
+ // Entries 260 - 600 http://www.wolframalpha.com/input/?i=TABLE[N[Bernoulli[i]%2C40]%2C+{i%2C258%2C600%2C2}]
+ -2.277640649601959593875058983506938037019e309L,
+ 3.945184036046326234163525556422667595884e312L,
+ -6.938525772130602106071724989641405550473e315L,
+ 1.238896367577564823729057820219210929986e319L,
+ -2.245542599169309759499987966025604480745e322L,
+ 4.131213176073842359732511639489669404266e325L,
+ -7.713581346815269584960928069762882771369e328L,
+ 1.461536066837669600638613788471335541313e332L,
+ -2.809904606225532896862935642992712059631e335L,
+ 5.480957121318876639512096994413992284327e338L,
+ -1.084573284087686110518125291186079616320e342L,
+ 2.176980775647663539729165173863716459962e345L,
+ -4.431998786117553751947439433256752608068e348L,
+ 9.150625657715535047417756278073770096073e351L,
+ -1.915867353003157351316577579148683133613e355L,
+ 4.067256303542212258698836003682016040629e358L,
+ -8.754223791037736616228150209910348734629e361L,
+ 1.910173688735533667244373747124109379826e365L,
+ -4.225001320265091714631115064713174404607e368L,
+ 9.471959352547827678466770796787503034505e371L,
+ -2.152149973279986829719817376756088198573e375L,
+ 4.955485775334221051344839716507812871361e378L,
+ -1.156225941759134696630956889716381968142e382L,
+ 2.733406597646137698610991926705098514017e385L,
+ -6.546868135325176947099912523279938546333e388L,
+ 1.588524912441221472814692121069821695547e392L,
+ -3.904354800861715180218598151050191841308e395L,
+ 9.719938686092045781827273411668132975319e398L,
+ -2.450763621049522051234479737511375679283e402L,
+ 6.257892098396815305085674126334317095277e405L,
+ -1.618113552083806592527989531636955084420e409L,
+ 4.236528795217618357348618613216833722648e412L,
+ -1.123047068199051008086174989124136878992e416L,
+ 3.013971787525654770217283559392286666886e419L,
+ -8.188437573221553030375681429202969070420e422L,
+ 2.251910591336716809153958146725775718707e426L,
+ -6.268411292043789823075314151509139413399e429L,
+ 1.765990845202322642693572112511312471527e433L,
+ -5.035154436231331651259071296731160882240e436L,
+ 1.452779356460483245253765356664402207266e440L,
+ -4.241490890130137339052414960684151515166e443L,
+ 1.252966001692427774088293833338841893293e447L,
+ -3.744830047478272947978103227876747240343e450L,
+ 1.132315806695710930595876001089232216024e454L,
+ -3.463510845942701805991786197773934662578e457L,
+ 1.071643382649675572086865465873916611537e461L,
+ -3.353824475439933688957233489984711465335e464L,
+ 1.061594257145875875963152734129803268488e468L,
+ -3.398420969215528955528654193586189805265e471L,
+ 1.100192502000434096206138068020551065890e475L,
+ -3.601686379213993374332690210094863486472e478L,
+ 1.192235170430164900533187239994513019475e482L,
+ -3.990342751779668381699052942504119409180e485L,
+ 1.350281800938769780891258894167663309221e489L,
+ -4.619325443466054312873093650888507562249e492L,
+ 1.597522243968586548227514639959727696694e496L,
+ -5.584753729092155108530929002119620487652e499L,
+ 1.973443623104646193229794524759543752089e503L,
+ -7.048295441989615807045620880311201930244e506L,
+ 2.544236702499719094591873151590280263560e510L,
+ -9.281551595258615205927443367289948150345e513L,
+ 3.421757163154453657766296828520235351572e517L,
+ -1.274733639384538364282697627345068947433e521L,
+ 4.798524805311016034711205886780460173566e524L,
+ -1.825116948422858388787806917284878870034e528L,
+ 7.013667442807288452441777981425055613982e531L,
+ -2.723003862685989740898815670978399383114e535L,
+ 1.068014853917260290630122222858884658850e539L,
+ -4.231650952273697842269381683768681118533e542L,
+ 1.693650052202594386658903598564772900388e546L,
+ -6.846944855806453360616258582310883597678e549L,
+ 2.795809132238082267120232174243715559601e553L,
+ -1.153012972808983269106716828311318981951e557L,
+ 4.802368854268746357511997492039592697149e560L,
+ -2.019995255271910836389761734035403905781e564L,
+ 8.580207235032617856059250643095019760968e567L,
+ -3.680247942263468164408192134916355198549e571L,
+ 1.593924457586765331397457407661306895942e575L,
+ -6.970267175232643679233530367569943057501e578L,
+ 3.077528087427698518703282907890556154309e582L,
+ -1.371846760052887888926055417297342106614e586L,
+ 6.173627360829553396851763207025505289166e589L,
+ -2.804703130495506384463249394043486916669e593L,
+ 1.286250900087150126167490951216207186092e597L,
+ -5.954394420063617872366818601092036543220e600L,
+ 2.782297785278756426177542270854984091406e604L,
+ -1.312214674935307746141207680066262384215e608L,
+ 6.246299145383554153167974732783934504370e611L,
+ -3.000812007679574430883792565577444226490e615L,
+ 1.454904877136007844493861746476079537075e619L,
+ -7.118558521873800304612781121044077357278e622L,
+ 3.514739820897817389472822276832677887997e626L,
+ -1.751137068816377401163011262831890828437e630L,
+ 8.803498091818800678575314081978951179602e633L,
+ -4.465612911700593572269200981612564161010e637L,
+ 2.285494565287530681465757798517033542888e641L,
+ -1.180145168917737098025683613598595411329e645L,
+ 6.147941849198393232663105284575149616925e648L,
+ -3.231069156963603593233679426198974663352e652L,
+ 1.713042725635435041806895849197608270935e656L,
+ -9.161761363270648920537613435771882898051e659L,
+ 4.942675965960539112005679080810117766825e663L,
+ -2.689684712697383518131267222872386600031e667L,
+ 1.476320014229917759615308193449511534656e671L,
+ -8.173037740864781506597184122049453514594e674L,
+ 4.563462313190521363235182420178784459580e678L,
+ -2.569790015236158475703055501886439298708e682L,
+ 1.459410219452119981958355737832022375085e686L,
+ -8.358304882556983795372406183642486436653e689L,
+ 4.827305091483557818593092377664570208355e693L,
+ -2.811394311081493166793414157061950132403e697L,
+ 1.651026863340675349245561261339568827739e701L,
+ -9.776578579336866764167878646459810047899e704L,
+ 5.837207965197521880181236529616560780535e708L,
+ -3.513938957938032127105389702846371181520e712L,
+ 2.132747371360190507595748444536911078788e716L,
+ -1.305047363239192640729466563372665311602e720L,
+ 8.050825342678337497636292798039996484780e723L,
+ -5.006884161223862543665524155681082112689e727L,
+ 3.139016066011452177570812014513491361235e731L,
+ -1.983829535212711378291469356666001365873e735L,
+ 1.263822427649676371257598052486237628698e739L,
+ -8.115678659900522918802121684491754629503e742L,
+ 5.252995164972075271667364371449050412435e746L,
+ -3.427038125662404660056511738625477058135e750L,
+ 2.253446011834352733279946306835940729858e754L,
+ -1.493407341897034717876962786798831719683e758L,
+ 9.974681322653365118752729509398728354442e761L,
+ -6.714230142773850863927710112350816379426e765L,
+ 4.554668668931723346600337564274944733530e769L,
+ -3.113635386023220127834102980385275379533e773L,
+ 2.144945411287666204679363498162954050208e777L,
+ -1.488982121181387164932397544378555256016e781L,
+ 1.041537218854627455352298173588983048748e785L,
+ -7.341073881786613676177562822942175683993e788L,
+ 5.213524272587199574980117351016322518428e792L,
+ -3.730592531776514409283897139216167197989e796L,
+ 2.689592876341877079083449497724049500175e800L,
+ -1.953643788231947582529884602972233135002e804L,
+ 1.429691073080500563348668321308878246277e808L,
+ -1.054059177095488639836063073070536825675e812L,
+ 7.828919160938693948399336431565350676613e815L,
+ -5.857884457184396382550955498026762014753e819L,
+ 4.415401588264172474136969345712659422380e823L,
+ -3.352573884181287635796498822858109969161e827L,
+ 2.564210385719224000156548240934108974447e831L,
+ -1.975534392116037602837941409848663077528e835L,
+ 1.533062123975940045180943006948008486466e839L,
+ -1.198306160488763291730059994812781226903e843L,
+ 9.434034267770711698676321369174735725321e846L,
+ -7.480619200038505368468483892246806488879e850L,
+ 5.974161898439971564124576801455052907638e854L,
+ -4.805125663714699771668630995361572639386e858L,
+ 3.892332138028039952403812726744593073776e862L,
+ -3.175276505779699340738548328810180869575e866L,
+ 2.608608681939322393581069188271626122519e870L,
+ -2.158148554392732439392868052394994052628e874L,
+ 1.797993483301448477700600221980862686033e878L,
+ -1.508407575089108597171576068862286462909e882L,
+ 1.274273406242459482708930389008701147244e886L,
+ -1.083950475353171986748233157909397370193e890L,
+ 9.284292630726328432038470356821265395331e893L,
+ -8.007012115449516364480417355063446317414e897L,
+ 6.952871948429568933888979915833266241471e901L,
+ -6.078828929473797621198666799700739891205e905L,
+ 5.350908089710964244671334224708057812633e909L,
+ -4.742168072503284973969982758434401589090e913L,
+ 4.231149239401967697257534662010605751136e917L,
+ -3.800684612827828851942743291026898158947e921L,
+ 3.436984796314246158361599955909956583986e925L,
+ -3.128930718993658356398482705317381808301e929L,
+ //
+ // 602-1300: http://www.wolframalpha.com/input/?i=TABLE[N[Bernoulli[i]%2C40]%2C+{i%2C602%2C1300%2C2}]
+ 2.867524740577223817164663595437919813239e933L, -2.645462974939090580963101220449509725942e937L, 2.456800827789169780295419018499543141869e941L, -2.296690549725790064673528302231294870532e945L, 2.161174697699793265715182091764676666457e949L, -2.047023224586087259305754002882269123194e953L, 1.951604806042481282712736234132803700277e957L, -1.872785206668284042110390583158639495143e961L, 1.808847160923282257302788929692654262867e965L, -1.758427529634609613399327744595257497188e969L, 1.720468488019528147087036246754294757647e973L, -1.694180279355332648057740852839804839425e977L, 1.679013685251183870616469618951463869496e981L, -1.674640861433092946269144173974414945664e985L, 1.680943600147858322148767806987527412112e989L, -1.698008433134805056489370119323402510305e993L, 1.726128304411348354183882648263448448633e997L, -1.765810838736918108045764015629875016219e1001L, 1.817793526882665071123822455897912718293e1005L, -1.883066459765807128944897377914669600374e1009L, 1.962903588035940537938222992228124233567e1013L, -2.058903881920696086033171142046100185783e1017L, 2.173044241735786946064676598703393618281e1021L, -2.307746591425236218893160658331303115253e1025L, 2.465962312241418731528973526597433097256e1029L, -2.651278087802503406316742676403301581549e1033L, 2.868048395658440423778896607880692085708e1037L, -3.121561373094393453726645989392054731637e1041L, 3.418246710091027042099932753084126095820e1045L, -3.765936717592482928796920675282930034018e1049L, 4.174194967165213973474293718362757753877e1053L, -4.654731142471753017867105249805137855862e1057L, 5.221926310090434518253178454907900079787e1061L, -5.893500145664015254409680930288710794031e1065L, 6.691361332576333738130720616841706994101e1069L, -7.642695184575063524608775697714741180954e1073L, 8.781359617440634128952082759434723165820e1077L, -1.014968338800868135594698909567734048618e1082L, 1.180079105471061498849752479044520598414e1086L, -1.380162016721660241308046692646452732446e1090L, 1.623685158291375662775444238282343536948e1094L, -1.921404880943289359290531906131400049399e1098L, 2.287040419533950152851434188305457266969e1102L, -2.738162880206032093123060939173765335255e1106L, 3.297371307848643161532227459901386725801e1110L, -3.993854689967542662299211323085023297602e1114L, 4.865474805885735467044047308902313673643e1118L, -5.961554732739027308247618738765152679497e1122L, 7.346627151757492821447573639763873833441e1126L, -9.105493288459908620636712748727395637965e1130L, 1.135007867626164861991621396462821975167e1135L, -1.422876214067403769204874786137232627418e1139L, 1.793912271573925309173135913914667878908e1143L, -2.274542916104231188526120123855259514144e1147L, 2.900273688809987694128857655036783261991e1151L, -3.719022795563122339874875448447744493398e1155L, 4.795753420982845153626611023078973364321e1159L, -6.218937220186281310109009529226561379773e1163L, 8.109611247999584815668395828940708619394e1167L, -1.063412316303440216539797215354141158589e1172L, 1.402214363674117662460496032135704328989e1176L, -1.859223235464558752766840772026058694872e1180L, 2.478828203789903637835992128856742276028e1184L, -3.323169416193176673655321536761413885767e1188L, 4.479640207312477092938541546776915956580e1192L, -6.071721672924085739424644485636889518799e1196L, 8.274698015123579607850404326757887762270e1200L, -1.133855131459773018024052539697784205966e1205L, 1.562146222050424344025824344480153248984e1209L, -2.163904570724750459592352173471446831752e1213L, 3.013703210722669908901286635073603018696e1217L, -4.219903244242308803914269531001720703294e1221L, 5.940703220571043642186808904696174833998e1225L, -8.408147464216029127243257448169774333631e1229L, 1.196419999747411909144144315499654470715e1234L, -1.711518922741148710381740436694440587059e1238L, 2.461434539630850545757453894977350505251e1242L, -3.558748530932574002484841810677232366801e1246L, 5.172525606281917297657859608800373729529e1250L, -7.557850217376323621984784308774476917753e1254L, 1.110141075986004209769735296234549704181e1259L, -1.639216556732622481406083885926912451281e1263L, 2.433138328152562628385514545400044125983e1267L, -3.630476645219033020888837165221286413171e1271L, 5.445289518636306992942604775585977779418e1275L, -8.209806424989072060381590985042272020067e1279L, 1.244209849774134691374848390346442737613e1284L, -1.895384488692308848372754844910263931874e1288L, 2.902272596647764894203369746806169285113e1292L, -4.466944174025026625137032739317650862593e1296L, 6.910485739507636504313238347702354354916e1300L, -1.074550085668784170644854815272144687769e1305L, 1.679419258904938802199084915274175753529e1309L, -2.638155207645646220849795321076977230763e1313L, 4.165284786632654168563096850610185378233e1317L, -6.609774274649031371770290191295685774584e1321L, 1.054194100570841329575393359295845860860e1326L, -1.689822316104196916970708778265725885275e1330L, 2.722340957904912685605914893019783431164e1334L, -4.407776313964403233676810178851005163725e1338L, 7.172436210641903635864868181569129834361e1342L, -1.172947440100495955246356688225986736990e1347L, 1.927745674072824377954824961348211728006e1351L, -3.184013467435655962214317208087993711563e1355L, 5.285045125125832341263897233405196808096e1359L, -8.815883582819232027207118521581424783107e1363L, 1.477818368424505276711779171224799759099e1368L, -2.489482576496570159333357550363134602876e1372L, 4.214292881345076419678976329218843808204e1376L, -7.169068531615459070909644981451297906220e1380L, 1.225513133750594558180516896275774441895e1385L, -2.105160827387119480607950260289853896637e1389L, 3.633787605672960549893307203363402915249e1393L, -6.302830804027849515239463308430185990705e1397L, 1.098521433860299633481449685364914115468e1402L, -1.923858597401607622723144320370279518600e1406L, 3.385512828549942051667348582951554570164e1410L, -5.986286250836771248147827011780631183980e1414L, 1.063572794668186370728928272374836554300e1419L, -1.898666684876492795233907174493757572290e1423L, 3.405627002840442789235393111726609930533e1427L, -6.137724140284450036591063946055819333244e1431L, 1.111411024660941507986132154479364267486e1436L, -2.022060876221034821890406900217875915949e1440L, 3.696248025817144690840539132103538834108e1444L, -6.788448439024998306316860676030442691610e1448L, 1.252615233049059554031883468823648511657e1453L, -2.322190433141265975888955985950824418729e1457L, 4.325200102353909846882217732999001735342e1461L, -8.093531903011880118699218269369570178812e1465L, 1.521558881878323790120983450270946857209e1470L, -2.873780311010933807686415826253380907421e1474L, 5.452903697278823304173192839252276211670e1478L, -1.039457922537509500320638240809547113575e1483L, 1.990610112724715126895008793014214505760e1487L, -3.829667853173777076954453401761025071562e1491L, 7.401624504283011888971231756333356050310e1495L, -1.437075122764477911733220492562365990710e1500L, 2.802940275035867428066581228962104019228e1504L, -5.491938363067613321364335249495394164430e1508L, 1.080961960603953462180593404647115933651e1513L, -2.137290931892412298654741768897581319007e1517L, 4.245031321673807283498263276791307370788e1521L, -8.469499523038763989328773224520912663309e1525L, 1.697421812794203793865032206191322699261e1530L, -3.417217332563937242285349373774004020539e1534L, 6.910378594841763785923780822895851271770e1538L, -1.403696282437585785557998429691459557649e1543L, 2.864060533055333035232343601021192111053e1547L, -5.869818290384811353182423286543086530728e1551L, 1.208359745327224593486268988808338456906e1556L, -2.498576742140453770373914215325521001990e1560L, 5.189311407347546310078739863704346083861e1564L, -1.082537954843916294257278789980768336964e1569L, 2.268238255751421312559806122980932952706e1573L, -4.773557403917983369065731568732198697502e1577L, 1.009019097334998841920279535262007639746e1582L, -2.142181266523235177327239693359275472557e1586L, 4.567814904130855969979178320003286614868e1590L, -9.782550516204803195398428611221899469345e1594L, 2.104180123097086948576304557651398411373e1599L, -4.545658958087323864004652894518442709646e1603L, 9.862563944609427542603740078470901803131e1607L, -2.149105846582226970866569209122813809019e1612L, 4.703235567543888152049628411354542509156e1616L, -1.033719212601584878353206879472796545848e1621L, 2.281767401903848796732740825793310514456e1625L, -5.058236070813950229238666252351966279306e1629L, 1.126112519657857205642546937554224492775e1634L, -2.517766761987679577706779689880657777343e1638L, 5.653225190181653388317503182908983211029e1642L, -1.274735955461074142223278576503188429497e1647L, 2.886578974679460464298863945016671299242e1651L, -6.564203307141426181809363135003467581753e1655L, 1.499036144473064593308260681782048262301e1660L, -3.437714715599902386917108442954580869236e1664L, 7.916830957072777234152907034541325149479e1668L, -1.830850567422571420661248197094782575285e1673L, 4.251778280827419894527511469762091846660e1677L, -9.915182507286989818033146623995507108134e1681L, 2.321878208636697663781227497233334385222e1686L, -5.459879022461660582811365437190884471726e1690L, 1.289222044549922720398543474297554204559e1695L, -3.056819658344217799458557578658863826289e1699L, 7.277891759142725294172926258364455941365e1703L, -1.739928293433385104144012025546489673795e1708L, 4.176797408823713136137404972612780406904e1712L, -1.006788178307821554781930741698052910780e1717L, 2.436754569909644399766538111317379484511e1721L, -5.921896599028498715774458493117079340155e1725L, 1.445045688171565118619109316933316429671e1730L, -3.540547766876069233350621578795319652040e1734L, 8.710114552028472554054293344204504325978e1738L, -2.151484527880464463303897113553085899101e1743L, 5.335928195512405709733771642389502809087e1747L, -1.328726408335015910030370523083559660016e1752L, 3.322090527232917400247098823651437597786e1756L, -8.339387326241218096865362177688582376376e1760L, 2.101842203781264395369771906884644062395e1765L, -5.318704469415522036482913743767085545209e1769L, 1.351288005941730688647540059088127991581e1774L, -3.446853546858473171100748720136784228698e1778L, 8.827284762030783576089954173424852998700e1782L, -2.269642226090373319660782216907175419317e1787L, 5.858820683661708553422363777419430816755e1791L, -1.518385813684321665045387969920683656625e1796L, 3.950661327164595923092260035122668890334e1800L, -1.031976516347387969958181456058243183780e1805L, 2.706317892325103782207094286049104555552e1809L, -7.125140422584701175967252533378906957380e1813L, 1.883260203116768075569432925204868418472e1818L, -4.997193687108743666000994570700725873035e1822L, 1.331182722092654526185433799891693838871e1827L, -3.559930289076558484535632566755216035553e1831L, 9.557281027056970446117541983785660301558e1835L, -2.575805002229372523547972911961335317502e1840L, 6.969058431277067406841032797913179025984e1844L, -1.892842481279278678390672746902260183506e1849L, 5.160964211693777744707760614147460787285e1853L, -1.412602588198037643242529860614298968137e1858L, 3.881313379962387603749693387037174052146e1862L, -1.070542170988009009334148472388319844527e1867L, 2.964094312414144330805731101996829908435e1871L, -8.238350132106899955856124602934281976453e1875L, 2.298504171050560756192352106062598639825e1880L, -6.437303944649223478093890316531995121228e1884L, 1.809727811843121957353712606428292269805e1889L, -5.107047553992257935533518628886728031061e1893L, 1.446674478990385642488446075734631327506e1898L, -4.113513327511444762766719175770513771122e1902L, 1.174067517257431444028448391638451935667e1907L, -3.363630086409895071362533854123306097827e1911L, 9.672868956071838221096869293070568259792e1915L, -2.792101741911955365960369780457612630184e1920L, 8.089710604557382430162031502761771390568e1924L, -2.352650988877130983061761312962677887796e1929L, 6.867549079740051556501575104006222995568e1933L, -2.012161201632998475706904405535757516336e1938L, 5.917489529279588702317256137229398357271e1942L, -1.746718667239329545125902248821502764273e1947L, 5.175069416058975040990816515838893249437e1951L, -1.538913401594651457295303469904084052963e1956L, 4.593185746210984655636051293374195150815e1960L, -1.375981868450401919299150690829612124045e1965L, 4.137207965217520410530508053863759216958e1969L, -1.248518564582257710069294326648626362439e1974L, 3.781575291117895093413381897917341286951e1978L, -1.149575999691408110085856948595444100435e1983L, 3.507413095836612229403470531176947165451e1987L, -1.074032838410645352804690949680310176413e1992L, 3.300857202456564870338466973024760446263e1996L, -1.018149578840803516349758843017979498322e2001L, 3.151876950233613792531594490714752800621e2005L, -9.792574827376149360558532022944033224780e2009L, 3.053456145978161645823454710737904504036e2014L, -9.555442346102849014299990542596620094035e2018L, 3.001037449298122384017009412541525703002e2023L, -9.459120112371096268275049056229023773120e2027L, 2.992168042152196502453442556462819104060e2032L, -9.498922680869041470681858599915282791899e2036L, 3.026307717971075309746179763189393755074e2041L, -9.676079238806159594565350708123427510151e2045L, 3.104778286352798464772361361434013339088e2050L, -9.997786802782252742109475924344598057966e2054L, 3.230847952724856366943939804248186203776e2059L, -1.047769651900498931701604323213605884945e2064L, 3.409958102134053489747140426163802214042e2068L, -1.113687894644055086152064258459886518528e2073L, 3.650114509271160332136458711252217684956e2077L, -1.200536387553969483433239131469825141412e2082L, 3.962482337718333099498977337189304099484e2086L, -1.312441206957064803437100929905979391106e2091L, 4.362246723746013772563799740886664288515e2095L, -1.454975881895253548422481637083633839534e2100L, 4.869831412214692119172895822285084162147e2104L, -1.635618419512383251104125916207188960680e2109L, 5.512611314145041257838234038980389596534e2113L, -1.864392957231340288547618808749072127289e2118L, 6.327317613106621547060670091824665547127e2122L, -2.154772001506498703267302897994526372056e2127L, 7.363426139490286496267931634843475368903e2131L, -2.524950643808031915843604894357998905460e2136L, 8.687956390288096215918373666581638675156e2140L, -2.999656978200020459428228924242615592768e2145L, 1.039231328851609224822335039430898644149e2150L, -3.612742437616019936358910410005123924796e2154L, 1.260211309932738404790711574105022002093e2159L, -4.410916378453971105434385837025433805752e2163L, 1.549140617923265948720013792673729394719e2168L, -5.459173749226782924959103886664322964926e2172L, 1.930343307630952098252884031069043541182e2177L, -6.848749229218425353808144618581305978045e2181L, 2.438117138001365487681440577590059588102e2186L, -8.708873656769794358508423272379627581292e2190L, 3.121268068338199458891764932384819739714e2195L, -1.122430216307539309816165910733145404999e2200L, 4.049900779207199370582177687160985635615e2204L, -1.466167983141158219266077836130256565915e2209L, 5.325678718693772500250292767751070974887e2213L, -1.940955845102272053048140384364058448998e2218L, 7.097467198361219669927211698104447309186e2222L, -2.603968771680987683436428778397387110896e2227L, 9.585403285394812946713320044815117440444e2231L, -3.540176030547640510648455468270569908446e2236L, 1.311827683984025111744358347783996339730e2241L, -4.877124229155333857009747836542843294702e2245L, 1.819213075760490882591173222316749809951e2250L, -6.808221630329265915405178596748950929642e2254L, 2.556299969544109052724772800143396857058e2259L, -9.629763347675306704861859899230073979116e2263L, 3.639508580119285595844040783082958425575e2268L, -1.380037493555816309137481185927387732499e2273L, 5.249980712165216709135893538080020409581e2277L, -2.003737844109055078145975651407367170529e2282L, 7.672522280806944397358668566379646540213e2286L, -2.947454993639165318799389781921184991045e2291L, 1.135966912801707623489383623092951142963e2296L, -4.392293711194501621873299212059053651432e2300L, 1.703813210168560937608104155973968112409e2305L, -6.630636743874062041158387022015853902938e2309L, 2.588742636486379690203698247275411406029e2314L, -1.013959594068423546627946242481463893979e2319L, 3.984265821528043268586235974854766821078e2323L, -1.570614519682157047612769672066387881154e2328L, 6.211297381339606877062824459742129064477e2332L, -2.464246931985476159686671650962783785426e2337L, 9.807833742601662212615240518855757197483e2341L, -3.916036434571217691317276306031837539092e2346L, 1.568566392975837368624727722120313955274e2351L, -6.302885887601142677858008037129298948063e2355L, 2.540704455306077495480843691828334210014e2360L, -1.027412480318234348899627142408950111875e2365L, 4.167823618450297116765978030480648316769e2369L, -1.696076602731914277275203926124423530377e2374L, 6.923904505633301788461482786634220738504e2378L, -2.835463065742506394026733592206185459035e2383L, 1.164828772275756526225951620927486307632e2388L, -4.800242878545012539781545966693324656699e2392L, 1.984381759611877246529319121941597679107e2397L, -8.228979942542641498511023600269641046627e2401L, 3.423130231367101727862739208673375060101e2406L, -1.428418168129733054582191895023094524495e2411L, 5.979153801634459282232521647160044877770e2415L, -2.510581926948409809562349588087762800160e2420L, 1.057443785053915411991029410076722022815e2425L, -4.467723713549428749678277264414266162837e2429L, 1.893474116528533144079731251913008472748e2434L, -8.049601965052954947260081891142509464888e2438L, 3.432648527503971149009691133946275281368e2443L, -1.468324699963694393989960228042259134294e2448L,
+ //
+ // 1302-1600: http://www.wolframalpha.com/input/?i=TABLE[N[Bernoulli[i]%2C40]%2C+{i%2C1302%2C1600%2C2}]
+ 6.300146502435743791500010801885493871234e2452L, -2.711520667146768856688291798851999580833e2457L, 1.170595555513900137297344452318266434006e2462L, -5.069095411973246242900074508988493530542e2466L, 2.201819284807954055092117706033113168896e2471L, -9.593088725189386197503123561368325167085e2475L, 4.192362385909155628936230811010649614060e2480L, -1.837725836941968309866675158105812946762e2485L, 8.080201101491972605313807752565294881374e2489L, -3.563536075527215702966392543784039539240e2494L, 1.576361051321107275181955665159661781175e2499L, -6.994292466180175594372663323941761853364e2503L, 3.112744353537336702834647901141392426258e2508L, -1.389481328370627358752727485697345194612e2513L, 6.221134636655213696041740685131223999953e2517L, -2.793779613656947577224654924852010601105e2522L, 1.258399062987759035354039924686781081603e2527L, -5.685208194704131918461885165870560583895e2531L, 2.576167857759537340210434756292816456179e2536L, -1.170846052338591953257169251219597581763e2541L, 5.337296787116189575571202979672747140313e2545L, -2.440264475369219459038748840841422948951e2550L, 1.119037151526195093932933161706501865175e2555L, -5.146858829220973887154576240993607686435e2559L, 2.374259791963193693837576781321391741634e2564L, -1.098501215269400934956638118646657823799e2569L, 5.097500369683616795005376807036889542869e2573L, -2.372446971688020647583535886090779018865e2578L, 1.107430282014636546248612381377039463753e2583L, -5.184597227131050012643138079903381280471e2587L, 2.434392040100910394476893838832599310265e2592L, -1.146412753331162872665743308094817095949e2597L, 5.414578104816988124950636101250217797539e2601L, -2.564835392810685332173156758121489913946e2606L, 1.218495070518549208066544111736985586178e2611L, -5.805713573821806672815019495319510297824e2615L, 2.774298194574319430697819781128985128618e2620L, -1.329580186505564627453485444017911980430e2625L, 6.390545858902318479863947547243743500916e2629L, -3.080502542499571035376377703435361520427e2634L, 1.489236104239976282318361008292980814533e2639L, -7.220413839991892382038608955317126799684e2643L, 3.510874916591640642524021216241607185085e2648L, -1.712070118580404599831061485055269100525e2653L, 8.372956919832386730490070625622785478703e2657L, -4.106629146981883685523102256292669054596e2662L, 2.019945438530802964718619732330776495740e2667L, -9.964133277392242111939720494354938982970e2671L, 4.929278642971447854669801547226335041410e2676L, -2.445509657169810919463982615395074704130e2681L, 1.216734421265677299127016883839223226884e2686L, -6.071008437677720186241562251151490713584e2690L, 3.037824949882992896564570441252792097027e2695L, -1.524402878612630565501569310883356490225e2700L, 7.671320530781999359200097739951316234193e2704L, -3.871436167706734376478728954716915204399e2709L, 1.959313530432202158587932399068682252335e2714L, -9.944063618400630821320953821427307024297e2718L, 5.061161998202463346818982228476199873781e2723L, -2.583219090831132705328958245740715185448e2728L, 1.322193991367293532684189527174543501836e2733L, -6.786569982732483290873213417465458376706e2737L, 3.493212334804776543395067018414547811062e2742L, -1.803090099978261928508495412750404640933e2747L, 9.333100843930216567894508007158644926767e2751L, -4.844499031405982604449146511179496492045e2756L, 2.521648090959971240812330574936006906830e2761L, -1.316227870932708474838173333385377250286e2766L, 6.889488826832738674261056521130795910494e2770L, -3.616184242864384509259984293501533623932e2775L, 1.903356124758119137116543283603627028779e2780L, -1.004601544584640657081847200643996069583e2785L, 5.317043885597842225603585588404817559596e2789L, -2.821938866752488868682751438901900485500e2794L, 1.501842023003449590337997900945924161741e2799L, -8.014908048137216649348740300633172710524e2803L, 4.289126235121619907138036129192558937445e2808L, -2.301619137231461344870820700320913118444e2813L, 1.238485136850053215006962645111854705210e2818L, -6.682503731149007943059244518074044280490e2822L, 3.615572393938012932030234169574978859655e2827L, -1.961565108627429629104703146282982075623e2832L, 1.067123259692924564435881096382837264046e2837L, -5.821179870182035246401397327057170726418e2841L, 3.184127229476322727732208017279268211356e2846L, -1.746429902183019597973436257300843998825e2851L, 9.604873565299766333876882842813498685054e2855L, -5.296759978724702692134960752308186890356e2860L, 2.928906353338652198977536576170287112391e2865L, -1.623961162577704769945821804737884742792e2870L, 9.028574047002736235613238355032484299017e2874L, -5.033087486357905828950503441308068892610e2879L, 2.813325650062267479031371852434194635210e2884L, -1.576791132296320840138263753339056345362e2889L, 8.861258343945925667272164531504265693289e2893L, -4.993236404321511029440212686547068244002e2898L, 2.821192993950901287717082243608730217471e2903L, -1.598254169674379493385730199445427966752e2908L, 9.078617590346932363947095804057608979359e2912L, -5.170742114456472142154347566092068443393e2917L, 2.952866185102528847516095880416675972086e2922L, -1.690794578626103552690094140317813413244e2927L, 9.707168799669516048238542260085175133847e2931L, -5.587884732306715493795271931175883605707e2936L, 3.225179489154957423492905957887744116530e2941L, -1.866424419669188178697802576490431604300e2946L, 1.082967626854618222657109354056973072044e2951L, -6.300392007169862865282706277272018077291e2955L, 3.675066377245428685118763485986517510658e2960L, -2.149348371085132073107516253339849053182e2965L, 1.260349351812619395000600434630904474324e2970L, -7.409963623771231302980906971935254993610e2974L, 4.367980758467862686643231700861155889684e2979L, -2.581566823350789671250829457603555544100e2984L, 1.529757357568342629912560827243282062227e2989L, -9.088595394263364554625061567617375176719e2993L, 5.413829169254585648363594604231030415354e2998L, -3.233288119606092759447005827969216281573e3003L, 1.936042437734875803183915765854038424658e3008L, -1.162289934202291715747729318797398221667e3013L, 6.995870350500567071550614251287615697508e3017L, -4.221776496490106417392945233048068288503e3022L, 2.554309239868912570382343877718991746122e3027L, -1.549440871550119801225143558087410562418e3032L, 9.423199525954784955533959981278992475051e3036L, -5.745689660772387668861183913170050552119e3041L, 3.512407521007240798565045328376471603253e3046L, -2.152708113797517364614914569890010876143e3051L, 1.322761289733739440340237168659770154654e3056L, -8.148777388506488753591136948542248584098e3060L, 5.032880858479326069741729004270784264612e3065L, -3.116396010103058126269735274818345780360e3070L, 1.934634831148214353514796782480703021435e3075L, -1.204077166243116651938489240924641810276e3080L, 7.513065583444964704795707060501161621868e3084L, -4.699873512563164914493150520500838535415e3089L, 2.947541197349762411713872934523813866703e3094L, -1.853262416286420077763886100673646141885e3099L, 1.168196427912100545575264493997591040800e3104L, -7.382362285873345348505276546404015842875e3108L, 4.677071041058096429847797962954927487730e3113L, -2.970642034084362431442183248944824506476e3118L, 1.891572688282564476274920103912259755482e3123L, -1.207509963440193713810418554061532113326e3128L, 7.727731208240101791845515599659441557781e3132L, -4.957988488048495669466804712012179891532e3137L, 3.188965862446236259925047956715566822864e3142L, -2.056286895821370106507670239256782411337e3147L, 1.329246918771714093479509313343886287414e3152L, -8.614188519577835653765633797787633659253e3156L,
+ //
+ // 1602-1900: http://www.wolframalpha.com/input/?i=TABLE[N[Bernoulli[i]%2C40]%2C+{i%2C1602%2C1900%2C2}]
+ 5.596396533621874175909933615343145642161e3161L, -3.644908483469388437457938883454376864180e3166L, 2.379838409026860469990569665632800095988e3171L, -1.557720925267669865362152155022069166772e3176L, 1.022143420270029721682551084917730373739e3181L, -6.723767358891570842116651998814252095792e3185L, 4.433950491570308179905446963723780229747e3190L, -2.931196854668917448553150023532223509373e3195L, 1.942557068752664549549945921392100172355e3200L, -1.290553202978622786891265558106235068695e3205L, 8.595082329732118303768775883557789195136e3209L, -5.738453265222970049867280061719670658457e3214L, 3.840687915100689856736926915331157331684e3219L, -2.576862441955523551149886625900059307506e3224L, 1.733166107320377310388765047659987844208e3229L, -1.168569552450178559412843683052610870569e3234L, 7.898289836694980777809433306209459851871e3238L, -5.351485909164216694400535493924387979018e3243L, 3.634772439350395177931952925644409735777e3248L, -2.474801048002975145046569303233576339695e3253L, 1.689126939254790850063878942448569759390e3258L, -1.155691524500722774057997965355407962525e3263L, 7.926435404542361405718288670391575676323e3267L, -5.449654814183048796524718620178906854846e3272L, 3.755898589900254795894812942275711835138e3277L, -2.594843902682143854622514329649211211808e3282L, 1.797048752397789969347915328338360264536e3287L, -1.247551415074438712713815166107969504456e3292L, 8.681719521514448143910215886388510318746e3296L, -6.056203898213120922016159444227958572276e3301L, 4.234882876331814099029781995617143573641e3306L, -2.968432911643338866295929748049749932906e3311L, 2.085723508930484816454740610260790948864e3316L, -1.469023169879432026361623513301566735138e3321L, 1.037150346505052892302077637883522696572e3326L, -7.339977067836656769144838365069396168014e3330L, 5.206985412168234130596004552956337839140e3335L, -3.702673773319239583641029108403509825141e3340L, 2.639251227995760315076225206168354089692e3345L, -1.885736353072698581595150856674914203383e3350L, 1.350563292338261784288559687678302458996e3355L, -9.695749980998301526113046898985991802000e3359L, 6.977167462628398202151721319169989304520e3364L, -5.032768280399753942925624560483352299263e3369L, 3.638844963651800168080623511900705036698e3374L, -2.637228631269251606169613775399022890118e3379L, 1.915836351653767108720464847696767898597e3384L, -1.395064293615007319328267865803567670760e3389L, 1.018249052614943190644465556486933211307e3394L, -7.449662162606857550867922631658930320805e3398L, 5.463119632208085241594107781601567713991e3403L, -4.015736541676989144201935890497836963875e3408L, 2.958754190183866660901503059509579790900e3413L, -2.185096074054288399312733179064098492511e3418L, 1.617517444557020250864919655301189186103e3423L, -1.200170662015511746748935675940010250555e3428L, 8.925888349899029449015791684428724952411e3432L, -6.653851763691885517669938275618991145962e3437L, 4.971722031098457895973348076474071155918e3442L, -3.723500582577984967442020337848702786829e3447L, 2.795153783541721373364976034391375710110e3452L, -2.103141577212720698169118819883801186873e3457L, 1.586129575320959267959148073466004084241e3462L, -1.198988457279648730711646682156242973137e3467L, 9.084402368157025658430300252246526602197e3471L, -6.898927494435965163817354296023108913714e3476L, 5.251332286149361587885046891266325872375e3481L, -4.006442950956739933884502808470603581850e3486L, 3.063718202820270282280659950794978994604e3491L, -2.348215284130973783732145823834807395920e3496L, 1.803952490148087317330011096671019781340e3501L, -1.389022326803437345760911068933754707688e3506L, 1.071986115818329525986099441493200866389e3511L, -8.292085224650940719705699485423856363908e3515L, 6.428829064452939640541475198655560890344e3520L, -4.995654440302797445368056643032307686314e3525L, 3.890847042582299188849273838681034339406e3530L, -3.037288555751484681537442833929275697351e3535L, 2.376385803695694695338601696534348875191e3540L, -1.863527130251861900692886008704804849076e3545L, 1.464674913498036269270793715104706378182e3550L, -1.153804954579033578659954846698233083197e3555L, 9.109783835348935092264268296199541780964e3559L, -7.208869193983001804305451104827153729326e3564L, 5.717530734277611949162917337810749919265e3569L, -4.544970302634007326980094771330550661605e3574L, 3.621042850825283032134228901678636353355e3579L, -2.891447067949778492831490654980043715471e3584L, 2.314060419397710657435821461707043283167e3589L, -1.856140759923563235273220981623595304434e3594L, 1.492185412981476596273279338314204171587e3599L, -1.202290032627175365810126250991853594801e3604L, 9.708881154579770196658265042625239421053e3608L, -7.857809850747029705680072304049448493252e3613L, 6.373898598298513400228819113197728735438e3618L, -5.181780406472117449048907989647202286666e3623L, 4.222036621953044040518942750638183171221e3628L, -3.447728386429130175025813550845575613047e3633L, 2.821701521717856346224159586852612710800e3638L, -2.314488376711998526455043944505424906920e3643L, 1.902671298033180765286213227393060711096e3648L, -1.567603736821312488140289549008391847440e3653L, 1.294408945316538946551785312385509945367e3658L, -1.071194533081615830960091702262923009420e3663L, 8.884351908108581551151252566466606126397e3667L, -7.384866682828103669170236267589653324531e3672L, 6.152023838008155718180876735217718355563e3677L, -5.136304310431705506236573876510219357975e3682L, 4.297736808124296434723193397876220759378e3687L, -3.603994887745884762510172194982172483480e3692L, 3.028884745605031552399167746007361297342e3697L, -2.551141302205187365552982635794121855138e3702L, 2.153467982869535549299173317536193051608e3707L, -1.821769476343602094059466497311600827296e3712L, 1.544537580582347892980177956984101211006e3717L, -1.312358705945937257247030754517293537539e3722L, 1.117518229297781388884979995402355617235e3727L, -9.536820860779441793021624381677086661097e3731L, 8.156400668831968026931547065507466530546e3736L, -6.990984948728184142718575396052260691181e3741L, 6.005124901126818071638224144541102727563e3746L, -5.169500241880947716732682089328427995109e3751L, 4.459815478235310026240134567325749844182e3756L, -3.855902253361684187081283218890336962427e3761L, 3.340988024176995223515640815937037040546e3766L, -2.901099226680215736735094376078800376829e3771L, 2.524573363444334459448089563912567842927e3776L, -2.201659455716348555524529213295341212492e3781L, 1.924190302190936448078364755844591374353e3786L, -1.685313186099770223843319514432495898517e3791L, 1.479268235966730475749985741048766689808e3796L, -1.301205702893883803117530921635013780575e3801L, 1.147035071153450453405384269242743907426e3806L, -1.013300250456366849150496776951686112298e3811L, 8.970761720605591762300958007557533865346e3815L, -7.958829781488943084496783248922217392838e3820L, 7.076146954685024795720193943027902028642e3825L, -6.304798526260409199660290516451546966159e3830L, 5.629519616664188107056583939722984509867e3835L, -5.037281594099054092767959480843344929292e3840L, 4.516946091316834843581919268794683123349e3845L, -4.058975118925834202620358386772092359951e3850L, 3.655187798978978909014603682039470653549e3855L, -3.298555903041546671060101785513812175322e3860L, 2.983031738662727912016882399515879119620e3865L, -2.703403043317732979516341931451317866898e3870L, 2.455170460800096241793872443768546335444e3875L, -2.234443928432490538417605502448376856290e3880L, 2.037854924078003280537856980560782325730e3885L, -1.862482033918775734840779765743099458137e3890L,
+ //
+ // 1902-2200: http://www.wolframalpha.com/input/?i=TABLE[N[Bernoulli[i]%2C40]%2C+{i%2C1902%2C2200%2C2}]
+ 1.705787724951999960095629912416210969679e3895L, -1.565564556110550991891247404758895970376e3900L, 1.439889351869832939488618785632174464789e3905L, -1.327084102784257406218693901793045990520e3910L, 1.225682557296027075027021534960026145706e3915L, -1.134401635488994148555787301654561211982e3920L, 1.052116934052356802920509999705307165985e3925L, -9.778417073593082219082361206542342793584e3929L, 9.107088061888562704837019028349522303725e3934L, -8.499551364633102138471246155980056936129e3939L, 7.949082681085658044610890152056533167407e3944L, -7.449748809722797718736397140511396011691e3949L, 6.996307824769340144608141799981589288378e3954L, -6.584122718472954006131003060359621706243e3959L, 6.209086595833487707192492087176843233407e3964L, -5.867557793863165391821489909125720982339e3969L, 5.556303538475260373917478405626416604297e3974L, -5.272450955936249442242634142613834212778e3979L, 5.013444428433789818228792126117223030641e3984L, -4.777008429684552423800736200488532033034e3989L, 4.561115100786341787876705283291018781137e3994L, -4.363955932181992701667719449097126840439e3999L, 4.183917007557000586305945495258591147615e4004L, -4.019557342177353010692923286760895584096e4009L, 3.869589913635745758786275231296652917580e4014L, -3.732865038934070181861017140563175000872e4019L, 3.608355799736107390800162778737339576843e4024L, -3.495145258697474565347261083975193776541e4029L, 3.392415245050326563747729613872524362741e4034L, -3.299436517958948801426629481782413630714e4039L, 3.215560142306355508598119430378551642857e4044L, -3.140209934146377815556058799557727461298e4049L, 3.072875852591406752692761744649563131272e4054L, -3.013108231854799187724018548255922550991e4059L, 2.960512761914376268185064129600549308882e4064L, -2.914746139139036596123006476633770383901e4069L, 2.875512319506974985103149834921665445532e4074L, -2.842559316984704569380036093537576068104e4079L, 2.815676498441436148701483904115879856704e4084L, -2.794692334326268275058539147656334465534e4089L, 2.779472571396106785963004020814493340829e4094L, -2.769918800191406321625251621260024635680e4099L, 2.765967395840433013288935879837390099329e4104L, -2.767588816244119880300161388073836623878e4109L, 2.774787246856347651152278076466043136230e4114L, -2.787600586224957950622601135620189837948e4119L, 2.806100771288225169339048358106052817280e4124L, -2.830394446218080573456394167711739786431e4129L, 2.860623983452244712039094143642843717029e4134L, -2.896968870550611723525738907034588104300e4139L, 2.939647481737606306044335918078617963078e4144L, -2.988919258547518526076380181812161398808e4149L, 3.045087329976721023952450383837883029431e4154L, -3.108501609077197464748958150625867523408e4159L, 3.179562410123820875787052833975010965963e4164L, -3.258724638491880104953913719767939138170e4169L, 3.346502614347964869115073881474258766546e4174L, -3.443475601364631413158991572423086599816e4179L, 3.550294123121350747300886840907918182129e4184L, -3.667687162886053419715985091863398517145e4189L, 3.796470357354794420044278000297864085607e4194L, -3.937555311976846882455930574021795626971e4199L, 4.091960185075595842547638450930710467324e4204L, -4.260821710519620959138720129506770036460e4209L, 4.445408854703156440576808070360934740837e4214L, -4.647138333645908068599900650548418672065e4219L, 4.867592250805288922190809906525766574205e4224L, -5.108538156515551259475573296900660666192e4229L, 5.371951876776035157276013631113314852508e4234L, -5.660043513521220243900043448456234873940e4239L, 5.975287081834808618140945840817834710330e4244L, -6.320454323372684034118816565375206053746e4249L, 6.698653321371992324876559665938996023646e4254L, -7.113372643219128807424340495235606473967e4259L, 7.568531854202750881338746432078817214052e4264L, -8.068539383842553693076672384509126681464e4269L, 8.618358887685935324188596304168259394311e4274L, -9.223585437012291673660319256730398171887e4279L, 9.890533091606747031464718533600572123091e4284L, -1.062633567277107015128545384570274268438e4290L, 1.143906286231591191271274413511275981288e4295L, -1.233785411712565904499340744089870916842e4300L, 1.333307331840530219050170916015276125870e4305L, -1.443648758235403286296065629219598769529e4310L, 1.566147425967471851736562867318748510088e4315L, -1.702326086290842780634120184324081017286e4320L, 1.853920350455786350409148418966087344063e4325L, -2.022911043115598592197907512410632615740e4330L, 2.211561842992792253055716743938240466613e4335L, -2.422463130294011318178080247305407476096e4340L, 2.658583129381772791030436640519847627789e4345L, -2.923327636881988941081365085520742216540e4350L, 3.220609866329557159104267531058019683271e4355L, -3.554932228621330128152149026066400241546e4360L, 3.931482212643167323798366327390058684499e4365L, -4.356244944221399578650235478583297389113e4370L, 4.836135498303121165971331625888490168138e4375L, -5.379154636371461359750682662639062606297e4380L, 5.994572359716861309678596804350346692501e4385L, -6.693144535124290060793936095397161934045e4390L, 7.487368894313509797084395689517008597061e4395L, -8.391787970609807810531578161564037339793e4400L, 9.423348062978921203475110312003096820035e4405L, -1.060182516651648405903017734022504884319e4411L, 1.195033105063952979885086754342706651656e4416L, -1.349591538868673992167798923586925758429e4421L, 1.527028315253291113905307092657539132480e4426L, -1.731065051510920640409442255224015234974e4431L, 1.966076741510092840076264635935585216200e4436L, -2.237214093245750681191361238831105906202e4441L, 2.550550094903891445719729187215253324232e4446L, -2.913255853313667303707651906277658164129e4451L, 3.333811847072394764285817140850092324169e4456L, -3.822262084288044913490118858492563410392e4461L, 4.390520310533864198186202368026630430120e4466L, -5.052739449335052080092114976206610871466e4471L, 5.825757966350870043117899492954521458799e4476L, -6.729639942938203582008846884575881320532e4481L, 7.788329466816396015493306357116312471970e4486L, -9.030444674469025073047417528762134025409e4491L, 1.049024263381993629167658236142000524752e4497L, -1.220879351508964912255081664072251573277e4502L, 1.423541151220109512749655991050110438471e4507L, -1.662940118618541616964708044356967429362e4512L, 1.946219185900482116137855064775635250366e4517L, -2.281995008842006909631764011781911322493e4522L, 2.680678198213108543648324254258111216040e4527L, -3.154866427472784086389609599207759103500e4532L, 3.719827710160801797530420206201570269720e4537L, -4.394095404360277919140027580071549980218e4542L, 5.200201854779615608741690339830306148442e4547L, -6.165584312943608652377791415603277251516e4552L, 7.323705248531382981433751104158852636445e4557L, -8.715439846124090647163930834760361817820e4562L, 1.039079696609215651011736087603304766850e4568L, -1.241105689556982425619608247473478857800e4573L, 1.485143079696380339521658550262280772546e4578L, -1.780437412164973637340821168154300094802e4583L, 2.138372099157518882088209435171770222745e4588L, -2.572985071149069551034276570909360759588e4593L, 3.101615379617643734762997559011097203354e4598L, -3.745713657616368229906151946770042703357e4603L, 4.531859496161940719835150033082561700677e4608L, -5.493040495326927998321538336584233566465e4613L, 6.670262730603009306595018122252730741798e4618L, -8.114581584793494903775255213273982440688e4623L, 9.889666561810883044159054730371102725871e4628L, -1.207504541653929734716275932570097623330e4634L, 1.477021377885843688233899471354959308782e4639L, -1.809984912147908767583043524070645821179e4644L,
+ //
+ // 2202-2320: http://www.wolframalpha.com/input/?i=TABLE[N[Bernoulli[i]%2C40]%2C+{i%2C2202%2C2320%2C2}]
+ 2.222043594325228980916360265527780300093e4649L, -2.732869701246338361699515268224049951411e4654L, 3.367233945421922463553518272642397177145e4659L, -4.156377225041273602431272489314020150392e4664L, 5.139764368092890466235162431795350591151e4669L, -6.367329693760865476879589228002216011370e4674L, 7.902356742934106007362514378717026407839e4679L, -9.825176966314431712897976595483070301406e4684L, 1.223792760178593282435724837135946867088e4690L, -1.527068151452750404853140815207477555192e4695L, 1.908935682572268829496101580401263597905e4700L, -2.390593888616966248780378941331847473699e4705L, 2.999171106576893833644521002894489856321e4710L, -3.769440655453736670024798444784356437578e4715L, 4.746047769851891438576002047529258107351e4720L, -5.986405469241447720766576164546767533359e4725L, 7.564466155536872051712519119999711534616e4730L, -9.575641408047918720040356745796976488951e4735L, 1.214322951835035451699619713803395497423e4741L, -1.542682591979864353012093794301924196234e4746L, 1.963334539793192183270983986567556358603e4751L, -2.503148969013901182572118121398034622584e4756L, 3.197076711250102964526567664729089847162e4761L, -4.090653552025822488578293526174572934858e4766L, 5.243302769651520536759521264615159906699e4771L, -6.732697170903775309261288127044088674182e4776L, 8.660529543801770516930589210020128142543e4781L, -1.116015823611149634592870112730519454113e4787L, 1.440675306432920129218036927923030695520e4792L, -1.863078034853256227415397798026969938881e4797L, 2.413595413458810442409656314019115041699e4802L, -3.132317029597258599678590012779717945144e4807L, 4.072246763371584312534474102756137619716e4812L, -5.303577511521827157146305369181950467569e4817L, 6.919417518688636032335131253584331645491e4822L, -9.043473312934241153732087612484569398979e4827L, 1.184037400265044213826044590639924237359e4833L, -1.552956685415800894409743993367334099777e4838L, 2.040404893052952221581694807126473204625e4843L, -2.685565763841580219033402331219206776210e4848L, 3.540927057361929050327811875290025248120e4853L, -4.676912607538885419407656762767991163574e4858L, 6.188165903566760647569323704623433330229e4863L, -8.202087471895029964699042637255411806373e4868L, 1.089045274355389654614196651761310970580e4874L, -1.448524684976553869119447042300206226148e4879L, 1.930028100376784839502387280956424581974e4884L, -2.576074799096023589462128312524664980682e4889L, 3.444369635011990347297134928452972402038e4894L, -4.613354441299253694113609154769978684993e4899L, 6.189834306866879018555349507257537840922e4904L, -8.319470760665157534580593571258276368233e4909L, 1.120124240070996761986102680587384813245e4915L, -1.510740451399746828351090108638980398124e4920L, 2.041108231091323198877509959371257503819e4925L, -2.762447751447012472733302936575873838539e4930L,
+#endif
+ }};
+
+ return bernoulli_data[n];
+}
+
+template <class T>
+inline T unchecked_bernoulli_imp(std::size_t n, const mpl::int_<4>& )
+{
+ //
+ // Special case added for multiprecision types that have no conversion from long long,
+ // there are very few such types, but mpfr_class is one.
+ //
+ static const boost::array<boost::int32_t, 1 + max_bernoulli_b2n<T>::value> numerators =
+ {{
+ boost::int32_t( +1LL),
+ boost::int32_t( +1LL),
+ boost::int32_t( -1LL),
+ boost::int32_t( +1LL),
+ boost::int32_t( -1LL),
+ boost::int32_t( +5LL),
+ boost::int32_t( -691LL),
+ boost::int32_t( +7LL),
+ boost::int32_t( -3617LL),
+ boost::int32_t( +43867LL),
+ boost::int32_t( -174611LL),
+ boost::int32_t( +854513LL),
+ }};
+
+ static const boost::array<boost::int32_t, 1 + max_bernoulli_b2n<T>::value> denominators =
+ {{
+ boost::int32_t( 1LL),
+ boost::int32_t( 6LL),
+ boost::int32_t( 30LL),
+ boost::int32_t( 42LL),
+ boost::int32_t( 30LL),
+ boost::int32_t( 66LL),
+ boost::int32_t( 2730LL),
+ boost::int32_t( 6LL),
+ boost::int32_t( 510LL),
+ boost::int32_t( 798LL),
+ boost::int32_t( 330LL),
+ boost::int32_t( 138LL),
+ }};
+ return T(numerators[n]) / T(denominators[n]);
+}
+
+} // namespace detail
+
+template<class T>
+inline T unchecked_bernoulli_b2n(const std::size_t n)
+{
+ typedef mpl::int_<detail::bernoulli_imp_variant<T>::value> tag_type;
+
+ return detail::unchecked_bernoulli_imp<T>(n, tag_type());
+}
+
+}} // namespaces
+
+#endif // BOOST_MATH_UNCHECKED_BERNOULLI_HPP
diff --git a/boost/math/special_functions/detail/unchecked_factorial.hpp b/boost/math/special_functions/detail/unchecked_factorial.hpp
index eb8927a..3c23d6e 100644
--- a/boost/math/special_functions/detail/unchecked_factorial.hpp
+++ b/boost/math/special_functions/detail/unchecked_factorial.hpp
@@ -15,7 +15,9 @@
#pragma warning(push) // Temporary until lexical cast fixed.
#pragma warning(disable: 4127 4701)
#endif
+#ifndef BOOST_MATH_NO_LEXICAL_CAST
#include <boost/lexical_cast.hpp>
+#endif
#ifdef BOOST_MSVC
#pragma warning(pop)
#endif
@@ -266,6 +268,196 @@ struct max_factorial<long double>
BOOST_STATIC_CONSTANT(unsigned, value = 170);
};
+#ifdef BOOST_MATH_USE_FLOAT128
+
+template <>
+inline BOOST_MATH_FLOAT128_TYPE unchecked_factorial<BOOST_MATH_FLOAT128_TYPE>(unsigned i)
+{
+ static const boost::array<BOOST_MATH_FLOAT128_TYPE, 171> factorials = { {
+ 1,
+ 1,
+ 2,
+ 6,
+ 24,
+ 120,
+ 720,
+ 5040,
+ 40320,
+ 362880.0Q,
+ 3628800.0Q,
+ 39916800.0Q,
+ 479001600.0Q,
+ 6227020800.0Q,
+ 87178291200.0Q,
+ 1307674368000.0Q,
+ 20922789888000.0Q,
+ 355687428096000.0Q,
+ 6402373705728000.0Q,
+ 121645100408832000.0Q,
+ 0.243290200817664e19Q,
+ 0.5109094217170944e20Q,
+ 0.112400072777760768e22Q,
+ 0.2585201673888497664e23Q,
+ 0.62044840173323943936e24Q,
+ 0.15511210043330985984e26Q,
+ 0.403291461126605635584e27Q,
+ 0.10888869450418352160768e29Q,
+ 0.304888344611713860501504e30Q,
+ 0.8841761993739701954543616e31Q,
+ 0.26525285981219105863630848e33Q,
+ 0.822283865417792281772556288e34Q,
+ 0.26313083693369353016721801216e36Q,
+ 0.868331761881188649551819440128e37Q,
+ 0.29523279903960414084761860964352e39Q,
+ 0.103331479663861449296666513375232e41Q,
+ 0.3719933267899012174679994481508352e42Q,
+ 0.137637530912263450463159795815809024e44Q,
+ 0.5230226174666011117600072241000742912e45Q,
+ 0.203978820811974433586402817399028973568e47Q,
+ 0.815915283247897734345611269596115894272e48Q,
+ 0.3345252661316380710817006205344075166515e50Q,
+ 0.1405006117752879898543142606244511569936e52Q,
+ 0.6041526306337383563735513206851399750726e53Q,
+ 0.265827157478844876804362581101461589032e55Q,
+ 0.1196222208654801945619631614956577150644e57Q,
+ 0.5502622159812088949850305428800254892962e58Q,
+ 0.2586232415111681806429643551536119799692e60Q,
+ 0.1241391559253607267086228904737337503852e62Q,
+ 0.6082818640342675608722521633212953768876e63Q,
+ 0.3041409320171337804361260816606476884438e65Q,
+ 0.1551118753287382280224243016469303211063e67Q,
+ 0.8065817517094387857166063685640376697529e68Q,
+ 0.427488328406002556429801375338939964969e70Q,
+ 0.2308436973392413804720927426830275810833e72Q,
+ 0.1269640335365827592596510084756651695958e74Q,
+ 0.7109985878048634518540456474637249497365e75Q,
+ 0.4052691950487721675568060190543232213498e77Q,
+ 0.2350561331282878571829474910515074683829e79Q,
+ 0.1386831185456898357379390197203894063459e81Q,
+ 0.8320987112741390144276341183223364380754e82Q,
+ 0.507580213877224798800856812176625227226e84Q,
+ 0.3146997326038793752565312235495076408801e86Q,
+ 0.1982608315404440064116146708361898137545e88Q,
+ 0.1268869321858841641034333893351614808029e90Q,
+ 0.8247650592082470666723170306785496252186e91Q,
+ 0.5443449390774430640037292402478427526443e93Q,
+ 0.3647111091818868528824985909660546442717e95Q,
+ 0.2480035542436830599600990418569171581047e97Q,
+ 0.1711224524281413113724683388812728390923e99Q,
+ 0.1197857166996989179607278372168909873646e101Q,
+ 0.8504785885678623175211676442399260102886e102Q,
+ 0.6123445837688608686152407038527467274078e104Q,
+ 0.4470115461512684340891257138125051110077e106Q,
+ 0.3307885441519386412259530282212537821457e108Q,
+ 0.2480914081139539809194647711659403366093e110Q,
+ 0.188549470166605025498793226086114655823e112Q,
+ 0.1451830920282858696340707840863082849837e114Q,
+ 0.1132428117820629783145752115873204622873e116Q,
+ 0.8946182130782975286851441715398316520698e117Q,
+ 0.7156945704626380229481153372318653216558e119Q,
+ 0.5797126020747367985879734231578109105412e121Q,
+ 0.4753643337012841748421382069894049466438e123Q,
+ 0.3945523969720658651189747118012061057144e125Q,
+ 0.3314240134565353266999387579130131288001e127Q,
+ 0.2817104114380550276949479442260611594801e129Q,
+ 0.2422709538367273238176552320344125971528e131Q,
+ 0.210775729837952771721360051869938959523e133Q,
+ 0.1854826422573984391147968456455462843802e135Q,
+ 0.1650795516090846108121691926245361930984e137Q,
+ 0.1485715964481761497309522733620825737886e139Q,
+ 0.1352001527678402962551665687594951421476e141Q,
+ 0.1243841405464130725547532432587355307758e143Q,
+ 0.1156772507081641574759205162306240436215e145Q,
+ 0.1087366156656743080273652852567866010042e147Q,
+ 0.103299784882390592625997020993947270954e149Q,
+ 0.9916779348709496892095714015418938011582e150Q,
+ 0.9619275968248211985332842594956369871234e152Q,
+ 0.942689044888324774562618574305724247381e154Q,
+ 0.9332621544394415268169923885626670049072e156Q,
+ 0.9332621544394415268169923885626670049072e158Q,
+ 0.9425947759838359420851623124482936749562e160Q,
+ 0.9614466715035126609268655586972595484554e162Q,
+ 0.990290071648618040754671525458177334909e164Q,
+ 0.1029901674514562762384858386476504428305e167Q,
+ 0.1081396758240290900504101305800329649721e169Q,
+ 0.1146280563734708354534347384148349428704e171Q,
+ 0.1226520203196137939351751701038733888713e173Q,
+ 0.132464181945182897449989183712183259981e175Q,
+ 0.1443859583202493582204882102462797533793e177Q,
+ 0.1588245541522742940425370312709077287172e179Q,
+ 0.1762952551090244663872161047107075788761e181Q,
+ 0.1974506857221074023536820372759924883413e183Q,
+ 0.2231192748659813646596607021218715118256e185Q,
+ 0.2543559733472187557120132004189335234812e187Q,
+ 0.2925093693493015690688151804817735520034e189Q,
+ 0.339310868445189820119825609358857320324e191Q,
+ 0.396993716080872089540195962949863064779e193Q,
+ 0.4684525849754290656574312362808384164393e195Q,
+ 0.5574585761207605881323431711741977155627e197Q,
+ 0.6689502913449127057588118054090372586753e199Q,
+ 0.8094298525273443739681622845449350829971e201Q,
+ 0.9875044200833601362411579871448208012564e203Q,
+ 0.1214630436702532967576624324188129585545e206Q,
+ 0.1506141741511140879795014161993280686076e208Q,
+ 0.1882677176888926099743767702491600857595e210Q,
+ 0.237217324288004688567714730513941708057e212Q,
+ 0.3012660018457659544809977077527059692324e214Q,
+ 0.3856204823625804217356770659234636406175e216Q,
+ 0.4974504222477287440390234150412680963966e218Q,
+ 0.6466855489220473672507304395536485253155e220Q,
+ 0.8471580690878820510984568758152795681634e222Q,
+ 0.1118248651196004307449963076076169029976e225Q,
+ 0.1487270706090685728908450891181304809868e227Q,
+ 0.1992942746161518876737324194182948445223e229Q,
+ 0.269047270731805048359538766214698040105e231Q,
+ 0.3659042881952548657689727220519893345429e233Q,
+ 0.5012888748274991661034926292112253883237e235Q,
+ 0.6917786472619488492228198283114910358867e237Q,
+ 0.9615723196941089004197195613529725398826e239Q,
+ 0.1346201247571752460587607385894161555836e242Q,
+ 0.1898143759076170969428526414110767793728e244Q,
+ 0.2695364137888162776588507508037290267094e246Q,
+ 0.3854370717180072770521565736493325081944e248Q,
+ 0.5550293832739304789551054660550388118e250Q,
+ 0.80479260574719919448490292577980627711e252Q,
+ 0.1174997204390910823947958271638517164581e255Q,
+ 0.1727245890454638911203498659308620231933e257Q,
+ 0.2556323917872865588581178015776757943262e259Q,
+ 0.380892263763056972698595524350736933546e261Q,
+ 0.571338395644585459047893286526105400319e263Q,
+ 0.8627209774233240431623188626544191544816e265Q,
+ 0.1311335885683452545606724671234717114812e268Q,
+ 0.2006343905095682394778288746989117185662e270Q,
+ 0.308976961384735088795856467036324046592e272Q,
+ 0.4789142901463393876335775239063022722176e274Q,
+ 0.7471062926282894447083809372938315446595e276Q,
+ 0.1172956879426414428192158071551315525115e279Q,
+ 0.1853271869493734796543609753051078529682e281Q,
+ 0.2946702272495038326504339507351214862195e283Q,
+ 0.4714723635992061322406943211761943779512e285Q,
+ 0.7590705053947218729075178570936729485014e287Q,
+ 0.1229694218739449434110178928491750176572e290Q,
+ 0.2004401576545302577599591653441552787813e292Q,
+ 0.3287218585534296227263330311644146572013e294Q,
+ 0.5423910666131588774984495014212841843822e296Q,
+ 0.9003691705778437366474261723593317460744e298Q,
+ 0.1503616514864999040201201707840084015944e301Q,
+ 0.2526075744973198387538018869171341146786e303Q,
+ 0.4269068009004705274939251888899566538069e305Q,
+ 0.7257415615307998967396728211129263114717e307Q,
+ } };
+
+ return factorials[i];
+}
+
+template <>
+struct max_factorial<BOOST_MATH_FLOAT128_TYPE>
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 170);
+};
+
+#endif
+
template <>
inline double unchecked_factorial<double>(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(double))
{
@@ -279,6 +471,8 @@ struct max_factorial<double>
value = ::boost::math::max_factorial<long double>::value);
};
+#ifndef BOOST_MATH_NO_LEXICAL_CAST
+
template <class T>
inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
{
@@ -403,6 +597,22 @@ struct max_factorial
BOOST_STATIC_CONSTANT(unsigned, value = 100);
};
+#else // BOOST_MATH_NO_LEXICAL_CAST
+
+template <class T>
+inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))
+{
+ return 1;
+}
+
+template <class T>
+struct max_factorial
+{
+ BOOST_STATIC_CONSTANT(unsigned, value = 0);
+};
+
+#endif
+
#ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION
template <class T>
const unsigned max_factorial<T>::value;