BUG: Use lop1p to improve numerical precision

Use log1p(-x) instead of log(1 - x)
Seperate legacy version from current

closes #17020

2 files changed, 56 insertions, 21 deletions

diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
index 3df819d9d..6b4deb925 100644
--- a/numpy/random/src/distributions/distributions.c
+++ b/numpy/random/src/distributions/distributions.c
@@ -906,15 +906,9 @@ double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
   }
 }
 
-/*
- * RAND_INT_TYPE is used to share integer generators with RandomState which
- * used long in place of int64_t. If changing a distribution that uses
- * RAND_INT_TYPE, then the original unmodified copy must be retained for
- * use in RandomState by copying to the legacy distributions source file.
- */
-RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p) {
+int64_t random_logseries(bitgen_t *bitgen_state, double p) {
   double q, r, U, V;
-  RAND_INT_TYPE result;
+  int64_t result;
 
   r = npy_log1p(-p);
 
@@ -926,7 +920,7 @@ RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p) {
     U = next_double(bitgen_state);
     q = 1.0 - exp(r * U);
     if (V <= q * q) {
-      result = (RAND_INT_TYPE)floor(1 + log(V) / log(q));
+      result = (int64_t)floor(1 + log(V) / log(q));
       if ((result < 1) || (V == 0.0)) {
         continue;
       } else {
@@ -940,6 +934,14 @@ RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p) {
   }
 }
 
+/*
+ * RAND_INT_TYPE is used to share integer generators with RandomState which
+ * used long in place of int64_t. If changing a distribution that uses
+ * RAND_INT_TYPE, then the original unmodified copy must be retained for
+ * use in RandomState by copying to the legacy distributions source file.
+ */
+
+/* Still used but both generator and mtrand via legacy_random_geometric */
 RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p) {
   double U;
   RAND_INT_TYPE X;
@@ -957,11 +959,11 @@ RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p) {
   return X;
 }
 
-RAND_INT_TYPE random_geometric_inversion(bitgen_t *bitgen_state, double p) {
-  return (RAND_INT_TYPE)ceil(npy_log1p(-next_double(bitgen_state)) / npy_log1p(-p));
+int64_t random_geometric_inversion(bitgen_t *bitgen_state, double p) {
+  return (int64_t)ceil(npy_log1p(-next_double(bitgen_state)) / npy_log1p(-p));
 }
 
-RAND_INT_TYPE random_geometric(bitgen_t *bitgen_state, double p) {
+int64_t random_geometric(bitgen_t *bitgen_state, double p) {
   if (p >= 0.333333333333333333333333) {
     return random_geometric_search(bitgen_state, p);
   } else {
diff --git a/numpy/random/src/legacy/legacy-distributions.c b/numpy/random/src/legacy/legacy-distributions.c
index 5b17401dd..bfea15e40 100644
--- a/numpy/random/src/legacy/legacy-distributions.c
+++ b/numpy/random/src/legacy/legacy-distributions.c
@@ -379,23 +379,28 @@ int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, int64_t good,
 }
 
 
-int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) {
-  return (int64_t)random_logseries(bitgen_state, p);
-}
-
- int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam) {
+int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam) {
   return (int64_t)random_poisson(bitgen_state, lam);
 }
 
- int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a) {
+int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a) {
   return (int64_t)random_zipf(bitgen_state, a);
 }
 
- int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p) {
-  return (int64_t)random_geometric(bitgen_state, p);
+
+static long legacy_geometric_inversion(bitgen_t *bitgen_state, double p) {
+  return (long)ceil(npy_log1p(-next_double(bitgen_state)) / log(1 - p));
+}
+
+int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p) {
+  if (p >= 0.333333333333333333333333) {
+    return (int64_t)random_geometric_search(bitgen_state, p);
+  } else {
+    return (int64_t)legacy_geometric_inversion(bitgen_state, p);
+  }
 }
 
- void legacy_random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n,
+void legacy_random_multinomial(bitgen_t *bitgen_state, RAND_INT_TYPE n,
                                RAND_INT_TYPE *mnix, double *pix, npy_intp d,
                                binomial_t *binomial) {
   return random_multinomial(bitgen_state, n, mnix, pix, d, binomial);
@@ -457,4 +462,32 @@ double legacy_vonmises(bitgen_t *bitgen_state, double mu, double kappa) {
 
     return mod;
   }
+}
+
+int64_t legacy_logseries(bitgen_t *bitgen_state, double p) {
+  double q, r, U, V;
+  long result;
+
+  r = log(1.0 - p);
+
+  while (1) {
+    V = next_double(bitgen_state);
+    if (V >= p) {
+      return 1;
+    }
+    U = next_double(bitgen_state);
+    q = 1.0 - exp(r * U);
+    if (V <= q * q) {
+      result = (long)floor(1 + log(V) / log(q));
+      if ((result < 1) || (V == 0.0)) {
+        continue;
+      } else {
+        return (int64_t)result;
+      }
+    }
+    if (V >= q) {
+      return 1;
+    }
+    return 2;
+  }
 }
\ No newline at end of file