Add logit, logistic_sigmoid, and use in logistic dist to support promotion policies

include/boost/math/distributions/logistic.hpp

@@ -19,6 +19,8 @@
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/policies/policy.hpp>
 #include <boost/math/policies/error_handling.hpp>
+#include <boost/math/special_functions/logit.hpp>
+#include <boost/math/special_functions/logistic_sigmoid.hpp>
 
 namespace boost { namespace math { 
 
@@ -145,12 +147,8 @@ namespace boost { namespace math {
           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 1 / (1 + exp(power)); 
+       RealType power = -(location - x) / scale;
+       return logistic_sigmoid(power, Policy());
     }
 
     template <class RealType, class Policy>
@@ -221,15 +219,8 @@ namespace boost { namespace math {
        {
           return policies::raise_overflow_error<RealType>(function,"probability argument is 1, must be >0 and <1",Policy());
        }
-       //Expressions to try
-       //return location+scale*log(p/(1-p));
-       //return location+scale*log1p((2*p-1)/(1-p));
-
-       //return location - scale*log( (1-p)/p);
-       //return location - scale*log1p((1-2*p)/p);
 
-       //return -scale*log(1/p-1) + location;
-       return location - scale * log((1 - p) / p);
+       return location - scale * -logit(p, Policy());
      } // RealType quantile(const logistic_distribution<RealType, Policy>& dist, const RealType& p)
 
     template <class RealType, class Policy>
@@ -259,12 +250,8 @@ namespace boost { namespace math {
        {
           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 1 / (1 + exp(power)); 
+       RealType power = -(x - location) / scale;
+       return logistic_sigmoid(power, Policy());
     } 
 
     template <class RealType, class Policy>
@@ -328,15 +315,8 @@ namespace boost { namespace math {
        {
           return policies::raise_overflow_error<RealType>(function,"probability argument is 0, but must be >0 and <1",Policy());
        }
-       //Expressions to try 
-       //return location+scale*log((1-q)/q);
-       return location + scale * log((1 - q) / q);
-
-       //return location-scale*log(q/(1-q));
-       //return location-scale*log1p((2*q-1)/(1-q));
 
-       //return location+scale*log(1/q-1);
-       //return location+scale*log1p(1/q-2);
+       return location + scale * -logit(q, Policy());
     } 
 
     template <class RealType, class Policy>

include/boost/math/special_functions/logistic_sigmoid.hpp

@@ -0,0 +1,46 @@
+// Copyright Matt Borland 2025.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SF_EXPIT_HPP
+#define BOOST_MATH_SF_EXPIT_HPP
+
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/tools/precision.hpp>
+#include <cmath>
+
+namespace boost {
+namespace math {
+
+template <typename RealType, typename Policy>
+RealType logistic_sigmoid(RealType x, const Policy&)
+{
+    BOOST_MATH_STD_USING
+
+    using promoted_real_type = typename policies::evaluation<RealType, Policy>::type;
+
+    if(-x >= tools::log_max_value<RealType>())
+    {
+        return static_cast<RealType>(0);
+    }
+    if(-x <= -tools::log_max_value<RealType>())
+    {
+        return static_cast<RealType>(1);
+    }
+
+    const auto res {static_cast<RealType>(1 / (1 + exp(static_cast<promoted_real_type>(-x))))};
+    return res;
+}
+
+template <typename RealType>
+RealType logistic_sigmoid(RealType x)
+{
+    return logistic_sigmoid(x, policies::policy<>());
+}
+
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_SF_EXPIT_HPP

include/boost/math/special_functions/logit.hpp

@@ -0,0 +1,56 @@
+// Copyright Matt Borland 2025.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SF_LOGIT_HPP
+#define BOOST_MATH_SF_LOGIT_HPP
+
+#include <boost/math/tools/config.hpp>
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/policies/error_handling.hpp>
+#include <cmath>
+#include <cfenv>
+
+namespace boost {
+namespace math {
+
+template <typename RealType, typename Policy>
+RealType logit(RealType p, const Policy&)
+{
+    BOOST_MATH_STD_USING
+    using std::atanh;
+
+    using promoted_real_type = typename policies::evaluation<RealType, Policy>::type;
+
+    if (p < tools::min_value<RealType>())
+    {
+        return -policies::raise_overflow_error<RealType>("logit", "sub-normals will overflow ln(x/(1-x))", Policy());
+    }
+
+    static const RealType crossover {RealType{1}/4};
+    const auto promoted_p {static_cast<promoted_real_type>(p)};
+    RealType result {};
+    if (p > crossover)
+    {
+        result = static_cast<RealType>(2 * atanh(2 * promoted_p - 1));
+    }
+    else
+    {
+        result = static_cast<RealType>(log(promoted_p / (1 - promoted_p)));
+    }
+
+    return result;
+}
+
+template <typename RealType>
+RealType logit(RealType p)
+{
+    return logit(p, policies::policy<>());
+}
+
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_SF_LOGIT_HPP

test/Jamfile.v2

@@ -590,6 +590,8 @@ test-suite special_fun :
    [ run test_sinc.cpp /boost/test//boost_unit_test_framework pch_light ]
    [ run test_fibonacci.cpp /boost/test//boost_unit_test_framework ]
    [ run test_prime.cpp /boost/test//boost_unit_test_framework ]
+   [ run test_logistic_sigmoid.cpp ]
+   [ run test_logit.cpp ]
 ;
 
 test-suite distribution_tests :

test/git_issue_1294.cpp

@@ -0,0 +1,33 @@
+// Copyright 2025 Matt Borland
+// 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)
+//
+// See: https://github.com/boostorg/math/issues/1294
+
+#include <boost/math/distributions/logistic.hpp>
+#include "math_unit_test.hpp"
+
+int main()
+{
+    using namespace boost::math::policies;
+    using boost::math::logistic_distribution;
+
+    typedef policy<
+        promote_float<true>,
+        promote_double<true>
+    > with_promotion;
+
+    constexpr double p = 2049.0/4096;
+    constexpr double ref = 9.76562577610225755e-04;
+
+    logistic_distribution<double, with_promotion> dist_promote;
+    const double x = quantile(dist_promote, p);
+
+    // Previously we had: 9.76562577610170027e-04
+    // Which is an ULP distance of 256
+    CHECK_ULP_CLOSE(x, ref, 1);
+
+    return boost::math::test::report_errors();
+}

test/math_unit_test.hpp

@@ -58,7 +58,7 @@ bool check_mollified_close(Real expected, Real computed, Real tol, std::string c
     }
     using std::max;
     using std::abs;
-    Real denom = (max)(abs(expected), Real(1));
+    Real denom = (max)(Real(abs(expected)), Real(1));
     Real mollified_relative_error = abs(expected - computed)/denom;
     if (mollified_relative_error > tol)
     {

test/test_logistic_sigmoid.cpp

@@ -0,0 +1,91 @@
+// Copyright Matt Borland 2025.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/math/special_functions/logistic_sigmoid.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include "math_unit_test.hpp"
+#include <array>
+#include <cfloat>
+#include <cfenv>
+
+#pragma STDC FENV_ACCESS ON
+
+template <typename RealType>
+void test()
+{
+    const std::array<RealType, 5> x_values = {
+        0,
+        1,
+        1000,
+        0.5,
+        0.75
+    };
+    const std::array<RealType, 5> y_values = {
+        static_cast<RealType>(1) / 2,
+        static_cast<RealType>(0.73105857863000487925115924182183627436514464016505651927636590791904045307),
+        static_cast<RealType>(1),
+        static_cast<RealType>(0.62245933120185456463890056574550847875327936530891016305943716265854500),
+        static_cast<RealType>(0.6791786991753929731596801157765790212342212482195760219829517436805)
+    };
+
+    for (std::size_t i = 0; i < x_values.size(); ++i)
+    {
+        const RealType test_value {boost::math::logistic_sigmoid(x_values[i])};
+        BOOST_MATH_IF_CONSTEXPR (std::is_same<RealType, float>::value || std::is_same<RealType, double>::value)
+        {
+            CHECK_ULP_CLOSE(test_value, y_values[i], 1);
+        }
+        else
+        {
+            RealType comparison_value = y_values[i];
+            CHECK_MOLLIFIED_CLOSE(test_value, comparison_value, static_cast<RealType>(1e-15));
+        }
+
+        bool fe {false};
+        if (std::fetestexcept(FE_OVERFLOW))
+        {
+            fe = true;                                  // LCOV_EXCL_LINE
+            std::cerr << "FE_OVERFLOW" << std::endl;    // LCOV_EXCL_LINE
+        }
+        if (std::fetestexcept(FE_UNDERFLOW))
+        {
+            fe = true;                                  // LCOV_EXCL_LINE
+            std::cerr << "FE_UNDERFLOW" << std::endl;   // LCOV_EXCL_LINE
+        }
+        if (std::fetestexcept(FE_DIVBYZERO))
+        {
+            fe = true;                                  // LCOV_EXCL_LINE
+            std::cerr << "FE_DIVBYZERO" << std::endl;   // LCOV_EXCL_LINE
+        }
+        if (std::fetestexcept(FE_INVALID))
+        {
+            fe = true;                                  // LCOV_EXCL_LINE
+            std::cerr << "FE_INVALID" << std::endl;     // LCOV_EXCL_LINE
+        }
+
+        CHECK_EQUAL(fe, false);
+    }
+}
+
+int main()
+{
+    std::feclearexcept(FE_ALL_EXCEPT);
+    test<float>();
+
+    std::feclearexcept(FE_ALL_EXCEPT);
+    test<double>();
+
+    std::feclearexcept(FE_ALL_EXCEPT);
+    test<long double>();
+
+    std::feclearexcept(FE_ALL_EXCEPT);
+    test<boost::multiprecision::cpp_bin_float_quad>();
+
+    test<boost::multiprecision::cpp_dec_float_50>();
+
+    return boost::math::test::report_errors();
+}