#include <math.h>
#define freq(x) (0.5*erfc(-0.707106781186547524*(x)))
#define rdfreq(x) (exp(0.5*(x)*(x))*2.50662827463100050242)

#define C0 2.515517
#define C1 0.802853
#define C2 0.010328
#define D0 1.0
#define D1 1.432788
#define D2 0.189269
#define D3 0.001308

/* Joel Heinrich    23 March 2005

   returns Gaussian deviate with mean 0, variance 1

   The user must supply a uniform random generator: gauss_gen's 
   argument is a pointer to a function that takes no arguments 
   and returns a double, where the double is a uniform random deviate.

   Uses the transformation method, inverting the integral of the Gaussian.
   Approximate inverse taken from Abramowitz and Stegun 26.2.23


*/


double gauss_gen( double (*ru)() ) {
  const double y = ru(), yy = (y>0.5) ? 1-y : y, t = sqrt(-2*log(yy));
  double x = (C0+t*(C1+t*C2))/(D0+t*(D1+t*(D2+t*D3))) - t; /* approximate */
  x -= (freq(x) - yy)*rdfreq(x);  /* two Newton iterations to refine */
  x -= (freq(x) - yy)*rdfreq(x);
  return (y>0.5) ? -x : x;
}