#include "misc.h"
#include "rnd.h"
#include <stdlib.h>
#include <math.h>
// Each invocation of rnd returns a new random normal deviate.
static R scale = 0x1.p-30;
R rnd(){static int f = 0;
  static R m;
  if(f) {f = 0; return m;}
  while (1){R x = scale*random()-1, y = scale*random()-1,
    s = x*x +y*y; 
    if(s<1) {R r = sqrt(- 2*log(s)/s);
    f = 1; m = x*r; return y*r;}}}