diff options
Diffstat (limited to 'boost/math/special_functions')
19 files changed, 540 insertions, 177 deletions
diff --git a/boost/math/special_functions/detail/bernoulli_details.hpp b/boost/math/special_functions/detail/bernoulli_details.hpp index c12a172b60..ce95034439 100644 --- a/boost/math/special_functions/detail/bernoulli_details.hpp +++ b/boost/math/special_functions/detail/bernoulli_details.hpp @@ -11,6 +11,7 @@ #include <boost/detail/lightweight_mutex.hpp> #include <boost/utility/enable_if.hpp> #include <boost/math/tools/toms748_solve.hpp> +#include <vector> #ifdef BOOST_HAS_THREADS @@ -271,6 +272,7 @@ struct fixed_vector : private std::allocator<T> T* begin()const { return m_data; } T* end()const { return m_data + m_used; } unsigned capacity()const { return m_capacity; } + void clear() { m_used = 0; } private: T* m_data; unsigned m_used, m_capacity; @@ -284,6 +286,7 @@ public: #if defined(BOOST_HAS_THREADS) && !defined(BOOST_MATH_NO_ATOMIC_INT) , m_counter(0) #endif + , m_current_precision(boost::math::tools::digits<T>()) {} typedef fixed_vector<T> container_type; @@ -436,6 +439,13 @@ public: // // Single threaded code, very simple: // + if(m_current_precision < boost::math::tools::digits<T>()) + { + bn.clear(); + tn.clear(); + m_intermediates.clear(); + m_current_precision = boost::math::tools::digits<T>(); + } if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); @@ -452,6 +462,13 @@ public: // 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(m_current_precision < boost::math::tools::digits<T>()) + { + bn.clear(); + tn.clear(); + m_intermediates.clear(); + m_current_precision = boost::math::tools::digits<T>(); + } if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); @@ -471,12 +488,22 @@ public: // // 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) + if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) + || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) { 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((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) + || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) { + if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()) + { + bn.clear(); + tn.clear(); + m_intermediates.clear(); + m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release); + m_current_precision = boost::math::tools::digits<T>(); + } if(start + n >= bn.size()) { std::size_t new_size = (std::min)((std::max)((std::max)(std::size_t(start + n), std::size_t(bn.size() + 20)), std::size_t(50)), std::size_t(bn.capacity())); @@ -539,6 +566,13 @@ public: // // Single threaded code, very simple: // + if(m_current_precision < boost::math::tools::digits<T>()) + { + bn.clear(); + tn.clear(); + m_intermediates.clear(); + m_current_precision = boost::math::tools::digits<T>(); + } 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())); @@ -563,6 +597,13 @@ public: // 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(m_current_precision < boost::math::tools::digits<T>()) + { + bn.clear(); + tn.clear(); + m_intermediates.clear(); + m_current_precision = boost::math::tools::digits<T>(); + } 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())); @@ -590,12 +631,22 @@ public: // // 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) + if((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) + || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) { 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((static_cast<std::size_t>(m_counter.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < start + n) + || (static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>())) { + if(static_cast<int>(m_current_precision.load(BOOST_MATH_ATOMIC_NS::memory_order_consume)) < boost::math::tools::digits<T>()) + { + bn.clear(); + tn.clear(); + m_intermediates.clear(); + m_counter.store(0, BOOST_MATH_ATOMIC_NS::memory_order_release); + m_current_precision = boost::math::tools::digits<T>(); + } 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())); @@ -627,18 +678,20 @@ private: // // The caches for Bernoulli and tangent numbers, once allocated, // these must NEVER EVER reallocate as it breaks our thread - // safety guarentees: + // safety guarantees: // fixed_vector<T> bn, tn; std::vector<T> m_intermediates; // The value at which we know overflow has already occurred for the Bn: std::size_t m_overflow_limit; #if !defined(BOOST_HAS_THREADS) + int m_current_precision; #elif defined(BOOST_MATH_NO_ATOMIC_INT) boost::detail::lightweight_mutex m_mutex; + int m_current_precision; #else boost::detail::lightweight_mutex m_mutex; - atomic_counter_type m_counter; + atomic_counter_type m_counter, m_current_precision; #endif }; diff --git a/boost/math/special_functions/detail/bessel_k0.hpp b/boost/math/special_functions/detail/bessel_k0.hpp index e209168e87..42041dc673 100644 --- a/boost/math/special_functions/detail/bessel_k0.hpp +++ b/boost/math/special_functions/detail/bessel_k0.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/tools/rational.hpp> @@ -148,5 +150,9 @@ T bessel_k0(T x, const Policy& pol) }}} // namespaces +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_BESSEL_K0_HPP diff --git a/boost/math/special_functions/detail/bessel_k1.hpp b/boost/math/special_functions/detail/bessel_k1.hpp index 0d17cd3057..c0c5478520 100644 --- a/boost/math/special_functions/detail/bessel_k1.hpp +++ b/boost/math/special_functions/detail/bessel_k1.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/tools/rational.hpp> @@ -144,5 +146,9 @@ T bessel_k1(T x, const Policy& pol) }}} // namespaces +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_BESSEL_K1_HPP diff --git a/boost/math/special_functions/detail/bessel_y0.hpp b/boost/math/special_functions/detail/bessel_y0.hpp index 533ab7c8a0..84b64aaa73 100644 --- a/boost/math/special_functions/detail/bessel_y0.hpp +++ b/boost/math/special_functions/detail/bessel_y0.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/special_functions/detail/bessel_j0.hpp> @@ -220,5 +222,9 @@ T bessel_y0(T x, const Policy& pol) }}} // namespaces +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_BESSEL_Y0_HPP diff --git a/boost/math/special_functions/detail/bessel_y1.hpp b/boost/math/special_functions/detail/bessel_y1.hpp index 8396f8fe11..7dc99c42cd 100644 --- a/boost/math/special_functions/detail/bessel_y1.hpp +++ b/boost/math/special_functions/detail/bessel_y1.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/special_functions/detail/bessel_j1.hpp> @@ -192,5 +194,9 @@ T bessel_y1(T x, const Policy& pol) }}} // namespaces +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_BESSEL_Y1_HPP diff --git a/boost/math/special_functions/detail/erf_inv.hpp b/boost/math/special_functions/detail/erf_inv.hpp index 35072d5155..4e48300521 100644 --- a/boost/math/special_functions/detail/erf_inv.hpp +++ b/boost/math/special_functions/detail/erf_inv.hpp @@ -8,6 +8,9 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4127) // Conditional expression is constant +#pragma warning(disable:4702) // Unreachable code: optimization warning #endif namespace boost{ namespace math{ @@ -334,31 +337,37 @@ struct erf_inv_initializer 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()); - // 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 std::numeric_limits<T>::digits is zero, we must not call + // our inituialization code here as the precision presumably + // varies at runtime, and will not have been set yet. + if(std::numeric_limits<T>::digits) + { + 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()); + // 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()); - // 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: + // 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(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()); + 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(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()); + 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{} }; @@ -532,5 +541,9 @@ inline typename tools::promote_args<T>::type erf_inv(T z) } // namespace math } // namespace boost +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_SF_ERF_INV_HPP diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp index b2dd7ea292..df1a047435 100644 --- a/boost/math/special_functions/detail/lanczos_sse2.hpp +++ b/boost/math/special_functions/detail/lanczos_sse2.hpp @@ -190,6 +190,22 @@ inline double lanczos13m53::lanczos_sum_expG_scaled<double>(const double& x) return t[0] / t[1]; } +#ifdef _MSC_VER + +BOOST_STATIC_ASSERT(sizeof(double) == sizeof(long double)); + +template <> +inline long double lanczos13m53::lanczos_sum<long double>(const long double& x) +{ + return lanczos_sum<double>(static_cast<double>(x)); +} +template <> +inline long double lanczos13m53::lanczos_sum_expG_scaled<long double>(const long double& x) +{ + return lanczos_sum_expG_scaled<double>(static_cast<double>(x)); +} +#endif + } // namespace lanczos } // namespace math } // namespace boost diff --git a/boost/math/special_functions/detail/polygamma.hpp b/boost/math/special_functions/detail/polygamma.hpp index 20a0292fb4..c0e4932907 100644 --- a/boost/math/special_functions/detail/polygamma.hpp +++ b/boost/math/special_functions/detail/polygamma.hpp @@ -11,7 +11,7 @@ #ifndef _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ #define _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ - #include <cmath> +#include <cmath> #include <limits> #include <boost/cstdint.hpp> #include <boost/math/policies/policy.hpp> @@ -27,7 +27,13 @@ #include <boost/static_assert.hpp> #include <boost/type_traits/is_convertible.hpp> - namespace boost { namespace math { namespace detail{ +#ifdef _MSC_VER +#pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) +#endif + +namespace boost { namespace math { namespace detail{ template<class T, class Policy> T polygamma_atinfinityplus(const int n, const T& x, const Policy& pol, const char* function) // for large values of x such as for x> 400 @@ -402,8 +408,18 @@ static boost::detail::lightweight_mutex m; boost::detail::lightweight_mutex::scoped_lock l(m); #endif + static int digits = tools::digits<T>(); static std::vector<std::vector<T> > table(1, std::vector<T>(1, T(-1))); + int current_digits = tools::digits<T>(); + + if(digits != current_digits) + { + // Oh my... our precision has changed! + table = std::vector<std::vector<T> >(1, std::vector<T>(1, T(-1))); + digits = current_digits; + } + int index = n - 1; if(index >= (int)table.size()) @@ -534,5 +550,9 @@ } } } // namespace boost::math::detail +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // _BOOST_POLYGAMMA_DETAIL_2013_07_30_HPP_ diff --git a/boost/math/special_functions/detail/unchecked_factorial.hpp b/boost/math/special_functions/detail/unchecked_factorial.hpp index 3c23d6e15a..ca9a752291 100644 --- a/boost/math/special_functions/detail/unchecked_factorial.hpp +++ b/boost/math/special_functions/detail/unchecked_factorial.hpp @@ -474,7 +474,29 @@ struct max_factorial<double> #ifndef BOOST_MATH_NO_LEXICAL_CAST template <class T> -inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) +struct unchecked_factorial_initializer +{ + struct init + { + init() + { + boost::math::unchecked_factorial<T>(3); + } + void force_instantiate()const {} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T> +const typename unchecked_factorial_initializer<T>::init unchecked_factorial_initializer<T>::initializer; + + +template <class T, int N> +inline T unchecked_factorial_imp(unsigned i, const mpl::int_<N>&) { BOOST_STATIC_ASSERT(!boost::is_integral<T>::value); // factorial<unsigned int>(n) is not implemented @@ -484,114 +506,255 @@ inline T unchecked_factorial(unsigned i BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n)); // See factorial documentation for more detail. + unchecked_factorial_initializer<T>::force_instantiate(); + static const boost::array<T, 101> factorials = {{ - boost::lexical_cast<T>("1"), - boost::lexical_cast<T>("1"), - boost::lexical_cast<T>("2"), - boost::lexical_cast<T>("6"), - boost::lexical_cast<T>("24"), - boost::lexical_cast<T>("120"), - boost::lexical_cast<T>("720"), - boost::lexical_cast<T>("5040"), - boost::lexical_cast<T>("40320"), - boost::lexical_cast<T>("362880"), - boost::lexical_cast<T>("3628800"), - boost::lexical_cast<T>("39916800"), - boost::lexical_cast<T>("479001600"), - boost::lexical_cast<T>("6227020800"), - boost::lexical_cast<T>("87178291200"), - boost::lexical_cast<T>("1307674368000"), - boost::lexical_cast<T>("20922789888000"), - boost::lexical_cast<T>("355687428096000"), - boost::lexical_cast<T>("6402373705728000"), - boost::lexical_cast<T>("121645100408832000"), - boost::lexical_cast<T>("2432902008176640000"), - boost::lexical_cast<T>("51090942171709440000"), - boost::lexical_cast<T>("1124000727777607680000"), - boost::lexical_cast<T>("25852016738884976640000"), - boost::lexical_cast<T>("620448401733239439360000"), - boost::lexical_cast<T>("15511210043330985984000000"), - boost::lexical_cast<T>("403291461126605635584000000"), - boost::lexical_cast<T>("10888869450418352160768000000"), - boost::lexical_cast<T>("304888344611713860501504000000"), - boost::lexical_cast<T>("8841761993739701954543616000000"), - boost::lexical_cast<T>("265252859812191058636308480000000"), - boost::lexical_cast<T>("8222838654177922817725562880000000"), - boost::lexical_cast<T>("263130836933693530167218012160000000"), - boost::lexical_cast<T>("8683317618811886495518194401280000000"), - boost::lexical_cast<T>("295232799039604140847618609643520000000"), - boost::lexical_cast<T>("10333147966386144929666651337523200000000"), - boost::lexical_cast<T>("371993326789901217467999448150835200000000"), - boost::lexical_cast<T>("13763753091226345046315979581580902400000000"), - boost::lexical_cast<T>("523022617466601111760007224100074291200000000"), - boost::lexical_cast<T>("20397882081197443358640281739902897356800000000"), - boost::lexical_cast<T>("815915283247897734345611269596115894272000000000"), - boost::lexical_cast<T>("33452526613163807108170062053440751665152000000000"), - boost::lexical_cast<T>("1405006117752879898543142606244511569936384000000000"), - boost::lexical_cast<T>("60415263063373835637355132068513997507264512000000000"), - boost::lexical_cast<T>("2658271574788448768043625811014615890319638528000000000"), - boost::lexical_cast<T>("119622220865480194561963161495657715064383733760000000000"), - boost::lexical_cast<T>("5502622159812088949850305428800254892961651752960000000000"), - boost::lexical_cast<T>("258623241511168180642964355153611979969197632389120000000000"), - boost::lexical_cast<T>("12413915592536072670862289047373375038521486354677760000000000"), - boost::lexical_cast<T>("608281864034267560872252163321295376887552831379210240000000000"), - boost::lexical_cast<T>("30414093201713378043612608166064768844377641568960512000000000000"), - boost::lexical_cast<T>("1551118753287382280224243016469303211063259720016986112000000000000"), - boost::lexical_cast<T>("80658175170943878571660636856403766975289505440883277824000000000000"), - boost::lexical_cast<T>("4274883284060025564298013753389399649690343788366813724672000000000000"), - boost::lexical_cast<T>("230843697339241380472092742683027581083278564571807941132288000000000000"), - boost::lexical_cast<T>("12696403353658275925965100847566516959580321051449436762275840000000000000"), - boost::lexical_cast<T>("710998587804863451854045647463724949736497978881168458687447040000000000000"), - boost::lexical_cast<T>("40526919504877216755680601905432322134980384796226602145184481280000000000000"), - boost::lexical_cast<T>("2350561331282878571829474910515074683828862318181142924420699914240000000000000"), - boost::lexical_cast<T>("138683118545689835737939019720389406345902876772687432540821294940160000000000000"), - boost::lexical_cast<T>("8320987112741390144276341183223364380754172606361245952449277696409600000000000000"), - boost::lexical_cast<T>("507580213877224798800856812176625227226004528988036003099405939480985600000000000000"), - boost::lexical_cast<T>("31469973260387937525653122354950764088012280797258232192163168247821107200000000000000"), - boost::lexical_cast<T>("1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000"), - boost::lexical_cast<T>("126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000"), - boost::lexical_cast<T>("8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000"), - boost::lexical_cast<T>("544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000"), - boost::lexical_cast<T>("36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000"), - boost::lexical_cast<T>("2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000"), - boost::lexical_cast<T>("171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000"), - boost::lexical_cast<T>("11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000"), - boost::lexical_cast<T>("850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000"), - boost::lexical_cast<T>("61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000"), - boost::lexical_cast<T>("4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000"), - boost::lexical_cast<T>("330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000"), - boost::lexical_cast<T>("24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000"), - boost::lexical_cast<T>("1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000"), - boost::lexical_cast<T>("145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000"), - boost::lexical_cast<T>("11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000"), - boost::lexical_cast<T>("894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000"), - boost::lexical_cast<T>("71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000"), - boost::lexical_cast<T>("5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000"), - boost::lexical_cast<T>("475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000"), - boost::lexical_cast<T>("39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000"), - boost::lexical_cast<T>("3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000"), - boost::lexical_cast<T>("281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000"), - boost::lexical_cast<T>("24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000"), - boost::lexical_cast<T>("2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000"), - boost::lexical_cast<T>("185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000"), - boost::lexical_cast<T>("16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000"), - boost::lexical_cast<T>("1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000"), - boost::lexical_cast<T>("135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000"), - boost::lexical_cast<T>("12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000"), - boost::lexical_cast<T>("1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000"), - boost::lexical_cast<T>("108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000"), - boost::lexical_cast<T>("10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000"), - boost::lexical_cast<T>("991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000"), - boost::lexical_cast<T>("96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000"), - boost::lexical_cast<T>("9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000"), - boost::lexical_cast<T>("933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000"), - boost::lexical_cast<T>("93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000"), + T(boost::math::tools::convert_from_string<T>("1")), + T(boost::math::tools::convert_from_string<T>("1")), + T(boost::math::tools::convert_from_string<T>("2")), + T(boost::math::tools::convert_from_string<T>("6")), + T(boost::math::tools::convert_from_string<T>("24")), + T(boost::math::tools::convert_from_string<T>("120")), + T(boost::math::tools::convert_from_string<T>("720")), + T(boost::math::tools::convert_from_string<T>("5040")), + T(boost::math::tools::convert_from_string<T>("40320")), + T(boost::math::tools::convert_from_string<T>("362880")), + T(boost::math::tools::convert_from_string<T>("3628800")), + T(boost::math::tools::convert_from_string<T>("39916800")), + T(boost::math::tools::convert_from_string<T>("479001600")), + T(boost::math::tools::convert_from_string<T>("6227020800")), + T(boost::math::tools::convert_from_string<T>("87178291200")), + T(boost::math::tools::convert_from_string<T>("1307674368000")), + T(boost::math::tools::convert_from_string<T>("20922789888000")), + T(boost::math::tools::convert_from_string<T>("355687428096000")), + T(boost::math::tools::convert_from_string<T>("6402373705728000")), + T(boost::math::tools::convert_from_string<T>("121645100408832000")), + T(boost::math::tools::convert_from_string<T>("2432902008176640000")), + T(boost::math::tools::convert_from_string<T>("51090942171709440000")), + T(boost::math::tools::convert_from_string<T>("1124000727777607680000")), + T(boost::math::tools::convert_from_string<T>("25852016738884976640000")), + T(boost::math::tools::convert_from_string<T>("620448401733239439360000")), + T(boost::math::tools::convert_from_string<T>("15511210043330985984000000")), + T(boost::math::tools::convert_from_string<T>("403291461126605635584000000")), + T(boost::math::tools::convert_from_string<T>("10888869450418352160768000000")), + T(boost::math::tools::convert_from_string<T>("304888344611713860501504000000")), + T(boost::math::tools::convert_from_string<T>("8841761993739701954543616000000")), + T(boost::math::tools::convert_from_string<T>("265252859812191058636308480000000")), + T(boost::math::tools::convert_from_string<T>("8222838654177922817725562880000000")), + T(boost::math::tools::convert_from_string<T>("263130836933693530167218012160000000")), + T(boost::math::tools::convert_from_string<T>("8683317618811886495518194401280000000")), + T(boost::math::tools::convert_from_string<T>("295232799039604140847618609643520000000")), + T(boost::math::tools::convert_from_string<T>("10333147966386144929666651337523200000000")), + T(boost::math::tools::convert_from_string<T>("371993326789901217467999448150835200000000")), + T(boost::math::tools::convert_from_string<T>("13763753091226345046315979581580902400000000")), + T(boost::math::tools::convert_from_string<T>("523022617466601111760007224100074291200000000")), + T(boost::math::tools::convert_from_string<T>("20397882081197443358640281739902897356800000000")), + T(boost::math::tools::convert_from_string<T>("815915283247897734345611269596115894272000000000")), + T(boost::math::tools::convert_from_string<T>("33452526613163807108170062053440751665152000000000")), + T(boost::math::tools::convert_from_string<T>("1405006117752879898543142606244511569936384000000000")), + T(boost::math::tools::convert_from_string<T>("60415263063373835637355132068513997507264512000000000")), + T(boost::math::tools::convert_from_string<T>("2658271574788448768043625811014615890319638528000000000")), + T(boost::math::tools::convert_from_string<T>("119622220865480194561963161495657715064383733760000000000")), + T(boost::math::tools::convert_from_string<T>("5502622159812088949850305428800254892961651752960000000000")), + T(boost::math::tools::convert_from_string<T>("258623241511168180642964355153611979969197632389120000000000")), + T(boost::math::tools::convert_from_string<T>("12413915592536072670862289047373375038521486354677760000000000")), + T(boost::math::tools::convert_from_string<T>("608281864034267560872252163321295376887552831379210240000000000")), + T(boost::math::tools::convert_from_string<T>("30414093201713378043612608166064768844377641568960512000000000000")), + T(boost::math::tools::convert_from_string<T>("1551118753287382280224243016469303211063259720016986112000000000000")), + T(boost::math::tools::convert_from_string<T>("80658175170943878571660636856403766975289505440883277824000000000000")), + T(boost::math::tools::convert_from_string<T>("4274883284060025564298013753389399649690343788366813724672000000000000")), + T(boost::math::tools::convert_from_string<T>("230843697339241380472092742683027581083278564571807941132288000000000000")), + T(boost::math::tools::convert_from_string<T>("12696403353658275925965100847566516959580321051449436762275840000000000000")), + T(boost::math::tools::convert_from_string<T>("710998587804863451854045647463724949736497978881168458687447040000000000000")), + T(boost::math::tools::convert_from_string<T>("40526919504877216755680601905432322134980384796226602145184481280000000000000")), + T(boost::math::tools::convert_from_string<T>("2350561331282878571829474910515074683828862318181142924420699914240000000000000")), + T(boost::math::tools::convert_from_string<T>("138683118545689835737939019720389406345902876772687432540821294940160000000000000")), + T(boost::math::tools::convert_from_string<T>("8320987112741390144276341183223364380754172606361245952449277696409600000000000000")), + T(boost::math::tools::convert_from_string<T>("507580213877224798800856812176625227226004528988036003099405939480985600000000000000")), + T(boost::math::tools::convert_from_string<T>("31469973260387937525653122354950764088012280797258232192163168247821107200000000000000")), + T(boost::math::tools::convert_from_string<T>("1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000")), + T(boost::math::tools::convert_from_string<T>("126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000")), + T(boost::math::tools::convert_from_string<T>("8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000")), + T(boost::math::tools::convert_from_string<T>("544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000")), + T(boost::math::tools::convert_from_string<T>("36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000")), + T(boost::math::tools::convert_from_string<T>("2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000")), + T(boost::math::tools::convert_from_string<T>("171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000")), + T(boost::math::tools::convert_from_string<T>("11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000")), + T(boost::math::tools::convert_from_string<T>("850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000")), + T(boost::math::tools::convert_from_string<T>("61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000")), + T(boost::math::tools::convert_from_string<T>("4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000")), + T(boost::math::tools::convert_from_string<T>("330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000")), + T(boost::math::tools::convert_from_string<T>("24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000")), + T(boost::math::tools::convert_from_string<T>("1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000")), + T(boost::math::tools::convert_from_string<T>("145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000")), + T(boost::math::tools::convert_from_string<T>("11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000")), + T(boost::math::tools::convert_from_string<T>("894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000")), + T(boost::math::tools::convert_from_string<T>("71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000")), + T(boost::math::tools::convert_from_string<T>("5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000")), + T(boost::math::tools::convert_from_string<T>("475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000")), + T(boost::math::tools::convert_from_string<T>("39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000")), + T(boost::math::tools::convert_from_string<T>("3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000")), + T(boost::math::tools::convert_from_string<T>("281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000")), + T(boost::math::tools::convert_from_string<T>("24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000")), + T(boost::math::tools::convert_from_string<T>("2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000")), + T(boost::math::tools::convert_from_string<T>("185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000")), + T(boost::math::tools::convert_from_string<T>("16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000")), + T(boost::math::tools::convert_from_string<T>("1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000")), + T(boost::math::tools::convert_from_string<T>("93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000")), }}; return factorials[i]; } template <class T> +inline T unchecked_factorial_imp(unsigned i, const mpl::int_<0>&) +{ + BOOST_STATIC_ASSERT(!boost::is_integral<T>::value); + // factorial<unsigned int>(n) is not implemented + // because it would overflow integral type T for too small n + // to be useful. Use instead a floating-point type, + // and convert to an unsigned type if essential, for example: + // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n)); + // See factorial documentation for more detail. +#ifdef BOOST_NO_CXX11_THREAD_LOCAL + unchecked_factorial_initializer<T>::force_instantiate(); +#endif + static const char* const factorial_strings[] = { + "1", + "1", + "2", + "6", + "24", + "120", + "720", + "5040", + "40320", + "362880", + "3628800", + "39916800", + "479001600", + "6227020800", + "87178291200", + "1307674368000", + "20922789888000", + "355687428096000", + "6402373705728000", + "121645100408832000", + "2432902008176640000", + "51090942171709440000", + "1124000727777607680000", + "25852016738884976640000", + "620448401733239439360000", + "15511210043330985984000000", + "403291461126605635584000000", + "10888869450418352160768000000", + "304888344611713860501504000000", + "8841761993739701954543616000000", + "265252859812191058636308480000000", + "8222838654177922817725562880000000", + "263130836933693530167218012160000000", + "8683317618811886495518194401280000000", + "295232799039604140847618609643520000000", + "10333147966386144929666651337523200000000", + "371993326789901217467999448150835200000000", + "13763753091226345046315979581580902400000000", + "523022617466601111760007224100074291200000000", + "20397882081197443358640281739902897356800000000", + "815915283247897734345611269596115894272000000000", + "33452526613163807108170062053440751665152000000000", + "1405006117752879898543142606244511569936384000000000", + "60415263063373835637355132068513997507264512000000000", + "2658271574788448768043625811014615890319638528000000000", + "119622220865480194561963161495657715064383733760000000000", + "5502622159812088949850305428800254892961651752960000000000", + "258623241511168180642964355153611979969197632389120000000000", + "12413915592536072670862289047373375038521486354677760000000000", + "608281864034267560872252163321295376887552831379210240000000000", + "30414093201713378043612608166064768844377641568960512000000000000", + "1551118753287382280224243016469303211063259720016986112000000000000", + "80658175170943878571660636856403766975289505440883277824000000000000", + "4274883284060025564298013753389399649690343788366813724672000000000000", + "230843697339241380472092742683027581083278564571807941132288000000000000", + "12696403353658275925965100847566516959580321051449436762275840000000000000", + "710998587804863451854045647463724949736497978881168458687447040000000000000", + "40526919504877216755680601905432322134980384796226602145184481280000000000000", + "2350561331282878571829474910515074683828862318181142924420699914240000000000000", + "138683118545689835737939019720389406345902876772687432540821294940160000000000000", + "8320987112741390144276341183223364380754172606361245952449277696409600000000000000", + "507580213877224798800856812176625227226004528988036003099405939480985600000000000000", + "31469973260387937525653122354950764088012280797258232192163168247821107200000000000000", + "1982608315404440064116146708361898137544773690227268628106279599612729753600000000000000", + "126886932185884164103433389335161480802865516174545192198801894375214704230400000000000000", + "8247650592082470666723170306785496252186258551345437492922123134388955774976000000000000000", + "544344939077443064003729240247842752644293064388798874532860126869671081148416000000000000000", + "36471110918188685288249859096605464427167635314049524593701628500267962436943872000000000000000", + "2480035542436830599600990418569171581047399201355367672371710738018221445712183296000000000000000", + "171122452428141311372468338881272839092270544893520369393648040923257279754140647424000000000000000", + "11978571669969891796072783721689098736458938142546425857555362864628009582789845319680000000000000000", + "850478588567862317521167644239926010288584608120796235886430763388588680378079017697280000000000000000", + "61234458376886086861524070385274672740778091784697328983823014963978384987221689274204160000000000000000", + "4470115461512684340891257138125051110076800700282905015819080092370422104067183317016903680000000000000000", + "330788544151938641225953028221253782145683251820934971170611926835411235700971565459250872320000000000000000", + "24809140811395398091946477116594033660926243886570122837795894512655842677572867409443815424000000000000000000", + "1885494701666050254987932260861146558230394535379329335672487982961844043495537923117729972224000000000000000000", + "145183092028285869634070784086308284983740379224208358846781574688061991349156420080065207861248000000000000000000", + "11324281178206297831457521158732046228731749579488251990048962825668835325234200766245086213177344000000000000000000", + "894618213078297528685144171539831652069808216779571907213868063227837990693501860533361810841010176000000000000000000", + "71569457046263802294811533723186532165584657342365752577109445058227039255480148842668944867280814080000000000000000000", + "5797126020747367985879734231578109105412357244731625958745865049716390179693892056256184534249745940480000000000000000000", + "475364333701284174842138206989404946643813294067993328617160934076743994734899148613007131808479167119360000000000000000000", + "39455239697206586511897471180120610571436503407643446275224357528369751562996629334879591940103770870906880000000000000000000", + "3314240134565353266999387579130131288000666286242049487118846032383059131291716864129885722968716753156177920000000000000000000", + "281710411438055027694947944226061159480056634330574206405101912752560026159795933451040286452340924018275123200000000000000000000", + "24227095383672732381765523203441259715284870552429381750838764496720162249742450276789464634901319465571660595200000000000000000000", + "2107757298379527717213600518699389595229783738061356212322972511214654115727593174080683423236414793504734471782400000000000000000000", + "185482642257398439114796845645546284380220968949399346684421580986889562184028199319100141244804501828416633516851200000000000000000000", + "16507955160908461081216919262453619309839666236496541854913520707833171034378509739399912570787600662729080382999756800000000000000000000", + "1485715964481761497309522733620825737885569961284688766942216863704985393094065876545992131370884059645617234469978112000000000000000000000", + "135200152767840296255166568759495142147586866476906677791741734597153670771559994765685283954750449427751168336768008192000000000000000000000", + "12438414054641307255475324325873553077577991715875414356840239582938137710983519518443046123837041347353107486982656753664000000000000000000000", + "1156772507081641574759205162306240436214753229576413535186142281213246807121467315215203289516844845303838996289387078090752000000000000000000000", + "108736615665674308027365285256786601004186803580182872307497374434045199869417927630229109214583415458560865651202385340530688000000000000000000000", + "10329978488239059262599702099394727095397746340117372869212250571234293987594703124871765375385424468563282236864226607350415360000000000000000000000", + "991677934870949689209571401541893801158183648651267795444376054838492222809091499987689476037000748982075094738965754305639874560000000000000000000000", + "96192759682482119853328425949563698712343813919172976158104477319333745612481875498805879175589072651261284189679678167647067832320000000000000000000000", + "9426890448883247745626185743057242473809693764078951663494238777294707070023223798882976159207729119823605850588608460429412647567360000000000000000000000", + "933262154439441526816992388562667004907159682643816214685929638952175999932299156089414639761565182862536979208272237582511852109168640000000000000000000000", + "93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000", + }; + + static BOOST_MATH_THREAD_LOCAL T factorials[sizeof(factorial_strings) / sizeof(factorial_strings[0])]; + static BOOST_MATH_THREAD_LOCAL int digits = 0; + + int current_digits = boost::math::tools::digits<T>(); + + if(digits != current_digits) + { + digits = current_digits; + for(unsigned k = 0; k < sizeof(factorials) / sizeof(factorials[0]); ++k) + factorials[k] = static_cast<T>(boost::math::tools::convert_from_string<T>(factorial_strings[k])); + } + + return factorials[i]; +} + +template <class T> +inline T unchecked_factorial(unsigned i) +{ + typedef typename boost::math::policies::precision<T, boost::math::policies::policy<> >::type tag_type; + return unchecked_factorial_imp<T>(i, tag_type()); +} + +template <class T> struct max_factorial { BOOST_STATIC_CONSTANT(unsigned, value = 100); diff --git a/boost/math/special_functions/digamma.hpp b/boost/math/special_functions/digamma.hpp index 718eaf9278..f6431c67a6 100644 --- a/boost/math/special_functions/digamma.hpp +++ b/boost/math/special_functions/digamma.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/special_functions/math_fwd.hpp> @@ -489,7 +491,7 @@ T digamma_imp(T x, const mpl::int_<0>* t, const Policy& pol) // limit = 10 at 50 bit precision and // limit = 250 at 1000 bit precision. // - T lim = 10 + (tools::digits<T>() - 50) * 240 / 950; + int lim = 10 + ((tools::digits<T>() - 50) * 240L) / 950; T two_x = ldexp(x, 1); if(x >= lim) { @@ -624,5 +626,10 @@ inline typename tools::promote_args<T>::type } // namespace math } // namespace boost + +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif diff --git a/boost/math/special_functions/ellint_1.hpp b/boost/math/special_functions/ellint_1.hpp index 62a0bf3fde..d1d9d72e30 100644 --- a/boost/math/special_functions/ellint_1.hpp +++ b/boost/math/special_functions/ellint_1.hpp @@ -106,7 +106,6 @@ T ellint_f_imp(T phi, T k, const Policy& pol) sinp *= sinp; T cosp = cos(rphi); cosp *= cosp; - T c = 1 / sinp; BOOST_MATH_INSTRUMENT_VARIABLE(sinp); BOOST_MATH_INSTRUMENT_VARIABLE(cosp); if(sinp > tools::min_value<T>()) @@ -115,6 +114,7 @@ T ellint_f_imp(T phi, T k, const Policy& pol) // Use http://dlmf.nist.gov/19.25#E5, note that // c-1 simplifies to cot^2(rphi) which avoid cancellation: // + T c = 1 / sinp; result = rphi == 0 ? static_cast<T>(0) : static_cast<T>(s * ellint_rf_imp(T(cosp / sinp), T(c - k * k), c, pol)); } else diff --git a/boost/math/special_functions/ellint_2.hpp b/boost/math/special_functions/ellint_2.hpp index 7c1eab45c2..9ee6b63821 100644 --- a/boost/math/special_functions/ellint_2.hpp +++ b/boost/math/special_functions/ellint_2.hpp @@ -94,10 +94,6 @@ T ellint_e_imp(T phi, T k, const Policy& pol) s = -1; rphi = constants::half_pi<T>() - rphi; } - T sinp = sin(rphi); - T cosp = cos(rphi); - T c = 1 / (sinp * sinp); - T cm1 = cosp * cosp / (sinp * sinp); // c - 1 T k2 = k * k; if(k2 > 1) { @@ -111,6 +107,10 @@ T ellint_e_imp(T phi, T k, const Policy& pol) else { // http://dlmf.nist.gov/19.25#E10 + T sinp = sin(rphi); + T cosp = cos(rphi); + T c = 1 / (sinp * sinp); + T cm1 = cosp * cosp / (sinp * sinp); // c - 1 result = s * ((1 - k2) * ellint_rf_imp(cm1, T(c - k2), c, pol) + k2 * (1 - k2) * ellint_rd(cm1, c, T(c - k2), pol) / 3 + k2 * sqrt(cm1 / (c * (c - k2)))); } if(m != 0) diff --git a/boost/math/special_functions/ellint_3.hpp b/boost/math/special_functions/ellint_3.hpp index b6670a196f..b8b36729cf 100644 --- a/boost/math/special_functions/ellint_3.hpp +++ b/boost/math/special_functions/ellint_3.hpp @@ -58,19 +58,19 @@ T ellint_pi_imp(T v, T phi, T k, T vc, const Policy& pol) T sphi = sin(fabs(phi)); T result = 0; - if(v > 1 / (sphi * sphi)) - { - // Complex result is a domain error: - return policies::raise_domain_error<T>(function, - "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol); - } - // Special cases first: if(v == 0) { // A&S 17.7.18 & 19 return (k == 0) ? phi : ellint_f_imp(phi, k, pol); } + if((v > 0) && (1 / v < (sphi * sphi))) + { + // Complex result is a domain error: + return policies::raise_domain_error<T>(function, + "Got v = %1%, but result is complex for v > 1 / sin^2(phi)", v, pol); + } + if(v == 1) { // http://functions.wolfram.com/08.06.03.0008.01 diff --git a/boost/math/special_functions/expint.hpp b/boost/math/special_functions/expint.hpp index c26420db9e..fc656f3355 100644 --- a/boost/math/special_functions/expint.hpp +++ b/boost/math/special_functions/expint.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/tools/precision.hpp> @@ -1667,6 +1669,10 @@ inline typename tools::promote_args<T>::type }} // namespaces +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_EXPINT_HPP diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp index 3a3191a807..a2b30f8551 100644 --- a/boost/math/special_functions/gamma.hpp +++ b/boost/math/special_functions/gamma.hpp @@ -33,6 +33,7 @@ #include <boost/math/special_functions/detail/unchecked_factorial.hpp> #include <boost/math/special_functions/detail/lgamma_small.hpp> #include <boost/math/special_functions/bernoulli.hpp> +#include <boost/math/special_functions/zeta.hpp> #include <boost/type_traits/is_convertible.hpp> #include <boost/assert.hpp> #include <boost/mpl/greater.hpp> @@ -495,6 +496,37 @@ T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) } template <class T, class Policy> +inline T log_gamma_near_1(const T& z, Policy const& pol) +{ + // + // This is for the multiprecision case where there is + // no lanczos support... + // + BOOST_MATH_STD_USING // ADL of std names + + BOOST_ASSERT(fabs(z) < 1); + + T result = -constants::euler<T>() * z; + + T power_term = z * z; + T term; + unsigned j = 0; + + do + { + term = boost::math::zeta<T>(j + 2, pol) * power_term / (j + 2); + if(j & 1) + result -= term; + else + result += term; + power_term *= z; + ++j; + } while(fabs(result) * tools::epsilon<T>() < fabs(term)); + + return result; +} + +template <class T, class Policy> T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign) { BOOST_MATH_STD_USING @@ -529,7 +561,15 @@ T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sig // Here we simply take the logarithm of tgamma(). This is somewhat // inefficient, but simple. The rationale is that the argument here // is relatively small and overflow is not expected to be likely. - if (z > -tools::root_epsilon<T>()) + if(fabs(z - 1) < 0.25) + { + return log_gamma_near_1(T(zz - 1), pol); + } + else if(fabs(z - 2) < 0.25) + { + return log_gamma_near_1(T(zz - 2), pol) + log(zz - 1); + } + else if (z > -tools::root_epsilon<T>()) { // Reflection formula may fail if z is very close to zero, let the series // expansion for tgamma close to zero do the work: @@ -634,7 +674,7 @@ T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l) mpl::greater<precision_type, mpl::int_<113> > >, typename mpl::if_< - is_same<Lanczos, lanczos::lanczos24m113>, + mpl::and_<is_same<Lanczos, lanczos::lanczos24m113>, mpl::greater<precision_type, mpl::int_<0> > >, mpl::int_<113>, mpl::int_<0> >::type, @@ -680,32 +720,16 @@ T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l) } template <class T, class Policy> -inline T tgammap1m1_imp(T dz, Policy const& pol, - const ::boost::math::lanczos::undefined_lanczos& l) +inline T tgammap1m1_imp(T z, Policy const& pol, + const ::boost::math::lanczos::undefined_lanczos&) { BOOST_MATH_STD_USING // ADL of std names - // - // There should be a better solution than this, but the - // algebra isn't easy for the general case.... - // Start by subracting 1 from tgamma: - // - T result = gamma_imp(T(1 + dz), pol, l) - 1; - BOOST_MATH_INSTRUMENT_CODE(result); - // - // Test the level of cancellation error observed: we loose one bit - // for each power of 2 the result is less than 1. If we would get - // more bits from our most precise lgamma rational approximation, - // then use that instead: - // - BOOST_MATH_INSTRUMENT_CODE((dz > -0.5)); - BOOST_MATH_INSTRUMENT_CODE((dz < 2)); - BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34)); - if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34)) + + if(fabs(z) < 0.55) { - result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113()); - BOOST_MATH_INSTRUMENT_CODE(result); + return boost::math::expm1(log_gamma_near_1(z, pol)); } - return result; + return boost::math::expm1(boost::math::lgamma(1 + z, pol)); } // @@ -1396,16 +1420,26 @@ T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& } T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>()); T result; - if(fabs(delta) < 10) + if(z + delta == z) { - result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol)); + if(fabs(delta) < 10) + result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol)); + else + result = 1; } else { - result = pow(zgh / (zgh + delta), z - constants::half<T>()); + if(fabs(delta) < 10) + { + result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol)); + } + else + { + result = pow(zgh / (zgh + delta), z - constants::half<T>()); + } + // Split the calculation up to avoid spurious overflow: + result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta)); } - // Split the calculation up to avoid spurious overflow: - result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta)); result *= pow(constants::e<T>() / (zgh + delta), delta); return result; } @@ -1663,7 +1697,14 @@ struct igamma_initializer template <int N> static void do_init(const mpl::int_<N>&) { - boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy()); + // If std::numeric_limits<T>::digits is zero, we must not call + // our inituialization code here as the precision presumably + // varies at runtime, and will not have been set yet. Plus the + // code requiring initialization isn't called when digits == 0. + if(std::numeric_limits<T>::digits) + { + boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy()); + } } static void do_init(const mpl::int_<53>&){} void force_instantiate()const{} diff --git a/boost/math/special_functions/lanczos.hpp b/boost/math/special_functions/lanczos.hpp index 0db21d3d16..c9f4f90dcf 100644 --- a/boost/math/special_functions/lanczos.hpp +++ b/boost/math/special_functions/lanczos.hpp @@ -1284,7 +1284,7 @@ struct lanczos } // namespace boost #if !defined(_CRAYC) && !defined(__CUDACC__) && (!defined(__GNUC__) || (__GNUC__ > 3) || ((__GNUC__ == 3) && (__GNUC_MINOR__ > 3))) -#if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) +#if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(__SSE2__) || defined(_M_AMD64) || defined(_M_X64) #include <boost/math/special_functions/detail/lanczos_sse2.hpp> #endif #endif diff --git a/boost/math/special_functions/log1p.hpp b/boost/math/special_functions/log1p.hpp index 62f5b8027c..50bdaeb385 100644 --- a/boost/math/special_functions/log1p.hpp +++ b/boost/math/special_functions/log1p.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/config/no_tr1/cmath.hpp> @@ -497,6 +499,10 @@ inline typename tools::promote_args<T>::type log1pmx(T x) } // namespace math } // namespace boost +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_LOG1P_INCLUDED diff --git a/boost/math/special_functions/powm1.hpp b/boost/math/special_functions/powm1.hpp index fe2fce35d8..37db8746c8 100644 --- a/boost/math/special_functions/powm1.hpp +++ b/boost/math/special_functions/powm1.hpp @@ -8,6 +8,8 @@ #ifdef _MSC_VER #pragma once +#pragma warning(push) +#pragma warning(disable:4702) // Unreachable code (release mode only warning) #endif #include <boost/math/special_functions/math_fwd.hpp> @@ -70,6 +72,10 @@ inline typename tools::promote_args<T1, T2>::type } // namespace math } // namespace boost +#ifdef _MSC_VER +#pragma warning(pop) +#endif + #endif // BOOST_MATH_POWM1 diff --git a/boost/math/special_functions/zeta.hpp b/boost/math/special_functions/zeta.hpp index 616bb0cccc..91b83c1b2e 100644 --- a/boost/math/special_functions/zeta.hpp +++ b/boost/math/special_functions/zeta.hpp @@ -878,11 +878,19 @@ T zeta_imp_odd_integer(int s, const T&, const Policy&, const mpl::true_&) template <class T, class Policy> T zeta_imp_odd_integer(int s, const T& sc, const Policy& pol, const mpl::false_&) { - static bool is_init = false; - static T results[50] = {}; + static BOOST_MATH_THREAD_LOCAL bool is_init = false; + static BOOST_MATH_THREAD_LOCAL T results[50] = {}; + static BOOST_MATH_THREAD_LOCAL int digits = tools::digits<T>(); + int current_digits = tools::digits<T>(); + if(digits != current_digits) + { + // Oh my precision has changed... + is_init = false; + } if(!is_init) { is_init = true; + digits = current_digits; for(unsigned k = 0; k < sizeof(results) / sizeof(results[0]); ++k) { T arg = k * 2 + 3; |