In C, when the integer part of a double exceeds the maximum int64, casting
the double to int64_t is undefined behavior.  On some platforms, the result
is 2**64-1 and on others it is -2**64.  That results in rng.geometric()
returning -2**64 when p was sufficiently small on some platforms.

The fix is to never invoke undefined behavior by ensuring we don't attempt
to cast very large double values to int64_t.

Closes gh-22100

The formula to convert a 32 bit random integer to a random float32,

    (next_uint32(bitgen_state) >> 9) * (1.0f / 8388608.0f)

shifts by one bit too many, resulting in uniform float32 samples always
having a 0 in the least significant bit. The formula is corrected to

    (next_uint32(bitgen_state) >> 8) * (1.0f / 16777216.0f)

Occurrences of the incorrect formula in numpy/random/tests/test_direct.py
were also corrected.

Closes gh-17478.

* Fix typos

Use exponentials rather than generating exponentials via inversion

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

closes #17020

This reduces the number of cases when floating point precision makes the argument to `log` become `0`

Avoid conditional compilation and move old version to legacy_vonmises
Small clean up
Additional comments

Apply vonmises fix only to Generator
Add tests for correctness

closes #17378
closes #17275