API: restructure and document numpy.random C-API (#14604)

* API: restructure and document numpy.random C-API

* DOC: fix bad reference

* API: ship, document, and start to test numpy.random C-API examples

* API, DOC, TST: fix tests, refactor documentation to include snippets

* BUILD: move public headers to numpy/core/include/numpy/random

* TST: ignore DeprecationWarnings in setuptools and numba

* DOC: document the C-API as used from Cython

30 files changed, 377 insertions, 187 deletions

diff --git a/doc/source/reference/random/bit_generators/index.rst b/doc/source/reference/random/bit_generators/index.rst
index 94d3d8a3c..315657172 100644
--- a/doc/source/reference/random/bit_generators/index.rst
+++ b/doc/source/reference/random/bit_generators/index.rst
@@ -19,7 +19,7 @@ The included BitGenerators are:
   and can be advanced by an arbitrary amount. See the documentation for
   :meth:`~.PCG64.advance`. PCG-64 has a period of :math:`2^{128}`. See the `PCG
   author's page`_ for more details about this class of PRNG.
-* MT19937 - The standard Python BitGenerator. Adds a `~mt19937.MT19937.jumped`
+* MT19937 - The standard Python BitGenerator. Adds a `MT19937.jumped`
   function that returns a new generator with state as-if :math:`2^{128}` draws have
   been made.
 * Philox - A counter-based generator capable of being advanced an
diff --git a/doc/source/reference/random/c-api.rst b/doc/source/reference/random/c-api.rst
new file mode 100644
index 000000000..3c901f3b4
--- /dev/null
+++ b/doc/source/reference/random/c-api.rst
@@ -0,0 +1,177 @@
+Cython API for random
+---------------------
+
+.. currentmodule:: numpy.random
+
+Typed versions of many of the `Generator` and `BitGenerator` can be accessed
+directly from Cython: the complete list is given below.
+
+The ``_bit_generator`` module is usable via::
+
+    cimport numpy.random._bit_generator
+
+It provides function pointers for quickly accessing the next bytes in the
+`BitGenerator`::
+
+    struct bitgen:
+        void *state
+        uint64_t (*next_uint64)(void *st) nogil
+        uint32_t (*next_uint32)(void *st) nogil
+        double (*next_double)(void *st) nogil
+        uint64_t (*next_raw)(void *st) nogil
+
+    ctypedef bitgen bitgen_t
+
+See `extending` for examples of using these functions.
+
+The ``_generator`` module is usable via::
+
+    cimport numpy.random._generator
+
+It provides low-level functions for various distributions. All the functions require a ``bitgen_t`` BitGenerator structure. The functions are named with the followig cconventions:
+
+- "standard" refers to the reference values for any parameters. For instance
+  "standard_uniform" means a uniform distribution on the interval ``0.0`` to
+  ``1.0``
+
+- "fill" functions will fill the provided ``out`` with ``cnt`` values.
+
+- The functions without "standard" in their name require additional parameters
+  to describe the distributions.
+
+.. c:function:: double random_standard_uniform(bitgen_t *bitgen_state)
+
+.. c:function:: void random_standard_uniform_fill(bitgen_t* bitgen_state, np.npy_intp cnt, double *out)
+
+.. c:function:: double random_standard_exponential(bitgen_t *bitgen_state)
+
+.. c:function:: void random_standard_exponential_fill(bitgen_t *bitgen_state, np.npy_intp cnt, double *out)
+
+.. c:function:: double random_standard_exponential_zig(bitgen_t *bitgen_state)
+
+.. c:function:: void random_standard_exponential_zig_fill(bitgen_t *bitgen_state, np.npy_intp cnt, double *out)
+
+.. c:function:: double random_standard_normal(bitgen_t* bitgen_state)
+
+.. c:function:: void random_standard_normal_fill(bitgen_t *bitgen_state, np.npy_intp count, double *out)
+
+.. c:function:: void random_standard_normal_fill_f(bitgen_t *bitgen_state, np.npy_intp count, float *out)
+
+.. c:function:: double random_standard_gamma(bitgen_t *bitgen_state, double shape)
+
+.. c:function:: float random_standard_uniform_f(bitgen_t *bitgen_state)
+
+.. c:function:: void random_standard_uniform_fill_f(bitgen_t* bitgen_state, np.npy_intp cnt, float *out)
+
+.. c:function:: float random_standard_exponential_f(bitgen_t *bitgen_state)
+
+.. c:function:: float random_standard_exponential_zig_f(bitgen_t *bitgen_state)
+
+.. c:function:: void random_standard_exponential_fill_f(bitgen_t *bitgen_state, np.npy_intp cnt, float *out)
+
+.. c:function:: void random_standard_exponential_zig_fill_f(bitgen_t *bitgen_state, np.npy_intp cnt, float *out)
+
+.. c:function:: float random_standard_normal_f(bitgen_t* bitgen_state)
+
+.. c:function:: float random_standard_gamma_f(bitgen_t *bitgen_state, float shape)
+
+.. c:function:: double random_normal(bitgen_t *bitgen_state, double loc, double scale)
+
+.. c:function:: double random_gamma(bitgen_t *bitgen_state, double shape, double scale)
+
+.. c:function:: float random_gamma_f(bitgen_t *bitgen_state, float shape, float scale)
+
+.. c:function:: double random_exponential(bitgen_t *bitgen_state, double scale)
+
+.. c:function:: double random_uniform(bitgen_t *bitgen_state, double lower, double range)
+.. c:function:: double random_beta(bitgen_t *bitgen_state, double a, double b)
+
+.. c:function:: double random_chisquare(bitgen_t *bitgen_state, double df)
+
+.. c:function:: double random_f(bitgen_t *bitgen_state, double dfnum, double dfden)
+
+.. c:function:: double random_standard_cauchy(bitgen_t *bitgen_state)
+
+.. c:function:: double random_pareto(bitgen_t *bitgen_state, double a)
+
+.. c:function:: double random_weibull(bitgen_t *bitgen_state, double a)
+
+.. c:function:: double random_power(bitgen_t *bitgen_state, double a)
+
+.. c:function:: double random_laplace(bitgen_t *bitgen_state, double loc, double scale)
+
+.. c:function:: double random_gumbel(bitgen_t *bitgen_state, double loc, double scale)
+
+.. c:function:: double random_logistic(bitgen_t *bitgen_state, double loc, double scale)
+
+.. c:function:: double random_lognormal(bitgen_t *bitgen_state, double mean, double sigma)
+
+.. c:function:: double random_rayleigh(bitgen_t *bitgen_state, double mode)
+
+.. c:function:: double random_standard_t(bitgen_t *bitgen_state, double df)
+
+.. c:function:: double random_noncentral_chisquare(bitgen_t *bitgen_state, double df,
+                                       double nonc)
+.. c:function:: double random_noncentral_f(bitgen_t *bitgen_state, double dfnum,
+                               double dfden, double nonc)
+.. c:function:: double random_wald(bitgen_t *bitgen_state, double mean, double scale)
+
+.. c:function:: double random_vonmises(bitgen_t *bitgen_state, double mu, double kappa)
+
+.. c:function:: double random_triangular(bitgen_t *bitgen_state, double left, double mode,
+                             double right)
+
+.. c:function:: int64_t random_poisson(bitgen_t *bitgen_state, double lam)
+
+.. c:function:: int64_t random_negative_binomial(bitgen_t *bitgen_state, double n, double p)
+
+.. c:function:: int64_t random_binomial(bitgen_t *bitgen_state, double p, int64_t n, binomial_t *binomial)
+
+.. c:function:: int64_t random_logseries(bitgen_t *bitgen_state, double p)
+
+.. c:function:: int64_t random_geometric_search(bitgen_t *bitgen_state, double p)
+
+.. c:function:: int64_t random_geometric_inversion(bitgen_t *bitgen_state, double p)
+
+.. c:function:: int64_t random_geometric(bitgen_t *bitgen_state, double p)
+
+.. c:function:: int64_t random_zipf(bitgen_t *bitgen_state, double a)
+
+.. c:function:: int64_t random_hypergeometric(bitgen_t *bitgen_state, int64_t good, int64_t bad,
+                                    int64_t sample)
+
+.. c:function:: uint64_t random_interval(bitgen_t *bitgen_state, uint64_t max)
+
+.. c:function:: void random_multinomial(bitgen_t *bitgen_state, int64_t n, int64_t *mnix,
+                            double *pix, np.npy_intp d, binomial_t *binomial)
+
+.. c:function:: int random_mvhg_count(bitgen_t *bitgen_state,
+                          int64_t total,
+                          size_t num_colors, int64_t *colors,
+                          int64_t nsample,
+                          size_t num_variates, int64_t *variates)
+
+.. c:function:: void random_mvhg_marginals(bitgen_t *bitgen_state,
+                               int64_t total,
+                               size_t num_colors, int64_t *colors,
+                               int64_t nsample,
+                               size_t num_variates, int64_t *variates)
+
+Generate a single integer
+
+.. c:function:: int64_t random_positive_int64(bitgen_t *bitgen_state)
+
+.. c:function:: int32_t random_positive_int32(bitgen_t *bitgen_state)
+
+.. c:function:: int64_t random_positive_int(bitgen_t *bitgen_state)
+
+.. c:function:: uint64_t random_uint(bitgen_t *bitgen_state)
+
+
+Generate random uint64 numbers in closed interval [off, off + rng].
+
+.. c:function:: uint64_t random_bounded_uint64(bitgen_t *bitgen_state,
+                                   uint64_t off, uint64_t rng,
+                                   uint64_t mask, bint use_masked)
+
+
diff --git a/doc/source/reference/random/examples/cython/extending.pyx b/doc/source/reference/random/examples/cython/extending.pyx
new file mode 100644
index 000000000..0cfbc146f
--- /dev/null
+++ b/doc/source/reference/random/examples/cython/extending.pyx
@@ -0,0 +1,4 @@
+extending.pyx
+-------------
+
+.. include:: ../../../../../../numpy/random/examples/extending.pyx
diff --git a/doc/source/reference/random/examples/cython/extending.pyx.rst b/doc/source/reference/random/examples/cython/extending.pyx.rst
new file mode 100644
index 000000000..bc31488d7
--- /dev/null
+++ b/doc/source/reference/random/examples/cython/extending.pyx.rst
@@ -0,0 +1,5 @@
+extending.pyx
+-------------
+
+.. literalinclude:: ../../../../../../numpy/random/examples/cython/extending.pyx
+    :language: cython
diff --git a/doc/source/reference/random/examples/cython/extending_distributions.pyx.rst b/doc/source/reference/random/examples/cython/extending_distributions.pyx.rst
new file mode 100644
index 000000000..a1bb01f45
--- /dev/null
+++ b/doc/source/reference/random/examples/cython/extending_distributions.pyx.rst
@@ -0,0 +1,5 @@
+extending_distributions.pyx
+---------------------------
+
+.. literalinclude:: ../../../../../../numpy/random/examples/cython/extending_distributions.pyx
+    :language: cython
diff --git a/doc/source/reference/random/examples/cython/index.rst b/doc/source/reference/random/examples/cython/index.rst
new file mode 100644
index 000000000..9c5da9559
--- /dev/null
+++ b/doc/source/reference/random/examples/cython/index.rst
@@ -0,0 +1,9 @@
+
+Extending `numpy.random` via Cython
+-----------------------------------
+
+
+.. toctree::
+    setup.py.rst
+    extending.pyx
+    extending_distributions.pyx
diff --git a/doc/source/reference/random/examples/cython/setup.py.rst b/doc/source/reference/random/examples/cython/setup.py.rst
new file mode 100644
index 000000000..381b45fd4
--- /dev/null
+++ b/doc/source/reference/random/examples/cython/setup.py.rst
@@ -0,0 +1,5 @@
+setup.py
+--------
+
+.. literalinclude:: ../../../../../../numpy/random/examples/cython/setup.py
+    :language: python
diff --git a/doc/source/reference/random/examples/numba.rst b/doc/source/reference/random/examples/numba.rst
new file mode 100644
index 000000000..a780afde7
--- /dev/null
+++ b/doc/source/reference/random/examples/numba.rst
@@ -0,0 +1,5 @@
+Extending via Numba
+-------------------
+
+.. literalinclude:: ../../../../../numpy/random/examples/numba/extending.py
+    :language: python
diff --git a/doc/source/reference/random/examples/numba_cffi.rst b/doc/source/reference/random/examples/numba_cffi.rst
new file mode 100644
index 000000000..ad4767a7a
--- /dev/null
+++ b/doc/source/reference/random/examples/numba_cffi.rst
@@ -0,0 +1,5 @@
+Extending via Numba and CFFI
+----------------------------
+
+.. literalinclude:: ../../../../../numpy/random/examples/numba/extending_distributions.py
+    :language: python
diff --git a/doc/source/reference/random/extending.rst b/doc/source/reference/random/extending.rst
index 22f9cb7e4..c63fb1a1b 100644
--- a/doc/source/reference/random/extending.rst
+++ b/doc/source/reference/random/extending.rst
@@ -12,128 +12,42 @@ Numba
 Numba can be used with either CTypes or CFFI.  The current iteration of the
 BitGenerators all export a small set of functions through both interfaces.
 
-This example shows how numba can be used to produce Box-Muller normals using
+This example shows how numba can be used to produce gaussian samples using
 a pure Python implementation which is then compiled.  The random numbers are
 provided by ``ctypes.next_double``.
 
-.. code-block:: python
-
-    from numpy.random import PCG64
-    import numpy as np
-    import numba as nb
-
-    x = PCG64()
-    f = x.ctypes.next_double
-    s = x.ctypes.state
-    state_addr = x.ctypes.state_address
-
-    def normals(n, state):
-        out = np.empty(n)
-        for i in range((n+1)//2):
-            x1 = 2.0*f(state) - 1.0
-            x2 = 2.0*f(state) - 1.0
-            r2 = x1*x1 + x2*x2
-            while r2 >= 1.0 or r2 == 0.0:
-                x1 = 2.0*f(state) - 1.0
-                x2 = 2.0*f(state) - 1.0
-                r2 = x1*x1 + x2*x2
-            g = np.sqrt(-2.0*np.log(r2)/r2)
-            out[2*i] = g*x1
-            if 2*i+1 < n:
-                out[2*i+1] = g*x2
-        return out
-
-    # Compile using Numba
-    print(normals(10, s).var())
-    # Warm up
-    normalsj = nb.jit(normals, nopython=True)
-    # Must use state address not state with numba
-    normalsj(1, state_addr)
-    %timeit normalsj(1000000, state_addr)
-    print('1,000,000 Box-Muller (numba/PCG64) randoms')
-    %timeit np.random.standard_normal(1000000)
-    print('1,000,000 Box-Muller (NumPy) randoms')
-
+.. literalinclude:: ../../../../numpy/random/examples/numba/extending.py
+    :language: python
+    :end-before: example 2
 
 Both CTypes and CFFI allow the more complicated distributions to be used
-directly in Numba after compiling the file distributions.c into a DLL or so.
-An example showing the use of a more complicated distribution is in the
-examples folder.
+directly in Numba after compiling the file distributions.c into a ``DLL`` or
+``so``.  An example showing the use of a more complicated distribution is in
+the `examples` section below.
 
-.. _randomgen_cython:
+.. _random_cython:
 
 Cython
 ======
 
 Cython can be used to unpack the ``PyCapsule`` provided by a BitGenerator.
-This example uses `~pcg64.PCG64` and
-``random_gauss_zig``, the Ziggurat-based generator for normals, to fill an
-array.  The usual caveats for writing high-performance code using Cython --
-removing bounds checks and wrap around, providing array alignment information
--- still apply.
-
-.. code-block:: cython
-
-    import numpy as np
-    cimport numpy as np
-    cimport cython
-    from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
-    from numpy.random.common cimport *
-    from numpy.random.distributions cimport random_gauss_zig
-    from numpy.random import PCG64
-
-
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    def normals_zig(Py_ssize_t n):
-        cdef Py_ssize_t i
-        cdef bitgen_t *rng
-        cdef const char *capsule_name = "BitGenerator"
-        cdef double[::1] random_values
-
-        x = PCG64()
-        capsule = x.capsule
-        if not PyCapsule_IsValid(capsule, capsule_name):
-            raise ValueError("Invalid pointer to anon_func_state")
-        rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name)
-        random_values = np.empty(n)
-        # Best practice is to release GIL and acquire the lock
-        with x.lock, nogil:
-            for i in range(n):
-                random_values[i] = random_gauss_zig(rng)
-        randoms = np.asarray(random_values)
-        return randoms
+This example uses `PCG64` and the example from above.  The usual caveats
+for writing high-performance code using Cython -- removing bounds checks and
+wrap around, providing array alignment information -- still apply.
+
+.. literalinclude:: ../../../../numpy/random/examples/cython/extending_distributions.pyx
+    :language: cython
+    :end-before: example 2
 
 The BitGenerator can also be directly accessed using the members of the basic
 RNG structure.
 
-.. code-block:: cython
-
-    @cython.boundscheck(False)
-    @cython.wraparound(False)
-    def uniforms(Py_ssize_t n):
-        cdef Py_ssize_t i
-        cdef bitgen_t *rng
-        cdef const char *capsule_name = "BitGenerator"
-        cdef double[::1] random_values
-
-        x = PCG64()
-        capsule = x.capsule
-        # Optional check that the capsule if from a BitGenerator
-        if not PyCapsule_IsValid(capsule, capsule_name):
-            raise ValueError("Invalid pointer to anon_func_state")
-        # Cast the pointer
-        rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name)
-        random_values = np.empty(n)
-        with x.lock, nogil:
-            for i in range(n):
-                # Call the function
-                random_values[i] = rng.next_double(rng.state)
-        randoms = np.asarray(random_values)
-        return randoms
+.. literalinclude:: ../../../../numpy/random/examples/cython/extending_distributions.pyx
+    :language: cython
+    :start-after: example 2
 
 These functions along with a minimal setup file are included in the
-examples folder.
+`examples` folder, ``numpy.random.examples``.
 
 New Basic RNGs
 ==============
@@ -163,3 +77,11 @@ the next 64-bit unsigned integer function if not needed. Functions inside
 .. code-block:: c
 
   bitgen_state->next_uint64(bitgen_state->state)
+
+Examples
+========
+
+.. toctree::
+    Numba <examples/numba>
+    CFFI + Numba <examples/numba_cffi> 
+    Cython <examples/cython/index>
diff --git a/doc/source/reference/random/index.rst b/doc/source/reference/random/index.rst
index 15c161244..d28646df9 100644
--- a/doc/source/reference/random/index.rst
+++ b/doc/source/reference/random/index.rst
@@ -32,7 +32,7 @@ instance's methods are imported into the numpy.random namespace, see
 Quick Start
 -----------
 
-By default, `~Generator` uses bits provided by `~pcg64.PCG64` which
+By default, `~Generator` uses bits provided by `PCG64` which
 has better statistical properties than the legacy mt19937 random
 number generator in `~.RandomState`.
 
@@ -155,7 +155,7 @@ What's New or Different
   (`~.PCG64.ctypes`) and CFFI (`~.PCG64.cffi`). This allows the bit generators
   to be used in numba.
 * The bit generators can be used in downstream projects via
-  :ref:`Cython <randomgen_cython>`.
+  :ref:`Cython <random_cython>`.
 * `~.Generator.integers` is now the canonical way to generate integer
   random numbers from a discrete uniform distribution. The ``rand`` and
   ``randn`` methods are only available through the legacy `~.RandomState`.
@@ -199,7 +199,8 @@ Features
    Multithreaded Generation <multithreading>
    new-or-different
    Comparing Performance <performance>
-   extending
+   c-api
+   Examples of using Numba, Cython, CFFI <extending>
 
 Original Source
 ~~~~~~~~~~~~~~~
diff --git a/numpy/random/include/bitgen.h b/numpy/core/include/numpy/random/bitgen.h
index 83c2858dd..83c2858dd 100644
--- a/numpy/random/include/bitgen.h
+++ b/numpy/core/include/numpy/random/bitgen.h
diff --git a/numpy/random/include/distributions.h b/numpy/core/include/numpy/random/distributions.h
index c02ea605e..e489e69b8 100644
--- a/numpy/random/include/distributions.h
+++ b/numpy/core/include/numpy/random/distributions.h
@@ -8,7 +8,7 @@
 #include <stdint.h>
 
 #include "numpy/npy_math.h"
-#include "include/bitgen.h"
+#include "numpy/random/bitgen.h"
 
 /*
  * RAND_INT_TYPE is used to share integer generators with RandomState which
diff --git a/numpy/random/_bit_generator.pxd b/numpy/random/_bit_generator.pxd
index 30fa4a27d..bd5e47a20 100644
--- a/numpy/random/_bit_generator.pxd
+++ b/numpy/random/_bit_generator.pxd
@@ -1,7 +1,7 @@
 cimport numpy as np
 from libc.stdint cimport uint32_t, uint64_t
 
-cdef extern from "include/bitgen.h":
+cdef extern from "numpy/random/bitgen.h":
     struct bitgen:
         void *state
         uint64_t (*next_uint64)(void *st) nogil
diff --git a/numpy/random/_bounded_integers.pyx.in b/numpy/random/_bounded_integers.pyx.in
index c0068dab2..7e19471e4 100644
--- a/numpy/random/_bounded_integers.pyx.in
+++ b/numpy/random/_bounded_integers.pyx.in
@@ -8,7 +8,7 @@ __all__ = []
 
 np.import_array()
 
-cdef extern from "include/distributions.h":
+cdef extern from "numpy/random/distributions.h":
     # Generate random numbers in closed interval [off, off + rng].
     uint64_t random_bounded_uint64(bitgen_t *bitgen_state,
                                    uint64_t off, uint64_t rng,
diff --git a/numpy/random/_generator.pyx b/numpy/random/_generator.pyx
index b6c222cc0..39d05b32b 100644
--- a/numpy/random/_generator.pyx
+++ b/numpy/random/_generator.pyx
@@ -28,7 +28,7 @@ from ._common cimport (POISSON_LAM_MAX, CONS_POSITIVE, CONS_NONE,
         )
 
 
-cdef extern from "include/distributions.h":
+cdef extern from "numpy/random/distributions.h":
 
     struct s_binomial_t:
         int has_binomial
diff --git a/numpy/random/examples/__init__.py b/numpy/random/examples/__init__.py
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/numpy/random/examples/__init__.py
diff --git a/numpy/random/examples/cython/extending.pyx b/numpy/random/examples/cython/extending.pyx
index a6a4ba4bf..2a866648d 100644
--- a/numpy/random/examples/cython/extending.pyx
+++ b/numpy/random/examples/cython/extending.pyx
@@ -8,7 +8,7 @@ import numpy as np
 cimport numpy as np
 cimport cython
 
-from numpy.random.common cimport bitgen_t
+from numpy.random._bit_generator cimport bitgen_t
 from numpy.random import PCG64
 
 np.import_array()
diff --git a/numpy/random/examples/cython/extending_distributions.pyx b/numpy/random/examples/cython/extending_distributions.pyx
index 3cefec97e..d17da45c1 100644
--- a/numpy/random/examples/cython/extending_distributions.pyx
+++ b/numpy/random/examples/cython/extending_distributions.pyx
@@ -1,21 +1,23 @@
 #!/usr/bin/env python
 #cython: language_level=3
 """
-This file shows how the distributions that are accessed through
-distributions.pxd can be used Cython code.
+This file shows how the to use a BitGenerator to create a distribution.
 """
 import numpy as np
 cimport numpy as np
 cimport cython
 from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
-from numpy.random.common cimport *
-from numpy.random.distributions cimport random_gauss_zig
+from numpy.random._bit_generator cimport bitgen_t
 from numpy.random import PCG64
 
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def normals_zig(Py_ssize_t n):
+def uniforms(Py_ssize_t n):
+    """ Create an array of `n` uniformly distributed doubles.
+    A 'real' distribution would want to process the values into
+    some non-uniform distribution
+    """
     cdef Py_ssize_t i
     cdef bitgen_t *rng
     cdef const char *capsule_name = "BitGenerator"
@@ -23,21 +25,25 @@ def normals_zig(Py_ssize_t n):
 
     x = PCG64()
     capsule = x.capsule
+    # Optional check that the capsule if from a BitGenerator
     if not PyCapsule_IsValid(capsule, capsule_name):
         raise ValueError("Invalid pointer to anon_func_state")
+    # Cast the pointer
     rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name)
-    random_values = np.empty(n)
-    # Best practice is to release GIL and acquire the lock
+    random_values = np.empty(n, dtype='float64')
     with x.lock, nogil:
         for i in range(n):
-            random_values[i] = random_gauss_zig(rng)
+            # Call the function
+            random_values[i] = rng.next_double(rng.state)
     randoms = np.asarray(random_values)
-    return randoms
 
+    return randoms
+ 
+# cython example 2
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def uniforms(Py_ssize_t n):
+def uint16_uniforms(Py_ssize_t n):
     cdef Py_ssize_t i
     cdef bitgen_t *rng
     cdef const char *capsule_name = "BitGenerator"
@@ -45,15 +51,13 @@ def uniforms(Py_ssize_t n):
 
     x = PCG64()
     capsule = x.capsule
-    # Optional check that the capsule if from a BitGenerator
     if not PyCapsule_IsValid(capsule, capsule_name):
         raise ValueError("Invalid pointer to anon_func_state")
-    # Cast the pointer
     rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name)
-    random_values = np.empty(n)
+    random_values = np.empty(n, dtype='uint32')
+    # Best practice is to release GIL and acquire the lock
     with x.lock, nogil:
         for i in range(n):
-            # Call the function
-            random_values[i] = rng.next_double(rng.state)
+            random_values[i] = rng.next_uint32(rng.state)
     randoms = np.asarray(random_values)
     return randoms
diff --git a/numpy/random/examples/cython/setup.py b/numpy/random/examples/cython/setup.py
index 69f057ed5..315527a2d 100644
--- a/numpy/random/examples/cython/setup.py
+++ b/numpy/random/examples/cython/setup.py
@@ -9,15 +9,20 @@ import numpy as np
 from distutils.core import setup
 from Cython.Build import cythonize
 from setuptools.extension import Extension
-from os.path import join
+from os.path import join, abspath, dirname
+
+curpath = abspath(dirname(__file__))
 
 extending = Extension("extending",
-                      sources=['extending.pyx'],
-                      include_dirs=[np.get_include()])
+                      sources=[join(curpath, 'extending.pyx')],
+                      include_dirs=[
+                            np.get_include(),
+                            join(curpath, '..', '..')
+                        ],
+                      )
 distributions = Extension("extending_distributions",
-                          sources=['extending_distributions.pyx',
-                                   join('..', '..', 'src',
-                                        'distributions', 'distributions.c')],
+                          sources=[join(curpath, 'extending_distributions.pyx'),
+                                  ],
                           include_dirs=[np.get_include()])
 
 extensions = [extending, distributions]
diff --git a/numpy/random/examples/numba/extending.py b/numpy/random/examples/numba/extending.py
index d41c2d76f..0d240596b 100644
--- a/numpy/random/examples/numba/extending.py
+++ b/numpy/random/examples/numba/extending.py
@@ -1,14 +1,57 @@
-import datetime as dt
-
 import numpy as np
 import numba as nb
 
 from numpy.random import PCG64
+from timeit import timeit
+
+bit_gen = PCG64()
+next_d = bit_gen.cffi.next_double
+state_addr = bit_gen.cffi.state_address
+
+def normals(n, state):
+    out = np.empty(n)
+    for i in range((n + 1) // 2):
+        x1 = 2.0 * next_d(state) - 1.0
+        x2 = 2.0 * next_d(state) - 1.0
+        r2 = x1 * x1 + x2 * x2
+        while r2 >= 1.0 or r2 == 0.0:
+            x1 = 2.0 * next_d(state) - 1.0
+            x2 = 2.0 * next_d(state) - 1.0
+            r2 = x1 * x1 + x2 * x2
+        f = np.sqrt(-2.0 * np.log(r2) / r2)
+        out[2 * i] = f * x1
+        if 2 * i + 1 < n:
+            out[2 * i + 1] = f * x2
+    return out
+
+# Compile using Numba
+normalsj = nb.jit(normals, nopython=True)
+# Must use state address not state with numba
+n = 10000
+
+def numbacall():
+    return normalsj(n, state_addr)
+
+rg = np.random.Generator(PCG64())
 
-x = PCG64()
-f = x.ctypes.next_uint32
-s = x.ctypes.state
+def numpycall():
+    return rg.normal(size=n)
 
+# Check that the functions work
+r1 = numbacall()
+r2 = numpycall()
+assert r1.shape == (n,)
+assert r1.shape == r2.shape
+
+t1 = timeit(numbacall, number=1000)
+print('{:.2f} secs for {} PCG64 (Numba/PCG64) gaussian randoms'.format(t1, n))
+t2 = timeit(numpycall, number=1000)
+print('{:.2f} secs for {} PCG64 (NumPy/PCG64) gaussian randoms'.format(t2, n))
+
+# example 2
+
+next_u32 = bit_gen.ctypes.next_uint32
+ctypes_state = bit_gen.ctypes.state
 
 @nb.jit(nopython=True)
 def bounded_uint(lb, ub, state):
@@ -19,14 +62,14 @@ def bounded_uint(lb, ub, state):
     mask |= mask >> 8
     mask |= mask >> 16
 
-    val = f(state) & mask
+    val = next_u32(state) & mask
     while val > delta:
-        val = f(state) & mask
+        val = next_u32(state) & mask
 
     return lb + val
 
 
-print(bounded_uint(323, 2394691, s.value))
+print(bounded_uint(323, 2394691, ctypes_state.value))
 
 
 @nb.jit(nopython=True)
@@ -36,42 +79,6 @@ def bounded_uints(lb, ub, n, state):
         out[i] = bounded_uint(lb, ub, state)
 
 
-bounded_uints(323, 2394691, 10000000, s.value)
-
-g = x.cffi.next_double
-cffi_state = x.cffi.state
-state_addr = x.cffi.state_address
-
-
-def normals(n, state):
-    out = np.empty(n)
-    for i in range((n + 1) // 2):
-        x1 = 2.0 * g(state) - 1.0
-        x2 = 2.0 * g(state) - 1.0
-        r2 = x1 * x1 + x2 * x2
-        while r2 >= 1.0 or r2 == 0.0:
-            x1 = 2.0 * g(state) - 1.0
-            x2 = 2.0 * g(state) - 1.0
-            r2 = x1 * x1 + x2 * x2
-        f = np.sqrt(-2.0 * np.log(r2) / r2)
-        out[2 * i] = f * x1
-        if 2 * i + 1 < n:
-            out[2 * i + 1] = f * x2
-    return out
+bounded_uints(323, 2394691, 10000000, ctypes_state.value)
 
 
-print(normals(10, cffi_state).var())
-# Warm up
-normalsj = nb.jit(normals, nopython=True)
-normalsj(1, state_addr)
-
-start = dt.datetime.now()
-normalsj(1000000, state_addr)
-ms = 1000 * (dt.datetime.now() - start).total_seconds()
-print('1,000,000 Polar-transform (numba/PCG64) randoms in '
-      '{ms:0.1f}ms'.format(ms=ms))
-
-start = dt.datetime.now()
-np.random.standard_normal(1000000)
-ms = 1000 * (dt.datetime.now() - start).total_seconds()
-print('1,000,000 Polar-transform (NumPy) randoms in {ms:0.1f}ms'.format(ms=ms))
diff --git a/numpy/random/include/legacy-distributions.h b/numpy/random/include/legacy-distributions.h
index 6a0fc7dc4..b8ba0841c 100644
--- a/numpy/random/include/legacy-distributions.h
+++ b/numpy/random/include/legacy-distributions.h
@@ -2,7 +2,7 @@
 #define _RANDOMDGEN__DISTRIBUTIONS_LEGACY_H_
 
 
-#include "distributions.h"
+#include "numpy/random/distributions.h"
 
 typedef struct aug_bitgen {
   bitgen_t *bit_generator;
diff --git a/numpy/random/mtrand.pyx b/numpy/random/mtrand.pyx
index b21607282..691a6e6e7 100644
--- a/numpy/random/mtrand.pyx
+++ b/numpy/random/mtrand.pyx
@@ -25,7 +25,7 @@ from ._common cimport (POISSON_LAM_MAX, CONS_POSITIVE, CONS_NONE,
             check_array_constraint, check_constraint, disc, discrete_broadcast_iii,
         )
 
-cdef extern from "include/distributions.h":
+cdef extern from "numpy/random/distributions.h":
     struct s_binomial_t:
         int has_binomial
         double psave
diff --git a/numpy/random/setup.py b/numpy/random/setup.py
index f9059d7d7..93d5144ea 100644
--- a/numpy/random/setup.py
+++ b/numpy/random/setup.py
@@ -34,6 +34,7 @@ def configuration(parent_package='', top_path=None):
 
     defs.append(('NPY_NO_DEPRECATED_API', 0))
     config.add_data_dir('tests')
+    config.add_data_dir('examples')
 
     EXTRA_LINK_ARGS = []
     # Math lib
@@ -92,6 +93,7 @@ def configuration(parent_package='', top_path=None):
                              depends=['%s.pyx' % gen, '%s.pxd' % gen,],
                              define_macros=defs,
                              )
+        config.add_data_files('{0}.pxd'.format(gen))
     other_srcs = [
         'src/distributions/logfactorial.c',
         'src/distributions/distributions.c',
@@ -110,6 +112,7 @@ def configuration(parent_package='', top_path=None):
                              depends=['%s.pyx' % gen],
                              define_macros=defs,
                              )
+    config.add_data_files('_bounded_inteters.pxd')
     config.add_extension('mtrand',
                          sources=['mtrand.c',
                                   'src/legacy/legacy-distributions.c',
diff --git a/numpy/random/src/distributions/distributions.c b/numpy/random/src/distributions/distributions.c
index ab8de8bcb..df3323408 100644
--- a/numpy/random/src/distributions/distributions.c
+++ b/numpy/random/src/distributions/distributions.c
@@ -1,4 +1,4 @@
-#include "include/distributions.h"
+#include "numpy/random/distributions.h"
 #include "ziggurat_constants.h"
 #include "logfactorial.h"
 
diff --git a/numpy/random/src/distributions/random_hypergeometric.c b/numpy/random/src/distributions/random_hypergeometric.c
index da5ea9c68..0da49bd62 100644
--- a/numpy/random/src/distributions/random_hypergeometric.c
+++ b/numpy/random/src/distributions/random_hypergeometric.c
@@ -1,4 +1,4 @@
-#include "include/distributions.h"
+#include "numpy/random/distributions.h"
 #include "logfactorial.h"
 #include <stdint.h>
 
diff --git a/numpy/random/src/distributions/random_mvhg_count.c b/numpy/random/src/distributions/random_mvhg_count.c
index 9c0cc045d..0c46ea417 100644
--- a/numpy/random/src/distributions/random_mvhg_count.c
+++ b/numpy/random/src/distributions/random_mvhg_count.c
@@ -2,7 +2,7 @@
 #include <stdlib.h>
 #include <stdbool.h>
 
-#include "include/distributions.h"
+#include "numpy/random/distributions.h"
 
 /*
  *  random_mvhg_count
diff --git a/numpy/random/src/distributions/random_mvhg_marginals.c b/numpy/random/src/distributions/random_mvhg_marginals.c
index 301a4acad..7e4c24988 100644
--- a/numpy/random/src/distributions/random_mvhg_marginals.c
+++ b/numpy/random/src/distributions/random_mvhg_marginals.c
@@ -3,7 +3,7 @@
 #include <stdbool.h>
 #include <math.h>
 
-#include "include/distributions.h"
+#include "numpy/random/distributions.h"
 #include "logfactorial.h"
 
 
diff --git a/numpy/random/tests/test_extending.py b/numpy/random/tests/test_extending.py
new file mode 100644
index 000000000..01e43b9e8
--- /dev/null
+++ b/numpy/random/tests/test_extending.py
@@ -0,0 +1,32 @@
+import os, sys
+import pytest
+import warnings
+
+try:
+    with warnings.catch_warnings(record=True) as w:
+        # numba issue gh-4733
+        warnings.filterwarnings('always', '', DeprecationWarning)
+        import numba
+    import cffi
+except ImportError:
+    numba = None
+
+def test_cython():
+    curdir = os.getcwd()
+    argv = sys.argv
+    examples = (os.path.dirname(__file__), '..', 'examples')
+    try:
+        os.chdir(os.path.join(*examples))
+        sys.argv = argv[:1] + ['build']
+        with warnings.catch_warnings(record=True) as w:
+            # setuptools issue gh-1885
+            warnings.filterwarnings('always', '', DeprecationWarning)
+            from numpy.random.examples.cython import setup
+    finally:
+        sys.argv = argv
+        os.chdir(curdir)
+
+@pytest.mark.skipif(numba is None, reason="requires numba")
+def test_numba():
+        from numpy.random.examples.numba import extending
+
diff --git a/numpy/tests/test_public_api.py b/numpy/tests/test_public_api.py
index c71d03432..ba3961bef 100644
--- a/numpy/tests/test_public_api.py
+++ b/numpy/tests/test_public_api.py
@@ -334,6 +334,7 @@ SKIP_LIST = [
     "numpy.core.cversions",
     "numpy.core.generate_numpy_api",
     "numpy.distutils.msvc9compiler",
+    'numpy.random.examples',
 ]
 
 
@@ -386,7 +387,7 @@ SKIP_LIST_2 = [
     'numpy.matlib.fft',
     'numpy.matlib.random',
     'numpy.matlib.ctypeslib',
-    'numpy.matlib.ma'
+    'numpy.matlib.ma',
 ]