diff --git a/include/boost/math/distributions/geometric.hpp b/include/boost/math/distributions/geometric.hpp
new file mode 100644
index 000000000..4aa51746f
--- /dev/null
+++ b/include/boost/math/distributions/geometric.hpp
@@ -0,0 +1,516 @@
+// boost\math\distributions\geometric.hpp
+
+// Copyright John Maddock 2010.
+// Copyright Paul A. Bristow 2010.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// geometric distribution is a discrete probability distribution.
+// It expresses the probability distribution of the number (k) of
+// events, occurrences, failures or arrivals before the first success.
+// supported on the set {0, 1, 2, 3...}
+
+// Note that the set includes zero (unlike some definitions that start at one).
+
+// The random variate k is the number of events, occurrences or arrivals.
+// k argument may be integral, signed, or unsigned, or floating point.
+// If necessary, it has already been promoted from an integral type.
+
+// Note that the geometric distribution
+// (like others including the binomial, geometric & Bernoulli)
+// is strictly defined as a discrete function:
+// only integral values of k are envisaged.
+// However because the method of calculation uses a continuous gamma function,
+// it is convenient to treat it as if a continous function,
+// and permit non-integral values of k.
+// To enforce the strict mathematical model, users should use floor or ceil functions
+// on k outside this function to ensure that k is integral.
+
+// See http://en.wikipedia.org/wiki/geometric_distribution
+// http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html
+// http://mathworld.wolfram.com/GeometricDistribution.html
+
+#ifndef BOOST_MATH_SPECIAL_GEOMETRIC_HPP
+#define BOOST_MATH_SPECIAL_GEOMETRIC_HPP
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x) == Ix(a, b).
+#include <boost/math/distributions/complement.hpp> // complement.
+#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks domain_error & logic_error.
+#include <boost/math/special_functions/fpclassify.hpp> // isnan.
+#include <boost/math/tools/roots.hpp> // for root finding.
+#include <boost/math/distributions/detail/inv_discrete_quantile.hpp>
+
+#include <boost/type_traits/is_floating_point.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/mpl/if.hpp>
+
+#include <limits> // using std::numeric_limits;
+#include <utility>
+
+#if defined (BOOST_MSVC)
+#  pragma warning(push)
+// This believed not now necessary, so commented out.
+//#  pragma warning(disable: 4702) // unreachable code.
+// in domain_error_imp in error_handling.
+#endif
+
+namespace boost
+{
+  namespace math
+  {
+    namespace geometric_detail
+    {
+      // Common error checking routines for geometric distribution function:
+      template <class RealType, class Policy>
+      inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol)
+      {
+        if( !(boost::math::isfinite)(p) || (p < 0) || (p > 1) )
+        {
+          *result = policies::raise_domain_error<RealType>(
+            function,
+            "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol);
+          return false;
+        }
+        return true;
+      }
+
+      template <class RealType, class Policy>
+      inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& pol)
+      {
+        return check_success_fraction(function, p, result, pol);
+      }
+
+      template <class RealType, class Policy>
+      inline bool check_dist_and_k(const char* function,  const RealType& p, RealType k, RealType* result, const Policy& pol)
+      {
+        if(check_dist(function, p, result, pol) == false)
+        {
+          return false;
+        }
+        if( !(boost::math::isfinite)(k) || (k < 0) )
+        { // Check k failures.
+          *result = policies::raise_domain_error<RealType>(
+            function,
+            "Number of failures argument is %1%, but must be >= 0 !", k, pol);
+          return false;
+        }
+        return true;
+      } // Check_dist_and_k
+
+      template <class RealType, class Policy>
+      inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& pol)
+      {
+        if(check_dist(function, p, result, pol) && detail::check_probability(function, prob, result, pol) == false)
+        {
+          return false;
+        }
+        return true;
+      } // check_dist_and_prob
+    } //  namespace geometric_detail
+
+    template <class RealType = double, class Policy = policies::policy<> >
+    class geometric_distribution
+    {
+    public:
+      typedef RealType value_type;
+      typedef Policy policy_type;
+
+      geometric_distribution(RealType p) : m_p(p)
+      { // Constructor stores success_fraction p.
+        RealType result;
+        geometric_detail::check_dist(
+          "geometric_distribution<%1%>::geometric_distribution",
+          m_p, // Check success_fraction 0 <= p <= 1.
+          &result, Policy());
+      } // geometric_distribution constructor.
+
+      // Private data getter class member functions.
+      RealType success_fraction() const
+      { // Probability of success as fraction in range 0 to 1.
+        return m_p;
+      }
+      RealType successes() const
+      { // Total number of successes r = 1 (for compatibility with negative binomial?).
+        return 1;
+      }
+
+      // Parameter estimation.
+      // (These are copies of negative_binomial distribution with successes = 1).
+      static RealType find_lower_bound_on_p(
+        RealType trials,
+        RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test.
+      {
+        static const char* function = "boost::math::geometric<%1%>::find_lower_bound_on_p";
+        RealType result;  // of error checks.
+        RealType successes = 1;
+        RealType failures = trials - successes;
+        if(false == detail::check_probability(function, alpha, &result, Policy())
+          && geometric_detail::check_dist_and_k(
+          function, RealType(0), failures, &result, Policy()))
+        {
+          return result;
+        }
+        // Use complement ibeta_inv function for lower bound.
+        // This is adapted from the corresponding binomial formula
+        // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
+        // This is a Clopper-Pearson interval, and may be overly conservative,
+        // see also "A Simple Improved Inferential Method for Some
+        // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY
+        // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf
+        //
+        return ibeta_inv(successes, failures + 1, alpha, static_cast<RealType*>(0), Policy());
+      } // find_lower_bound_on_p
+
+      static RealType find_upper_bound_on_p(
+        RealType trials,
+        RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test.
+      {
+        static const char* function = "boost::math::geometric<%1%>::find_upper_bound_on_p";
+        RealType result;  // of error checks.
+        RealType successes = 1;
+        RealType failures = trials - successes;
+        if(false == geometric_detail::check_dist_and_k(
+          function, RealType(0), failures, &result, Policy())
+          && detail::check_probability(function, alpha, &result, Policy()))
+        {
+          return result;
+        }
+        if(failures == 0)
+        {
+           return 1;
+        }// Use complement ibetac_inv function for upper bound.
+        // Note adjusted failures value: *not* failures+1 as usual.
+        // This is adapted from the corresponding binomial formula
+        // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
+        // This is a Clopper-Pearson interval, and may be overly conservative,
+        // see also "A Simple Improved Inferential Method for Some
+        // Discrete Distributions" Yong CAI and K. Krishnamoorthy
+        // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf
+        //
+        return ibetac_inv(successes, failures, alpha, static_cast<RealType*>(0), Policy());
+      } // find_upper_bound_on_p
+
+      // Estimate number of trials :
+      // "How many trials do I need to be P% sure of seeing k or fewer failures?"
+
+      static RealType find_minimum_number_of_trials(
+        RealType k,     // number of failures (k >= 0).
+        RealType p,     // success fraction 0 <= p <= 1.
+        RealType alpha) // risk level threshold 0 <= alpha <= 1.
+      {
+        static const char* function = "boost::math::geometric<%1%>::find_minimum_number_of_trials";
+        // Error checks:
+        RealType result;
+        if(false == geometric_detail::check_dist_and_k(
+          function, p, k, &result, Policy())
+          && detail::check_probability(function, alpha, &result, Policy()))
+        {
+          return result;
+        }
+        result = ibeta_inva(k + 1, p, alpha, Policy());  // returns n - k
+        return result + k;
+      } // RealType find_number_of_failures
+
+      static RealType find_maximum_number_of_trials(
+        RealType k,     // number of failures (k >= 0).
+        RealType p,     // success fraction 0 <= p <= 1.
+        RealType alpha) // risk level threshold 0 <= alpha <= 1.
+      {
+        static const char* function = "boost::math::geometric<%1%>::find_maximum_number_of_trials";
+        // Error checks:
+        RealType result;
+        if(false == geometric_detail::check_dist_and_k(
+          function, p, k, &result, Policy())
+          &&  detail::check_probability(function, alpha, &result, Policy()))
+        { 
+          return result;
+        }
+        result = ibetac_inva(k + 1, p, alpha, Policy());  // returns n - k
+        return result + k;
+      } // RealType find_number_of_trials complemented
+
+    private:
+      //RealType m_r; // successes fixed at unity.
+      RealType m_p; // success_fraction
+    }; // template <class RealType, class Policy> class geometric_distribution
+
+    typedef geometric_distribution<double> geometric; // Reserved name of type double.
+
+    template <class RealType, class Policy>
+    inline const std::pair<RealType, RealType> range(const geometric_distribution<RealType, Policy>& /* dist */)
+    { // Range of permissible values for random variable k.
+       using boost::math::tools::max_value;
+       return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer?
+    }
+
+    template <class RealType, class Policy>
+    inline const std::pair<RealType, RealType> support(const geometric_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>(static_cast<RealType>(0),  max_value<RealType>()); // max_integer?
+    }
+
+    template <class RealType, class Policy>
+    inline RealType mean(const geometric_distribution<RealType, Policy>& dist)
+    { // Mean of geometric distribution = (1-p)/p.
+      return (1 - dist.success_fraction() ) / dist.success_fraction();
+    } // mean
+
+    // median implemented via quantile(half) in derived accessors.
+
+    template <class RealType, class Policy>
+    inline RealType mode(const geometric_distribution<RealType, Policy>&)
+    { // Mode of geometric distribution = zero.
+      BOOST_MATH_STD_USING // ADL of std functions.
+      return 0;
+    } // mode
+    
+    template <class RealType, class Policy>
+    inline RealType variance(const geometric_distribution<RealType, Policy>& dist)
+    { // Variance of Binomial distribution = (1-p) / p^2.
+      return  (1 - dist.success_fraction())
+        / (dist.success_fraction() * dist.success_fraction());
+    } // variance
+
+    template <class RealType, class Policy>
+    inline RealType skewness(const geometric_distribution<RealType, Policy>& dist)
+    { // skewness of geometric distribution = 2-p / (sqrt(r(1-p))
+      BOOST_MATH_STD_USING // ADL of std functions.
+      RealType p = dist.success_fraction();
+      return (2 - p) / sqrt(1 - p);
+    } // skewness
+
+    template <class RealType, class Policy>
+    inline RealType kurtosis(const geometric_distribution<RealType, Policy>& dist)
+    { // kurtosis of geometric distribution
+      // http://en.wikipedia.org/wiki/geometric is kurtosis_excess so add 3
+      RealType p = dist.success_fraction();
+      return 3 + (p*p - 6*p + 6) / (1 - p);
+    } // kurtosis
+
+     template <class RealType, class Policy>
+    inline RealType kurtosis_excess(const geometric_distribution<RealType, Policy>& dist)
+    { // kurtosis excess of geometric distribution
+      // http://mathworld.wolfram.com/Kurtosis.html table of kurtosis_excess
+      RealType p = dist.success_fraction();
+      return (p*p - 6*p + 6) / (1 - p);
+    } // kurtosis_excess
+
+    // RealType standard_deviation(const geometric_distribution<RealType, Policy>& dist)
+    // standard_deviation provided by derived accessors.
+    // RealType hazard(const geometric_distribution<RealType, Policy>& dist)
+    // hazard of geometric distribution provided by derived accessors.
+    // RealType chf(const geometric_distribution<RealType, Policy>& dist)
+    // chf of geometric distribution provided by derived accessors.
+
+    template <class RealType, class Policy>
+    inline RealType pdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)
+    { // Probability Density/Mass Function.
+      BOOST_FPU_EXCEPTION_GUARD
+      BOOST_MATH_STD_USING  // For ADL of math functions.
+      static const char* function = "boost::math::pdf(const geometric_distribution<%1%>&, %1%)";
+
+      RealType p = dist.success_fraction();
+      RealType result;
+      if(false == geometric_detail::check_dist_and_k(
+        function,
+        p,
+        k,
+        &result, Policy()))
+      {
+        return result;
+      }
+      if (k == 0)
+      {
+        return p; // success_fraction
+      }
+      RealType q = 1 - p;  // Inaccurate for small p?
+      // So try to avoid inaccuracy for large or small p.
+      // but has little effect > last significant bit.
+      //cout << "p *  pow(q, k) " << result << endl; // seems best whatever p
+      //cout << "exp(p * k * log1p(-p)) " << p * exp(k * log1p(-p)) << endl;
+      //if (p < 0.5)
+      //{
+      //  result = p *  pow(q, k);
+      //}
+      //else
+      //{
+      //  result = p * exp(k * log1p(-p));
+      //}
+      result = p * pow(q, k);
+      return result;
+    } // geometric_pdf
+
+    template <class RealType, class Policy>
+    inline RealType cdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)
+    { // Cumulative Distribution Function of geometric.
+      static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";
+
+      // k argument may be integral, signed, or unsigned, or floating point.
+      // If necessary, it has already been promoted from an integral type.
+      RealType p = dist.success_fraction();
+      // Error check:
+      RealType result;
+      if(false == geometric_detail::check_dist_and_k(
+        function,
+        p,
+        k,
+        &result, Policy()))
+      {
+        return result;
+      }
+      if(k == 0)
+      {
+        return p; // success_fraction
+      }
+      //RealType q = 1 - p;  // Bad for small p
+      //RealType probability = 1 - std::pow(q, k+1);
+
+      RealType z = log1p(-p) * (k+1);
+      RealType probability = -expm1(z);
+
+      return probability;
+    } // cdf Cumulative Distribution Function geometric.
+
+      template <class RealType, class Policy>
+      inline RealType cdf(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)
+      { // Complemented Cumulative Distribution Function geometric.
+
+      static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";
+      // k argument may be integral, signed, or unsigned, or floating point.
+      // If necessary, it has already been promoted from an integral type.
+      RealType const& k = c.param;
+      geometric_distribution<RealType, Policy> const& dist = c.dist;
+      RealType p = dist.success_fraction();
+      // Error check:
+      RealType result;
+      if(false == geometric_detail::check_dist_and_k(
+        function,
+        p,
+        k,
+        &result, Policy()))
+      {
+        return result;
+      }
+      RealType z = log1p(-p) * (k+1);
+      RealType probability = exp(z);
+      return probability;
+    } // cdf Complemented Cumulative Distribution Function geometric.
+
+    template <class RealType, class Policy>
+    inline RealType quantile(const geometric_distribution<RealType, Policy>& dist, const RealType& x)
+    { // Quantile, percentile/100 or Percent Point geometric function.
+      // Return the number of expected failures k for a given probability p.
+
+      // Inverse cumulative Distribution Function or Quantile (percentile / 100) of geometric Probability.
+      // k argument may be integral, signed, or unsigned, or floating point.
+
+      static const char* function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)";
+      BOOST_MATH_STD_USING // ADL of std functions.
+
+      RealType success_fraction = dist.success_fraction();
+      // Check dist and x.
+      RealType result;
+      if(false == geometric_detail::check_dist_and_prob
+        (function, success_fraction, x, &result, Policy()))
+      {
+        return result;
+      }
+
+      // Special cases.
+      if (x == 1)
+      {  // Would need +infinity failures for total confidence.
+        result = policies::raise_overflow_error<RealType>(
+            function,
+            "Probability argument is 1, which implies infinite failures !", Policy());
+        return result;
+       // usually means return +std::numeric_limits<RealType>::infinity();
+       // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR
+      }
+      if (x == 0)
+      { // No failures are expected if P = 0.
+        return 0; // Total trials will be just dist.successes.
+      }
+      // if (P <= pow(dist.success_fraction(), 1))
+      if (x <= success_fraction)
+      { // p <= pdf(dist, 0) == cdf(dist, 0)
+        return 0;
+      }
+      if (x == 1)
+      {
+        return 0;
+      }
+   
+      // log(1-x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small
+      result = log1p(-x) / log1p(-success_fraction) -1;
+      // Subtract a few epsilons here too?
+      // to make sure it doesn't slip over, so ceil would be one too many.
+      return result;
+    } // RealType quantile(const geometric_distribution dist, p)
+
+    template <class RealType, class Policy>
+    inline RealType quantile(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)
+    {  // Quantile or Percent Point Binomial function.
+       // Return the number of expected failures k for a given
+       // complement of the probability Q = 1 - P.
+       static const char* function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)";
+
+       // Error checks:
+       RealType x = c.param;
+       const geometric_distribution<RealType, Policy>& dist = c.dist;
+       RealType success_fraction = dist.success_fraction();
+       RealType result;
+       if(false == geometric_detail::check_dist_and_prob(
+          function,
+          success_fraction,
+          x,
+          &result, Policy()))
+       {
+          return result;
+       }
+
+       // Special cases:
+       if(x == 1)
+       {  // There may actually be no answer to this question,
+          // since the probability of zero failures may be non-zero,
+          return 0; // but zero is the best we can do:
+       }
+       if (-x <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy()))
+       {  // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
+          return 0; //
+       }
+       if(x == 0)
+       {  // Probability 1 - Q  == 1 so infinite failures to achieve certainty.
+          // Would need +infinity failures for total confidence.
+          result = policies::raise_overflow_error<RealType>(
+             function,
+             "Probability argument complement is 0, which implies infinite failures !", Policy());
+          return result;
+          // usually means return +std::numeric_limits<RealType>::infinity();
+          // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR
+       }
+       // log(x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small
+       result = log(x) / log1p(-success_fraction) -1;
+      return result;
+
+    } // 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>
+
+#if defined (BOOST_MSVC)
+# pragma warning(pop)
+#endif
+
+#endif // BOOST_MATH_SPECIAL_GEOMETRIC_HPP
diff --git a/include/boost/math/distributions/inverse_gaussian.hpp b/include/boost/math/distributions/inverse_gaussian.hpp
new file mode 100644
index 000000000..05a60178f
--- /dev/null
+++ b/include/boost/math/distributions/inverse_gaussian.hpp
@@ -0,0 +1,578 @@
+//  Copyright John Maddock 2010.
+//  Copyright Paul A. Bristow 2010.
+
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_STATS_INVERSE_GAUSSIAN_HPP
+#define BOOST_STATS_INVERSE_GAUSSIAN_HPP
+
+// http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution
+// http://mathworld.wolfram.com/InverseGaussianDistribution.html
+
+// The normal-inverse Gaussian distribution
+// also called the Wald distribution (some sources limit this to when mean = 1).
+
+// It is the continuous probability distribution
+// that is defined as the normal variance-mean mixture where the mixing density is the 
+// inverse Gaussian distribution. The tails of the distribution decrease more slowly
+// than the normal distribution. It is therefore suitable to model phenomena
+// where numerically large values are more probable than is the case for the normal distribution.
+
+// The Inverse Gaussian distribution was first studied in relationship to Brownian motion.
+// In 1956 M.C.K. Tweedie used the name 'Inverse Gaussian' because there is an inverse 
+// relationship between the time to cover a unit distance and distance covered in unit time.
+
+// Examples are returns from financial assets and turbulent wind speeds. 
+// The normal-inverse Gaussian distributions form
+// a subclass of the generalised hyperbolic distributions.
+
+// See also
+
+// http://en.wikipedia.org/wiki/Normal_distribution
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
+// Also:
+// Weisstein, Eric W. "Normal Distribution."
+// From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/NormalDistribution.html
+
+// http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.
+// ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/
+
+// http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/inverse_gaussian.html
+// R package for dinverse_gaussian, ...
+
+// http://www.statsci.org/s/inverse_gaussian.s  and http://www.statsci.org/s/inverse_gaussian.html
+
+//#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/erf.hpp> // for erf/erfc.
+#include <boost/math/distributions/complement.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/normal.hpp>
+#include <boost/math/distributions/gamma.hpp> // for gamma function
+// using boost::math::gamma_p;
+
+#include <boost/math/tr1.hpp>
+//using std::tr1::tuple;
+//using std::tr1::make_tuple;
+#include <boost/math/tools/roots.hpp>
+//using boost::math::tools::newton_raphson_iterate;
+
+#include <utility>
+
+namespace boost{ namespace math{
+
+template <class RealType = double, class Policy = policies::policy<> >
+class inverse_gaussian_distribution
+{
+public:
+   typedef RealType value_type;
+   typedef Policy policy_type;
+
+   inverse_gaussian_distribution(RealType mean = 1, RealType scale = 1)
+      : m_mean(mean), m_scale(scale)
+   { // Default is a 1,1 inverse_gaussian distribution.
+     static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution";
+
+     RealType result;
+     detail::check_scale(function, scale, &result, Policy());
+     detail::check_location(function, mean, &result, Policy());
+   }
+
+   RealType mean()const
+   { // alias for location.
+      return m_mean; // aka mu
+   }
+
+   // Synonyms, provided to allow generic use of find_location and find_scale.
+   RealType location()const
+   { // location, aka mu.
+      return m_mean;
+   }
+   RealType scale()const
+   { // scale, aka lambda.
+      return m_scale;
+   }
+
+   RealType shape()const
+   { // shape, aka phi = lambda/mu.
+      return m_scale / m_mean;
+   }
+
+private:
+   //
+   // Data members:
+   //
+   RealType m_mean;  // distribution mean or location, aka mu.
+   RealType m_scale;    // distribution standard deviation or scale, aka lambda.
+}; // class normal_distribution
+
+typedef inverse_gaussian_distribution<double> inverse_gaussian;
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> range(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
+{ // Range of permissible values for random variable x, zero to max.
+   using boost::math::tools::max_value;
+   return std::pair<RealType, RealType>(static_cast<RealType>(0.), max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> support(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
+{ // Range of supported values for random variable x, zero to max.
+  // 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>(static_cast<RealType>(0.),  max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline RealType pdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Probability Density Function
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType scale = dist.scale();
+   RealType mean = dist.mean();
+   RealType result;
+   static const char* function = "boost::math::pdf(const inverse_gaussian_distribution<%1%>&, %1%)";
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_location(function, mean, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_positive_x(function, x, &result, Policy()))
+   {
+      return std::numeric_limits<RealType>::quiet_NaN();
+   }
+
+    if (x == 0)
+    {
+      return 0; // Convenient, even if not defined mathematically.
+    }
+
+   result =
+     sqrt(scale / (constants::two_pi<RealType>() * x * x * x))
+    * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
+   return result;
+} // pdf
+
+template <class RealType, class Policy>
+inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Cumulative Density Function.
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType scale = dist.scale();
+   RealType mean = dist.mean();
+   static const char* function = "boost::math::cdf(const inverse_gaussian_distribution<%1%>&, %1%)";
+   RealType result;
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_location(function, mean, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_positive_x(function, x, &result, Policy()))
+   {
+     return result;
+   }
+   if (x == 0)
+   {
+     return 0; // Convenient, even if not defined mathematically.
+   }
+   // Problem with this formula for large scale > 1000 or small x, 
+   //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two<RealType>(), Policy()) + 1)
+   //  + exp(2 * scale / mean) / 2 
+   //  * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));
+   // so use normal distribution version:
+   // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution.
+
+   normal_distribution<RealType> n01;
+
+   RealType n0 = sqrt(scale / x);
+   n0 *= ((x / mean) -1);
+   RealType n1 = cdf(n01, n0);
+   RealType expfactor = exp(2 * scale / mean);
+   RealType n3 = - sqrt(scale / x);
+   n3 *= (x / mean) + 1;
+   // cout << "((x / mean) +1) = " << n3 << endl;
+   RealType n4 = cdf(n01, n3);
+   //cout << "phi1 = " << n1 << ", exp(2 * scale / mean) = " << n2 << ", exp * phi2 = " << n4 * n2 << endl;
+
+   result = n1 + expfactor * n4;
+
+   if(false)
+   { // Output some diagnostic values.
+     cout <<"_\n    cdf===========================" << endl;
+     cout << "sqrt(scale / x)*((x / mean) -1) = " << n0 << endl;
+     cout << "cdf(n01, n1) = " << n1 << endl;
+     cout << "exp(2 * scale / mean) = " << expfactor << endl;
+     cout << " - sqrt(scale / x)*((x / mean) +1) = " << n3 << endl;
+     cout << "cdf(n01, - sqrt(scale / x)*((x / mean) +1)) = " << n4 << endl;
+     cout << "exp * cdf_2 = " << n4 * expfactor << endl;
+     cout << "cdf_1 + exp * cdf_2 = " << result << endl;
+   }
+   return result;
+} // cdf
+
+template <class RealType>
+struct inverse_gaussian_quantile_functor
+{ 
+
+  inverse_gaussian_quantile_functor(const boost::math::inverse_gaussian_distribution<RealType> dist, RealType const& p)
+    : distribution(dist), prob(p)
+  {
+  }
+  boost::math::tuple<RealType, RealType> operator()(RealType const& x)
+  {
+    RealType c = cdf(distribution, x);
+    RealType fx = c - prob;  // Difference cdf - value - to minimize.
+    RealType dx = pdf(distribution, x); // pdf is 1st derivative.
+    if(false)
+    {
+     cout << "cdf(dist, " << x << ") = " << c << ", diff = " << fx << ", dx " << dx << endl;
+    }
+    // return both function evaluation difference f(x) and 1st derivative f'(x).
+    return std::tr1::make_tuple(fx, dx);
+  }
+  private:
+  const boost::math::inverse_gaussian_distribution<RealType> distribution;
+  RealType prob; 
+};
+
+template <class RealType>
+struct inverse_gaussian_quantile_complement_functor
+{ 
+    inverse_gaussian_quantile_complement_functor(const boost::math::inverse_gaussian_distribution<RealType> dist, RealType const& p)
+    : distribution(dist), prob(p)
+  {
+  }
+  boost::math::tuple<RealType, RealType> operator()(RealType const& x)
+  {
+    RealType c = cdf(complement(distribution, x));
+    RealType fx = c - prob;  // Difference cdf - value - to minimize.
+    RealType dx = -pdf(distribution, x); // pdf is 1st derivative.
+    if(true)
+    {
+      std::streamsize precision = cout.precision(); // Save
+      if (false)
+      {
+        cout << setprecision(numeric_limits<RealType>::max_digits10) 
+         << "cdf((complement(dist, " << x << ")) = " << c 
+         << setprecision(4)
+         << ", diff = " << fx << ", dx " << dx << endl;
+        cout.precision(precision);  // restore
+      }
+    }
+    // return both function evaluation difference f(x) and 1st derivative f'(x).
+    return std::tr1::make_tuple(fx, dx);
+  }
+  private:
+  const boost::math::inverse_gaussian_distribution<RealType> distribution;
+  RealType prob; 
+};
+
+namespace detail
+{
+  template <class RealType>
+  inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1)
+  { // guess at random variate value x for inverse gaussian quantile.
+      using boost::math::policies::policy;
+      // Error type.
+      using boost::math::policies::overflow_error;
+      // Action.
+      using boost::math::policies::ignore_error;
+
+      typedef policy<
+        overflow_error<ignore_error> // Ignore overflow (return infinity)
+      > no_overthrow_policy;
+
+    RealType x; // result is guess at random variate value x.
+    RealType phi = lambda / mu;
+    if (phi > 2.)
+    { // Big phi, so starting to look like normal Gaussian distribution.
+      //    x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu);
+      // Whitmore, G.A. and Yalovsky, M.
+      // A normalising logarithmic transformation for inverse Gaussian random variables,
+      // Technometrics 20-2, 207-208 (1978), but using expression from
+      // V Seshadri, Inverse Gaussian distribution (1998) ISBN 0387 98618 9, page 6.
+      //x = qnorm(p, 0, 1, true, false);
+      //x /= sqrt(phi);
+      //x = x - 1. / (2 * phi);
+      //x = mu * exp(x);
+      // x = mu * exp(qnorm(p, 0, 1, true, false) / sqrt(phi) - 1/(2 * phi));
+      normal_distribution<RealType, no_overthrow_policy> n01;
+      x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi));
+      // Might add about 0.006 to this to get closer?
+      //RealType pi = 3.1459;
+      //cout << "1 / sqrt(8 * pi * phi) = " << 1 / sqrt(8 * pi * phi) << endl;
+      // Bagshaw guess is:
+      // RealType U = quantile(n01, p); // U <- qnorm (p)
+      // RealType r1 = 1 + U / sqrt (phi) + U * U / (2 * phi) + U * U * U / (8 * phi * sqrt(phi));
+    }
+    else
+    { // phi < 2 so much less symmetrical with long tail,
+      // so use gamma distribution as an approximation.
+      using boost::math::gamma_distribution;
+
+      // Define the distribution, using gamma_nooverflow:
+      typedef gamma_distribution<RealType, no_overthrow_policy> gamma_nooverflow;
+
+      gamma_distribution<RealType, no_overthrow_policy> g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
+
+      // gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
+      // R qgamma(0.2, 0.5, 1)  0.0320923
+      RealType qg = quantile(complement(g, p));
+      //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false);
+      //cout << "quantile(complement(g, p)) = " << qg << ", qgamma(1.- p, 0.5, 1.0, true, false); = " << qg1 << endl; // 49.014664823030209 1.#INF
+      x = lambda / (qg * 2);
+      // 
+      if (x > mu/2) // x > mu /2?
+      { // x too large for the gamma approximation to work well.
+        //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807
+        RealType q = quantile(g, p);
+        // cout << "quantile((g, p)) = " << q << endl;// 49.014664823030209 
+          //<< ", qgamma(1.- p, 0.5, 1.0); = " << x << endl; // 1.#INF
+       // x = mu * exp(q * static_cast<RealType>(0.1));  // Said to improve at high p
+       // x = mu * x;  // Improves at high p?
+        x = mu * exp(q / sqrt(phi) - 1/(2 * phi));
+      }
+    }
+    return x;
+  }  // guess_ig
+} // namespace detail
+
+template <class RealType, class Policy>
+inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& p)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions.
+   // No closed form exists so guess and use Newton Raphson iteration.
+
+   RealType mean = dist.mean();
+   RealType scale = dist.scale();
+   static const char* function = "boost::math::quantile(const inverse_gaussian_distribution<%1%>&, %1%)";
+
+   RealType result;
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+      return result;
+   if(false == detail::check_location(function, mean, &result, Policy()))
+      return result;
+   if(false == detail::check_probability(function, p, &result, Policy()))
+      return result;
+   if (p == 0)
+   {
+     return 0; // Convenient, even if not defined mathematically?
+   }
+   if (p == 1)
+   { // Might not return infinity?
+     return numeric_limits<RealType>::infinity();
+   }
+  //RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1);
+
+  RealType guess = detail::guess_ig(p, dist.mean(), dist.scale());
+  using boost::math::tools::max_value;
+
+  RealType min = 0.; // Minimum possible value is bottom of range of distribution.
+  RealType max = max_value<RealType>();// Maximum possible value is top of range. 
+  // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
+  // digits used to control how accurate to try to make the result.
+  // To allow user to control accuracy versus speed,
+  int get_digits = policies::digits<RealType, Policy>();// get digits from policy, 
+  boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
+  if(false)
+  {
+  cout << "Probability " << p << ", guess " << guess 
+    << ", min " << min << ", max " << max 
+    //<< ", std::numeric_limits<" << typeid(RealType).name()  << ">::digits = " << digits 
+     << ", accuracy " << get_digits << " bits."
+     << ", max iterations set by policy " << m
+     << endl;
+  }
+  using boost::math::tools::newton_raphson_iterate;
+  result =
+    newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType>(dist, p), guess, min, max, get_digits, m);
+   //cout << m << " iterations." << endl;
+   return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType cdf(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType scale = c.dist.scale();
+   RealType mean = c.dist.mean();
+   RealType x = c.param;
+   static const char* function = "boost::math::cdf(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
+   // infinite arguments not supported.
+   //if((boost::math::isinf)(x))
+   //{
+   //  if(x < 0) return 1; // cdf complement -infinity is unity.
+   //  return 0; // cdf complement +infinity is zero
+   //}
+   // These produce MSVC 4127 warnings, so the above used instead.
+   //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
+   //{ // cdf complement +infinity is zero.
+   //  return 0;
+   //}
+   //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
+   //{ // cdf complement -infinity is unity.
+   //  return 1;
+   //}
+   RealType result;
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+      return result;
+   if(false == detail::check_location(function, mean, &result, Policy()))
+      return result;
+   if(false == detail::check_x(function, x, &result, Policy()))
+      return result;
+
+   normal_distribution<RealType> n01;
+   RealType n0 = sqrt(scale / x);
+   n0 *= ((x / mean) -1);
+   RealType cdf_1 = cdf(complement(n01, n0));
+
+   RealType expfactor = exp(2 * scale / mean);
+   RealType n3 = - sqrt(scale / x);
+   n3 *= (x / mean) + 1;
+
+   //RealType n5 = +sqrt(scale/x) * ((x /mean) + 1); // note now positive sign.
+   RealType n6 = cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1)));
+   RealType n4 = cdf(n01, n3); // = 
+   result = cdf_1 - expfactor * n6; 
+   if(false)
+   {
+     cout <<"_\n      cdf(complement ===========================" << endl;
+     cout << "sqrt(scale / x)*((x / mean) -1) = " << n0 << endl;
+     cout << "cdf(complement(n01, n1)) = " << cdf_1 << endl;
+     cout << "-sqrt(scale / x) * ((x / mean) +1) = " << n3 << endl;
+     cout << "exp(2 * scale / mean) = " << expfactor << endl;
+     cout << "cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1))) = " << n6 << endl;
+     cout << "cdf((n01, ) exp(2 * scale / mean) * (x / mean) + 1) = " << n4 << endl;
+     cout << "exp * cdf_2 = " << result << endl;
+   }
+   //cout << "cdf(complement) result = " << result << endl;
+   return result;
+} // cdf complement
+
+template <class RealType, class Policy>
+inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType scale = c.dist.scale();
+   RealType mean = c.dist.mean();
+   static const char* function = "boost::math::quantile(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
+   RealType result;
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+      return result;
+   if(false == detail::check_location(function, mean, &result, Policy()))
+      return result;
+   RealType q = c.param;
+   if(false == detail::check_probability(function, q, &result, Policy()))
+      return result;
+
+   RealType guess = detail::guess_ig(q, mean, scale);
+   // Complement.
+   using boost::math::tools::max_value;
+
+  RealType min = 0.; // Minimum possible value is bottom of range of distribution.
+  RealType max = max_value<RealType>();// Maximum possible value is top of range. 
+  // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
+  // digits used to control how accurate to try to make the result.
+  int get_digits = policies::digits<RealType, Policy>();
+  boost::uintmax_t m = policies::get_max_root_iterations<Policy>();
+  if(false)
+  {
+    cout << "Probability " << q << ", guess at x = " << guess 
+    //<< ", min " << min << ", max " << max 
+    ////<< ", std::numeric_limits<" << typeid(RealType).name()  << ">::digits = " << digits 
+    // << ", accuracy " << get_digits << " bits."
+    // << ", max iterations set by policy " << m
+     << endl;
+  }
+  using boost::math::tools::newton_raphson_iterate;
+  result =
+    newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType>(c.dist, q), guess, min, max, get_digits, m);
+   //cout << m << " iterations." << endl;
+   return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType mean(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{ // aka mu
+   return dist.mean();
+}
+
+template <class RealType, class Policy>
+inline RealType scale(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{ // aka lambda
+   return dist.scale();
+}
+
+template <class RealType, class Policy>
+inline RealType shape(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{ // aka phi
+   return dist.shape();
+}
+
+template <class RealType, class Policy>
+inline RealType standard_deviation(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType mean = dist.mean();
+  RealType result = sqrt(mean * mean * mean / scale);
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType mode(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale)) 
+      - 3 * mean / (2 * scale));
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType skewness(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = 3 * sqrt(mean/scale);
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = 15 * mean / scale -3;
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis_excess(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = 15 * mean / scale;
+  return result;
+}
+
+} // 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_INVERSE_GAUSSIAN_HPP
+
+
diff --git a/include/boost/math/distributions/inverse_normal.hpp b/include/boost/math/distributions/inverse_normal.hpp
new file mode 100644
index 000000000..731d568de
--- /dev/null
+++ b/include/boost/math/distributions/inverse_normal.hpp
@@ -0,0 +1,361 @@
+//  Copyright John Maddock 2010.
+//  Copyright Paul A. Bristow 2010.
+
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_STATS_INVERSE_NORMAL_HPP
+#define BOOST_STATS_INVERSE_NORMAL_HPP
+
+// http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution
+// http://mathworld.wolfram.com/InverseGaussianDistribution.html
+
+// The normal-inverse Gaussian distribution (also called the Wald distribution when mean = 1)
+// is the continuous probability distribution
+// that is defined as the normal variance-mean mixture where the mixing density is the 
+// inverse Gaussian distribution. The tails of the distribution decrease more slowly
+// than the normal distribution. It is therefore suitable to model phenomena
+// where numerically large values are more probable than is the case for the normal distribution.
+
+// Examples are returns from financial assets and turbulent wind speeds. 
+// The normal-inverse Gaussian distributions form
+// a subclass of the generalised hyperbolic distributions.
+
+// See also
+
+// http://en.wikipedia.org/wiki/Normal_distribution
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
+// Also:
+// Weisstein, Eric W. "Normal Distribution."
+// From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/NormalDistribution.html
+
+// http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.
+// ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/
+
+// http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/invGauss.html
+// R package for dinvGauss, ...
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/erf.hpp> // for erf/erfc.
+#include <boost/math/distributions/complement.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/normal.hpp>
+
+#include <utility>
+
+namespace boost{ namespace math{
+
+template <class RealType = double, class Policy = policies::policy<> >
+class inverse_normal_distribution
+{
+public:
+   typedef RealType value_type;
+   typedef Policy policy_type;
+
+   inverse_normal_distribution(RealType mean = 1, RealType sd = 1)
+      : m_mean(mean), m_sd(sd)
+   { // Default is a 1,1 inverse_normal distribution.
+     static const char* function = "boost::math::inverse_normal_distribution<%1%>::inverse_normal_distribution";
+
+     RealType result;
+     detail::check_scale(function, sd, &result, Policy());
+     detail::check_location(function, mean, &result, Policy());
+   }
+
+   RealType mean()const
+   { // alias for location.
+      return m_mean; // aka mu
+   }
+
+   RealType standard_deviation()const
+   { // alias for scale.
+      return m_sd; // aka lambda.
+   }
+
+   // Synonyms, provided to allow generic use of find_location and find_scale.
+   RealType location()const
+   { // location, aka mu.
+      return m_mean;
+   }
+   RealType scale()const
+   { // scale, aka lambda.
+      return m_sd;
+   }
+
+private:
+   //
+   // Data members:
+   //
+   RealType m_mean;  // distribution mean or location, aka mu.
+   RealType m_sd;    // distribution standard deviation or scale, aka lambda.
+}; // class normal_distribution
+
+typedef inverse_normal_distribution<double> inverse_normal;
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> range(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{ // Range of permissible values for random variable x, zero to max.
+   using boost::math::tools::max_value;
+   return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> support(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{ // Range of supported values for random variable x, zero to max.
+   // 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>(static_cast<RealType>(0),  max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline RealType pdf(const inverse_normal_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Probability Density Function
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType scale = dist.scale();
+   RealType mean = dist.mean();
+   RealType result;
+   static const char* function = "boost::math::pdf(const inverse_normal_distribution<%1%>&, %1%)";
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_location(function, mean, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_x_gt0(function, x, &result, Policy()))
+   {
+      return numeric_limits<RealType>::quiet_NaN();
+   }
+
+   //result =
+   //  sqrt(scale / (2 * constants::pi<RealType>() * x * x * x))
+   // * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
+
+   result =
+     sqrt(scale / (constants::two_pi<RealType>() * x * x * x))
+    * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
+   return result;
+} // pdf
+
+template <class RealType, class Policy>
+inline RealType cdf(const inverse_normal_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Cumulative Density Function.
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType scale = dist.scale();
+   RealType mean = dist.mean();
+   static const char* function = "boost::math::cdf(const inverse_normal_distribution<%1%>&, %1%)";
+   RealType result;
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_location(function, mean, &result, Policy()))
+   {
+      return result;
+   }
+   if(false == detail::check_x_gt0(function, x, &result, Policy()))
+   {
+     return result;
+   }
+
+   //result = 0.5 * (erf(sqrt(scale / x) * (x / mean -1) / sqrt(2.L), Policy()) + 1)
+   //  + exp(2 * scale / mean) / 2 
+   //  * (1 - erf(sqrt(scale / x) * (x / mean +1) / sqrt(2.L), Policy()));
+
+   result = 0.5 * (erf(sqrt(scale / x) * (x / mean - 1) / constants::root_two<RealType>(), Policy()) + 1)
+     + exp(2 * scale / mean) / 2 
+     * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));
+
+   return result;
+} // cdf
+
+template <class RealType, class Policy>
+inline RealType quantile(const inverse_normal_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType mean = dist.mean();
+   RealType scale = dist.scale();
+   static const char* function = "boost::math::quantile(const inverse_normal_distribution<%1%>&, %1%)";
+
+   RealType result;
+   if(false == detail::check_scale(function, scale, &result, Policy()))
+      return result;
+   if(false == detail::check_location(function, mean, &result, Policy()))
+      return result;
+   if(false == detail::check_probability(function, x, &result, Policy()))
+      return result;
+
+   cout << "x " << x << endl;
+   RealType a = sqrt(scale / x); // a scale = lambda/x
+   RealType b = x / mean; // b = x/mu
+
+   // pnorm q, mean, sd, lower.tail = true;
+
+  //double q=1.0-pnorm(+a*(b-1.0), 0, 1, true, false);
+  //double p=    pnorm(-a*(b+1.0), 0, 1, true, false);
+
+   //boost::math::normal_distribution<RealType> norm01;
+   using boost::math::normal;
+   normal norm01;
+
+   double qx = a * (b - 1.0);
+   RealType q = 1 - ((qx <= 0) ? 0 : cdf(norm01, qx));
+   cout << "a = " << a << ", b = " <<  b << ", qx = " << qx << ", pnorm= " << pnorm01(qx) << ", cdf= "  << cdf(norm01, qx) << " q = " << q << endl;
+
+   //cout << "1 - pnorm01(qx) " << 1.0 - pnorm01(qx) << endl;
+
+
+
+   RealType px =  -a * (b + 1.0);
+   RealType p =  pnorm01(px);
+   RealType cdfpx = (px <= 0) ? 0 : cdf(norm01, px);
+   cout << "-a*(b+1.0) == px = " << px <<", pnorm01(p) = " << p << ", cdfpx = " << cdfpx << endl;
+
+   if (p == 0)
+   {
+     result = q;
+   }
+   else
+   {
+     RealType r2 = 2 * scale / mean;
+     if (r2 >= numeric_limits<RealType>::max() )
+     {
+       result = numeric_limits<RealType>::quiet_NaN();
+     }
+     else
+     {
+       result = q - exp(r2) * p;
+     }
+   }
+   return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType cdf(const complemented2_type<inverse_normal_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType sd = c.dist.standard_deviation();
+   RealType mean = c.dist.mean();
+   RealType x = c.param;
+   static const char* function = "boost::math::cdf(const complement(inverse_normal_distribution<%1%>&), %1%)";
+
+   if((boost::math::isinf)(x))
+   {
+     if(x < 0) return 1; // cdf complement -infinity is unity.
+     return 0; // cdf complement +infinity is zero
+   }
+   // These produce MSVC 4127 warnings, so the above used instead.
+   //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
+   //{ // cdf complement +infinity is zero.
+   //  return 0;
+   //}
+   //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
+   //{ // cdf complement -infinity is unity.
+   //  return 1;
+   //}
+   RealType result;
+   if(false == detail::check_scale(function, sd, &result, Policy()))
+      return result;
+   if(false == detail::check_location(function, mean, &result, Policy()))
+      return result;
+   if(false == detail::check_x(function, x, &result, Policy()))
+      return result;
+
+   RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
+   result = boost::math::erfc(diff, Policy()) / 2;
+   return result;
+} // cdf complement
+
+template <class RealType, class Policy>
+inline RealType quantile(const complemented2_type<inverse_normal_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   RealType sd = c.dist.standard_deviation();
+   RealType mean = c.dist.mean();
+   static const char* function = "boost::math::quantile(const complement(inverse_normal_distribution<%1%>&), %1%)";
+   RealType result;
+   if(false == detail::check_scale(function, sd, &result, Policy()))
+      return result;
+   if(false == detail::check_location(function, mean, &result, Policy()))
+      return result;
+   RealType q = c.param;
+   if(false == detail::check_probability(function, q, &result, Policy()))
+      return result;
+   result = boost::math::erfc_inv(2 * q, Policy());
+   result *= sd * constants::root_two<RealType>();
+   result += mean;
+   return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType mean(const inverse_normal_distribution<RealType, Policy>& dist)
+{
+   return dist.mean();
+}
+
+template <class RealType, class Policy>
+inline RealType standard_deviation(const inverse_normal_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType mean = dist.mean();
+  RealType result = sqrt(mean * mean * mean / scale)
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType mode(const inverse_normal_distribution<RealType, Policy>& dist)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale)) 
+      - 3 * mean / (2 * scale));
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType skewness(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{
+   RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = 3 * sqrt(mean/scale);
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = 12 * mean / scale ;
+  return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis_excess(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{
+  RealType scale = dist.scale();
+  RealType  mean = dist.mean();
+  RealType result = 15 * mean / scale;
+  return result;
+}
+
+} // 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_INVERSE_NORMAL_HPP
+
+
diff --git a/include/boost/math/distributions/inverse_uniform.hpp b/include/boost/math/distributions/inverse_uniform.hpp
new file mode 100644
index 000000000..7ab097a6f
--- /dev/null
+++ b/include/boost/math/distributions/inverse_uniform.hpp
@@ -0,0 +1,400 @@
+//  Copyright John Maddock 2010.
+//  Copyright Paul A. Bristow 2010.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_STATS_INVERSE_UNIFORM_HPP
+#define BOOST_STATS_INVERSE_UNIFORM_HPP
+
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm
+// http://mathworld.wolfram.com/UniformDistribution.html
+// http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/UniformDistribution.html
+// http://en.wikipedia.org/wiki/Uniform_distribution_%28continuous%29
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/complement.hpp>
+
+#include <utility>
+
+namespace boost{ namespace math
+{
+  namespace detail
+  {
+    template <class RealType, class Policy>
+    inline bool check_inverse_uniform_lower(
+      const char* function,
+      RealType lower,
+      RealType* result, const Policy& pol)
+    {
+      if((boost::math::isfinite)(lower))
+      { // any finite value is OK.
+        return true;
+      }
+      else
+      { // Not finite.
+        *result = policies::raise_domain_error<RealType>(
+          function,
+          "Lower parameter is %1%, but must be >= 0!", lower, pol);
+        return false;
+      }
+    } // bool check_inverse_uniform_lower(
+
+    template <class RealType, class Policy>
+    inline bool check_inverse_uniform_upper(
+      const char* function,
+      RealType upper,
+      RealType* result, const Policy& pol)
+    {
+      if((boost::math::isfinite)(upper))
+      { // Any finite value is OK.
+        return true;
+      }
+      else
+      { // Not finite.
+        *result = policies::raise_domain_error<RealType>(
+          function,
+          "Upper parameter is %1%, but must be finite!", upper, pol);
+        return false;
+      }
+    } // bool check_inverse_uniform_upper(
+
+    template <class RealType, class Policy>
+    inline bool check_inverse_uniform_x(
+      const char* function,
+      RealType const& x,
+      RealType* result, const Policy& pol)
+    {
+      if((boost::math::isfinite)(x))
+      { // Any finite value - if < lower or >upper will return NaN
+        return true;
+      }
+      else
+      { // Not finite..
+        *result = policies::raise_domain_error<RealType>(
+          function,
+          "y parameter is %1%, but must be finite!", x, pol);
+        return false;
+      }
+    } // bool check_inverse_uniform_x
+
+    template <class RealType, class Policy>
+    inline bool check_inverse_uniform(
+      const char* function,
+      RealType lower,
+      RealType upper,
+      RealType* result, const Policy& pol)
+    {
+      if((check_inverse_uniform_lower(function, lower, result, pol) == false)
+        || (check_inverse_uniform_upper(function, upper, result, pol) == false))
+      {
+        return false;
+      }
+      else if (lower >= upper) // If lower == upper then 1 / (upper-lower) = 1/0 = +infinity!
+      { // upper and lower have been checked before, so must be lower >= upper.
+        *result = policies::raise_domain_error<RealType>(
+          function,
+          "lower parameter is %1%, but must be less than upper!", lower, pol);
+        return false;
+      }
+      else
+      { // All OK,
+        return true;
+      }
+    } // bool check_inverse_uniform(
+
+  } // namespace detail
+
+  template <class RealType = double, class Policy = policies::policy<> >
+  class inverse_uniform_distribution
+  {
+  public:
+    typedef RealType value_type;
+    typedef Policy policy_type;
+
+    inverse_uniform_distribution(RealType lower = 0, RealType upper = 1) // Constructor.
+      : m_lower(lower), m_upper(upper) // Default is standard uniform distribution.
+    {
+      RealType result;
+      detail::check_inverse_uniform(
+        "boost::math::inverse_uniform_distribution<%1%>::inverse_uniform_distribution",
+        lower, upper, &result, Policy());
+    }
+    // Accessor functions.
+    RealType lower()const
+    {
+      return m_lower;
+    }
+
+    RealType upper()const
+    {
+      return m_upper;
+    }
+  private:
+    // Data members:
+    RealType m_lower;  // distribution lower aka a.
+    RealType m_upper;  // distribution upper aka b.
+  }; // class inverse_uniform_distribution
+
+  typedef inverse_uniform_distribution<double> inverse_uniform;
+
+  template <class RealType, class Policy>
+  inline const std::pair<RealType, RealType> range(const inverse_uniform_distribution<RealType, Policy>& /* dist */)
+  { // Range of permissible values for random variable x.
+     using boost::math::tools::max_value;
+     return std::pair<RealType, RealType>(dist.lower(), dist.upper()); // 0 to 1.
+     // Note RealType infinity is NOT permitted, only max_value.
+  }
+
+  template <class RealType, class Policy>
+  inline const std::pair<RealType, RealType> support(const inverse_uniform_distribution<RealType, Policy>& dist)
+  { // Range of supported values for random variable x.
+    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
+     using boost::math::tools::max_value;
+     return std::pair<RealType, RealType>(dist.lower(),  dist.upper());
+  }
+
+  template <class RealType, class Policy>
+  inline RealType pdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform(
+      "boost::math::pdf(const inverse_uniform_distribution<%1%>&, %1%)",
+      lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    if(false == detail::check_inverse_uniform_x(
+      "boost::math::pdf(const inverse_uniform_distribution<%1%>&, %1%)", x, &result, Policy()))
+    {
+      return result;
+    }
+    // Undefined (singularity) outside lower to upper.
+    if((x < lower) || (x > upper) )
+    {
+      return std::numeric_limits<RealType>::quiet_NaN();
+    }
+    else
+    {
+      return 1 / (upper - lower);
+    }
+  } // RealType pdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+
+  template <class RealType, class Policy>
+  inline RealType cdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform(
+      "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+      lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    if(false == detail::check_inverse_uniform_x(
+      "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+      x, &result, Policy()))
+    {
+      return result;
+    }
+    // Undefined (singularity) outside 0 to 1.
+    if (x < 0)
+    {
+      return std::numeric_limits<RealType>::quiet_NaN();
+    }
+    if (x > 1)
+    {
+      return std::numeric_limits<RealType>::quiet_NaN();
+    }
+    return x * (upper - lower) + lower; // lower <= x <= upper
+  } // RealType cdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+
+  template <class RealType, class Policy>
+  inline RealType quantile(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& p)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks
+    if(false == detail::check_inverse_uniform(
+      "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+      lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    if(false == detail::check_probability(
+      "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+      p, &result, Policy()))
+    {
+      return result;
+    }
+    if(p == 0)
+    {
+      return lower;
+    }
+    if(p == 1)
+    {
+      return upper;
+    }
+    return p * (upper - lower) + lower;
+  } // RealType quantile(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& p)
+
+  template <class RealType, class Policy>
+  inline RealType cdf(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+  {
+    RealType lower = c.dist.lower();
+    RealType upper = c.dist.upper();
+    RealType x = c.param;
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform(
+      "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+      lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    if(false == detail::check_inverse_uniform_x(
+      "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+      x, &result, Policy()))
+    {
+      return result;
+    }
+    if (x < lower)
+    {
+      return 0;
+    }
+    if (x > upper)
+    {
+      return 1;
+    }
+    return (upper - x) / (upper - lower);
+  } // RealType cdf(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+
+  template <class RealType, class Policy>
+  inline RealType quantile(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+  {
+    RealType lower = c.dist.lower();
+    RealType upper = c.dist.upper();
+    RealType q = c.param;
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform(
+      "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+      lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    if(false == detail::check_probability(
+      "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+      q, &result, Policy()))
+      if(q == 0)
+      {
+        return lower;
+      }
+      if(q == 1)
+      {
+        return upper;
+      }
+      return -q * (upper - lower) + upper;
+  } // RealType quantile(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+
+  template <class RealType, class Policy>
+  inline RealType mean(const inverse_uniform_distribution<RealType, Policy>& dist)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result;  // of checks.
+    if(false == detail::check_inverse_uniform(
+      "boost::math::mean(const inverse_uniform_distribution<%1%>&)",
+      lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    return (lower + upper ) / 2;
+  } // RealType mean(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+  template <class RealType, class Policy>
+  inline RealType variance(const inverse_uniform_distribution<RealType, Policy>& dist)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform("boost::math::variance(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    return (upper - lower) * ( upper - lower) / 12;
+    // for standard inverse_uniform = 0.833333333333333333333333333333333333333333;
+  } // RealType variance(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+  template <class RealType, class Policy>
+  inline RealType mode(const inverse_uniform_distribution<RealType, Policy>& dist)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform("boost::math::mode(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    result = lower; // Any value [lower, upper] but arbitrarily choose lower.
+    return result;
+  }
+
+  template <class RealType, class Policy>
+  inline RealType median(const inverse_uniform_distribution<RealType, Policy>& dist)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform("boost::math::median(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    return (lower + upper) / 2; //
+  }
+  template <class RealType, class Policy>
+  inline RealType skewness(const inverse_uniform_distribution<RealType, Policy>& dist)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result; // of checks.
+    if(false == detail::check_inverse_uniform("boost::math::skewness(const inverse_uniform_distribution<%1%>&)",lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    return 0;
+  } // RealType skewness(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+  template <class RealType, class Policy>
+  inline RealType kurtosis_excess(const inverse_uniform_distribution<RealType, Policy>& dist)
+  {
+    RealType lower = dist.lower();
+    RealType upper = dist.upper();
+    RealType result;  // of checks.
+    if(false == detail::check_inverse_uniform("boost::math::kurtosis_execess(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+    {
+      return result;
+    }
+    return static_cast<RealType>(-6)/5; //  -6/5 = -1.2;
+  } // RealType kurtosis_excess(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+  template <class RealType, class Policy>
+  inline RealType kurtosis(const inverse_uniform_distribution<RealType, Policy>& 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_INVERSE_UNIFORM_HPP
+
+
+