#include <stdlib.h>
#include <math.h>
double rand1(void) { return(drand48()*2.-1.); }
double randg(void)      /* after Box-Muller method */
{   static int k=0;
    static double u, v, r, f;
    while (k == 0) {
        u = rand1(); v = rand1(); r = u*u+v*v;
        if ((r<1.) && (r>0.)) { f = sqrt(-2.*log(r)/r); k = 1; return(u*f); }
    }
    k = 0; return(v*f);
}
float rand1_(void) { return((float)rand1()); }
float randg_(void) { return((float)randg()); }