Added gamma distribution and tests.

Fixed a lot of gcc-linux issues: mostly missing "using namespace std" declarations. Added new concept check to test for these kinds of failures at compile time. [SVN r3344]
2026-01-19 04:22:09 +00:00 · 2006-11-03 18:42:23 +00:00
parent 041d476768
commit feacc055d3
19 changed files with 1238 additions and 23 deletions
--- a/include/boost/math/concepts/distributions.hpp
+++ b/include/boost/math/concepts/distributions.hpp
@@ -97,7 +97,7 @@ struct DistributionConcept
      function_requires<CopyConstructibleConcept<Distribution> >();
      function_requires<AssignableConcept<Distribution> >();

-      typedef Distribution::value_type value_type;
+      typedef typename Distribution::value_type value_type;

      const Distribution& dist = DistributionConcept<Distribution>::get_object();

--- a/include/boost/math/concepts/std_real_concept.hpp
+++ b/include/boost/math/concepts/std_real_concept.hpp
@@ -0,0 +1,344 @@
+//  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)
+
+// std_real_concept is an archetype for built-in Real types.
+
+// The main purpose in providing this type, is to verify
+// that std lib functions are found via a using declaration
+// bringing those functions into the current scope, and not
+// just because they happen to be in global scope.
+//
+// If ::pow is found rather than std::pow say, then the code
+// will silently compile, but trunctaion of long doubles to
+// double will cause a significant loss of precision.
+// A template instantiated with std_real_concept will *only*
+// compile if it std::whatever is in scope.
+
+#include <boost/config.hpp>
+#include <boost/limits.hpp>
+#include <boost/math/tools/real_cast.hpp>
+#include <boost/math/tools/precision.hpp>
+
+#include <ostream>
+#include <istream>
+#include <cmath>
+
+#ifndef BOOST_MATH_STD_REAL_CONCEPT_HPP
+#define BOOST_MATH_STD_REAL_CONCEPT_HPP
+
+namespace boost{ namespace math{
+
+namespace concepts
+{
+
+class std_real_concept
+{
+public:
+   // Constructors:
+   std_real_concept() : m_value(0){}
+   std_real_concept(char c) : m_value(c){}
+#ifndef BOOST_NO_INTRINSIC_WCHAR_T
+   std_real_concept(wchar_t c) : m_value(c){}
+#endif
+   std_real_concept(unsigned char c) : m_value(c){}
+   std_real_concept(signed char c) : m_value(c){}
+   std_real_concept(unsigned short c) : m_value(c){}
+   std_real_concept(short c) : m_value(c){}
+   std_real_concept(unsigned int c) : m_value(c){}
+   std_real_concept(int c) : m_value(c){}
+   std_real_concept(unsigned long c) : m_value(c){}
+   std_real_concept(long c) : m_value(c){}
+#ifdef BOOST_HAS_LONG_LONG
+   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)){}
+#endif
+   std_real_concept(float c) : m_value(c){}
+   std_real_concept(double c) : m_value(c){}
+   std_real_concept(long double c) : m_value(c){}
+
+   // Assignment:
+   std_real_concept& operator=(char c) { m_value = c; return *this; }
+   std_real_concept& operator=(unsigned char c) { m_value = c; return *this; }
+   std_real_concept& operator=(signed char c) { m_value = c; return *this; }
+#ifndef BOOST_NO_INTRINSIC_WCHAR_T
+   std_real_concept& operator=(wchar_t c) { m_value = c; return *this; }
+#endif
+   std_real_concept& operator=(short c) { m_value = c; return *this; }
+   std_real_concept& operator=(unsigned short c) { m_value = c; return *this; }
+   std_real_concept& operator=(int c) { m_value = c; return *this; }
+   std_real_concept& operator=(unsigned int c) { m_value = c; return *this; }
+   std_real_concept& operator=(long c) { m_value = c; return *this; }
+   std_real_concept& operator=(unsigned long c) { m_value = c; return *this; }
+#ifdef BOOST_HAS_LONG_LONG
+   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; }
+#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; }
+
+   // Member arithmetic:
+   std_real_concept& operator+=(const std_real_concept& other)
+   { m_value += other.value(); return *this; }
+   std_real_concept& operator-=(const std_real_concept& other)
+   { m_value -= other.value(); return *this; }
+   std_real_concept& operator*=(const std_real_concept& other)
+   { m_value *= other.value(); return *this; }
+   std_real_concept& operator/=(const std_real_concept& other)
+   { m_value /= other.value(); return *this; }
+   std_real_concept operator-()const
+   { return -m_value; }
+   std_real_concept const& operator+()const
+   { return *this; }
+
+private:
+   long double m_value;
+};
+
+// Non-member arithmetic:
+inline std_real_concept operator+(const std_real_concept& a, const std_real_concept& b)
+{
+   std_real_concept result(a);
+   result += b;
+   return result;
+}
+inline std_real_concept operator-(const std_real_concept& a, const std_real_concept& b)
+{
+   std_real_concept result(a);
+   result -= b;
+   return result;
+}
+inline std_real_concept operator*(const std_real_concept& a, const std_real_concept& b)
+{
+   std_real_concept result(a);
+   result *= b;
+   return result;
+}
+inline std_real_concept operator/(const std_real_concept& a, const std_real_concept& b)
+{
+   std_real_concept result(a);
+   result /= b;
+   return result;
+}
+
+// Comparison:
+inline bool operator == (const std_real_concept& a, const std_real_concept& b)
+{ return a.value() == b.value(); }
+inline bool operator != (const std_real_concept& a, const std_real_concept& b)
+{ return a.value() != b.value();}
+inline bool operator < (const std_real_concept& a, const std_real_concept& b)
+{ return a.value() < b.value(); }
+inline bool operator <= (const std_real_concept& a, const std_real_concept& b)
+{ return a.value() <= b.value(); }
+inline bool operator > (const std_real_concept& a, const std_real_concept& b)
+{ return a.value() > b.value(); }
+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
+
+namespace std{
+
+// Non-member functions:
+inline boost::math::concepts::std_real_concept acos(boost::math::concepts::std_real_concept a)
+{ return std::acos(a.value()); }
+inline boost::math::concepts::std_real_concept cos(boost::math::concepts::std_real_concept a)
+{ return std::cos(a.value()); }
+inline boost::math::concepts::std_real_concept asin(boost::math::concepts::std_real_concept a)
+{ return std::asin(a.value()); }
+inline boost::math::concepts::std_real_concept atan(boost::math::concepts::std_real_concept a)
+{ return std::atan(a.value()); }
+inline boost::math::concepts::std_real_concept atan2(boost::math::concepts::std_real_concept a, boost::math::concepts::std_real_concept b)
+{ return std::atan2(a.value(), b.value()); }
+inline boost::math::concepts::std_real_concept ceil(boost::math::concepts::std_real_concept a)
+{ return std::ceil(a.value()); }
+inline boost::math::concepts::std_real_concept fmod(boost::math::concepts::std_real_concept a, boost::math::concepts::std_real_concept b)
+{ return std::fmod(a.value(), b.value()); }
+inline boost::math::concepts::std_real_concept cosh(boost::math::concepts::std_real_concept a)
+{ return std::cosh(a.value()); }
+inline boost::math::concepts::std_real_concept exp(boost::math::concepts::std_real_concept a)
+{ return std::exp(a.value()); }
+inline boost::math::concepts::std_real_concept fabs(boost::math::concepts::std_real_concept a)
+{ return std::fabs(a.value()); }
+inline boost::math::concepts::std_real_concept abs(boost::math::concepts::std_real_concept a)
+{ return std::abs(a.value()); }
+inline boost::math::concepts::std_real_concept floor(boost::math::concepts::std_real_concept a)
+{ 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);
+   *ipart = ip;
+   return result;
+}
+inline boost::math::concepts::std_real_concept frexp(boost::math::concepts::std_real_concept a, int* expon)
+{ return std::frexp(a.value(), expon); }
+inline boost::math::concepts::std_real_concept ldexp(boost::math::concepts::std_real_concept a, int expon)
+{ return std::ldexp(a.value(), expon); }
+inline boost::math::concepts::std_real_concept log(boost::math::concepts::std_real_concept a)
+{ return std::log(a.value()); }
+inline boost::math::concepts::std_real_concept log10(boost::math::concepts::std_real_concept a)
+{ return std::log10(a.value()); }
+inline boost::math::concepts::std_real_concept tan(boost::math::concepts::std_real_concept a)
+{ return std::tan(a.value()); }
+inline boost::math::concepts::std_real_concept pow(boost::math::concepts::std_real_concept a, boost::math::concepts::std_real_concept b)
+{ return std::pow(a.value(), b.value()); }
+inline boost::math::concepts::std_real_concept pow(boost::math::concepts::std_real_concept a, int b)
+{ return std::pow(a.value(), b); }
+inline boost::math::concepts::std_real_concept sin(boost::math::concepts::std_real_concept a)
+{ return std::sin(a.value()); }
+inline boost::math::concepts::std_real_concept sinh(boost::math::concepts::std_real_concept a)
+{ return std::sinh(a.value()); }
+inline boost::math::concepts::std_real_concept sqrt(boost::math::concepts::std_real_concept a)
+{ return std::sqrt(a.value()); }
+inline boost::math::concepts::std_real_concept tanh(boost::math::concepts::std_real_concept a)
+{ return std::tanh(a.value()); }
+
+} // namespace std
+
+namespace boost{ namespace math{ namespace concepts{
+
+// Streaming:
+template <class charT, class traits>
+inline std::basic_ostream<charT, traits>& operator<<(std::basic_ostream<charT, traits>& os, const std_real_concept& a)
+{
+   return os << a.value();
+}
+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;
+   is >> v;
+   a = v;
+   return is;
+}
+
+} // namespace concepts
+
+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)
+{
+   return static_cast<unsigned int>(r.value());
+}
+
+template <>
+inline int real_cast<int, concepts::std_real_concept>(concepts::std_real_concept r)
+{
+   return static_cast<int>(r.value());
+}
+
+template <>
+inline long real_cast<long, concepts::std_real_concept>(concepts::std_real_concept r)
+{
+   return static_cast<long>(r.value());
+}
+
+// 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)
+{
+   return static_cast<float>(r.value());
+}
+template <>
+inline double real_cast<double, concepts::std_real_concept>(concepts::std_real_concept r)
+{
+   return static_cast<double>(r.value());
+}
+template <>
+inline long double real_cast<long double, concepts::std_real_concept>(concepts::std_real_concept r)
+{
+   return r.value();
+}
+
+template <>
+inline concepts::std_real_concept max_value<concepts::std_real_concept>(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+{
+   return max_value<long double>();
+}
+
+template <>
+inline concepts::std_real_concept min_value<concepts::std_real_concept>(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+{
+   return min_value<long double>();
+}
+
+template <>
+inline concepts::std_real_concept log_max_value<concepts::std_real_concept>(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+{
+   return log_max_value<long double>();
+}
+
+template <>
+inline concepts::std_real_concept log_min_value<concepts::std_real_concept>(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+{
+   return log_min_value<long double>();
+}
+
+template <>
+inline concepts::std_real_concept epsilon(BOOST_EXPLICIT_TEMPLATE_TYPE_SPEC(concepts::std_real_concept))
+{
+   return std::numeric_limits<long double>::epsilon();
+}
+
+template <>
+inline int digits<concepts::std_real_concept>(BOOST_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 std::numeric_limits<long double>::digits;
+}
+
+} // namespace tools
+
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_STD_REAL_CONCEPT_HPP
+
--- a/include/boost/math/distributions/extreme_value.hpp
+++ b/include/boost/math/distributions/extreme_value.hpp
@@ -157,7 +157,7 @@ RealType quantile(const complemented2_type<extreme_value_distribution<RealType>,
   if(q == 1)
      return -tools::overflow_error<RealType>(BOOST_CURRENT_FUNCTION, 0);

-   result = a - log(-log1p(-q)) * b;
+   result = a - log(-boost::math::log1p(-q)) * b;

   return result;
 }
--- a/include/boost/math/distributions/fisher_f.hpp
+++ b/include/boost/math/distributions/fisher_f.hpp
@@ -282,6 +282,7 @@ inline RealType mode(const fisher_f_distribution<RealType>& dist)
 template <class RealType>
 inline RealType skewness(const fisher_f_distribution<RealType>& dist)
 {
+   using namespace std; // ADL of std names
   // See http://mathworld.wolfram.com/F-Distribution.html
   RealType df1 = dist.degrees_of_freedom1();
   RealType df2 = dist.degrees_of_freedom2();
--- a/include/boost/math/distributions/gamma.hpp
+++ b/include/boost/math/distributions/gamma.hpp
@@ -0,0 +1,300 @@
+//  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_STATS_WEIBULL_HPP
+#define BOOST_STATS_WEIBULL_HPP
+
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm
+// http://mathworld.wolfram.com/WeibullDistribution.html
+
+#include <boost/math/special_functions/gamma.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/complement.hpp>
+
+namespace boost{ namespace math{
+
+namespace detail{
+
+template <class RealType>
+bool check_gamma_shape(
+      const char* function,
+      RealType shape,
+      RealType* result)
+{
+   if((shape <= 0) || !(boost::math::isfinite)(shape))
+   {
+      *result = tools::domain_error<RealType>(
+         function,
+         "Shape parameter is %1%, but must be > 0 !", shape);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType>
+bool check_gamma_x(
+      const char* function,
+      RealType const& x,
+      RealType* result)
+{
+   if((x < 0) || !(boost::math::isfinite)(x))
+   {
+      *result = tools::domain_error<RealType>(
+         function,
+         "Random variate is %1% but must be >= 0 !", x);
+      return false;
+   }
+   return true;
+}
+
+template <class RealType>
+inline bool check_gamma(
+      const char* function,
+      RealType scale,
+      RealType shape,
+      RealType* result)
+{
+   return check_scale(function, scale, result) && check_gamma_shape(function, shape, result);
+}
+
+} // namespace detail
+
+template <class RealType = double>
+class gamma_distribution
+{
+public:
+   typedef RealType value_type;
+
+   gamma_distribution(RealType shape, RealType scale = 1)
+      : m_shape(shape), m_scale(scale)
+   {
+      RealType result;
+      detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result);
+   }
+
+   RealType shape()const
+   {
+      return m_shape;
+   }
+
+   RealType scale()const
+   {
+      return m_scale;
+   }
+private:
+   //
+   // Data members:
+   //
+   RealType m_shape;     // distribution shape
+   RealType m_scale;     // distribution scale
+};
+
+template <class RealType>
+RealType pdf(const gamma_distribution<RealType>& dist, const RealType& x)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+   if(false == detail::check_gamma_x(BOOST_CURRENT_FUNCTION, x, &result))
+      return result;
+
+   if(x == 0)
+      return 0;
+
+   result = gamma_P_derivative(shape, x / scale) / scale;
+
+   return result;
+}
+
+template <class RealType>
+inline RealType cdf(const gamma_distribution<RealType>& dist, const RealType& x)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+   if(false == detail::check_gamma_x(BOOST_CURRENT_FUNCTION, x, &result))
+      return result;
+
+   result = boost::math::gamma_P(shape, x / scale);
+
+   return result;
+}
+
+template <class RealType>
+RealType quantile(const gamma_distribution<RealType>& dist, const RealType& p)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+   if(false == detail::check_probability(BOOST_CURRENT_FUNCTION, p, &result))
+      return result;
+
+   if(p == 1)
+      return tools::overflow_error<RealType>(BOOST_CURRENT_FUNCTION, 0);
+
+   result = gamma_P_inv(shape, p) * scale;
+
+   return result;
+}
+
+template <class RealType>
+RealType cdf(const complemented2_type<gamma_distribution<RealType>, RealType>& c)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = c.dist.shape();
+   RealType scale = c.dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+   if(false == detail::check_gamma_x(BOOST_CURRENT_FUNCTION, c.param, &result))
+      return result;
+
+   result = gamma_Q(shape, c.param / scale);
+
+   return result;
+}
+
+template <class RealType>
+RealType quantile(const complemented2_type<gamma_distribution<RealType>, RealType>& c)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = c.dist.shape();
+   RealType scale = c.dist.scale();
+   RealType q = c.param;
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+   if(false == detail::check_probability(BOOST_CURRENT_FUNCTION, q, &result))
+      return result;
+
+   if(q == 0)
+      return tools::overflow_error<RealType>(BOOST_CURRENT_FUNCTION, 0);
+
+   result = gamma_Q_inv(shape, q) * scale;
+
+   return result;
+}
+
+template <class RealType>
+inline RealType mean(const gamma_distribution<RealType>& dist)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+
+   result = shape * scale;
+   return result;
+}
+
+template <class RealType>
+RealType variance(const gamma_distribution<RealType>& dist)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+
+   result = shape * scale * scale;
+   return result;
+}
+
+template <class RealType>
+inline RealType mode(const gamma_distribution<RealType>& dist)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+
+   if(shape < 1)
+      return tools::domain_error<RealType>(
+         BOOST_CURRENT_FUNCTION,
+         "The mode of the gamma distribution is only defined for values of the shape parameter >= 1, but got %1%.",
+         shape);
+
+   result = (shape - 1) * scale;
+   return result;
+}
+
+template <class RealType>
+inline RealType skewness(const gamma_distribution<RealType>& dist)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+
+   result = 2 / sqrt(shape);
+   return result;
+}
+
+template <class RealType>
+inline RealType kurtosis_excess(const gamma_distribution<RealType>& dist)
+{
+   using namespace std;  // for ADL of std functions
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result;
+   if(false == detail::check_gamma(BOOST_CURRENT_FUNCTION, scale, shape, &result))
+      return result;
+
+   result = 6 / shape;
+   return result;
+}
+
+template <class RealType>
+inline RealType kurtosis(const gamma_distribution<RealType>& dist)
+{
+   return kurtosis_excess(dist) + 3;
+}
+
+} // namespace math
+} // namespace boost
+
+// This include must be at the end, *after* the accessors
+// for this distribution have been defined, in order to
+// keep compilers that support two-phase lookup happy.
+#include <boost/math/distributions/detail/derived_accessors.hpp>
+
+#endif // BOOST_STATS_WEIBULL_HPP
+
--- a/include/boost/math/distributions/weibull.hpp
+++ b/include/boost/math/distributions/weibull.hpp
@@ -153,7 +153,7 @@ RealType quantile(const weibull_distribution<RealType>& dist, const RealType& p)
   if(p == 1)
      return tools::overflow_error<RealType>(BOOST_CURRENT_FUNCTION, 0);

-   result = scale * pow(-log1p(-p), 1 / shape);
+   result = scale * pow(-boost::math::log1p(-p), 1 / shape);

   return result;
 }
--- a/include/boost/math/special_functions/beta.hpp
+++ b/include/boost/math/special_functions/beta.hpp
@@ -780,6 +780,7 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool
 template <class T>
 T binomial_ccdf(T n, T k, T x, T y)
 {
+   using namespace std; // ADL of std names
   T result = pow(x, n);
   T term = result;
   for(unsigned i = tools::real_cast<unsigned>(n - 1); i > k; --i)
--- a/include/boost/math/special_functions/detail/ibeta_inverse.hpp
+++ b/include/boost/math/special_functions/detail/ibeta_inverse.hpp
@@ -23,6 +23,8 @@ struct temme_root_finder

   std::tr1::tuple<T, T> operator()(T x)
   {
+      using namespace std; // ADL of std names
+
      T y = 1 - x;
      if(y == 0)
      {
@@ -51,6 +53,8 @@ private:
 template <class T>
 T temme_method_1_ibeta_inverse(T a, T b, T z)
 {
+   using namespace std; // ADL of std names
+
   const T r2 = sqrt(T(2));
   //
   // get the first approximation for eta from the inverse
@@ -127,6 +131,8 @@ T temme_method_1_ibeta_inverse(T a, T b, T z)
 template <class T>
 T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r, T theta)
 {
+   using namespace std; // ADL of std names
+
   //
   // Get first estimate for eta, see Eq 3.9 and 3.10,
   // but note there is a typo in Eq 3.10:
@@ -302,6 +308,8 @@ T temme_method_2_ibeta_inverse(T /*a*/, T /*b*/, T z, T r, T theta)
 template <class T>
 T temme_method_3_ibeta_inverse(T a, T b, T p, T q)
 {
+   using namespace std; // ADL of std names
+
   //
   // Begin by getting an initial approximation for the quantity
   // eta from the dominant part of the incomplete beta:
@@ -404,6 +412,8 @@ struct ibeta_roots

   std::tr1::tuple<T, T, T> operator()(T x)
   {
+      using namespace std; // ADL of std names
+
      BOOST_FPU_EXCEPTION_GUARD
      T f = ibeta_imp(a, b, x, L(), invert, true) - target;
      T f1 = invert ?
--- a/include/boost/math/special_functions/detail/igamma_large.hpp
+++ b/include/boost/math/special_functions/detail/igamma_large.hpp
@@ -70,6 +70,7 @@ inline T igamma_temme_large(T, T, mpl::int_<0> const *)
 template <class T>
 T igamma_temme_large(T a, T x, mpl::int_<64> const *)
 {
+   using namespace std; // ADL of std functions
   T sigma = (x - a) / a;
   T phi = -log1pmx(sigma);
   T y = a * phi;
@@ -263,7 +264,7 @@ T igamma_temme_large(T a, T x, mpl::int_<64> const *)
   if(x < a)
      result = -result;

-   result += erfc(sqrt(y)) / 2;
+   result += boost::math::erfc(sqrt(y)) / 2;

   return result;
 }
@@ -274,6 +275,7 @@ T igamma_temme_large(T a, T x, mpl::int_<64> const *)
 template <class T>
 T igamma_temme_large(T a, T x, mpl::int_<53> const *)
 {
+   using namespace std; // ADL of std functions
   T sigma = (x - a) / a;
   T phi = -log1pmx(sigma);
   T y = a * phi;
@@ -404,7 +406,7 @@ T igamma_temme_large(T a, T x, mpl::int_<53> const *)
   if(x < a)
      result = -result;

-   result += erfc(sqrt(y)) / 2;
+   result += boost::math::erfc(sqrt(y)) / 2;

   return result;
 }
@@ -415,6 +417,7 @@ T igamma_temme_large(T a, T x, mpl::int_<53> const *)
 template <class T>
 T igamma_temme_large(T a, T x, mpl::int_<24> const *)
 {
+   using namespace std; // ADL of std functions
   T sigma = (x - a) / a;
   T phi = -log1pmx(sigma);
   T y = a * phi;
@@ -456,7 +459,7 @@ T igamma_temme_large(T a, T x, mpl::int_<24> const *)
   if(x < a)
      result = -result;

-   result += erfc(sqrt(y)) / 2;
+   result += boost::math::erfc(sqrt(y)) / 2;

   return result;
 }
@@ -470,6 +473,7 @@ T igamma_temme_large(T a, T x, mpl::int_<24> const *)
 template <class T>
 T igamma_temme_large(T a, T x, mpl::int_<113> const *)
 {
+   using namespace std; // ADL of std functions
   T sigma = (x - a) / a;
   T phi = -log1pmx(sigma);
   T y = a * phi;
@@ -752,7 +756,7 @@ T igamma_temme_large(T a, T x, mpl::int_<113> const *)
   if(x < a)
      result = -result;

-   result += erfc(sqrt(y)) / 2;
+   result += boost::math::erfc(sqrt(y)) / 2;

   return result;
 }
--- a/include/boost/math/special_functions/factorials.hpp
+++ b/include/boost/math/special_functions/factorials.hpp
@@ -126,6 +126,7 @@ T rising_factorial(T x, int n)
 template <class T>
 inline T falling_factorial(T x, unsigned n)
 {
+   using namespace std; // ADL of std names
   if(x == 0)
      return 0;
   if(x < 0)
--- a/include/boost/math/special_functions/gamma.hpp
+++ b/include/boost/math/special_functions/gamma.hpp
@@ -958,6 +958,7 @@ template <class T, class Tag>
 T tgammap1m1_imp(T dz, Tag const& tag,
                 const ::boost::math::lanczos::undefined_lanczos& l)
 {
+   using namespace std; // ADL of std names
   //
   // There should be a better solution than this, but the
   // algebra isn't easy for the general case....
@@ -1270,7 +1271,7 @@ T finite_half_gamma_Q(T a, T x)
   // Calculates normalised Q when a is a half-integer:
   //
   using namespace std;
-   T e = erfc(sqrt(x));
+   T e = boost::math::erfc(sqrt(x));
   if((e != 0) && (a > 1))
   {
      T term = exp(-x) / sqrt(constants::pi<T>() * x);
--- a/include/boost/math/special_functions/lanczos.hpp
+++ b/include/boost/math/special_functions/lanczos.hpp
@@ -1163,7 +1163,25 @@ template <class T>
 struct lanczos_traits
 {
   typedef T value_type;
-   typedef undefined_lanczos evaluation_type;
+   // typedef undefined_lanczos evaluation_type;
+   typedef typename mpl::if_c<
+      (std::numeric_limits<T>::is_specialized == 0)
+      || (std::numeric_limits<T>::digits > 113),
+      undefined_lanczos,
+      typename mpl::if_c<
+         std::numeric_limits<T>::digits <= 24,
+         lanczos6m24,
+         typename mpl::if_c<
+            std::numeric_limits<T>::digits <= 53,
+            lanczos13m53,
+            typename mpl::if_c<
+               std::numeric_limits<T>::digits <= 64,
+               lanczos17m64,
+               lanczos24m113
+            >::type
+         >::type
+      >::type
+   >::type evaluation_type;
 };

 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
--- a/include/boost/math/tools/fraction.hpp
+++ b/include/boost/math/tools/fraction.hpp
@@ -66,11 +66,6 @@ namespace detail

 } // namespace detail

-// Note all before are also in namespace detail too.
-
-// TODO move } // namespace detail to end??
-
-
 //
 // continued_fraction_b
 // Evaluates:
@@ -88,6 +83,8 @@ namespace detail
 template <class Gen>
 typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits)
 {
+   using namespace std; // ADL of std names
+
   typedef detail::fraction_traits<Gen> traits;
   typedef typename traits::result_type result_type;
   typedef typename traits::value_type value_type;
@@ -123,6 +120,8 @@ typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g,
 template <class Gen>
 typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g, int bits, boost::uintmax_t& max_terms)
 {
+   using namespace std; // ADL of std names
+
   typedef detail::fraction_traits<Gen> traits;
   typedef typename traits::result_type result_type;
   typedef typename traits::value_type value_type;
@@ -176,6 +175,8 @@ typename detail::fraction_traits<Gen>::result_type continued_fraction_b(Gen& g,
 template <class Gen>
 typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits)
 {
+   using namespace std; // ADL of std names
+
   typedef detail::fraction_traits<Gen> traits;
   typedef typename traits::result_type result_type;
   typedef typename traits::value_type value_type;
@@ -212,6 +213,8 @@ typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g,
 template <class Gen>
 typename detail::fraction_traits<Gen>::result_type continued_fraction_a(Gen& g, int bits, boost::uintmax_t& max_terms)
 {
+   using namespace std; // ADL of std names
+
   typedef detail::fraction_traits<Gen> traits;
   typedef typename traits::result_type result_type;
   typedef typename traits::value_type value_type;
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -29,6 +29,7 @@ run test_extreme_value.cpp ;
 run test_factorials.cpp ;
 run test_fisher_f.cpp ;
 run test_gamma.cpp ;
+run test_gamma_dist.cpp ;
 run test_ibeta.cpp ;
 run test_ibeta_inv.cpp ;
 run test_ibeta_inv_ab.cpp ;
@@ -37,6 +38,7 @@ run test_igamma_inv.cpp ;
 run test_igamma_inva.cpp ;
 run test_lognormal.cpp ;
 run test_minima.cpp ;
+run test_negative_binomial.cpp ;
 run test_normal.cpp ;
 run test_poisson.cpp ;
 run test_promotion.cpp ;
@@ -49,6 +51,9 @@ run test_toms748_solve.cpp ;
 run test_weibull.cpp ;

 compile distribution_concept_check.cpp ;
+compile std_real_concept_check.cpp ;
+
+



--- a/test/powm1_sqrtp1m1_test.cpp
+++ b/test/powm1_sqrtp1m1_test.cpp
@@ -48,7 +48,7 @@ void expected_results()
      ".*",                          // platform
      ".*",                          // test type(s)
      ".*",                          // test data group
-      ".*", 3, 2);                   // test function
+      ".*", 5, 2);                   // test function

   //
   // Finish off by printing out the compiler/stdlib/platform names,
--- a/test/std_real_concept_check.cpp
+++ b/test/std_real_concept_check.cpp
@@ -0,0 +1,280 @@
+//  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)
+
+#include <boost/math/concepts/std_real_concept.hpp>
+#include <boost/math/concepts/distributions.hpp>
+
+#include <boost/math/distributions/normal.hpp>
+#include <boost/math/distributions/students_t.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/extreme_value.hpp>
+#include <boost/math/distributions/fisher_f.hpp>
+#include <boost/math/distributions/weibull.hpp>
+#include <boost/math/distributions/lognormal.hpp>
+#include <boost/math/special_functions/digamma.hpp>
+#include <boost/math/special_functions/cbrt.hpp>
+
+//
+// The purpose of this test is to verify that our code compiles
+// cleanly with a type whose std lib functions are in namespace
+// std and can *not* be found by ADL.  This verifies that we're
+// not finding std lib functions that are in the global namespace
+// for example calling ::pow(double) rather than std::pow(long double).
+// This is a silent error that does the wrong thing at runtime, and
+// of course we can't call std::pow() directly because we want
+// the functions to be found by ADL when that's appropriate.
+//
+// Furthermore our code does different things internally depending
+// on numeric_limits<>::digits, so there are some macros that can
+// be defined that cause our concept-archetype to emulate various
+// floating point types:
+//
+// EMULATE32: 32-bit float
+// EMULATE64: 64-bit double
+// EMULATE80: 80-bit long double
+// EMULATE128: 128-bit long double
+//
+// In order to ensure total code coverage this file must be
+// compiled with each of the above macros in turn, and then 
+// without any of the above as well!
+//
+
+#define NULL_MACRO /**/
+#ifdef EMULATE32
+namespace std{
+template<>
+struct numeric_limits<boost::math::concepts::std_real_concept>
+{
+   static const bool is_specialized = true;
+   static boost::math::concepts::std_real_concept min NULL_MACRO() throw();
+   static boost::math::concepts::std_real_concept max NULL_MACRO() throw();
+   static const int digits = 24;
+   static const int digits10 = 6;
+   static const bool is_signed = true;
+   static const bool is_integer = false;
+   static const bool is_exact = false;
+   static const int radix = 2;
+   static boost::math::concepts::std_real_concept epsilon() throw();
+   static boost::math::concepts::std_real_concept round_error() throw();
+   static const int min_exponent = -125;
+   static const int min_exponent10 = -37;
+   static const int max_exponent = 128;
+   static const int max_exponent10 = 38;
+   static const bool has_infinity = true;
+   static const bool has_quiet_NaN = true;
+   static const bool has_signaling_NaN = true;
+   static const float_denorm_style has_denorm = denorm_absent;
+   static const bool has_denorm_loss = false;
+   static boost::math::concepts::std_real_concept infinity() throw();
+   static boost::math::concepts::std_real_concept quiet_NaN() throw();
+   static boost::math::concepts::std_real_concept signaling_NaN() throw();
+   static boost::math::concepts::std_real_concept denorm_min() throw();
+   static const bool is_iec559 = true;
+   static const bool is_bounded = false;
+   static const bool is_modulo = false;
+   static const bool traps = false;
+   static const bool tinyness_before = false;
+   static const float_round_style round_style = round_toward_zero;
+};
+}
+#endif
+#ifdef EMULATE64
+namespace std{
+template<>
+struct numeric_limits<boost::math::concepts::std_real_concept>
+{
+   static const bool is_specialized = true;
+   static boost::math::concepts::std_real_concept min NULL_MACRO() throw();
+   static boost::math::concepts::std_real_concept max NULL_MACRO() throw();
+   static const int digits = 53;
+   static const int digits10 = 15;
+   static const bool is_signed = true;
+   static const bool is_integer = false;
+   static const bool is_exact = false;
+   static const int radix = 2;
+   static boost::math::concepts::std_real_concept epsilon() throw();
+   static boost::math::concepts::std_real_concept round_error() throw();
+   static const int min_exponent = -1021;
+   static const int min_exponent10 = -307;
+   static const int max_exponent = 1024;
+   static const int max_exponent10 = 308;
+   static const bool has_infinity = true;
+   static const bool has_quiet_NaN = true;
+   static const bool has_signaling_NaN = true;
+   static const float_denorm_style has_denorm = denorm_absent;
+   static const bool has_denorm_loss = false;
+   static boost::math::concepts::std_real_concept infinity() throw();
+   static boost::math::concepts::std_real_concept quiet_NaN() throw();
+   static boost::math::concepts::std_real_concept signaling_NaN() throw();
+   static boost::math::concepts::std_real_concept denorm_min() throw();
+   static const bool is_iec559 = true;
+   static const bool is_bounded = false;
+   static const bool is_modulo = false;
+   static const bool traps = false;
+   static const bool tinyness_before = false;
+   static const float_round_style round_style = round_toward_zero;
+};
+}
+#endif
+#ifdef EMULATE80
+namespace std{
+template<>
+struct numeric_limits<boost::math::concepts::std_real_concept>
+{
+   static const bool is_specialized = true;
+   static boost::math::concepts::std_real_concept min NULL_MACRO() throw();
+   static boost::math::concepts::std_real_concept max NULL_MACRO() throw();
+   static const int digits = 64;
+   static const int digits10 = 18;
+   static const bool is_signed = true;
+   static const bool is_integer = false;
+   static const bool is_exact = false;
+   static const int radix = 2;
+   static boost::math::concepts::std_real_concept epsilon() throw();
+   static boost::math::concepts::std_real_concept round_error() throw();
+   static const int min_exponent = -16381;
+   static const int min_exponent10 = -4931;
+   static const int max_exponent = 16384;
+   static const int max_exponent10 = 4932;
+   static const bool has_infinity = true;
+   static const bool has_quiet_NaN = true;
+   static const bool has_signaling_NaN = true;
+   static const float_denorm_style has_denorm = denorm_absent;
+   static const bool has_denorm_loss = false;
+   static boost::math::concepts::std_real_concept infinity() throw();
+   static boost::math::concepts::std_real_concept quiet_NaN() throw();
+   static boost::math::concepts::std_real_concept signaling_NaN() throw();
+   static boost::math::concepts::std_real_concept denorm_min() throw();
+   static const bool is_iec559 = true;
+   static const bool is_bounded = false;
+   static const bool is_modulo = false;
+   static const bool traps = false;
+   static const bool tinyness_before = false;
+   static const float_round_style round_style = round_toward_zero;
+};
+}
+#endif
+#ifdef EMULATE128
+namespace std{
+template<>
+struct numeric_limits<boost::math::concepts::std_real_concept>
+{
+   static const bool is_specialized = true;
+   static boost::math::concepts::std_real_concept min NULL_MACRO() throw();
+   static boost::math::concepts::std_real_concept max NULL_MACRO() throw();
+   static const int digits = 113;
+   static const int digits10 = 33;
+   static const bool is_signed = true;
+   static const bool is_integer = false;
+   static const bool is_exact = false;
+   static const int radix = 2;
+   static boost::math::concepts::std_real_concept epsilon() throw();
+   static boost::math::concepts::std_real_concept round_error() throw();
+   static const int min_exponent = -16381;
+   static const int min_exponent10 = -4931;
+   static const int max_exponent = 16384;
+   static const int max_exponent10 = 4932;
+   static const bool has_infinity = true;
+   static const bool has_quiet_NaN = true;
+   static const bool has_signaling_NaN = true;
+   static const float_denorm_style has_denorm = denorm_absent;
+   static const bool has_denorm_loss = false;
+   static boost::math::concepts::std_real_concept infinity() throw();
+   static boost::math::concepts::std_real_concept quiet_NaN() throw();
+   static boost::math::concepts::std_real_concept signaling_NaN() throw();
+   static boost::math::concepts::std_real_concept denorm_min() throw();
+   static const bool is_iec559 = true;
+   static const bool is_bounded = false;
+   static const bool is_modulo = false;
+   static const bool traps = false;
+   static const bool tinyness_before = false;
+   static const float_round_style round_style = round_toward_zero;
+};
+}
+#endif
+
+template <class RealType>
+void instantiate(RealType)
+{
+   using namespace boost;
+   using namespace boost::math;
+   using namespace boost::math::concepts;
+
+   function_requires<DistributionConcept<normal_distribution<RealType> > >();
+   function_requires<DistributionConcept<binomial_distribution<RealType> > >();
+   function_requires<DistributionConcept<cauchy_distribution<RealType> > >();
+   function_requires<DistributionConcept<chi_squared_distribution<RealType> > >();
+   function_requires<DistributionConcept<exponential_distribution<RealType> > >();
+   function_requires<DistributionConcept<extreme_value_distribution<RealType> > >();
+   function_requires<DistributionConcept<fisher_f_distribution<RealType> > >();
+   function_requires<DistributionConcept<students_t_distribution<RealType> > >();
+   function_requires<DistributionConcept<weibull_distribution<RealType> > >();
+   function_requires<DistributionConcept<lognormal_distribution<RealType> > >();
+}
+
+void test_functions()
+{
+   int i;
+   boost::math::concepts::std_real_concept v1(0.5), v2(0.5), v3(0.5);
+   boost::math::tgamma(v1);
+   boost::math::tgamma1pm1(v1);
+   boost::math::lgamma(v1);
+   boost::math::lgamma(v1, &i);
+   boost::math::digamma(v1);
+   boost::math::tgamma_ratio(v1, v2);
+   boost::math::tgamma_delta_ratio(v1, v2);
+   boost::math::factorial<boost::math::concepts::std_real_concept>(i);
+   boost::math::unchecked_factorial<boost::math::concepts::std_real_concept>(i);
+   i = boost::math::max_factorial<boost::math::concepts::std_real_concept>::value;
+   boost::math::double_factorial<boost::math::concepts::std_real_concept>(i);
+   boost::math::rising_factorial(v1, i);
+   boost::math::falling_factorial(v1, i);
+   boost::math::tgamma(v1, v2);
+   boost::math::tgamma_lower(v1, v2);
+   boost::math::gamma_P(v1, v2);
+   boost::math::gamma_Q(v1, v2);
+   boost::math::gamma_P_inv(v1, v2);
+   boost::math::gamma_Q_inv(v1, v2);
+   boost::math::gamma_P_inva(v1, v2);
+   boost::math::gamma_Q_inva(v1, v2);
+   boost::math::erf(v1);
+   boost::math::erfc(v1);
+   boost::math::erf_inv(v1);
+   boost::math::erfc_inv(v1);
+   boost::math::beta(v1, v2);
+   boost::math::beta(v1, v2, v3);
+   boost::math::betac(v1, v2, v3);
+   boost::math::ibeta(v1, v2, v3);
+   boost::math::ibetac(v1, v2, v3);
+   boost::math::ibeta_inv(v1, v2, v3);
+   boost::math::ibetac_inv(v1, v2, v3);
+   boost::math::ibeta_inva(v1, v2, v3);
+   boost::math::ibetac_inva(v1, v2, v3);
+   boost::math::ibeta_invb(v1, v2, v3);
+   boost::math::ibetac_invb(v1, v2, v3);
+   boost::math::gamma_P_derivative(v2, v3);
+   boost::math::ibeta_derivative(v1, v2, v3);
+   boost::math::fpclassify(v1);
+   boost::math::isfinite(v1);
+   boost::math::isnormal(v1);
+   boost::math::isnan(v1);
+   boost::math::isinf(v1);
+   boost::math::log1p(v1);
+   boost::math::expm1(v1);
+   boost::math::cbrt(v1);
+   boost::math::sqrt1pm1(v1);
+   boost::math::powm1(v1, v2);
+
+}
+
+
+int main()
+{
+   instantiate(boost::math::concepts::std_real_concept(0));
+}
+
--- a/test/test_binomial_coeff.cpp
+++ b/test/test_binomial_coeff.cpp
@@ -60,7 +60,7 @@ void expected_results()
      ".*",                          // platform
      largest_type,                  // test type(s)
      ".*large.*",                   // test data group
-      ".*", 50, 20);                 // test function
+      ".*", 70, 20);                 // test function
   add_expected_result(
      ".*",                          // compiler
      ".*",                          // stdlib
@@ -74,7 +74,7 @@ void expected_results()
      ".*",                          // platform
      "real_concept",                // test type(s)
      ".*",                          // test data group
-      ".*", 50, 30);                 // test function
+      ".*", 100, 30);                 // test function
   add_expected_result(
      ".*",                          // compiler
      ".*",                          // stdlib
--- a/test/test_factorials.cpp
+++ b/test/test_factorials.cpp
@@ -139,13 +139,13 @@ void test_spots(T)
      static_cast<T>(-9.76168312768123676601980433377916854311706629232503473758698e26L), tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::rising_factorial(static_cast<T>(-30.25), -21),
-      static_cast<T>(-1.50079704000923674318934280259377728203516775215430875839823e-34), tolerance);
+      static_cast<T>(-1.50079704000923674318934280259377728203516775215430875839823e-34L), tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::rising_factorial(static_cast<T>(-30.25), 5),
      static_cast<T>(-1.78799177197265625000000e7L), tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::rising_factorial(static_cast<T>(-30.25), -5),
-      static_cast<T>(-2.47177487004482195012362027432181137141899692171397467859150e-8), tolerance);
+      static_cast<T>(-2.47177487004482195012362027432181137141899692171397467859150e-8L), tolerance);
   BOOST_CHECK_CLOSE(
      ::boost::math::rising_factorial(static_cast<T>(-30.25), 6),
      static_cast<T>(4.5146792242309570312500000e8L), tolerance);
@@ -208,13 +208,13 @@ void test_spots(T)
   if(boost::math::tools::digits<T>() > 50)
   {
      BOOST_CHECK_CLOSE(
-         ::boost::math::falling_factorial(static_cast<T>(-30.75), 30),
+         ::boost::math::falling_factorial(static_cast<T>(-30.75L), 30),
         static_cast<T>(naive_falling_factorial(-30.75L, 30)),
-         tolerance);
+         tolerance * 3);
      BOOST_CHECK_CLOSE(
-         ::boost::math::falling_factorial(static_cast<T>(-30.75), 27),
+         ::boost::math::falling_factorial(static_cast<T>(-30.75L), 27),
         static_cast<T>(naive_falling_factorial(-30.75L, 27)),
-         tolerance);
+         tolerance * 3);
   }
   BOOST_CHECK_CLOSE(
      ::boost::math::falling_factorial(static_cast<T>(-12.0), 6),
--- a/test/test_gamma_dist.cpp
+++ b/test/test_gamma_dist.cpp
@@ -0,0 +1,247 @@
+// 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)
+
+// test_gamma_dist.cpp
+
+// http://en.wikipedia.org/wiki/Gamma_distribution
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm
+// Also:
+// Weisstein, Eric W. "Gamma Distribution." 
+// From MathWorld--A Wolfram Web Resource. 
+// http://mathworld.wolfram.com/GammaDistribution.html
+
+
+#define BOOST_MATH_THROW_ON_DOMAIN_ERROR
+
+#ifdef _MSC_VER
+#  pragma warning(disable: 4127) // conditional expression is constant.
+#  pragma warning(disable: 4100) // unreferenced formal parameter.
+#  pragma warning(disable: 4512) // assignment operator could not be generated.
+#  pragma warning(disable: 4510) // default constructor could not be generated.
+#  pragma warning(disable: 4610) // can never be instantiated - user defined constructor required.
+#  if !(defined _SCL_SECURE_NO_DEPRECATE) || (_SCL_SECURE_NO_DEPRECATE == 0)
+#    pragma warning(disable: 4996) // 'std::char_traits<char>::copy' was declared deprecated.
+     // #define _SCL_SECURE_NO_DEPRECATE = 1 // avoid C4996 warning.
+#  endif
+//#  pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data.
+#endif
+
+#include <boost/math/concepts/real_concept.hpp> // for real_concept
+#include <boost/test/included/test_exec_monitor.hpp> // Boost.Test
+#include <boost/test/floating_point_comparison.hpp>
+
+#include <boost/math/distributions/gamma.hpp>
+	 using boost::math::gamma_distribution;
+#include <boost/math/tools/test.hpp> 
+
+#include <iostream>
+	using std::cout;
+	using std::endl;
+	using std::setprecision;
+#include <limits>
+  using std::numeric_limits;
+
+template <class RealType>
+RealType NaivePDF(RealType shape, RealType scale, RealType x)
+{
+   // Deliberately naive PDF calculator again which
+   // we'll compare our pdf function.  However some
+   // published values to compare against would be better....
+   using namespace std;
+   RealType result = log(x) * (shape - 1) - x / scale - boost::math::lgamma(shape) - log(scale) * shape;
+   return exp(result);
+}
+
+template <class RealType>
+void check_gamma(RealType shape, RealType scale, RealType x, RealType p, RealType q, RealType tol)
+{
+   BOOST_CHECK_CLOSE(
+      ::boost::math::cdf(
+         gamma_distribution<RealType>(shape, scale),    // distribution.
+         x),                                            // random variable.
+         p,                                             // probability.
+			tol);                                          // %tolerance.
+   BOOST_CHECK_CLOSE(
+      ::boost::math::cdf(
+         complement(
+            gamma_distribution<RealType>(shape, scale), // distribution.
+            x)),                                        // random variable.
+         q,                                             // probability complement.
+			tol);                                          // %tolerance.
+   if(p < 0.999)
+   {
+      BOOST_CHECK_CLOSE(
+         ::boost::math::quantile(
+            gamma_distribution<RealType>(shape, scale),    // distribution.
+            p),                                            // probability.
+            x,                                             // random variable.
+			   tol);                                          // %tolerance.
+   }
+   if(q < 0.999)
+   {
+      BOOST_CHECK_CLOSE(
+         ::boost::math::quantile(
+            complement(
+               gamma_distribution<RealType>(shape, scale), // distribution.
+               q)),                                        // probability complement.
+            x,                                             // random variable.
+			   tol);                                          // %tolerance.
+   }
+   // PDF:
+   BOOST_CHECK_CLOSE(
+      boost::math::pdf(
+         gamma_distribution<RealType>(shape, scale),    // distribution.
+         x),                                            // random variable.
+         NaivePDF(shape, scale, x),                     // PDF
+		   tol);                                          // %tolerance.
+}
+
+template <class RealType>
+void test_spots(RealType T)
+{
+   // Basic santity checks
+   //
+	// 15 decimal places expressed as a persentage.
+   // 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; 
+	cout << "Tolerance for type " << typeid(T).name()  << " is " << tolerance << " %" << endl;
+
+   check_gamma(
+      static_cast<RealType>(0.5),
+      static_cast<RealType>(1),
+      static_cast<RealType>(0.5),
+      static_cast<RealType>(0.682689492137085),
+      static_cast<RealType>(1-0.682689492137085),
+      tolerance);
+   check_gamma(
+      static_cast<RealType>(2),
+      static_cast<RealType>(1),
+      static_cast<RealType>(0.5),
+      static_cast<RealType>(0.090204010431050),
+      static_cast<RealType>(1-0.090204010431050),
+      tolerance);
+   check_gamma(
+      static_cast<RealType>(40),
+      static_cast<RealType>(1),
+      static_cast<RealType>(10),
+      static_cast<RealType>(7.34163631456064E-13),
+      static_cast<RealType>(1-7.34163631456064E-13),
+      tolerance);
+
+   //
+   // Some more test data generated by the online
+   // calculator at http://espse.ed.psu.edu/edpsych/faculty/rhale/hale/507Mat/statlets/free/pdist.htm
+   // This has the advantage of supporting the scale parameter as well
+   // as shape, but has only a few digits accuracy, and produces
+   // some deeply suspect values if the shape parameter is < 1
+   // (it doesn't agree with MathCAD or this implementation).
+   // To be fair the incomplete gamma is tricky to get right in this area...
+   //
+	tolerance = 1e-5f * 100; // 5 decimal places as a persentage
+	cout << "Tolerance for type " << typeid(T).name()  << " is " << tolerance << " %" << endl;
+
+   check_gamma(
+      static_cast<RealType>(2),
+      static_cast<RealType>(1)/5,
+      static_cast<RealType>(0.1),
+      static_cast<RealType>(0.090204),
+      static_cast<RealType>(1-0.090204),
+      tolerance);
+   check_gamma(
+      static_cast<RealType>(2),
+      static_cast<RealType>(1)/5,
+      static_cast<RealType>(0.5),
+      static_cast<RealType>(1-0.287298),
+      static_cast<RealType>(0.287298),
+      tolerance);
+   check_gamma(
+      static_cast<RealType>(3),
+      static_cast<RealType>(2),
+      static_cast<RealType>(1),
+      static_cast<RealType>(0.014388),
+      static_cast<RealType>(1-0.014388),
+      tolerance * 10); // one less decimal place in the test value
+   check_gamma(
+      static_cast<RealType>(3),
+      static_cast<RealType>(2),
+      static_cast<RealType>(5),
+      static_cast<RealType>(0.456187),
+      static_cast<RealType>(1-0.456187),
+      tolerance);
+
+
+    RealType tol2 = boost::math::tools::epsilon<RealType>() * 5;
+    gamma_distribution<RealType> dist(8, 3);
+    RealType x = static_cast<RealType>(0.125);
+    using namespace std; // ADL of std names.
+    // mean:
+    BOOST_CHECK_CLOSE(
+       mean(dist)
+       , static_cast<RealType>(8*3), tol2);
+    // variance:
+    BOOST_CHECK_CLOSE(
+       variance(dist)
+       , static_cast<RealType>(8*3*3), tol2);
+    // std deviation:
+    BOOST_CHECK_CLOSE(
+       standard_deviation(dist)
+       , sqrt(static_cast<RealType>(8*3*3)), 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);
+    // mode:
+    BOOST_CHECK_CLOSE(
+       mode(dist)
+       , static_cast<RealType>(7 * 3), tol2);
+    // skewness:
+    BOOST_CHECK_CLOSE(
+       skewness(dist)
+       , 2 / sqrt(static_cast<RealType>(8)), tol2);
+    // kertosis:
+    BOOST_CHECK_CLOSE(
+       kurtosis(dist)
+       , 3 + 6 / static_cast<RealType>(8), tol2);
+    // kertosis excess:
+    BOOST_CHECK_CLOSE(
+       kurtosis_excess(dist)
+       , 6 / static_cast<RealType>(8), tol2);
+
+} // template <class RealType>void test_spots(RealType)
+
+int test_main(int, char* [])
+{
+   // Basic sanity-check spot values.
+	// (Parameter value, arbitrarily zero, only communicates the floating point type).
+  test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
+  test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+  test_spots(0.0L); // Test long double.
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582))
+  test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
+#endif
+#else
+   std::cout << "<note>The long double tests have been disabled on this platform "
+      "either because the long double overloads of the usual math functions are "
+      "not available at all, or because they are too inaccurate for these tests "
+      "to pass.</note>" << std::cout;
+#endif
+
+   return 0;
+} // int test_main(int, char* [])
+
+