randdist.nim

import random, math

template gaussImpl(randCall) =
  const K = sqrt(2 / E)
  var
    a = 0.0
    b = 0.0
  while true:
    a = randCall
    b = (2.0 * randCall - 1.0) * K
    if  b * b <= -4.0 * a * a * ln(a): break
  result = mu + sigma * (b / a)

proc gauss*(r: var Rand; mu = 0.0; sigma = 1.0): float =
  ## Returns a Gaussian random variate,
  ## with mean ``mu`` and standard deviation ``sigma``
  ## using the given state.
  # Ratio of uniforms method for normal
  # http://www2.econ.osaka-u.ac.jp/~tanizaki/class/2013/econome3/13.pdf
  gaussImpl(rand(r, 1.0))

proc gauss*(mu = 0.0, sigma = 1.0): float =
  ## Returns a Gaussian random variate,
  ## with mean ``mu`` and standard deviation ``sigma``.
  ##
  ## If `randomize<#randomize>`_ has not been called, the order of outcomes
  ## from this proc will always be the same.
  ##
  ## This proc uses the default random number generator. Thus, it is **not**
  ## thread-safe.
  gaussImpl(rand(1.0))

echo gauss()