From f4781481d28e7a20280b2f833338ae143e52533a Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sun, 29 Jan 2023 10:45:12 -0800
Subject: [PATCH 01/19] Add generic logcdf

---
 .../boost/math/distributions/detail/derived_accessors.hpp  | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/include/boost/math/distributions/detail/derived_accessors.hpp b/include/boost/math/distributions/detail/derived_accessors.hpp
index d278e78776..23419b021d 100644
--- a/include/boost/math/distributions/detail/derived_accessors.hpp
+++ b/include/boost/math/distributions/detail/derived_accessors.hpp
@@ -122,6 +122,13 @@ inline typename Distribution::value_type cdf(const Distribution& dist, const Rea
    typedef typename Distribution::value_type value_type;
    return cdf(dist, static_cast<value_type>(x));
 }
+template <class Distribution, class Realtype>
+inline typename Distribution::value_type logcdf(const Distribution& dist, const Realtype& x)
+{
+   using std::log;
+   using value_type = typename Distribution::value_type;
+   return log(cdf(dist, static_cast<value_type>(x)));
+}
 template <class Distribution, class RealType>
 inline typename Distribution::value_type quantile(const Distribution& dist, const RealType& x)
 {

From 29d8fb19c3c0207fabfcc8b999f41818077778b2 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sun, 5 Feb 2023 17:31:21 -0800
Subject: [PATCH 02/19] Add logcdf specialization to extreme value distribution

---
 .../math/distributions/extreme_value.hpp      | 26 +++++++++++++++++++
 test/test_extreme_value.cpp                   | 17 +++++++++++-
 2 files changed, 42 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp
index b2e8bea966..ac06ec37c2 100644
--- a/include/boost/math/distributions/extreme_value.hpp
+++ b/include/boost/math/distributions/extreme_value.hpp
@@ -174,6 +174,32 @@ inline RealType cdf(const extreme_value_distribution<RealType, Policy>& dist, co
    return result;
 } // cdf
 
+template <class RealType, class Policy>
+inline RealType logcdf(const extreme_value_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
+
+   if((boost::math::isinf)(x))
+      return x < 0 ? 0.0f : 1.0f;
+   RealType a = dist.location();
+   RealType b = dist.scale();
+   RealType result = 0;
+   if(0 == detail::verify_scale_b(function, b, &result, Policy()))
+      return result;
+   if(0 == detail::check_finite(function, a, &result, Policy()))
+      return result;
+   if(0 == detail::check_finite(function, a, &result, Policy()))
+      return result;
+   if(0 == detail::check_x("boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)", x, &result, Policy()))
+      return result;
+
+   result = -exp((a-x)/b);
+
+   return result;
+} // logcdf
+
 template <class RealType, class Policy>
 RealType quantile(const extreme_value_distribution<RealType, Policy>& dist, const RealType& p)
 {
diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp
index 2493f466df..7e83067429 100644
--- a/test/test_extreme_value.cpp
+++ b/test/test_extreme_value.cpp
@@ -23,7 +23,7 @@
    using std::setprecision;
 
 template <class RealType>
-void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealType tolerance)
+void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealType logp, RealType tolerance)
 {
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
@@ -37,6 +37,12 @@ void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealT
          x)),
          q,
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         extreme_value_distribution<RealType>(a, b),      
+         x),
+         logp,
+         tolerance); // %
    if((p < 0.999) && (p > 0))
    {
       BOOST_CHECK_CLOSE(
@@ -78,6 +84,7 @@ void test_spots(RealType)
       static_cast<RealType>(0.125), // x
       static_cast<RealType>(0.27692033409990891617007608217222L), // p
       static_cast<RealType>(0.72307966590009108382992391782778L), //q
+      static_cast<RealType>(-1.2840254166877414840734205680624364583362808652814L), // Log(p)
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
@@ -85,6 +92,7 @@ void test_spots(RealType)
       static_cast<RealType>(-5), // x
       static_cast<RealType>(1.6087601139887776413169427645933e-7L), // p
       static_cast<RealType>(0.99999983912398860112223586830572L), //q
+      static_cast<RealType>(-15.6426318841881716102126980461566588450380350341076L), // Log(p)
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
@@ -92,6 +100,7 @@ void test_spots(RealType)
       static_cast<RealType>(0.75), // x
       static_cast<RealType>(0.69220062755534635386542199718279L), // p
       static_cast<RealType>(0.30779937244465364613457800281721), //q
+      static_cast<RealType>(-0.36787944117144232159552377016146086744581113103177L), // Log(p)
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
@@ -99,6 +108,7 @@ void test_spots(RealType)
       static_cast<RealType>(5), // x
       static_cast<RealType>(0.99999998477002037126351248727041L), // p
       static_cast<RealType>(1.5229979628736487512729586276294e-8L), //q
+      static_cast<RealType>(-1.52299797447126284361366292335174318621748e-8L), // Log(p)
       tolerance);
 
    BOOST_CHECK_CLOSE(
@@ -195,6 +205,9 @@ void test_spots(RealType)
    BOOST_MATH_CHECK_THROW(
        cdf(extreme_value_distribution<RealType>(0, -1), RealType(1)),
        std::domain_error);
+   BOOST_MATH_CHECK_THROW(
+       logcdf(extreme_value_distribution<RealType>(0, -1), RealType(1)),
+       std::domain_error);
    BOOST_MATH_CHECK_THROW(
        quantile(dist, RealType(-1)),
        std::domain_error);
@@ -211,6 +224,8 @@ void test_spots(RealType)
       BOOST_CHECK_EQUAL(cdf(extreme_value_distribution<RealType>(), inf), 1);
       BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution<RealType>(), -inf)), 1);
       BOOST_CHECK_EQUAL(cdf(complement(extreme_value_distribution<RealType>(), inf)), 0);
+      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), -inf), 0);
+      BOOST_CHECK_EQUAL(logcdf(extreme_value_distribution<RealType>(), inf), 1);
    }
    //
    // Bug reports:

From 04501b8b2fd8c892c6f5ed11bb35c936bdb4a452 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sun, 5 Feb 2023 20:04:50 -0800
Subject: [PATCH 03/19] Add logcdf specialization for laplace distribution

---
 include/boost/math/distributions/laplace.hpp | 44 ++++++++++++++++++
 test/test_laplace.cpp                        | 48 +++++++++++++++++++-
 2 files changed, 91 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp
index cc922a879a..ad14d9e95e 100644
--- a/include/boost/math/distributions/laplace.hpp
+++ b/include/boost/math/distributions/laplace.hpp
@@ -226,6 +226,50 @@ inline RealType cdf(const laplace_distribution<RealType, Policy>& dist, const Re
    return result;
 } // cdf
 
+template <class RealType, class Policy>
+inline RealType logcdf(const laplace_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING  // For ADL of std functions.
+
+   RealType result = 0;
+   // Checking function argument.
+   const char* function = "boost::math::logcdf(const laplace_distribution<%1%>&, %1%)";
+   // Check scale and location.
+   if (false == dist.check_parameters(function, &result)) 
+   {
+      return result;
+   }
+
+   // Special cdf values:
+   if((boost::math::isinf)(x))
+   {
+      if(x < 0) 
+      {
+         return 0; // -infinity.
+      }
+      return 1; // + infinity.
+   }
+
+   if (false == detail::check_x(function, x, &result, Policy())) 
+   {
+      return result;
+   }
+
+   // General cdf  values
+   RealType scale( dist.scale() );
+   RealType location( dist.location() );
+
+   if (x < location)
+   {
+      result = ((x - location) / scale) - boost::math::constants::ln_two<RealType>();
+   }
+   else
+   {
+      result = log(2 - exp((location - x) / scale)) - boost::math::constants::ln_two<RealType>();
+   }
+
+   return result;
+} // logcdf
 
 template <class RealType, class Policy>
 inline RealType quantile(const laplace_distribution<RealType, Policy>& dist, const RealType& p)
diff --git a/test/test_laplace.cpp b/test/test_laplace.cpp
index 716c738736..fc1d175593 100644
--- a/test/test_laplace.cpp
+++ b/test/test_laplace.cpp
@@ -225,6 +225,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.067667641618306345946999747486242201703815773119812L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(-2.L)),
+      //static_cast<RealType>(0.06766764161831L),
+      static_cast<RealType>(-2.6931471805599453094172321214581765680755001319215L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
       //static_cast<RealType>(0.18393972058572L),
@@ -239,7 +245,14 @@ void test_pdf_cdf_ocatave()
 
    BOOST_CHECK_CLOSE(
       cdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
-      static_cast<RealType>(0.18393972058572L),
+      //static_cast<RealType>(0.18393972058572L),
+      static_cast<RealType>(0.18393972058572116079776188508073043372290556554506L),
+   tolerance);
+
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
+      //log(static_cast<RealType>(0.18393972058572L)),
+      static_cast<RealType>(-1.6931471805599453094172321214581765680755001342016L),
    tolerance);
 
    BOOST_CHECK_CLOSE(
@@ -254,6 +267,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.30326532985631671180189976749559022672095906778368L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(-0.5L)),
+      //static_cast<RealType>(0.30326532985632L),
+      static_cast<RealType>(-1.1931471805599453094172321214581765680755001342281L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(0.0L)),
       static_cast<RealType>(0.5L),
@@ -269,6 +288,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.5L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(0.0L)),
+      //log(static_cast<RealType>(0.5L)),
+      static_cast<RealType>(-0.69314718055994530941723212145817656807550013436026L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(0.5L)),
       //static_cast<RealType>(0.30326532985632L),
@@ -286,6 +311,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.69673467014368328819810023250440977327904093221632L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(0.5L)),
+      // static_cast<RealType>(0.69673467014368L),
+      static_cast<RealType>(-0.36135061480875908299180460835979393734587166765723L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(1.0L)),
     //  static_cast<RealType>(0.18393972058572L),
@@ -304,6 +335,14 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.81606027941427883920223811491926956627709443427977L),
    tolerance);
 
+   // Tolerance of naive implementation requires tol * 1e6
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(1.00000000000000L)),
+     // static_cast<RealType>(0.81606027941428L),
+      static_cast<RealType>(-0.2032670549151953327040710243903843274153282048153L),
+      //log(static_cast<RealType>(0.81606027941427883920223811491926956627709443427977L)),
+   tolerance*1e5);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(2.0L)),
     //  static_cast<RealType>(0.06766764161831L),
@@ -322,6 +361,13 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(0.93233235838169365405300025251375779829618422688019L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(laplace_distribution<RealType>(), static_cast<RealType>(2.0L)),
+   //   static_cast<RealType>(0.93233235838169L),
+      static_cast<RealType>(-0.07006592016028141016812825150355188141186366978811L),
+      //log(static_cast<RealType>(0.93233235838169365405300025251375779829618422688019L)),
+   tolerance*1e6);
+
    check_out_of_range<laplace_distribution<RealType> >(0, 1);
    BOOST_MATH_CHECK_THROW(laplace_distribution<RealType>(0, 0), std::domain_error);
    BOOST_MATH_CHECK_THROW(laplace_distribution<RealType>(0, -1), std::domain_error);

From c73c96ce955d81b36c8ee586ffa5983f6535bca2 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 6 Feb 2023 10:40:53 -0800
Subject: [PATCH 04/19] Add logcdf specialization to logistic distribution

---
 include/boost/math/distributions/logistic.hpp | 43 +++++++++++++++
 test/test_logistic_dist.cpp                   | 52 ++++++++++++++-----
 2 files changed, 81 insertions(+), 14 deletions(-)

diff --git a/include/boost/math/distributions/logistic.hpp b/include/boost/math/distributions/logistic.hpp
index 3c7d1289d1..776c4072ff 100644
--- a/include/boost/math/distributions/logistic.hpp
+++ b/include/boost/math/distributions/logistic.hpp
@@ -146,6 +146,49 @@ namespace boost { namespace math {
        if(power < -tools::log_max_value<RealType>())
           return 1;
        return 1 / (1 + exp(power)); 
+    }
+
+    template <class RealType, class Policy>
+    inline RealType logcdf(const logistic_distribution<RealType, Policy>& dist, const RealType& x)
+    {
+       RealType scale = dist.scale();
+       RealType location = dist.location();
+       RealType result = 0; // of checks.
+       static const char* function = "boost::math::cdf(const logistic_distribution<%1%>&, %1%)";
+       if(false == detail::check_scale(function, scale, &result, Policy()))
+       {
+          return result;
+       }
+       if(false == detail::check_location(function, location, &result, Policy()))
+       {
+          return result;
+       }
+
+       if((boost::math::isinf)(x))
+       {
+          if(x < 0) 
+          {
+            return 0; // -infinity
+          }
+          return 1; // + infinity
+       }
+
+       if(false == detail::check_x(function, x, &result, Policy()))
+       {
+          return result;
+       }
+       BOOST_MATH_STD_USING
+       RealType power = (location - x) / scale;
+       if(power > tools::log_max_value<RealType>())
+       {
+          return 0;
+       }
+       if(power < -tools::log_max_value<RealType>())
+       {
+          return 1;
+       }
+
+       return -log(1 + exp(power));
     } 
     
     template <class RealType, class Policy>
diff --git a/test/test_logistic_dist.cpp b/test/test_logistic_dist.cpp
index cb6673315a..93cd95fbbd 100644
--- a/test/test_logistic_dist.cpp
+++ b/test/test_logistic_dist.cpp
@@ -32,7 +32,7 @@
    using std::setprecision;
 
 template <class RealType>
-void test_spot(RealType location, RealType scale, RealType x, RealType p, RealType q, RealType tolerance)
+void test_spot(RealType location, RealType scale, RealType x, RealType p, RealType q, RealType logp, RealType tolerance, RealType logtolerance)
 {
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
@@ -40,6 +40,12 @@ void test_spot(RealType location, RealType scale, RealType x, RealType p, RealTy
       x),
       p,
       tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+      logistic_distribution<RealType>(location,scale),      
+      x),
+      logp,
+      logtolerance); // %
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
       complement(logistic_distribution<RealType>(location,scale),      
@@ -93,23 +99,30 @@ void test_spots(RealType T)
       static_cast<RealType>(0.1L), // x
       static_cast<RealType>(0.141851064900487789594278108470953L), // p
       static_cast<RealType>(0.858148935099512210405721891529047L), //q
+      static_cast<RealType>(-1.95297761052607413459937686381969L), // log(p)
+      tolerance,
       tolerance);
-
+   
    test_spot(
       static_cast<RealType>(5), // location
       static_cast<RealType>(2), // scale
       static_cast<RealType>(3.123123123L),//x 
       static_cast<RealType>(0.281215878622547904873088053477813L), // p
       static_cast<RealType>(0.718784121377452095126911946522187L), //q
+      static_cast<RealType>(-1.26863265327659162571930803425160L),
+      tolerance,
       tolerance);
-
+   
    test_spot(
       static_cast<RealType>(1.2345L), // location
       static_cast<RealType>(0.12345L), // scale
       static_cast<RealType>(3.123123123L),//x
       static_cast<RealType>(0.999999773084685079723328282229357L), // p
       static_cast<RealType>(2.26915314920276671717770643005212e-7L), //q
-      tolerance);
+      //static_cast<RealType>(-2.2691534066556063906815875163479345e-7L), // Log(p) Edge case
+      log(static_cast<RealType>(0.999999773084685079723328282229357L)),
+      tolerance,
+      tolerance * 100000);
 
 
    //High probability
@@ -119,7 +132,9 @@ void test_spots(RealType T)
       static_cast<RealType>(10), // x
       static_cast<RealType>(0.99999998477002048723965105559179L), // p  
       static_cast<RealType>(1.5229979512760348944408208801237e-8L), //q
-      tolerance);
+      log(static_cast<RealType>(0.99999998477002048723965105559179L)),
+      tolerance,
+      tolerance * 100000);
 
    //negative x
    test_spot(
@@ -128,18 +143,20 @@ void test_spots(RealType T)
       static_cast<RealType>(-0.1L), // scale
       static_cast<RealType>(0.0724264853615177178439235061476928L), // p
       static_cast<RealType>(0.927573514638482282156076493852307L), //q
+      static_cast<RealType>(-2.6251832265757900570194868863670116L), // log(p)
+      tolerance,
       tolerance);
 
-
    test_spot(
       static_cast<RealType>(5), // location
       static_cast<RealType>(2), // scale
       static_cast<RealType>(-20), // x
       static_cast<RealType>(3.72663928418656138608800947863869e-6L), // p
       static_cast<RealType>(0.999996273360715813438613911990521L), //q
+      static_cast<RealType>(-12.500003726646228123990312715291263L),
+      tolerance,
       tolerance);
 
-
    // Test value to check cancellation error in straight/complemented quantile.
    // The subtraction in the formula location-scale*log term introduces catastrophic
    // cancellation error if location and scale*log term are close.
@@ -150,7 +167,9 @@ void test_spots(RealType T)
       static_cast<RealType>(-0.00125796420642514024493852425918807L),// x
       static_cast<RealType>(0.7L), // p
       static_cast<RealType>(0.3L), //q
-      80*tolerance);   
+      static_cast<RealType>(-0.35667494393873237891263871124118448L), // Log(p)
+      80*tolerance,
+      tolerance);   
 
    test_spot(
       static_cast<RealType>(1.2345L), // location
@@ -158,17 +177,19 @@ void test_spots(RealType T)
       static_cast<RealType>(0.0012345L), // x
       static_cast<RealType>(0.0000458541039469413343331170952855318L), // p
       static_cast<RealType>(0.999954145896053058665666882904714L), //q
-      80*tolerance);
-
-
-
+      static_cast<RealType>(-9.9900458551552785044234447480622748L), // Log(p)
+      80*tolerance,
+      tolerance);
+   
    test_spot(
       static_cast<RealType>(5L), // location
       static_cast<RealType>(2L), // scale
       static_cast<RealType>(0.0012345L), // x
       static_cast<RealType>(0.0759014628704232983512906076564256L), // p
       static_cast<RealType>(0.924098537129576701648709392343574L), //q
-      80*tolerance);
+      static_cast<RealType>(-2.5783193211111713937468666265459702L), // Log(p)
+      80*tolerance,
+      tolerance);
 
    //negative location
    test_spot(
@@ -177,7 +198,10 @@ void test_spots(RealType T)
       static_cast<RealType>(3), // x
       static_cast<RealType>(0.999999999999999999999999984171276L), // p
       static_cast<RealType>(1.58287236765203121622150720373972e-26L), //q
-      tolerance);
+      log(static_cast<RealType>(0.999999999999999999999999984171276L)), // Log(p)
+      tolerance,
+      tolerance * 100000);
+
    //PDF Testing
    BOOST_CHECK_CLOSE(
       ::boost::math::pdf(

From 58006b2250635ebdc53e0f51acb70293a5139527 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 6 Feb 2023 11:12:51 -0800
Subject: [PATCH 05/19] Add generic logcdf for complemented two type

---
 .../boost/math/distributions/detail/derived_accessors.hpp | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/include/boost/math/distributions/detail/derived_accessors.hpp b/include/boost/math/distributions/detail/derived_accessors.hpp
index 23419b021d..eb76409a1c 100644
--- a/include/boost/math/distributions/detail/derived_accessors.hpp
+++ b/include/boost/math/distributions/detail/derived_accessors.hpp
@@ -150,6 +150,14 @@ inline typename Distribution::value_type cdf(const complemented2_type<Distributi
    return cdf(complement(c.dist, static_cast<value_type>(c.param)));
 }
 
+template <class Distribution, class RealType>
+inline typename Distribution::value_type logcdf(const complemented2_type<Distribution, RealType>& c)
+{
+   using std::log;
+   typedef typename Distribution::value_type value_type;
+   return log(cdf(complement(c.dist, static_cast<value_type>(c.param))));
+}
+
 template <class Distribution, class RealType>
 inline typename Distribution::value_type quantile(const complemented2_type<Distribution, RealType>& c)
 {

From f7a4adc0ccb820d0506404f2b8b00ca5b27f0191 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 6 Feb 2023 11:24:09 -0800
Subject: [PATCH 06/19] Add logcdf specialization to weibull distribution

---
 include/boost/math/distributions/weibull.hpp | 42 ++++++++++++++++++++
 test/test_weibull.cpp                        | 13 ++++++
 2 files changed, 55 insertions(+)

diff --git a/include/boost/math/distributions/weibull.hpp b/include/boost/math/distributions/weibull.hpp
index 724cce6ed8..ca4bbd7b53 100644
--- a/include/boost/math/distributions/weibull.hpp
+++ b/include/boost/math/distributions/weibull.hpp
@@ -212,6 +212,27 @@ inline RealType cdf(const weibull_distribution<RealType, Policy>& dist, const Re
    return result;
 }
 
+template <class RealType, class Policy>
+inline RealType logcdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   static const char* function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)";
+
+   RealType shape = dist.shape();
+   RealType scale = dist.scale();
+
+   RealType result = 0;
+   if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
+      return result;
+   if(false == detail::check_weibull_x(function, x, &result, Policy()))
+      return result;
+
+   result = log1p(-exp(-pow(x / scale, shape)), Policy());
+
+   return result;
+}
+
 template <class RealType, class Policy>
 inline RealType quantile(const weibull_distribution<RealType, Policy>& dist, const RealType& p)
 {
@@ -257,6 +278,27 @@ inline RealType cdf(const complemented2_type<weibull_distribution<RealType, Poli
    return result;
 }
 
+template <class RealType, class Policy>
+inline RealType logcdf(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING  // for ADL of std functions
+
+   static const char* function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)";
+
+   RealType shape = c.dist.shape();
+   RealType scale = c.dist.scale();
+
+   RealType result = 0;
+   if(false == detail::check_weibull(function, scale, shape, &result, Policy()))
+      return result;
+   if(false == detail::check_weibull_x(function, c.param, &result, Policy()))
+      return result;
+
+   result = -pow(c.param / scale, shape);
+
+   return result;
+}
+
 template <class RealType, class Policy>
 inline RealType quantile(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)
 {
diff --git a/test/test_weibull.cpp b/test/test_weibull.cpp
index 6cacf5402a..4b31a7f0b7 100644
--- a/test/test_weibull.cpp
+++ b/test/test_weibull.cpp
@@ -41,6 +41,12 @@ void check_weibull(RealType shape, RealType scale, RealType x, RealType p, RealT
          x),                                            // random variable.
          p,                                             // probability.
          tol);                                          // %tolerance.
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         weibull_distribution<RealType>(shape, scale),       // distribution.
+         x),                                            // random variable.
+         log(p),                                             // probability.
+         tol);   
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
          complement(
@@ -48,6 +54,13 @@ void check_weibull(RealType shape, RealType scale, RealType x, RealType p, RealT
             x)),                                        // random variable.
          q,                                             // probability complement.
          tol);                                          // %tolerance.
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         complement(
+            weibull_distribution<RealType>(shape, scale),    // distribution.
+            x)),                                        // random variable.
+         log(q),                                             // probability complement.
+         tol);   
    BOOST_CHECK_CLOSE(
       ::boost::math::quantile(
          weibull_distribution<RealType>(shape, scale),       // distribution.

From b02674d50fb6e4d8689a6b33a1706b24e2cca59d Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 6 Feb 2023 18:02:59 -0800
Subject: [PATCH 07/19] Improve logistic distribution logcdf

---
 include/boost/math/distributions/logistic.hpp | 2 +-
 test/test_logistic_dist.cpp                   | 7 +++----
 2 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/include/boost/math/distributions/logistic.hpp b/include/boost/math/distributions/logistic.hpp
index 776c4072ff..32076eacf4 100644
--- a/include/boost/math/distributions/logistic.hpp
+++ b/include/boost/math/distributions/logistic.hpp
@@ -188,7 +188,7 @@ namespace boost { namespace math {
           return 1;
        }
 
-       return -log(1 + exp(power));
+       return -log1p(exp(power));
     } 
     
     template <class RealType, class Policy>
diff --git a/test/test_logistic_dist.cpp b/test/test_logistic_dist.cpp
index 93cd95fbbd..737e868cc9 100644
--- a/test/test_logistic_dist.cpp
+++ b/test/test_logistic_dist.cpp
@@ -119,8 +119,7 @@ void test_spots(RealType T)
       static_cast<RealType>(3.123123123L),//x
       static_cast<RealType>(0.999999773084685079723328282229357L), // p
       static_cast<RealType>(2.26915314920276671717770643005212e-7L), //q
-      //static_cast<RealType>(-2.2691534066556063906815875163479345e-7L), // Log(p) Edge case
-      log(static_cast<RealType>(0.999999773084685079723328282229357L)),
+      static_cast<RealType>(-2.2691534066556063906815875163479345e-7L), // Log(p) Edge case
       tolerance,
       tolerance * 100000);
 
@@ -132,7 +131,7 @@ void test_spots(RealType T)
       static_cast<RealType>(10), // x
       static_cast<RealType>(0.99999998477002048723965105559179L), // p  
       static_cast<RealType>(1.5229979512760348944408208801237e-8L), //q
-      log(static_cast<RealType>(0.99999998477002048723965105559179L)),
+      static_cast<RealType>(-1.5229979628736488101501003766470705e-8),
       tolerance,
       tolerance * 100000);
 
@@ -198,7 +197,7 @@ void test_spots(RealType T)
       static_cast<RealType>(3), // x
       static_cast<RealType>(0.999999999999999999999999984171276L), // p
       static_cast<RealType>(1.58287236765203121622150720373972e-26L), //q
-      log(static_cast<RealType>(0.999999999999999999999999984171276L)), // Log(p)
+      static_cast<RealType>(-1.58287236765203121622150720373972e-26L), // Log(p)
       tolerance,
       tolerance * 100000);
 

From c0f89544c3a19f465e17d0644f17fc3e899d9801 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Mon, 6 Feb 2023 18:14:49 -0800
Subject: [PATCH 08/19] Add logistic distribution logcdf for complement two
 type

---
 include/boost/math/distributions/logistic.hpp | 36 +++++++++++++++
 test/test_logistic_dist.cpp                   | 46 ++++++++++++++++---
 2 files changed, 75 insertions(+), 7 deletions(-)

diff --git a/include/boost/math/distributions/logistic.hpp b/include/boost/math/distributions/logistic.hpp
index 32076eacf4..d12de48c59 100644
--- a/include/boost/math/distributions/logistic.hpp
+++ b/include/boost/math/distributions/logistic.hpp
@@ -262,6 +262,42 @@ namespace boost { namespace math {
        return 1 / (1 + exp(power)); 
     } 
 
+    template <class RealType, class Policy>
+    inline RealType logcdf(const complemented2_type<logistic_distribution<RealType, Policy>, RealType>& c)
+    {
+       BOOST_MATH_STD_USING
+       RealType location = c.dist.location();
+       RealType scale = c.dist.scale();
+       RealType x = c.param;
+       static const char* function = "boost::math::cdf(const complement(logistic_distribution<%1%>&), %1%)";
+
+       RealType result = 0;
+       if(false == detail::check_scale(function, scale, &result, Policy()))
+       {
+          return result;
+       }
+       if(false == detail::check_location(function, location, &result, Policy()))
+       {
+          return result;
+       }
+       if((boost::math::isinf)(x))
+       {
+          if(x < 0) return 1; // cdf complement -infinity is unity.
+          return 0; // cdf complement +infinity is zero.
+       }
+       if(false == detail::check_x(function, x, &result, Policy()))
+       {
+          return result;
+       }
+       RealType power = (x - location) / scale;
+       if(power > tools::log_max_value<RealType>())
+          return 0;
+       if(power < -tools::log_max_value<RealType>())
+          return 1;
+       
+       return -log1p(exp(power));
+    }  
+
     template <class RealType, class Policy>
     inline RealType quantile(const complemented2_type<logistic_distribution<RealType, Policy>, RealType>& c)
     {
diff --git a/test/test_logistic_dist.cpp b/test/test_logistic_dist.cpp
index 737e868cc9..178f92d34c 100644
--- a/test/test_logistic_dist.cpp
+++ b/test/test_logistic_dist.cpp
@@ -26,14 +26,20 @@
 #include <boost/test/tools/floating_point_comparison.hpp>
 #include "test_out_of_range.hpp"
 
+#include <type_traits>
 #include <iostream>
    using std::cout;
    using std::endl;
    using std::setprecision;
 
 template <class RealType>
-void test_spot(RealType location, RealType scale, RealType x, RealType p, RealType q, RealType logp, RealType tolerance, RealType logtolerance)
+void test_spot(RealType location, RealType scale, RealType x, RealType p, RealType q, RealType logp, RealType logq, RealType tolerance, RealType logtolerance, RealType logtoleranceq)
 {
+   if (std::is_same<RealType, long double>::value)
+   {
+      logtoleranceq *= 100;
+   }
+   
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
       logistic_distribution<RealType>(location,scale),      
@@ -52,6 +58,12 @@ void test_spot(RealType location, RealType scale, RealType x, RealType p, RealTy
       x)),
       q,
       tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+      complement(logistic_distribution<RealType>(location,scale),      
+      x)),
+      log(q),
+      logtoleranceq); // %
    if(p < 0.999)
    {
       BOOST_CHECK_CLOSE(
@@ -100,6 +112,8 @@ void test_spots(RealType T)
       static_cast<RealType>(0.141851064900487789594278108470953L), // p
       static_cast<RealType>(0.858148935099512210405721891529047L), //q
       static_cast<RealType>(-1.95297761052607413459937686381969L), // log(p)
+      static_cast<RealType>(-0.15297761052607413459937686381969L), // log(q)
+      tolerance,
       tolerance,
       tolerance);
    
@@ -109,7 +123,9 @@ void test_spots(RealType T)
       static_cast<RealType>(3.123123123L),//x 
       static_cast<RealType>(0.281215878622547904873088053477813L), // p
       static_cast<RealType>(0.718784121377452095126911946522187L), //q
-      static_cast<RealType>(-1.26863265327659162571930803425160L),
+      static_cast<RealType>(-1.26863265327659162571930803425160L), // log(p)
+      static_cast<RealType>(-0.33019421477659162571930803425160L), // log(q)
+      tolerance,
       tolerance,
       tolerance);
    
@@ -120,7 +136,9 @@ void test_spots(RealType T)
       static_cast<RealType>(0.999999773084685079723328282229357L), // p
       static_cast<RealType>(2.26915314920276671717770643005212e-7L), //q
       static_cast<RealType>(-2.2691534066556063906815875163479345e-7L), // Log(p) Edge case
+      static_cast<RealType>(-15.298688951095170556204624487357L), // log(q)
       tolerance,
+      tolerance * 100000,
       tolerance * 100000);
 
 
@@ -131,8 +149,10 @@ void test_spots(RealType T)
       static_cast<RealType>(10), // x
       static_cast<RealType>(0.99999998477002048723965105559179L), // p  
       static_cast<RealType>(1.5229979512760348944408208801237e-8L), //q
-      static_cast<RealType>(-1.5229979628736488101501003766470705e-8),
+      static_cast<RealType>(-1.5229979628736488101501003766470705e-8), // log(p)
+      static_cast<RealType>(-18.000000015229979628736488101501L), // Log(q)
       tolerance,
+      tolerance * 100000,
       tolerance * 100000);
 
    //negative x
@@ -143,6 +163,8 @@ void test_spots(RealType T)
       static_cast<RealType>(0.0724264853615177178439235061476928L), // p
       static_cast<RealType>(0.927573514638482282156076493852307L), //q
       static_cast<RealType>(-2.6251832265757900570194868863670116L), // log(p)
+      static_cast<RealType>(-0.075183226575790057019486886367011L), // log(q)
+      tolerance,
       tolerance,
       tolerance);
 
@@ -152,9 +174,11 @@ void test_spots(RealType T)
       static_cast<RealType>(-20), // x
       static_cast<RealType>(3.72663928418656138608800947863869e-6L), // p
       static_cast<RealType>(0.999996273360715813438613911990521L), //q
-      static_cast<RealType>(-12.500003726646228123990312715291263L),
+      static_cast<RealType>(-12.500003726646228123990312715291263L), // log(p)
+      static_cast<RealType>(-3.7266462281239903127152912633031e-6L), // Log(q)
       tolerance,
-      tolerance);
+      tolerance,
+      tolerance * 100000);
 
    // Test value to check cancellation error in straight/complemented quantile.
    // The subtraction in the formula location-scale*log term introduces catastrophic
@@ -167,7 +191,9 @@ void test_spots(RealType T)
       static_cast<RealType>(0.7L), // p
       static_cast<RealType>(0.3L), //q
       static_cast<RealType>(-0.35667494393873237891263871124118448L), // Log(p)
+      static_cast<RealType>(-1.2039728043259359926227462177618L), // log(q)
       80*tolerance,
+      tolerance,
       tolerance);   
 
    test_spot(
@@ -177,8 +203,10 @@ void test_spots(RealType T)
       static_cast<RealType>(0.0000458541039469413343331170952855318L), // p
       static_cast<RealType>(0.999954145896053058665666882904714L), //q
       static_cast<RealType>(-9.9900458551552785044234447480622748L), // Log(p)
+      static_cast<RealType>(-0.000045855155278504423444748062274193L), // log(q)
       80*tolerance,
-      tolerance);
+      tolerance,
+      tolerance * 100000);
    
    test_spot(
       static_cast<RealType>(5L), // location
@@ -187,7 +215,9 @@ void test_spots(RealType T)
       static_cast<RealType>(0.0759014628704232983512906076564256L), // p
       static_cast<RealType>(0.924098537129576701648709392343574L), //q
       static_cast<RealType>(-2.5783193211111713937468666265459702L), // Log(p)
+      static_cast<RealType>(-0.078936571111171393746866626545970L), // log(q)
       80*tolerance,
+      tolerance,
       tolerance);
 
    //negative location
@@ -198,8 +228,10 @@ void test_spots(RealType T)
       static_cast<RealType>(0.999999999999999999999999984171276L), // p
       static_cast<RealType>(1.58287236765203121622150720373972e-26L), //q
       static_cast<RealType>(-1.58287236765203121622150720373972e-26L), // Log(p)
+      static_cast<RealType>(-59.407971267074894017899199262178L), // log(q)
       tolerance,
-      tolerance * 100000);
+      tolerance * 100000,
+      tolerance * 10000000);
 
    //PDF Testing
    BOOST_CHECK_CLOSE(

From 6159c05a2be40338d1b55dd4359530ec2672f0be Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Tue, 7 Feb 2023 10:13:22 -0800
Subject: [PATCH 09/19] Improve laplace logcdf

---
 include/boost/math/distributions/laplace.hpp | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp
index ad14d9e95e..0a03ff6d9e 100644
--- a/include/boost/math/distributions/laplace.hpp
+++ b/include/boost/math/distributions/laplace.hpp
@@ -17,6 +17,7 @@
 #ifndef BOOST_STATS_LAPLACE_HPP
 #define BOOST_STATS_LAPLACE_HPP
 
+#include <boost/math/special_functions/log1p.hpp>
 #include <boost/math/distributions/detail/common_error_handling.hpp>
 #include <boost/math/distributions/complement.hpp>
 #include <boost/math/constants/constants.hpp>
@@ -265,7 +266,7 @@ inline RealType logcdf(const laplace_distribution<RealType, Policy>& dist, const
    }
    else
    {
-      result = log(2 - exp((location - x) / scale)) - boost::math::constants::ln_two<RealType>();
+      result = log1p(-exp((location - x) / scale) / 2);
    }
 
    return result;

From a77f0fa196221b8277cfc50a65dd68d7a1b0b804 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sat, 11 Feb 2023 16:42:20 -0800
Subject: [PATCH 10/19] Add logcdf specialization to exponential distribution

---
 .../boost/math/distributions/exponential.hpp  | 39 ++++++++++++
 test/test_exponential_dist.cpp                | 62 ++++++++++++++++++-
 2 files changed, 98 insertions(+), 3 deletions(-)

diff --git a/include/boost/math/distributions/exponential.hpp b/include/boost/math/distributions/exponential.hpp
index 5214575a64..164e01f205 100644
--- a/include/boost/math/distributions/exponential.hpp
+++ b/include/boost/math/distributions/exponential.hpp
@@ -163,6 +163,24 @@ inline RealType cdf(const exponential_distribution<RealType, Policy>& dist, cons
    return result;
 } // cdf
 
+template <class RealType, class Policy>
+inline RealType logcdf(const exponential_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   static const char* function = "boost::math::logcdf(const exponential_distribution<%1%>&, %1%)";
+
+   RealType result = 0;
+   RealType lambda = dist.lambda();
+   if(0 == detail::verify_lambda(function, lambda, &result, Policy()))
+      return result;
+   if(0 == detail::verify_exp_x(function, x, &result, Policy()))
+      return result;
+   result = boost::math::log1p(-exp(-x * lambda), Policy());
+
+   return result;
+} // cdf
+
 template <class RealType, class Policy>
 inline RealType quantile(const exponential_distribution<RealType, Policy>& dist, const RealType& p)
 {
@@ -207,6 +225,27 @@ inline RealType cdf(const complemented2_type<exponential_distribution<RealType,
    return result;
 }
 
+template <class RealType, class Policy>
+inline RealType logcdf(const complemented2_type<exponential_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   static const char* function = "boost::math::logcdf(const exponential_distribution<%1%>&, %1%)";
+
+   RealType result = 0;
+   RealType lambda = c.dist.lambda();
+   if(0 == detail::verify_lambda(function, lambda, &result, Policy()))
+      return result;
+   if(0 == detail::verify_exp_x(function, c.param, &result, Policy()))
+      return result;
+   // Workaround for VC11/12 bug:
+   if (c.param >= tools::max_value<RealType>())
+      return 0;
+   result = -c.param * lambda;
+
+   return result;
+}
+
 template <class RealType, class Policy>
 inline RealType quantile(const complemented2_type<exponential_distribution<RealType, Policy>, RealType>& c)
 {
diff --git a/test/test_exponential_dist.cpp b/test/test_exponential_dist.cpp
index 12d9fcad19..c659b91c2b 100644
--- a/test/test_exponential_dist.cpp
+++ b/test/test_exponential_dist.cpp
@@ -18,14 +18,22 @@
 #include <boost/test/tools/floating_point_comparison.hpp>
 #include "test_out_of_range.hpp"
 
+#include <cmath>
+#include <type_traits>
 #include <iostream>
    using std::cout;
    using std::endl;
    using std::setprecision;
+   using std::log;
 
 template <class RealType>
-void test_spot(RealType l, RealType x, RealType p, RealType q, RealType tolerance)
+void test_spot(RealType l, RealType x, RealType p, RealType q, RealType logp, RealType logq, RealType tolerance, RealType logtolerance)
 {
+   BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value)
+   {
+      logtolerance *= 100;
+   }
+   
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
          exponential_distribution<RealType>(l),
@@ -38,6 +46,18 @@ void test_spot(RealType l, RealType x, RealType p, RealType q, RealType toleranc
          x)),
          q,
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         exponential_distribution<RealType>(l),
+         x),
+         logp,
+         logtolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         complement(exponential_distribution<RealType>(l),
+         x)),
+         logq,
+         logtolerance); // %
    if(p < 0.999)
    {
       BOOST_CHECK_CLOSE(
@@ -82,24 +102,36 @@ void test_spots(RealType T)
       static_cast<RealType>(0.125), // x
       static_cast<RealType>(0.060586937186524213880289175377695L), // p
       static_cast<RealType>(0.93941306281347578611971082462231L), //q
+      static_cast<RealType>(log(static_cast<RealType>(0.060586937186524213880289175377695L))),
+      static_cast<RealType>(log(static_cast<RealType>(0.93941306281347578611971082462231L))),
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // lambda
       static_cast<RealType>(5), // x
       static_cast<RealType>(0.91791500137610120483047132553284L), // p
       static_cast<RealType>(0.08208499862389879516952867446716L), //q
+      static_cast<RealType>(log(static_cast<RealType>(0.91791500137610120483047132553284L))),
+      static_cast<RealType>(log(static_cast<RealType>(0.08208499862389879516952867446716L))),
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(2), // lambda
       static_cast<RealType>(0.125), // x
       static_cast<RealType>(0.22119921692859513175482973302168L), // p
       static_cast<RealType>(0.77880078307140486824517026697832L), //q
+      static_cast<RealType>(log(static_cast<RealType>(0.22119921692859513175482973302168L))), // log(p)
+      static_cast<RealType>(log(static_cast<RealType>(0.77880078307140486824517026697832L))), // log(q)
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(2), // lambda
       static_cast<RealType>(5), // x
       static_cast<RealType>(0.99995460007023751514846440848444L), // p
       static_cast<RealType>(4.5399929762484851535591515560551e-5L), //q
+      static_cast<RealType>(-0.00004540096037048920950444635987890882815054L),
+      static_cast<RealType>(-10.0000000000000000000000000000000000000000000L),
+      tolerance,
       tolerance);
 
    //
@@ -110,49 +142,73 @@ void test_spots(RealType T)
       static_cast<RealType>(1), // x
       static_cast<RealType>(6.321205588285580E-001L), // p
       static_cast<RealType>(1-6.321205588285580E-001L), //q
+      static_cast<RealType>(-0.4586751453870818910216436450673297018770L), // log(p)
+      static_cast<RealType>(-1.0000000000000000000000000000000000000000L),
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(2), // lambda
       static_cast<RealType>(1), // x
       static_cast<RealType>(8.646647167633870E-001L), // p
       static_cast<RealType>(1-8.646647167633870E-001L), //q
+      static_cast<RealType>(-0.1454134578688590569726481500994740599617L), //log(p)
+      static_cast<RealType>(-2.000000000000000000000000000000000000000L),
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(1), // lambda
       static_cast<RealType>(0.5), // x
       static_cast<RealType>(3.934693402873670E-001L), // p
       static_cast<RealType>(1-3.934693402873670E-001L), //q
+      static_cast<RealType>(-0.9327521295671885718946410001485004516325L), 
+      static_cast<RealType>(-0.5000000000000000000000000000000000000000L),
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(0.1), // lambda
       static_cast<RealType>(1), // x
       static_cast<RealType>(9.516258196404040E-002L), // p
       static_cast<RealType>(1-9.516258196404040E-002L), //q
+      static_cast<RealType>(-2.352168461044090808919497187062612584712L),
+      static_cast<RealType>(-0.100000000000000000000000000000000000000L),
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(10), // lambda
       static_cast<RealType>(1), // x
       static_cast<RealType>(9.999546000702380E-001L), // p
       static_cast<RealType>(1-9.999546000702380E-001L), //q
-      tolerance*10000); // we loose four digits to cancellation
+      static_cast<RealType>(-0.00004540096037048920950444635987890882815054L),
+      static_cast<RealType>(-10.0000000000000000000000000000000000000000000L),
+      tolerance*10000, // we loose four digits to cancellation
+      tolerance); 
    test_spot(
       static_cast<RealType>(0.1), // lambda
       static_cast<RealType>(10), // x
       static_cast<RealType>(6.321205588285580E-001L), // p
       static_cast<RealType>(1-6.321205588285580E-001L), //q
+      static_cast<RealType>(-0.4586751453870818910216436450673297018770L), // log(p),
+      static_cast<RealType>(-1.0000000000000000000000000000000000000000L), // log(q)
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(1), // lambda
       static_cast<RealType>(0.01), // x
       static_cast<RealType>(9.950166250831950E-003L), // p
       static_cast<RealType>(1-9.950166250831950E-003L), //q
+      static_cast<RealType>(-4.610166019324896918080084954100938122002L), // log(p),
+      static_cast<RealType>(-0.010000000000000000000000000000000000000L), // log(q)
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(1), // lambda
       static_cast<RealType>(0.0001), // x
       static_cast<RealType>(9.999500016666250E-005L), // p
       static_cast<RealType>(1-9.999500016666250E-005L), //q
-      tolerance);
+      static_cast<RealType>(-9.210390371559516069440021374287500922116L), // log(p),
+      static_cast<RealType>(-0.000100000000000000000000000000000000000L), // log(q)
+      tolerance,
+      std::is_same<RealType, float>::value ? tolerance * 10 : tolerance);
    /*
    // This test data appears to be erroneous, MathCad appears
    // to suffer from cancellation error as x -> 0

From e387f9c429be63e397509923427701c4af74e182 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sat, 11 Feb 2023 17:28:27 -0800
Subject: [PATCH 11/19] Add logcdf specialization to geometric distribution

---
 .../boost/math/distributions/geometric.hpp    | 62 ++++++++++++++++++-
 test/test_geometric.cpp                       | 52 +++++++++++++---
 2 files changed, 103 insertions(+), 11 deletions(-)

diff --git a/include/boost/math/distributions/geometric.hpp b/include/boost/math/distributions/geometric.hpp
index b091bf1c63..baff0e249f 100644
--- a/include/boost/math/distributions/geometric.hpp
+++ b/include/boost/math/distributions/geometric.hpp
@@ -43,9 +43,11 @@
 #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/math/special_functions/log1p.hpp>
 
 #include <limits> // using std::numeric_limits;
 #include <utility>
+#include <cmath>
 
 #if defined (BOOST_MSVC)
 #  pragma warning(push)
@@ -378,9 +380,39 @@ namespace boost
       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.
+    template <class RealType, class Policy>
+    inline RealType logcdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)
+    { // Cumulative Distribution Function of geometric.
+      using std::pow;
+      static const char* function = "boost::math::logcdf(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 = 0;
+      if(false == geometric_detail::check_dist_and_k(
+        function,
+        p,
+        k,
+        &result, Policy()))
+      {
+        return -std::numeric_limits<RealType>::infinity();
+      }
+      if(k == 0)
+      {
+        return log(p); // success_fraction
+      }
+      //RealType q = 1 - p;  // Bad for small p
+      //RealType probability = 1 - std::pow(q, k+1);
+
+      RealType z = boost::math::log1p(-p, Policy()) * (k + 1);
+      return log1p(-exp(z), Policy());
+    } // logcdf 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.
       BOOST_MATH_STD_USING
       static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";
       // k argument may be integral, signed, or unsigned, or floating point.
@@ -403,6 +435,30 @@ namespace boost
       return probability;
     } // cdf Complemented Cumulative Distribution Function geometric.
 
+    template <class RealType, class Policy>
+    inline RealType logcdf(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)
+    { // Complemented Cumulative Distribution Function geometric.
+      BOOST_MATH_STD_USING
+      static const char* function = "boost::math::logcdf(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 = 0;
+      if(false == geometric_detail::check_dist_and_k(
+        function,
+        p,
+        k,
+        &result, Policy()))
+      {
+        return -std::numeric_limits<RealType>::infinity();
+      }
+
+      return boost::math::log1p(-p, Policy()) * (k+1);
+    } // logcdf 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.
diff --git a/test/test_geometric.cpp b/test/test_geometric.cpp
index 85a6e64543..ccc3673228 100644
--- a/test/test_geometric.cpp
+++ b/test/test_geometric.cpp
@@ -48,18 +48,31 @@ using std::setprecision;
 using std::showpoint;
 #include <limits>
 using std::numeric_limits;
+#include <cmath>
+using std::log;
+using std::abs;
+#include <type_traits>
 
 template <class RealType>
 void test_spot( // Test a single spot value against 'known good' values.
-               RealType k,    // Number of failures.
-               RealType p,    // Probability of success_fraction.
-               RealType P,    // CDF probability.
-               RealType Q,    // Complement of CDF.
-               RealType tol)  // Test tolerance.
+               RealType k,       // Number of failures.
+               RealType p,       // Probability of success_fraction.
+               RealType P,       // CDF probability.
+               RealType Q,       // Complement of CDF.
+               RealType logP,    // Logcdf probability
+               RealType logQ,    // Complement of logcdf
+               RealType tol,     // Test tolerance
+               RealType logtol)  // Logcdf Test tolerance.
 {
+   BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value || std::is_same<RealType, real_concept>::value)
+   {
+     logtol *= 100;
+   }
+   
    boost::math::geometric_distribution<RealType> g(p);
    BOOST_CHECK_EQUAL(p, g.success_fraction());
    BOOST_CHECK_CLOSE_FRACTION(cdf(g, k), P, tol);
+   BOOST_CHECK_CLOSE_FRACTION(logcdf(g, k), logP, logtol);
 
   if((P < 0.99) && (Q < 0.99))
   {
@@ -68,6 +81,8 @@ void test_spot( // Test a single spot value against 'known good' values.
     //
     BOOST_CHECK_CLOSE_FRACTION(
       cdf(complement(g, k)), Q, tol);
+    BOOST_CHECK_CLOSE_FRACTION(
+      logcdf(complement(g, k)), logQ, logtol);
     if(k != 0)
     {
       BOOST_CHECK_CLOSE_FRACTION(
@@ -222,6 +237,9 @@ void test_spots(RealType)
   static_cast<RealType>(0.5), // Probability of success as fraction, p
   static_cast<RealType>(0.875L), // Probability of result (CDF), P
   static_cast<RealType>(0.125L),  // complement CCDF Q = 1 - P
+  static_cast<RealType>(-0.1335313926245226231463436209313499745894L),
+  static_cast<RealType>(-2.079441541679835928251696364374529704227L),
+  tolerance,
   tolerance);
 
   test_spot( //
@@ -229,6 +247,9 @@ void test_spots(RealType)
   static_cast<RealType>(0.25), // Probability of success as fraction, p
   static_cast<RealType>(0.25),   // Probability of result (CDF), P
   static_cast<RealType>(0.75),   // Q = 1 - P
+  static_cast<RealType>(-1.386294361119890618834464242916353136151L),
+  static_cast<RealType>(-0.2876820724517809274392190059938274315035L),
+  tolerance,
   tolerance);
 
   test_spot(
@@ -239,6 +260,9 @@ void test_spots(RealType)
   static_cast<RealType>(0.25),  // Probability of success, p
   static_cast<RealType>(0.95776486396789551L),  // Probability of result (CDF), P
   static_cast<RealType>(0.042235136032104499L), // Q = 1 - P
+  static_cast<RealType>(-0.04315297584768019483875419429616349387993L),
+  static_cast<RealType>(-3.164502796969590201831409065932101746539L),
+  tolerance,
   tolerance);
 
   test_spot(  //
@@ -248,6 +272,9 @@ void test_spots(RealType)
   static_cast<RealType>(0.25),     // Probability of success, p
   static_cast<RealType>(0.99999957525875771),  // Probability of result (CDF), P
   static_cast<RealType>(4.2474124232020353e-07),   // Q = 1 - P
+  static_cast<RealType>(-4.247413325227902241937783772756893037512e-7L),
+  static_cast<RealType>(-14.67178569504082729940016930568519898711L),
+  tolerance,
   tolerance);
   /*
   // This causes failures in find_upper_bound_on_p p is small branch.
@@ -266,6 +293,9 @@ void test_spots(RealType)
   static_cast<RealType>(0.99),    // Probability of success, p
   static_cast<RealType>(1), // Probability of result (CDF), P
   static_cast<RealType>(1.0000000000000364e-102),   // Q = 1 - P
+  static_cast<RealType>(-1.0000000000000364e-102L),
+  static_cast<RealType>(-std::numeric_limits<RealType>::infinity()),
+  tolerance,
   tolerance);
 
   test_spot(  // > formatC(pgeom(1,0.99, TRUE), digits=17) [1] "0.99990000000000001"
@@ -274,7 +304,10 @@ void test_spots(RealType)
   static_cast<RealType>(0.99),                    // Probability of success, p
   static_cast<RealType>(0.9999),     // Probability of result (CDF), P
   static_cast<RealType>(0.0001),   // Q = 1 - P
-  tolerance);
+  static_cast<RealType>(-0.0001000050003333583353335000142869643968354L),
+  static_cast<RealType>(-9.210340371976182736071965818737456830404L),
+  tolerance,
+  tolerance * 100);
 
 if(std::numeric_limits<RealType>::is_specialized)
 { // An extreme value test that is more accurate than using negative binomial.
@@ -284,8 +317,11 @@ if(std::numeric_limits<RealType>::is_specialized)
   static_cast<RealType>(10000L), // Number of failures, k
   static_cast<RealType>(0.001L),                    // Probability of success, p
   static_cast<RealType>(0.99995487182736897L),     // Probability of result (CDF), P
-  static_cast<RealType>(4.5128172631071587e-05L),   // Q = 1 - P
-  tolerance); //
+  static_cast<RealType>(4.5128172631071587e-05L),   // Q = 1 - P,
+  static_cast<RealType>(-0.00004512919093769043386238651458397312570531L),
+  static_cast<RealType>(-10.00600383616891853492996552293751795172L),
+  tolerance,
+  tolerance * 100); //
   } // numeric_limit is specialized
  // End of single spot tests using RealType
 

From 614445505af70cfb8a7dbc54160251f53bbbedac Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sat, 11 Feb 2023 17:49:29 -0800
Subject: [PATCH 12/19] Add logcdf specialization to pareto distribution

---
 include/boost/math/distributions/pareto.hpp | 48 +++++++++++++++++++++
 test/test_pareto.cpp                        | 22 ++++++++++
 2 files changed, 70 insertions(+)

diff --git a/include/boost/math/distributions/pareto.hpp b/include/boost/math/distributions/pareto.hpp
index 45d43587f0..778310d10e 100644
--- a/include/boost/math/distributions/pareto.hpp
+++ b/include/boost/math/distributions/pareto.hpp
@@ -1,5 +1,6 @@
 //  Copyright John Maddock 2007.
 //  Copyright Paul A. Bristow 2007, 2009
+//  Copyright Matt Borland 2023.
 //  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)
@@ -20,8 +21,11 @@
 #include <boost/math/distributions/complement.hpp>
 #include <boost/math/distributions/detail/common_error_handling.hpp>
 #include <boost/math/special_functions/powm1.hpp>
+#include <boost/math/special_functions/log1p.hpp>
 
 #include <utility> // for BOOST_CURRENT_VALUE?
+#include <limits>
+#include <cmath>
 
 namespace boost
 {
@@ -225,6 +229,28 @@ namespace boost
       return result;
     } // cdf
 
+    template <class RealType, class Policy>
+    inline RealType logcdf(const pareto_distribution<RealType, Policy>& dist, const RealType& x)
+    {
+      BOOST_MATH_STD_USING  // for ADL of std function pow.
+      static const char* function = "boost::math::logcdf(const pareto_distribution<%1%>&, %1%)";
+      RealType scale = dist.scale();
+      RealType shape = dist.shape();
+      RealType result = 0;
+
+      if(false == (detail::check_pareto_x(function, x, &result, Policy())
+         && detail::check_pareto(function, scale, shape, &result, Policy())))
+         return result;
+
+      if (x <= scale)
+      { // regardless of shape, cdf is zero.
+        return -std::numeric_limits<RealType>::infinity();
+      }
+
+      result = log1p(-pow(scale/x, shape), Policy());
+      return result;
+    } // logcdf
+
     template <class RealType, class Policy>
     inline RealType quantile(const pareto_distribution<RealType, Policy>& dist, const RealType& p)
     {
@@ -274,6 +300,28 @@ namespace boost
        return result;
     } // cdf complement
 
+    template <class RealType, class Policy>
+    inline RealType logcdf(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
+    {
+       BOOST_MATH_STD_USING  // for ADL of std function pow.
+       static const char* function = "boost::math::logcdf(const pareto_distribution<%1%>&, %1%)";
+       RealType result = 0;
+       RealType x = c.param;
+       RealType scale = c.dist.scale();
+       RealType shape = c.dist.shape();
+       if(false == (detail::check_pareto_x(function, x, &result, Policy())
+           && detail::check_pareto(function, scale, shape, &result, Policy())))
+         return result;
+
+       if (x <= scale)
+       { // regardless of shape, cdf is zero, and complement is unity.
+         return 0;
+       }
+       result = log(pow((scale/x), shape));
+
+       return result;
+    } // logcdf complement
+
     template <class RealType, class Policy>
     inline RealType quantile(const complemented2_type<pareto_distribution<RealType, Policy>, RealType>& c)
     {
diff --git a/test/test_pareto.cpp b/test/test_pareto.cpp
index 792787d279..35a5bb0098 100644
--- a/test/test_pareto.cpp
+++ b/test/test_pareto.cpp
@@ -1,5 +1,6 @@
 // Copyright Paul A. Bristow 2007, 2009.
 // Copyright John Maddock 2006.
+// Copyright Matt Borland 2023.
 // Use, modification and distribution are subject to the
 // Boost Software License, Version 1.0.
 // (See accompanying file LICENSE_1_0.txt
@@ -40,10 +41,18 @@
    using std::setprecision;
 #include <limits>
   using std::numeric_limits;
+#include <type_traits>
 
   template <class RealType>
   void check_pareto(RealType scale, RealType shape, RealType x, RealType p, RealType q, RealType tol)
   {
+    RealType logtol = tol * 10;
+    BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value || 
+                        std::is_same<RealType, boost::math::concepts::real_concept>::value)
+    {
+      logtol *= 100;
+    }
+    
     BOOST_CHECK_CLOSE_FRACTION(
       ::boost::math::cdf(
       pareto_distribution<RealType>(scale, shape),   // distribution.
@@ -57,6 +66,19 @@
       x)),                                           // random variable.
       q,                                             // probability complement.
       tol);                                          // tolerance eps.
+    BOOST_CHECK_CLOSE_FRACTION(
+      ::boost::math::logcdf(
+      pareto_distribution<RealType>(scale, shape),   // distribution.
+      x),                                            // random variable.
+      log(p),                                             // probability.
+      logtol);                                          // tolerance eps.
+    BOOST_CHECK_CLOSE_FRACTION(
+      ::boost::math::logcdf(
+      complement(
+      pareto_distribution<RealType>(scale, shape),   // distribution.
+      x)),                                           // random variable.
+      log(q),                                             // probability complement.
+      logtol);                                          // tolerance eps.
     BOOST_CHECK_CLOSE_FRACTION(
       ::boost::math::quantile(
       pareto_distribution<RealType>(scale, shape),   // distribution.

From 7b678a90010c0fb2ae4d86d47680e324765ad50f Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sat, 11 Feb 2023 18:07:27 -0800
Subject: [PATCH 13/19] Add logcdf specialization to rayleigh distribution

---
 include/boost/math/distributions/rayleigh.hpp | 47 +++++++++++++++++++
 test/test_rayleigh.cpp                        | 22 +++++++++
 2 files changed, 69 insertions(+)

diff --git a/include/boost/math/distributions/rayleigh.hpp b/include/boost/math/distributions/rayleigh.hpp
index cbbf39876d..4e741313c8 100644
--- a/include/boost/math/distributions/rayleigh.hpp
+++ b/include/boost/math/distributions/rayleigh.hpp
@@ -1,4 +1,5 @@
 //  Copyright Paul A. Bristow 2007.
+//  Copyright Matt Borland 2023.
 //  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)
@@ -19,6 +20,7 @@
 #endif
 
 #include <utility>
+#include <limits>
 #include <cmath>
 
 namespace boost{ namespace math{
@@ -168,6 +170,26 @@ inline RealType cdf(const rayleigh_distribution<RealType, Policy>& dist, const R
    return result;
 } // cdf
 
+template <class RealType, class Policy>
+inline RealType logcdf(const rayleigh_distribution<RealType, Policy>& dist, const RealType& x)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   RealType result = 0;
+   RealType sigma = dist.sigma();
+   static const char* function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";
+   if(false == detail::verify_sigma(function, sigma, &result, Policy()))
+   {
+      return -std::numeric_limits<RealType>::infinity();
+   }
+   if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
+   {
+      return -std::numeric_limits<RealType>::infinity();
+   }
+   result = log1p(-exp(-x * x / ( 2 * sigma * sigma)), Policy());   
+   return result;
+} // logcdf
+
 template <class RealType, class Policy>
 inline RealType quantile(const rayleigh_distribution<RealType, Policy>& dist, const RealType& p)
 {
@@ -218,6 +240,31 @@ inline RealType cdf(const complemented2_type<rayleigh_distribution<RealType, Pol
    return result;
 } // cdf complement
 
+template <class RealType, class Policy>
+inline RealType logcdf(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   RealType result = 0;
+   RealType sigma = c.dist.sigma();
+   static const char* function = "boost::math::logcdf(const rayleigh_distribution<%1%>&, %1%)";
+   if(false == detail::verify_sigma(function, sigma, &result, Policy()))
+   {
+      return -std::numeric_limits<RealType>::infinity();
+   }
+   RealType x = c.param;
+   if(false == detail::verify_rayleigh_x(function, x, &result, Policy()))
+   {
+      return -std::numeric_limits<RealType>::infinity();
+   }
+   RealType ea = x * x / (2 * sigma * sigma);
+   // Fix for VC11/12 x64 bug in exp(float):
+   if (ea >= tools::max_value<RealType>())
+      return 0;
+   result = -ea;
+   return result;
+} // logcdf complement
+
 template <class RealType, class Policy>
 inline RealType quantile(const complemented2_type<rayleigh_distribution<RealType, Policy>, RealType>& c)
 {
diff --git a/test/test_rayleigh.cpp b/test/test_rayleigh.cpp
index 6aa2d2fb83..de92dfa848 100644
--- a/test/test_rayleigh.cpp
+++ b/test/test_rayleigh.cpp
@@ -1,4 +1,5 @@
 // Copyright John Maddock 2006.
+// Copyright Matt Borland 2023.
 
 // Use, modification and distribution are subject to the
 // Boost Software License, Version 1.0.
@@ -28,10 +29,19 @@
    using std::setprecision;
 #include <cmath>
    using std::log;
+#include <type_traits>
 
 template <class RealType>
 void test_spot(RealType s, RealType x, RealType p, RealType q, RealType tolerance)
 {
+   RealType logtolerance = tolerance;
+
+   BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value || 
+                       std::is_same<RealType, boost::math::concepts::real_concept>::value)
+   {
+      logtolerance *= 100;
+   }
+   
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
          rayleigh_distribution<RealType>(s),
@@ -44,6 +54,18 @@ void test_spot(RealType s, RealType x, RealType p, RealType q, RealType toleranc
          x)),
          q,
          tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         rayleigh_distribution<RealType>(s),
+         x),
+         log(p),
+         tolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         complement(rayleigh_distribution<RealType>(s),
+         x)),
+         log(q),
+         tolerance); // %
    // Special extra tests for p and q near to unity.
    if(p < 0.999)
    {

From fb958d384c2723006bf6e7e8956c053ac5703486 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sat, 11 Feb 2023 18:08:29 -0800
Subject: [PATCH 14/19] Relax tolerance on edge case

---
 test/test_geometric.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_geometric.cpp b/test/test_geometric.cpp
index ccc3673228..928a2aa0ed 100644
--- a/test/test_geometric.cpp
+++ b/test/test_geometric.cpp
@@ -296,7 +296,7 @@ void test_spots(RealType)
   static_cast<RealType>(-1.0000000000000364e-102L),
   static_cast<RealType>(-std::numeric_limits<RealType>::infinity()),
   tolerance,
-  tolerance);
+  tolerance * 100);
 
   test_spot(  // > formatC(pgeom(1,0.99, TRUE), digits=17) [1] "0.99990000000000001"
     // > formatC(pgeom(1,0.99, FALSE), digits=17) [1] "0.00010000000000000009"

From 1503eb3a8e68422e422178a8ec9e213d9456b54c Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sun, 12 Feb 2023 08:14:58 -0800
Subject: [PATCH 15/19] Change tolerance of long doubles on ARM

---
 test/test_logistic_dist.cpp | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/test/test_logistic_dist.cpp b/test/test_logistic_dist.cpp
index 178f92d34c..3a36cfff6d 100644
--- a/test/test_logistic_dist.cpp
+++ b/test/test_logistic_dist.cpp
@@ -94,6 +94,14 @@ void test_spots(RealType T)
    RealType tolerance = (std::max)(
       static_cast<RealType>(1e-33L),
       boost::math::tools::epsilon<RealType>());
+
+   #if defined(__arm__) || defined(__aarch64__)
+   if (std::is_same<RealType, long double>::value || std::is_same<RealType, boost::math::concepts::real_concept>::value)
+   {
+      tolerance = std::numeric_limits<double>::epsilon();
+   } 
+   #endif
+
    cout<<"Absolute tolerance:"<<tolerance<<endl;
 
    tolerance *= 50 * 100; 

From 56758d7ac239788729eea26cae414ea30a07185b Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sun, 12 Feb 2023 08:16:28 -0800
Subject: [PATCH 16/19] Long double/real concept fixes

---
 test/test_logistic_dist.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_logistic_dist.cpp b/test/test_logistic_dist.cpp
index 3a36cfff6d..90e11d2f5c 100644
--- a/test/test_logistic_dist.cpp
+++ b/test/test_logistic_dist.cpp
@@ -35,7 +35,7 @@
 template <class RealType>
 void test_spot(RealType location, RealType scale, RealType x, RealType p, RealType q, RealType logp, RealType logq, RealType tolerance, RealType logtolerance, RealType logtoleranceq)
 {
-   if (std::is_same<RealType, long double>::value)
+   if (std::is_same<RealType, long double>::value || std::is_same<RealType, boost::math::concepts::real_concept>::value)
    {
       logtoleranceq *= 100;
    }
@@ -157,7 +157,7 @@ void test_spots(RealType T)
       static_cast<RealType>(10), // x
       static_cast<RealType>(0.99999998477002048723965105559179L), // p  
       static_cast<RealType>(1.5229979512760348944408208801237e-8L), //q
-      static_cast<RealType>(-1.5229979628736488101501003766470705e-8), // log(p)
+      static_cast<RealType>(-1.5229979628736488101501003766470705e-8L), // log(p)
       static_cast<RealType>(-18.000000015229979628736488101501L), // Log(q)
       tolerance,
       tolerance * 100000,

From 0b0af30f495aa909a58cdf102bfa0630b58a7192 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Sun, 12 Feb 2023 10:27:52 -0800
Subject: [PATCH 17/19] Change tolerance for non-x86 platforms

---
 test/test_logistic_dist.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_logistic_dist.cpp b/test/test_logistic_dist.cpp
index 90e11d2f5c..3bb092ce7e 100644
--- a/test/test_logistic_dist.cpp
+++ b/test/test_logistic_dist.cpp
@@ -95,7 +95,7 @@ void test_spots(RealType T)
       static_cast<RealType>(1e-33L),
       boost::math::tools::epsilon<RealType>());
 
-   #if defined(__arm__) || defined(__aarch64__)
+   #if !(defined(__amd64__) || defined(__x86_64__) || defined(_M_AMD64))
    if (std::is_same<RealType, long double>::value || std::is_same<RealType, boost::math::concepts::real_concept>::value)
    {
       tolerance = std::numeric_limits<double>::epsilon();

From fcf05006cba8d91922d7e6b8c970aab802a39adf Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Tue, 14 Feb 2023 08:43:47 -0800
Subject: [PATCH 18/19] Add complement logcdf for extreme value distribution

---
 .../math/distributions/extreme_value.hpp      | 24 +++++++++++++
 test/test_extreme_value.cpp                   | 35 +++++++++++++++----
 2 files changed, 53 insertions(+), 6 deletions(-)

diff --git a/include/boost/math/distributions/extreme_value.hpp b/include/boost/math/distributions/extreme_value.hpp
index ac06ec37c2..1bde2743c0 100644
--- a/include/boost/math/distributions/extreme_value.hpp
+++ b/include/boost/math/distributions/extreme_value.hpp
@@ -251,6 +251,30 @@ inline RealType cdf(const complemented2_type<extreme_value_distribution<RealType
    return result;
 }
 
+template <class RealType, class Policy>
+inline RealType logcdf(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
+{
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   static const char* function = "boost::math::logcdf(const extreme_value_distribution<%1%>&, %1%)";
+
+   if((boost::math::isinf)(c.param))
+      return c.param < 0 ? 1.0f : 0.0f;
+   RealType a = c.dist.location();
+   RealType b = c.dist.scale();
+   RealType result = 0;
+   if(0 == detail::verify_scale_b(function, b, &result, Policy()))
+      return result;
+   if(0 == detail::check_finite(function, a, &result, Policy()))
+      return result;
+   if(0 == detail::check_x(function, c.param, &result, Policy()))
+      return result;
+
+   result = log1p(-exp(-exp((a-c.param)/b)), Policy());
+
+   return result;
+}
+
 template <class RealType, class Policy>
 RealType quantile(const complemented2_type<extreme_value_distribution<RealType, Policy>, RealType>& c)
 {
diff --git a/test/test_extreme_value.cpp b/test/test_extreme_value.cpp
index 7e83067429..fd8d928630 100644
--- a/test/test_extreme_value.cpp
+++ b/test/test_extreme_value.cpp
@@ -21,10 +21,16 @@
    using std::cout;
    using std::endl;
    using std::setprecision;
+#include <type_traits>
 
 template <class RealType>
-void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealType logp, RealType tolerance)
+void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealType logp, RealType logq, RealType tolerance, RealType logtolerance)
 {
+   BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value || std::is_same<RealType, boost::math::concepts::real_concept>::value)
+   {
+      logtolerance *= 100;
+   }
+   
    BOOST_CHECK_CLOSE(
       ::boost::math::cdf(
          extreme_value_distribution<RealType>(a, b),      
@@ -42,7 +48,13 @@ void test_spot(RealType a, RealType b, RealType x, RealType p, RealType q, RealT
          extreme_value_distribution<RealType>(a, b),      
          x),
          logp,
-         tolerance); // %
+         logtolerance); // %
+   BOOST_CHECK_CLOSE(
+      ::boost::math::logcdf(
+         complement(extreme_value_distribution<RealType>(a, b),      
+         x)),
+         logq,
+         logtolerance); // %
    if((p < 0.999) && (p > 0))
    {
       BOOST_CHECK_CLOSE(
@@ -85,6 +97,8 @@ void test_spots(RealType)
       static_cast<RealType>(0.27692033409990891617007608217222L), // p
       static_cast<RealType>(0.72307966590009108382992391782778L), //q
       static_cast<RealType>(-1.2840254166877414840734205680624364583362808652814L), // Log(p)
+      static_cast<RealType>(-0.324235874926689525622193916272L), // Log(q)
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
@@ -93,24 +107,33 @@ void test_spots(RealType)
       static_cast<RealType>(1.6087601139887776413169427645933e-7L), // p
       static_cast<RealType>(0.99999983912398860112223586830572L), //q
       static_cast<RealType>(-15.6426318841881716102126980461566588450380350341076L), // Log(p)
+      static_cast<RealType>(-1.60876024339424673820018469895e-7), // Log(q)
+      tolerance,
       tolerance);
    test_spot(
       static_cast<RealType>(0.5), // a
       static_cast<RealType>(0.25), // b
       static_cast<RealType>(0.75), // x
       static_cast<RealType>(0.69220062755534635386542199718279L), // p
-      static_cast<RealType>(0.30779937244465364613457800281721), //q
+      static_cast<RealType>(0.30779937244465364613457800281721L), //q
       static_cast<RealType>(-0.36787944117144232159552377016146086744581113103177L), // Log(p)
+      static_cast<RealType>(-1.17830709642071784241681100298L), // Log(q)
+      tolerance,
       tolerance);
-   test_spot(
+   // Edge case throws overflow exception in complement logcdf for float type
+   BOOST_IF_CONSTEXPR (!std::is_same<RealType, float>::value)
+   {  
+      test_spot(
       static_cast<RealType>(0.5), // a
       static_cast<RealType>(0.25), // b
       static_cast<RealType>(5), // x
       static_cast<RealType>(0.99999998477002037126351248727041L), // p
       static_cast<RealType>(1.5229979628736487512729586276294e-8L), //q
       static_cast<RealType>(-1.52299797447126284361366292335174318621748e-8L), // Log(p)
-      tolerance);
-
+      static_cast<RealType>(-18.0000000076149898626916357587L), // Log(q)
+      tolerance,
+      tolerance * 10000);
+   }
    BOOST_CHECK_CLOSE(
       ::boost::math::pdf(
          extreme_value_distribution<RealType>(0.5, 2),      

From 3886633231130663bb5afff5a88a6004c9151c2f Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Tue, 14 Feb 2023 09:39:28 -0800
Subject: [PATCH 19/19] Add complement logcdf for laplace distribution

---
 include/boost/math/distributions/laplace.hpp | 40 +++++++++++++++
 test/test_laplace.cpp                        | 51 +++++++++++++++++++-
 2 files changed, 90 insertions(+), 1 deletion(-)

diff --git a/include/boost/math/distributions/laplace.hpp b/include/boost/math/distributions/laplace.hpp
index 0a03ff6d9e..81ae8fed9d 100644
--- a/include/boost/math/distributions/laplace.hpp
+++ b/include/boost/math/distributions/laplace.hpp
@@ -347,6 +347,46 @@ inline RealType cdf(const complemented2_type<laplace_distribution<RealType, Poli
    return result;
 } // cdf complement
 
+template <class RealType, class Policy>
+inline RealType logcdf(const complemented2_type<laplace_distribution<RealType, Policy>, RealType>& c)
+{
+   // Calculate complement of logcdf.
+   BOOST_MATH_STD_USING // for ADL of std functions
+
+   RealType scale = c.dist.scale();
+   RealType location = c.dist.location();
+   RealType x = c.param;
+   RealType result = 0;
+
+   // Checking function argument.
+   const char* function = "boost::math::logcdf(const complemented2_type<laplace_distribution<%1%>, %1%>&)";
+
+   // Check scale and location.
+    if (false == c.dist.check_parameters(function, &result)) return result;
+
+   // Special cdf values.
+   if((boost::math::isinf)(x))
+   {
+     if(x < 0) 
+     { 
+       return 1; // cdf complement -infinity is unity.
+     }
+
+     return 0; // cdf complement +infinity is zero.
+   }
+   if(false == detail::check_x(function, x, &result, Policy()))return result;
+
+   // Cdf interval value.
+   if (-x < -location)
+   {
+      result = (-x+location)/scale - boost::math::constants::ln_two<RealType>();
+   }
+   else
+   {
+      result = log1p(-exp( (-location+x)/scale )/2, Policy());
+   }
+   return result;
+} // cdf complement
 
 template <class RealType, class Policy>
 inline RealType quantile(const complemented2_type<laplace_distribution<RealType, Policy>, RealType>& c)
diff --git a/test/test_laplace.cpp b/test/test_laplace.cpp
index fc1d175593..df79c8a584 100644
--- a/test/test_laplace.cpp
+++ b/test/test_laplace.cpp
@@ -69,6 +69,7 @@ Test 8: test_extreme_function_arguments()
 using boost::math::laplace_distribution;
 #include <cmath>
 using std::log;
+#include <type_traits>
 
 /*
 #include <iostream>
@@ -206,6 +207,11 @@ template <class RealType>
 void test_pdf_cdf_ocatave()
 {
    RealType tolerance(1e-10f);
+   RealType logtolerance = std::numeric_limits<RealType>::epsilon() * 500;
+   BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value || std::is_same<RealType, boost::math::concepts::real_concept>::value)
+   {
+      logtolerance *= 100;
+   }
    
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(-2.L)),
@@ -231,6 +237,11 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(-2.6931471805599453094172321214581765680755001319215L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(-2.L))),
+      static_cast<RealType>(-0.0700659201602814101681282515035518814118L),
+   logtolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(-1.L)),
       //static_cast<RealType>(0.18393972058572L),
@@ -255,6 +266,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(-1.6931471805599453094172321214581765680755001342016L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(-1.L))),
+      //log(static_cast<RealType>(0.18393972058572L)),
+      static_cast<RealType>(-0.203267054915195332704071024390384327415L),
+   logtolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(-0.5L)),
    // static_cast<RealType>(0.30326532985632L),
@@ -273,6 +290,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(-1.1931471805599453094172321214581765680755001342281L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(-0.5L))),
+      //static_cast<RealType>(0.30326532985632L),
+      static_cast<RealType>(-0.3613506148087590829918046083597939373459L),
+   logtolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(0.0L)),
       static_cast<RealType>(0.5L),
@@ -294,6 +317,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(-0.69314718055994530941723212145817656807550013436026L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(0.0L))),
+      //log(static_cast<RealType>(0.5L)),
+      static_cast<RealType>(-0.69314718055994530941723212145817656807550013436026L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(0.5L)),
       //static_cast<RealType>(0.30326532985632L),
@@ -317,6 +346,12 @@ void test_pdf_cdf_ocatave()
       static_cast<RealType>(-0.36135061480875908299180460835979393734587166765723L),
    tolerance);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(0.5L))),
+      // static_cast<RealType>(0.69673467014368L),
+      static_cast<RealType>(-1.193147180559945309417232121458176568075L),
+   tolerance);
+
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(1.0L)),
     //  static_cast<RealType>(0.18393972058572L),
@@ -341,7 +376,14 @@ void test_pdf_cdf_ocatave()
      // static_cast<RealType>(0.81606027941428L),
       static_cast<RealType>(-0.2032670549151953327040710243903843274153282048153L),
       //log(static_cast<RealType>(0.81606027941427883920223811491926956627709443427977L)),
-   tolerance*1e5);
+   logtolerance);
+
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(1.00000000000000L))),
+     // static_cast<RealType>(0.81606027941428L),
+      static_cast<RealType>(-1.693147180559945309417232121458176568076L),
+      //log(static_cast<RealType>(0.81606027941427883920223811491926956627709443427977L)),
+   logtolerance);
 
    BOOST_CHECK_CLOSE(
       pdf(laplace_distribution<RealType>(), static_cast<RealType>(2.0L)),
@@ -368,6 +410,13 @@ void test_pdf_cdf_ocatave()
       //log(static_cast<RealType>(0.93233235838169365405300025251375779829618422688019L)),
    tolerance*1e6);
 
+   BOOST_CHECK_CLOSE(
+      logcdf(complement(laplace_distribution<RealType>(), static_cast<RealType>(2.0L))),
+   //   static_cast<RealType>(0.93233235838169L),
+      static_cast<RealType>(-2.693147180559945309417232121458176568075L),
+      //log(static_cast<RealType>(0.93233235838169365405300025251375779829618422688019L)),
+   logtolerance);
+
    check_out_of_range<laplace_distribution<RealType> >(0, 1);
    BOOST_MATH_CHECK_THROW(laplace_distribution<RealType>(0, 0), std::domain_error);
    BOOST_MATH_CHECK_THROW(laplace_distribution<RealType>(0, -1), std::domain_error);