content: -- Algorithms taken from Dodge and Jerse's Computer Music: Synthesis,
  -- Composition, and Performance, Chapter 11. 
  
  module System.Random.Distributions (
    -- * Random Distributions
    linear, exponential, bilExp, gaussian, cauchy, poisson, frequency
  
    -- * Utility Functions
    , rands
  
    ) where
  
  import System.Random
  
  {- | Given a random number generator, generates a linearly distributed
  random variable between 0 and 1.  Returns the random value together
  with a new random number generator.  The probability density function
  is given by
  
  > f(x) = 2(1-x)  0 <= x <= 1
  >      = 0       otherwise
  
  -}
  linear :: (RandomGen g, Floating a, Random a, Ord a) => g -> (a,g)
  linear g0 = 
      let (r1, g1) = randomR (0, 1) g0
          (r2, g2) = randomR (0, 1) g1
      in (min r1 r2, g2)
  
  {- | Takes a random number generator and produces another one
   that avoids generating the given number.
  -}
  avoid :: (Random a, Eq a, RandomGen g) => a -> (g -> (a,g)) -> g -> (a,g)
  avoid x f g = if r == x then avoid x f g' else (r,g')
      where (r,g') = f g
  
  {- | Generates an exponentially distributed random variable given a
  spread parameter lambda.  A larger spread increases the probability of
  generating a small number.  The mean of the distribution is
  1/lambda.  The range of the generated number is [0,inf] although
  the chance of getting a very large number is very small.
  
  The probability density function is given by
  
  > f(x) = lambda e^(-lambda * x)
  -}
  exponential :: (RandomGen g, Floating a, Random a, Eq a) => 
              a  -- ^ horizontal spread of the function.
           -> g  -- ^ a random number generator.
           -> (a,g)
  exponential lambda g0 = (-log r1 / lambda, g1)
      where (r1, g1) = avoid 0 random g0
  
  {- | Generates a random number with a bilateral exponential distribution.  
  Similar to exponential, but the mean of the distribution is 0 and
  50% of the results fall between (-1/lambda, 1/lambda).
  
  -}
  bilExp :: (Floating a, Ord a, Random a, RandomGen g) =>
            a    -- ^ horizontal spread of the function.
         -> g    -- ^ a random number generator.
         -> (a,g)
  bilExp lambda g0 = 
      let (r', g1) = avoid 0 random g0
          r = 2 * r'
          u = if r > 1 then 2 - r else r
      in (signum (1 - r) * log u / lambda, g1)
  
  {- | Generates a random number with a Gaussian distribution.  
  -}
  gaussian :: (Floating a, Random a, RandomGen g) =>
              a     -- ^ standard deviation.
           -> a     -- ^ mean.
           -> g     -- ^ a random number generator.
           -> (a,g)
  gaussian stddev center g0 = 
      let n = 12
          s = sum $ take n $ randoms g0
      in (stddev * (s - fromIntegral n / 2) + center, fst (split g0))
  
  {- | Generates a Cauchy-distributed random variable.  
  The distribution is symmetric with a mean of 0.  
  
  -}
  cauchy :: (Floating a, Random a, RandomGen g, Eq a) =>
            a   -- ^ alpha (density).
         -> g   -- ^ a random number generator.
         -> (a,g)
  cauchy density g0 = (density * tan (u * pi), g1)
      where (u, g1) = avoid 0.5 random g0
  
  {- | Generates a Poisson-distributed random variable.
  The given parameter lambda is the mean of the distribution.  
  If lambda is an integer, the probability that the result j=lambda-1
  will be as great as that of j=lambda.  The Poisson distribution
  is discrete. The returned value will be a non-negative
  integer.
  
  -}
  poisson :: (Num t, Ord a, Floating a, RandomGen g, Random a) =>
             a -> g -> (t, g)
  poisson lambda g0 = (k 0 us, g1)
      where v   = exp (-lambda)
            us  = scanl1 (*) (randoms g0)
            g1  = fst (split g0)
            k n (u:us)
                | u >= v     = k (n+1) us
                | otherwise  = n
            k _ [] = error "System.Random.Distributions.poisson: randoms did not return an infinite list"
  
  {- | Given a list of weight-value pairs, generates a value randomly picked
  from the list, weighting the probability of choosing each value by the 
  weight given.
  
  -}
  frequency :: (Floating w, Ord w, Random w, RandomGen g) 
               => [(w, a)] -> g -> (a,g)
  frequency xs g0 = (pick r xs, g1)
      where (r, g1) = randomR (0, tot) g0
            tot = sum (map fst xs)
            pick n ((w,a):xs) 
               | n <= w    = a
               | otherwise = pick (n-w) xs
            pick _ [] = error "System.Random.Distributions.frequency: The impossible happened"
  
  
  {- | Given a function generating a random number variable and a random
  number generator, produces an infinite list of random values 
  generated from the given function.
  
  -}
  rands :: (RandomGen g, Random a) => 
           (g -> (a,g)) -> g -> [a]
  rands f g = x : rands f g' where (x,g') = f g
  
  