diff --git a/include/xsf/legendre.h b/include/xsf/legendre.h
index 23fb9e8f3..54105ee4f 100644
--- a/include/xsf/legendre.h
+++ b/include/xsf/legendre.h
@@ -301,7 +301,11 @@ struct assoc_legendre_p_recurrence_n<T, assoc_legendre_norm_policy> {
 template <typename NormPolicy, typename T>
 void assoc_legendre_p_pm1(NormPolicy norm, int n, int m, T z, int branch_cut, T &res) {
     if (m == 0) {
-        res = T(1);
+        if (real(z) >= 0) {
+            res = T(1);
+        } else {
+            res = T(std::pow(-1, n));
+        }
     } else {
         res = T(0);
     }
@@ -310,7 +314,11 @@ void assoc_legendre_p_pm1(NormPolicy norm, int n, int m, T z, int branch_cut, T
 template <typename NormPolicy, typename T, size_t Order>
 void assoc_legendre_p_pm1(NormPolicy norm, int n, int m, dual<T, Order> z, int branch_cut, dual<T, Order> &res) {
     if (m == 0) {
-        res[0] = T(1);
+        if (real(z[0]) >= 0) {
+            res[0] = T(1);
+        } else {
+            res[0] = T(std::pow(-1, n));
+        }
     } else {
         res[0] = T(0);
     }
diff --git a/tests/xsf_tests/test_assoc_legendre_p.cpp b/tests/xsf_tests/test_assoc_legendre_p.cpp
new file mode 100644
index 000000000..3378ab663
--- /dev/null
+++ b/tests/xsf_tests/test_assoc_legendre_p.cpp
@@ -0,0 +1,23 @@
+#include "../testing_utils.h"
+#include <tuple>
+#include <xsf/legendre.h>
+
+TEST_CASE("assoc_legendre_p scipy/gh-23101", "[assoc_legendre_p][xsf_tests]") {
+    using test_case = std::tuple<int, int, double, double, double>;
+    // Reference values were computed with the Python library mpmath.
+    auto [n, m, z, ref, rtol] = GENERATE(
+        test_case{1, 0, -1.0, -1.0, 1e-15}, test_case{1, 0, std::nextafter(-1.0, -2.0), -1.0000000000000002, 1e-15},
+        test_case{1, 0, std::nextafter(-1.0, 0.0), -0.9999999999999999, 1e-15}, test_case{2, 0, -1.0, 1.0, 1e-15},
+        test_case{2, 0, std::nextafter(-1.0, -2.0), 1.0000000000000007, 1e-15},
+        test_case{2, 0, std::nextafter(-1.0, 0.0), 0.9999999999999997, 1e-15}, test_case{3, 0, -1.0, -1.0, 1e-15},
+        test_case{3, 0, std::nextafter(-1.0, -2.0), -1.0000000000000013, 1e-15},
+        test_case{3, 0, std::nextafter(-1.0, 0.0), -0.9999999999999993, 1e-15}, test_case{4, 0, -1.0, 1.0, 1e-15},
+        test_case{4, 0, std::nextafter(-1.0, -2.0), 1.0000000000000022, 4e-15},
+        test_case{4, 0, std::nextafter(-1.0, 0.0), 0.9999999999999989, 1e-15}
+    );
+    const double w = xsf::assoc_legendre_p(xsf::assoc_legendre_unnorm, n, m, z, 2);
+    const auto rel_error = xsf::extended_relative_error(w, ref);
+
+    CAPTURE(n, m, z, w, ref, rtol, rel_error);
+    REQUIRE(rel_error <= rtol);
+}