From 82548699a56ef3e6108f650e76d68a27e3802bae Mon Sep 17 00:00:00 2001 From: John Maddock Date: Thu, 26 Oct 2006 11:51:59 +0000 Subject: [PATCH] Improved and generally tidied up sqrt1pm1 and powm1. Renamed sqrtp1m1 as sqrt1pm1 to match log1p etc. Added docs for misc power functions, plus weibull distribution. Updated tests to match the above changes. [SVN r3321] --- doc/dist_reference.qbk | 1 + doc/distributions/weibull.qbk | 90 +++++++++++ doc/math.qbk | 6 + doc/powers.qbk | 149 ++++++++++++++++++ .../boost/math/special_functions/gamma.hpp | 6 +- .../boost/math/special_functions/math_fwd.hpp | 2 +- .../boost/math/special_functions/powm1.hpp | 76 ++------- .../boost/math/special_functions/sqrt1pm1.hpp | 34 ++++ .../boost/math/special_functions/sqrtp1m1.hpp | 58 ------- test/powm1_sqrtp1m1_test.cpp | 8 +- test/test_gamma.cpp | 36 ++--- 11 files changed, 316 insertions(+), 150 deletions(-) create mode 100644 doc/distributions/weibull.qbk create mode 100644 doc/powers.qbk create mode 100644 include/boost/math/special_functions/sqrt1pm1.hpp delete mode 100644 include/boost/math/special_functions/sqrtp1m1.hpp diff --git a/doc/dist_reference.qbk b/doc/dist_reference.qbk index 66f86283b..e29579c10 100644 --- a/doc/dist_reference.qbk +++ b/doc/dist_reference.qbk @@ -40,6 +40,7 @@ [include distributions/normal.qbk] [include distributions/poisson.qbk] [include distributions/students_t.qbk] +[include distributions/weibull.qbk] [endsect][/section:dists Distributions] diff --git a/doc/distributions/weibull.qbk b/doc/distributions/weibull.qbk new file mode 100644 index 000000000..42d060816 --- /dev/null +++ b/doc/distributions/weibull.qbk @@ -0,0 +1,90 @@ +[section:weibull Weibull Distribution] + +[template alpha '''α'''] +[template beta '''β'''] +[template GAMMA '''Γ'''] + + +``#include `` + + namespace boost{ namespace math{ + + template + class weibull_distribution; + + typedef weibull_distribution weibull; + + template + class weibull_distribution + { + public: + typedef RealType value_type; + // Construct: + weibull_distribution(RealType shape, RealType scale = 1) + // Accessors: + RealType shape()const; + RealType scale()const; + }; + + }} // namespaces + +The Weibull distribution is a continuous distribution with the probablity +density function: + +f(x; [alpha], [beta]) = ([alpha]\/[beta]) * (x \/ [beta])[super [alpha] - 1] * e[super -(x\/[beta])[super [alpha]]] + +For shape parameter [alpha], and scale parameter [beta]. + +The Weibull distribution is often used in the field of failure analysis, in +particular it can mimic distributions where the failure rate varies over time. + +[h4 Member Functions] + + weibull_distribution(RealType shape, RealType scale = 1); + +Constructs a weibull distribution with shape /shape/ and +scale /scale/. + + RealType shape()const; + +Returns the /shape/ parameter of this distribution. + + RealType scale()const; + +Returns the /scale/ parameter of this distribution. + +[h4 Non-member Accessors] + +All the [link math_toolkit.dist.dist_ref.nmp usual non-member accessor functions] that are generic to all +distributions are supported: __usual_accessors. + +[h4 Accuracy] + +The weibull distribution is implemented in terms of the +standard library log and exp functions +and as such should have very low error rates. + +[h4 Implementation] + + +In the following table [alpha] is the shape parameter of the distribution, +[beta] is it's scale parameter, /x/ is the random variate, /p/ is the probability +and /q = 1-p/. + +[table +[[Function][Implementation Notes]] +[[pdf][Using the relation: pdf = [alpha][beta][super -[alpha] ]x[super [alpha] - 1] e[super -(x/beta)[super alpha]] ]] +[[cdf][Using the relation: p = -__expm1(-(x\/[beta])[super [alpha]]) ]] +[[cdf complement][Using the relation: q = e[super -(x\/[beta])[super [alpha]]] ]] +[[quantile][Using the relation: x = [beta] * (-__log1p(-p))[super 1\/[alpha]] ]] +[[quantile from the complement][Using the relation: x = [beta] * (-log(q))[super 1\/[alpha]] ]] +[[mean][[beta] * [GAMMA](1 + 1\/[alpha]) ]] +[[variance][[beta][super 2]([GAMMA](1 + 2\/[alpha]) - [GAMMA][super 2](1 + 1\/[alpha])) ]] +[[mode][[beta](([alpha] - 1) \/ [alpha])[super 1\/[alpha]] ]] +[[skewness][Refer to [@http://mathworld.wolfram.com/WeibullDistribution.html Weisstein, Eric W. "Weibull Distribution." From MathWorld--A Wolfram Web Resource.] ]] +[[kurtosis][Refer to [@http://mathworld.wolfram.com/WeibullDistribution.html Weisstein, Eric W. "Weibull Distribution." From MathWorld--A Wolfram Web Resource.] ]] +[[kurtosis excess][Refer to [@http://mathworld.wolfram.com/WeibullDistribution.html Weisstein, Eric W. "Weibull Distribution." From MathWorld--A Wolfram Web Resource.] ]] +] + +[endsect][/section:normal_dist Normal] + diff --git a/doc/math.qbk b/doc/math.qbk index 5e4a5ee68..6fe018032 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -101,6 +101,11 @@ and use the function's name as the link text] [def __isnan [link math_toolkit.special.fpclass isnan]] [def __isinf [link math_toolkit.special.fpclass isinf]] [def __isnormal [link math_toolkit.special.fpclass isnormal]] +[def __expm1 [link math_toolkit.special.powers expm1]] +[def __log1p [link math_toolkit.special.powers log1p]] +[def __cbrt [link math_toolkit.special.powers cbrt]] +[def __sqrt1pm1 [link math_toolkit.special.powers sqrt1pm1]] +[def __powm1 [link math_toolkit.special.powers powm1]] @@ -181,6 +186,7 @@ some typical applications in statistics. [include ibeta_inv.qbk] [include beta_gamma_derivatives.qbk] [include fpclassify.qbk] +[include powers.qbk] [include error_handling.qbk] [endsect][/section:special Special Functions] diff --git a/doc/powers.qbk b/doc/powers.qbk new file mode 100644 index 000000000..5d373ac9f --- /dev/null +++ b/doc/powers.qbk @@ -0,0 +1,149 @@ +[section:powers Logs, Powers, Roots and Exponentials] + +[caution __caution ] + +[h4 Synopsis] + +`` +#include +`` + + namespace boost{ namespace math{ + + template + T log1p(T x); + + }} // namespaces + +`` +#include +`` + + namespace boost{ namespace math{ + + template + T expm1(T x); + + }} // namespaces + +`` +#include +`` + + namespace boost{ namespace math{ + + template + T cbrt(T x); + + }} // namespaces + +`` +#include +`` + + namespace boost{ namespace math{ + + template + T sqrt1pm1(T x); + + }} // namespaces + +`` +#include +`` + + namespace boost{ namespace math{ + + template + T powm1(T x, T y); + + }} // namespaces + +[h4 Description] + + template + T log1p(T x); + +Returns the natural logarithm of `x+1`. + +There are many situations where it is desirable to compute `log(x+1)`. +However, for small `x` then `x+1` suffers from catastrophic cancellation errors +so that `x+1 == 1` and `log(x+1) == 0`, when in fact for very small x, the +best approximation to `log(x+1)` would be `x`. `log1p` calculates the best +approximation to `log(1+x)` using a Taylor series expansion for accuracy +(less than __te). +Alternatively note that there are faster methods available, +for example using the equivalence: + + log(1+x) == (log(1+x) * x) / ((1-x) - 1) + +However, experience has shown that these methods tend to fail quite spectacularly +once the compiler's optimizations are turned on, consequently they are +used only when known not to break with a particular compiler. +In contrast, the series expansion method seems to be reasonably +immune to optimizer-induced errors. + +Finally when BOOST_HAS_LOG1P is defined then the `float/double/long double` +specializations of this template simply forward to the platform's +native implementation of this function. + + template + T expm1(T x); + +Returns e[super x] - 1. + +For small x, then __ex is very close to 1, as a result calculating __exm1 results +in catastrophic cancellation errors when x is small. `expm1` calculates __exm1 using +rational approximations (for up to 128-bit long doubles), otherwise via +a series expansion when x is small (giving an accuracy of less than __te). + +Finally when BOOST_HAS_EXPM1 is defined then the `float/double/long double` +specializations of this template simply forward to the platform's +native implementation of this function. + + template + T cbrt(T x); + +Returns the cubed root of x. + +Implemented using Halley iteration. + + template + T sqrt1pm1(T x); + +Returns `sqrt(1+x) - 1`. + +This function is useful when you need the difference between sqrt(x) and 1, when +x is itself close to 1. + +Implemented in terms of `log1p` and `expm1`. + + template + T powm1(T x, T y); + +Returns x[super y ] - 1. + +There are two domains where this is useful: when y is very small, or when +x is close to 1. + +Implemented in terms of `expm1`. + +[h4 Accuracy] + +For built in floating point types `expm1`, `log1p` and `cbrt` +should have approximately 1 epsilon accuracy, +the other functions about double that. + +[h4 Testing] + +A mixture of spot test sanity checks, and random high precision test values +calculated using NTL::RR at 1000-bit precision. + +[endsect][/section:beta_function The Beta Function] +[/ + Copyright 2006 John Maddock and Paul A. Bristow. + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + diff --git a/include/boost/math/special_functions/gamma.hpp b/include/boost/math/special_functions/gamma.hpp index 999ceca55..e0a05d23a 100644 --- a/include/boost/math/special_functions/gamma.hpp +++ b/include/boost/math/special_functions/gamma.hpp @@ -27,7 +27,7 @@ #include #include #include -#include +#include #include #include #include @@ -910,7 +910,7 @@ T tgammap1m1_imp(T dz, const L&) T zgh = (L::g() + T(0.5) + dz) / boost::math::constants::e(); T A = boost::math::powm1(zgh, dz); - T B = boost::math::sqrtp1m1(dz / (L::g() + T(0.5))); + T B = boost::math::sqrt1pm1(dz / (L::g() + T(0.5))); T C = L::lanczos_sum_near_1(dz); T Ap1 = pow(zgh, dz); T Bp1 = sqrt(1 + (dz / (L::g() + T(0.5))) ); @@ -1544,7 +1544,7 @@ inline T lgamma(T x) } template -inline T tgammap1m1(T z) +inline T tgamma1pm1(T z) { BOOST_FPU_EXCEPTION_GUARD typedef typename lanczos::lanczos_traits::type>::value_type value_type; diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp index f7c8a9d6a..8c1deb706 100644 --- a/include/boost/math/special_functions/math_fwd.hpp +++ b/include/boost/math/special_functions/math_fwd.hpp @@ -171,7 +171,7 @@ namespace boost // sqrt template - T sqrtp1m1(const T&); + T sqrt1pm1(const T&); } // namespace math } // namespace boost diff --git a/include/boost/math/special_functions/powm1.hpp b/include/boost/math/special_functions/powm1.hpp index ec10c4700..cbf4ca191 100644 --- a/include/boost/math/special_functions/powm1.hpp +++ b/include/boost/math/special_functions/powm1.hpp @@ -6,81 +6,25 @@ #ifndef BOOST_MATH_POWM1 #define BOOST_MATH_POWM1 -#include -#include +#include +#include #include -// -// This algorithm computes (x^y)-1, it's only called for small -// values of x and y, in fact x and y < 1, any other argument -// probably is not supported, and will not yield accurate results. -// Probably best not to use this unless you know that the conditions -// above are satisfied. -// - -namespace boost{ namespace math{ namespace detail{ +namespace boost{ namespace math{ template -struct powm1_series -{ - typedef T result_type; - powm1_series(T z, T a_) : k(1)/*, a(a_)*/ - { - using namespace std; - result = lz = log(z) * a_; - } - T operator()() - { - T r = result; - result *= lz / ++k; - return r; - } -private: - T result, lz/*, a*/; - int k; -}; - -template -struct powp1m1_series -{ - typedef T result_type; - - powp1m1_series(T z_, T a_) : k(1), result(z_*a_), z(z_), a(a_){} - T operator()() - { - T r = result; - result *= (a-k)*z; - result /= ++k; - return r; - } -private: - int k; - T result, z, a; -}; - -} // namespace detail - -template -T powm1(const T a, const T z) +inline T powm1(const T a, const T z) { using namespace std; - T result = pow(a, z) - 1; - - if(fabs(result) < T(0.5)) + if((fabs(a) < 1) || (fabs(z) < 1)) { - if(fabs(a-1) < fabs(z)) - { - detail::powp1m1_series gen(a-1, z); - result = tools::kahan_sum_series(gen, ::boost::math::tools::digits()); - } - else - { - detail::powm1_series gen(a, z); - result = tools::kahan_sum_series(gen, ::boost::math::tools::digits()); - } + T p = log(a) * z; + if(fabs(p) < 2) + return boost::math::expm1(p); + // otherwise fall though: } - return result; + return pow(a, z) - 1; } } // namespace math diff --git a/include/boost/math/special_functions/sqrt1pm1.hpp b/include/boost/math/special_functions/sqrt1pm1.hpp new file mode 100644 index 000000000..291a0def0 --- /dev/null +++ b/include/boost/math/special_functions/sqrt1pm1.hpp @@ -0,0 +1,34 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SQRT1PM1 +#define BOOST_MATH_SQRT1PM1 + +#include +#include + +// +// This algorithm computes sqrt(1+x)-1 for small x: +// + +namespace boost{ namespace math{ + +template +inline T sqrt1pm1(const T& val) +{ + using namespace std; + + if(fabs(val) > 0.75) + return sqrt(1 + val) - 1; + return boost::math::expm1(boost::math::log1p(val) / 2); +} + +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_SQRT1PM1 + + + diff --git a/include/boost/math/special_functions/sqrtp1m1.hpp b/include/boost/math/special_functions/sqrtp1m1.hpp deleted file mode 100644 index 881685ce3..000000000 --- a/include/boost/math/special_functions/sqrtp1m1.hpp +++ /dev/null @@ -1,58 +0,0 @@ -// (C) Copyright John Maddock 2006. -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. (See accompanying file -// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_MATH_SQRTP1M1 -#define BOOST_MATH_SQRTP1M1 - -#include -#include - -// -// This algorithm computes sqrt(1+x)-1 for small x only: -// - -namespace boost{ namespace math{ - -namespace detail -{ - - template - struct sqrtp1m1_series - { - typedef T result_type; - sqrtp1m1_series(T z_) : result(z_/2), z(z_), k(1){} - T operator()() - { - T r = result; - result *= z * (k - T(0.5)); - ++k; - result /= -k; - return r; - } - private: - T result, z; - int k; - }; - -} // namespace detail - -template -T sqrtp1m1(const T& val) -{ - using namespace std; - - if(fabs(val) > 0.75) - return sqrt(1 + val) - 1; - detail::sqrtp1m1_series gen(val); - return tools::kahan_sum_series(gen, ::boost::math::tools::digits()); -} - -} // namespace math -} // namespace boost - -#endif // BOOST_MATH_SQRTP1M1 - - - diff --git a/test/powm1_sqrtp1m1_test.cpp b/test/powm1_sqrtp1m1_test.cpp index 413cecdab..64470fed4 100644 --- a/test/powm1_sqrtp1m1_test.cpp +++ b/test/powm1_sqrtp1m1_test.cpp @@ -6,7 +6,7 @@ #include #include #include -#include +#include #include #include #include @@ -19,7 +19,7 @@ // DESCRIPTION: // ~~~~~~~~~~~~ // -// This file tests the functions powm1 and sqrtp1m1. +// This file tests the functions powm1 and sqrt1pm1. // The accuracy tests // use values generated with NTL::RR at 1000-bit precision // and our generic versions of these functions. @@ -1617,7 +1617,7 @@ void test_powm1_sqrtp1m1(T, const char* type_name) using namespace std; typedef T (*func_t)(const T&); - func_t f = &boost::math::sqrtp1m1; + func_t f = &boost::math::sqrt1pm1; boost::math::tools::test_result result = boost::math::tools::test( sqrtp1m1_data, @@ -1626,7 +1626,7 @@ void test_powm1_sqrtp1m1(T, const char* type_name) std::cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n" "Test results for type " << type_name << std::endl << std::endl; - handle_test_result(result, sqrtp1m1_data[result.worst()], result.worst(), type_name, "boost::math::sqrtp1m1", "sqrtp1m1"); + handle_test_result(result, sqrtp1m1_data[result.worst()], result.worst(), type_name, "boost::math::sqrt1pm1", "sqrt1pm1"); typedef T (*func2_t)(T, T); func2_t f2 = &boost::math::powm1; diff --git a/test/test_gamma.cpp b/test/test_gamma.cpp index fd35bd548..dcf0dc206 100644 --- a/test/test_gamma.cpp +++ b/test/test_gamma.cpp @@ -23,7 +23,7 @@ // ~~~~~~~~~~~~ // // This file tests the functions tgamma and lgamma, and the -// function tgammap1m1. There are two sets of tests, spot +// function tgamma1pm1. There are two sets of tests, spot // tests which compare our results with selected values computed // using the online special function calculator at // functions.wolfram.com, while the bulk of the accuracy tests @@ -100,8 +100,8 @@ void expected_results() ".*", // stdlib "linux", // platform largest_type, // test type(s) - "tgammap1m1.*", // test data group - "boost::math::tgammap1m1", 50, 15); // test function + "tgamma1pm1.*", // test data group + "boost::math::tgamma1pm1", 50, 15); // test function add_expected_result( ".*", // compiler ".*", // stdlib @@ -121,8 +121,8 @@ void expected_results() ".*", // stdlib "linux", // platform "real_concept", // test type(s) - "tgammap1m1.*", // test data group - "boost::math::tgammap1m1", 40, 10); // test function + "tgamma1pm1.*", // test data group + "boost::math::tgamma1pm1", 40, 10); // test function // // HP-UX results: // @@ -159,8 +159,8 @@ void expected_results() ".*", // stdlib "HP-UX", // platform largest_type, // test type(s) - "tgammap1m1.*", // test data group - "boost::math::tgammap1m1", 200, 70); // test function + "tgamma1pm1.*", // test data group + "boost::math::tgamma1pm1", 200, 70); // test function add_expected_result( ".*", // compiler ".*", // stdlib @@ -173,8 +173,8 @@ void expected_results() ".*", // stdlib "HP-UX", // platform "real_concept", // test type(s) - "tgammap1m1.*", // test data group - "boost::math::tgammap1m1", 200, 60); // test function + "tgamma1pm1.*", // test data group + "boost::math::tgamma1pm1", 200, 60); // test function // // Catch all cases come last: @@ -219,8 +219,8 @@ void expected_results() ".*", // stdlib ".*", // platform largest_type, // test type(s) - "tgammap1m1.*", // test data group - "boost::math::tgammap1m1", 30, 9); // test function + "tgamma1pm1.*", // test data group + "boost::math::tgamma1pm1", 30, 9); // test function add_expected_result( ".*", // compiler @@ -255,8 +255,8 @@ void expected_results() ".*", // stdlib ".*", // platform "real_concept", // test type(s) - "tgammap1m1.*", // test data group - "boost::math::tgammap1m1", 20, 5); // test function + "tgamma1pm1.*", // test data group + "boost::math::tgamma1pm1", 20, 5); // test function // // Finish off by printing out the compiler/stdlib/platform names, @@ -328,7 +328,7 @@ void do_test_gammap1m1(const T& data, const char* type_name, const char* test_na typedef typename row_type::value_type value_type; typedef value_type (*pg)(value_type); - pg funcp = boost::math::tgammap1m1; + pg funcp = boost::math::tgamma1pm1; boost::math::tools::test_result result; @@ -336,13 +336,13 @@ void do_test_gammap1m1(const T& data, const char* type_name, const char* test_na << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"; // - // test tgammap1m1 against data: + // test tgamma1pm1 against data: // result = boost::math::tools::test( data, boost::lambda::bind(funcp, boost::lambda::ret(boost::lambda::_1[0])), boost::lambda::ret(boost::lambda::_1[1])); - handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::tgammap1m1", test_name); + handle_test_result(result, data[result.worst()], result.worst(), type_name, "boost::math::tgamma1pm1", test_name); std::cout << std::endl; } @@ -389,9 +389,9 @@ void test_gamma(T, const char* name) do_test_gamma(near_m55, name, "near -55"); // - // And now tgammap1m1 which computes gamma(1+dz)-1: + // And now tgamma1pm1 which computes gamma(1+dz)-1: // - do_test_gammap1m1(gammap1m1_data, name, "tgammap1m1(dz)"); + do_test_gammap1m1(gammap1m1_data, name, "tgamma1pm1(dz)"); } template