…uments.

The changes in this pull request restrict the arguments ngood and nbad
of numpy.random.hypergeometric to be less than 2**30.  This "fix" (in the
sense of the old joke "Doctor, it hurts when I do *this*") makes the
following problems impossible to encounter:

* Python hangs with these calls (see numpygh-11443):

      np.random.hypergeometric(2**62-1, 2**62-1, 26, size=2)
      np.random.hypergeometric(2**62-1, 2**62-1, 11, size=12)

* Also reported in numpygh-11443: the distribution of the variates of

      np.random.hypergeometric(ngood, nbad, nsample)

  is not correct when nsample >= 10 and ngood + nbad is sufficiently large.
  I don't have a rigorous bound, but a histogram of the p-values of 1500
  runs of a G-test with size=5000000 in each run starts to look suspicious
  when ngood + nbad is approximately 2**37.  When the sum is greater than
  2**38, the distribution is clearly wrong, and, as reported in the pull
  request numpygh-9834, the following call returns all zeros:

      np.random.hypergeometric(2**55, 2**55, 10, size=20)

* The code does not check for integer overflow of ngood + nbad, so the
  following call

      np.random.hypergeometric(2**62, 2**62, 1)

  on a system where the default integer is 64 bit results in
  `ValueError: ngood + nbad < nsample`.

By restricting the two arguments to be less than 2**30, the values are
well within the space of inputs for which I have no evidence of the
function generating an incorrect distribution, and the overflow problem
is avoided for both 32-bit and 64-bit systems.

I replaced the test for hypergeometric(2**40 - 2, 2**40 - 2, 2**40 - 2) > 0
with hypegeometric(2**30 - 1, 2**30 - 1, 2**30 - 1) > 0.  The test was a
regression test for a bug in which the function would enter an infinite
loop when the arguments were sufficiently large.  The regression test for
the fix included the call hypergeometric(2**40 - 2, 2**40 - 2, 2**40 - 2)
on systems with 64 bit integers.  This call is now disallowed, so I add a
call with the maximum allowed values of ngood and nbad.

Summary of the changes:

* Move the functions random_hypergeometric_hyp, random_hypergeometric_hrua
  and random_hypergeometric from distributions.c to legacy-distributions.c.
  These are now the legacy implementation of hypergeometric.
* Add the files logfactorial.c and logfactorial.h,
  containing the function logfactorial(int64_t k).
* Add the files random_hypergeometric.c and random_hypergeometric.h,
  containing the function random_hypergeometric (the new implementation
  of the hypergeometric distribution).  See more details below.
* Fix two tests in numpy/random/tests/test_generator_mt19937.py that
  used values returned by the hypergeometric distribution.  The
  new implementation changes the stream, so those tests needed to
  be updated.
* Remove another test obviated by an added constraint on the arguments
  of hypergeometric.

Details of the rewrite:

If you carefully step through the old function rk_hypergeometric_hyp(),
you'll see that the end result is basically the same as the new function
hypergeometric_sample(), but the new function accomplishes the result with
just integers.  The floating point calculations in the old code caused
problems when the arguments were extremely large (explained in more detail
in the unmerged pull request numpy#9834).

The new version of hypergeometric_hrua() is a new translation of
Stadlober's ratio-of-uniforms algorithm for the hypergeometric
distribution.  It fixes a mistake in the old implementation that made the
method less efficient than it could be (see the details in the unmerged
pull request numpy#11138), and uses a faster
function for computing log(k!).

The HRUA algorithm suffers from loss of floating point precision when
the arguments are *extremely* large (see the comments in github issue
11443).  To avoid these problems, the arguments `ngood` and `nbad` of
hypergeometric must be less than 10**9.  This constraint obviates an
existing regression test that was run on systems with 64 bit long
integers, so that test was removed.