Merge pull request #1297 from boostorg/1296

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

doc/math.qbk

@@ -339,6 +339,10 @@ and use the function's name as the link text.]
 [def __hermite [link math_toolkit.sf_poly.hermite hermite]]
 [def __cardinal_b_splines [link math_toolkit.sf_poly.cardinal_b_splines cardinal_b_splines]]
 
+[/logistic functions]
+[def __logit [link math_toolkit.logistic.logit logit]]
+[def __logistic_sigmoid [link math_toolkit.logistic.logistic_sigmoid logistic_sigmoid]]
+
 [/Misc]
 [def __expint [link math_toolkit.expint.expint_i expint]]
 [def __spherical_harmonic [link math_toolkit.sf_poly.sph_harm spherical_harmonic]]
@@ -666,6 +670,11 @@ and as a CD ISBN 0-9504833-2-X  978-0-9504833-2-0, Classification 519.2-dc22.
 [include sf/jacobi.qbk]
 [endsect] [/section:sf_poly Polynomials]
 
+[section:logistic Logistic Functions]
+[include sf/logit.qbk]
+[include sf/logistic_sigmoid.qbk]
+[endsect] [/section:logistic Logistic Functions]
+
 [section:bessel Bessel Functions]
 [include sf/bessel_introduction.qbk]
 [include sf/bessel_jy.qbk]

doc/sf/logistic_sigmoid.qbk

@@ -0,0 +1,30 @@
+[/
+  Copyright Matt Borland 2025
+  Distributed under 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).
+]
+
+[section:logistic_sigmoid logistic_sigmoid]
+
+[h4 Synopsis]
+
+    #include <boost/math/special_functions/logistic_sigmoid.hpp>
+
+    namespace boost { namespace math {
+
+    template <typename RealType, typename Policy>
+    RealType logistic_sigmoid(RealType x, const Policy&);
+
+    template <typename RealType>
+    RealType logistic_sigmoid(RealType x)
+
+    }} // namespaces
+
+[h4 Description]
+
+Returns the [@https://en.wikipedia.org/wiki/Logistic_function logistic sigmoid function]
+This function is broadly useful, and is used to compute the CDF of the logistic distribution.
+It is also sometimes referred to as expit as it is the inverse of the logit function.
+
+[endsect]

doc/sf/logit.qbk

@@ -0,0 +1,30 @@
+[/
+  Copyright Matt Borland 2025
+  Distributed under 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).
+]
+
+[section:logit logit]
+
+[h4 Synopsis]
+
+    #include <boost/math/special_functions/logit.hpp>
+
+    namespace boost { namespace math {
+
+    template <typename RealType, typename Policy>
+    RealType logit(RealType x, const Policy&);
+
+    template <typename RealType>
+    RealType logit(RealType x)
+
+    }} // namespaces
+
+[h4 Description]
+
+Returns the [@https://en.wikipedia.org/wiki/Logit logit function]
+This function is broadly useful, and is used to compute the Quantile of the logistic distribution.
+The inverse of this function is the logistic_sigmoid.
+
+[endsect]

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 { 
 
@@ -144,13 +146,10 @@ 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)); 
+
+       using promoted_real_type = typename policies::evaluation<RealType, Policy>::type;
+       promoted_real_type power = (static_cast<promoted_real_type>(x) - static_cast<promoted_real_type>(location)) / static_cast<promoted_real_type>(scale);
+       return logistic_sigmoid(power, policies::make_forwarding_policy_t<Policy>());
     }
 
     template <class RealType, class Policy>
@@ -199,7 +198,6 @@ namespace boost { namespace math {
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType quantile(const logistic_distribution<RealType, Policy>& dist, const RealType& p)
     {
-       BOOST_MATH_STD_USING
        RealType location = dist.location();
        RealType scale = dist.scale();
 
@@ -221,21 +219,13 @@ 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>
     BOOST_MATH_GPU_ENABLED inline RealType cdf(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;
@@ -259,12 +249,10 @@ 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)); 
+
+       using promoted_real_type = typename policies::evaluation<RealType, Policy>::type;
+       promoted_real_type power = (static_cast<promoted_real_type>(location) - static_cast<promoted_real_type>(x)) / static_cast<promoted_real_type>(scale);
+       return logistic_sigmoid(power, policies::make_forwarding_policy_t<Policy>());
     } 
 
     template <class RealType, class Policy>
@@ -306,7 +294,6 @@ namespace boost { namespace math {
     template <class RealType, class Policy>
     BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<logistic_distribution<RealType, Policy>, RealType>& c)
     {
-       BOOST_MATH_STD_USING
        RealType scale = c.dist.scale();
        RealType location = c.dist.location();
        constexpr auto function = "boost::math::quantile(const complement(logistic_distribution<%1%>&), %1%)";
@@ -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/policies/policy.hpp

@@ -1003,6 +1003,25 @@ struct is_noexcept_error_policy
       && (t8::value != throw_on_error) && (t8::value != user_error));
 };
 
+// Generate a forwarding policy to stop further promotion from occurring
+// For example if a special function for float promotes to double, we don't want the next
+// function in the call chain to then promote to long double
+template <typename Policy>
+struct make_forwarding_policy
+{
+   using type =
+          typename policies::normalise<
+             Policy,
+             policies::promote_float<false>,
+             policies::promote_double<false>,
+             policies::discrete_quantile<>,
+             policies::assert_undefined<>
+          >::type;
+};
+
+template <typename Policy>
+using make_forwarding_policy_t = typename make_forwarding_policy<Policy>::type;
+
 }}} // namespaces
 
 #endif // BOOST_MATH_POLICY_HPP

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();
+}

test/test_logit.cpp

@@ -0,0 +1,107 @@
+// 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/logit.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include "math_unit_test.hpp"
+#include <array>
+#include <limits>
+#include <cfenv>
+
+#pragma STDC FENV_ACCESS ON
+
+template <typename RealType>
+void test()
+{
+    const std::array<RealType, 5> x_values = {
+        0.01,
+        0.24,
+        0.5,
+        0.75,
+        0.995,
+    };
+    const std::array<RealType, 5> y_values = {
+        RealType{-4.595119850134589926852434051810180709116687969582916078687956376405},
+        RealType{-1.15267950993838545919655007350715126451438856911612411268258589327840479},
+        RealType{0},
+        RealType{1.09861228866810969139524523692252570464749055782274945173469433363749429},
+        RealType{5.2933048247244923954101212918685372018911052805694724989064609879440992}
+    };
+
+    for (std::size_t i = 0; i < x_values.size(); ++i)
+    {
+        const RealType test_value {boost::math::logit(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], 5);
+        }
+        else
+        {
+            CHECK_MOLLIFIED_CLOSE(test_value, y_values[i], 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);
+    }
+
+    #if defined(_CPPUNWIND) || defined(__EXCEPTIONS)
+
+    BOOST_MATH_IF_CONSTEXPR (std::is_arithmetic<RealType>::value)
+    {
+        bool thrown {false};
+        try
+        {
+            boost::math::logit(std::numeric_limits<RealType>::denorm_min());
+        }
+        catch (...)
+        {
+            thrown = true;
+        }
+
+        CHECK_EQUAL(thrown, true);
+    }
+
+    #endif // Exceptional environments
+}
+
+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>();
+
+    return boost::math::test::report_errors();
+}