diff options
Diffstat (limited to 'boost/math')
51 files changed, 4489 insertions, 135 deletions
diff --git a/boost/math/bindings/detail/big_lanczos.hpp b/boost/math/bindings/detail/big_lanczos.hpp index a11f1b8196..314b4f88c4 100644 --- a/boost/math/bindings/detail/big_lanczos.hpp +++ b/boost/math/bindings/detail/big_lanczos.hpp @@ -32,6 +32,7 @@ struct lanczos22UDT : public mpl::int_<120> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos22UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[22] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 46198410803245094237463011094.12173081986)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 43735859291852324413622037436.321513777)), @@ -86,6 +87,7 @@ struct lanczos22UDT : public mpl::int_<120> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos22UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[22] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 6939996264376682180.277485395074954356211)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 6570067992110214451.87201438870245659384)), @@ -141,6 +143,7 @@ struct lanczos22UDT : public mpl::int_<120> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos22UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[21] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 8.318998691953337183034781139546384476554)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, -63.15415991415959158214140353299240638675)), @@ -219,6 +222,7 @@ struct lanczos31UDT template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos31UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[31] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 0.2579646553333513328235723061836959833277e46)), static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 0.2444796504337453845497419271639377138264e46)), @@ -291,6 +295,7 @@ struct lanczos31UDT template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos31UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[31] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 30137154810677525966583148469478.52374216)), static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 28561746428637727032849890123131.36314653)), @@ -364,6 +369,7 @@ struct lanczos31UDT template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos31UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[30] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 11.80038544942943603508206880307972596807)), static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, -130.6355975335626214564236363322099481079)), @@ -407,6 +413,7 @@ struct lanczos31UDT template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos31UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[30] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 147.9979641587472136175636384176549713358)), static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, -1638.404318611773924210055619836375434296)), @@ -461,6 +468,7 @@ struct lanczos61UDT template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos61UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() using namespace boost; static const T d[61] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 2.50662827463100050241576528481104525300698674060993831662992357634229365460784197494659584)), @@ -536,6 +544,7 @@ struct lanczos61UDT template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos61UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() using namespace boost; static const T d[61] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 0.901751806425638853077358552989167785490911341809902155556127108480303870921448984935411583e-27)), @@ -611,6 +620,7 @@ struct lanczos61UDT template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos61UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() using namespace boost; static const T d[60] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 23.2463658527729692390378860713647146932236940604550445351214987229819352880524561852919518)), @@ -685,6 +695,7 @@ struct lanczos61UDT template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos61UDT, T>::force_instantiate(); // Ensure our constants get initialized before main() using namespace boost; static const T d[60] = { static_cast<T>(BOOST_MATH_HUGE_CONSTANT(T, 150, 557.56438192770795764344217888434355281097193198928944200046501607026919782564033547346298)), diff --git a/boost/math/bindings/e_float.hpp b/boost/math/bindings/e_float.hpp index 959bff109b..d7f9a73508 100644 --- a/boost/math/bindings/e_float.hpp +++ b/boost/math/bindings/e_float.hpp @@ -318,6 +318,11 @@ inline e_float atan(const e_float& v) return ::ef::atan(v.value()); } +inline e_float atan2(const e_float& v, const e_float& u) +{ + return ::ef::atan2(v.value(), u.value()); +} + inline e_float ldexp(const e_float& v, int e) { return v.value() * ::ef::pow2(e); diff --git a/boost/math/bindings/mpfr.hpp b/boost/math/bindings/mpfr.hpp index 95be0ee2ba..e5c6ba071d 100644 --- a/boost/math/bindings/mpfr.hpp +++ b/boost/math/bindings/mpfr.hpp @@ -860,7 +860,11 @@ inline mpfr_class bessel_i1(mpfr_class x) } // namespace detail -}} +} + +template<> struct is_convertible<long double, mpfr_class> : public mpl::false_{}; + +} #endif // BOOST_MATH_MPLFR_BINDINGS_HPP diff --git a/boost/math/bindings/rr.hpp b/boost/math/bindings/rr.hpp index 58c22b8891..6ec79f953d 100644 --- a/boost/math/bindings/rr.hpp +++ b/boost/math/bindings/rr.hpp @@ -763,6 +763,17 @@ namespace ntl{ NTL::RR::precision()); } + inline RR atan2(RR y, RR x) + { + if(x > 0) + return atan(y / x); + if(x < 0) + { + return y < 0 ? atan(y / x) - boost::math::constants::pi<RR>() : atan(y / x) + boost::math::constants::pi<RR>(); + } + return y < 0 ? -boost::math::constants::half_pi<RR>() : boost::math::constants::half_pi<RR>() ; + } + inline RR sinh(RR z) { return (expm1(z.value()) - expm1(-z.value())) / 2; diff --git a/boost/math/concepts/real_concept.hpp b/boost/math/concepts/real_concept.hpp index d9000754c2..1ed2c1df00 100644 --- a/boost/math/concepts/real_concept.hpp +++ b/boost/math/concepts/real_concept.hpp @@ -28,6 +28,7 @@ #include <boost/math/special_functions/round.hpp> #include <boost/math/special_functions/trunc.hpp> #include <boost/math/special_functions/modf.hpp> +#include <boost/math/tools/big_constant.hpp> #include <boost/math/tools/precision.hpp> #include <boost/math/policies/policy.hpp> #if defined(__SGI_STL_PORT) @@ -324,6 +325,12 @@ namespace tools { template <> +inline concepts::real_concept make_big_value<concepts::real_concept>(long double val, const char* , mpl::false_ const&, mpl::false_ const&) +{ + return val; // Can't use lexical_cast here, sometimes it fails.... +} + +template <> inline concepts::real_concept max_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept)) { return max_value<concepts::real_concept_base_type>(); diff --git a/boost/math/concepts/real_type_concept.hpp b/boost/math/concepts/real_type_concept.hpp index f6ec016211..97da780129 100644 --- a/boost/math/concepts/real_type_concept.hpp +++ b/boost/math/concepts/real_type_concept.hpp @@ -6,6 +6,7 @@ #ifndef BOOST_MATH_REAL_TYPE_CONCEPT_HPP #define BOOST_MATH_REAL_TYPE_CONCEPT_HPP +#include <boost/config.hpp> #ifdef BOOST_MSVC #pragma warning(push) #pragma warning(disable: 4100) diff --git a/boost/math/concepts/std_real_concept.hpp b/boost/math/concepts/std_real_concept.hpp index 645579573b..4c4eb6ac32 100644 --- a/boost/math/concepts/std_real_concept.hpp +++ b/boost/math/concepts/std_real_concept.hpp @@ -323,12 +323,19 @@ inline std::basic_istream<charT, traits>& operator>>(std::basic_istream<charT, t }} #include <boost/math/tools/precision.hpp> +#include <boost/math/tools/big_constant.hpp> namespace boost{ namespace math{ namespace tools { template <> +inline concepts::std_real_concept make_big_value<concepts::std_real_concept>(long double val, const char* , mpl::false_ const&, mpl::false_ const&) +{ + return val; // Can't use lexical_cast here, sometimes it fails.... +} + +template <> inline concepts::std_real_concept max_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept)) { return max_value<concepts::std_real_concept_base_type>(); diff --git a/boost/math/constants/calculate_constants.hpp b/boost/math/constants/calculate_constants.hpp new file mode 100644 index 0000000000..0b78929e71 --- /dev/null +++ b/boost/math/constants/calculate_constants.hpp @@ -0,0 +1,959 @@ +// Copyright John Maddock 2010, 2012. +// Copyright Paul A. Bristow 2011, 2012. + +// 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_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED +#define BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED + +#include <boost/math/special_functions/trunc.hpp> + +namespace boost{ namespace math{ namespace constants{ namespace detail{ + +template <class T> +template<int N> +inline T constant_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + + return ldexp(acos(T(0)), 1); + + /* + // Although this code works well, it's usually more accurate to just call acos + // and access the number types own representation of PI which is usually calculated + // at slightly higher precision... + + T result; + T a = 1; + T b; + T A(a); + T B = 0.5f; + T D = 0.25f; + + T lim; + lim = boost::math::tools::epsilon<T>(); + + unsigned k = 1; + + do + { + result = A + B; + result = ldexp(result, -2); + b = sqrt(B); + a += b; + a = ldexp(a, -1); + A = a * a; + B = A - result; + B = ldexp(B, 1); + result = A - B; + bool neg = boost::math::sign(result) < 0; + if(neg) + result = -result; + if(result <= lim) + break; + if(neg) + result = -result; + result = ldexp(result, k - 1); + D -= result; + ++k; + lim = ldexp(lim, 1); + } + while(true); + + result = B / D; + return result; + */ +} + +template <class T> +template<int N> +inline T constant_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return 2 * pi<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> // 2 / pi +template<int N> +inline T constant_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return 2 / pi<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> // sqrt(2/pi) +template <int N> +inline T constant_root_two_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return sqrt((2 / pi<T, policies::policy<policies::digits2<N> > >())); +} + +template <class T> +template<int N> +inline T constant_one_div_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return 1 / two_pi<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> +template<int N> +inline T constant_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(pi<T, policies::policy<policies::digits2<N> > >()); +} + +template <class T> +template<int N> +inline T constant_root_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(pi<T, policies::policy<policies::digits2<N> > >() / 2); +} + +template <class T> +template<int N> +inline T constant_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(two_pi<T, policies::policy<policies::digits2<N> > >()); +} + +template <class T> +template<int N> +inline T constant_root_ln_four<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(log(static_cast<T>(4))); +} + +template <class T> +template<int N> +inline T constant_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + // + // Although we can clearly calculate this from first principles, this hooks into + // T's own notion of e, which hopefully will more accurate than one calculated to + // a few epsilon: + // + BOOST_MATH_STD_USING + return exp(static_cast<T>(1)); +} + +template <class T> +template<int N> +inline T constant_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return static_cast<T>(1) / static_cast<T>(2); +} + +template <class T> +template<int M> +inline T constant_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<M>)) +{ + BOOST_MATH_STD_USING + // + // This is the method described in: + // "Some New Algorithms for High-Precision Computation of Euler's Constant" + // Richard P Brent and Edwin M McMillan. + // Mathematics of Comnputation, Volume 34, Number 149, Jan 1980, pages 305-312. + // See equation 17 with p = 2. + // + T n = 3 + (M ? (std::min)(M, tools::digits<T>()) : tools::digits<T>()) / 4; + T lim = M ? ldexp(T(1), (std::min)(M, tools::digits<T>())) : tools::epsilon<T>(); + T lnn = log(n); + T term = 1; + T N = -lnn; + T D = 1; + T Hk = 0; + T one = 1; + + for(unsigned k = 1;; ++k) + { + term *= n * n; + term /= k * k; + Hk += one / k; + N += term * (Hk - lnn); + D += term; + + if(term < D * lim) + break; + } + return N / D; +} + +template <class T> +template<int N> +inline T constant_euler_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return euler<T, policies::policy<policies::digits2<N> > >() + * euler<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> +template<int N> +inline T constant_one_div_euler<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(1) + / euler<T, policies::policy<policies::digits2<N> > >(); +} + + +template <class T> +template<int N> +inline T constant_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(static_cast<T>(2)); +} + + +template <class T> +template<int N> +inline T constant_root_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(static_cast<T>(3)); +} + +template <class T> +template<int N> +inline T constant_half_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(static_cast<T>(2)) / 2; +} + +template <class T> +template<int N> +inline T constant_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + // + // Although there are good ways to calculate this from scratch, this hooks into + // T's own notion of log(2) which will hopefully be accurate to the full precision + // of T: + // + BOOST_MATH_STD_USING + return log(static_cast<T>(2)); +} + +template <class T> +template<int N> +inline T constant_ln_ten<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return log(static_cast<T>(10)); +} + +template <class T> +template<int N> +inline T constant_ln_ln_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return log(log(static_cast<T>(2))); +} + +template <class T> +template<int N> +inline T constant_third<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(1) / static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_twothirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(2) / static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_two_thirds<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(2) / static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_three_quarters<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(3) / static_cast<T>(4); +} + +template <class T> +template<int N> +inline T constant_pi_minus_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_four_minus_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return static_cast<T>(4) - pi<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> +template<int N> +inline T constant_pow23_four_minus_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pow(four_minus_pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1.5)); +} + +template <class T> +template<int N> +inline T constant_exp_minus_half<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return exp(static_cast<T>(-0.5)); +} + +// Pi +template <class T> +template<int N> +inline T constant_one_div_root_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return static_cast<T>(1) / root_two<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> +template<int N> +inline T constant_one_div_root_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return static_cast<T>(1) / root_pi<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> +template<int N> +inline T constant_one_div_root_two_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + return static_cast<T>(1) / root_two_pi<T, policies::policy<policies::digits2<N> > >(); +} + +template <class T> +template<int N> +inline T constant_root_one_div_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(static_cast<T>(1) / pi<T, policies::policy<policies::digits2<N> > >()); +} + + +template <class T> +template<int N> +inline T constant_four_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(4) / static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_half_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(2); +} + + +template <class T> +template<int N> +inline T constant_third_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_sixth_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() / static_cast<T>(6); +} + +template <class T> +template<int N> +inline T constant_two_thirds_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(2) / static_cast<T>(3); +} + +template <class T> +template<int N> +inline T constant_three_quarters_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() * static_cast<T>(3) / static_cast<T>(4); +} + +template <class T> +template<int N> +inline T constant_pi_pow_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pow(pi<T, policies::policy<policies::digits2<N> > >(), e<T, policies::policy<policies::digits2<N> > >()); // +} + +template <class T> +template<int N> +inline T constant_pi_sqr<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >() ; // +} + +template <class T> +template<int N> +inline T constant_pi_sqr_div_six<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >() + / static_cast<T>(6); // +} + + +template <class T> +template<int N> +inline T constant_pi_cubed<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >() + ; // +} + +template <class T> +template<int N> +inline T constant_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3)); +} + +template <class T> +template<int N> +inline T constant_one_div_cbrt_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(1) + / pow(pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(1)/ static_cast<T>(3)); +} + +// Euler's e + +template <class T> +template<int N> +inline T constant_e_pow_pi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pow(e<T, policies::policy<policies::digits2<N> > >(), pi<T, policies::policy<policies::digits2<N> > >()); // +} + +template <class T> +template<int N> +inline T constant_root_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sqrt(e<T, policies::policy<policies::digits2<N> > >()); +} + +template <class T> +template<int N> +inline T constant_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return log10(e<T, policies::policy<policies::digits2<N> > >()); +} + +template <class T> +template<int N> +inline T constant_one_div_log10_e<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(1) / + log10(e<T, policies::policy<policies::digits2<N> > >()); +} + +// Trigonometric + +template <class T> +template<int N> +inline T constant_degree<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return pi<T, policies::policy<policies::digits2<N> > >() + / static_cast<T>(180) + ; // +} + +template <class T> +template<int N> +inline T constant_radian<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(180) + / pi<T, policies::policy<policies::digits2<N> > >() + ; // +} + +template <class T> +template<int N> +inline T constant_sin_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sin(static_cast<T>(1)) ; // +} + +template <class T> +template<int N> +inline T constant_cos_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return cos(static_cast<T>(1)) ; // +} + +template <class T> +template<int N> +inline T constant_sinh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return sinh(static_cast<T>(1)) ; // +} + +template <class T> +template<int N> +inline T constant_cosh_one<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return cosh(static_cast<T>(1)) ; // +} + +template <class T> +template<int N> +inline T constant_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return (static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ; // +} + +template <class T> +template<int N> +inline T constant_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + //return log(phi<T, policies::policy<policies::digits2<N> > >()); // ??? + return log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ); +} +template <class T> +template<int N> +inline T constant_one_div_ln_phi<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + return static_cast<T>(1) / + log((static_cast<T>(1) + sqrt(static_cast<T>(5)) )/static_cast<T>(2) ); +} + +// Zeta + +template <class T> +template<int N> +inline T constant_zeta_two<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + BOOST_MATH_STD_USING + + return pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >() + /static_cast<T>(6); +} + +template <class T> +template<int N> +inline T constant_zeta_three<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + // http://mathworld.wolfram.com/AperysConstant.html + // http://en.wikipedia.org/wiki/Mathematical_constant + + // http://oeis.org/A002117/constant + //T zeta3("1.20205690315959428539973816151144999076" + // "4986292340498881792271555341838205786313" + // "09018645587360933525814619915"); + + //"1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915" A002117 + // 1.202056903159594285399738161511449990, 76498629234049888179227155534183820578631309018645587360933525814619915780, +00); + //"1.2020569031595942 double + // http://www.spaennare.se/SSPROG/ssnum.pdf // section 11, Algorithmfor Apery’s constant zeta(3). + // Programs to Calculate some Mathematical Constants to Large Precision, Document Version 1.50 + + // by Stefan Spannare September 19, 2007 + // zeta(3) = 1/64 * sum + BOOST_MATH_STD_USING + T n_fact=static_cast<T>(1); // build n! for n = 0. + T sum = static_cast<double>(77); // Start with n = 0 case. + // for n = 0, (77/1) /64 = 1.203125 + //double lim = std::numeric_limits<double>::epsilon(); + T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>(); + for(unsigned int n = 1; n < 40; ++n) + { // three to five decimal digits per term, so 40 should be plenty for 100 decimal digits. + //cout << "n = " << n << endl; + n_fact *= n; // n! + T n_fact_p10 = n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact * n_fact; // (n!)^10 + T num = ((205 * n * n) + (250 * n) + 77) * n_fact_p10; // 205n^2 + 250n + 77 + // int nn = (2 * n + 1); + // T d = factorial(nn); // inline factorial. + T d = 1; + for(unsigned int i = 1; i <= (n+n + 1); ++i) // (2n + 1) + { + d *= i; + } + T den = d * d * d * d * d; // [(2n+1)!]^5 + //cout << "den = " << den << endl; + T term = num/den; + if (n % 2 != 0) + { //term *= -1; + sum -= term; + } + else + { + sum += term; + } + //cout << "term = " << term << endl; + //cout << "sum/64 = " << sum/64 << endl; + if(abs(term) < lim) + { + break; + } + } + return sum / 64; +} + +template <class T> +template<int N> +inline T constant_catalan<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ // http://oeis.org/A006752/constant + //T c("0.915965594177219015054603514932384110774" + //"149374281672134266498119621763019776254769479356512926115106248574"); + + // 9.159655941772190150546035149323841107, 74149374281672134266498119621763019776254769479356512926115106248574422619, -01); + + // This is equation (entry) 31 from + // http://www-2.cs.cmu.edu/~adamchik/articles/catalan/catalan.htm + // See also http://www.mpfr.org/algorithms.pdf + BOOST_MATH_STD_USING + T k_fact = 1; + T tk_fact = 1; + T sum = 1; + T term; + T lim = N ? ldexp(T(1), 1 - (std::min)(N, tools::digits<T>())) : tools::epsilon<T>(); + + for(unsigned k = 1;; ++k) + { + k_fact *= k; + tk_fact *= (2 * k) * (2 * k - 1); + term = k_fact * k_fact / (tk_fact * (2 * k + 1) * (2 * k + 1)); + sum += term; + if(term < lim) + { + break; + } + } + return boost::math::constants::pi<T, boost::math::policies::policy<> >() + * log(2 + boost::math::constants::root_three<T, boost::math::policies::policy<> >()) + / 8 + + 3 * sum / 8; +} + +namespace khinchin_detail{ + +template <class T> +T zeta_polynomial_series(T s, T sc, int digits) +{ + BOOST_MATH_STD_USING + // + // This is algorithm 3 from: + // + // "An Efficient Algorithm for the Riemann Zeta Function", P. Borwein, + // Canadian Mathematical Society, Conference Proceedings, 2000. + // See: http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P155.pdf + // + BOOST_MATH_STD_USING + int n = (digits * 19) / 53; + T sum = 0; + T two_n = ldexp(T(1), n); + int ej_sign = 1; + for(int j = 0; j < n; ++j) + { + sum += ej_sign * -two_n / pow(T(j + 1), s); + ej_sign = -ej_sign; + } + T ej_sum = 1; + T ej_term = 1; + for(int j = n; j <= 2 * n - 1; ++j) + { + sum += ej_sign * (ej_sum - two_n) / pow(T(j + 1), s); + ej_sign = -ej_sign; + ej_term *= 2 * n - j; + ej_term /= j - n + 1; + ej_sum += ej_term; + } + return -sum / (two_n * (1 - pow(T(2), sc))); +} + +template <class T> +T khinchin(int digits) +{ + BOOST_MATH_STD_USING + T sum = 0; + T term; + T lim = ldexp(T(1), 1-digits); + T factor = 0; + unsigned last_k = 1; + T num = 1; + for(unsigned n = 1;; ++n) + { + for(unsigned k = last_k; k <= 2 * n - 1; ++k) + { + factor += num / k; + num = -num; + } + last_k = 2 * n; + term = (zeta_polynomial_series(T(2 * n), T(1 - T(2 * n)), digits) - 1) * factor / n; + sum += term; + if(term < lim) + break; + } + return exp(sum / boost::math::constants::ln_two<T, boost::math::policies::policy<> >()); +} + +} + +template <class T> +template<int N> +inline T constant_khinchin<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>(); + return khinchin_detail::khinchin<T>(n); +} + +template <class T> +template<int N> +inline T constant_extreme_value_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ // from e_float constants.cpp + // Mathematica: N[12 Sqrt[6] Zeta[3]/Pi^3, 1101] + BOOST_MATH_STD_USING + T ev(12 * sqrt(static_cast<T>(6)) * zeta_three<T, policies::policy<policies::digits2<N> > >() + / pi_cubed<T, policies::policy<policies::digits2<N> > >() ); + +//T ev( +//"1.1395470994046486574927930193898461120875997958365518247216557100852480077060706857071875468869385150" +//"1894272048688553376986765366075828644841024041679714157616857834895702411080704529137366329462558680" +//"2015498788776135705587959418756809080074611906006528647805347822929577145038743873949415294942796280" +//"0895597703063466053535550338267721294164578901640163603544404938283861127819804918174973533694090594" +//"3094963822672055237678432023017824416203652657301470473548274848068762500300316769691474974950757965" +//"8640779777748741897542093874605477776538884083378029488863880220988107155275203245233994097178778984" +//"3488995668362387892097897322246698071290011857605809901090220903955815127463328974447572119951192970" +//"3684453635456559086126406960279692862247058250100678008419431185138019869693206366891639436908462809" +//"9756051372711251054914491837034685476095423926553367264355374652153595857163724698198860485357368964" +//"3807049634423621246870868566707915720704996296083373077647528285782964567312903914752617978405994377" +//"9064157147206717895272199736902453130842229559980076472936976287378945035706933650987259357729800315"); + + return ev; +} + +namespace detail{ +// +// Calculation of the Glaisher constant depends upon calculating the +// derivative of the zeta function at 2, we can then use the relation: +// zeta'(2) = 1/6 pi^2 [euler + ln(2pi)-12ln(A)] +// To get the constant A. +// See equation 45 at http://mathworld.wolfram.com/RiemannZetaFunction.html. +// +// The derivative of the zeta function is computed by direct differentiation +// of the relation: +// (1-2^(1-s))zeta(s) = SUM(n=0, INF){ (-n)^n / (n+1)^s } +// Which gives us 2 slowly converging but alternating sums to compute, +// for this we use Algorithm 1 from "Convergent Acceleration of Alternating Series", +// Henri Cohen, Fernando Rodriguez Villegas and Don Zagier, Experimental Mathematics 9:1 (1999). +// See http://www.math.utexas.edu/users/villegas/publications/conv-accel.pdf +// +template <class T> +T zeta_series_derivative_2(unsigned digits) +{ + // Derivative of the series part, evaluated at 2: + BOOST_MATH_STD_USING + int n = digits * 301 * 13 / 10000; + boost::math::itrunc((std::numeric_limits<T>::digits10 + 1) * 1.3); + T d = pow(3 + sqrt(T(8)), n); + d = (d + 1 / d) / 2; + T b = -1; + T c = -d; + T s = 0; + for(int k = 0; k < n; ++k) + { + T a = -log(T(k+1)) / ((k+1) * (k+1)); + c = b - c; + s = s + c * a; + b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1)); + } + return s / d; +} + +template <class T> +T zeta_series_2(unsigned digits) +{ + // Series part of zeta at 2: + BOOST_MATH_STD_USING + int n = digits * 301 * 13 / 10000; + T d = pow(3 + sqrt(T(8)), n); + d = (d + 1 / d) / 2; + T b = -1; + T c = -d; + T s = 0; + for(int k = 0; k < n; ++k) + { + T a = T(1) / ((k + 1) * (k + 1)); + c = b - c; + s = s + c * a; + b = (k + n) * (k - n) * b / ((k + T(0.5f)) * (k + 1)); + } + return s / d; +} + +template <class T> +inline T zeta_series_lead_2() +{ + // lead part at 2: + return 2; +} + +template <class T> +inline T zeta_series_derivative_lead_2() +{ + // derivative of lead part at 2: + return -2 * boost::math::constants::ln_two<T>(); +} + +template <class T> +inline T zeta_derivative_2(unsigned n) +{ + // zeta derivative at 2: + return zeta_series_derivative_2<T>(n) * zeta_series_lead_2<T>() + + zeta_series_derivative_lead_2<T>() * zeta_series_2<T>(n); +} + +} // namespace detail + +template <class T> +template<int N> +inline T constant_glaisher<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ + + BOOST_MATH_STD_USING + typedef policies::policy<policies::digits2<N> > forwarding_policy; + int n = N ? (std::min)(N, tools::digits<T>()) : tools::digits<T>(); + T v = detail::zeta_derivative_2<T>(n); + v *= 6; + v /= boost::math::constants::pi<T, forwarding_policy>() * boost::math::constants::pi<T, forwarding_policy>(); + v -= boost::math::constants::euler<T, forwarding_policy>(); + v -= log(2 * boost::math::constants::pi<T, forwarding_policy>()); + v /= -12; + return exp(v); + + /* + // from http://mpmath.googlecode.com/svn/data/glaisher.txt + // 20,000 digits of the Glaisher-Kinkelin constant A = exp(1/2 - zeta'(-1)) + // Computed using A = exp((6 (-zeta'(2))/pi^2 + log 2 pi + gamma)/12) + // with Euler-Maclaurin summation for zeta'(2). + T g( + "1.282427129100622636875342568869791727767688927325001192063740021740406308858826" + "46112973649195820237439420646120399000748933157791362775280404159072573861727522" + "14334327143439787335067915257366856907876561146686449997784962754518174312394652" + "76128213808180219264516851546143919901083573730703504903888123418813674978133050" + "93770833682222494115874837348064399978830070125567001286994157705432053927585405" + "81731588155481762970384743250467775147374600031616023046613296342991558095879293" + "36343887288701988953460725233184702489001091776941712153569193674967261270398013" + "52652668868978218897401729375840750167472114895288815996668743164513890306962645" + "59870469543740253099606800842447417554061490189444139386196089129682173528798629" + "88434220366989900606980888785849587494085307347117090132667567503310523405221054" + "14176776156308191919997185237047761312315374135304725819814797451761027540834943" + "14384965234139453373065832325673954957601692256427736926358821692159870775858274" + "69575162841550648585890834128227556209547002918593263079373376942077522290940187"); + + return g; + */ +} + +template <class T> +template<int N> +inline T constant_rayleigh_skewness<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ // From e_float + // 1100 digits of the Rayleigh distribution skewness + // Mathematica: N[2 Sqrt[Pi] (Pi - 3)/((4 - Pi)^(3/2)), 1100] + + BOOST_MATH_STD_USING + T rs(2 * root_pi<T, policies::policy<policies::digits2<N> > >() + * pi_minus_three<T, policies::policy<policies::digits2<N> > >() + / pow(four_minus_pi<T, policies::policy<policies::digits2<N> > >(), static_cast<T>(3./2)) + ); + // 6.31110657818937138191899351544227779844042203134719497658094585692926819617473725459905027032537306794400047264, + + //"0.6311106578189371381918993515442277798440422031347194976580945856929268196174737254599050270325373067" + //"9440004726436754739597525250317640394102954301685809920213808351450851396781817932734836994829371322" + //"5797376021347531983451654130317032832308462278373358624120822253764532674177325950686466133508511968" + //"2389168716630349407238090652663422922072397393006683401992961569208109477307776249225072042971818671" + //"4058887072693437217879039875871765635655476241624825389439481561152126886932506682176611183750503553" + //"1218982627032068396407180216351425758181396562859085306247387212297187006230007438534686340210168288" + //"8956816965453815849613622117088096547521391672977226658826566757207615552041767516828171274858145957" + //"6137539156656005855905288420585194082284972984285863898582313048515484073396332610565441264220790791" + //"0194897267890422924599776483890102027823328602965235306539844007677157873140562950510028206251529523" + //"7428049693650605954398446899724157486062545281504433364675815915402937209673727753199567661561209251" + //"4695589950526053470201635372590001578503476490223746511106018091907936826431407434894024396366284848"); ; + return rs; +} + +template <class T> +template<int N> +inline T constant_rayleigh_kurtosis_excess<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ // - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2) + // Might provide provide and calculate this using pi_minus_four. + BOOST_MATH_STD_USING + return - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >()) + - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) ) + / + ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)) + * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))) + ); +} + +template <class T> +template<int N> +inline T constant_rayleigh_kurtosis<T>::compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>)) +{ // 3 - (6 Pi^2 - 24 Pi + 16)/((Pi - 4)^2) + // Might provide provide and calculate this using pi_minus_four. + BOOST_MATH_STD_USING + return static_cast<T>(3) - (((static_cast<T>(6) * pi<T, policies::policy<policies::digits2<N> > >() + * pi<T, policies::policy<policies::digits2<N> > >()) + - (static_cast<T>(24) * pi<T, policies::policy<policies::digits2<N> > >()) + static_cast<T>(16) ) + / + ((pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4)) + * (pi<T, policies::policy<policies::digits2<N> > >() - static_cast<T>(4))) + ); +} + +}}}} // namespaces + +#endif // BOOST_MATH_CALCULATE_CONSTANTS_CONSTANTS_INCLUDED diff --git a/boost/math/constants/constants.hpp b/boost/math/constants/constants.hpp index 78aa42eff9..bb2260e7a7 100644 --- a/boost/math/constants/constants.hpp +++ b/boost/math/constants/constants.hpp @@ -1,5 +1,5 @@ -// Copyright John Maddock 2005-2006. -// Copyright Paul A. Bristow 2006-2010. +// Copyright John Maddock 2005-2006, 2011. +// Copyright Paul A. Bristow 2006-2011. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -8,6 +8,8 @@ #define BOOST_MATH_CONSTANTS_CONSTANTS_INCLUDED #include <boost/math/tools/config.hpp> +#include <boost/math/policies/policy.hpp> +#include <boost/math/tools/precision.hpp> #ifdef BOOST_MSVC #pragma warning(push) #pragma warning(disable: 4127 4701) @@ -16,6 +18,10 @@ #ifdef BOOST_MSVC #pragma warning(pop) #endif +#include <boost/mpl/if.hpp> +#include <boost/mpl/and.hpp> +#include <boost/mpl/int.hpp> +#include <boost/type_traits/is_convertible.hpp> namespace boost{ namespace math { @@ -35,43 +41,252 @@ namespace boost{ namespace math // (This is necessary because you can't use a numeric constant // since even a long double might not have enough digits). + enum construction_method + { + construct_from_float = 1, + construct_from_double = 2, + construct_from_long_double = 3, + construct_from_string = 4 + }; - #define BOOST_DEFINE_MATH_CONSTANT(name, x, y, exp)\ - template <class T> inline T name(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T))\ + // + // Max number of binary digits in the string representations + // of our constants: + // + BOOST_STATIC_CONSTANT(int, max_string_digits = (101 * 1000L) / 301L); + + template <class Real, class Policy> + struct construction_traits + { + private: + typedef typename policies::precision<Real, Policy>::type t1; + typedef typename policies::precision<float, Policy>::type t2; + typedef typename policies::precision<double, Policy>::type t3; + typedef typename policies::precision<long double, Policy>::type t4; + public: + typedef typename mpl::if_< + mpl::and_<boost::is_convertible<float, Real>, mpl::bool_< t1::value <= t2::value>, mpl::bool_<0 != t1::value> >, + mpl::int_<construct_from_float>, + typename mpl::if_< + mpl::and_<boost::is_convertible<double, Real>, mpl::bool_< t1::value <= t3::value>, mpl::bool_<0 != t1::value> >, + mpl::int_<construct_from_double>, + typename mpl::if_< + mpl::and_<boost::is_convertible<long double, Real>, mpl::bool_< t1::value <= t4::value>, mpl::bool_<0 != t1::value> >, + mpl::int_<construct_from_long_double>, + typename mpl::if_< + mpl::and_<mpl::bool_< t1::value <= max_string_digits>, mpl::bool_<0 != t1::value> >, + mpl::int_<construct_from_string>, + mpl::int_<t1::value> + >::type + >::type + >::type + >::type type; + }; + +#ifdef BOOST_HAS_THREADS +#define BOOST_MATH_CONSTANT_THREAD_HELPER(name, prefix) \ + boost::once_flag f = BOOST_ONCE_INIT;\ + boost::call_once(f, &BOOST_JOIN(BOOST_JOIN(string_, get_), name)<T>); +#else +#define BOOST_MATH_CONSTANT_THREAD_HELPER(name, prefix) +#endif + + namespace detail{ + + template <class Real> + Real convert_from_string(const char* p, const mpl::false_&) + { + return boost::lexical_cast<Real>(p); + } + template <class Real> + const char* convert_from_string(const char* p, const mpl::true_&) + { + return p; + } + + template <class T, T (*F)(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(T))> + struct constant_initializer + { + static void force_instantiate() + { + init.force_instantiate(); + } + private: + struct initializer + { + initializer() + { + F( + #ifdef BOOST_NO_EXPLICIT_FUNCTION_TEMPLATE_ARGUMENTS + 0 + #endif + ); + } + void force_instantiate()const{} + }; + static const initializer init; + }; + + template <class T, T (*F)(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(T))> + typename constant_initializer<T, F>::initializer const constant_initializer<T, F>::init; + + template <class T, int N, T (*F)(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>) BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))> + struct constant_initializer2 + { + static void force_instantiate() + { + init.force_instantiate(); + } + private: + struct initializer + { + initializer() + { + F( + #ifdef BOOST_NO_EXPLICIT_FUNCTION_TEMPLATE_ARGUMENTS + mpl::int_<N>() , 0 + #endif + ); + } + void force_instantiate()const{} + }; + static const initializer init; + }; + + template <class T, int N, T (*F)(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>) BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T))> + typename constant_initializer2<T, N, F>::initializer const constant_initializer2<T, N, F>::init; + + } + + #define BOOST_DEFINE_MATH_CONSTANT(name, x, y)\ + namespace detail{\ + template <class T> struct BOOST_JOIN(constant_, name){\ + private:\ + /* The default implementations come next: */ \ + static inline T get_from_string()\ {\ - static const T result = ::boost::lexical_cast<T>(BOOST_STRINGIZE(BOOST_JOIN(BOOST_JOIN(x, y), BOOST_JOIN(e, exp))));\ + static const T result = convert_from_string<T>(y, boost::is_convertible<const char*, T>());\ return result;\ }\ - template <> inline float name<float>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(float))\ - { return BOOST_JOIN(BOOST_JOIN(x, BOOST_JOIN(e, exp)), F); }\ - template <> inline double name<double>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(double))\ - { return BOOST_JOIN(x, BOOST_JOIN(e, exp)); }\ - template <> inline long double name<long double>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(long double))\ - { return BOOST_JOIN(BOOST_JOIN(x, BOOST_JOIN(e, exp)), L); } - - BOOST_DEFINE_MATH_CONSTANT(pi, 3.141592653589793238462643383279502884197169399375105820974944, 59230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196, 0) - BOOST_DEFINE_MATH_CONSTANT(two_pi, 6.2831853071795864769252867665590057683943388015061, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(root_pi, 1.7724538509055160272981674833411451827975, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(root_half_pi, 1.253314137315500251207882642405522626503, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(root_two_pi, 2.506628274631000502415765284811045253007, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(root_ln_four, 1.1774100225154746910115693264596996377473856893858205385225257565000, 2658854698492680841813836877081, 0) - BOOST_DEFINE_MATH_CONSTANT(e, 2.7182818284590452353602874713526624977572470936999595749669676, 27724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011, 0) - BOOST_DEFINE_MATH_CONSTANT(half, 0.5, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(euler, 0.577215664901532860606512090082402431042159335939923598805, 76723488486, 0) - BOOST_DEFINE_MATH_CONSTANT(root_two, 1.414213562373095048801688724209698078569671875376948073, 17667973799073247846210703885038753432764157273501384623091229702492483605585073721264412149709993583141322266592750559275579995050115278206, 0) - BOOST_DEFINE_MATH_CONSTANT(half_root_two, 0.70710678118654752440084436210484903928483593756084, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(ln_two, 0.693147180559945309417232121458176568075500134360255254, 120680009493393621969694715605863326996418687, 0) - BOOST_DEFINE_MATH_CONSTANT(ln_ln_two, -0.36651292058166432701243915823266946945426344783710526305367771367056, 16153193527385494558228566989083583025230453648347655663425171940646634, 0) - BOOST_DEFINE_MATH_CONSTANT(third, 0.3333333333333333333333333333333333333333333333333333333333333333333333, 3333333333333333333333333333333333333333333333333333333333333333333333333, 0) - BOOST_DEFINE_MATH_CONSTANT(twothirds, 0.66666666666666666666666666666666666666666666666666666666666666666666, 66666666666666666666666666666666666666666666666666666666666666666666667, 0) //deprecated, use two_thirds. - BOOST_DEFINE_MATH_CONSTANT(two_thirds, 0.66666666666666666666666666666666666666666666666666666666666666666666, 66666666666666666666666666666666666666666666666666666666666666666666667, 0) - BOOST_DEFINE_MATH_CONSTANT(pi_minus_three, 0.141592653589793238462643383279502884197169399375105820974944, 59230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196, 0) - BOOST_DEFINE_MATH_CONSTANT(four_minus_pi, 0.85840734641020676153735661672049711580283060062489417902505540769218359, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(pow23_four_minus_pi, 0.79531676737159754434839533505680658072763917332771320544530223438582161, 0, 0) - BOOST_DEFINE_MATH_CONSTANT(exp_minus_half, 0.6065306597126334236037995349911804534419181354871869556828921587350565194137, 484239986476115079894560, 0) - } // namespace constants + /* This one is for very high precision that is none the less known at compile time: */ \ + template <int N> static T compute(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpl::int_<N>));\ + /* public getters come next */\ + public:\ + static inline T get(const mpl::int_<construct_from_string>&)\ + {\ + constant_initializer<T, & BOOST_JOIN(constant_, name)<T>::get_from_string >::force_instantiate();\ + return get_from_string();\ + }\ + static inline BOOST_CONSTEXPR T get(const mpl::int_<construct_from_float>)\ + { return BOOST_JOIN(x, F); }\ + static inline BOOST_CONSTEXPR T get(const mpl::int_<construct_from_double>&)\ + { return x; }\ + static inline BOOST_CONSTEXPR T get(const mpl::int_<construct_from_long_double>&)\ + { return BOOST_JOIN(x, L); }\ + template <int N> static inline T get(const mpl::int_<N>& n)\ + {\ + constant_initializer2<T, N, & BOOST_JOIN(constant_, name)<T>::template compute<N> >::force_instantiate();\ + return compute<N>(); \ + }\ + /* This one is for true arbitary precision, which may well vary at runtime: */ \ + static inline T get(const mpl::int_<0>&)\ + { return tools::digits<T>() > max_string_digits ? compute<0>() : get(mpl::int_<construct_from_string>()); }\ + }; /* end of struct */\ + } /* namespace detail */ \ + \ + \ + /* The actual forwarding function: */ \ + template <class T, class Policy> inline BOOST_CONSTEXPR T name(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T) BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(Policy))\ + { return detail:: BOOST_JOIN(constant_, name)<T>::get(typename construction_traits<T, Policy>::type()); }\ + template <class T> inline BOOST_CONSTEXPR T name(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T))\ + { return name<T, boost::math::policies::policy<> >(); }\ + \ + \ + /* Now the namespace specific versions: */ \ + } namespace float_constants{ static const float name = BOOST_JOIN(x, F); }\ + namespace double_constants{ static const double name = x; } \ + namespace long_double_constants{ static const long double name = BOOST_JOIN(x, L); }\ + namespace constants{ + + BOOST_DEFINE_MATH_CONSTANT(half, 5.000000000000000000000000000000000000e-01, "5.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01") + BOOST_DEFINE_MATH_CONSTANT(third, 3.333333333333333333333333333333333333e-01, "3.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333e-01") + BOOST_DEFINE_MATH_CONSTANT(twothirds, 6.666666666666666666666666666666666666e-01, "6.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01") + BOOST_DEFINE_MATH_CONSTANT(two_thirds, 6.666666666666666666666666666666666666e-01, "6.66666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667e-01") + BOOST_DEFINE_MATH_CONSTANT(three_quarters, 7.500000000000000000000000000000000000e-01, "7.50000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000e-01") + BOOST_DEFINE_MATH_CONSTANT(root_two, 1.414213562373095048801688724209698078e+00, "1.41421356237309504880168872420969807856967187537694807317667973799073247846210703885038753432764157273501384623e+00") + BOOST_DEFINE_MATH_CONSTANT(root_three, 1.732050807568877293527446341505872366e+00, "1.73205080756887729352744634150587236694280525381038062805580697945193301690880003708114618675724857567562614142e+00") + BOOST_DEFINE_MATH_CONSTANT(half_root_two, 7.071067811865475244008443621048490392e-01, "7.07106781186547524400844362104849039284835937688474036588339868995366239231053519425193767163820786367506923115e-01") + BOOST_DEFINE_MATH_CONSTANT(ln_two, 6.931471805599453094172321214581765680e-01, "6.93147180559945309417232121458176568075500134360255254120680009493393621969694715605863326996418687542001481021e-01") + BOOST_DEFINE_MATH_CONSTANT(ln_ln_two, -3.665129205816643270124391582326694694e-01, "-3.66512920581664327012439158232669469454263447837105263053677713670561615319352738549455822856698908358302523045e-01") + BOOST_DEFINE_MATH_CONSTANT(root_ln_four, 1.177410022515474691011569326459699637e+00, "1.17741002251547469101156932645969963774738568938582053852252575650002658854698492680841813836877081106747157858e+00") + BOOST_DEFINE_MATH_CONSTANT(one_div_root_two, 7.071067811865475244008443621048490392e-01, "7.07106781186547524400844362104849039284835937688474036588339868995366239231053519425193767163820786367506923115e-01") + BOOST_DEFINE_MATH_CONSTANT(pi, 3.141592653589793238462643383279502884e+00, "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651e+00") + BOOST_DEFINE_MATH_CONSTANT(half_pi, 1.570796326794896619231321691639751442e+00, "1.57079632679489661923132169163975144209858469968755291048747229615390820314310449931401741267105853399107404326e+00") + BOOST_DEFINE_MATH_CONSTANT(third_pi, 1.047197551196597746154214461093167628e+00, "1.04719755119659774615421446109316762806572313312503527365831486410260546876206966620934494178070568932738269550e+00") + BOOST_DEFINE_MATH_CONSTANT(sixth_pi, 5.235987755982988730771072305465838140e-01, "5.23598775598298873077107230546583814032861566562517636829157432051302734381034833104672470890352844663691347752e-01") + BOOST_DEFINE_MATH_CONSTANT(two_pi, 6.283185307179586476925286766559005768e+00, "6.28318530717958647692528676655900576839433879875021164194988918461563281257241799725606965068423413596429617303e+00") + BOOST_DEFINE_MATH_CONSTANT(two_thirds_pi, 2.094395102393195492308428922186335256e+00, "2.09439510239319549230842892218633525613144626625007054731662972820521093752413933241868988356141137865476539101e+00") + BOOST_DEFINE_MATH_CONSTANT(three_quarters_pi, 2.356194490192344928846982537459627163e+00, "2.35619449019234492884698253745962716314787704953132936573120844423086230471465674897102611900658780098661106488e+00") + BOOST_DEFINE_MATH_CONSTANT(four_thirds_pi, 4.188790204786390984616857844372670512e+00, "4.18879020478639098461685784437267051226289253250014109463325945641042187504827866483737976712282275730953078202e+00") + BOOST_DEFINE_MATH_CONSTANT(one_div_two_pi, 1.591549430918953357688837633725143620e-01, "1.59154943091895335768883763372514362034459645740456448747667344058896797634226535090113802766253085956072842727e-01") + BOOST_DEFINE_MATH_CONSTANT(one_div_root_two_pi, 3.989422804014326779399460599343818684e-01, "3.98942280401432677939946059934381868475858631164934657665925829670657925899301838501252333907306936430302558863e-01") + BOOST_DEFINE_MATH_CONSTANT(root_pi, 1.772453850905516027298167483341145182e+00, "1.77245385090551602729816748334114518279754945612238712821380778985291128459103218137495065673854466541622682362e+00") + BOOST_DEFINE_MATH_CONSTANT(root_half_pi, 1.253314137315500251207882642405522626e+00, "1.25331413731550025120788264240552262650349337030496915831496178817114682730392098747329791918902863305800498633e+00") + BOOST_DEFINE_MATH_CONSTANT(root_two_pi, 2.506628274631000502415765284811045253e+00, "2.50662827463100050241576528481104525300698674060993831662992357634229365460784197494659583837805726611600997267e+00") + BOOST_DEFINE_MATH_CONSTANT(one_div_root_pi, 5.641895835477562869480794515607725858e-01, "5.64189583547756286948079451560772585844050629328998856844085721710642468441493414486743660202107363443028347906e-01") + BOOST_DEFINE_MATH_CONSTANT(root_one_div_pi, 5.641895835477562869480794515607725858e-01, "5.64189583547756286948079451560772585844050629328998856844085721710642468441493414486743660202107363443028347906e-01") + BOOST_DEFINE_MATH_CONSTANT(pi_minus_three, 1.415926535897932384626433832795028841e-01, "1.41592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513e-01") + BOOST_DEFINE_MATH_CONSTANT(four_minus_pi, 8.584073464102067615373566167204971158e-01, "8.58407346410206761537356616720497115802830600624894179025055407692183593713791001371965174657882932017851913487e-01") + BOOST_DEFINE_MATH_CONSTANT(pow23_four_minus_pi, 7.953167673715975443483953350568065807e-01, "7.95316767371597544348395335056806580727639173327713205445302234388856268267518187590758006888600828436839800178e-01") + BOOST_DEFINE_MATH_CONSTANT(pi_pow_e, 2.245915771836104547342715220454373502e+01, "2.24591577183610454734271522045437350275893151339966922492030025540669260403991179123185197527271430315314500731e+01") + BOOST_DEFINE_MATH_CONSTANT(pi_sqr, 9.869604401089358618834490999876151135e+00, "9.86960440108935861883449099987615113531369940724079062641334937622004482241920524300177340371855223182402591377e+00") + BOOST_DEFINE_MATH_CONSTANT(pi_sqr_div_six, 1.644934066848226436472415166646025189e+00, "1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890061975870530400431896e+00") + BOOST_DEFINE_MATH_CONSTANT(pi_cubed, 3.100627668029982017547631506710139520e+01, "3.10062766802998201754763150671013952022252885658851076941445381038063949174657060375667010326028861930301219616e+01") + BOOST_DEFINE_MATH_CONSTANT(cbrt_pi, 1.464591887561523263020142527263790391e+00, "1.46459188756152326302014252726379039173859685562793717435725593713839364979828626614568206782035382089750397002e+00") + BOOST_DEFINE_MATH_CONSTANT(one_div_cbrt_pi, 6.827840632552956814670208331581645981e-01, "6.82784063255295681467020833158164598108367515632448804042681583118899226433403918237673501922595519865685577274e-01") + BOOST_DEFINE_MATH_CONSTANT(e, 2.718281828459045235360287471352662497e+00, "2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138217852516642742746639193e+00") + BOOST_DEFINE_MATH_CONSTANT(exp_minus_half, 6.065306597126334236037995349911804534e-01, "6.06530659712633423603799534991180453441918135487186955682892158735056519413748423998647611507989456026423789794e-01") + BOOST_DEFINE_MATH_CONSTANT(e_pow_pi, 2.314069263277926900572908636794854738e+01, "2.31406926327792690057290863679485473802661062426002119934450464095243423506904527835169719970675492196759527048e+01") + BOOST_DEFINE_MATH_CONSTANT(root_e, 1.648721270700128146848650787814163571e+00, "1.64872127070012814684865078781416357165377610071014801157507931164066102119421560863277652005636664300286663776e+00") + BOOST_DEFINE_MATH_CONSTANT(log10_e, 4.342944819032518276511289189166050822e-01, "4.34294481903251827651128918916605082294397005803666566114453783165864649208870774729224949338431748318706106745e-01") + BOOST_DEFINE_MATH_CONSTANT(one_div_log10_e, 2.302585092994045684017991454684364207e+00, "2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404e+00") + BOOST_DEFINE_MATH_CONSTANT(ln_ten, 2.302585092994045684017991454684364207e+00, "2.30258509299404568401799145468436420760110148862877297603332790096757260967735248023599720508959829834196778404e+00") + BOOST_DEFINE_MATH_CONSTANT(degree, 1.745329251994329576923690768488612713e-02, "1.74532925199432957692369076848861271344287188854172545609719144017100911460344944368224156963450948221230449251e-02") + BOOST_DEFINE_MATH_CONSTANT(radian, 5.729577951308232087679815481410517033e+01, "5.72957795130823208767981548141051703324054724665643215491602438612028471483215526324409689958511109441862233816e+01") + BOOST_DEFINE_MATH_CONSTANT(sin_one, 8.414709848078965066525023216302989996e-01, "8.41470984807896506652502321630298999622563060798371065672751709991910404391239668948639743543052695854349037908e-01") + BOOST_DEFINE_MATH_CONSTANT(cos_one, 5.403023058681397174009366074429766037e-01, "5.40302305868139717400936607442976603732310420617922227670097255381100394774471764517951856087183089343571731160e-01") + BOOST_DEFINE_MATH_CONSTANT(sinh_one, 1.175201193643801456882381850595600815e+00, "1.17520119364380145688238185059560081515571798133409587022956541301330756730432389560711745208962339184041953333e+00") + BOOST_DEFINE_MATH_CONSTANT(cosh_one, 1.543080634815243778477905620757061682e+00, "1.54308063481524377847790562075706168260152911236586370473740221471076906304922369896426472643554303558704685860e+00") + BOOST_DEFINE_MATH_CONSTANT(phi, 1.618033988749894848204586834365638117e+00, "1.61803398874989484820458683436563811772030917980576286213544862270526046281890244970720720418939113748475408808e+00") + BOOST_DEFINE_MATH_CONSTANT(ln_phi, 4.812118250596034474977589134243684231e-01, "4.81211825059603447497758913424368423135184334385660519661018168840163867608221774412009429122723474997231839958e-01") + BOOST_DEFINE_MATH_CONSTANT(one_div_ln_phi, 2.078086921235027537601322606117795767e+00, "2.07808692123502753760132260611779576774219226778328348027813992191974386928553540901445615414453604821933918634e+00") + BOOST_DEFINE_MATH_CONSTANT(euler, 5.772156649015328606065120900824024310e-01, "5.77215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447250e-01") + BOOST_DEFINE_MATH_CONSTANT(one_div_euler, 1.732454714600633473583025315860829681e+00, "1.73245471460063347358302531586082968115577655226680502204843613287065531408655243008832840219409928068072365714e+00") + BOOST_DEFINE_MATH_CONSTANT(euler_sqr, 3.331779238077186743183761363552442266e-01, "3.33177923807718674318376136355244226659417140249629743150833338002265793695756669661263268631715977303039565603e-01") + BOOST_DEFINE_MATH_CONSTANT(zeta_two, 1.644934066848226436472415166646025189e+00, "1.64493406684822643647241516664602518921894990120679843773555822937000747040320087383362890061975870530400431896e+00") + BOOST_DEFINE_MATH_CONSTANT(zeta_three, 1.202056903159594285399738161511449990e+00, "1.20205690315959428539973816151144999076498629234049888179227155534183820578631309018645587360933525814619915780e+00") + BOOST_DEFINE_MATH_CONSTANT(catalan, 9.159655941772190150546035149323841107e-01, "9.15965594177219015054603514932384110774149374281672134266498119621763019776254769479356512926115106248574422619e-01") + BOOST_DEFINE_MATH_CONSTANT(glaisher, 1.282427129100622636875342568869791727e+00, "1.28242712910062263687534256886979172776768892732500119206374002174040630885882646112973649195820237439420646120e+00") + BOOST_DEFINE_MATH_CONSTANT(khinchin, 2.685452001065306445309714835481795693e+00, "2.68545200106530644530971483548179569382038229399446295305115234555721885953715200280114117493184769799515346591e+00") + BOOST_DEFINE_MATH_CONSTANT(extreme_value_skewness, 1.139547099404648657492793019389846112e+00, "1.13954709940464865749279301938984611208759979583655182472165571008524800770607068570718754688693851501894272049e+00") + BOOST_DEFINE_MATH_CONSTANT(rayleigh_skewness, 6.311106578189371381918993515442277798e-01, "6.31110657818937138191899351544227779844042203134719497658094585692926819617473725459905027032537306794400047264e-01") + BOOST_DEFINE_MATH_CONSTANT(rayleigh_kurtosis, 3.245089300687638062848660410619754415e+00, "3.24508930068763806284866041061975441541706673178920936177133764493367904540874159051490619368679348977426462633e+00") + BOOST_DEFINE_MATH_CONSTANT(rayleigh_kurtosis_excess, 2.450893006876380628486604106197544154e-01, "2.45089300687638062848660410619754415417066731789209361771337644933679045408741590514906193686793489774264626328e-01") + + BOOST_DEFINE_MATH_CONSTANT(two_div_pi, 6.366197723675813430755350534900574481e-01, "6.36619772367581343075535053490057448137838582961825794990669376235587190536906140360455211065012343824291370907e-01") + BOOST_DEFINE_MATH_CONSTANT(root_two_div_pi, 7.978845608028653558798921198687637369e-01, "7.97884560802865355879892119868763736951717262329869315331851659341315851798603677002504667814613872860605117725e-01") + + +} // namespace constants } // namespace math } // namespace boost +// +// We deliberately include this *after* all the declarations above, +// that way the calculation routines can call on other constants above: +// +#include <boost/math/constants/calculate_constants.hpp> + #endif // BOOST_MATH_CONSTANTS_CONSTANTS_INCLUDED + diff --git a/boost/math/constants/generate.hpp b/boost/math/constants/generate.hpp new file mode 100644 index 0000000000..dfb15633a5 --- /dev/null +++ b/boost/math/constants/generate.hpp @@ -0,0 +1,76 @@ +// Copyright John Maddock 2010. +// 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_CONSTANTS_GENERATE_INCLUDED +#define BOOST_MATH_CONSTANTS_GENERATE_INCLUDED + +#include <boost/math/constants/constants.hpp> +#include <boost/regex.hpp> +#include <iostream> +#include <iomanip> +#include <sstream> + +#ifdef USE_MPFR +#include <boost/math/bindings/mpfr.hpp> +#elif defined(USE_MPREAL) +#include <boost/math/bindings/mpreal.hpp> +#elif defined(USE_CPP_FLOAT) +#include <boost/multiprecision/cpp_float.hpp> +#else +#include <boost/math/bindings/rr.hpp> +#endif + +namespace boost{ namespace math{ namespace constants{ + +#ifdef USE_MPFR +typedef mpfr_class generator_type; +#elif defined(USE_MPREAL) +typedef mpfr::mpreal generator_type; +#elif defined(USE_CPP_FLOAT) +typedef boost::multiprecision::mp_number<boost::multiprecision::cpp_float<500> > generator_type; +#else +typedef ntl::RR generator_type; +#endif + +inline void print_constant(const char* name, generator_type(*f)(const mpl::int_<0>&)) +{ +#ifdef USE_MPFR + mpfr_class::set_dprec(((200 + 1) * 1000L) / 301L); +#elif defined(USE_MPREAL) + mpfr::mpreal::set_default_prec(((200 + 1) * 1000L) / 301L); +#elif defined(USE_CPP_FLOAT) + // Nothing to do, precision is already set. +#else + ntl::RR::SetPrecision(((200 + 1) * 1000L) / 301L); + ntl::RR::SetOutputPrecision(102); +#endif + generator_type value = f(boost::mpl::int_<0>()); + std::stringstream os; + os << std::setprecision(110) << std::scientific; + os << value; + std::string s = os.str(); + static const regex e("([+-]?\\d+(?:\\.\\d{0,36})?)(\\d*)(?:e([+-]?\\d+))?"); + smatch what; + if(regex_match(s, what, e)) + { + std::cout << + "BOOST_DEFINE_MATH_CONSTANT(" << name << ", " + << what[1] << "e" << (what[3].length() ? what[3].str() : std::string("0")) << ", " + << "\"" << what[1] << what[2] << "e" << (what[3].length() ? what[3].str() : std::string("0")) + << "\");" << std::endl; + } + else + { + std::cout << "Format of numeric constant was not recognised!!" << std::endl; + } +} + +#define BOOST_CONSTANTS_GENERATE(name) \ + boost::math::constants::print_constant(#name, \ + & boost::math::constants::detail::BOOST_JOIN(constant_, name)<boost::math::constants::generator_type>::get) + +}}} // namespaces + +#endif // BOOST_MATH_CONSTANTS_GENERATE_INCLUDED diff --git a/boost/math/constants/info.hpp b/boost/math/constants/info.hpp new file mode 100644 index 0000000000..e3e6a517d4 --- /dev/null +++ b/boost/math/constants/info.hpp @@ -0,0 +1,163 @@ +// Copyright John Maddock 2010. +// 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) + +#ifdef _MSC_VER +# pragma once +#endif + +#ifndef BOOST_MATH_CONSTANTS_INFO_INCLUDED +#define BOOST_MATH_CONSTANTS_INFO_INCLUDED + +#include <boost/math/constants/constants.hpp> +#include <iostream> +#include <iomanip> +#include <typeinfo> + +namespace boost{ namespace math{ namespace constants{ + + namespace detail{ + + template <class T> + const char* nameof(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) + { + return typeid(T).name(); + } + template <> + const char* nameof<float>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(float)) + { + return "float"; + } + template <> + const char* nameof<double>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(double)) + { + return "double"; + } + template <> + const char* nameof<long double>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(long double)) + { + return "long double"; + } + + } + +template <class T, class Policy> +void print_info_on_type(std::ostream& os = std::cout BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T) BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(Policy)) +{ + using detail::nameof; +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4127) +#endif + os << + "Information on the Implementation and Handling of \n" + "Mathematical Constants for Type " << nameof<T>() << + "\n\n" + "Checking for std::numeric_limits<" << nameof<T>() << "> specialisation: " << + (std::numeric_limits<T>::is_specialized ? "yes" : "no") << std::endl; + if(std::numeric_limits<T>::is_specialized) + { + os << + "std::numeric_limits<" << nameof<T>() << ">::digits reports that the radix is " << std::numeric_limits<T>::radix << ".\n"; + if (std::numeric_limits<T>::radix == 2) + { + os << + "std::numeric_limits<" << nameof<T>() << ">::digits reports that the precision is \n" << std::numeric_limits<T>::digits << " binary digits.\n"; + } + else if (std::numeric_limits<T>::radix == 10) + { + os << + "std::numeric_limits<" << nameof<T>() << ">::digits reports that the precision is \n" << std::numeric_limits<T>::digits10 << " decimal digits.\n"; + os << + "std::numeric_limits<" << nameof<T>() << ">::digits reports that the precision is \n" + << std::numeric_limits<T>::digits * 1000L /301L << " binary digits.\n"; // divide by log2(10) - about 3 bits per decimal digit. + } + else + { + os << "Unknown radix = " << std::numeric_limits<T>::radix << "\n"; + } + } + typedef typename boost::math::policies::precision<T, Policy>::type precision_type; + if(precision_type::value) + { + if (std::numeric_limits<T>::radix == 2) + { + os << + "boost::math::policies::precision<" << nameof<T>() << ", " << nameof<Policy>() << " reports that the compile time precision is \n" << precision_type::value << " binary digits.\n"; + } + else if (std::numeric_limits<T>::radix == 10) + { + os << + "boost::math::policies::precision<" << nameof<T>() << ", " << nameof<Policy>() << " reports that the compile time precision is \n" << precision_type::value << " binary digits.\n"; + } + else + { + os << "Unknown radix = " << std::numeric_limits<T>::radix << "\n"; + } + } + else + { + os << + "boost::math::policies::precision<" << nameof<T>() << ", Policy> \n" + "reports that there is no compile type precision available.\n" + "boost::math::tools::digits<" << nameof<T>() << ">() \n" + "reports that the current runtime precision is \n" << + boost::math::tools::digits<T>() << " binary digits.\n"; + } + + typedef typename construction_traits<T, Policy>::type construction_type; + + switch(construction_type::value) + { + case 0: + os << + "No compile time precision is available, the construction method \n" + "will be decided at runtime and results will not be cached \n" + "- this may lead to poor runtime performance.\n" + "Current runtime precision indicates that\n"; + if(boost::math::tools::digits<T>() > max_string_digits) + { + os << "the constant will be recalculated on each call.\n"; + } + else + { + os << "the constant will be constructed from a string on each call.\n"; + } + break; + case 1: + os << + "The constant will be constructed from a float.\n"; + break; + case 2: + os << + "The constant will be constructed from a double.\n"; + break; + case 3: + os << + "The constant will be constructed from a long double.\n"; + break; + case 4: + os << + "The constant will be constructed from a string (and the result cached).\n"; + break; + default: + os << + "The constant will be calculated (and the result cached).\n"; + break; + } + os << std::endl; +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif +} + +template <class T> +void print_info_on_type(std::ostream& os = std::cout BOOST_MATH_APPEND_EXPLICIT_TEMPLATE_TYPE_SPEC(T)) +{ + print_info_on_type<T, boost::math::policies::policy<> >(os); +} + +}}} // namespaces + +#endif // BOOST_MATH_CONSTANTS_INFO_INCLUDED diff --git a/boost/math/distributions/detail/inv_discrete_quantile.hpp b/boost/math/distributions/detail/inv_discrete_quantile.hpp index 4d28e52936..9397e7c7c2 100644 --- a/boost/math/distributions/detail/inv_discrete_quantile.hpp +++ b/boost/math/distributions/detail/inv_discrete_quantile.hpp @@ -124,6 +124,8 @@ typename Dist::value_type --count; if(fb == 0) return b; + if(a == b) + return b; // can't go any higher! } else { @@ -135,6 +137,8 @@ typename Dist::value_type --count; if(fa == 0) return a; + if(a == b) + return a; // We can't go any lower than this! } } } @@ -208,7 +212,7 @@ typename Dist::value_type // Zero is to the right of x2, so walk upwards // until we find it: // - while((boost::math::sign)(fb) == (boost::math::sign)(fa)) + while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) { if(count == 0) policies::raise_evaluation_error(function, "Unable to bracket root, last nearest value was %1%", b, policy_type()); @@ -228,7 +232,7 @@ typename Dist::value_type // Zero is to the left of a, so walk downwards // until we find it: // - while((boost::math::sign)(fb) == (boost::math::sign)(fa)) + while(((boost::math::sign)(fb) == (boost::math::sign)(fa)) && (a != b)) { if(fabs(a) < tools::min_value<value_type>()) { @@ -255,6 +259,8 @@ typename Dist::value_type return a; if(fb == 0) return b; + if(a == b) + return b; // Ran out of bounds trying to bracket - there is no answer! // // Adjust bounds so that if we're looking for an integer // result, then both ends round the same way: diff --git a/boost/math/distributions/skew_normal.hpp b/boost/math/distributions/skew_normal.hpp new file mode 100644 index 0000000000..58526defdb --- /dev/null +++ b/boost/math/distributions/skew_normal.hpp @@ -0,0 +1,716 @@ +// (C) Benjamin Sobotta 2012 + +// 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_STATS_SKEW_NORMAL_HPP +#define BOOST_STATS_SKEW_NORMAL_HPP + +// http://en.wikipedia.org/wiki/Skew_normal_distribution +// http://azzalini.stat.unipd.it/SN/ +// Also: +// Azzalini, A. (1985). "A class of distributions which includes the normal ones". +// Scand. J. Statist. 12: 171-178. + +#include <boost/math/distributions/fwd.hpp> // TODO add skew_normal distribution to fwd.hpp! +#include <boost/math/special_functions/owens_t.hpp> // Owen's T function +#include <boost/math/distributions/complement.hpp> +#include <boost/math/distributions/normal.hpp> +#include <boost/math/distributions/detail/common_error_handling.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/tools/tuple.hpp> +#include <boost/math/tools/roots.hpp> // Newton-Raphson +#include <boost/assert.hpp> +#include <boost/math/distributions/detail/generic_mode.hpp> // pdf max finder. + +#include <utility> +#include <algorithm> // std::lower_bound, std::distance + +namespace boost{ namespace math{ + + namespace detail + { + template <class RealType, class Policy> + inline bool check_skew_normal_shape( + const char* function, + RealType shape, + RealType* result, + const Policy& pol) + { + if(!(boost::math::isfinite)(shape)) + { + *result = + policies::raise_domain_error<RealType>(function, + "Shape parameter is %1%, but must be finite!", + shape, pol); + return false; + } + return true; + } + + } // namespace detail + + template <class RealType = double, class Policy = policies::policy<> > + class skew_normal_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + skew_normal_distribution(RealType location = 0, RealType scale = 1, RealType shape = 0) + : location_(location), scale_(scale), shape_(shape) + { // Default is a 'standard' normal distribution N01. (shape=0 results in the normal distribution with no skew) + static const char* function = "boost::math::skew_normal_distribution<%1%>::skew_normal_distribution"; + + RealType result; + detail::check_scale(function, scale, &result, Policy()); + detail::check_location(function, location, &result, Policy()); + detail::check_skew_normal_shape(function, shape, &result, Policy()); + } + + RealType location()const + { + return location_; + } + + RealType scale()const + { + return scale_; + } + + RealType shape()const + { + return shape_; + } + + + private: + // + // Data members: + // + RealType location_; // distribution location. + RealType scale_; // distribution scale. + RealType shape_; // distribution shape. + }; // class skew_normal_distribution + + typedef skew_normal_distribution<double> skew_normal; + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const skew_normal_distribution<RealType, Policy>& /*dist*/) + { // Range of permissible values for random variable x. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const skew_normal_distribution<RealType, Policy>& /*dist*/) + { // Range of supported values for random variable x. + // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. + + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value. + } + + template <class RealType, class Policy> + inline RealType pdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::pdf(const skew_normal_distribution<%1%>&, %1%)"; + if((boost::math::isinf)(x)) + { + return 0; // pdf + and - infinity is zero. + } + // Below produces MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity()) + //{ // pdf + and - infinity is zero. + // return 0; + //} + + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, location, &result, Policy())) + { + return result; + } + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + { + return result; + } + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType transformed_x = (x-location)/scale; + + normal_distribution<RealType, Policy> std_normal; + + result = pdf(std_normal, transformed_x) * cdf(std_normal, shape*transformed_x) * 2 / scale; + + return result; + } // pdf + + template <class RealType, class Policy> + inline RealType cdf(const skew_normal_distribution<RealType, Policy>& dist, const RealType& x) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::cdf(const skew_normal_distribution<%1%>&, %1%)"; + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + { + return result; + } + if(false == detail::check_location(function, location, &result, Policy())) + { + return result; + } + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + { + return result; + } + if((boost::math::isinf)(x)) + { + if(x < 0) return 0; // -infinity + return 1; // + infinity + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + //{ // cdf +infinity is unity. + // return 1; + //} + //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + //{ // cdf -infinity is zero. + // return 0; + //} + if(false == detail::check_x(function, x, &result, Policy())) + { + return result; + } + + const RealType transformed_x = (x-location)/scale; + + normal_distribution<RealType, Policy> std_normal; + + result = cdf(std_normal, transformed_x) - owens_t(transformed_x, shape)*static_cast<RealType>(2); + + return result; + } // cdf + + template <class RealType, class Policy> + inline RealType cdf(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c) + { + const RealType scale = c.dist.scale(); + const RealType location = c.dist.location(); + const RealType shape = c.dist.shape(); + const RealType x = c.param; + + static const char* function = "boost::math::cdf(const complement(skew_normal_distribution<%1%>&), %1%)"; + + if((boost::math::isinf)(x)) + { + if(x < 0) return 1; // cdf complement -infinity is unity. + return 0; // cdf complement +infinity is zero + } + // These produce MSVC 4127 warnings, so the above used instead. + //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity()) + //{ // cdf complement +infinity is zero. + // return 0; + //} + //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity()) + //{ // cdf complement -infinity is unity. + // return 1; + //} + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + if(false == detail::check_x(function, x, &result, Policy())) + return result; + + const RealType transformed_x = (x-location)/scale; + + normal_distribution<RealType, Policy> std_normal; + + result = cdf(complement(std_normal, transformed_x)) + owens_t(transformed_x, shape)*static_cast<RealType>(2); + return result; + } // cdf complement + + template <class RealType, class Policy> + inline RealType location(const skew_normal_distribution<RealType, Policy>& dist) + { + return dist.location(); + } + + template <class RealType, class Policy> + inline RealType scale(const skew_normal_distribution<RealType, Policy>& dist) + { + return dist.scale(); + } + + template <class RealType, class Policy> + inline RealType shape(const skew_normal_distribution<RealType, Policy>& dist) + { + return dist.shape(); + } + + template <class RealType, class Policy> + inline RealType mean(const skew_normal_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // for ADL of std functions + + using namespace boost::math::constants; + + //const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape()); + + //return dist.location() + dist.scale() * delta * root_two_div_pi<RealType>(); + + return dist.location() + dist.scale() * dist.shape() / sqrt(pi<RealType>()+pi<RealType>()*dist.shape()*dist.shape()) * root_two<RealType>(); + } + + template <class RealType, class Policy> + inline RealType variance(const skew_normal_distribution<RealType, Policy>& dist) + { + using namespace boost::math::constants; + + const RealType delta2 = static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())); + //const RealType inv_delta2 = static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape()); + + RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()*delta2); + //RealType variance = dist.scale()*dist.scale()*(static_cast<RealType>(1)-two_div_pi<RealType>()/inv_delta2); + + return variance; + } + + namespace detail + { + /* + TODO No closed expression for mode, so use max of pdf. + */ + + template <class RealType, class Policy> + inline RealType mode_fallback(const skew_normal_distribution<RealType, Policy>& dist) + { // mode. + static const char* function = "mode(skew_normal_distribution<%1%> const&)"; + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + RealType result; + if(!detail::check_scale( + function, + scale, &result, Policy()) + || + !detail::check_skew_normal_shape( + function, + shape, + &result, + Policy())) + return result; + + if( shape == 0 ) + { + return location; + } + + if( shape < 0 ) + { + skew_normal_distribution<RealType, Policy> D(0, 1, -shape); + result = mode_fallback(D); + result = location-scale*result; + return result; + } + + BOOST_MATH_STD_USING + + // 21 elements + static const RealType shapes[] = { + 0.0, + 1.000000000000000e-004, + 2.069138081114790e-004, + 4.281332398719396e-004, + 8.858667904100824e-004, + 1.832980710832436e-003, + 3.792690190732250e-003, + 7.847599703514606e-003, + 1.623776739188722e-002, + 3.359818286283781e-002, + 6.951927961775606e-002, + 1.438449888287663e-001, + 2.976351441631319e-001, + 6.158482110660261e-001, + 1.274274985703135e+000, + 2.636650898730361e+000, + 5.455594781168514e+000, + 1.128837891684688e+001, + 2.335721469090121e+001, + 4.832930238571753e+001, + 1.000000000000000e+002}; + + // 21 elements + static const RealType guess[] = { + 0.0, + 5.000050000525391e-005, + 1.500015000148736e-004, + 3.500035000350010e-004, + 7.500075000752560e-004, + 1.450014500145258e-003, + 3.050030500305390e-003, + 6.250062500624765e-003, + 1.295012950129504e-002, + 2.675026750267495e-002, + 5.525055250552491e-002, + 1.132511325113255e-001, + 2.249522495224952e-001, + 3.992539925399257e-001, + 5.353553535535358e-001, + 4.954549545495457e-001, + 3.524535245352451e-001, + 2.182521825218249e-001, + 1.256512565125654e-001, + 6.945069450694508e-002, + 3.735037350373460e-002 + }; + + const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape); + + typedef typename std::iterator_traits<RealType*>::difference_type diff_type; + + const diff_type d = std::distance(shapes, result_ptr); + + BOOST_ASSERT(d > static_cast<diff_type>(0)); + + // refine + if(d < static_cast<diff_type>(21)) // shape smaller 100 + { + result = guess[d-static_cast<diff_type>(1)] + + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)]) + * (shape-shapes[d-static_cast<diff_type>(1)]); + } + else // shape greater 100 + { + result = 1e-4; + } + + skew_normal_distribution<RealType, Policy> helper(0, 1, shape); + + result = detail::generic_find_mode_01(helper, result, function); + + result = result*scale + location; + + return result; + } // mode_fallback + + + /* + * TODO No closed expression for mode, so use f'(x) = 0 + */ + template <class RealType, class Policy> + struct skew_normal_mode_functor + { + skew_normal_mode_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist) + : distribution(dist) + { + } + + boost::math::tuple<RealType, RealType> operator()(RealType const& x) + { + normal_distribution<RealType, Policy> std_normal; + const RealType shape = distribution.shape(); + const RealType pdf_x = pdf(distribution, x); + const RealType normpdf_x = pdf(std_normal, x); + const RealType normpdf_ax = pdf(std_normal, x*shape); + RealType fx = static_cast<RealType>(2)*shape*normpdf_ax*normpdf_x - x*pdf_x; + RealType dx = static_cast<RealType>(2)*shape*x*normpdf_x*normpdf_ax*(static_cast<RealType>(1) + shape*shape) + pdf_x + x*fx; + // return both function evaluation difference f(x) and 1st derivative f'(x). + return boost::math::make_tuple(fx, -dx); + } + private: + const boost::math::skew_normal_distribution<RealType, Policy> distribution; + }; + + } // namespace detail + + template <class RealType, class Policy> + inline RealType mode(const skew_normal_distribution<RealType, Policy>& dist) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::mode(const skew_normal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + + if( shape == 0 ) + { + return location; + } + + if( shape < 0 ) + { + skew_normal_distribution<RealType, Policy> D(0, 1, -shape); + result = mode(D); + result = location-scale*result; + return result; + } + + // 21 elements + static const RealType shapes[] = { + 0.0, + 1.000000000000000e-004, + 2.069138081114790e-004, + 4.281332398719396e-004, + 8.858667904100824e-004, + 1.832980710832436e-003, + 3.792690190732250e-003, + 7.847599703514606e-003, + 1.623776739188722e-002, + 3.359818286283781e-002, + 6.951927961775606e-002, + 1.438449888287663e-001, + 2.976351441631319e-001, + 6.158482110660261e-001, + 1.274274985703135e+000, + 2.636650898730361e+000, + 5.455594781168514e+000, + 1.128837891684688e+001, + 2.335721469090121e+001, + 4.832930238571753e+001, + 1.000000000000000e+002}; + + // 21 elements + static const RealType guess[] = { + 0.0, + 5.000050000525391e-005, + 1.500015000148736e-004, + 3.500035000350010e-004, + 7.500075000752560e-004, + 1.450014500145258e-003, + 3.050030500305390e-003, + 6.250062500624765e-003, + 1.295012950129504e-002, + 2.675026750267495e-002, + 5.525055250552491e-002, + 1.132511325113255e-001, + 2.249522495224952e-001, + 3.992539925399257e-001, + 5.353553535535358e-001, + 4.954549545495457e-001, + 3.524535245352451e-001, + 2.182521825218249e-001, + 1.256512565125654e-001, + 6.945069450694508e-002, + 3.735037350373460e-002 + }; + + const RealType* result_ptr = std::lower_bound(shapes, shapes+21, shape); + + typedef typename std::iterator_traits<RealType*>::difference_type diff_type; + + const diff_type d = std::distance(shapes, result_ptr); + + BOOST_ASSERT(d > static_cast<diff_type>(0)); + + // TODO: make the search bounds smarter, depending on the shape parameter + RealType search_min = 0; // below zero was caught above + RealType search_max = 0.55; // will never go above 0.55 + + // refine + if(d < static_cast<diff_type>(21)) // shape smaller 100 + { + // it is safe to assume that d > 0, because shape==0.0 is caught earlier + result = guess[d-static_cast<diff_type>(1)] + + (guess[d]-guess[d-static_cast<diff_type>(1)])/(shapes[d]-shapes[d-static_cast<diff_type>(1)]) + * (shape-shapes[d-static_cast<diff_type>(1)]); + } + else // shape greater 100 + { + result = 1e-4; + search_max = guess[19]; // set 19 instead of 20 to have a safety margin because the table may not be exact @ shape=100 + } + + const int get_digits = policies::digits<RealType, Policy>();// get digits from policy, + boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations. + + skew_normal_distribution<RealType, Policy> helper(0, 1, shape); + + result = tools::newton_raphson_iterate(detail::skew_normal_mode_functor<RealType, Policy>(helper), result, + search_min, search_max, get_digits, m); + + result = result*scale + location; + + return result; + } + + + + template <class RealType, class Policy> + inline RealType skewness(const skew_normal_distribution<RealType, Policy>& dist) + { + BOOST_MATH_STD_USING // for ADL of std functions + using namespace boost::math::constants; + + static const RealType factor = four_minus_pi<RealType>()/static_cast<RealType>(2); + const RealType delta = dist.shape() / sqrt(static_cast<RealType>(1)+dist.shape()*dist.shape()); + + return factor * pow(root_two_div_pi<RealType>() * delta, 3) / + pow(static_cast<RealType>(1)-two_div_pi<RealType>()*delta*delta, static_cast<RealType>(1.5)); + } + + template <class RealType, class Policy> + inline RealType kurtosis(const skew_normal_distribution<RealType, Policy>& dist) + { + return kurtosis_excess(dist)+static_cast<RealType>(3); + } + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const skew_normal_distribution<RealType, Policy>& dist) + { + using namespace boost::math::constants; + + static const RealType factor = pi_minus_three<RealType>()*static_cast<RealType>(2); + + const RealType delta2 = static_cast<RealType>(1) / (static_cast<RealType>(1)+static_cast<RealType>(1)/(dist.shape()*dist.shape())); + + const RealType x = static_cast<RealType>(1)-two_div_pi<RealType>()*delta2; + const RealType y = two_div_pi<RealType>() * delta2; + + return factor * y*y / (x*x); + } + + namespace detail + { + + template <class RealType, class Policy> + struct skew_normal_quantile_functor + { + skew_normal_quantile_functor(const boost::math::skew_normal_distribution<RealType, Policy> dist, RealType const& p) + : distribution(dist), prob(p) + { + } + + boost::math::tuple<RealType, RealType> operator()(RealType const& x) + { + RealType c = cdf(distribution, x); + RealType fx = c - prob; // Difference cdf - value - to minimize. + RealType dx = pdf(distribution, x); // pdf is 1st derivative. + // return both function evaluation difference f(x) and 1st derivative f'(x). + return boost::math::make_tuple(fx, dx); + } + private: + const boost::math::skew_normal_distribution<RealType, Policy> distribution; + RealType prob; + }; + + } // namespace detail + + template <class RealType, class Policy> + inline RealType quantile(const skew_normal_distribution<RealType, Policy>& dist, const RealType& p) + { + const RealType scale = dist.scale(); + const RealType location = dist.location(); + const RealType shape = dist.shape(); + + static const char* function = "boost::math::quantile(const skew_normal_distribution<%1%>&, %1%)"; + + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + if(false == detail::check_probability(function, p, &result, Policy())) + return result; + + // compute initial guess via Cornish-Fisher expansion + RealType x = -boost::math::erfc_inv(2 * p, Policy()) * constants::root_two<RealType>(); + + // avoid unnecessary computations if there is no skew + if(shape != 0) + { + const RealType skew = skewness(dist); + const RealType exk = kurtosis_excess(dist); + + x = x + (x*x-static_cast<RealType>(1))*skew/static_cast<RealType>(6) + + x*(x*x-static_cast<RealType>(3))*exk/static_cast<RealType>(24) + - x*(static_cast<RealType>(2)*x*x-static_cast<RealType>(5))*skew*skew/static_cast<RealType>(36); + } // if(shape != 0) + + result = standard_deviation(dist)*x+mean(dist); + + // handle special case of non-skew normal distribution + if(shape == 0) + return result; + + // refine the result by numerically searching the root of (p-cdf) + + const RealType search_min = range(dist).first; + const RealType search_max = range(dist).second; + + const int get_digits = policies::digits<RealType, Policy>();// get digits from policy, + boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations. + + result = tools::newton_raphson_iterate(detail::skew_normal_quantile_functor<RealType, Policy>(dist, p), result, + search_min, search_max, get_digits, m); + + return result; + } // quantile + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<skew_normal_distribution<RealType, Policy>, RealType>& c) + { + const RealType scale = c.dist.scale(); + const RealType location = c.dist.location(); + const RealType shape = c.dist.shape(); + + static const char* function = "boost::math::quantile(const complement(skew_normal_distribution<%1%>&), %1%)"; + RealType result = 0; + if(false == detail::check_scale(function, scale, &result, Policy())) + return result; + if(false == detail::check_location(function, location, &result, Policy())) + return result; + if(false == detail::check_skew_normal_shape(function, shape, &result, Policy())) + return result; + RealType q = c.param; + if(false == detail::check_probability(function, q, &result, Policy())) + return result; + + skew_normal_distribution<RealType, Policy> D(-location, scale, -shape); + + result = -quantile(D, q); + + return result; + } // quantile + + +} // namespace math +} // namespace boost + +// This include must be at the end, *after* the accessors +// for this distribution have been defined, in order to +// keep compilers that support two-phase lookup happy. +#include <boost/math/distributions/detail/derived_accessors.hpp> + +#endif // BOOST_STATS_SKEW_NORMAL_HPP + + diff --git a/boost/math/distributions/triangular.hpp b/boost/math/distributions/triangular.hpp index ee607e1776..735d20235c 100644 --- a/boost/math/distributions/triangular.hpp +++ b/boost/math/distributions/triangular.hpp @@ -470,7 +470,7 @@ namespace boost{ namespace math RealType mode = dist.mode(); RealType upper = dist.upper(); RealType result = 0; // of checks. - if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy())) + if(false == boost::math::detail::check_triangular(function,lower, mode, upper, &result, Policy())) { return result; } diff --git a/boost/math/policies/error_handling.hpp b/boost/math/policies/error_handling.hpp index 8e1b96d6f5..d4306c6006 100644 --- a/boost/math/policies/error_handling.hpp +++ b/boost/math/policies/error_handling.hpp @@ -12,6 +12,7 @@ #include <iomanip> #include <string> #include <cerrno> +#include <complex> #include <boost/config/no_tr1/cmath.hpp> #include <stdexcept> #include <boost/math/tools/config.hpp> @@ -576,6 +577,15 @@ inline bool check_overflow(T val, R* result, const char* function, const Policy& return false; } template <class R, class T, class Policy> +inline bool check_overflow(std::complex<T> val, R* result, const char* function, const Policy& pol) +{ + typedef typename R::value_type r_type; + r_type re, im; + bool r = check_overflow<r_type>(val.real(), &re, function, pol) || check_overflow<r_type>(val.imag(), &im, function, pol); + *result = R(re, im); + return r; +} +template <class R, class T, class Policy> inline bool check_underflow(T val, R* result, const char* function, const Policy& pol) { if((val != 0) && (static_cast<R>(val) == 0)) @@ -586,6 +596,15 @@ inline bool check_underflow(T val, R* result, const char* function, const Policy return false; } template <class R, class T, class Policy> +inline bool check_underflow(std::complex<T> val, R* result, const char* function, const Policy& pol) +{ + typedef typename R::value_type r_type; + r_type re, im; + bool r = check_underflow<r_type>(val.real(), &re, function, pol) || check_underflow<r_type>(val.imag(), &im, function, pol); + *result = R(re, im); + return r; +} +template <class R, class T, class Policy> inline bool check_denorm(T val, R* result, const char* function, const Policy& pol) { BOOST_MATH_STD_USING @@ -596,14 +615,29 @@ inline bool check_denorm(T val, R* result, const char* function, const Policy& p } return false; } +template <class R, class T, class Policy> +inline bool check_denorm(std::complex<T> val, R* result, const char* function, const Policy& pol) +{ + typedef typename R::value_type r_type; + r_type re, im; + bool r = check_denorm<r_type>(val.real(), &re, function, pol) || check_denorm<r_type>(val.imag(), &im, function, pol); + *result = R(re, im); + return r; +} // Default instantiations with ignore_error policy. template <class R, class T> inline bool check_overflow(T /* val */, R* /* result */, const char* /* function */, const overflow_error<ignore_error>&){ return false; } template <class R, class T> +inline bool check_overflow(std::complex<T> /* val */, R* /* result */, const char* /* function */, const overflow_error<ignore_error>&){ return false; } +template <class R, class T> inline bool check_underflow(T /* val */, R* /* result */, const char* /* function */, const underflow_error<ignore_error>&){ return false; } template <class R, class T> +inline bool check_underflow(std::complex<T> /* val */, R* /* result */, const char* /* function */, const underflow_error<ignore_error>&){ return false; } +template <class R, class T> inline bool check_denorm(T /* val */, R* /* result*/, const char* /* function */, const denorm_error<ignore_error>&){ return false; } +template <class R, class T> +inline bool check_denorm(std::complex<T> /* val */, R* /* result*/, const char* /* function */, const denorm_error<ignore_error>&){ return false; } } // namespace detail diff --git a/boost/math/special_functions.hpp b/boost/math/special_functions.hpp index f31a0694c9..00fe866f4a 100644 --- a/boost/math/special_functions.hpp +++ b/boost/math/special_functions.hpp @@ -55,5 +55,7 @@ #include <boost/math/special_functions/trunc.hpp> #include <boost/math/special_functions/pow.hpp> #include <boost/math/special_functions/next.hpp> +#include <boost/math/special_functions/owens_t.hpp> +#include <boost/math/special_functions/hankel.hpp> #endif // BOOST_MATH_SPECIAL_FUNCTIONS_HPP diff --git a/boost/math/special_functions/detail/bessel_i0.hpp b/boost/math/special_functions/detail/bessel_i0.hpp index 2c129facc7..7dc65d1a1b 100644 --- a/boost/math/special_functions/detail/bessel_i0.hpp +++ b/boost/math/special_functions/detail/bessel_i0.hpp @@ -21,8 +21,38 @@ namespace boost { namespace math { namespace detail{ template <typename T> +T bessel_i0(T x); + +template <class T> +struct bessel_i0_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_i0(T(1)); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T> +const typename bessel_i0_initializer<T>::init bessel_i0_initializer<T>::initializer; + +template <typename T> T bessel_i0(T x) { + bessel_i0_initializer<T>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2335582639474375249e+15)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -5.5050369673018427753e+14)), diff --git a/boost/math/special_functions/detail/bessel_i1.hpp b/boost/math/special_functions/detail/bessel_i1.hpp index aa4596cfd4..47f1b79883 100644 --- a/boost/math/special_functions/detail/bessel_i1.hpp +++ b/boost/math/special_functions/detail/bessel_i1.hpp @@ -21,8 +21,39 @@ namespace boost { namespace math { namespace detail{ template <typename T> +T bessel_i1(T x); + +template <class T> +struct bessel_i1_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_i1(T(1)); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T> +const typename bessel_i1_initializer<T>::init bessel_i1_initializer<T>::initializer; + +template <typename T> T bessel_i1(T x) { + + bessel_i1_initializer<T>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4577180278143463643e+15)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.7732037840791591320e+14)), diff --git a/boost/math/special_functions/detail/bessel_j0.hpp b/boost/math/special_functions/detail/bessel_j0.hpp index ee25d46f61..a07052d73e 100644 --- a/boost/math/special_functions/detail/bessel_j0.hpp +++ b/boost/math/special_functions/detail/bessel_j0.hpp @@ -22,8 +22,38 @@ namespace boost { namespace math { namespace detail{ template <typename T> +T bessel_j0(T x); + +template <class T> +struct bessel_j0_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_j0(T(1)); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T> +const typename bessel_j0_initializer<T>::init bessel_j0_initializer<T>::initializer; + +template <typename T> T bessel_j0(T x) { + bessel_j0_initializer<T>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -4.1298668500990866786e+11)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.7282507878605942706e+10)), diff --git a/boost/math/special_functions/detail/bessel_j1.hpp b/boost/math/special_functions/detail/bessel_j1.hpp index 3db2503ff6..09d862c240 100644 --- a/boost/math/special_functions/detail/bessel_j1.hpp +++ b/boost/math/special_functions/detail/bessel_j1.hpp @@ -22,8 +22,38 @@ namespace boost { namespace math{ namespace detail{ template <typename T> +T bessel_j1(T x); + +template <class T> +struct bessel_j1_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_j1(T(1)); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T> +const typename bessel_j1_initializer<T>::init bessel_j1_initializer<T>::initializer; + +template <typename T> T bessel_j1(T x) { + bessel_j1_initializer<T>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -1.4258509801366645672e+11)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 6.6781041261492395835e+09)), diff --git a/boost/math/special_functions/detail/bessel_jy.hpp b/boost/math/special_functions/detail/bessel_jy.hpp index 19f951ab70..d60dda2d41 100644 --- a/boost/math/special_functions/detail/bessel_jy.hpp +++ b/boost/math/special_functions/detail/bessel_jy.hpp @@ -215,8 +215,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) Cr = br + a / temp; Dr = br; Di = bi; - //std::cout << "C = " << Cr << " " << Ci << std::endl; - //std::cout << "D = " << Dr << " " << Di << std::endl; if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; } if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; } temp = Dr * Dr + Di * Di; @@ -227,7 +225,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) temp = fr; fr = temp * delta_r - fi * delta_i; fi = temp * delta_i + fi * delta_r; - //std::cout << fr << " " << fi << std::endl; for (k = 2; k < policies::get_max_series_iterations<Policy>(); k++) { a = k - 0.5f; @@ -239,8 +236,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) Ci = bi - a * Ci / temp; Dr = br + a * Dr; Di = bi + a * Di; - //std::cout << "C = " << Cr << " " << Ci << std::endl; - //std::cout << "D = " << Dr << " " << Di << std::endl; if (fabs(Cr) + fabs(Ci) < tiny) { Cr = tiny; } if (fabs(Dr) + fabs(Di) < tiny) { Dr = tiny; } temp = Dr * Dr + Di * Di; @@ -253,7 +248,6 @@ int CF2_jy(T v, T x, T* p, T* q, const Policy& pol) fi = temp * delta_i + fi * delta_r; if (fabs(delta_r - 1) + fabs(delta_i) < tolerance) break; - //std::cout << fr << " " << fi << std::endl; } policies::check_series_iterations<T>("boost::math::bessel_jy<%1%>(%1%,%1%) in CF2_jy", k, pol); *p = fr; @@ -491,6 +485,16 @@ int bessel_jy(T v, T x, T* J, T* Y, int kind, const Policy& pol) 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 diff --git a/boost/math/special_functions/detail/bessel_k0.hpp b/boost/math/special_functions/detail/bessel_k0.hpp index 81407dab10..e209168e87 100644 --- a/boost/math/special_functions/detail/bessel_k0.hpp +++ b/boost/math/special_functions/detail/bessel_k0.hpp @@ -22,10 +22,40 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> +T bessel_k0(T x, const Policy&); + +template <class T, class Policy> +struct bessel_k0_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_k0(T(1), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename bessel_k0_initializer<T, Policy>::init bessel_k0_initializer<T, Policy>::initializer; + +template <typename T, typename Policy> T bessel_k0(T x, const Policy& pol) { BOOST_MATH_INSTRUMENT_CODE(x); + bessel_k0_initializer<T, Policy>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2.4708152720399552679e+03)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.9169059852270512312e+03)), diff --git a/boost/math/special_functions/detail/bessel_k1.hpp b/boost/math/special_functions/detail/bessel_k1.hpp index 225209f7ba..0d17cd3057 100644 --- a/boost/math/special_functions/detail/bessel_k1.hpp +++ b/boost/math/special_functions/detail/bessel_k1.hpp @@ -22,8 +22,38 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> +T bessel_k1(T x, const Policy&); + +template <class T, class Policy> +struct bessel_k1_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_k1(T(1), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename bessel_k1_initializer<T, Policy>::init bessel_k1_initializer<T, Policy>::initializer; + +template <typename T, typename Policy> T bessel_k1(T x, const Policy& pol) { + bessel_k1_initializer<T, Policy>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -2.2149374878243304548e+06)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 7.1938920065420586101e+05)), diff --git a/boost/math/special_functions/detail/bessel_y0.hpp b/boost/math/special_functions/detail/bessel_y0.hpp index e23f861bf0..289bda5f18 100644 --- a/boost/math/special_functions/detail/bessel_y0.hpp +++ b/boost/math/special_functions/detail/bessel_y0.hpp @@ -24,8 +24,38 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> +T bessel_y0(T x, const Policy&); + +template <class T, class Policy> +struct bessel_y0_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_y0(T(1), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename bessel_y0_initializer<T, Policy>::init bessel_y0_initializer<T, Policy>::initializer; + +template <typename T, typename Policy> T bessel_y0(T x, const Policy& pol) { + bessel_y0_initializer<T, Policy>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1.0723538782003176831e+11)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -8.3716255451260504098e+09)), diff --git a/boost/math/special_functions/detail/bessel_y1.hpp b/boost/math/special_functions/detail/bessel_y1.hpp index b85e7011ea..caf09ffd26 100644 --- a/boost/math/special_functions/detail/bessel_y1.hpp +++ b/boost/math/special_functions/detail/bessel_y1.hpp @@ -24,8 +24,38 @@ namespace boost { namespace math { namespace detail{ template <typename T, typename Policy> +T bessel_y1(T x, const Policy&); + +template <class T, class Policy> +struct bessel_y1_initializer +{ + struct init + { + init() + { + do_init(); + } + static void do_init() + { + bessel_y1(T(1), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename bessel_y1_initializer<T, Policy>::init bessel_y1_initializer<T, Policy>::initializer; + +template <typename T, typename Policy> T bessel_y1(T x, const Policy& pol) { + bessel_y1_initializer<T, Policy>::force_instantiate(); + static const T P1[] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.0535726612579544093e+13)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 5.4708611716525426053e+12)), diff --git a/boost/math/special_functions/detail/erf_inv.hpp b/boost/math/special_functions/detail/erf_inv.hpp index f2f625f991..d51db9d52f 100644 --- a/boost/math/special_functions/detail/erf_inv.hpp +++ b/boost/math/special_functions/detail/erf_inv.hpp @@ -322,12 +322,47 @@ T erf_inv_imp(const T& p, const T& q, const Policy& pol, const boost::mpl::int_< return result; } +template <class T, class Policy> +struct erf_inv_initializer +{ + struct init + { + init() + { + do_init(); + } + 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) + 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) + 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) + boost::math::erfc_inv(static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 1e-900)), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename erf_inv_initializer<T, Policy>::init erf_inv_initializer<T, Policy>::initializer; + } // namespace detail template <class T, class Policy> typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol) { typedef typename tools::promote_args<T>::type result_type; + // // Begin by testing for domain errors, and other special cases: // @@ -378,6 +413,8 @@ typename tools::promote_args<T>::type erfc_inv(T z, const Policy& pol) policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; + detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate(); + // // And get the result, negating where required: // @@ -389,6 +426,7 @@ template <class T, class Policy> typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol) { typedef typename tools::promote_args<T>::type result_type; + // // Begin by testing for domain errors, and other special cases: // @@ -445,6 +483,8 @@ typename tools::promote_args<T>::type erf_inv(T z, const Policy& pol) // precision internally if it's appropriate: // typedef typename policies::evaluation<result_type, Policy>::type eval_type; + + detail::erf_inv_initializer<eval_type, forwarding_policy>::force_instantiate(); // // And get the result, negating where required: // diff --git a/boost/math/special_functions/detail/iconv.hpp b/boost/math/special_functions/detail/iconv.hpp index 8916eaed1d..4256ffcc88 100644 --- a/boost/math/special_functions/detail/iconv.hpp +++ b/boost/math/special_functions/detail/iconv.hpp @@ -25,7 +25,7 @@ template <class T, class Policy> inline int iconv_imp(T v, Policy const& pol, mpl::false_ const&) { BOOST_MATH_STD_USING - return iround(v); + return iround(v, pol); } template <class T, class Policy> diff --git a/boost/math/special_functions/detail/igamma_large.hpp b/boost/math/special_functions/detail/igamma_large.hpp index f9a810c489..eb3d4ba93e 100644 --- a/boost/math/special_functions/detail/igamma_large.hpp +++ b/boost/math/special_functions/detail/igamma_large.hpp @@ -759,7 +759,6 @@ T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<113> const *) return result; } - } // namespace detail } // namespace math } // namespace math diff --git a/boost/math/special_functions/detail/lanczos_sse2.hpp b/boost/math/special_functions/detail/lanczos_sse2.hpp index 6a3f3e5347..f8846bf376 100644 --- a/boost/math/special_functions/detail/lanczos_sse2.hpp +++ b/boost/math/special_functions/detail/lanczos_sse2.hpp @@ -13,7 +13,7 @@ #include <emmintrin.h> #if defined(__GNUC__) || defined(__PGI) -#define ALIGN16 __attribute__((aligned(16))) +#define ALIGN16 __attribute__((__aligned__(16))) #else #define ALIGN16 __declspec(align(16)) #endif @@ -194,8 +194,11 @@ inline double lanczos13m53::lanczos_sum_expG_scaled<double>(const double& x) } // namespace math } // namespace boost +#undef ALIGN16 + #endif // BOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS + diff --git a/boost/math/special_functions/detail/lgamma_small.hpp b/boost/math/special_functions/detail/lgamma_small.hpp index 526a573583..ec28ed2adf 100644 --- a/boost/math/special_functions/detail/lgamma_small.hpp +++ b/boost/math/special_functions/detail/lgamma_small.hpp @@ -15,6 +15,14 @@ namespace boost{ namespace math{ namespace detail{ // +// These need forward declaring to keep GCC happy: +// +template <class T, class Policy, class Lanczos> +T gamma_imp(T z, const Policy& pol, const Lanczos& l); +template <class T, class Policy> +T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l); + +// // lgamma for small arguments: // template <class T, class Policy, class Lanczos> diff --git a/boost/math/special_functions/digamma.hpp b/boost/math/special_functions/digamma.hpp index 0dbddc77e2..1268b64dc9 100644 --- a/boost/math/special_functions/digamma.hpp +++ b/boost/math/special_functions/digamma.hpp @@ -407,6 +407,31 @@ T digamma_imp(T x, const Tag* t, const Policy& pol) return result; } +// +// Initializer: ensure all our constants are initialized prior to the first call of main: +// +template <class T, class Policy> +struct digamma_initializer +{ + struct init + { + init() + { + boost::math::digamma(T(1.5), Policy()); + boost::math::digamma(T(500), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename digamma_initializer<T, Policy>::init digamma_initializer<T, Policy>::initializer; + } // namespace detail template <class T, class Policy> @@ -433,6 +458,9 @@ inline typename tools::promote_args<T>::type >::type >::type tag_type; + // Force initialization of constants: + detail::digamma_initializer<result_type, Policy>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, Policy>(detail::digamma_imp( static_cast<value_type>(x), static_cast<const tag_type*>(0), pol), "boost::math::digamma<%1%>(%1%)"); diff --git a/boost/math/special_functions/erf.hpp b/boost/math/special_functions/erf.hpp index 1abb59177f..e67332a61a 100644 --- a/boost/math/special_functions/erf.hpp +++ b/boost/math/special_functions/erf.hpp @@ -978,6 +978,59 @@ T erf_imp(T z, bool invert, const Policy& pol, const mpl::int_<113>& t) return result; } // template <class T, class Lanczos>T erf_imp(T z, bool invert, const Lanczos& l, const mpl::int_<113>& t) +template <class T, class Policy, class tag> +struct erf_initializer +{ + struct init + { + init() + { + do_init(tag()); + } + static void do_init(const mpl::int_<0>&){} + static void do_init(const mpl::int_<53>&) + { + boost::math::erf(static_cast<T>(1e-12), Policy()); + boost::math::erf(static_cast<T>(0.25), Policy()); + boost::math::erf(static_cast<T>(1.25), Policy()); + boost::math::erf(static_cast<T>(2.25), Policy()); + boost::math::erf(static_cast<T>(4.25), Policy()); + boost::math::erf(static_cast<T>(5.25), Policy()); + } + static void do_init(const mpl::int_<64>&) + { + boost::math::erf(static_cast<T>(1e-12), Policy()); + boost::math::erf(static_cast<T>(0.25), Policy()); + boost::math::erf(static_cast<T>(1.25), Policy()); + boost::math::erf(static_cast<T>(2.25), Policy()); + boost::math::erf(static_cast<T>(4.25), Policy()); + boost::math::erf(static_cast<T>(5.25), Policy()); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::erf(static_cast<T>(1e-22), Policy()); + boost::math::erf(static_cast<T>(0.25), Policy()); + boost::math::erf(static_cast<T>(1.25), Policy()); + boost::math::erf(static_cast<T>(2.125), Policy()); + boost::math::erf(static_cast<T>(2.75), Policy()); + boost::math::erf(static_cast<T>(3.25), Policy()); + boost::math::erf(static_cast<T>(5.25), Policy()); + boost::math::erf(static_cast<T>(7.25), Policy()); + boost::math::erf(static_cast<T>(11.25), Policy()); + boost::math::erf(static_cast<T>(12.5), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy, class tag> +const typename erf_initializer<T, Policy, tag>::init erf_initializer<T, Policy, tag>::initializer; + } // namespace detail template <class T, class Policy> @@ -1017,6 +1070,8 @@ inline typename tools::promote_args<T>::type erf(T z, const Policy& /* pol */) BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name()); + detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp( static_cast<value_type>(z), false, @@ -1061,6 +1116,8 @@ inline typename tools::promote_args<T>::type erfc(T z, const Policy& /* pol */) BOOST_MATH_INSTRUMENT_CODE("tag_type = " << typeid(tag_type).name()); + detail::erf_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); // Force constants to be initialized before main + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::erf_imp( static_cast<value_type>(z), true, diff --git a/boost/math/special_functions/expint.hpp b/boost/math/special_functions/expint.hpp index 06bd78ff48..1c86d282fa 100644 --- a/boost/math/special_functions/expint.hpp +++ b/boost/math/special_functions/expint.hpp @@ -1473,6 +1473,94 @@ T expint_i_imp(T z, const Policy& pol, const mpl::int_<113>& tag) return result; } +template <class T, class Policy, class tag> +struct expint_i_initializer +{ + struct init + { + init() + { + do_init(tag()); + } + static void do_init(const mpl::int_<0>&){} + static void do_init(const mpl::int_<53>&) + { + boost::math::expint(T(5)); + boost::math::expint(T(7)); + boost::math::expint(T(18)); + boost::math::expint(T(38)); + boost::math::expint(T(45)); + } + static void do_init(const mpl::int_<64>&) + { + boost::math::expint(T(5)); + boost::math::expint(T(7)); + boost::math::expint(T(18)); + boost::math::expint(T(38)); + boost::math::expint(T(45)); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::expint(T(5)); + boost::math::expint(T(7)); + boost::math::expint(T(17)); + boost::math::expint(T(25)); + boost::math::expint(T(40)); + boost::math::expint(T(50)); + boost::math::expint(T(80)); + boost::math::expint(T(200)); + boost::math::expint(T(220)); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy, class tag> +const typename expint_i_initializer<T, Policy, tag>::init expint_i_initializer<T, Policy, tag>::initializer; + +template <class T, class Policy, class tag> +struct expint_1_initializer +{ + struct init + { + init() + { + do_init(tag()); + } + static void do_init(const mpl::int_<0>&){} + static void do_init(const mpl::int_<53>&) + { + boost::math::expint(1, T(0.5)); + boost::math::expint(1, T(2)); + } + static void do_init(const mpl::int_<64>&) + { + boost::math::expint(1, T(0.5)); + boost::math::expint(1, T(2)); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::expint(1, T(0.5)); + boost::math::expint(1, T(2)); + boost::math::expint(1, T(6)); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy, class tag> +const typename expint_1_initializer<T, Policy, tag>::init expint_1_initializer<T, Policy, tag>::initializer; + template <class T, class Policy> inline typename tools::promote_args<T>::type expint_forwarder(T z, const Policy& /*pol*/, mpl::true_ const&) @@ -1504,6 +1592,8 @@ inline typename tools::promote_args<T>::type >::type >::type tag_type; + expint_i_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_i_imp( static_cast<value_type>(z), forwarding_policy(), @@ -1550,6 +1640,8 @@ inline typename tools::promote_args<T>::type >::type >::type tag_type; + detail::expint_1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expint_imp( n, static_cast<value_type>(z), diff --git a/boost/math/special_functions/expm1.hpp b/boost/math/special_functions/expm1.hpp index 345220fcee..9ff2541fb1 100644 --- a/boost/math/special_functions/expm1.hpp +++ b/boost/math/special_functions/expm1.hpp @@ -65,34 +65,37 @@ namespace detail expm1_series& operator=(const expm1_series&); }; -template <class T, bool b = boost::is_pod<T>::value> -struct expm1_init_on_startup +template <class T, class Policy, class tag> +struct expm1_initializer { struct init { init() { - boost::math::expm1(T(0.5f)); + do_init(tag()); } - void do_nothing()const{} + template <int N> + static void do_init(const mpl::int_<N>&){} + static void do_init(const mpl::int_<64>&) + { + expm1(T(0.5)); + } + static void do_init(const mpl::int_<113>&) + { + expm1(T(0.5)); + } + void force_instantiate()const{} }; - - static void do_nothing() + static const init initializer; + static void force_instantiate() { - initializer.do_nothing(); + initializer.force_instantiate(); } - - static const init initializer; }; -template <class T, bool b> -const typename expm1_init_on_startup<T, b>::init expm1_init_on_startup<T, b>::initializer; +template <class T, class Policy, class tag> +const typename expm1_initializer<T, Policy, tag>::init expm1_initializer<T, Policy, tag>::initializer; -template <class T> -struct expm1_init_on_startup<T, true> -{ - static void do_nothing(){} -}; // // Algorithm expm1 is part of C99, but is not yet provided by many compilers. // @@ -133,8 +136,6 @@ T expm1_imp(T x, const mpl::int_<53>&, const P& pol) { BOOST_MATH_STD_USING - expm1_init_on_startup<T>::do_nothing(); - T a = fabs(x); if(a > T(0.5L)) { @@ -162,8 +163,6 @@ T expm1_imp(T x, const mpl::int_<64>&, const P& pol) { BOOST_MATH_STD_USING - expm1_init_on_startup<T>::do_nothing(); - T a = fabs(x); if(a > T(0.5L)) { @@ -207,8 +206,6 @@ T expm1_imp(T x, const mpl::int_<113>&, const P& pol) { BOOST_MATH_STD_USING - expm1_init_on_startup<T>::do_nothing(); - T a = fabs(x); if(a > T(0.5L)) { @@ -287,6 +284,8 @@ inline typename tools::promote_args<T>::type expm1(T x, const Policy& /* pol */) >::type >::type tag_type; + detail::expm1_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::expm1_imp( static_cast<value_type>(x), tag_type(), forwarding_policy()), "boost::math::expm1<%1%>(%1%)"); diff --git a/boost/math/special_functions/gamma.hpp b/boost/math/special_functions/gamma.hpp index 1ae965f18c..86d15b7f2a 100644 --- a/boost/math/special_functions/gamma.hpp +++ b/boost/math/special_functions/gamma.hpp @@ -1258,6 +1258,101 @@ inline typename tools::promote_args<T>::type return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)"); } +template <class T, class Policy> +struct igamma_initializer +{ + struct init + { + init() + { + 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::less_equal<precision_type, mpl::int_<53> >, + mpl::int_<53>, + typename mpl::if_< + mpl::less_equal<precision_type, mpl::int_<64> >, + mpl::int_<64>, + mpl::int_<113> + >::type + >::type + >::type tag_type; + + do_init(tag_type()); + } + template <int N> + static void do_init(const mpl::int_<N>&) + { + 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{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer; + +template <class T, class Policy> +struct lgamma_initializer +{ + struct init + { + init() + { + typedef typename policies::precision<T, Policy>::type precision_type; + typedef typename mpl::if_< + mpl::and_< + mpl::less_equal<precision_type, mpl::int_<64> >, + mpl::greater<precision_type, mpl::int_<0> > + >, + mpl::int_<64>, + typename mpl::if_< + mpl::and_< + mpl::less_equal<precision_type, mpl::int_<113> >, + mpl::greater<precision_type, mpl::int_<0> > + >, + mpl::int_<113>, mpl::int_<0> >::type + >::type tag_type; + do_init(tag_type()); + } + static void do_init(const mpl::int_<64>&) + { + boost::math::lgamma(static_cast<T>(2.5), Policy()); + boost::math::lgamma(static_cast<T>(1.25), Policy()); + boost::math::lgamma(static_cast<T>(1.75), Policy()); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::lgamma(static_cast<T>(2.5), Policy()); + boost::math::lgamma(static_cast<T>(1.25), Policy()); + boost::math::lgamma(static_cast<T>(1.5), Policy()); + boost::math::lgamma(static_cast<T>(1.75), Policy()); + } + static void do_init(const mpl::int_<0>&) + { + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy> +const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer; + template <class T1, class T2, class Policy> inline typename tools::promote_args<T1, T2>::type tgamma(T1 a, T2 z, const Policy&, const mpl::false_) @@ -1272,6 +1367,9 @@ inline typename tools::promote_args<T1, T2>::type policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; + + igamma_initializer<value_type, forwarding_policy>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>( detail::gamma_incomplete_imp(static_cast<value_type>(a), static_cast<value_type>(z), false, true, @@ -1285,6 +1383,7 @@ inline typename tools::promote_args<T1, T2>::type return tgamma(a, z, policies::policy<>(), tag); } + } // namespace detail template <class T> @@ -1308,6 +1407,9 @@ inline typename tools::promote_args<T>::type policies::promote_double<false>, policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; + + detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)"); } @@ -1395,6 +1497,8 @@ inline typename tools::promote_args<T1, T2>::type policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; + detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>( detail::gamma_incomplete_imp(static_cast<value_type>(a), static_cast<value_type>(z), false, false, @@ -1424,6 +1528,8 @@ inline typename tools::promote_args<T1, T2>::type policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; + detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>( detail::gamma_incomplete_imp(static_cast<value_type>(a), static_cast<value_type>(z), true, true, @@ -1453,6 +1559,8 @@ inline typename tools::promote_args<T1, T2>::type policies::discrete_quantile<>, policies::assert_undefined<> >::type forwarding_policy; + detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>( detail::gamma_incomplete_imp(static_cast<value_type>(a), static_cast<value_type>(z), true, false, diff --git a/boost/math/special_functions/hankel.hpp b/boost/math/special_functions/hankel.hpp new file mode 100644 index 0000000000..bc3fc2d742 --- /dev/null +++ b/boost/math/special_functions/hankel.hpp @@ -0,0 +1,178 @@ +// Copyright John Maddock 2012. +// 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_HANKEL_HPP +#define BOOST_MATH_HANKEL_HPP + +#include <boost/math/special_functions/bessel.hpp> + +namespace boost{ namespace math{ + +namespace detail{ + +template <class T, class Policy> +std::complex<T> hankel_imp(T v, T x, const bessel_no_int_tag&, const Policy& pol, int sign) +{ + BOOST_MATH_STD_USING + static const char* function = "boost::math::cyl_hankel_1<%1%>(%1%,%1%)"; + + if(x < 0) + { + bool isint_v = floor(v) == v; + T j, y; + bessel_jy(v, -x, &j, &y, need_j | need_y, pol); + std::complex<T> cx(x), cv(v); + std::complex<T> j_result, y_result; + if(isint_v) + { + int s = (iround(j) & 1) ? -1 : 1; + j_result = j * s; + y_result = T(s) * (y - (2 / constants::pi<T>()) * (log(-x) - log(cx)) * j); + } + else + { + j_result = pow(cx, v) * pow(-cx, -v) * j; + T p1 = pow(-x, v); + std::complex<T> p2 = pow(cx, v); + y_result = p1 * y / p2 + + (p2 / p1 - p1 / p2) * j / tan(constants::pi<T>() * v); + } + // multiply y_result by i: + y_result = std::complex<T>(-sign * y_result.imag(), sign * y_result.real()); + return j_result + y_result; + } + + if(x == 0) + { + if(v == 0) + { + // J is 1, Y is -INF + return std::complex<T>(1, sign * -policies::raise_overflow_error<T>(function, 0, pol)); + } + else + { + // At least one of J and Y is complex infinity: + return std::complex<T>(policies::raise_overflow_error<T>(function, 0, pol), sign * policies::raise_overflow_error<T>(function, 0, pol)); + } + } + + T j, y; + bessel_jy(v, x, &j, &y, need_j | need_y, pol); + return std::complex<T>(j, sign * y); +} + +template <class T, class Policy> +std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign); + +template <class T, class Policy> +inline std::complex<T> hankel_imp(T v, T x, const bessel_maybe_int_tag&, const Policy& pol, int sign) +{ + BOOST_MATH_STD_USING // ADL of std names. + int ival = detail::iconv(v, pol); + if(0 == v - ival) + { + return hankel_imp(ival, x, bessel_int_tag(), pol, sign); + } + return hankel_imp(v, x, bessel_no_int_tag(), pol, sign); +} + +template <class T, class Policy> +inline std::complex<T> hankel_imp(int v, T x, const bessel_int_tag&, const Policy& pol, int sign) +{ + BOOST_MATH_STD_USING + if((std::abs(v < 200)) && (x > 0)) + return std::complex<T>(bessel_jn(v, x, pol), sign * bessel_yn(v, x, pol)); + return hankel_imp(static_cast<T>(v), x, bessel_no_int_tag(), pol, sign); +} + +template <class T, class Policy> +inline std::complex<T> sph_hankel_imp(T v, T x, const Policy& pol, int sign) +{ + BOOST_MATH_STD_USING + return constants::root_half_pi<T>() * hankel_imp(v + 0.5f, x, bessel_no_int_tag(), pol, sign) / sqrt(std::complex<T>(x)); +} + +} // namespace detail + +template <class T1, class T2, class Policy> +inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; + typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, 1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)"); +} + +template <class T1, class T2> +inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x) +{ + return cyl_hankel_1(v, x, policies::policy<>()); +} + +template <class T1, class T2, class Policy> +inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; + typedef typename detail::bessel_traits<T1, T2, Policy>::optimisation_tag tag_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::hankel_imp<value_type>(v, static_cast<value_type>(x), tag_type(), pol, -1), "boost::math::cyl_hankel_1<%1%>(%1%,%1%)"); +} + +template <class T1, class T2> +inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x) +{ + return cyl_hankel_2(v, x, policies::policy<>()); +} + +template <class T1, class T2, class Policy> +inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy& pol) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), 1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)"); +} + +template <class T1, class T2> +inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x) +{ + return sph_hankel_1(v, x, policies::policy<>()); +} + +template <class T1, class T2, class Policy> +inline std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy& pol) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename detail::bessel_traits<T1, T2, Policy>::result_type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float<false>, + policies::promote_double<false>, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + + return policies::checked_narrowing_cast<std::complex<result_type>, Policy>(detail::sph_hankel_imp<value_type>(static_cast<value_type>(v), static_cast<value_type>(x), forwarding_policy(), -1), "boost::math::sph_hankel_1<%1%>(%1%,%1%)"); +} + +template <class T1, class T2> +inline std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x) +{ + return sph_hankel_2(v, x, policies::policy<>()); +} + +}} // namespaces + +#endif // BOOST_MATH_HANKEL_HPP
\ No newline at end of file diff --git a/boost/math/special_functions/lanczos.hpp b/boost/math/special_functions/lanczos.hpp index 20ff969359..ed891549f1 100644 --- a/boost/math/special_functions/lanczos.hpp +++ b/boost/math/special_functions/lanczos.hpp @@ -34,6 +34,35 @@ namespace boost{ namespace math{ namespace lanczos{ // http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at // http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision. // +// Begin with a small helper to force initialization of constants prior +// to main. This makes the constant initialization thread safe, even +// when called with a user-defined number type. +// +template <class Lanczos, class T> +struct lanczos_initializer +{ + struct init + { + init() + { + T t(1); + Lanczos::lanczos_sum(t); + Lanczos::lanczos_sum_expG_scaled(t); + Lanczos::lanczos_sum_near_1(t); + Lanczos::lanczos_sum_near_2(t); + Lanczos::g(); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; +template <class Lanczos, class T> +typename lanczos_initializer<Lanczos, T>::init const lanczos_initializer<Lanczos, T>::initializer; +// // Lanczos Coefficients for N=6 G=5.581 // Max experimental error (with arbitary precision arithmetic) 9.516e-12 // Generated with compiler: Microsoft Visual C++ version 8.0 on Win32 at Mar 23 2006 @@ -47,6 +76,7 @@ struct lanczos6 : public mpl::int_<35> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[6] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 8706.349592549009182288174442774377925882)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 8523.650341121874633477483696775067709735)), @@ -69,6 +99,7 @@ struct lanczos6 : public mpl::int_<35> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[6] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 32.81244541029783471623665933780748627823)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 32.12388941444332003446077108933558534361)), @@ -92,6 +123,7 @@ struct lanczos6 : public mpl::int_<35> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[5] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 2.044879010930422922760429926121241330235)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, -2.751366405578505366591317846728753993668)), @@ -110,6 +142,7 @@ struct lanczos6 : public mpl::int_<35> template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos6, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[5] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, 5.748142489536043490764289256167080091892)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 35, -7.734074268282457156081021756682138251825)), @@ -143,6 +176,7 @@ struct lanczos11 : public mpl::int_<60> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[11] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 38474670393.31776828316099004518914832218)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 36857665043.51950660081971227404959150474)), @@ -175,6 +209,7 @@ struct lanczos11 : public mpl::int_<60> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[11] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 709811.662581657956893540610814842699825)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 679979.847415722640161734319823103390728)), @@ -208,6 +243,7 @@ struct lanczos11 : public mpl::int_<60> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[10] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 4.005853070677940377969080796551266387954)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, -13.17044315127646469834125159673527183164)), @@ -231,6 +267,7 @@ struct lanczos11 : public mpl::int_<60> template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos11, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[10] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, 19.05889633808148715159575716844556056056)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 60, -62.66183664701721716960978577959655644762)), @@ -269,6 +306,7 @@ struct lanczos13 : public mpl::int_<72> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[13] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 44012138428004.60895436261759919070125699)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 41590453358593.20051581730723108131357995)), @@ -305,6 +343,7 @@ struct lanczos13 : public mpl::int_<72> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[13] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 86091529.53418537217994842267760536134841)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 81354505.17858011242874285785316135398567)), @@ -342,6 +381,7 @@ struct lanczos13 : public mpl::int_<72> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[12] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 4.832115561461656947793029596285626840312)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, -19.86441536140337740383120735104359034688)), @@ -367,6 +407,7 @@ struct lanczos13 : public mpl::int_<72> template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos13, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[12] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, 26.96979819614830698367887026728396466395)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, -110.8705424709385114023884328797900204863)), @@ -382,7 +423,7 @@ struct lanczos13 : public mpl::int_<72> static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 72, -0.9685385411006641478305219367315965391289e-9)), }; T result = 0; - T z = z = 2; + T z = dz + 2; for(unsigned k = 1; k <= sizeof(d)/sizeof(d[0]); ++k) { result += (-d[k-1]*dz)/(z + k*z + k*k - 1); @@ -407,6 +448,7 @@ struct lanczos22 : public mpl::int_<120> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[22] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 46198410803245094237463011094.12173081986)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 43735859291852324413622037436.321513777)), @@ -461,6 +503,7 @@ struct lanczos22 : public mpl::int_<120> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[22] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 6939996264376682180.277485395074954356211)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 6570067992110214451.87201438870245659384)), @@ -516,6 +559,7 @@ struct lanczos22 : public mpl::int_<120> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[21] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 8.318998691953337183034781139546384476554)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, -63.15415991415959158214140353299240638675)), @@ -550,6 +594,7 @@ struct lanczos22 : public mpl::int_<120> template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos22, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[21] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, 75.39272007105208086018421070699575462226)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 120, -572.3481967049935412452681346759966390319)), @@ -830,6 +875,7 @@ struct lanczos17m64 : public mpl::int_<64> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[17] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 553681095419291969.2230556393350368550504)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 731918863887667017.2511276782146694632234)), @@ -874,6 +920,7 @@ struct lanczos17m64 : public mpl::int_<64> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[17] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 2715894658327.717377557655133124376674911)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 3590179526097.912105038525528721129550434)), @@ -919,6 +966,7 @@ struct lanczos17m64 : public mpl::int_<64> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[16] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 4.493645054286536365763334986866616581265)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -16.95716370392468543800733966378143997694)), @@ -948,6 +996,7 @@ struct lanczos17m64 : public mpl::int_<64> template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos17m64, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[16] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, 23.56409085052261327114594781581930373708)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 64, -88.92116338946308797946237246006238652361)), @@ -991,6 +1040,7 @@ struct lanczos24m113 : public mpl::int_<113> template <class T> static T lanczos_sum(const T& z) { + lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[24] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 2029889364934367661624137213253.22102954656825019111612712252027267955023987678816620961507)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 2338599599286656537526273232565.2727349714338768161421882478417543004440597874814359063158)), @@ -1049,6 +1099,7 @@ struct lanczos24m113 : public mpl::int_<113> template <class T> static T lanczos_sum_expG_scaled(const T& z) { + lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T num[24] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 3035162425359883494754.02878223286972654682199012688209026810841953293372712802258398358538)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 3496756894406430103600.16057175075063458536101374170860226963245118484234495645518505519827)), @@ -1108,6 +1159,7 @@ struct lanczos24m113 : public mpl::int_<113> template<class T> static T lanczos_sum_near_1(const T& dz) { + lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[23] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 7.4734083002469026177867421609938203388868806387315406134072298925733950040583068760685908)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, -50.4225805042247530267317342133388132970816607563062253708655085754357843064134941138154171)), @@ -1144,6 +1196,7 @@ struct lanczos24m113 : public mpl::int_<113> template<class T> static T lanczos_sum_near_2(const T& dz) { + lanczos_initializer<lanczos24m113, T>::force_instantiate(); // Ensure our constants get initialized before main() static const T d[23] = { static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, 61.4165001061101455341808888883960361969557848005400286332291451422461117307237198559485365)), static_cast<T>(BOOST_MATH_BIG_CONSTANT(T, 113, -414.372973678657049667308134761613915623353625332248315105320470271523320700386200587519147)), diff --git a/boost/math/special_functions/log1p.hpp b/boost/math/special_functions/log1p.hpp index 9bae7165e4..989bdc21b6 100644 --- a/boost/math/special_functions/log1p.hpp +++ b/boost/math/special_functions/log1p.hpp @@ -258,6 +258,34 @@ T log1p_imp(T const& x, const Policy& pol, const mpl::int_<24>&) return result; } +template <class T, class Policy, class tag> +struct log1p_initializer +{ + struct init + { + init() + { + do_init(tag()); + } + template <int N> + static void do_init(const mpl::int_<N>&){} + static void do_init(const mpl::int_<64>&) + { + boost::math::log1p(static_cast<T>(0.25), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy, class tag> +const typename log1p_initializer<T, Policy, tag>::init log1p_initializer<T, Policy, tag>::initializer; + + } // namespace detail template <class T, class Policy> @@ -286,6 +314,9 @@ inline typename tools::promote_args<T>::type log1p(T x, const Policy&) >::type >::type >::type tag_type; + + detail::log1p_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>( detail::log1p_imp(static_cast<value_type>(x), forwarding_policy(), tag_type()), "boost::math::log1p<%1%>(%1%)"); } diff --git a/boost/math/special_functions/math_fwd.hpp b/boost/math/special_functions/math_fwd.hpp index 14364a3d5c..982cdf7ca3 100644 --- a/boost/math/special_functions/math_fwd.hpp +++ b/boost/math/special_functions/math_fwd.hpp @@ -28,6 +28,7 @@ #include <boost/math/policies/policy.hpp> #include <boost/mpl/comparison.hpp> #include <boost/config/no_tr1/complex.hpp> +#include <complex> #define BOOST_NO_MACRO_EXPAND /**/ @@ -614,6 +615,30 @@ namespace boost template <class T> typename detail::bessel_traits<T, T, policies::policy<> >::result_type sph_neumann(unsigned v, T x); + template <class T1, class T2, class Policy> + std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_1(T1 v, T2 x, const Policy& pol); + + template <class T1, class T2> + std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_1(T1 v, T2 x); + + template <class T1, class T2, class Policy> + std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> cyl_hankel_2(T1 v, T2 x, const Policy& pol); + + template <class T1, class T2> + std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> cyl_hankel_2(T1 v, T2 x); + + template <class T1, class T2, class Policy> + std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_1(T1 v, T2 x, const Policy& pol); + + template <class T1, class T2> + std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_1(T1 v, T2 x); + + template <class T1, class T2, class Policy> + std::complex<typename detail::bessel_traits<T1, T2, Policy>::result_type> sph_hankel_2(T1 v, T2 x, const Policy& pol); + + template <class T1, class T2> + std::complex<typename detail::bessel_traits<T1, T2, policies::policy<> >::result_type> sph_hankel_2(T1 v, T2 x); + template <class T, class Policy> typename tools::promote_args<T>::type sin_pi(T x, const Policy&); @@ -681,6 +706,13 @@ namespace boost template <class T, class Policy> typename tools::promote_args<T>::type zeta(T s, const Policy&); + // Owen's T function: + template <class T1, class T2, class Policy> + typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol); + + template <class T1, class T2> + typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a); + template <class T> typename tools::promote_args<T>::type zeta(T s); @@ -1063,6 +1095,26 @@ namespace boost template <class T> T float_next(const T& a){ return boost::math::float_next(a, Policy()); }\ template <class T> T float_prior(const T& a){ return boost::math::float_prior(a, Policy()); }\ template <class T> T float_distance(const T& a, const T& b){ return boost::math::float_distance(a, b, Policy()); }\ + \ + template <class RT1, class RT2>\ + inline typename boost::math::tools::promote_args<RT1, RT2>::type owens_t(RT1 a, RT2 z){ return boost::math::owens_t(a, z, Policy()); }\ + \ + template <class T1, class T2>\ + inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> cyl_hankel_1(T1 v, T2 x)\ + { return boost::math::cyl_hankel_1(v, x, Policy()); }\ + \ + template <class T1, class T2>\ + inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> cyl_hankel_2(T1 v, T2 x)\ + { return boost::math::cyl_hankel_2(v, x, Policy()); }\ + \ + template <class T1, class T2>\ + inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> sph_hankel_1(T1 v, T2 x)\ + { return boost::math::sph_hankel_1(v, x, Policy()); }\ + \ + template <class T1, class T2>\ + inline std::complex<typename boost::math::detail::bessel_traits<T1, T2, Policy >::result_type> sph_hankel_2(T1 v, T2 x)\ + { return boost::math::sph_hankel_2(v, x, Policy()); }\ + #endif // BOOST_MATH_SPECIAL_MATH_FWD_HPP diff --git a/boost/math/special_functions/nonfinite_num_facets.hpp b/boost/math/special_functions/nonfinite_num_facets.hpp index 9fa61481b5..84d3f1070a 100644 --- a/boost/math/special_functions/nonfinite_num_facets.hpp +++ b/boost/math/special_functions/nonfinite_num_facets.hpp @@ -1,8 +1,9 @@ #ifndef BOOST_MATH_NONFINITE_NUM_FACETS_HPP #define BOOST_MATH_NONFINITE_NUM_FACETS_HPP -// Copyright (c) 2006 Johan Rade -// Copyright 2011 Paul A. Bristow (comments) +// Copyright 2006 Johan Rade +// Copyright 2012 K R Walker +// Copyright 2011, 2012 Paul A. Bristow // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt @@ -28,10 +29,9 @@ #include <boost/math/special_functions/sign.hpp> #ifdef _MSC_VER -# pragma warning(push) -# pragma warning(disable : 4127) // conditional expression is constant. -# pragma warning(disable : 4706) // assignment within conditional expression. -# pragma warning(disable : 4224) // formal parameter 'version' was previously defined as a type. +# pragma warning(push) +# pragma warning(disable : 4127) // conditional expression is constant. +# pragma warning(disable : 4706) // assignment within conditional expression. #endif namespace boost { @@ -65,16 +65,14 @@ namespace boost { protected: virtual OutputIterator do_put( - OutputIterator it, std::ios_base& iosb, - CharType fill, double val) const + OutputIterator it, std::ios_base& iosb, CharType fill, double val) const { put_and_reset_width(it, iosb, fill, val); return it; } virtual OutputIterator do_put( - OutputIterator it, std::ios_base& iosb, - CharType fill, long double val) const + OutputIterator it, std::ios_base& iosb, CharType fill, long double val) const { put_and_reset_width(it, iosb, fill, val); return it; @@ -93,86 +91,136 @@ namespace boost { OutputIterator& it, std::ios_base& iosb, CharType fill, ValType val) const { - switch((boost::math::fpclassify)(val)) { + static const CharType prefix_plus[2] = { '+', '\0' }; + static const CharType prefix_minus[2] = { '-', '\0' }; + static const CharType body_inf[4] = { 'i', 'n', 'f', '\0' }; + static const CharType body_nan[4] = { 'n', 'a', 'n', '\0' }; + static const CharType* null_string = 0; + + switch((boost::math::fpclassify)(val)) + { case FP_INFINITE: if(flags_ & trap_infinity) + { throw std::ios_base::failure("Infinity"); + } else if((boost::math::signbit)(val)) - put_num_and_fill(it, iosb, "-", "inf", fill); + { // negative infinity. + put_num_and_fill(it, iosb, prefix_minus, body_inf, fill, val); + } else if(iosb.flags() & std::ios_base::showpos) - put_num_and_fill(it, iosb, "+", "inf", fill); + { // Explicit "+inf" wanted. + put_num_and_fill(it, iosb, prefix_plus, body_inf, fill, val); + } else - put_num_and_fill(it, iosb, "", "inf", fill); + { // just "inf" wanted. + put_num_and_fill(it, iosb, null_string, body_inf, fill, val); + } break; case FP_NAN: if(flags_ & trap_nan) + { throw std::ios_base::failure("NaN"); + } else if((boost::math::signbit)(val)) - put_num_and_fill(it, iosb, "-", "nan", fill); + { // negative so "-nan". + put_num_and_fill(it, iosb, prefix_minus, body_nan, fill, val); + } else if(iosb.flags() & std::ios_base::showpos) - put_num_and_fill(it, iosb, "+", "nan", fill); + { // explicit "+nan" wanted. + put_num_and_fill(it, iosb, prefix_plus, body_nan, fill, val); + } else - put_num_and_fill(it, iosb, "", "nan", fill); + { // Just "nan". + put_num_and_fill(it, iosb, null_string, body_nan, fill, val); + } break; case FP_ZERO: - if(flags_ & signed_zero) { - if((boost::math::signbit)(val)) - put_num_and_fill(it, iosb, "-", "0", fill); - else if(iosb.flags() & std::ios_base::showpos) - put_num_and_fill(it, iosb, "+", "0", fill); - else - put_num_and_fill(it, iosb, "", "0", fill); + if((flags_ & signed_zero) && ((boost::math::signbit)(val))) + { // Flag set to distinguish between positive and negative zero. + // But string "0" should have stuff after decimal point if setprecision and/or exp format. + + std::basic_ostringstream<CharType> zeros; // Needs to be CharType version. + + // Copy flags, fill, width and precision. + zeros.flags(iosb.flags()); + zeros.unsetf(std::ios::showpos); // Ignore showpos because must be negative. + zeros.precision(iosb.precision()); + //zeros.width is set by put_num_and_fill + zeros.fill(static_cast<char>(fill)); + zeros << ValType(0); + put_num_and_fill(it, iosb, prefix_minus, zeros.str().c_str(), fill, val); } else - put_num_and_fill(it, iosb, "", "0", fill); + { // Output the platform default for positive and negative zero. + put_num_and_fill(it, iosb, null_string, null_string, fill, val); + } break; - default: - it = std::num_put<CharType, OutputIterator>::do_put( - it, iosb, fill, val); + default: // Normal non-zero finite value. + it = std::num_put<CharType, OutputIterator>::do_put(it, iosb, fill, val); break; } } + template<class ValType> void put_num_and_fill( - OutputIterator& it, std::ios_base& iosb, const char* prefix, - const char* body, CharType fill) const + OutputIterator& it, std::ios_base& iosb, const CharType* prefix, + const CharType* body, CharType fill, ValType val) const { - int width = (int)std::strlen(prefix) + (int)std::strlen(body); - std::ios_base::fmtflags adjust - = iosb.flags() & std::ios_base::adjustfield; + int prefix_length = prefix ? (int)std::char_traits<CharType>::length(prefix) : 0; + int body_length = body ? (int)std::char_traits<CharType>::length(body) : 0; + int width = prefix_length + body_length; + std::ios_base::fmtflags adjust = iosb.flags() & std::ios_base::adjustfield; const std::ctype<CharType>& ct = std::use_facet<std::ctype<CharType> >(iosb.getloc()); - if(adjust != std::ios_base::internal && adjust != std::ios_base::left) - put_fill(it, iosb, fill, width); + if(body || prefix) + { // adjust == std::ios_base::right, so leading fill needed. + if(adjust != std::ios_base::internal && adjust != std::ios_base::left) + put_fill(it, iosb, fill, width); + } - while(*prefix) - *it = ct.widen(*(prefix++)); + if(prefix) + { // Adjust width for prefix. + while(*prefix) + *it = *(prefix++); + iosb.width( iosb.width() - prefix_length ); + width -= prefix_length; + } - if(adjust == std::ios_base::internal) - put_fill(it, iosb, fill, width); + if(body) + { // + if(adjust == std::ios_base::internal) + { // Put fill between sign and digits. + put_fill(it, iosb, fill, width); + } + if(iosb.flags() & std::ios_base::uppercase) + { + while(*body) + *it = ct.toupper(*(body++)); + } + else + { + while(*body) + *it = *(body++); + } - if(iosb.flags() & std::ios_base::uppercase) { - while(*body) - *it = ct.toupper(ct.widen(*(body++))); + if(adjust == std::ios_base::left) + put_fill(it, iosb, fill, width); } - else { - while(*body) - *it = ct.widen(*(body++)); + else + { + it = std::num_put<CharType, OutputIterator>::do_put(it, iosb, fill, val); } - - if(adjust == std::ios_base::left) - put_fill(it, iosb, fill, width); } void put_fill( - OutputIterator& it, std::ios_base& iosb, - CharType fill, int width) const - { + OutputIterator& it, std::ios_base& iosb, CharType fill, int width) const + { // Insert fill chars. for(std::streamsize i = iosb.width() - static_cast<std::streamsize>(width); i > 0; --i) *it = fill; } @@ -540,4 +588,5 @@ namespace boost { # pragma warning(pop) #endif -#endif +#endif // BOOST_MATH_NONFINITE_NUM_FACETS_HPP + diff --git a/boost/math/special_functions/owens_t.hpp b/boost/math/special_functions/owens_t.hpp new file mode 100644 index 0000000000..98d6380c39 --- /dev/null +++ b/boost/math/special_functions/owens_t.hpp @@ -0,0 +1,1061 @@ +// (C) Benjamin Sobotta 2012 + +// 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_OWENS_T_HPP +#define BOOST_OWENS_T_HPP + +// Reference: +// Mike Patefield, David Tandy +// FAST AND ACCURATE CALCULATION OF OWEN'S T-FUNCTION +// Journal of Statistical Software, 5 (5), 1-25 + +#ifdef _MSC_VER +# pragma once +#endif + +#include <boost/config/no_tr1/cmath.hpp> +#include <boost/math/special_functions/erf.hpp> +#include <boost/math/special_functions/expm1.hpp> +#include <boost/throw_exception.hpp> +#include <boost/assert.hpp> +#include <boost/math/constants/constants.hpp> +#include <boost/math/tools/big_constant.hpp> + +#include <stdexcept> + +namespace boost +{ + namespace math + { + namespace detail + { + // owens_t_znorm1(x) = P(-oo<Z<=x)-0.5 with Z being normally distributed. + template<typename RealType> + inline RealType owens_t_znorm1(const RealType x) + { + using namespace boost::math::constants; + return erf(x*one_div_root_two<RealType>())*half<RealType>(); + } // RealType owens_t_znorm1(const RealType x) + + // owens_t_znorm2(x) = P(x<=Z<oo) with Z being normally distributed. + template<typename RealType> + inline RealType owens_t_znorm2(const RealType x) + { + using namespace boost::math::constants; + return erfc(x*one_div_root_two<RealType>())*half<RealType>(); + } // RealType owens_t_znorm2(const RealType x) + + // Auxiliary function, it computes an array key that is used to determine + // the specific computation method for Owen's T and the order thereof + // used in owens_t_dispatch. + template<typename RealType> + inline unsigned short owens_t_compute_code(const RealType h, const RealType a) + { + static const RealType hrange[] = + {0.02, 0.06, 0.09, 0.125, 0.26, 0.4, 0.6, 1.6, 1.7, 2.33, 2.4, 3.36, 3.4, 4.8}; + + static const RealType arange[] = {0.025, 0.09, 0.15, 0.36, 0.5, 0.9, 0.99999}; + /* + original select array from paper: + 1, 1, 2,13,13,13,13,13,13,13,13,16,16,16, 9 + 1, 2, 2, 3, 3, 5, 5,14,14,15,15,16,16,16, 9 + 2, 2, 3, 3, 3, 5, 5,15,15,15,15,16,16,16,10 + 2, 2, 3, 5, 5, 5, 5, 7, 7,16,16,16,16,16,10 + 2, 3, 3, 5, 5, 6, 6, 8, 8,17,17,17,12,12,11 + 2, 3, 5, 5, 5, 6, 6, 8, 8,17,17,17,12,12,12 + 2, 3, 4, 4, 6, 6, 8, 8,17,17,17,17,17,12,12 + 2, 3, 4, 4, 6, 6,18,18,18,18,17,17,17,12,12 + */ + // subtract one because the array is written in FORTRAN in mind - in C arrays start @ zero + static const unsigned short select[] = + { + 0, 0 , 1 , 12 ,12 , 12 , 12 , 12 , 12 , 12 , 12 , 15 , 15 , 15 , 8, + 0 , 1 , 1 , 2 , 2 , 4 , 4 , 13 , 13 , 14 , 14 , 15 , 15 , 15 , 8, + 1 , 1 , 2 , 2 , 2 , 4 , 4 , 14 , 14 , 14 , 14 , 15 , 15 , 15 , 9, + 1 , 1 , 2 , 4 , 4 , 4 , 4 , 6 , 6 , 15 , 15 , 15 , 15 , 15 , 9, + 1 , 2 , 2 , 4 , 4 , 5 , 5 , 7 , 7 , 16 ,16 , 16 , 11 , 11 , 10, + 1 , 2 , 4 , 4 , 4 , 5 , 5 , 7 , 7 , 16 , 16 , 16 , 11 , 11 , 11, + 1 , 2 , 3 , 3 , 5 , 5 , 7 , 7 , 16 , 16 , 16 , 16 , 16 , 11 , 11, + 1 , 2 , 3 , 3 , 5 , 5 , 17 , 17 , 17 , 17 , 16 , 16 , 16 , 11 , 11 + }; + + unsigned short ihint = 14, iaint = 7; + for(unsigned short i = 0; i != 14; i++) + { + if( h <= hrange[i] ) + { + ihint = i; + break; + } + } // for(unsigned short i = 0; i != 14; i++) + + for(unsigned short i = 0; i != 7; i++) + { + if( a <= arange[i] ) + { + iaint = i; + break; + } + } // for(unsigned short i = 0; i != 7; i++) + + // interprete select array as 8x15 matrix + return select[iaint*15 + ihint]; + + } // unsigned short owens_t_compute_code(const RealType h, const RealType a) + + template<typename RealType> + inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<53>&) + { + static const unsigned short ord[] = {2, 3, 4, 5, 7, 10, 12, 18, 10, 20, 30, 0, 4, 7, 8, 20, 0, 0}; // 18 entries + + BOOST_ASSERT(icode<18); + + return ord[icode]; + } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<53> const&) + + template<typename RealType> + inline unsigned short owens_t_get_order_imp(const unsigned short icode, RealType, const mpl::int_<64>&) + { + // method ================>>> {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6} + static const unsigned short ord[] = {3, 4, 5, 6, 8, 11, 13, 19, 10, 20, 30, 0, 7, 10, 11, 23, 0, 0}; // 18 entries + + BOOST_ASSERT(icode<18); + + return ord[icode]; + } // unsigned short owens_t_get_order(const unsigned short icode, RealType, mpl::int<64> const&) + + template<typename RealType, typename Policy> + inline unsigned short owens_t_get_order(const unsigned short icode, RealType r, const Policy&) + { + typedef typename policies::precision<RealType, Policy>::type precision_type; + typedef typename mpl::if_< + mpl::or_< + mpl::less_equal<precision_type, mpl::int_<0> >, + mpl::greater<precision_type, mpl::int_<53> > + >, + mpl::int_<64>, + mpl::int_<53> + >::type tag_type; + + return owens_t_get_order_imp(icode, r, tag_type()); + } + + // compute the value of Owen's T function with method T1 from the reference paper + template<typename RealType> + inline RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const RealType hs = -h*h*half<RealType>(); + const RealType dhs = exp( hs ); + const RealType as = a*a; + + unsigned short j=1; + RealType jj = 1; + RealType aj = a * one_div_two_pi<RealType>(); + RealType dj = expm1( hs ); + RealType gj = hs*dhs; + + RealType val = atan( a ) * one_div_two_pi<RealType>(); + + while( true ) + { + val += dj*aj/jj; + + if( m <= j ) + break; + + j++; + jj += static_cast<RealType>(2); + aj *= as; + dj = gj - dj; + gj *= hs / static_cast<RealType>(j); + } // while( true ) + + return val; + } // RealType owens_t_T1(const RealType h, const RealType a, const unsigned short m) + + // compute the value of Owen's T function with method T2 from the reference paper + template<typename RealType, class Policy> + inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::false_&) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const unsigned short maxii = m+m+1; + const RealType hs = h*h; + const RealType as = -a*a; + const RealType y = static_cast<RealType>(1) / hs; + + unsigned short ii = 1; + RealType val = 0; + RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>(); + RealType z = owens_t_znorm1(ah)/h; + + while( true ) + { + val += z; + if( maxii <= ii ) + { + val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>(); + break; + } // if( maxii <= ii ) + z = y * ( vi - static_cast<RealType>(ii) * z ); + vi *= as; + ii += 2; + } // while( true ) + + return val; + } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah) + + // compute the value of Owen's T function with method T3 from the reference paper + template<typename RealType> + inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<53>&) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const unsigned short m = 20; + + static const RealType c2[] = + { + 0.99999999999999987510, + -0.99999999999988796462, 0.99999999998290743652, + -0.99999999896282500134, 0.99999996660459362918, + -0.99999933986272476760, 0.99999125611136965852, + -0.99991777624463387686, 0.99942835555870132569, + -0.99697311720723000295, 0.98751448037275303682, + -0.95915857980572882813, 0.89246305511006708555, + -0.76893425990463999675, 0.58893528468484693250, + -0.38380345160440256652, 0.20317601701045299653, + -0.82813631607004984866E-01, 0.24167984735759576523E-01, + -0.44676566663971825242E-02, 0.39141169402373836468E-03 + }; + + const RealType as = a*a; + const RealType hs = h*h; + const RealType y = static_cast<RealType>(1)/hs; + + RealType ii = 1; + unsigned short i = 0; + RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>(); + RealType zi = owens_t_znorm1(ah)/h; + RealType val = 0; + + while( true ) + { + BOOST_ASSERT(i < 21); + val += zi*c2[i]; + if( m <= i ) // if( m < i+1 ) + { + val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>(); + break; + } // if( m < i ) + zi = y * (ii*zi - vi); + vi *= as; + ii += 2; + i++; + } // while( true ) + + return val; + } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah) + + // compute the value of Owen's T function with method T3 from the reference paper + template<class RealType> + inline RealType owens_t_T3_imp(const RealType h, const RealType a, const RealType ah, const mpl::int_<64>&) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const unsigned short m = 30; + + static const RealType c2[] = + { + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999999999729978162447266851932041876728736094298092917625009873), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99999999999999999999467056379678391810626533251885323416799874878563998732905968), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999999999824849349313270659391127814689133077036298754586814091034842536), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999999997703859616213643405880166422891953033591551179153879839440241685), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99999999999998394883415238173334565554173013941245103172035286759201504179038147), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999999993063616095509371081203145247992197457263066869044528823599399470977), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999999797336340409464429599229870590160411238245275855903767652432017766116267), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999999999574958412069046680119051639753412378037565521359444170241346845522403274), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999999933226234193375324943920160947158239076786103108097456617750134812033362048), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.9999999188923242461073033481053037468263536806742737922476636768006622772762168467), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.9999992195143483674402853783549420883055129680082932629160081128947764415749728967), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.999993935137206712830997921913316971472227199741857386575097250553105958772041501), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.99996135597690552745362392866517133091672395614263398912807169603795088421057688716), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.99979556366513946026406788969630293820987757758641211293079784585126692672425362469), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.999092789629617100153486251423850590051366661947344315423226082520411961968929483), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.996593837411918202119308620432614600338157335862888580671450938858935084316004769854), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.98910017138386127038463510314625339359073956513420458166238478926511821146316469589567), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.970078558040693314521331982203762771512160168582494513347846407314584943870399016019), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.92911438683263187495758525500033707204091967947532160289872782771388170647150321633673), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.8542058695956156057286980736842905011429254735181323743367879525470479126968822863), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.73796526033030091233118357742803709382964420335559408722681794195743240930748630755), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.58523469882837394570128599003785154144164680587615878645171632791404210655891158), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.415997776145676306165661663581868460503874205343014196580122174949645271353372263), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.2588210875241943574388730510317252236407805082485246378222935376279663808416534365), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.1375535825163892648504646951500265585055789019410617565727090346559210218472356689), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.0607952766325955730493900985022020434830339794955745989150270485056436844239206648), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0216337683299871528059836483840390514275488679530797294557060229266785853764115), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -0.00593405693455186729876995814181203900550014220428843483927218267309209471516256), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 0.0011743414818332946510474576182739210553333860106811865963485870668929503649964142), + BOOST_MATH_BIG_CONSTANT(RealType, 260, -1.489155613350368934073453260689881330166342484405529981510694514036264969925132e-4), + BOOST_MATH_BIG_CONSTANT(RealType, 260, 9.072354320794357587710929507988814669454281514268844884841547607134260303118208e-6) + }; + + const RealType as = a*a; + const RealType hs = h*h; + const RealType y = 1 / hs; + + RealType ii = 1; + unsigned short i = 0; + RealType vi = a * exp( -ah*ah*half<RealType>() ) * one_div_root_two_pi<RealType>(); + RealType zi = owens_t_znorm1(ah)/h; + RealType val = 0; + + while( true ) + { + BOOST_ASSERT(i < 31); + val += zi*c2[i]; + if( m <= i ) // if( m < i+1 ) + { + val *= exp( -hs*half<RealType>() ) * one_div_root_two_pi<RealType>(); + break; + } // if( m < i ) + zi = y * (ii*zi - vi); + vi *= as; + ii += 2; + i++; + } // while( true ) + + return val; + } // RealType owens_t_T3(const RealType h, const RealType a, const RealType ah) + + template<class RealType, class Policy> + inline RealType owens_t_T3(const RealType h, const RealType a, const RealType ah, const Policy&) + { + typedef typename policies::precision<RealType, Policy>::type precision_type; + typedef typename mpl::if_< + mpl::or_< + mpl::less_equal<precision_type, mpl::int_<0> >, + mpl::greater<precision_type, mpl::int_<53> > + >, + mpl::int_<64>, + mpl::int_<53> + >::type tag_type; + + return owens_t_T3_imp(h, a, ah, tag_type()); + } + + // compute the value of Owen's T function with method T4 from the reference paper + template<typename RealType> + inline RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const unsigned short maxii = m+m+1; + const RealType hs = h*h; + const RealType as = -a*a; + + unsigned short ii = 1; + RealType ai = a * exp( -hs*(static_cast<RealType>(1)-as)*half<RealType>() ) * one_div_two_pi<RealType>(); + RealType yi = 1; + RealType val = 0; + + while( true ) + { + val += ai*yi; + if( maxii <= ii ) + break; + ii += 2; + yi = (static_cast<RealType>(1)-hs*yi) / static_cast<RealType>(ii); + ai *= as; + } // while( true ) + + return val; + } // RealType owens_t_T4(const RealType h, const RealType a, const unsigned short m) + + // compute the value of Owen's T function with method T5 from the reference paper + template<typename RealType> + inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<53>&) + { + BOOST_MATH_STD_USING + /* + NOTICE: + - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre + polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre + quadrature, because T5(h,a,m) contains only x^2 terms. + - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor + of 1/(2*pi) according to T5(h,a,m). + */ + + const unsigned short m = 13; + static const RealType pts[] = {0.35082039676451715489E-02, + 0.31279042338030753740E-01, 0.85266826283219451090E-01, + 0.16245071730812277011, 0.25851196049125434828, + 0.36807553840697533536, 0.48501092905604697475, + 0.60277514152618576821, 0.71477884217753226516, + 0.81475510988760098605, 0.89711029755948965867, + 0.95723808085944261843, 0.99178832974629703586}; + static const RealType wts[] = { 0.18831438115323502887E-01, + 0.18567086243977649478E-01, 0.18042093461223385584E-01, + 0.17263829606398753364E-01, 0.16243219975989856730E-01, + 0.14994592034116704829E-01, 0.13535474469662088392E-01, + 0.11886351605820165233E-01, 0.10070377242777431897E-01, + 0.81130545742299586629E-02, 0.60419009528470238773E-02, + 0.38862217010742057883E-02, 0.16793031084546090448E-02}; + + const RealType as = a*a; + const RealType hs = -h*h*boost::math::constants::half<RealType>(); + + RealType val = 0; + for(unsigned short i = 0; i < m; ++i) + { + BOOST_ASSERT(i < 13); + const RealType r = static_cast<RealType>(1) + as*pts[i]; + val += wts[i] * exp( hs*r ) / r; + } // for(unsigned short i = 0; i < m; ++i) + + return val*a; + } // RealType owens_t_T5(const RealType h, const RealType a) + + // compute the value of Owen's T function with method T5 from the reference paper + template<typename RealType> + inline RealType owens_t_T5_imp(const RealType h, const RealType a, const mpl::int_<64>&) + { + BOOST_MATH_STD_USING + /* + NOTICE: + - The pts[] array contains the squares (!) of the abscissas, i.e. the roots of the Legendre + polynomial P_n(x), instead of the plain roots as required in Gauss-Legendre + quadrature, because T5(h,a,m) contains only x^2 terms. + - The wts[] array contains the weights for Gauss-Legendre quadrature scaled with a factor + of 1/(2*pi) according to T5(h,a,m). + */ + + const unsigned short m = 19; + static const RealType pts[] = { + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0016634282895983227941), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.014904509242697054183), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.04103478879005817919), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.079359853513391511008), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.1288612130237615133), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.18822336642448518856), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.25586876186122962384), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.32999972011807857222), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.40864620815774761438), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.48971819306044782365), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.57106118513245543894), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.6505134942981533829), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.72596367859928091618), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.79540665919549865924), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.85699701386308739244), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.90909804422384697594), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.95032536436570154409), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.97958418733152273717), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.99610366384229088321) + }; + static const RealType wts[] = { + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012975111395684900835), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012888764187499150078), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012716644398857307844), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012459897461364705691), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.012120231988292330388), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011699908404856841158), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.011201723906897224448), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.010628993848522759853), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0099855296835573320047), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0092756136096132857933), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0085039700881139589055), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0076757344408814561254), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0067964187616556459109), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.005871875456524750363), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0049082589542498110071), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0039119870792519721409), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0028897090921170700834), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.0018483371329504443947), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 0.00079623320100438873578) + }; + + const RealType as = a*a; + const RealType hs = -h*h*boost::math::constants::half<RealType>(); + + RealType val = 0; + for(unsigned short i = 0; i < m; ++i) + { + BOOST_ASSERT(i < 19); + const RealType r = 1 + as*pts[i]; + val += wts[i] * exp( hs*r ) / r; + } // for(unsigned short i = 0; i < m; ++i) + + return val*a; + } // RealType owens_t_T5(const RealType h, const RealType a) + + template<class RealType, class Policy> + inline RealType owens_t_T5(const RealType h, const RealType a, const Policy&) + { + typedef typename policies::precision<RealType, Policy>::type precision_type; + typedef typename mpl::if_< + mpl::or_< + mpl::less_equal<precision_type, mpl::int_<0> >, + mpl::greater<precision_type, mpl::int_<53> > + >, + mpl::int_<64>, + mpl::int_<53> + >::type tag_type; + + return owens_t_T5_imp(h, a, tag_type()); + } + + + // compute the value of Owen's T function with method T6 from the reference paper + template<typename RealType> + inline RealType owens_t_T6(const RealType h, const RealType a) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const RealType normh = owens_t_znorm2( h ); + const RealType y = static_cast<RealType>(1) - a; + const RealType r = atan2(y, static_cast<RealType>(1 + a) ); + + RealType val = normh * ( static_cast<RealType>(1) - normh ) * half<RealType>(); + + if( r != 0 ) + val -= r * exp( -y*h*h*half<RealType>()/r ) * one_div_two_pi<RealType>(); + + return val; + } // RealType owens_t_T6(const RealType h, const RealType a, const unsigned short m) + + template <class T, class Policy> + std::pair<T, T> owens_t_T1_accelerated(T h, T a, const Policy& pol) + { + // + // This is the same series as T1, but: + // * The Taylor series for atan has been combined with that for T1, + // reducing but not eliminating cancellation error. + // * The resulting alternating series is then accelerated using method 1 + // from H. Cohen, F. Rodriguez Villegas, D. Zagier, + // "Convergence acceleration of alternating series", Bonn, (1991). + // + BOOST_MATH_STD_USING + static const char* function = "boost::math::owens_t<%1%>(%1%, %1%)"; + T half_h_h = h * h / 2; + T a_pow = a; + T aa = a * a; + T exp_term = exp(-h * h / 2); + T one_minus_dj_sum = exp_term; + T sum = a_pow * exp_term; + T dj_pow = exp_term; + T term = sum; + T abs_err; + int j = 1; + + // + // Normally with this form of series acceleration we can calculate + // up front how many terms will be required - based on the assumption + // that each term decreases in size by a factor of 3. However, + // that assumption does not apply here, as the underlying T1 series can + // go quite strongly divergent in the early terms, before strongly + // converging later. Various "guestimates" have been tried to take account + // of this, but they don't always work.... so instead set "n" to the + // largest value that won't cause overflow later, and abort iteration + // when the last accelerated term was small enough... + // + int n; + try + { + n = itrunc(T(tools::log_max_value<T>() / 6)); + } + catch(...) + { + n = (std::numeric_limits<int>::max)(); + } + n = (std::min)(n, 1500); + T d = pow(3 + sqrt(T(8)), n); + d = (d + 1 / d) / 2; + T b = -1; + T c = -d; + c = b - c; + sum *= c; + b = -n * n * b * 2; + abs_err = ldexp(fabs(sum), -tools::digits<T>()); + + while(j < n) + { + a_pow *= aa; + dj_pow *= half_h_h / j; + one_minus_dj_sum += dj_pow; + term = one_minus_dj_sum * a_pow / (2 * j + 1); + c = b - c; + sum += c * term; + abs_err += ldexp(std::max(T(fabs(sum)), T(fabs(c*term))), -tools::digits<T>()); + b = (j + n) * (j - n) * b / ((j + T(0.5)) * (j + 1)); + ++j; + // + // Include an escape route to prevent calculating too many terms: + // + if((j > 10) && (fabs(sum * tools::epsilon<T>()) > fabs(c * term))) + break; + } + abs_err += fabs(c * term); + if(sum < 0) // sum must always be positive, if it's negative something really bad has happend: + policies::raise_evaluation_error(function, 0, T(0), pol); + return std::pair<T, T>((sum / d) / boost::math::constants::two_pi<T>(), abs_err / sum); + } + + template<typename RealType, class Policy> + inline RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah, const Policy&, const mpl::true_&) + { + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const unsigned short maxii = m+m+1; + const RealType hs = h*h; + const RealType as = -a*a; + const RealType y = static_cast<RealType>(1) / hs; + + unsigned short ii = 1; + RealType val = 0; + RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>(); + RealType z = owens_t_znorm1(ah)/h; + RealType last_z = fabs(z); + RealType lim = policies::get_epsilon<RealType, Policy>(); + + while( true ) + { + val += z; + // + // This series stops converging after a while, so put a limit + // on how far we go before returning our best guess: + // + if((fabs(lim * val) > fabs(z)) || ((ii > maxii) && (fabs(z) > last_z)) || (z == 0)) + { + val *= exp( -hs*half<RealType>() ) / root_two_pi<RealType>(); + break; + } // if( maxii <= ii ) + last_z = fabs(z); + z = y * ( vi - static_cast<RealType>(ii) * z ); + vi *= as; + ii += 2; + } // while( true ) + + return val; + } // RealType owens_t_T2(const RealType h, const RealType a, const unsigned short m, const RealType ah) + + template<typename RealType, class Policy> + inline std::pair<RealType, RealType> owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&) + { + // + // This is the same series as T2, but with acceleration applied. + // Note that we have to be *very* careful to check that nothing bad + // has happened during evaluation - this series will go divergent + // and/or fail to alternate at a drop of a hat! :-( + // + BOOST_MATH_STD_USING + using namespace boost::math::constants; + + const RealType hs = h*h; + const RealType as = -a*a; + const RealType y = static_cast<RealType>(1) / hs; + + unsigned short ii = 1; + RealType val = 0; + RealType vi = a * exp( -ah*ah*half<RealType>() ) / root_two_pi<RealType>(); + RealType z = boost::math::detail::owens_t_znorm1(ah)/h; + RealType last_z = fabs(z); + + // + // Normally with this form of series acceleration we can calculate + // up front how many terms will be required - based on the assumption + // that each term decreases in size by a factor of 3. However, + // that assumption does not apply here, as the underlying T1 series can + // go quite strongly divergent in the early terms, before strongly + // converging later. Various "guestimates" have been tried to take account + // of this, but they don't always work.... so instead set "n" to the + // largest value that won't cause overflow later, and abort iteration + // when the last accelerated term was small enough... + // + int n; + try + { + n = itrunc(RealType(tools::log_max_value<RealType>() / 6)); + } + catch(...) + { + n = (std::numeric_limits<int>::max)(); + } + n = (std::min)(n, 1500); + RealType d = pow(3 + sqrt(RealType(8)), n); + d = (d + 1 / d) / 2; + RealType b = -1; + RealType c = -d; + int s = 1; + + for(int k = 0; k < n; ++k) + { + // + // Check for both convergence and whether the series has gone bad: + // + if( + (fabs(z) > last_z) // Series has gone divergent, abort + || (fabs(val) * tools::epsilon<RealType>() > fabs(c * s * z)) // Convergence! + || (z * s < 0) // Series has stopped alternating - all bets are off - abort. + ) + { + break; + } + c = b - c; + val += c * s * z; + b = (k + n) * (k - n) * b / ((k + RealType(0.5)) * (k + 1)); + last_z = fabs(z); + s = -s; + z = y * ( vi - static_cast<RealType>(ii) * z ); + vi *= as; + ii += 2; + } // while( true ) + RealType err = fabs(c * z) / val; + return std::pair<RealType, RealType>(val * exp( -hs*half<RealType>() ) / (d * root_two_pi<RealType>()), err); + } // RealType owens_t_T2_accelerated(const RealType h, const RealType a, const RealType ah, const Policy&) + + template<typename RealType, typename Policy> + inline RealType T4_mp(const RealType h, const RealType a, const Policy& pol) + { + BOOST_MATH_STD_USING + + const RealType hs = h*h; + const RealType as = -a*a; + + unsigned short ii = 1; + RealType ai = constants::one_div_two_pi<RealType>() * a * exp( -0.5*hs*(1.0-as) ); + RealType yi = 1.0; + RealType val = 0.0; + + RealType lim = boost::math::policies::get_epsilon<RealType, Policy>(); + + while( true ) + { + RealType term = ai*yi; + val += term; + if((yi != 0) && (fabs(val * lim) > fabs(term))) + break; + ii += 2; + yi = (1.0-hs*yi) / static_cast<RealType>(ii); + ai *= as; + if(ii > (std::min)(1500, (int)policies::get_max_series_iterations<Policy>())) + policies::raise_evaluation_error("boost::math::owens_t<%1%>", 0, val, pol); + } // while( true ) + + return val; + } // arg_type owens_t_T4(const arg_type h, const arg_type a, const unsigned short m) + + + // This routine dispatches the call to one of six subroutines, depending on the values + // of h and a. + // preconditions: h >= 0, 0<=a<=1, ah=a*h + // + // Note there are different versions for different precisions.... + template<typename RealType, typename Policy> + inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, mpl::int_<64> const&) + { + // Simple main case for 64-bit precision or less, this is as per the Patefield-Tandy paper: + BOOST_MATH_STD_USING + // + // Handle some special cases first, these are from + // page 1077 of Owen's original paper: + // + if(h == 0) + { + return atan(a) * constants::one_div_two_pi<RealType>(); + } + if(a == 0) + { + return 0; + } + if(a == 1) + { + return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2; + } + if(a >= tools::max_value<RealType>()) + { + return owens_t_znorm2(RealType(fabs(h))); + } + RealType val = 0; // avoid compiler warnings, 0 will be overwritten in any case + const unsigned short icode = owens_t_compute_code(h, a); + const unsigned short m = owens_t_get_order(icode, val /* just a dummy for the type */, pol); + static const unsigned short meth[] = {1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 4, 4, 4, 4, 5, 6}; // 18 entries + + // determine the appropriate method, T1 ... T6 + switch( meth[icode] ) + { + case 1: // T1 + val = owens_t_T1(h,a,m); + break; + case 2: // T2 + typedef typename policies::precision<RealType, Policy>::type precision_type; + typedef mpl::bool_<(precision_type::value == 0) || (precision_type::value > 64)> tag_type; + val = owens_t_T2(h, a, m, ah, pol, tag_type()); + break; + case 3: // T3 + val = owens_t_T3(h,a,ah, pol); + break; + case 4: // T4 + val = owens_t_T4(h,a,m); + break; + case 5: // T5 + val = owens_t_T5(h,a, pol); + break; + case 6: // T6 + val = owens_t_T6(h,a); + break; + default: + BOOST_THROW_EXCEPTION(std::logic_error("selection routine in Owen's T function failed")); + } + return val; + } + + template<typename RealType, typename Policy> + inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<65>&) + { + // Arbitrary precision version: + BOOST_MATH_STD_USING + // + // Handle some special cases first, these are from + // page 1077 of Owen's original paper: + // + if(h == 0) + { + return atan(a) * constants::one_div_two_pi<RealType>(); + } + if(a == 0) + { + return 0; + } + if(a == 1) + { + return owens_t_znorm2(RealType(-h)) * owens_t_znorm2(h) / 2; + } + if(a >= tools::max_value<RealType>()) + { + return owens_t_znorm2(RealType(fabs(h))); + } + // Attempt arbitrary precision code, this will throw if it goes wrong: + typedef typename boost::math::policies::normalise<Policy, boost::math::policies::evaluation_error<> >::type forwarding_policy; + std::pair<RealType, RealType> p1(0, tools::max_value<RealType>()), p2(0, tools::max_value<RealType>()); + RealType target_precision = policies::get_epsilon<RealType, Policy>() * 1000; + bool have_t1(false), have_t2(false); + if(ah < 3) + { + try + { + have_t1 = true; + p1 = owens_t_T1_accelerated(h, a, forwarding_policy()); + if(p1.second < target_precision) + return p1.first; + } + catch(const boost::math::evaluation_error&){} // T1 may fail and throw, that's OK + } + if(ah > 1) + { + try + { + have_t2 = true; + p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy()); + if(p2.second < target_precision) + return p2.first; + } + catch(const boost::math::evaluation_error&){} // T2 may fail and throw, that's OK + } + // + // If we haven't tried T1 yet, do it now - sometimes it succeeds and the number of iterations + // is fairly low compared to T4. + // + if(!have_t1) + { + try + { + have_t1 = true; + p1 = owens_t_T1_accelerated(h, a, forwarding_policy()); + if(p1.second < target_precision) + return p1.first; + } + catch(const boost::math::evaluation_error&){} // T1 may fail and throw, that's OK + } + // + // If we haven't tried T2 yet, do it now - sometimes it succeeds and the number of iterations + // is fairly low compared to T4. + // + if(!have_t2) + { + try + { + have_t2 = true; + p2 = owens_t_T2_accelerated(h, a, ah, forwarding_policy()); + if(p2.second < target_precision) + return p2.first; + } + catch(const boost::math::evaluation_error&){} // T2 may fail and throw, that's OK + } + // + // OK, nothing left to do but try the most expensive option which is T4, + // this is often slow to converge, but when it does converge it tends to + // be accurate: + try + { + return T4_mp(h, a, pol); + } + catch(const boost::math::evaluation_error&){} // T4 may fail and throw, that's OK + // + // Now look back at the results from T1 and T2 and see if either gave better + // results than we could get from the 64-bit precision versions. + // + if((std::min)(p1.second, p2.second) < 1e-20) + { + return p1.second < p2.second ? p1.first : p2.first; + } + // + // We give up - no arbitrary precision versions succeeded! + // + return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>()); + } // RealType owens_t_dispatch(RealType h, RealType a, RealType ah) + template<typename RealType, typename Policy> + inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol, const mpl::int_<0>&) + { + // We don't know what the precision is until runtime: + if(tools::digits<RealType>() <= 64) + return owens_t_dispatch(h, a, ah, pol, mpl::int_<64>()); + return owens_t_dispatch(h, a, ah, pol, mpl::int_<65>()); + } + template<typename RealType, typename Policy> + inline RealType owens_t_dispatch(const RealType h, const RealType a, const RealType ah, const Policy& pol) + { + // Figure out the precision and forward to the correct version: + typedef typename policies::precision<RealType, Policy>::type precision_type; + typedef typename mpl::if_c< + precision_type::value == 0, + mpl::int_<0>, + typename mpl::if_c< + precision_type::value <= 64, + mpl::int_<64>, + mpl::int_<65> + >::type + >::type tag_type; + return owens_t_dispatch(h, a, ah, pol, tag_type()); + } + // compute Owen's T function, T(h,a), for arbitrary values of h and a + template<typename RealType, class Policy> + inline RealType owens_t(RealType h, RealType a, const Policy& pol) + { + BOOST_MATH_STD_USING + // exploit that T(-h,a) == T(h,a) + h = fabs(h); + + // Use equation (2) in the paper to remap the arguments + // such that h>=0 and 0<=a<=1 for the call of the actual + // computation routine. + + const RealType fabs_a = fabs(a); + const RealType fabs_ah = fabs_a*h; + + RealType val = 0.0; // avoid compiler warnings, 0.0 will be overwritten in any case + + if(fabs_a <= 1) + { + val = owens_t_dispatch(h, fabs_a, fabs_ah, pol); + } // if(fabs_a <= 1.0) + else + { + if( h <= 0.67 ) + { + const RealType normh = owens_t_znorm1(h); + const RealType normah = owens_t_znorm1(fabs_ah); + val = static_cast<RealType>(1)/static_cast<RealType>(4) - normh*normah - + owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol); + } // if( h <= 0.67 ) + else + { + const RealType normh = detail::owens_t_znorm2(h); + const RealType normah = detail::owens_t_znorm2(fabs_ah); + val = constants::half<RealType>()*(normh+normah) - normh*normah - + owens_t_dispatch(fabs_ah, static_cast<RealType>(1 / fabs_a), h, pol); + } // else [if( h <= 0.67 )] + } // else [if(fabs_a <= 1)] + + // exploit that T(h,-a) == -T(h,a) + if(a < 0) + { + return -val; + } // if(a < 0) + + return val; + } // RealType owens_t(RealType h, RealType a) + + template <class T, class Policy, class tag> + struct owens_t_initializer + { + struct init + { + init() + { + do_init(tag()); + } + template <int N> + static void do_init(const mpl::int_<N>&){} + static void do_init(const mpl::int_<64>&) + { + boost::math::owens_t(static_cast<T>(7), static_cast<T>(0.96875), Policy()); + boost::math::owens_t(static_cast<T>(2), static_cast<T>(0.5), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } + }; + + template <class T, class Policy, class tag> + const typename owens_t_initializer<T, Policy, tag>::init owens_t_initializer<T, Policy, tag>::initializer; + + } // namespace detail + + template <class T1, class T2, class Policy> + inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a, const Policy& pol) + { + typedef typename tools::promote_args<T1, T2>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + typedef typename policies::precision<value_type, Policy>::type precision_type; + typedef typename mpl::if_c< + precision_type::value == 0, + mpl::int_<0>, + typename mpl::if_c< + precision_type::value <= 64, + mpl::int_<64>, + mpl::int_<65> + >::type + >::type tag_type; + + detail::owens_t_initializer<result_type, Policy, tag_type>::force_instantiate(); + + return policies::checked_narrowing_cast<result_type, Policy>(detail::owens_t(static_cast<value_type>(h), static_cast<value_type>(a), pol), "boost::math::owens_t<%1%>(%1%,%1%)"); + } + + template <class T1, class T2> + inline typename tools::promote_args<T1, T2>::type owens_t(T1 h, T2 a) + { + return owens_t(h, a, policies::policy<>()); + } + + + } // namespace math +} // namespace boost + +#endif +// EOF diff --git a/boost/math/special_functions/round.hpp b/boost/math/special_functions/round.hpp index f11c6aeb1d..2b4497e198 100644 --- a/boost/math/special_functions/round.hpp +++ b/boost/math/special_functions/round.hpp @@ -76,7 +76,7 @@ inline boost::long_long_type llround(const T& v, const Policy& pol) BOOST_MATH_STD_USING T r = boost::math::round(v, pol); if((r > (std::numeric_limits<boost::long_long_type>::max)()) || (r < (std::numeric_limits<boost::long_long_type>::min)())) - return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", 0, v, 0LL, pol)); + return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::llround<%1%>(%1%)", 0, v, static_cast<boost::long_long_type>(0), pol)); return static_cast<boost::long_long_type>(r); } template <class T> diff --git a/boost/math/special_functions/trunc.hpp b/boost/math/special_functions/trunc.hpp index 520ae89f5d..7346afe6d1 100644 --- a/boost/math/special_functions/trunc.hpp +++ b/boost/math/special_functions/trunc.hpp @@ -76,7 +76,7 @@ inline boost::long_long_type lltrunc(const T& v, const Policy& pol) BOOST_MATH_STD_USING T r = boost::math::trunc(v, pol); if((r > (std::numeric_limits<boost::long_long_type>::max)()) || (r < (std::numeric_limits<boost::long_long_type>::min)())) - return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::lltrunc<%1%>(%1%)", 0, v, 0LL, pol)); + return static_cast<boost::long_long_type>(policies::raise_rounding_error("boost::math::lltrunc<%1%>(%1%)", 0, v, static_cast<boost::long_long_type>(0), pol)); return static_cast<boost::long_long_type>(r); } template <class T> diff --git a/boost/math/special_functions/zeta.hpp b/boost/math/special_functions/zeta.hpp index 4ed5f6a705..011182718e 100644 --- a/boost/math/special_functions/zeta.hpp +++ b/boost/math/special_functions/zeta.hpp @@ -542,7 +542,7 @@ T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<64>&) } template <class T, class Policy> -T zeta_imp_prec(T s, T sc, const Policy& pol, const mpl::int_<113>&) +T zeta_imp_prec(T s, T sc, const Policy&, const mpl::int_<113>&) { BOOST_MATH_STD_USING T result; @@ -896,6 +896,49 @@ T zeta_imp(T s, T sc, const Policy& pol, const Tag& tag) return result; } +template <class T, class Policy, class tag> +struct zeta_initializer +{ + struct init + { + init() + { + do_init(tag()); + } + static void do_init(const mpl::int_<0>&){} + static void do_init(const mpl::int_<53>&){} + static void do_init(const mpl::int_<64>&) + { + boost::math::zeta(static_cast<T>(0.5), Policy()); + boost::math::zeta(static_cast<T>(1.5), Policy()); + boost::math::zeta(static_cast<T>(3.5), Policy()); + boost::math::zeta(static_cast<T>(6.5), Policy()); + boost::math::zeta(static_cast<T>(14.5), Policy()); + boost::math::zeta(static_cast<T>(40.5), Policy()); + } + static void do_init(const mpl::int_<113>&) + { + boost::math::zeta(static_cast<T>(0.5), Policy()); + boost::math::zeta(static_cast<T>(1.5), Policy()); + boost::math::zeta(static_cast<T>(3.5), Policy()); + boost::math::zeta(static_cast<T>(5.5), Policy()); + boost::math::zeta(static_cast<T>(9.5), Policy()); + boost::math::zeta(static_cast<T>(16.5), Policy()); + boost::math::zeta(static_cast<T>(25), Policy()); + boost::math::zeta(static_cast<T>(70), Policy()); + } + void force_instantiate()const{} + }; + static const init initializer; + static void force_instantiate() + { + initializer.force_instantiate(); + } +}; + +template <class T, class Policy, class tag> +const typename zeta_initializer<T, Policy, tag>::init zeta_initializer<T, Policy, tag>::initializer; + } // detail template <class T, class Policy> @@ -929,6 +972,8 @@ inline typename tools::promote_args<T>::type zeta(T s, const Policy&) >::type tag_type; //typedef mpl::int_<0> tag_type; + detail::zeta_initializer<value_type, forwarding_policy, tag_type>::force_instantiate(); + return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::zeta_imp( static_cast<value_type>(s), static_cast<value_type>(1 - static_cast<value_type>(s)), diff --git a/boost/math/tools/big_constant.hpp b/boost/math/tools/big_constant.hpp index 45f13ceef9..119063164a 100644 --- a/boost/math/tools/big_constant.hpp +++ b/boost/math/tools/big_constant.hpp @@ -11,15 +11,17 @@ #include <boost/lexical_cast.hpp> #include <boost/type_traits/is_convertible.hpp> -namespace boost{ namespace math{ namespace tools{ +namespace boost{ namespace math{ + +namespace tools{ template <class T> -inline T make_big_value(long double v, const char*, mpl::true_ const&, mpl::false_ const&) +inline BOOST_CONSTEXPR_OR_CONST T make_big_value(long double v, const char*, mpl::true_ const&, mpl::false_ const&) { return static_cast<T>(v); } template <class T> -inline T make_big_value(long double v, const char*, mpl::true_ const&, mpl::true_ const&) +inline BOOST_CONSTEXPR_OR_CONST T make_big_value(long double v, const char*, mpl::true_ const&, mpl::true_ const&) { return static_cast<T>(v); } @@ -29,7 +31,7 @@ inline T make_big_value(long double, const char* s, mpl::false_ const&, mpl::fal return boost::lexical_cast<T>(s); } template <class T> -inline const char* make_big_value(long double, const char* s, mpl::false_ const&, mpl::true_ const&) +inline BOOST_CONSTEXPR const char* make_big_value(long double, const char* s, mpl::false_ const&, mpl::true_ const&) { return s; } @@ -38,7 +40,15 @@ inline const char* make_big_value(long double, const char* s, mpl::false_ const& // For constants which might fit in a long double (if it's big enough): // #define BOOST_MATH_BIG_CONSTANT(T, D, x)\ - boost::math::tools::make_big_value<T>(BOOST_JOIN(x, L), BOOST_STRINGIZE(x), mpl::bool_<D <= std::numeric_limits<long double>::digits>(), boost::is_convertible<const char*, T>()) + boost::math::tools::make_big_value<T>(\ + BOOST_JOIN(x, L), \ + BOOST_STRINGIZE(x), \ + mpl::bool_< (is_convertible<long double, T>::value) && \ + ((D <= std::numeric_limits<long double>::digits) \ + || is_floating_point<T>::value \ + || (std::numeric_limits<T>::is_specialized && \ + (std::numeric_limits<T>::digits10 <= std::numeric_limits<long double>::digits10))) >(), \ + boost::is_convertible<const char*, T>()) // // For constants too huge for any conceivable long double (and which generate compiler errors if we try and declare them as such): // diff --git a/boost/math/tools/config.hpp b/boost/math/tools/config.hpp index 96f5d81a43..b1fcd13856 100644 --- a/boost/math/tools/config.hpp +++ b/boost/math/tools/config.hpp @@ -16,6 +16,7 @@ #include <algorithm> // for min and max #include <boost/config/no_tr1/cmath.hpp> #include <climits> +#include <cfloat> #if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) # include <math.h> #endif @@ -24,7 +25,8 @@ #include <boost/math/special_functions/detail/round_fwd.hpp> #if (defined(__CYGWIN__) || defined(__FreeBSD__) || defined(__NetBSD__) \ - || (defined(__hppa) && !defined(__OpenBSD__)) || defined(__NO_LONG_DOUBLE_MATH)) && !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) + || (defined(__hppa) && !defined(__OpenBSD__)) || (defined(__NO_LONG_DOUBLE_MATH) && (DBL_MANT_DIG != LDBL_MANT_DIG))) \ + && !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) # define BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS #endif #if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) @@ -38,6 +40,13 @@ # define BOOST_MATH_CONTROL_FP _control87(MCW_EM,MCW_EM) # include <float.h> #endif +#ifdef __IBMCPP__ +// +// For reasons I don't unserstand, the tests with IMB's compiler all +// pass at long double precision, but fail with real_concept, those tests +// are disabled for now. (JM 2012). +# define BOOST_MATH_NO_REAL_CONCEPT_TESTS +#endif #if (defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__)) && ((LDBL_MANT_DIG == 106) || (__LDBL_MANT_DIG__ == 106)) && !defined(BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS) // // Darwin's rather strange "double double" is rather hard to diff --git a/boost/math/tools/precision.hpp b/boost/math/tools/precision.hpp index 73c969a862..8cdcd4eb87 100644 --- a/boost/math/tools/precision.hpp +++ b/boost/math/tools/precision.hpp @@ -195,7 +195,7 @@ inline T log_max_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) (std::numeric_limits<T>::max_exponent == 128 || std::numeric_limits<T>::max_exponent == 1024 || std::numeric_limits<T>::max_exponent == 16384), - mpl::int_<std::numeric_limits<T>::max_exponent>, + mpl::int_<(std::numeric_limits<T>::max_exponent > INT_MAX ? INT_MAX : static_cast<int>(std::numeric_limits<T>::max_exponent))>, mpl::int_<0> >::type tag_type; BOOST_STATIC_ASSERT( ::std::numeric_limits<T>::is_specialized); @@ -217,7 +217,7 @@ inline T log_min_value(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE(T)) (std::numeric_limits<T>::max_exponent == 128 || std::numeric_limits<T>::max_exponent == 1024 || std::numeric_limits<T>::max_exponent == 16384), - mpl::int_<std::numeric_limits<T>::max_exponent>, + mpl::int_<(std::numeric_limits<T>::max_exponent > INT_MAX ? INT_MAX : static_cast<int>(std::numeric_limits<T>::max_exponent))>, mpl::int_<0> >::type tag_type; diff --git a/boost/math/tools/test.hpp b/boost/math/tools/test.hpp index 90c46e6997..af51d9e3e1 100644 --- a/boost/math/tools/test.hpp +++ b/boost/math/tools/test.hpp @@ -255,6 +255,74 @@ test_result<typename calculate_result_type<A>::value_type> test(const A& a, F1 t return result; } +template <class Real, class A, class F1, class F2> +test_result<Real> test_hetero(const A& a, F1 test_func, F2 expect_func) +{ + typedef typename A::value_type row_type; + typedef Real value_type; + + test_result<value_type> result; + + for(unsigned i = 0; i < a.size(); ++i) + { + const row_type& row = a[i]; + value_type point; + try + { + point = test_func(row); + } + catch(const std::underflow_error&) + { + point = 0; + } + catch(const std::overflow_error&) + { + point = std::numeric_limits<value_type>::has_infinity ? + std::numeric_limits<value_type>::infinity() + : tools::max_value<value_type>(); + } + catch(const std::exception& e) + { + std::cerr << e.what() << std::endl; + print_row(row); + BOOST_ERROR("Unexpected exception."); + // so we don't get further errors: + point = expect_func(row); + } + value_type expected = expect_func(row); + value_type err = relative_error(point, expected); +#ifdef BOOST_INSTRUMENT + if(err != 0) + { + std::cout << row[0] << " " << err; + if(std::numeric_limits<value_type>::is_specialized) + { + std::cout << " (" << err / std::numeric_limits<value_type>::epsilon() << "eps)"; + } + std::cout << std::endl; + } +#endif + if(!(boost::math::isfinite)(point) && (boost::math::isfinite)(expected)) + { + std::cout << "CAUTION: Found non-finite result, when a finite value was expected at entry " << i << "\n"; + std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; + print_row(row); + BOOST_ERROR("Unexpected non-finite result"); + } + if(err > 0.5) + { + std::cout << "CAUTION: Gross error found at entry " << i << ".\n"; + std::cout << "Found: " << point << " Expected " << expected << " Error: " << err << std::endl; + print_row(row); + BOOST_ERROR("Gross error"); + } + result.add(err); + if((result.max)() == err) + result.set_worst(i); + } + return result; +} + } // namespace tools } // namespace math } // namespace boost diff --git a/boost/math/tools/test_data.hpp b/boost/math/tools/test_data.hpp index 4c28d7569d..56bae051b4 100644 --- a/boost/math/tools/test_data.hpp +++ b/boost/math/tools/test_data.hpp @@ -726,7 +726,7 @@ std::ostream& write_code(std::ostream& os, if(a == b) return os; - os << "#define SC_(x) static_cast<T>(BOOST_JOIN(x, L))\n" + os << "#ifndef SC_\n# define SC_(x) static_cast<T>(BOOST_JOIN(x, L))\n#endif\n" " static const boost::array<boost::array<T, " << a->size() << ">, " << data.size() << "> " << name << " = {{\n"; @@ -749,7 +749,7 @@ std::ostream& write_code(std::ostream& os, os << " }"; ++a; } - os << "\n }};\n#undef SC_\n\n"; + os << "\n }};\n//#undef SC_\n\n"; return os; } diff --git a/boost/math/tools/tuple.hpp b/boost/math/tools/tuple.hpp index f61f3e0594..2f297360a2 100644 --- a/boost/math/tools/tuple.hpp +++ b/boost/math/tools/tuple.hpp @@ -9,7 +9,7 @@ #include <boost/tr1/detail/config.hpp> // for BOOST_HAS_TR1_TUPLE -#ifndef BOOST_NO_0X_HDR_TUPLE +#ifndef BOOST_NO_CXX11_HDR_TUPLE #include <tuple> diff --git a/boost/math/tr1.hpp b/boost/math/tr1.hpp index 6a3d078d8a..df8ab0ef47 100644 --- a/boost/math/tr1.hpp +++ b/boost/math/tr1.hpp @@ -104,6 +104,7 @@ namespace boost{ namespace math{ namespace tr1{ extern "C"{ # include <boost/config/auto_link.hpp> #endif +#if !(defined(BOOST_INTEL) && defined(__APPLE__)) && !(defined(__FLT_EVAL_METHOD__) && !defined(__cplusplus)) #ifndef FLT_EVAL_METHOD typedef float float_t; typedef double double_t; @@ -117,6 +118,7 @@ typedef double double_t; typedef long double float_t; typedef long double double_t; #endif +#endif // C99 Functions: double BOOST_MATH_TR1_DECL boost_acosh BOOST_PREVENT_MACRO_SUBSTITUTION(double x) BOOST_MATH_C99_THROW_SPEC; |