BUG: Use lop1p to improve numerical precision

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

closes #17020

7 files changed, 74 insertions, 36 deletions

diff --git a/numpy/core/include/numpy/random/distributions.h b/numpy/core/include/numpy/random/distributions.h
index 3ffacc8f9..c58024605 100644
--- a/numpy/core/include/numpy/random/distributions.h
+++ b/numpy/core/include/numpy/random/distributions.h
@@ -123,8 +123,9 @@ DECLDIR RAND_INT_TYPE random_negative_binomial(bitgen_t *bitgen_state, double n,
 DECLDIR int64_t random_binomial(bitgen_t *bitgen_state, double p,
                                 int64_t n, binomial_t *binomial);
 
-DECLDIR RAND_INT_TYPE random_logseries(bitgen_t *bitgen_state, double p);
-DECLDIR RAND_INT_TYPE random_geometric(bitgen_t *bitgen_state, double p);
+DECLDIR int64_t random_logseries(bitgen_t *bitgen_state, double p);
+DECLDIR int64_t random_geometric(bitgen_t *bitgen_state, double p);
+DECLDIR RAND_INT_TYPE random_geometric_search(bitgen_t *bitgen_state, double p);
 DECLDIR RAND_INT_TYPE random_zipf(bitgen_t *bitgen_state, double a);
 DECLDIR int64_t random_hypergeometric(bitgen_t *bitgen_state,
                                       int64_t good, int64_t bad, int64_t sample);
diff --git a/numpy/random/_examples/cython/setup.py b/numpy/random/_examples/cython/setup.py
index 83f06fde8..7e0dd3e05 100644
--- a/numpy/random/_examples/cython/setup.py
+++ b/numpy/random/_examples/cython/setup.py
@@ -4,19 +4,20 @@ Build the Cython demonstrations of low-level access to NumPy random
 
 Usage: python setup.py build_ext -i
 """
+from distutils.core import setup
+from os.path import dirname, join, abspath
 
 import numpy as np
-from distutils.core import setup
 from Cython.Build import cythonize
+from numpy.distutils.misc_util import get_info
 from setuptools.extension import Extension
-from os.path import join, dirname
 
 path = dirname(__file__)
 src_dir = join(dirname(path), '..', 'src')
 defs = [('NPY_NO_DEPRECATED_API', 0)]
 inc_path = np.get_include()
-# not so nice. We need the random/lib library from numpy
-lib_path = join(np.get_include(), '..', '..', 'random', 'lib')
+lib_path = [abspath(join(np.get_include(), '..', '..', 'random', 'lib'))]
+lib_path += get_info('npymath')['library_dirs']
 
 extending = Extension("extending",
                       sources=[join('.', 'extending.pyx')],
@@ -29,10 +30,10 @@ extending = Extension("extending",
 distributions = Extension("extending_distributions",
                           sources=[join('.', 'extending_distributions.pyx')],
                           include_dirs=[inc_path],
-                          library_dirs=[lib_path],
-                          libraries=['npyrandom'],
+                          library_dirs=lib_path,
+                          libraries=['npyrandom', 'npymath'],
                           define_macros=defs,
-                         )
+                          )
 
 extensions = [extending, distributions]
 
diff --git a/numpy/random/include/legacy-distributions.h b/numpy/random/include/legacy-distributions.h
index f7ccd2cb5..3d882b73b 100644
--- a/numpy/random/include/legacy-distributions.h
+++ b/numpy/random/include/legacy-distributions.h
@@ -39,7 +39,7 @@ extern int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n,
 extern int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state,
                                             int64_t good, int64_t bad,
                                             int64_t sample);
-extern int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p);
+extern int64_t legacy_logseries(bitgen_t *bitgen_state, double p);
 extern int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam);
 extern int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a);
 extern int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p);
diff --git a/numpy/random/mtrand.pyx b/numpy/random/mtrand.pyx
index 4e12f8e59..d2fa7a074 100644
--- a/numpy/random/mtrand.pyx
+++ b/numpy/random/mtrand.pyx
@@ -89,7 +89,7 @@ cdef extern from "include/legacy-distributions.h":
                                    int64_t n, binomial_t *binomial) nogil
     int64_t legacy_negative_binomial(aug_bitgen_t *aug_state, double n, double p) nogil
     int64_t legacy_random_hypergeometric(bitgen_t *bitgen_state, int64_t good, int64_t bad, int64_t sample) nogil
-    int64_t legacy_random_logseries(bitgen_t *bitgen_state, double p) nogil
+    int64_t legacy_logseries(bitgen_t *bitgen_state, double p) nogil
     int64_t legacy_random_poisson(bitgen_t *bitgen_state, double lam) nogil
     int64_t legacy_random_zipf(bitgen_t *bitgen_state, double a) nogil
     int64_t legacy_random_geometric(bitgen_t *bitgen_state, double p) nogil
@@ -3954,7 +3954,7 @@ cdef class RandomState:
         >>> plt.show()
 
         """
-        out = disc(&legacy_random_logseries, &self._bitgen, size, self.lock, 1, 0,
+        out = disc(&legacy_logseries, &self._bitgen, size, self.lock, 1, 0,
                    p, 'p', CONS_BOUNDED_0_1,
                    0.0, '', CONS_NONE,
                    0.0, '', CONS_NONE)
diff --git a/numpy/random/setup.py b/numpy/random/setup.py
index bfd08e469..dce9a101e 100644
--- a/numpy/random/setup.py
+++ b/numpy/random/setup.py
@@ -23,7 +23,7 @@ def configuration(parent_package='', top_path=None):
 
     # enable unix large file support on 32 bit systems
     # (64 bit off_t, lseek -> lseek64 etc.)
-    if sys.platform[:3] == "aix":
+    if sys.platform[:3] == 'aix':
         defs = [('_LARGE_FILES', None)]
     else:
         defs = [('_FILE_OFFSET_BITS', '64'),
@@ -116,7 +116,7 @@ def configuration(parent_package='', top_path=None):
         # gen.pyx, src/distributions/distributions.c
         config.add_extension(gen,
                              sources=[f'{gen}.c'],
-                             libraries=EXTRA_LIBRARIES,
+                             libraries=EXTRA_LIBRARIES + ['npymath'],
                              extra_compile_args=EXTRA_COMPILE_ARGS,
                              include_dirs=['.', 'src'],
                              extra_link_args=EXTRA_LINK_ARGS,
@@ -124,13 +124,14 @@ def configuration(parent_package='', top_path=None):
                              define_macros=defs,
                              )
     config.add_data_files('_bounded_integers.pxd')
+    mtrand_libs = ['m', 'npymath'] if os.name != 'nt' else ['npymath']
     config.add_extension('mtrand',
                          sources=['mtrand.c',
                                   'src/legacy/legacy-distributions.c',
                                   'src/distributions/distributions.c',
                                  ],
                          include_dirs=['.', 'src', 'src/legacy'],
-                         libraries=['m'] if os.name != 'nt' else [],
+                         libraries=mtrand_libs,
                          extra_compile_args=EXTRA_COMPILE_ARGS,
                          extra_link_args=EXTRA_LINK_ARGS,
                          depends=depends + ['mtrand.pyx'],
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