Specifically the 'Logarithmic distribution' link is fixed.
For others the links or article names are updated.

Although the arguments are specified as 'high'
and 'low', it is possible to pass in values for
'low' and 'high' where 'low' >= 'high' and still
obtain well-defined behavior. The documentation
has been expanded to reflect this fact, with a
note to discourage passing in arguments satisfying
'low' > 'high'.

Fixes bug in np.random methods that would return scalars
when passed one-element array inputs. This is because
one-element ndarrays can be cast to integers / floats, which
is what functions like PyFloat_AsDouble do before converting
to the intended data type.

This commit changes the check used to determine whether the
inputs are purely scalar by converting all inputs to arrays
and checking if the resulting shape is an empty tuple (scalar)
or not (array).

Closes gh-4263.

Clarified the output size depending on whether
scalar or non-scalar inputs are passed to
functions in mtrand.pyx that can broadcast
their arguments.

Added a whole new suite of tests to ensure that
functions in mtrand.pyx which are broadcastable
actually broadcast their arguments properly.

Refactored methods that broadcast arguments
together by finding additional common ground
between code in the if...else branches that
involved a size parameter being passed in.

Apparently gcc only specializes one branch (the last one) so I went for
another 33% performance increase (matching #6776) in what's likely the
most common use case.

Only for 1d-ndarrays exactly, as subtypes (e.g. masked arrays) may not
allow direct shuffle of the underlying buffer (in fact, the old
implementation destroyed the underlying values of masked arrays while
shuffling).

Also handles struct-containing-object 1d ndarrays properly.

See #6776 for an earlier, less general (but even faster: ~6x)
improvement attempt, #5514 for the original issue.