Re: GSL 1.6 random numbers (gauss.c and rng.c)

["Charles Karney" <ckarney at sarnoff dot com>]

 > From: James Theiler <jt@lanl.gov>
 > Subject: Re: GSL 1.6 random numbers (gauss.c and rng.c)
 > Date: Wed, 26 Jan 2005 10:52:03 -0700 (MST)
 > 
 > ] Now, I've got some comments on
 > ] 
 > ]     gsl_ran_gaussian_ratio_method
 > ]     gsl_rng_uniform_int
 > ] 
 > ] in GSL 1.6
 > ] 
 > ] (a) [MINOR NIT] In gsl_ran_gaussian_ratio_method, the last digit of the
 > ]     constant for sqrt(8/e) is wrong.  Worse still, the number you have
 > ]     is too small.  Rounding the result UP still gives the correct
 > ]     results, however (at the expense of a few extra rejects).  So you
 > ]     might replace
 > ] 
 > ]       x = 1.71552776992141359295 * (v - 0.5) / u;
 > ]     by
 > ]       x = 1.715527769921413593 * (v - 0.5) / u;
 > ]     or
 > ]       x = 1.71552777 * (v - 0.5) / u;
 > ] 
 > ]     In my implementation I round up after additional digits:
 > ] 
 > ]       1.71552776992141359296037928256
 > 
 > Hmmm, if x were used solely as the rejection criterion, then I would
 > be inclined to agree that there is no harm to accuracy (and utterly
 > negligible harm to efficiency) to use a rounded-up value for the
 > constant.
 > 
 > But it's a little more complicated since x is also used as part of the
 > returned variate.  You may still be right that it doesn't matter,
 > maybe everything just scales up, but it's not obvious to me.

If you check the description of the algorithm in Knuth or the original
paper, you will see that you are just trying to sample a point uniformly
in a funny shaped region.  This is done by choosing a point in a
rectangle enclosing the region and rejecting the point unless it's
inside the region.  The rectangle [0,1] x [-sqrt(2/e), sqrt(2/e)] snugly
encloses the region.  However [0,1] x [-1,1] would also "work" at some
cost in efficiency and precision.  Leva uses 1.7156 for his multiplier.

 > ] (b) [EFFICIENCY] Leva gives quadratic tests to perform in
 > ]     gsl_ran_gaussian_ratio_method which result in log being evaluated
 > ]     0.012 per normal deviate instead of 1.369 times.  (Knuth's tests
 > ]     result in 0.232 logs/deviate).  Depending on the relative speed of
 > ]     log vs other arithmetic operations, this might speed up
 > ]     gsl_ran_gaussian_ratio_method significantly.  The reference is
 > ] 
 > ]       J. L. Leva, A Fast Normal Random Number Generator,
 > ]       ACM TOMS 18, 449-453 and 454-455 (1992).
 > 
 > Worth looking into. In my experience, logs are pretty fast on machines
 > I use (eg, Pentiums), so that the actual gains in reducing their
 > number are often not as substantial as one might hope.  But I haven't
 > done any benchmarks lately, and besides, if this is a better
 > algorithm, we should at least look into it.

I had mixed results on a Pentium with gcc.  My previous implemention
used Knuth's bounds.  Compiling with -ffast-math, Leva gives a 2.5%
speedup.  Without -ffast-math, I got a 18% speedup.  (As I recall, the
Knuthian bounds gave a 25% speedup over the raw ratio method.)

Note that the coefficients in Leva's quadratic bounds are only given to
a few sig. figs.  This is OK -- see discussion above about the the
sqrt(8/e) const.  The only important thing is that the FINAL exit test x
* x <= -4.0 * log (u) is accurate.

One other comment on both gsl_ran_gaussian_ratio_method and
gsl_ran_gaussian.  I think you would be better off using
gsl_rng_uniform_pos consistently in both routines.  This has the
property that (gsl_rng_uniform_pos - 0.5) is symmetrically distributed
about zero.  With gsl_rng_uniform, the gaussian routines return negative
numbers slightly more often than positive ones.  (A tiny effect -- but
better to fix it now than to worry about simulation results 5 years down
the road.)

 > ] (c) [GENERALIZE?] In gsl_rng_uniform_int you quit if n > range.  However
 > ]     within that test you could allow the case n-1 == range.  Thus
 > ] 
 > ]       if (n > range) 
 > ] 	{
 > ] 	  if (n - 1 == range)
 > ] 	    return ((r->type->get) (r->state)) - offset;
 > ] 	  GSL_ERROR_VAL ("n exceeds maximum value of generator",
 > ] 			    GSL_EINVAL, 0) ;
 > ] 	}
 > ] 
 > ]     (Of course if range = maxint - 1, n would not be representable...)
 > ] 
 > 
 > I have to admit, I haven't looked at this particular function before.  
 > It is confusing to me; it looks like range is defined in a funny way;
 > eg if min=0, max=99, shouldn't the range be 100 ? Or, 
 > if we are going to define range=max-min, then we 
 > should define the scale differently
 > 
 >       unsigned long int scale = (range + 1)/n;
 > 
 > where the "+ 1" is added to what we have now.  Then the rest proceeds
 > correctly and even more efficiently.
 > 
 > ...except when r->type->max is the maximum unsigned long (ie,
 > 0xffffffffUL), which is the case for a lot of the generators.  Then
 > the "+ 1" causes an overflow.  

I suspect that this was merely to avoid overflow.  Here's my independent
C++ implemention of this functionality where I encountered the same
problem.  The comments tell the story...

Here MM = 2^30 and UInt30() is Knuth's recommended RNG (returning a
result in [0, MM))

// A random number in [0, 2^32).
unsigned long Random::Integer() {
    // xor 2 30-bit numbers to give 32-bit result
    return (UInt30() << 2) ^ UInt30();
}

// Return random number in [0,n)
unsigned long Random::Integer(const unsigned long n) throw() {
    if (n == 0)
        // Special case, equivalent to Integer(2^32).  (Otherwise n == 0
        // is meaningless.)
        return Integer();
    else if (n == 1)
        // Saves using up a random number.
        return 0;
    else {
        // Do we need 32-bit randoms as our raw material?
        const bool big = n > MM;
        // The maximum random that we need.  We'd like the first number
        // to be 2^32, but we can't deal with such large numbers.  The
        // following will give the correct results at a slight cost in
        // efficiency.  (The case n = 2^31 will need two trips through
        // the loop instead of one.)
        const unsigned long q = big ? 0xffffffffUL : MM;
        const unsigned long r = q / n; // r > 0
        const unsigned long m = r * n; // 0 < m <= q
        unsigned long u;        // Find a random number in [0,m).
        do {
            // For small n, this is executed once (since m is nearly q).
            // Worst cases are n = MM/2 + 1 and 2*MM, in which case the
            // loop is executed slightly less than twice on average.
            u = big ? Integer() : UInt30();
        } while (u >= m);
        // Since m is r * n, u/r is unbiased random number in [0,n).  An
        // unbiased alternative is to use u % n.  However, division is
        // preferred, as it uses the "more random" leading bits in u.
        // This also makes Bool() and Integer(2) equivalent.
        return u / r;
    }
}

 > 
 > -- 
 > James Theiler            Space and Remote Sensing Sciences
 > MS-B244, ISR-2, LANL     Los Alamos National Laboratory
 > Los Alamos, NM 87545     http://nis-www.lanl.gov/~jt

One final remark.  knuthran.c is out of date!  See

    http://www-cs-faculty.stanford.edu/~knuth/news02.html#rng
    http://www-cs-faculty.stanford.edu/~knuth/programs/rng.c

There are two significant changes:

(1) The generator is "warmed up" more thoroughly.  (This avoids problems
    with correlations between random streams with adjacent seeds.)

(2) Only 100 out of every 1009 random numbers are used.  (This is needed
    so that the "birthday test" is satisfied.)

-- 
Charles Karney <ckarney@sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

URL: http://charles.karney.info
Tel: +1 609 734 2312
Fax: +1 609 734 2323