GSL 1.6: patch for gauss.c

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

OK, here's a patch to gauss.c (v 1.6) to implement Leva's bounds on the
ratio method for normal deviates.  It also uses gsl_rng_uniform_pos for
both Gaussian methods to ensure that the resulting normal deviates are
symmetric.

I checked Leva's paper on the timing...

    Polar:                333
    Ratio:                411
    Ratio+Knuth's bounds: 197
    Ratio+Leva's bounds:  166
    Polar (GSL):          666

These are times (us) per normal deviate on a 1992 era computer and I've
extrapolated to Polar (GSL) which only uses one of the pair of normal
deviates that the polar method produces.

Because modern computers implement log pretty efficiently, the benefits
of using the bounds are not as great as shown here.  But I think that
you will find that the ratio method should comfortably beat the polar
method

--- randist/gauss.c.orig	2003-07-25 11:18:13.000000000 -0400
+++ randist/gauss.c	2005-01-26 17:25:26.000000000 -0500
@@ -28,8 +28,8 @@
  * deviates.  We don't produce two, because then we'd have to save one
  * in a static variable for the next call, and that would screws up
  * re-entrant or threaded code, so we only produce one.  This makes
- * the Ratio method suddenly more appealing.  There are further tests
- * one can make if the log() is slow.  See Knuth for details */
+ * the Ratio method suddenly more appealing.  We include Leva's tests
+ * in the Ratio method to avoid computing the log() very often. */
 
 /* Both methods pass the statistical tests; but the polar method
  * seems to be a touch faster on my home Pentium, EVEN though we
@@ -47,8 +47,9 @@
     {
       /* choose x,y in uniform square (-1,-1) to (+1,+1) */
 
-      x = -1 + 2 * gsl_rng_uniform (r);
-      y = -1 + 2 * gsl_rng_uniform (r);
+      /* Use gsl_rng_uniform_pos so that square is centered on origin. */
+      x = -1 + 2 * gsl_rng_uniform_pos (r);
+      y = -1 + 2 * gsl_rng_uniform_pos (r);
 
       /* see if it is in the unit circle */
       r2 = x * x + y * y;
@@ -61,26 +62,38 @@
 
 /* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130 */
 /* K+M, ACM Trans Math Software 3 (1977) 257-260. */
+/* Leva, ACM Trans Math Software 18 (1992) 449-453 and 454-455. */
+/* This is a straighforward conversion of Leva's Fortran code. */
 
 double
 gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma)
 {
-  double u, v, x;
+  double u, v, x, y, Q;
+  const double s =  0.449871;
+  const double t = -0.386595;
+  const double a =  0.19600;
+  const double b =  0.25472;
+  const double r1 = 0.27597;
+  const double r2 = 0.27846;
 
-  do
+  do                 /* This loop is executed 1.369 times on average  */
     {
-      v = gsl_rng_uniform (r);
-      do
-        {
-          u = gsl_rng_uniform (r);
-        }
-      while (u == 0);
-      /* Const 1.715... = sqrt(8/e) */
-      x = 1.71552776992141359295 * (v - 0.5) / u;
+      /* Generate a point P = (u, v) uniform in a rectangle */
+      u = gsl_rng_uniform_pos (r); /* NB avoid u = 0 */
+      /* Const 1.7156 > sqrt(8/e) -- but not by too much*/
+      v = 1.7156 * (gsl_rng_uniform_pos (r) - 0.5); /* v is centered */
+      x = u - s;
+      y = abs(v) - t;
+      /* Compute quadratic form Q */
+      Q = x*x + y * (a*y - b*x);
     }
-  while (x * x > -4.0 * log (u));
+  while ( Q >= r1 &&            /* accept P if Q < r1 */
+          ( Q > r2 ||           /* reject P if Q > r2 */
+            /* Accept if v^2 <= -4 u^2 log(u) */
+            /* This final test is executed 0.012 times on average. */
+            v*v > -4 * u*u * log(u) ) );
 
-  return sigma * x;
+  return sigma * v/u;           /* Return slope */
 }
 
 double

-- 
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