MAINT: Get full precision for 32 bit floating point random values.

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

    (((rng)->next_uint32((rng)->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

    (((rng)->next_uint32((rng)->state) >> 8) * (1.0f / 16777216.0f))

See numpy/numpy#20314 for more details.

src/pgm_macros.h

@@ -52,7 +52,7 @@
  * adapted from a private <numpy/random/distributions.h> function of a similar name
  */
 #define next_float(rng) \
-    (((rng)->next_uint32((rng)->state) >> 9) * (1.0f / 8388608.0f))
+    (((rng)->next_uint32((rng)->state) >> 8) * (1.0f / 16777216.0f))
 
 /*
  * Generate a random double precision float in the range [0, 1). This macros is