🚀 go-pugleaf

paul@music.sie.arizona.edu (Paul J. Sanchez) writes:

>There are a couple of morals to this story:
>a) If you want to use Box-Muller, use something other than a LCG; and

Or, if you want to use a LCG then don't use the Box-Muller?

Could someone comment on the algorithm in Algorithm 488, CACM, Vol.
17, #12, p704, which draws from N(0,1).  

C source for it is here:

#define REAL float
/* Algorithm 488, CACM, Vol. 17, #12, p704.

Except on the first call ranstdn2 returns a pseudo-random number
having a gaussian (i.e.  normal) distribution with zero mean and unit
standard deviation.
Thus, the density is f(x) = exp(-0.5*x**2)/sqrt(2.0*pi).
The first call initializes ranstdn2 and returns zero
(no hassles -- it's a legit draw from N(0,1)).

It calls a function ranpm(), and it is assumed that successive calls
to ranpm(0) give independent pseudo- random numbers distributed
uniformly on (0, 1), possibly including 0 (but not 1).  the method
used was suggested by Von Neumann, and improved by Forsythe, Ahrens,
Dieter and Brent.  on the average there are 1.37746 calls of ranpm() for
each call of ranstdn().

Warning - dimension and data statements below are machine-dependent.
The dimension of d must be at least the number of bits in the fraction
of a REAL number.  Thus, on most machines the data statement
below can be truncated.

How is d calculated:
If the integral of sqrt(2.0/pi)*exp(-0.5*x**2) from a(i) to infinity
is 2**(-i), then d(i) = a(i) - a(i-1).   Thus it seems that we could
do better than the constants in this code by calculating d[] in
double precision using a modern pnorm1.

On a ELC using gcc 1.4 -O, it gives 2400 draws a CPU second.
*/

REAL ranstdn2(void)
{
    /* Initialized data */
    static REAL d[60] = { .67448975f,.47585963f,.383771164f,.328611323f,
	    .291142827f,.263684322f,.242508452f,.225567444f,.211634166f,
	    .199924267f,.189910758f,.181225181f,.1736014f,.166841909f,
	    .160796729f,.155349717f,.150409384f,.145902577f,.141770033f,
	    .137963174f,.134441762f,.13117215f,.128125965f,.12527909f,
	    .122610883f,.12010356f,.117741707f,.115511892f,.113402349f,
	    .11140272f,.109503852f,.107697617f,.105976772f,.104334841f,
	    .102766012f,.101265052f,.099827234f,.098448282f,.097124309f,
	    .095851778f,.094627461f,.093448407f,.092311909f,.091215482f,
	    .090156838f,.089133867f,.088144619f,.087187293f,.086260215f,
	    .085361834f,.084490706f,.083645487f,.082824924f,.082027847f,
	    .081253162f,.080499844f,.079766932f,.079053527f,.078358781f,
	    .077681899f };
    static REAL u = 0;

    /* Local variables */
    static REAL a;
    static int i;
    static REAL v, w;

/* U MUST BE PRESERVED BETWEEN CALLS. */
/* INITIALIZE DISPLACEMENT A AND COUNTER I. */
    a = 0;
    i = 0;
/* INCREMENT COUNTER AND DISPLACEMENT IF LEADING BIT */
/* OF U IS ONE. */
L10:
    u += u;
    if (u < 1) goto L20;
    u -= 1;
    ++i;
    a -= d[i - 1];
    goto L10;
/* FORM W UNIFORM ON 0 .LE. W .LT. D(I+1) FROM U. */
L20:
    w = d[i] * u;
/* FORM V = 0.5*((W-A)**2 - A**2). NOTE THAT 0 .LE. V .LT. LOG(2). */
    v = w * (w * .5f - a);
/* GENERATE NEW UNIFORM U. */
L30:
    u = ranpm();
/* ACCEPT W AS A RANDOM SAMPLE IF V .LE. U. */
    if (v <= u) goto L40;
/* GENERATE RANDOM V. */
    v = ranpm();
/* LOOP IF U .GT. V. */
    if (u > v) goto L30;
/* REJECT W AND FORM A NEW UNIFORM U FROM V AND U. */
    u = (v - u) / (1 - u);
    goto L20;
/* FORM NEW U (TO BE USED ON NEXT CALL) FROM U AND V. */
L40:
    u = (u - v) / (1 - v);
/* USE FIRST BIT OF U FOR SIGN, RETURN NORMAL VARIATE. */
    u += u;
    if (u < 1) return (a-w);
    else {
        u -= 1; return (w-a);
    }
} /* ranstdn2() */

#ifdef TESTING
#include <stdlib.h>
int main(int argc, char **argv)
{
    int N, i; 
    if (argc != 3) {
        fprintf(stderr, "Usage: %s [12] N\n", argv[0]);
        return 1;
    }
    N = atoi(argv[2]);
    if (atoi(argv[1]) == 2) {
        (void) ranstdn2(); /* swallow the 0 */
        for (i=0; i<N; i++) printf("%.6g\n", ranstdn2());
    } else
      for (i=0; i<N; i++) printf("%.6g\n", ranstdn());
    return 0;
}
#endif /* TESTING */

-- 
Ajay Shah, (213)749-8133, ajayshah@usc.edu