Currently, `gsl_rng_uniform is generator dependent and is in most cases implemented by dividing the output of `gsl_rng_get()` by `gsl_rng_max() + 1`, as stated in the docs. This has an obvious problem, the resolution of the generated floating points is solely dependent on the generators range.
None of the currently implemented generators have a range larger than 2^32 and most are smaller, which means that we can never even generate a random double to double precision. This issue is a lot bigger as many random distributions use `gsl_rng_uniform` internally and/or only return up to single-precision floating-point precision, as is the case with `gsl_ran_gaussian_ziggurat`. But if we want to generate a random double-precision floating-point number in the range [0,1), how many bits would be required? First, let's consider that we can't generate random floats uniformly using division by an integer: Integer: |----|----|----|----|----|----|----|----|----|----|----|----| | | / | / / / / \ __/ \ __/ | Float: ||||-|-|-|--|--|--|----|----|----|--------|--------|--------| There is no 1 to 1 mapping between an integer and a floating-point number. Multiple integer values are mapped to the same floating-point value and/or smaller floating-point values are skipped. (see the bottom of http://prng.di.unimi.it/ and the github link below) One method for generating random floats is to divide by the biggest number that doesn't map multiple integer values to the same floating-point number. That would be 2^53 for doubles. This is a good compromise between accuracy and performance and should be the default for random double generation. Hence, I propose adding a floating-point distribution function gsl_rand_uniform (and maybe gsl_rand_uniformf) that guarantee this precision for all generators. To make implementing this easier an additional generic interface could be added to each RNG: `void gsl_rng_write(gsl_rng *r, void *ptr, size_t n)`, that writes to `n` bytes of random data to `ptr`. There is another method for generating random floating-point values, that includes all representable numbers inside [a,b]. I won't explain it here, but I've written an implementation of this in https://github.com/camel-cdr/cauldron/blob/main/cauldron/random.h#L1353 and this could be added as something like `double gsl_rand_uniform_dense(gsl_rng *r, double min, double max)`. I'd like some feedback on these proposed changes, so I can write a patch in the future. Best regards, Olaf