|
|
Boost-Commit : |
From: john_at_[hidden]
Date: 2008-01-14 04:27:52
Author: johnmaddock
Date: 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
New Revision: 42744
URL: http://svn.boost.org/trac/boost/changeset/42744
Log:
Change concepts classes so they can be used when there is no long double support.
Added first cut of the non-central chi squared distribution.
Removed almost all occurrences of real_cast: replaced with calls to the truncation/rounding functions instead.
Added:
sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/test_nccs_hooks.hpp (contents, props changed)
Text files modified:
sandbox/math_toolkit/boost/math/concepts/real_concept.hpp | 169 ++++++++++++++++++---------------------
sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp | 137 ++++++++++++++------------------
sandbox/math_toolkit/boost/math/distributions.hpp | 1
sandbox/math_toolkit/boost/math/distributions/detail/common_error_handling.hpp | 17 ++++
sandbox/math_toolkit/boost/math/distributions/poisson.hpp | 3
sandbox/math_toolkit/boost/math/special_functions/bessel.hpp | 16 ++-
sandbox/math_toolkit/boost/math/special_functions/beta.hpp | 7
sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp | 4
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp | 5
sandbox/math_toolkit/boost/math/special_functions/factorials.hpp | 2
sandbox/math_toolkit/boost/math/special_functions/gamma.hpp | 8
sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp | 4
sandbox/math_toolkit/boost/math/tools/polynomial.hpp | 5
sandbox/math_toolkit/libs/math/performance/distributions.cpp | 16 ---
sandbox/math_toolkit/libs/math/test/functor.hpp | 12 +-
sandbox/math_toolkit/libs/math/test/test_bessel_i.cpp | 3
sandbox/math_toolkit/libs/math/test/test_bessel_j.cpp | 2
sandbox/math_toolkit/libs/math/test/test_bessel_k.cpp | 2
sandbox/math_toolkit/libs/math/test/test_bessel_y.cpp | 2
sandbox/math_toolkit/libs/math/test/test_binomial_coeff.cpp | 5
sandbox/math_toolkit/libs/math/test/test_expint.cpp | 5
sandbox/math_toolkit/libs/math/test/test_gamma_dist.cpp | 2
24 files changed, 210 insertions(+), 221 deletions(-)
Modified: sandbox/math_toolkit/boost/math/concepts/real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/real_concept.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -22,7 +22,9 @@
#include <boost/config.hpp>
#include <boost/limits.hpp>
-#include <boost/math/tools/real_cast.hpp>
+#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/precision.hpp>
#include <boost/math/policies/policy.hpp>
#include <ostream>
@@ -38,6 +40,12 @@
namespace concepts
{
+#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ typedef double real_concept_base_type;
+#else
+ typedef long double real_concept_base_type;
+#endif
+
class real_concept
{
public:
@@ -56,11 +64,11 @@
real_concept(unsigned long c) : m_value(c){}
real_concept(long c) : m_value(c){}
#if defined(BOOST_HAS_LONG_LONG) || defined(__DECCXX) || defined(__SUNPRO_CC)
- real_concept(unsigned long long c) : m_value(static_cast<long double>(c)){}
- real_concept(long long c) : m_value(static_cast<long double>(c)){}
+ real_concept(unsigned long long c) : m_value(static_cast<real_concept_base_type>(c)){}
+ real_concept(long long c) : m_value(static_cast<real_concept_base_type>(c)){}
#elif defined(BOOST_HAS_MS_INT64)
- real_concept(unsigned __int64 c) : m_value(static_cast<long double>(c)){}
- real_concept(__int64 c) : m_value(static_cast<long double>(c)){}
+ real_concept(unsigned __int64 c) : m_value(static_cast<real_concept_base_type>(c)){}
+ real_concept(__int64 c) : m_value(static_cast<real_concept_base_type>(c)){}
#endif
real_concept(float c) : m_value(c){}
real_concept(double c) : m_value(c){}
@@ -80,15 +88,15 @@
real_concept& operator=(long c) { m_value = c; return *this; }
real_concept& operator=(unsigned long c) { m_value = c; return *this; }
#ifdef BOOST_HAS_LONG_LONG
- real_concept& operator=(long long c) { m_value = static_cast<long double>(c); return *this; }
- real_concept& operator=(unsigned long long c) { m_value = static_cast<long double>(c); return *this; }
+ real_concept& operator=(long long c) { m_value = static_cast<real_concept_base_type>(c); return *this; }
+ real_concept& operator=(unsigned long long c) { m_value = static_cast<real_concept_base_type>(c); return *this; }
#endif
real_concept& operator=(float c) { m_value = c; return *this; }
real_concept& operator=(double c) { m_value = c; return *this; }
real_concept& operator=(long double c) { m_value = c; return *this; }
// Access:
- long double value()const{ return m_value; }
+ real_concept_base_type value()const{ return m_value; }
// Member arithmetic:
real_concept& operator+=(const real_concept& other)
@@ -105,7 +113,7 @@
{ return *this; }
private:
- long double m_value;
+ real_concept_base_type m_value;
};
// Non-member arithmetic:
@@ -148,40 +156,6 @@
inline bool operator >= (const real_concept& a, const real_concept& b)
{ return a.value() >= b.value(); }
-#if 0
-// Non-member mixed compare:
-template <class T>
-inline bool operator == (const T& a, const real_concept& b)
-{
- return a == b.value();
-}
-template <class T>
-inline bool operator != (const T& a, const real_concept& b)
-{
- return a != b.value();
-}
-template <class T>
-inline bool operator < (const T& a, const real_concept& b)
-{
- return a < b.value();
-}
-template <class T>
-inline bool operator > (const T& a, const real_concept& b)
-{
- return a > b.value();
-}
-template <class T>
-inline bool operator <= (const T& a, const real_concept& b)
-{
- return a <= b.value();
-}
-template <class T>
-inline bool operator >= (const T& a, const real_concept& b)
-{
- return a >= b.value();
-}
-#endif // Non-member mixed compare:
-
// Non-member functions:
inline real_concept acos(real_concept a)
{ return std::acos(a.value()); }
@@ -196,8 +170,13 @@
inline real_concept ceil(real_concept a)
{ return std::ceil(a.value()); }
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+// I've seen std::fmod(long double) crash on some platforms
+// so use fmodl instead:
inline real_concept fmod(real_concept a, real_concept b)
{ return fmodl(a.value(), b.value()); }
+#else
+inline real_concept fmod(real_concept a, real_concept b)
+{ return std::fmod(a.value(), b.value()); }
#endif
inline real_concept cosh(real_concept a)
{ return std::cosh(a.value()); }
@@ -211,8 +190,8 @@
{ return std::floor(a.value()); }
inline real_concept modf(real_concept a, real_concept* ipart)
{
- long double ip;
- long double result = std::modf(a.value(), &ip);
+ real_concept_base_type ip;
+ real_concept_base_type result = std::modf(a.value(), &ip);
*ipart = ip;
return result;
}
@@ -233,7 +212,7 @@
{ return std::pow(a.value(), b); }
#else
inline real_concept pow(real_concept a, int b)
-{ return std::pow(a.value(), static_cast<long double>(b)); }
+{ return std::pow(a.value(), static_cast<real_concept_base_type>(b)); }
#endif
inline real_concept sin(real_concept a)
{ return std::sin(a.value()); }
@@ -264,13 +243,13 @@
return is;
#elif defined(__SGI_STL_PORT)
std::string s;
- long double d;
+ real_concept_base_type d;
is >> s;
std::sscanf(s.c_str(), "%Lf", &d);
a = d;
return is;
#else
- long double v;
+ real_concept_base_type v;
is >> v;
a = v;
return is;
@@ -281,94 +260,102 @@
namespace tools
{
-// real_cast converts from T to integer and narrower floating-point types.
-
-// Convert from T to integer types.
template <>
-inline unsigned int real_cast<unsigned int, concepts::real_concept>(concepts::real_concept r)
+inline concepts::real_concept max_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
{
- return static_cast<unsigned int>(r.value());
+ return max_value<concepts::real_concept_base_type>();
}
template <>
-inline int real_cast<int, concepts::real_concept>(concepts::real_concept r)
+inline concepts::real_concept min_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
{
- return static_cast<int>(r.value());
+ return min_value<concepts::real_concept_base_type>();
}
template <>
-inline long real_cast<long, concepts::real_concept>(concepts::real_concept r)
+inline concepts::real_concept log_max_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
{
- return static_cast<long>(r.value());
+ return log_max_value<concepts::real_concept_base_type>();
}
-// Converts from T to narrower floating-point types, float, double & long double.
-
template <>
-inline float real_cast<float, concepts::real_concept>(concepts::real_concept r)
+inline concepts::real_concept log_min_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
{
- return static_cast<float>(r.value());
+ return log_min_value<concepts::real_concept_base_type>();
}
+
template <>
-inline double real_cast<double, concepts::real_concept>(concepts::real_concept r)
+inline concepts::real_concept epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
{
- return static_cast<double>(r.value());
+#ifdef __SUNPRO_CC
+ return std::numeric_limits<concepts::real_concept_base_type>::epsilon();
+#else
+ return tools::epsilon<concepts::real_concept_base_type>();
+#endif
}
+
template <>
-inline long double real_cast<long double, concepts::real_concept>(concepts::real_concept r)
-{
- return r.value();
+inline int digits<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+{
+ // Assume number of significand bits is same as real_concept_base_type,
+ // unless std::numeric_limits<T>::is_specialized to provide digits.
+ return tools::digits<concepts::real_concept_base_type>();
+ // Note that if numeric_limits real concept is NOT specialized to provide digits10
+ // (or max_digits10) then the default precision of 6 decimal digits will be used
+ // by Boost test (giving misleading error messages like
+ // "difference between {9.79796} and {9.79796} exceeds 5.42101e-19%"
+ // and by Boost lexical cast and serialization causing loss of accuracy.
}
+} // namespace tools
+
+//
+// Conversion and truncation routines:
+//
template <>
-inline concepts::real_concept max_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+inline int iround<concepts::real_concept>(const concepts::real_concept& v)
{
- return max_value<long double>();
+ return iround(v.value());
}
template <>
-inline concepts::real_concept min_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+inline long lround<concepts::real_concept>(const concepts::real_concept& v)
{
- return min_value<long double>();
+ return lround(v.value());
}
+#ifdef BOOST_HAS_LONG_LONG
+
template <>
-inline concepts::real_concept log_max_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+inline long long llround<concepts::real_concept>(const concepts::real_concept& v)
{
- return log_max_value<long double>();
+ return llround(v.value());
}
+#endif
+
template <>
-inline concepts::real_concept log_min_value<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+inline int itrunc<concepts::real_concept>(const concepts::real_concept& v)
{
- return log_min_value<long double>();
+ return itrunc(v.value());
}
template <>
-inline concepts::real_concept epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
+inline long ltrunc<concepts::real_concept>(const concepts::real_concept& v)
{
-#ifdef __SUNPRO_CC
- return std::numeric_limits<long double>::epsilon();
-#else
- return tools::epsilon<long double>();
-#endif
+ return ltrunc(v.value());
}
+#ifdef BOOST_HAS_LONG_LONG
+
template <>
-inline int digits<concepts::real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::real_concept))
-{
- // Assume number of significand bits is same as long double,
- // unless std::numeric_limits<T>::is_specialized to provide digits.
- return tools::digits<long double>();
- // Note that if numeric_limits real concept is NOT specialized to provide digits10
- // (or max_digits10) then the default precision of 6 decimal digits will be used
- // by Boost test (giving misleading error messages like
- // "difference between {9.79796} and {9.79796} exceeds 5.42101e-19%"
- // and by Boost lexical cast and serialization causing loss of accuracy.
+inline long long lltrunc<concepts::real_concept>(const concepts::real_concept& v)
+{
+ return lltrunc(v.value());
}
-} // namespace tools
+#endif
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp (original)
+++ sandbox/math_toolkit/boost/math/concepts/std_real_concept.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -18,9 +18,11 @@
#include <boost/config.hpp>
#include <boost/limits.hpp>
-#include <boost/math/tools/real_cast.hpp>
#include <boost/math/tools/precision.hpp>
#include <boost/math/policies/policy.hpp>
+#include <boost/math/special_functions/round.hpp>
+#include <boost/math/special_functions/trunc.hpp>
+#include <boost/math/special_functions/modf.hpp>
#include <ostream>
#include <istream>
@@ -35,6 +37,12 @@
namespace concepts
{
+#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ typedef double std_real_concept_base_type;
+#else
+ typedef long double std_real_concept_base_type;
+#endif
+
class std_real_concept
{
public:
@@ -53,8 +61,8 @@
std_real_concept(unsigned long c) : m_value(c){}
std_real_concept(long c) : m_value(c){}
#if defined(BOOST_HAS_LONG_LONG) || defined(__DECCXX) || defined(__SUNPRO_CC)
- std_real_concept(unsigned long long c) : m_value(static_cast<long double>(c)){}
- std_real_concept(long long c) : m_value(static_cast<long double>(c)){}
+ std_real_concept(unsigned long long c) : m_value(static_cast<std_real_concept_base_type>(c)){}
+ std_real_concept(long long c) : m_value(static_cast<std_real_concept_base_type>(c)){}
#endif
std_real_concept(float c) : m_value(c){}
std_real_concept(double c) : m_value(c){}
@@ -74,15 +82,15 @@
std_real_concept& operator=(long c) { m_value = c; return *this; }
std_real_concept& operator=(unsigned long c) { m_value = c; return *this; }
#if defined(BOOST_HAS_LONG_LONG) || defined(__DECCXX) || defined(__SUNPRO_CC)
- std_real_concept& operator=(long long c) { m_value = static_cast<long double>(c); return *this; }
- std_real_concept& operator=(unsigned long long c) { m_value = static_cast<long double>(c); return *this; }
+ std_real_concept& operator=(long long c) { m_value = static_cast<std_real_concept_base_type>(c); return *this; }
+ std_real_concept& operator=(unsigned long long c) { m_value = static_cast<std_real_concept_base_type>(c); return *this; }
#endif
std_real_concept& operator=(float c) { m_value = c; return *this; }
std_real_concept& operator=(double c) { m_value = c; return *this; }
std_real_concept& operator=(long double c) { m_value = c; return *this; }
// Access:
- long double value()const{ return m_value; }
+ std_real_concept_base_type value()const{ return m_value; }
// Member arithmetic:
std_real_concept& operator+=(const std_real_concept& other)
@@ -99,7 +107,7 @@
{ return *this; }
private:
- long double m_value;
+ std_real_concept_base_type m_value;
};
// Non-member arithmetic:
@@ -142,40 +150,6 @@
inline bool operator >= (const std_real_concept& a, const std_real_concept& b)
{ return a.value() >= b.value(); }
-#if 0
-// Non-member mixed compare:
-template <class T>
-inline bool operator == (const T& a, const std_real_concept& b)
-{
- return a == b.value();
-}
-template <class T>
-inline bool operator != (const T& a, const std_real_concept& b)
-{
- return a != b.value();
-}
-template <class T>
-inline bool operator < (const T& a, const std_real_concept& b)
-{
- return a < b.value();
-}
-template <class T>
-inline bool operator > (const T& a, const std_real_concept& b)
-{
- return a > b.value();
-}
-template <class T>
-inline bool operator <= (const T& a, const std_real_concept& b)
-{
- return a <= b.value();
-}
-template <class T>
-inline bool operator >= (const T& a, const std_real_concept& b)
-{
- return a >= b.value();
-}
-#endif // Non-member mixed compare:
-
} // namespace concepts
} // namespace math
} // namespace boost
@@ -214,8 +188,8 @@
{ return std::floor(a.value()); }
inline boost::math::concepts::std_real_concept modf(boost::math::concepts::std_real_concept a, boost::math::concepts::std_real_concept* ipart)
{
- long double ip;
- long double result = std::modf(a.value(), &ip);
+ boost::math::concepts::std_real_concept_base_type ip;
+ boost::math::concepts::std_real_concept_base_type result = std::modf(a.value(), &ip);
*ipart = ip;
return result;
}
@@ -255,7 +229,7 @@
template <class charT, class traits>
inline std::basic_istream<charT, traits>& operator>>(std::basic_istream<charT, traits>& is, std_real_concept& a)
{
- long double v;
+ std_real_concept_base_type v;
is >> v;
a = v;
return is;
@@ -265,84 +239,93 @@
namespace tools
{
-// real_cast converts from T to integer and narrower floating-point types.
-
-// Convert from T to integer types.
template <>
-inline unsigned int real_cast<unsigned int, concepts::std_real_concept>(concepts::std_real_concept r)
+inline concepts::std_real_concept max_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
{
- return static_cast<unsigned int>(r.value());
+ return max_value<concepts::std_real_concept_base_type>();
}
template <>
-inline int real_cast<int, concepts::std_real_concept>(concepts::std_real_concept r)
+inline concepts::std_real_concept min_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
{
- return static_cast<int>(r.value());
+ return min_value<concepts::std_real_concept_base_type>();
}
template <>
-inline long real_cast<long, concepts::std_real_concept>(concepts::std_real_concept r)
+inline concepts::std_real_concept log_max_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
{
- return static_cast<long>(r.value());
+ return log_max_value<concepts::std_real_concept_base_type>();
}
-// Converts from T to narrower floating-point types, float, double & long double.
-
template <>
-inline float real_cast<float, concepts::std_real_concept>(concepts::std_real_concept r)
+inline concepts::std_real_concept log_min_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
{
- return static_cast<float>(r.value());
+ return log_min_value<concepts::std_real_concept_base_type>();
}
+
template <>
-inline double real_cast<double, concepts::std_real_concept>(concepts::std_real_concept r)
+inline concepts::std_real_concept epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
{
- return static_cast<double>(r.value());
+ return tools::epsilon<concepts::std_real_concept_base_type>();
}
+
template <>
-inline long double real_cast<long double, concepts::std_real_concept>(concepts::std_real_concept r)
-{
- return r.value();
+inline int digits<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+{ // Assume number of significand bits is same as std_real_concept_base_type,
+ // unless std::numeric_limits<T>::is_specialized to provide digits.
+ return digits<concepts::std_real_concept_base_type>();
}
+} // namespace tools
+
+//
+// Conversion and truncation routines:
+//
template <>
-inline concepts::std_real_concept max_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+inline int iround<concepts::std_real_concept>(const concepts::std_real_concept& v)
{
- return max_value<long double>();
+ return iround(v.value());
}
template <>
-inline concepts::std_real_concept min_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+inline long lround<concepts::std_real_concept>(const concepts::std_real_concept& v)
{
- return min_value<long double>();
+ return lround(v.value());
}
+#ifdef BOOST_HAS_LONG_LONG
+
template <>
-inline concepts::std_real_concept log_max_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+inline long long llround<concepts::std_real_concept>(const concepts::std_real_concept& v)
{
- return log_max_value<long double>();
+ return llround(v.value());
}
+#endif
+
template <>
-inline concepts::std_real_concept log_min_value<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+inline int itrunc<concepts::std_real_concept>(const concepts::std_real_concept& v)
{
- return log_min_value<long double>();
+ return itrunc(v.value());
}
template <>
-inline concepts::std_real_concept epsilon(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+inline long ltrunc<concepts::std_real_concept>(const concepts::std_real_concept& v)
{
- return tools::epsilon<long double>();
+ return ltrunc(v.value());
}
+#ifdef BOOST_HAS_LONG_LONG
+
template <>
-inline int digits<concepts::std_real_concept>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
-{ // Assume number of significand bits is same as long double,
- // unless std::numeric_limits<T>::is_specialized to provide digits.
- return digits<long double>();
+inline long long lltrunc<concepts::std_real_concept>(const concepts::std_real_concept& v)
+{
+ return lltrunc(v.value());
}
-} // namespace tools
+#endif
+
} // namespace math
} // namespace boost
Modified: sandbox/math_toolkit/boost/math/distributions.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -24,6 +24,7 @@
#include <boost/math/distributions/gamma.hpp>
#include <boost/math/distributions/lognormal.hpp>
#include <boost/math/distributions/negative_binomial.hpp>
+#include <boost/math/distributions/non_central_chi_squared.hpp>
#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/pareto.hpp>
#include <boost/math/distributions/poisson.hpp>
Modified: sandbox/math_toolkit/boost/math/distributions/detail/common_error_handling.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/detail/common_error_handling.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/detail/common_error_handling.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -96,6 +96,23 @@
// leaving this test to catch any NaNs. see Normal and cauchy for example.
}
+template <class RealType, class Policy>
+inline bool check_non_centrality(
+ const char* function,
+ RealType ncp,
+ RealType* result,
+ const Policy& pol)
+{
+ if((ncp < 0) || !(boost::math::isfinite)(ncp))
+ { // Assume scale == 0 is NOT valid for any distribution.
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "Non centrality parameter is %1%, but must be > 0 !", ncp, pol);
+ return false;
+ }
+ return true;
+}
+
} // namespace detail
} // namespace math
} // namespace boost
Added: sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,644 @@
+// boost\math\distributions\non_central_chi_squared.hpp
+
+// Copyright John Maddock 2008.
+
+// 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_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
+#define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
+#include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
+#include <boost/math/special_functions/round.hpp> // for iround
+#include <boost/math/distributions/complement.hpp> // complements
+#include <boost/math/distributions/chi_squared.hpp> // central distribution
+#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
+#include <boost/math/special_functions/fpclassify.hpp> // isnan.
+#include <boost/math/tools/roots.hpp> // for root finding.
+
+namespace boost
+{
+ namespace math
+ {
+
+ template <class RealType, class Policy>
+ class non_central_chi_squared_distribution;
+
+ namespace detail{
+
+ template <class T, class Policy>
+ T non_central_chi_square_q(T x, T f, T theta, const Policy& pol)
+ {
+ //
+ // Computes the complement of the Non-Central Chi-Square
+ // Distribution CDF by summing a weighted sum of complements
+ // of the central-distributions. The weighting factor is
+ // a Poisson Distribution.
+ //
+ // This is an application of the technique described in:
+ //
+ // Computing discrete mixtures of continuous
+ // distributions: noncentral chisquare, noncentral t
+ // and the distribution of the square of the sample
+ // multiple correlation coeficient.
+ // D. Benton, K. Krishnamoorthy.
+ // Computational Statistics & Data Analysis 43 (2003) 249 – 267
+ //
+ BOOST_MATH_STD_USING
+ using namespace boost::math;
+
+ //
+ // Initialize the variables we'll be using:
+ //
+ T lambda = theta / 2;
+ T del = f / 2;
+ T y = x / 2;
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
+ T sum = 0;
+ //
+ // k is the starting location for iteration, we'll
+ // move both forwards and backwards from this point.
+ // k is chosen as the peek of the Poisson weights, which
+ // will occur *before* the largest term.
+ //
+ int k = iround(lambda);
+ // Forwards and backwards Poisson weights:
+ T poisf = boost::math::gamma_p_derivative(1 + k, lambda, pol);
+ T poisb = poisf * k / lambda;
+ // Initial forwards central chi squared term:
+ T gamf = boost::math::gamma_q(del + k, y, pol);
+ // Forwards and backwards recursion terms on the central chi squared:
+ T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
+ T xtermb = xtermf * (del + k) / y;
+ // Initial backwards central chi squared term:
+ T gamb = gamf - xtermb;
+
+ //
+ // Forwards iteration first, this is the
+ // stable direction for the gamma function
+ // recurrences:
+ //
+ for(int i = k; i < max_iter; ++i)
+ {
+ T term = poisf * gamf;
+ sum += term;
+ poisf *= lambda / (i + 1);
+ gamf += xtermf;
+ xtermf *= y / (del + i + 1);
+ if((sum == 0) || (term / sum < errtol))
+ break;
+ }
+ //
+ // Now backwards iteration: the gamma
+ // function recurrences are unstable in this
+ // direction, we rely on the terms deminishing in size
+ // faster than we introduce cancellation errors.
+ // For this reason it's very important that we start
+ // *before* the largest term so that backwards iteration
+ // is strictly converging.
+ //
+ for(int i = k - 1; i >= 0; --i)
+ {
+ T term = poisb * gamb;
+ sum += term;
+ poisb *= i / lambda;
+ xtermb *= (del + i) / y;
+ gamb -= xtermb;
+ if((sum == 0) || (term / sum < errtol))
+ break;
+ }
+
+ return sum;
+ }
+
+ template <class T, class Policy>
+ T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol)
+ {
+ //
+ // This is an implementation of:
+ //
+ // Algorithm AS 275:
+ // Computing the Non-Central #2 Distribution Function
+ // Cherng G. Ding
+ // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
+ //
+ // This uses a stable forward iteration to sum the
+ // CDF, unfortunately this can not be used for large
+ // values of the non-centrality parameter because:
+ // * The first term may underfow to zero.
+ // * We may need an extra-ordinary number of terms
+ // before we reach the first *significant* term.
+ //
+ BOOST_MATH_STD_USING
+ using namespace boost::math;
+ T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
+ T lambda = theta / 2;
+ T vk = exp(-lambda);
+ T uk = vk;
+ T sum = tk * vk;
+ if(sum == 0)
+ return sum;
+
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
+
+ for(int i = 1; i < max_iter; ++i)
+ {
+ tk = tk * x / (f + 2 * i);
+ uk = uk * lambda / i;
+ vk = vk + uk;
+ T term = vk * tk;
+ sum += term;
+ if(term / sum < errtol)
+ break;
+ }
+ return sum;
+ }
+
+
+ template <class T, class Policy>
+ T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol)
+ {
+ //
+ // This is taken more or less directly from:
+ //
+ // Computing discrete mixtures of continuous
+ // distributions: noncentral chisquare, noncentral t
+ // and the distribution of the square of the sample
+ // multiple correlation coeficient.
+ // D. Benton, K. Krishnamoorthy.
+ // Computational Statistics & Data Analysis 43 (2003) 249 – 267
+ //
+ // We're summing a Poisson weighting term multiplied by
+ // a central chi squared distribution.
+ //
+ using namespace boost::math;
+ BOOST_MATH_STD_USING
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
+ T errorf, errorb;
+
+ T x = y / 2;
+ T del = lambda / 2;
+ //
+ // Starting location for the iteration, we'll iterate
+ // both forwards and backwards from this point. The
+ // location chosen is the maximum of the Poisson weight
+ // function, which ocurrs *after* the largest term in the
+ // sum.
+ //
+ int k = iround(del);
+ T a = n / 2 + k;
+ // Central chi squared term for forward iteration:
+ T gamkf = boost::math::gamma_p(a, x, pol);
+
+ if(lambda == 0)
+ return gamkf;
+ // Central chi squared term for backward iteration:
+ T gamkb = gamkf;
+ // Forwards Poisson weight:
+ T poiskf = gamma_p_derivative(k+1, del, pol);
+ // Backwards Poisson weight:
+ T poiskb = poiskf;
+ // Forwards gamma function recursion term:
+ T xtermf = boost::math::gamma_p_derivative(a, x, pol);
+ // Backwards gamma function recursion term:
+ T xtermb = xtermf * x / a;
+ T sum = poiskf * gamkf;
+ if(sum == 0)
+ return sum;
+ int i = 1;
+ //
+ // Backwards recursion first, this is the stable
+ // direction for gamma function recurrences:
+ //
+ while(i <= k)
+ {
+ xtermb *= (a - i + 1) / x;
+ gamkb += xtermb;
+ poiskb = poiskb * (k - i + 1) / del;
+ errorb = gamkb * poiskb;
+ sum += errorb;
+ if(errorb / sum < errtol)
+ break;
+ ++i;
+ }
+ i = 1;
+ //
+ // Now forwards recursion, the gamma function
+ // recurrence relation is unstable in this direction,
+ // so we rely on the magnitude of successive terms
+ // decreasing faster than we introduce cancellation error.
+ // For this reason it's vital that k is chosen to be *after*
+ // the largest term, so that successive forward iterations
+ // are strictly (and rapidly) converging.
+ //
+ do
+ {
+ xtermf = xtermf * x / (a + i - 1);
+ gamkf = gamkf - xtermf;
+ poiskf = poiskf * del / (k + i);
+ errorf = poiskf * gamkf;
+ sum += errorf;
+ ++i;
+ }while((errorf / sum > errtol) && (i < max_iter));
+
+ return sum;
+ }
+
+ template <class RealType, class Policy>
+ inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
+ {
+ typedef typename policies::evaluation<RealType, 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;
+
+ value_type result;
+ if(l == 0)
+ result = cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x);
+ else if(x > k + l)
+ {
+ // Complement is the smaller of the two:
+ result = detail::non_central_chi_square_q(
+ static_cast<value_type>(x),
+ static_cast<value_type>(k),
+ static_cast<value_type>(l), forwarding_policy());
+ invert = !invert;
+ }
+ else if(l < 200)
+ {
+ // For small values of the non-centrality parameter
+ // we can use Ding's method:
+ result = detail::non_central_chi_square_p_ding(
+ static_cast<value_type>(x),
+ static_cast<value_type>(k),
+ static_cast<value_type>(l), forwarding_policy());
+ }
+ else
+ {
+ // For largers values of the non-centrality
+ // parameter Ding's method will consume an
+ // extra-ordinary number of terms, and worse
+ // may return zero when the result is in fact
+ // finite, use Krishnamoorthy's method instead:
+ result = detail::non_central_chi_square_p(
+ static_cast<value_type>(x),
+ static_cast<value_type>(k),
+ static_cast<value_type>(l), forwarding_policy());
+ }
+ if(invert)
+ result = 1 - result;
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ result,
+ "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
+ }
+
+ template <class T, class Policy>
+ struct nccs_quantile_functor
+ {
+ nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
+ : dist(d), target(t), comp(c) {}
+
+ T operator()(const T& x)
+ {
+ return comp ?
+ target - cdf(complement(dist, x))
+ : cdf(dist, x) - target;
+ }
+
+ private:
+ non_central_chi_squared_distribution<T,Policy> dist;
+ T target;
+ bool comp;
+ };
+
+ template <class RealType, class Policy>
+ RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
+ {
+ static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
+ typedef typename policies::evaluation<RealType, 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;
+
+ value_type k = dist.degrees_of_freedom();
+ value_type l = dist.non_centrality();
+ value_type r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy())
+ ||
+ !detail::check_probability(
+ function,
+ static_cast<value_type>(p),
+ &r,
+ Policy()))
+ return (RealType)r;
+ //
+ // Special cases first:
+ //
+ if(p == 0)
+ return comp
+ ? policies::raise_overflow_error<RealType>(function, 0, Policy())
+ : 0;
+ if(p == 1)
+ return !comp
+ ? policies::raise_overflow_error<RealType>(function, 0, Policy())
+ : 0;
+
+ value_type b = (l * l) / (k + 3 * l);
+ value_type c = (k + 3 * l) / (k + 2 * l);
+ value_type ff = (k + 2 * l) / (c * c);
+ value_type guess;
+ if(comp)
+ guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
+ else
+ guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
+
+ if(guess < 0)
+ guess = tools::min_value<value_type>();
+
+ detail::nccs_quantile_functor<value_type, Policy>
+ f(non_central_chi_squared_distribution<value_type, Policy>(k, l), p, comp);
+ tools::eps_tolerance<value_type> tol(policies::digits<RealType, Policy>());
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+ std::pair<value_type, value_type> ir = tools::bracket_and_solve_root(
+ f, guess, value_type(2), true, tol, max_iter, Policy());
+ value_type result = ir.first + (ir.second - ir.first) / 2;
+ if(max_iter == policies::get_max_root_iterations<Policy>())
+ {
+ policies::raise_evaluation_error<value_type>(function, "Unable to locate solution in a reasonable time:"
+ " either there is no answer to quantile of the non central chi squared distribution"
+ " or the answer is infinite. Current best guess is %1%", result, Policy());
+ }
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ result,
+ function);
+ }
+
+ }
+
+ template <class RealType = double, class Policy = policies::policy<> >
+ class non_central_chi_squared_distribution
+ {
+ public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
+ {
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
+ RealType r;
+ detail::check_df(
+ function,
+ df, &r, Policy());
+ detail::check_non_centrality(
+ function,
+ ncp,
+ &r,
+ Policy());
+ } // non_central_chi_squared_distribution constructor.
+
+ RealType degrees_of_freedom() const
+ { // Private data getter function.
+ return df;
+ }
+ RealType non_centrality() const
+ { // Private data getter function.
+ return ncp;
+ }
+ private:
+ // Data member, initialized by constructor.
+ RealType df; // degrees of freedom.
+ RealType ncp; // non-centrality parameter
+ }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
+
+ typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
+
+ // Non-member functions to give properties of the distribution.
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
+ { // Range of permissible values for random variable k.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(0, max_value<RealType>()); // Max integer?
+ }
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
+ { // Range of supported values for random variable k.
+ // 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>(0, max_value<RealType>());
+ }
+
+ template <class RealType, class Policy>
+ inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ { // Mean of poisson distribution = lambda.
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ return k + l;
+ } // mean
+
+ template <class RealType, class Policy>
+ inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ { // mode.
+ // TODO!
+ return 0;
+ }
+
+ template <class RealType, class Policy>
+ inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ { // variance.
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ return 2 * (2 * l + k);
+ }
+
+ // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ // standard_deviation provided by derived accessors.
+
+ template <class RealType, class Policy>
+ inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ { // skewness = sqrt(l).
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ BOOST_MATH_STD_USING
+ return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
+ }
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ {
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
+ } // kurtosis_excess
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
+ {
+ return kurtosis_excess(dist) + 3;
+ }
+
+ template <class RealType, class Policy>
+ RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
+ { // Probability Density/Mass Function.
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::pdf(%1%)";
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+
+ if(l == 0)
+ return pdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x);
+
+ r = -(x + l) / 2 + log(x / l) * (k / 4 - 0.5f);
+
+ return 0.5f * exp(r)
+ * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), Policy());
+ } // pdf
+
+ template <class RealType, class Policy>
+ RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
+ {
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+
+ return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
+ } // cdf
+
+ template <class RealType, class Policy>
+ RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
+ { // Complemented Cumulative Distribution Function
+ const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
+ non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
+ RealType x = c.param;
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+
+ return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
+ } // ccdf
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
+ { // Quantile (or Percent Point) function.
+ return detail::nccs_quantile(dist, p, false);
+ } // quantile
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
+ { // Quantile (or Percent Point) function.
+ return detail::nccs_quantile(c.dist, c.param, true);
+ } // quantile complement.
+
+ } // 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_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
+
+
+
Modified: sandbox/math_toolkit/boost/math/distributions/poisson.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/poisson.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/poisson.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -38,6 +38,7 @@
#include <boost/math/distributions/fwd.hpp>
#include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
+#include <boost/math/special_functions/trunc.hpp> // for incomplete gamma. gamma_q
#include <boost/math/distributions/complement.hpp> // complements
#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
#include <boost/math/special_functions/fpclassify.hpp> // isnan.
@@ -340,7 +341,7 @@
&& k < max_factorial<RealType>::value)
{ // k is small enough (for float 34, double 170 ...) to use factorial(k).
return exp(-mean) * pow(mean, k) /
- unchecked_factorial<RealType>(tools::real_cast<unsigned int>(floork));
+ unchecked_factorial<RealType>(itrunc(floork));
}
else
{ // Need to use log(factorial(k)) = lgamma(k+1)
Modified: sandbox/math_toolkit/boost/math/special_functions/bessel.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/bessel.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/bessel.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -24,6 +24,8 @@
#include <boost/math/special_functions/sin_pi.hpp>
#include <boost/math/special_functions/cos_pi.hpp>
#include <boost/math/special_functions/sinc.hpp>
+#include <boost/math/special_functions/trunc.hpp>
+#include <boost/math/special_functions/round.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/promotion.hpp>
@@ -127,7 +129,7 @@
if(floor(v) == v)
{
T r = cyl_bessel_j_imp(v, -x, t, pol);
- if(tools::real_cast<int>(v) & 1)
+ if(iround(v) & 1)
r = -r;
return r;
}
@@ -163,7 +165,7 @@
if(fabs(x) > asymptotic_bessel_j_limit<T>(v, tag_type()))
return asymptotic_bessel_j_large_x_2(v, x);
else
- return bessel_jn(tools::real_cast<int>(v), x, pol);
+ return bessel_jn(iround(v), x, pol);
}
return cyl_bessel_j_imp(v, x, bessel_no_int_tag(), pol);
}
@@ -226,7 +228,7 @@
if(floor(v) == v)
{
T r = cyl_bessel_i_imp(v, -x, pol);
- if(tools::real_cast<int>(v) & 1)
+ if(iround(v) & 1)
r = -r;
return r;
}
@@ -290,7 +292,7 @@
BOOST_MATH_STD_USING
if((floor(v) == v))
{
- return bessel_kn(tools::real_cast<int>(v), x, pol);
+ return bessel_kn(itrunc(v), x, pol);
}
return cyl_bessel_k_imp(v, x, bessel_no_int_tag(), pol);
}
@@ -335,12 +337,12 @@
if((fabs(x) > asymptotic_bessel_y_limit<T>(tag_type())) && (fabs(x) > 5 * abs(v)))
{
T r = asymptotic_bessel_y_large_x_2(static_cast<T>(abs(v)), x);
- if((v < 0) && (tools::real_cast<int>(v) & 1))
+ if((v < 0) && (itrunc(v) & 1))
r = -r;
return r;
}
else
- return bessel_yn(tools::real_cast<int>(v), x, pol);
+ return bessel_yn(itrunc(v), x, pol);
}
return cyl_neumann_imp<T>(v, x, bessel_no_int_tag(), pol);
}
@@ -358,7 +360,7 @@
return r;
}
else
- return bessel_yn(tools::real_cast<int>(v), x, pol);
+ return bessel_yn(v, x, pol);
}
template <class T, class Policy>
Modified: sandbox/math_toolkit/boost/math/special_functions/beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/beta.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -17,6 +17,7 @@
#include <boost/math/special_functions/erf.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/static_assert.hpp>
#include <cmath>
@@ -821,7 +822,7 @@
BOOST_MATH_STD_USING // ADL of std names
T result = pow(x, n);
T term = result;
- for(unsigned i = tools::real_cast<unsigned>(n - 1); i > k; --i)
+ for(unsigned i = itrunc(n - 1); i > k; --i)
{
term *= ((i + 1) * y) / ((n - i) * x) ;
result += term;
@@ -1059,7 +1060,7 @@
else if(a > 15)
{
// sidestep so we can use the series representation:
- int n = static_cast<int>(boost::math::tools::real_cast<long double>(floor(b)));
+ int n = itrunc(floor(b));
if(n == b)
--n;
T bbar = b - n;
@@ -1081,7 +1082,7 @@
// the formula here for the non-normalised case is tricky to figure
// out (for me!!), and requires two pochhammer calculations rather
// than one, so leave it for now....
- int n = static_cast<int>(boost::math::tools::real_cast<long double>(floor(b)));
+ int n = itrunc(floor(b));
T bbar = b - n;
if(bbar <= 0)
{
Modified: sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/cos_pi.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -12,7 +12,7 @@
#include <cmath>
#include <boost/math/tools/config.hpp>
-#include <boost/math/tools/real_cast.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/constants/constants.hpp>
namespace boost{ namespace math{
@@ -31,7 +31,7 @@
}
T rem = floor(x);
- if(tools::real_cast<int>(rem) & 1)
+ if(itrunc(rem) & 1)
invert = !invert;
rem = x - rem;
if(rem > 0.5f)
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_ik.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -234,7 +234,7 @@
v = -v; // v is non-negative from here
kind |= need_k;
}
- n = tools::real_cast<unsigned>(v + 0.5f);
+ n = iround(v);
u = v - n; // -1/2 <= u < 1/2
BOOST_MATH_INSTRUMENT_VARIABLE(n);
BOOST_MATH_INSTRUMENT_VARIABLE(u);
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/bessel_jy.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -225,7 +225,7 @@
v = -v; // v is non-negative from here
kind = need_j|need_y; // need both for reflection formula
}
- n = real_cast<unsigned>(v + 0.5L);
+ n = iround(v);
u = v - n; // -1/2 <= u < 1/2
if (x == 0)
Modified: sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/detail/t_distribution_inv.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -12,6 +12,7 @@
#endif
#include <boost/math/special_functions/cbrt.hpp>
+#include <boost/math/special_functions/round.hpp>
namespace boost{ namespace math{ namespace detail{
@@ -204,7 +205,7 @@
//
T tolerance = ldexp(1.0f, (2 * policies::digits<T, Policy>()) / 3);
- switch(boost::math::tools::real_cast<int>(df))
+ switch(itrunc(df))
{
case 1:
{
@@ -368,7 +369,7 @@
// where we use Shaw's tail series.
// The crossover point is roughly exponential in -df:
//
- T crossover = ldexp(1.0f, tools::real_cast<int>(df / -0.654f));
+ T crossover = ldexp(1.0f, iround(df / -0.654f));
if(u > crossover)
{
result = boost::math::detail::inverse_students_t_hill(df, u, pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/factorials.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/factorials.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/factorials.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -163,7 +163,7 @@
// handle it, split the product up into three parts:
//
T xp1 = x + 1;
- unsigned n2 = tools::real_cast<unsigned>(floor(xp1));
+ unsigned n2 = itrunc(floor(xp1));
if(n2 == xp1)
return 0;
T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
Modified: sandbox/math_toolkit/boost/math/special_functions/gamma.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/gamma.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/gamma.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -27,12 +27,12 @@
#include <boost/math/tools/series.hpp>
#include <boost/math/tools/fraction.hpp>
#include <boost/math/tools/precision.hpp>
-#include <boost/math/tools/real_cast.hpp>
#include <boost/math/tools/promotion.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/special_functions/log1p.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/special_functions/powm1.hpp>
#include <boost/math/special_functions/sqrt1pm1.hpp>
#include <boost/math/special_functions/lanczos.hpp>
@@ -165,7 +165,7 @@
}
if((floor(z) == z) && (z < max_factorial<T>::value))
{
- result *= unchecked_factorial<T>(tools::real_cast<unsigned>(z) - 1);
+ result *= unchecked_factorial<T>(itrunc(z) - 1);
}
else
{
@@ -370,7 +370,7 @@
BOOST_MATH_INSTRUMENT_CODE(prefix);
if((floor(z) == z) && (z < max_factorial<T>::value))
{
- prefix *= unchecked_factorial<T>(tools::real_cast<unsigned>(z) - 1);
+ prefix *= unchecked_factorial<T>(itrunc(z) - 1);
}
else
{
@@ -1095,7 +1095,7 @@
//
if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
{
- return unchecked_factorial<T>(tools::real_cast<unsigned>(z) - 1) / unchecked_factorial<T>(tools::real_cast<unsigned>(z + delta) - 1);
+ return unchecked_factorial<T>((unsigned)itrunc(z) - 1) / unchecked_factorial<T>((unsigned)itrunc(z + delta) - 1);
}
}
if(fabs(delta) < 20)
Modified: sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp (original)
+++ sandbox/math_toolkit/boost/math/special_functions/sin_pi.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -12,7 +12,7 @@
#include <cmath>
#include <boost/math/tools/config.hpp>
-#include <boost/math/tools/real_cast.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/constants/constants.hpp>
namespace boost{ namespace math{
@@ -36,7 +36,7 @@
invert = false;
T rem = floor(x);
- if(tools::real_cast<int>(rem) & 1)
+ if(itrunc(rem) & 1)
invert = !invert;
rem = x - rem;
if(rem > 0.5f)
Modified: sandbox/math_toolkit/boost/math/tools/polynomial.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/polynomial.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/polynomial.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -13,6 +13,7 @@
#include <boost/assert.hpp>
#include <boost/math/tools/rational.hpp>
#include <boost/math/tools/real_cast.hpp>
+#include <boost/math/special_functions/binomial.hpp>
#include <vector>
@@ -46,8 +47,8 @@
Seq polynomial_to_chebyshev(const Seq& s)
{
// Converts a Polynomial into Chebyshev form:
- typedef Seq::value_type value_type;
- typedef Seq::difference_type difference_type;
+ typedef typename Seq::value_type value_type;
+ typedef typename Seq::difference_type difference_type;
Seq result(s);
difference_type order = s.size() - 1;
difference_type even_order = order & 1 ? order - 1 : order;
Modified: sandbox/math_toolkit/libs/math/performance/distributions.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/performance/distributions.cpp (original)
+++ sandbox/math_toolkit/libs/math/performance/distributions.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -7,19 +7,7 @@
#include "performance_measure.hpp"
-#include <boost/math/distributions/beta.hpp>
-#include <boost/math/distributions/binomial.hpp>
-#include <boost/math/distributions/cauchy.hpp>
-#include <boost/math/distributions/chi_squared.hpp>
-#include <boost/math/distributions/exponential.hpp>
-#include <boost/math/distributions/fisher_f.hpp>
-#include <boost/math/distributions/gamma.hpp>
-#include <boost/math/distributions/lognormal.hpp>
-#include <boost/math/distributions/negative_binomial.hpp>
-#include <boost/math/distributions/normal.hpp>
-#include <boost/math/distributions/poisson.hpp>
-#include <boost/math/distributions/students_t.hpp>
-#include <boost/math/distributions/weibull.hpp>
+#include <boost/math/distributions.hpp>
double probabilities[] = {
1e-5,
@@ -201,6 +189,7 @@
BOOST_MATH_DISTRIBUTION1_TEST(poisson, real_values, int_values, probabilities)
BOOST_MATH_DISTRIBUTION1_TEST(students_t, int_values, real_values, probabilities)
BOOST_MATH_DISTRIBUTION2_TEST(weibull, real_values, real_values, real_values, probabilities)
+BOOST_MATH_DISTRIBUTION2_TEST(non_central_chi_squared, int_values, int_values, real_values, probabilities)
#ifdef TEST_R
@@ -341,6 +330,7 @@
BOOST_MATH_R_DISTRIBUTION1_TEST(pois, real_values, int_values, probabilities)
BOOST_MATH_R_DISTRIBUTION1_TEST(t, int_values, real_values, probabilities)
BOOST_MATH_R_DISTRIBUTION2_TEST(weibull, real_values, real_values, real_values, probabilities)
+BOOST_MATH_R_DISTRIBUTION2_TEST(nchisq, int_values, int_values, real_values, probabilities)
#endif
Modified: sandbox/math_toolkit/libs/math/test/functor.hpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/functor.hpp (original)
+++ sandbox/math_toolkit/libs/math/test/functor.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -3,6 +3,8 @@
// Boost Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+#include <boost/math/special_functions/trunc.hpp>
+
#ifndef BOOST_MATH_TEST_FUNCTOR_HPP
#define BOOST_MATH_TEST_FUNCTOR_HPP
@@ -121,7 +123,7 @@
template <class S>
typename S::value_type operator()(const S& row)
{
- return f(boost::math::tools::real_cast<int>(row[m_i]), row[m_j]);
+ return f(boost::math::itrunc(row[m_i]), row[m_j]);
}
private:
@@ -144,8 +146,8 @@
typename S::value_type operator()(const S& row)
{
return f(
- boost::math::tools::real_cast<int>(row[m_i]),
- boost::math::tools::real_cast<int>(row[m_j]),
+ boost::math::itrunc(row[m_i]),
+ boost::math::itrunc(row[m_j]),
row[m_k]);
}
@@ -169,8 +171,8 @@
typename S::value_type operator()(const S& row)
{
return f(
- boost::math::tools::real_cast<int>(row[m_i]),
- boost::math::tools::real_cast<int>(row[m_j]),
+ boost::math::itrunc(row[m_i]),
+ boost::math::itrunc(row[m_j]),
row[m_k],
row[m_l]);
}
Modified: sandbox/math_toolkit/libs/math/test/test_bessel_i.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_bessel_i.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_bessel_i.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -8,6 +8,7 @@
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/bessel.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/type_traits/is_floating_point.hpp>
#include <boost/array.hpp>
#include "functor.hpp"
@@ -103,7 +104,7 @@
{
return static_cast<T>(
boost::math::cyl_bessel_i(
- boost::math::tools::real_cast<int>(v), x));
+ boost::math::itrunc(v), x));
}
template <class T>
Modified: sandbox/math_toolkit/libs/math/test/test_bessel_j.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_bessel_j.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_bessel_j.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -330,7 +330,7 @@
template <class T>
T cyl_bessel_j_int_wrapper(T v, T x)
{
- return static_cast<T>(boost::math::cyl_bessel_j(boost::math::tools::real_cast<int>(v), x));
+ return static_cast<T>(boost::math::cyl_bessel_j(boost::math::itrunc(v), x));
}
Modified: sandbox/math_toolkit/libs/math/test/test_bessel_k.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_bessel_k.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_bessel_k.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -100,7 +100,7 @@
{
return static_cast<T>(
boost::math::cyl_bessel_k(
- boost::math::tools::real_cast<int>(v), x));
+ boost::math::itrunc(v), x));
}
template <class T>
Modified: sandbox/math_toolkit/libs/math/test/test_bessel_y.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_bessel_y.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_bessel_y.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -289,7 +289,7 @@
template <class T>
T cyl_neumann_int_wrapper(T v, T x)
{
- return static_cast<T>(boost::math::cyl_neumann(boost::math::tools::real_cast<int>(v), x));
+ return static_cast<T>(boost::math::cyl_neumann(boost::math::itrunc(v), x));
}
template <class T>
Modified: sandbox/math_toolkit/libs/math/test/test_binomial_coeff.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_binomial_coeff.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_binomial_coeff.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -7,6 +7,7 @@
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/special_functions/binomial.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/math/tools/test.hpp>
#include "functor.hpp"
#include <boost/array.hpp>
@@ -94,8 +95,8 @@
T binomial_wrapper(T n, T k)
{
return boost::math::binomial_coefficient<T>(
- boost::math::tools::real_cast<unsigned>(n),
- boost::math::tools::real_cast<unsigned>(k));
+ boost::math::itrunc(n),
+ boost::math::itrunc(k));
}
template <class T>
Modified: sandbox/math_toolkit/libs/math/test/test_expint.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_expint.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_expint.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -5,6 +5,7 @@
#include <boost/math/concepts/real_concept.hpp>
#include <boost/math/special_functions/expint.hpp>
+#include <boost/math/special_functions/trunc.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <boost/math/tools/stats.hpp>
@@ -106,7 +107,7 @@
T expint_wrapper(T n, T z)
{
return boost::math::expint(
- boost::math::tools::real_cast<unsigned>(n), z);
+ boost::math::itrunc(n), z);
}
#ifdef TEST_OTHER
@@ -114,7 +115,7 @@
T other_expint_wrapper(T n, T z)
{
return other::expint(
- boost::math::tools::real_cast<unsigned>(n), z);
+ boost::math::itrunc(n), z);
}
#endif
template <class T>
Modified: sandbox/math_toolkit/libs/math/test/test_gamma_dist.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_gamma_dist.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_gamma_dist.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -94,7 +94,7 @@
// The first tests use values generated by MathCAD,
// and should be accurate to around double precision.
//
- RealType tolerance = (std::max)(5e-14f, boost::math::tools::real_cast<float>(std::numeric_limits<RealType>::epsilon() * 20)) * 100;
+ RealType tolerance = (std::max)(RealType(5e-14f), std::numeric_limits<RealType>::epsilon() * 20) * 100;
cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl;
check_gamma(
Added: sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,379 @@
+// test_nc_chi_squared.cpp
+
+// Copyright John Maddock 2008.
+
+// 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)
+
+#include <boost/math/concepts/real_concept.hpp> // for real_concept
+#include <boost/math/distributions/non_central_chi_squared.hpp> // for chi_squared_distribution
+#include <boost/test/included/test_exec_monitor.hpp> // for test_main
+#include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
+
+#include "functor.hpp"
+#include "handle_test_result.hpp"
+#include "test_nccs_hooks.hpp"
+
+#include <iostream>
+using std::cout;
+using std::endl;
+#include <limits>
+using std::numeric_limits;
+
+#define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
+ {\
+ unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
+ BOOST_CHECK_CLOSE(a, b, prec); \
+ if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
+ {\
+ std::cerr << "Failure was at row " << i << std::endl;\
+ std::cerr << std::setprecision(35); \
+ std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
+ std::cerr << " , " << data[i][3] << " , " << data[i][4] << " , " << data[i][5] << " } " << std::endl;\
+ }\
+ }
+
+template <class RealType>
+RealType naive_pdf(RealType v, RealType lam, RealType x)
+{
+ // Formula direct from
+ // http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
+ // with no simplification:
+ RealType r = -(x+lam)/2 + log(x) * (v-1)/2 + log(lam) / 2;
+ r -= log(lam * x) * v/4;
+ r = exp(r) / 2;
+ r *= boost::math::cyl_bessel_i(v/2 - 1, sqrt(lam * x));
+ return r;
+}
+
+template <class RealType>
+void test_spot(
+ RealType df, // Degrees of freedom
+ RealType ncp, // non-centrality param
+ RealType cs, // Chi Square statistic
+ RealType P, // CDF
+ RealType Q, // Complement of CDF
+ RealType tol) // Test tolerance
+{
+ boost::math::non_central_chi_squared_distribution<RealType> dist(df, ncp);
+ BOOST_CHECK_CLOSE(
+ cdf(dist, cs), P, tol);
+ BOOST_CHECK_CLOSE(
+ pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), ncp, cs), tol);
+ if((P < 0.99) && (Q < 0.99))
+ {
+ //
+ // We can only check this if P is not too close to 1,
+ // so that we can guarentee Q is free of error:
+ //
+ BOOST_CHECK_CLOSE(
+ cdf(complement(dist, cs)), Q, tol);
+ BOOST_CHECK_CLOSE(
+ quantile(dist, P), cs, tol);
+ BOOST_CHECK_CLOSE(
+ quantile(complement(dist, Q)), cs, tol);
+ }
+}
+
+template <class RealType> // Any floating-point type RealType.
+void test_spots(RealType)
+{
+ RealType tolerance = (std::max)(
+ boost::math::tools::epsilon<RealType>(),
+ (RealType)boost::math::tools::epsilon<double>() * 5) * 100;
+ //
+ // At float precision we need to up the tolerance, since
+ // the input values are rounded off to inexact quantities
+ // the results get thrown off by a noticeable amount.
+ //
+ if(boost::math::tools::digits<RealType>() < 50)
+ tolerance *= 50;
+
+ cout << "Tolerance = " << tolerance << "%." << endl;
+
+ using boost::math::chi_squared_distribution;
+ using ::boost::math::chi_squared;
+ using ::boost::math::cdf;
+ using ::boost::math::pdf;
+ //
+ // Test against the data from Table 6 of:
+ //
+ // "Self-Validating Computations of Probabilities for Selected
+ // Central and Noncentral Univariate Probability Functions."
+ // Morgan C. Wang; William J. Kennedy
+ // Journal of the American Statistical Association,
+ // Vol. 89, No. 427. (Sep., 1994), pp. 878-887.
+ //
+ test_spot(
+ static_cast<RealType>(1), // degrees of freedom
+ static_cast<RealType>(6), // non centrality
+ static_cast<RealType>(0.00393), // Chi Squared statistic
+ static_cast<RealType>(0.2498463724258039e-2), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.2498463724258039e-2), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(5), // degrees of freedom
+ static_cast<RealType>(1), // non centrality
+ static_cast<RealType>(9.23636), // Chi Squared statistic
+ static_cast<RealType>(0.8272918751175548), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.8272918751175548), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(11), // degrees of freedom
+ static_cast<RealType>(21), // non centrality
+ static_cast<RealType>(24.72497), // Chi Squared statistic
+ static_cast<RealType>(0.2539481822183126), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.2539481822183126), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(31), // degrees of freedom
+ static_cast<RealType>(6), // non centrality
+ static_cast<RealType>(44.98534), // Chi Squared statistic
+ static_cast<RealType>(0.8125198785064969), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.8125198785064969), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(51), // degrees of freedom
+ static_cast<RealType>(1), // non centrality
+ static_cast<RealType>(38.56038), // Chi Squared statistic
+ static_cast<RealType>(0.8519497361859118e-1), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.8519497361859118e-1), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(100), // degrees of freedom
+ static_cast<RealType>(16), // non centrality
+ static_cast<RealType>(82.35814), // Chi Squared statistic
+ static_cast<RealType>(0.1184348822747824e-1), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.1184348822747824e-1), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(300), // degrees of freedom
+ static_cast<RealType>(16), // non centrality
+ static_cast<RealType>(331.78852), // Chi Squared statistic
+ static_cast<RealType>(0.7355956710306709), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.7355956710306709), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(500), // degrees of freedom
+ static_cast<RealType>(21), // non centrality
+ static_cast<RealType>(459.92612), // Chi Squared statistic
+ static_cast<RealType>(0.2797023600800060e-1), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.2797023600800060e-1), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(1), // degrees of freedom
+ static_cast<RealType>(1), // non centrality
+ static_cast<RealType>(0.00016), // Chi Squared statistic
+ static_cast<RealType>(0.6121428929881423e-2), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.6121428929881423e-2), // Q = 1 - P
+ tolerance);
+ test_spot(
+ static_cast<RealType>(1), // degrees of freedom
+ static_cast<RealType>(1), // non centrality
+ static_cast<RealType>(0.00393), // Chi Squared statistic
+ static_cast<RealType>(0.3033814229753780e-1), // Probability of result (CDF), P
+ static_cast<RealType>(1-0.3033814229753780e-1), // Q = 1 - P
+ tolerance);
+
+ RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage
+ boost::math::non_central_chi_squared_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12));
+ RealType x = 7;
+ using namespace std; // ADL of std names.
+ // mean:
+ BOOST_CHECK_CLOSE(
+ mean(dist)
+ , static_cast<RealType>(8+12), tol2);
+ // variance:
+ BOOST_CHECK_CLOSE(
+ variance(dist)
+ , static_cast<RealType>(64), tol2);
+ // std deviation:
+ BOOST_CHECK_CLOSE(
+ standard_deviation(dist)
+ , static_cast<RealType>(8), tol2);
+ // hazard:
+ BOOST_CHECK_CLOSE(
+ hazard(dist, x)
+ , pdf(dist, x) / cdf(complement(dist, x)), tol2);
+ // cumulative hazard:
+ BOOST_CHECK_CLOSE(
+ chf(dist, x)
+ , -log(cdf(complement(dist, x))), tol2);
+ // coefficient_of_variation:
+ BOOST_CHECK_CLOSE(
+ coefficient_of_variation(dist)
+ , standard_deviation(dist) / mean(dist), tol2);
+#if 0
+ // mode:
+ BOOST_CHECK_CLOSE(
+ mode(dist)
+ , static_cast<RealType>(6), tol2);
+#endif
+ BOOST_CHECK_CLOSE(
+ median(dist),
+ quantile(
+ non_central_chi_squared_distribution<RealType>(
+ static_cast<RealType>(8),
+ static_cast<RealType>(12)),
+ static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
+ // skewness:
+ BOOST_CHECK_CLOSE(
+ skewness(dist)
+ , static_cast<RealType>(0.6875), tol2);
+ // kurtosis:
+ BOOST_CHECK_CLOSE(
+ kurtosis(dist)
+ , static_cast<RealType>(3.65625), tol2);
+ // kurtosis excess:
+ BOOST_CHECK_CLOSE(
+ kurtosis_excess(dist)
+ , static_cast<RealType>(0.65625), tol2);
+} // template <class RealType>void test_spots(RealType)
+
+template <class T>
+T nccs_cdf(T df, T nc, T x)
+{
+ return cdf(boost::math::non_central_chi_squared_distribution<T>(df, nc), x);
+}
+
+template <class T>
+T nccs_ccdf(T df, T nc, T x)
+{
+ return cdf(complement(boost::math::non_central_chi_squared_distribution<T>(df, nc), x));
+}
+
+template <typename T>
+void do_test_nc_chi_squared(T& data, const char* type_name, const char* test)
+{
+ typedef typename T::value_type row_type;
+ typedef typename row_type::value_type value_type;
+
+ std::cout << "Testing: " << test << std::endl;
+
+ value_type (*fp1)(value_type, value_type, value_type) = nccs_cdf;
+ boost::math::tools::test_result<value_type> result;
+
+ result = boost::math::tools::test(
+ data,
+ bind_func(fp1, 0, 1, 2),
+ extract_result(3));
+ handle_test_result(result, data[result.worst()], result.worst(),
+ type_name, "CDF", test);
+
+ fp1 = nccs_ccdf;
+ result = boost::math::tools::test(
+ data,
+ bind_func(fp1, 0, 1, 2),
+ extract_result(4));
+ handle_test_result(result, data[result.worst()], result.worst(),
+ type_name, "CCDF", test);
+
+#ifdef TEST_OTHER
+ fp1 = other::nccs_cdf;
+ result = boost::math::tools::test(
+ data,
+ bind_func(fp1, 0, 1, 2),
+ extract_result(3));
+ handle_test_result(result, data[result.worst()], result.worst(),
+ type_name, "other::CDF", test);
+#endif
+
+ std::cout << std::endl;
+
+}
+
+template <typename T>
+void quantile_sanity_check(T& data, const char* type_name, const char* test)
+{
+ typedef typename T::value_type row_type;
+ typedef typename row_type::value_type value_type;
+
+ std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
+
+ //
+ // These sanity checks test for a round trip accuracy of one half
+ // of the bits in T, unless T is type float, in which case we check
+ // for just one decimal digit. The problem here is the sensitivity
+ // of the functions, not their accuracy. This test data was generated
+ // for the forward functions, which means that when it is used as
+ // the input to the inverses then it is necessarily inexact. This rounding
+ // of the input is what makes the data unsuitable for use as an accuracy check,
+ // and also demonstrates that you can't in general round-trip these functions.
+ // It is however a useful sanity check.
+ //
+ value_type precision = static_cast<value_type>(ldexp(1.0, 1-boost::math::policies::digits<value_type, boost::math::policies::policy<> >()/2)) * 100;
+ if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
+ precision = 1; // 1% or two decimal digits, all we can hope for when the input is truncated to float
+
+ for(unsigned i = 0; i < data.size(); ++i)
+ {
+ if(data[i][3] == 0)
+ {
+ BOOST_CHECK(0 == quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]));
+ }
+ else if(data[i][3] < 0.9999f)
+ {
+ value_type p = quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]);
+ value_type pt = data[i][2];
+ BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
+ }
+ if(data[i][4] == 0)
+ {
+ BOOST_CHECK(0 == quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3])));
+ }
+ else if(data[i][4] < 0.9999f)
+ {
+ value_type p = quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][4]));
+ value_type pt = data[i][2];
+ BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
+ }
+ }
+}
+
+template <typename T>
+void test_accuracy(T, const char* type_name)
+{
+#if 0
+#include "nccs.ipp"
+ do_test_nc_chi_squared(nccs, type_name, "Non Central Chi Squared, medium parameters");
+ quantile_sanity_check(nccs, type_name, "Non Central Chi Squared, medium parameters");
+#endif
+#include "nccs_big.ipp"
+ do_test_nc_chi_squared(nccs_big, type_name, "Non Central Chi Squared, large parameters");
+ quantile_sanity_check(nccs_big, type_name, "Non Central Chi Squared, large parameters");
+}
+
+int test_main(int, char* [])
+{
+ BOOST_MATH_CONTROL_FP;
+ // Basic sanity-check spot values.
+
+ // (Parameter value, arbitrarily zero, only communicates the floating point type).
+ test_spots(0.0F); // Test float.
+ test_spots(0.0); // Test double.
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ test_spots(0.0L); // Test long double.
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
+#endif
+#endif
+ test_accuracy(0.0F, "float"); // Test float.
+ test_accuracy(0.0, "double"); // Test double.
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ test_accuracy(0.0L, "long double"); // Test long double.
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ test_accuracy(boost::math::concepts::real_concept(0.), "real_concept"); // Test real concept.
+#endif
+#endif
+ return 0;
+} // int test_main(int, char* [])
+
+/*
+
+
+*/
+
+
+
Added: sandbox/math_toolkit/libs/math/test/test_nccs_hooks.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/test_nccs_hooks.hpp 2008-01-14 04:27:44 EST (Mon, 14 Jan 2008)
@@ -0,0 +1,61 @@
+// (C) Copyright John Maddock 2008.
+// 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_TEST_NCCS_OTHER_HOOKS_HPP
+#define BOOST_MATH_TEST_NCCS_OTHER_HOOKS_HPP
+
+
+#ifdef TEST_R
+#define MATHLIB_STANDALONE
+#include <rmath.h>
+namespace other{
+inline float nccs_cdf(float df, float nc, float x)
+{
+ return (float)pnchisq(x, df, nc, 1, 0);
+}
+inline double nccs_cdf(double df, double nc, double x)
+{
+ return pnchisq(x, df, nc, 1, 0);
+}
+inline long double nccs_cdf(long double df, long double nc, long double x)
+{
+ return pnchisq((double)x, (double)df, (double)nc, 1, 0);
+}
+}
+#define TEST_OTHER
+#endif
+
+#ifdef TEST_CDFLIB
+#include <cdflib.h>
+namespace other{
+inline double nccs_cdf(double df, double nc, double x)
+{
+ int kind(1), status(0);
+ double p, q, bound(0);
+ cdfchn(&kind, &p, &q, &x, &df, &nc, &status, &bound);
+ return p;
+}
+inline float nccs_cdf(float df, float nc, float x)
+{
+ return (double)nccs_cdf((double)df, (double)nc, (double)x);
+}
+inline long double nccs_cdf(long double df, long double nc, long double x)
+{
+ return nccs_cdf((double)df, (double)nc, (double)x);
+}
+}
+#define TEST_OTHER
+#endif
+
+#ifdef TEST_OTHER
+namespace other{
+ boost::math::concepts::real_concept nccs_cdf(boost::math::concepts::real_concept, boost::math::concepts::real_concept, boost::math::concepts::real_concept){ return 0; }
+}
+#endif
+
+
+#endif
+
+
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk