# This function computes random deviates from a standard normal
# distribution (zero mean, unit variance).  Called with no arguments,
# it returns a scalar.  Otherwise, the argument specifies the shape
# of the vector or matrix returned; see the `readnum' function for
# more information.

# The same seed is used in both `rand' and `randn'; it may be
# specified with the `srand' function.  If you don't call `srand',
# the seed is based on the system's clock.

# This implementation uses Knuth's polar algorithm.

randn = function (s)
{
  local (u; v; p; n; i; m; r; j; t);

  s = vector (s);
  n = 1;
  for (i in seq (s.ne)) { n = n*s[i]; }

  r = fill (0; 0.0);

  while (m = n - r.ne)
  {
    t = 1.3 * m + 5;
    u = 2.0 * rand (t) - 1;
    v = 2.0 * rand (t) - 1;
    p = u^2 + v^2;
    j = find (1; p < 1.0);
    if (j.ne > m) { j = j[1:m]; }
    p = p[j];
    u = u[j];
    r = r, sqrt (-2.0*log(p) / p) @ u;
  }

  return fill (s; r);
};