GSL 1.6: updated patch for gauss.c

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

Here's an updated patch for gauss.c (GSL 1.6).  It includes the
information on the timing, corrects the abs/fabs bug, and reverts to
gsl_rng_uniform for the "v" coordinate.  I had earlier argued that
having v in [-0.5, 0.5) resulted in an asymmetric result.  This is
bogus.  The endpoint v=-0.5 gets rejected in the while clause, so
there's no need to add the penalty of calling gsl_rng_uniform_pos.
(Note however, that this argument does NOT extend to the polar method.
There you really do need to avoid the edges of the square to guarantee
that the result is symmetric about zero.)  Finally I use
1-gsl_rng_uniform for u, since I presume that this will be more
efficient that calling gsl_rng_uniform_pos.

--- randist/gauss.c.orig	2003-07-25 11:18:13.000000000 -0400
+++ randist/gauss.c	2005-02-03 10:33:42.000000000 -0500
@@ -28,12 +28,17 @@
  * 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 */
-
-/* Both methods pass the statistical tests; but the polar method
- * seems to be a touch faster on my home Pentium, EVEN though we
- * are only using half of the available random deviates!
+ * the Ratio method suddenly more appealing.
+ *
+ * [Added by Charles Karney] We use Leva's implementation of the Ratio
+ * method which avoids calling log() nearly all the time and makes the
+ * Ratio method faster than the Polar method (when it produces just one
+ * result per call).  Timing per call (gcc -O2 on 866MHz Pentium,
+ * average over 10^8 calls)
+ *
+ *   Polar method: 660 ns
+ *   Ratio method: 368 ns
+ *
  */
 
 /* Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */
@@ -47,8 +52,11 @@
     {
       /* 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.
+       * If gsl_rng_uniform is used negative results would be slightly
+       * more likely than positive. */
+      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;
@@ -59,28 +67,54 @@
   return sigma * y * sqrt (-2.0 * log (r2) / r2);
 }
 
-/* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130 */
-/* K+M, ACM Trans Math Software 3 (1977) 257-260. */
+/* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130.
+ * K+M, ACM Trans Math Software 3 (1977) 257-260.
+ *
+ * [Added by Charles Karney] This is an implementation of Leva'a
+ * modifications to the original K+M method; see:
+ * J. L. Leva, ACM Trans Math Software 18 (1992) 449-453 and 454-455. */
 
 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;   /* Constants from Leva */
+  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 enclosing
+       * the K+M region v^2 <= - 4 u^2 log(u). */
+
+      /* u in (0, 1] to avoid singularity at u = 0 */
+      u = 1 - gsl_rng_uniform (r);
+
+      /* v is in the asymmetric interval [-0.5, 0.5).  However v = -0.5
+       * is rejected in the last part of the while clause.  The
+       * resulting normal deviate is strictly symmetric about 0
+       * (provided that v is symmetric once v = -0.5 is excluded). */
+      v = gsl_rng_uniform (r) - 0.5;
+
+      /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
+       * much (for efficiency). */
+      v *= 1.7156;
+
+      /* Compute Leva's quadratic form Q */
+      x = u - s;
+      y = fabs(v) - t;
+      Q = x*x + y * (a*y - b*x);
     }
-  while (x * x > -4.0 * log (u));
+  while ( Q >= r1 &&            /* accept P if Q < r1 (Leva) */
+          ( Q > r2 ||           /* reject P if Q > r2 (Leva) */
+            /* Accept if v^2 <= -4 u^2 log(u) (K+M) */
+            /* 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