[Help-gsl] Random number generation

[Victor Stinner]

Hi,

I'm writing a library for pseudo random number generation called Hasard. The 
goal is to get a simple API with portable code. Someone told me that GSL 
already exists. I read the GSL API and source code, but I dislike that 
gsl_rng_uniform_int() function.

This function sets an error if user range [0; n] if bigger than the generator 
range [0; r->type->max -r->type->min]. The code is very very simple and 
should be very fast, but it doesn't work is user needs a big number. Eg. [0; 
2^40-1] whereas the generator uses range [0; 2^32-1].

---

I wrote a function using the GMP library trying to get an uniform distribution 
for any generator and user range. Pseudo code:
  # Random tick in [0; tick_max]
  # User: [0; n]
  base = tick_max + 1
  result = 0
  quotient = 0
  nb_digits = 10
  for i in 1..nb_digits:
     result = (result * base) + random_tick()
     quotient = (quotient * base) + tick_max
  # convert range [0; quotient] to range [0; n]
  result = result * (n + 1) / (quotient + 1)
  result = floor(result)

Problems:
 - nb_digits has fixed value, which should be wrong in some cases
 - GMP is not fast: it's not really a problem
 - i'm not sure that the distribution is uniform :-/

Questions:

- How can I compute nb_digits?

- Do you know other algorithm to generate a random number with 
  a (proved) uniform distribution (working with any range size)?

---

I think that quotient % (n + 1) have to be zero to get a perfect uniform 
distribution. It can takes "many loops" (or will be unlimited) until this is 
true.

Example A:
  tick_max = n + 1 
  base = tick_max + 1
  quotient % n will be: 1, 3, 7, ..., (2^k - 1)

Example B:
  tick_max = n - 1 
  base = tick_max + 1
  quotient % n will be: n-1, n-1, ..., n-1

In example B, quotient%n will never reach zero. In base A, it depends on n.

---

My library under development:
   http://haypo.hachoir.org/trac/wiki/hasard

Victor