MAINT: special: Add log1mexp for log(1 - exp(x)) (#22439)

* MAINT: special: add private function _log1mexp
* TST: special: add tests for log1mexp
* MAINT: stats: replace Python _log1mexp with ufunc

---------

Co-authored-by: Matt Haberland <[email protected]>

Author: steppi

include/xsf/log_exp.h

@@ -3,6 +3,7 @@
 #include <cmath>
 
 #include "config.h"
+#include "error.h"
 
 namespace xsf {
 
@@ -65,4 +66,22 @@ T log_expit(T x) {
     return -std::log1p(std::exp(-x));
 };
 
+
+/* Compute log(1 - exp(x)). */
+template <typename T>
+T log1mexp(T x) {
+    if (x > 0) {
+	set_error("_log1mexp", SF_ERROR_DOMAIN, NULL);
+	return std::numeric_limits<T>::quiet_NaN();
+    }
+    if (x == 0) {
+	set_error("_log1mexp", SF_ERROR_SINGULAR, NULL);
+	return -std::numeric_limits<T>::infinity();
+    }
+    if (x < -1) {
+	return std::log1p(-std::exp(x));
+    }
+    return std::log(-std::expm1(x));
+}
+
 } // namespace xsf