# Algol68 Normal Random Deviate Generator #
MODE RNDS = STRUCT(BOOL f, REAL v);
PROC rnd = (REF RNDS s)REAL: (f OF s | f OF s := 1<0; v OF s
  | REAL x = nextrandom*2 - 1, y = nextrandom*2 - 1;
  REAL r2 = x*x + y*y;
    (r2 > 1.0 | rnd(s) | REAL q = sqrt(-2*ln(r2)/r2);
      s := (0<1, x*q); y*q));

#(RNDS rnds := (1<0, 0);
TO 20 DO print((rnd(rnds), newline)) OD)

(RNDS rnds := (1<0, 0); REAL s1 := 0, s2 := 0;
TO 1000 DO REAL d = rnd(rnds); s1 +:= d; s2 +:= d*d OD;
print((s1, s2)))#