[Help-gsl] Random number generator maximum values (and 64-bit)

[Joseph Wakeling]

Hello all,

A couple of days ago I installed a 64-bit Linux distro instead of a
32-bit one, and discovered an interesting bug in some C++ code I'd written.

The bug comes down to this little bit of code:

    unsigned long int M_ul = floor(log2(ULONG_MAX-1))-1;
    unsigned long int M_rng = floor(log2(r->type->max-1));
    unsigned long int M = min(M_ul,M_rng);
    unsigned long int N = 1<<M;

... where min() is from the <algorithm> header and floor() is from <cmath>.

The aim here is that N should be the largest power of 2 such that

  -- it is within the range of the random number generator, and

  -- 2*N is within the range of the unsigned long int implementation
     of the machine.

This is important because later on in the code I will call,

    gsl_rng_uniform_int(rng,N);

... which is where the code, compiled on my 64-bit system, falls over
with an error,

    gsl: ../gsl/gsl_rng.h:200: ERROR: invalid n, either 0 or exceeds
    maximum value of generator

Further investigation reveals that despite the bit of code to make sure
N is within the desired limits, N is being set to 2^64, which is WAY
beyond the maximum range of the RNG.  (I'm using mt19937; incidentally,
the reason the above code got written in the first place is because once
I used ranlux, which has a maximum value less than the range of a 32-bit
unsigned long int.)

More weirdly still, M = min(M_ul,M_rng) is being set to the correct
value of 31.

I suspect the solution to this problem is not really gsl-related, but
thought I'd post the query here as a likely source of insight and
interest ... :-)

Thanks & best wishes,

    -- Joe