boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Improved and generally tidied up sqrt1pm1 and powm1.

Browse Source 
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]
This commit is contained in:
John Maddock
2006-10-26 11:51:59 +00:00
parent 527bf0d227
commit 82548699a5
11 changed files with 316 additions and 150 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/dist_reference.qbk
View File

@@ -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]

90
doc/distributions/weibull.qbk Normal file
View File

@@ -0,0 +1,90 @@
[section:weibull Weibull Distribution]
[template alpha '''α''']
[template beta '''β''']
[template GAMMA '''Γ''']
``#include <boost/math/distributions/weibull.hpp>``
namespace boost{ namespace math{
template <class RealType>
class weibull_distribution;
typedef weibull_distribution<double> weibull;
template <class RealType>
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]

6
doc/math.qbk
View File

@@ -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]

149
doc/powers.qbk Normal file
View File

@@ -0,0 +1,149 @@
[section:powers Logs, Powers, Roots and Exponentials]
[caution __caution ]
[h4 Synopsis]
``
#include <boost/math/special_functions/log1p.hpp>
``
namespace boost{ namespace math{
template <class T>
T log1p(T x);
}} // namespaces
``
#include <boost/math/special_functions/expm1.hpp>
``
namespace boost{ namespace math{
template <class T>
T expm1(T x);
}} // namespaces
``
#include <boost/math/special_functions/cbrt.hpp>
``
namespace boost{ namespace math{
template <class T>
T cbrt(T x);
}} // namespaces
``
#include <boost/math/special_functions/sqrt1pm1.hpp>
``
namespace boost{ namespace math{
template <class T>
T sqrt1pm1(T x);
}} // namespaces
``
#include <boost/math/special_functions/powm1.hpp>
``
namespace boost{ namespace math{
template <class T>
T powm1(T x, T y);
}} // namespaces
[h4 Description]
template <class T>
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 <class T>
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 <class T>
T cbrt(T x);
Returns the cubed root of x.
Implemented using Halley iteration.
template <class T>
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 <class T>
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).
]

6
include/boost/math/special_functions/gamma.hpp
View File

@@ -27,7 +27,7 @@
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/powm1.hpp>
#include <boost/math/special_functions/sqrtp1m1.hpp>
#include <boost/math/special_functions/sqrt1pm1.hpp>
#include <boost/math/special_functions/lanczos.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/math/special_functions/detail/igamma_large.hpp>
@@ -910,7 +910,7 @@ T tgammap1m1_imp(T dz, const L&)
T zgh = (L::g() + T(0.5) + dz) / boost::math::constants::e<T>();
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 <class T>
inline T tgammap1m1(T z)
inline T tgamma1pm1(T z)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename lanczos::lanczos_traits<typename remove_cv<T>::type>::value_type value_type;

2
include/boost/math/special_functions/math_fwd.hpp
View File

@@ -171,7 +171,7 @@ namespace boost
// sqrt
template <class T>
T sqrtp1m1(const T&);
T sqrt1pm1(const T&);
} // namespace math
} // namespace boost

76
include/boost/math/special_functions/powm1.hpp
View File

@@ -6,81 +6,25 @@
#ifndef BOOST_MATH_POWM1
#define BOOST_MATH_POWM1
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
#include <boost/assert.hpp>
//
// 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 <class T>
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 <class T>
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 <class T>
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<T> gen(a-1, z);
result = tools::kahan_sum_series(gen, ::boost::math::tools::digits<T>());
}
else
{
detail::powm1_series<T> gen(a, z);
result = tools::kahan_sum_series(gen, ::boost::math::tools::digits<T>());
}
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

34
include/boost/math/special_functions/sqrt1pm1.hpp Normal file
View File

@@ -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 <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
//
// This algorithm computes sqrt(1+x)-1 for small x:
//
namespace boost{ namespace math{
template <class T>
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

58
include/boost/math/special_functions/sqrtp1m1.hpp
View File

@@ -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 <boost/math/tools/series.hpp>
#include <boost/math/tools/precision.hpp>
//
// This algorithm computes sqrt(1+x)-1 for small x only:
//
namespace boost{ namespace math{
namespace detail
{
template <class T>
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 <class T>
T sqrtp1m1(const T& val)
{
using namespace std;
if(fabs(val) > 0.75)
return sqrt(1 + val) - 1;
detail::sqrtp1m1_series<T> gen(val);
return tools::kahan_sum_series(gen, ::boost::math::tools::digits<T>());
}
} // namespace math
} // namespace boost
#endif // BOOST_MATH_SQRTP1M1

8
test/powm1_sqrtp1m1_test.cpp
View File

@@ -6,7 +6,7 @@
#include <boost/math/concepts/real_concept.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/sqrtp1m1.hpp>
#include <boost/math/special_functions/sqrt1pm1.hpp>
#include <boost/math/special_functions/powm1.hpp>
#include <boost/math/tools/test.hpp>
#include <boost/lambda/lambda.hpp>
@@ -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<T>;
func_t f = &boost::math::sqrt1pm1<T>;
boost::math::tools::test_result<T> 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<T>;

36
test/test_gamma.cpp
View File

@@ -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<value_type> 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<value_type>(boost::lambda::_1[0])),
boost::lambda::ret<value_type>(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 <class T>