/* rannor.c: random # ~ N(0,1)  */
/***********************************************************/
#include <math.h>
#include "ran.h"

double RanNor(void)
  /* return N(0,1) by the polar method, Knuth v2c3, page 117 */
{
  static int get_new_pair=1;
  static double x2;
  double a, s, v1, v2;

  if  (get_new_pair) {       /* generate a new pair (x1, x2) */
    do {
      v1=2.0*RanUni()-1.0;
      v2=2.0*RanUni()-1.0;
      s = v1*v1 + v2*v2;
    } while (s >= 1.0);
    a = sqrt(-2.0*log(s)/s);
    x2 = v2*a;         /* save x2 */
    get_new_pair = 0;
    return v1*a;       /* use x1 */
  }
  else {
    get_new_pair=1;
    return x2;  /* use x2 */
  }
}