Most compilers should optimize it, but it doesn't hurt to inline and has a better name
now.

Assert that an invalid value (2**n-1 for n = 8, 16, 32, 64) has not
been passed to the Lemire function.

When the input to Generator.integers was 2**32, the value 2**32-1
was being passed as the `rng` argument to the 32-bit Lemire method,
but that method requires `rng` be strictly less then 2**32-1.

The fix was to handle 2**32-1 by calling next_uint32 directly.
This also works for the legacy code without changing the stream
of random integers from `randint`.

Closes gh-16066.

* 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

When an 8 or 16 bit dtype was given to the integers() method of the
Generator class, the resulting sample was biased.  The problem was
the lines of the form

    const uint8_t threshold = -rng_excl % rng_excl;

in the implementations of Lemire's method, in the C file
distributions.c.  The intent was to compute
    (UINT8_MAX+1 - rng_excl) % rng_excl
However, when the type of rng_excl has integer conversion rank lower
than a C int (which is almost certainly the case for the 8 and 16
bit types), the terms in the expression -rng_excl % rng_excl are
promoted to int, and the result of the calculation is always 0.

The fix is to make the expression explicit, and write it as

    const uint8_t threshold = (UINT8_MAX - rng) % rng_excl;

rng is used, because rng_excl is simply rng + 1; by using rng, we
we only need the constant UINT#_MAX, without the extra +1.

For consistency, I made the same change for all the data types
(8, 16, 32 and 64 bit).

Closes gh-14774.

The new method

  multivariate_hypergeometric(self, object colors, object nsample,
                              size=None, method='marginals')

of the class numpy.random.Generator implements the multivariate
hypergeometric distribution; see
  https://en.wikipedia.org/wiki/Hypergeometric_distribution,
specifically the section "Multivariate hypergeometric distribution".

Two algorithms are implemented.  The user selects which algorithm
to use with the `method` parameter. The default, `method='marginals'`,
is based on repeated calls of the univariate hypergeometric
distribution function.  The other algorithm, selected with
`method='count'`, is a brute-force method that allocates an
internal array of length ``sum(colors)``.  It should only be used
when that value is small, but it can be much faster than the
"marginals" algorithm in that case.

The C implementations of the two methods are in the files
random_mvhg_count.c and random_mvhg_marginals.c in
numpy/random/src/distributions.