From 46514f3f6b4711e4909640635f70fa67d3dddcd8 Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Sun, 3 Feb 2008 16:38:17 +0000
Subject: [PATCH] Eliminate specials of rng instances *rng-mt19937*,
 *rng-cmrg*, *rng-default*, in favor of the more approprate letm-bound
 objects.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3281 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
---
 linear-least-squares.lisp      |  16 +++---
 monte-carlo.lisp               |  53 ++++++-----------
 random/bernoulli.lisp          |  30 +++++-----
 random/beta.lisp               |  32 ++++++-----
 random/binomial.lisp           |  46 +++++++--------
 random/cauchy.lisp             |  49 ++++++++--------
 random/chi-squared.lisp        |  49 ++++++++--------
 random/dirichlet.lisp          |  30 +++++-----
 random/discrete.lisp           |  42 ++++++--------
 random/exponential-power.lisp  |  48 +++++++---------
 random/exponential.lisp        |  50 ++++++++--------
 random/fdist.lisp              |  59 ++++++++++---------
 random/flat.lisp               |  49 ++++++++--------
 random/gamma.lisp              |  54 +++++++++--------
 random/gaussian-bivariate.lisp |  37 ++++++------
 random/gaussian-tail.lisp      |  42 ++++++--------
 random/gaussian.lisp           | 102 ++++++++++++++++-----------------
 random/generators.lisp         |   6 +-
 random/geometric.lisp          |  47 +++++++--------
 random/gumbel1.lisp            |  50 ++++++++--------
 random/gumbel2.lisp            |  59 +++++++++----------
 random/hypergeometric.lisp     |  57 ++++++++----------
 random/landau.lisp             |  24 +++-----
 random/laplace.lisp            |  49 ++++++++--------
 random/levy.lisp               |  69 ++++++++++------------
 random/logarithmic.lisp        |  27 ++++-----
 random/logistic.lisp           |  47 ++++++++-------
 random/lognormal.lisp          |  54 +++++++++--------
 random/multinomial.lisp        |  76 +++++++++++-------------
 random/negative-binomial.lisp  |  84 +++++++++++++--------------
 random/pareto.lisp             |  49 ++++++++--------
 random/poisson.lisp            |  24 ++++----
 random/rayleigh-tail.lisp      |  18 +++---
 random/rayleigh.lisp           |  57 +++++++++---------
 random/spherical-vector.lisp   |  53 ++++++++---------
 random/tdist.lisp              |  53 +++++++++--------
 random/weibull.lisp            |  47 ++++++++-------
 37 files changed, 809 insertions(+), 929 deletions(-)

diff --git a/linear-least-squares.lisp b/linear-least-squares.lisp
index 2eca08ce..ae3867b8 100644
--- a/linear-least-squares.lisp
+++ b/linear-least-squares.lisp
@@ -1,6 +1,6 @@
 ;; Linear least squares, or linear regression
 ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp>
-;; Time-stamp: <2008-02-02 21:01:35EST linear-least-squares.lisp>
+;; Time-stamp: <2008-02-03 11:31:47EST linear-least-squares.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -290,13 +290,13 @@
 (defun mv-linear-least-squares-data ()
   "Generate data for second example in Section 36.5 of the GSL
    manual."
-  (rng-set *rng-mt19937* 0)
-  (loop for x from 1/10 below 2 by 1/10
-	for xd = (coerce x 'double-float)
-	for y0 = (exp xd)
-	for sigma = (* 0.1d0 y0)
-	collect
-	(list xd (+ y0 (gaussian *rng-mt19937* sigma)) sigma)))
+  (letm ((rng (random-number-generator *mt19937* 0)))
+    (loop for x from 1/10 below 2 by 1/10
+	  for xd = (coerce x 'double-float)
+	  for y0 = (exp xd)
+	  for sigma = (* 0.1d0 y0)
+	  collect
+	  (list xd (+ y0 (gaussian rng sigma)) sigma))))
 
 (defun mv-linear-least-squares-example (data)
   "Second example in Section 36.5 of the GSL manual."
diff --git a/monte-carlo.lisp b/monte-carlo.lisp
index 758304e3..768ce204 100644
--- a/monte-carlo.lisp
+++ b/monte-carlo.lisp
@@ -1,6 +1,6 @@
 ;; Monte Carlo Integration
 ;; Liam Healy Sat Feb  3 2007 - 17:42
-;; Time-stamp: <2008-02-02 21:11:35EST monte-carlo.lisp>
+;; Time-stamp: <2008-02-03 11:35:59EST monte-carlo.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -288,46 +288,25 @@
 (def-mc-function monte-carlo-g 3)
 
 (defun random-walk-plain-example (&optional (nsamples 500000))
-  (letm ((ws (monte-carlo-plain 3)))
-    (with-data (lower vector-double 3)
-      (with-data (upper vector-double 3)
-	(setf (data lower) #(0.0d0 0.0d0 0.0d0)
-	      (data upper) (vector pi pi pi))
-	(rng-set *rng-mt19937* 0)
-	(monte-carlo-integrate-plain
-	 monte-carlo-g
-	 lower upper
-	 nsamples
-	 *rng-mt19937*
-	 ws)))))
+  (letm ((ws (monte-carlo-plain 3))
+	 (lower (vector-double #(0.0d0 0.0d0 0.0d0)))
+	 (upper (vector-double (vector pi pi pi)))
+	 (rng (random-number-generator *mt19937* 0)))
+    (monte-carlo-integrate-plain monte-carlo-g lower upper nsamples rng ws)))))
 
 (defun random-walk-miser-example (&optional (nsamples 500000))
-  (letm ((ws (monte-carlo-miser 3)))
-    (with-data (lower vector-double 3)
-      (with-data (upper vector-double 3)
-	(setf (data lower) #(0.0d0 0.0d0 0.0d0)
-	      (data upper) (vector pi pi pi))
-	(rng-set *rng-mt19937* 0)
-	(monte-carlo-integrate-miser
-	 monte-carlo-g
-	 lower upper
-	 nsamples
-	 *rng-mt19937*
-	 ws)))))
+  (letm ((ws (monte-carlo-miser 3))
+	 (lower (vector-double #(0.0d0 0.0d0 0.0d0)))
+	 (upper (vector-double (vector pi pi pi)))
+	 (rng (random-number-generator *mt19937* 0)))
+    (monte-carlo-integrate-miser monte-carlo-g lower upper nsamples rng ws)))))
 
 (defun random-walk-vegas-example (&optional (nsamples 500000))
-  (letm ((ws (monte-carlo-vegas 3)))
-    (with-data (lower vector-double 3)
-      (with-data (upper vector-double 3)
-	(setf (data lower) #(0.0d0 0.0d0 0.0d0)
-	      (data upper) (vector pi pi pi))
-	(rng-set *rng-mt19937* 0)
-	(monte-carlo-integrate-vegas
-	 monte-carlo-g
-	 lower upper
-	 nsamples
-	 *rng-mt19937*
-	 ws)))))
+  (letm ((ws (monte-carlo-vegas 3))
+	 (lower (vector-double #(0.0d0 0.0d0 0.0d0)))
+	 (upper (vector-double (vector pi pi pi)))
+	 (rng (random-number-generator *mt19937* 0)))
+    (monte-carlo-integrate-vegas monte-carlo-g lower upper nsamples rng ws)))))
 
 (lisp-unit:define-test monte-carlo
   (lisp-unit:assert-first-fp-equal
diff --git a/random/bernoulli.lisp b/random/bernoulli.lisp
index 69a80d79..ed9e8a0f 100644
--- a/random/bernoulli.lisp
+++ b/random/bernoulli.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        bernoulli.lisp                          
-; description: Bernoulli distribution                  
-; date:        Sat Nov 25 2006 - 16:59
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 17:01
-;********************************************************
-;;; $Id: $
+;; Bernoulli distribution
+;; Liam Healy, Sat Nov 25 2006 - 16:59
+;; Time-stamp: <2008-02-03 11:16:35EST bernoulli.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,9 +9,9 @@
   "gsl_ran_bernoulli"
   (((generator generator) :pointer) (p :double))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "Returns either 0 or 1, the result of a Bernoulli trial
-   with probability @var{p}.  The probability distribution for
+   with probability p.  The probability distribution for
    a Bernoulli trial is
    p(0) = 1 - p
    p(1) = p.")
@@ -23,19 +19,19 @@
 (defun-gsl bernoulli-pdf (k p)
   "gsl_ran_bernoulli_pdf" ((k :uint) (p :double))
   :c-return :double
-  :documentation "The probability @math{p(k)} of obtaining
-  @var{k} from a Bernoulli distribution with probability parameter
-  @var{p}, using the formula given in #'bernoulli.")
+  :documentation			; FDL
+  "The probability p(k) of obtaining
+  k from a Bernoulli distribution with probability parameter
+  p, using the formula given in #'bernoulli.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test bernoulli
   (lisp-unit:assert-equal
    '(0 1 1 0 1 1 0 0 0 0 0)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(bernoulli *rng-mt19937* 0.5d0))))
+	   collect
+	   (bernoulli rng 0.5d0))))
   (lisp-unit:assert-first-fp-equal
    "0.500000000000d+00"
    (bernoulli-pdf 0 0.5d0)))
diff --git a/random/beta.lisp b/random/beta.lisp
index 2e506106..d71a1747 100644
--- a/random/beta.lisp
+++ b/random/beta.lisp
@@ -1,6 +1,6 @@
 ;; Beta distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2008-02-02 23:04:13EST beta.lisp>
+;; Time-stamp: <2008-02-03 09:52:01EST beta.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -10,42 +10,46 @@
   "gsl_ran_beta"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the beta distribution.  The distribution function is
    p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx
-   @math{0 <= x <= 1}.")
+   0 <= x <= 1.")
 
 (defun-gsl beta-pdf (x a b)
   "gsl_ran_beta_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a beta distribution with parameters @var{a} and @var{b}, using the
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a beta distribution with parameters a and b, using the
    formula given in #'beta.")
 
 (defun-gsl beta-P (x a b)
   "gsl_cdf_beta_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the beta distribution with parameters a and b.")
 
 (defun-gsl beta-Q (x a b)
   "gsl_cdf_beta_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the beta distribution with parameters a and b.")
 
 (defun-gsl beta-Pinv (P a b)
   "gsl_cdf_beta_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the beta distribution with parameters a and b.")
 
 (defun-gsl beta-Qinv (Q a b)
   "gsl_cdf_beta_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the beta distribution with parameters a and b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test beta
diff --git a/random/binomial.lisp b/random/binomial.lisp
index d0a03177..f214d026 100644
--- a/random/binomial.lisp
+++ b/random/binomial.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        binomial.lisp                          
-; description: Binomial distribution
-; date:        Sat Nov 25 2006 - 16:00
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 18:08
-;********************************************************
-;;; $Id: $
+;; Binomial distribution
+;; Liam Healy, Sat Nov 25 2006 - 16:00
+;; Time-stamp: <2008-02-03 11:25:08EST binomial.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,44 +9,44 @@
   "gsl_ran_binomial"
   (((generator generator) :pointer) (p :double) (n :uint))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "A random integer from the binomial distribution,
-  the number of successes in @var{n} independent trials with probability
-  @var{p}.  The probability distribution for binomial variates is,
+  the number of successes in n independent trials with probability
+  p.  The probability distribution for binomial variates is,
   p(k) = {n! \over k! (n-k)!} p^k (1-p)^{n-k}
-  @math{0 <= k <= n}.")
+  0 <= k <= n.")
 
 (defun-gsl binomial-pdf (k p n)
   "gsl_ran_binomial_pdf" ((k :uint) (p :double) (n :uint))
   :c-return :double
-  :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a binomial distribution with parameters @var{p} and @var{n}, using
+  :documentation			; FDL
+  "The probability p(k) of obtaining k
+   from a binomial distribution with parameters p and n, using
    the formula given in #'binomial.")
 
 (defun-gsl binomial-P (k p n)
   "gsl_cdf_binomial_P" ((k :uint) (p :double) (n :uint))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(k)} for the Binomial distribution
-  with parameters @var{p} and @var{n}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(k) for the Binomial distribution with parameters p and n.")
 
 (defun-gsl binomial-Q (k p n)
   "gsl_cdf_binomial_Q" ((k :uint) (p :double) (n :uint))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(k)} for the Binomial distribution
-   with parameters @var{p} and @var{n}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+   Q(k) for the Binomial distribution
+   with parameters p and n.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test binomial
   (lisp-unit:assert-equal
    '(11 3 4 8 4 5 8 6 5 6 6)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(binomial *rng-mt19937* 0.4d0 12))))
+	   collect
+	   (binomial rng 0.4d0 12))))
   (lisp-unit:assert-first-fp-equal
    "0.227030335488d+00"
    (binomial-pdf 5 0.4d0 12))
diff --git a/random/cauchy.lisp b/random/cauchy.lisp
index bb64c1bf..e05cdfbc 100644
--- a/random/cauchy.lisp
+++ b/random/cauchy.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        cauchy.lisp                          
-; description: Cauchy distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Sep 30 2006 - 19:15
-;********************************************************
-;;; $Id: $
+;; Cauchy distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 11:02:19EST cauchy.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,45 +9,49 @@
   "gsl_ran_cauchy"
   (((generator generator) :pointer) (a :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Cauchy distribution with
-   scale parameter @var{a}.  The probability distribution for Cauchy
+   scale parameter a.  The probability distribution for Cauchy
    random variates is,
    p(x) dx = {1 \over a\pi (1 + (x/a)^2) } dx
-   for @math{x} in the range @math{-\infty} to @math{+\infty}.  The Cauchy
+   for x in the range -\infty to +\infty.  The Cauchy
    distribution is also known as the Lorentz distribution.")
 
 (defun-gsl cauchy-pdf (x a)
   "gsl_ran_cauchy_pdf" ((x :double) (a :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Cauchy distribution with scale parameter @var{a}, using the formula
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Cauchy distribution with scale parameter a, using the formula
    given for #'cauchy.")
 
 (defun-gsl cauchy-P (x a)
   "gsl_cdf_cauchy_P" ((x :double) (a :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the Cauchy distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the Cauchy distribution with scale parameter a.")
 
 (defun-gsl cauchy-Q (x a)
   "gsl_cdf_cauchy_Q" ((x :double) (a :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the Cauchy distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the Cauchy distribution with scale parameter a.")
 
 (defun-gsl cauchy-Pinv (P a)
   "gsl_cdf_cauchy_Pinv" ((P :double) (a :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the Cauchy distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the Cauchy distribution with scale parameter a.")
 
 (defun-gsl cauchy-Qinv (Q a)
   "gsl_cdf_cauchy_Qinv" ((Q :double) (a :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the Cauchy distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the Cauchy distribution with scale parameter a.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test cauchy
@@ -61,11 +61,10 @@
      "-0.134390491844d+01" "-0.103643632829d+02" "-0.790709314248d+02"
      "-0.106520710880d+02" "-0.939394824349d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (cauchy *rng-mt19937* 10.0d0)))))
+	    (cauchy rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.318309886184d-01"
    (cauchy-pdf 0.0d0 10.0d0))
diff --git a/random/chi-squared.lisp b/random/chi-squared.lisp
index 7bc51cb9..1ae331d4 100644
--- a/random/chi-squared.lisp
+++ b/random/chi-squared.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        chi-squared.lisp                          
-; description: Chi-squared distribution                  
-; date:        Sat Oct  7 2006 - 16:13
-; author:      Liam M. Healy                             
-; modified:    Sun Oct  8 2006 - 16:41
-;********************************************************
-;;; $Id: $
+;; Chi-squared distribution
+;; Liam Healy, Sat Oct  7 2006 - 16:13
+;; Time-stamp: <2008-02-03 10:23:21EST chi-squared.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,43 +9,47 @@
   "gsl_ran_chisq"
   (((generator generator) :pointer) (nu :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the chi-squared distribution
-  with @var{nu} degrees of freedom. The distribution function is
+  with nu degrees of freedom. The distribution function is
   p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx
-  @math{x >= 0}. ")
+  x >= 0. ")
 
 (defun-gsl chi-squared-pdf (x nu)
   "gsl_ran_chisq_pdf" ((x :double) (nu :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a chi-squared distribution with @var{nu} degrees of freedom, using
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a chi-squared distribution with nu degrees of freedom, using
    the formula given in #'chi-squared.")
 
 (defun-gsl chi-squared-P (x nu)
   "gsl_cdf_chisq_P" ((x :double) (nu :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the chi-squared distribution with nu degrees of freedom.")
 
 (defun-gsl chi-squared-Q (x nu)
   "gsl_cdf_chisq_Q" ((x :double) (nu :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the chi-squared distribution with nu degrees of freedom.")
 
 (defun-gsl chi-squared-Pinv (P nu)
   "gsl_cdf_chisq_Pinv" ((P :double) (nu :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the chi-squared distribution with nu degrees of freedom.")
 
 (defun-gsl chi-squared-Qinv (Q nu)
   "gsl_cdf_chisq_Qinv" ((Q :double) (nu :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+  :documentation 			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the chi-squared distribution with nu degrees of freedom.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test chi-squared
@@ -59,11 +59,10 @@
      "0.137566397216d+02" "0.113594580454d+02" "0.160682825906d+02"
      "0.129995142525d+02" "0.111473474244d+02")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (chi-squared *rng-mt19937* 10.0d0)))))
+	    (chi-squared rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.439391289468d+00"
    (chi-squared-pdf 0.5d0 1.0d0))
diff --git a/random/dirichlet.lisp b/random/dirichlet.lisp
index 618c46f2..86bec3be 100644
--- a/random/dirichlet.lisp
+++ b/random/dirichlet.lisp
@@ -1,6 +1,6 @@
 ;; Dirichlet distribution
 ;; Liam Healy, Sun Oct 29 2006
-;; Time-stamp: <2008-02-02 23:00:40EST dirichlet.lisp>
+;; Time-stamp: <2008-02-03 09:49:17EST dirichlet.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -13,20 +13,20 @@
    ;; theta had better be at least as long as alpha, or they'll be trouble
    ((gsl-array theta) :pointer))
   :c-return :void
-  :documentation
+  :documentation			; FDL
   "An array of K=(length alpha) random variates from a Dirichlet
   distribution of order K-1.  The distribution function is
   p(\theta_1,\ldots,\theta_K) \, d\theta_1 \cdots d\theta_K = 
         {1 \over Z} \prod_{i=1}^{K} \theta_i^{\alpha_i - 1} 
           \; \delta(1 -\sum_{i=1}^K \theta_i) d\theta_1 \cdots d\theta_K
-  @math{theta_i >= 0} and @math{alpha_i >= 0}.
-  The delta function ensures that @math{\sum \theta_i = 1}.
-  The normalization factor @math{Z} is
+  theta_i >= 0 and alpha_i >= 0.
+  The delta function ensures that \sum \theta_i = 1.
+  The normalization factor Z is
   Z = {\prod_{i=1}^K \Gamma(\alpha_i) \over \Gamma( \sum_{i=1}^K \alpha_i)}
-  The random variates are generated by sampling @var{K} values 
-  from gamma distributions with parameters @math{a=alpha_i, b=1}, 
+  The random variates are generated by sampling K values 
+  from gamma distributions with parameters a=alpha_i, b=1, 
   and renormalizing. 
-  See A.M. Law, W.D. Kelton, @cite{Simulation Modeling and Analysis}
+  See A.M. Law, W.D. Kelton, \"Simulation Modeling and Analysis\"
   (1991).")
 
 (defun-gsl dirichlet-pdf (alpha theta)
@@ -36,10 +36,10 @@
    ;; theta had better be at least as long as alpha, or they'll be trouble
    ((gsl-array theta) :pointer))
   :c-return :double
-  :documentation
-  "The probability density @math{p(\theta_1, ... , \theta_K)}
-   at @var{theta}[@var{K}] for a Dirichlet distribution with parameters 
-   @var{alpha}[@var{K}], using the formula given for #'dirichlet.")
+  :documentation			; FDL
+  "The probability density p(\theta_1, ... , \theta_K)
+   at theta[K] for a Dirichlet distribution with parameters 
+   alpha[K], using the formula given for #'dirichlet.")
 
 (defun-gsl dirichlet-log-pdf (alpha theta)
   "gsl_ran_dirichlet_lnpdf"
@@ -48,11 +48,11 @@
    ;; theta had better be at least as long as alpha, or they'll be trouble
    ((gsl-array theta) :pointer))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "The logarithm of the probability density 
-   @math{p(\theta_1, ... , \theta_K)}
+   p(\theta_1, ... , \theta_K)
    for a Dirichlet distribution with parameters 
-   @var{alpha}[@var{K}].")
+   alpha[K].")
 
 ;;; Examples and unit test
 (lisp-unit:define-test dirichlet
diff --git a/random/discrete.lisp b/random/discrete.lisp
index 064a6cea..763ca6ba 100644
--- a/random/discrete.lisp
+++ b/random/discrete.lisp
@@ -1,6 +1,6 @@
 ;; Discrete random variables
 ;; Liam Healy, Sat Nov 11 2006 - 21:51
-;; Time-stamp: <2008-02-02 20:59:58EST discrete.lisp>
+;; Time-stamp: <2008-02-03 10:58:16EST discrete.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -18,8 +18,8 @@
   (((dim0 probabilities) :size) ((gsl-array probabilities) :pointer))
   :c-return :pointer
   :export nil
-  :index (letm discrete)
-  :documentation
+  :index (letm discrete-random)
+  :documentation			; FDL
   "A pointer to a structure that contains the lookup
   table for the discrete random number generator.  The array probabilities contains
   the probabilities of the discrete events; these array elements must all be 
@@ -31,22 +31,15 @@
   "gsl_ran_discrete_free" ((table :pointer))
   :c-return :void
   :export nil
-  :index (letm discrete)
-  :documentation
+  :index (letm discrete-random)
+  :documentation			; FDL
   "De-allocates the lookup table created by #'discrete-preprocess.")
 
-(export 'with-discrete-table)
-(defmacro with-discrete-table ((probabilities table) &body body)
-  `(let ((,table (discrete-preprocess ,probabilities)))
-     (unwind-protect
-	  (progn ,@body)
-       (discrete-free ,table))))
-
 (defun-gsl discrete (generator table)
   "gsl_ran_discrete"
   (((generator generator) :pointer) (table :pointer))
   :c-return :size
-  :documentation
+  :documentation			; FDL
   "Generate discrete random numbers after running #'discrete-preprocess;
    the argument 'table is the value returned by #'discrete-preprocess.")
 
@@ -54,29 +47,28 @@
   "gsl_ran_discrete_pdf"
   ((k :size) (table :pointer))
   :c-return :double
-  :documentation
-  "The probability @math{P[k]} of observing the variable @var{k}.
-   Since @math{P[k]} is not stored as part of the lookup table, it must be
-   recomputed; this computation takes @math{O(K)}, so if @var{K} is large
-   and you care about the original array @math{P[k]} used to create the
-   lookup table, then you should just keep this original array @math{P[k]}
+  :documentation			; FDL
+  "The probability P[k] of observing the variable k.
+   Since P[k] is not stored as part of the lookup table, it must be
+   recomputed; this computation takes O(K), so if K is large
+   and you care about the original array P[k] used to create the
+   lookup table, then you should just keep this original array P[k]
    around.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test discrete
   (lisp-unit:assert-equal
    '(1 0 1 1 0 1 1 2 1 2 2)
-   (with-data (probabilities vector-double 3)
-     (setf (data probabilities) #(0.25d0 0.5d0 0.25d0))
+   ;; Must have two letms because the vector value is not set until
+   ;; the body, but the discrete-random needs that set value.
+   (letm ((probabilities (vector-double #(0.25d0 0.5d0 0.25d0))))
      (letm ((table (discrete-random probabilities))
 	    (rng (random-number-generator *mt19937* 0)))
        (loop for i from 0 to 10
 	     collect
-	     (discrete *rng-mt19937* table))))
+	     (discrete rng table))))
    (lisp-unit:assert-first-fp-equal
     "0.500000000000d+00"
-    (with-data (probabilities vector-double 3)
-      (setf (data probabilities)
-	    #(0.25d0 0.5d0 0.25d0))
+    (letm ((probabilities (vector-double #(0.25d0 0.5d0 0.25d0))))
       (letm ((table (discrete-random probabilities)))
 	(discrete-pdf 1 table))))))
diff --git a/random/exponential-power.lisp b/random/exponential-power.lisp
index 86001d83..4d69f708 100644
--- a/random/exponential-power.lisp
+++ b/random/exponential-power.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        exponential-power.lisp                          
-; description: Exponential power distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Sep 30 2006 - 19:03
-;********************************************************
-;;; $Id: $
+;; Exponential power distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 11:12:27EST exponential-power.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,37 +9,38 @@
   "gsl_ran_exppow"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the exponential power distribution
-   with scale parameter @var{a} and exponent @var{b}.  The distribution is
+   with scale parameter a and exponent b.  The distribution is
    p(x) dx = {1 \over 2 a \Gamma(1+1/b)} \exp(-|x/a|^b) dx
-   for @math{x >= 0}.  For @math{b = 1} this reduces to the Laplace
-   distribution.  For @math{b = 2} it has the same form as a gaussian
-   distribution, but with @c{$a = \sqrt{2} \sigma$}
-   @math{a = \sqrt@{2@} \sigma}.")
+   for x >= 0.  For b = 1 this reduces to the Laplace
+   distribution.  For b = 2 it has the same form as a gaussian
+   distribution, but with a = \sqrt{2} \sigma.")
 
 (defun-gsl exponential-power-pdf (x a b)
   "gsl_ran_exppow_pdf" 
   ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for an exponential power distribution with scale parameter @var{a}
-  and exponent @var{b}, using the formula given for #'exponential-power.")
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for an exponential power distribution with scale parameter a
+   and exponent b, using the formula given for #'exponential-power.")
 
 (defun-gsl exponential-power-P (x a b)
   "gsl_cdf_exppow_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-  @math{P(x)}, for the exponential power distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+  P(x), for the exponential power distribution with
+  parameters a and b.")
 
 (defun-gsl exponential-power-Q (x a b)
   "gsl_cdf_exppow_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions @math{Q(x)}
+  :documentation			; FDL
+  "The cumulative distribution functions Q(x)
   for the exponential power distribution with
-  parameters @var{a} and @var{b}.")
+  parameters a and b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test exponential-power
@@ -53,11 +50,10 @@
      "-0.169473362899d+01" "-0.480323610806d+00" "-0.276417363499d-01"
      "0.631839185605d+00" "-0.124788752274d-01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (exponential-power *rng-mt19937* 1.0d0 2.0d0)))))
+	    (exponential-power rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.564189583548d+00"
    (exponential-power-pdf 0.0d0 1.0d0 2.0d0))
diff --git a/random/exponential.lisp b/random/exponential.lisp
index 13b5b248..b8582069 100644
--- a/random/exponential.lisp
+++ b/random/exponential.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        exponential.lisp                          
-; description: Exponential distribution                  
-; date:        Sat Sep  2 2006 - 19:04                   
-; author:      Liam M. Healy                             
-; modified:    Mon Sep 11 2006 - 11:58
-;********************************************************
-;;; $Id: $
+;; Exponential distribution
+;; Liam Healy, Sat Sep  2 2006 - 19:04
+;; Time-stamp: <2008-02-03 10:42:55EST exponential.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,43 +9,47 @@
   "gsl_ran_exponential"
   (((generator generator) :pointer) (mu :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the exponential distribution
-   with mean @var{mu}. The distribution is
+   with mean mu. The distribution is
    p(x) dx = {1 \over \mu} \exp(-x/\mu) dx
-   @math{x >= 0}.")
+   x >= 0.")
 
 (defun-gsl exponential-pdf (x mu)
   "gsl_ran_exponential_pdf" ((x :double) (mu :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-  for an exponential distribution with mean @var{mu}, using the formula
+  :documentation			; FDL
+  "The probability density p(x) at x
+  for an exponential distribution with mean mu, using the formula
   given for exponential.")
 
 (defun-gsl exponential-P (x mu)
   "gsl_cdf_exponential_P" ((x :double) (mu :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-   @math{P(x)} for the exponential distribution with mean @var{mu}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+   P(x) for the exponential distribution with mean mu.")
 
 (defun-gsl exponential-Q (x mu)
   "gsl_cdf_exponential_Q" ((x :double) (mu :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-   @math{Q(x)} for the exponential distribution with mean @var{mu}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+   Q(x) for the exponential distribution with mean mu.")
 
 (defun-gsl exponential-Pinv (P mu)
   "gsl_cdf_exponential_Pinv" ((P :double) (mu :double))
   :c-return :double
-  :documentation "The inverse cumulative distribution function
-   @math{P(x)} for the exponential distribution with mean @var{mu}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function
+   P(x) for the exponential distribution with mean mu.")
 
 (defun-gsl exponential-Qinv (Q mu)
   "gsl_cdf_exponential_Qinv" ((Q :double) (mu :double))
   :c-return :double
-  :documentation "The inverse cumulative distribution function
-   @math{Q(x)} for the exponential distribution with mean @var{mu}.")
+  :documentation 			; FDL
+  "The inverse cumulative distribution function
+   Q(x) for the exponential distribution with mean mu.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test exponential
@@ -59,11 +59,10 @@
      "0.434536244968d+00" "0.295303920905d+01" "0.616105293907d+01"
      "0.301168633354d+01" "0.274510798194d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (exponential *rng-mt19937* 10.0d0)))))
+	    (exponential rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.100000000000d+00"
    (exponential-pdf 0.0d0 10.0d0))
@@ -79,3 +78,4 @@
   (lisp-unit:assert-first-fp-equal
    "0.100000000000d+01"
    (exponential-qinv 0.6065306597126334d0 2.0d0)))
+ 
diff --git a/random/fdist.lisp b/random/fdist.lisp
index 979b6988..a15dee02 100644
--- a/random/fdist.lisp
+++ b/random/fdist.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        fdist.lisp                          
-; description: Fdist distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sun Oct  8 2006 - 16:31
-;********************************************************
-;;; $Id: $
+;; Fdist distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 10:17:28EST fdist.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,51 +9,55 @@
   "gsl_ran_fdist"
   (((generator generator) :pointer) (nu1 :double) (nu2 :double))
   :c-return :double
-  :documentation
-  "A random variate from the F-distribution with degrees of freedom @var{nu1}
-   and @var{nu2}. The distribution function is
+  :documentation			; FDL
+  "A random variate from the F-distribution with degrees of freedom nu1
+   and nu2.  The distribution function is
    p(x) dx = 
    { \Gamma((\nu_1 + \nu_2)/2)
         \over \Gamma(\nu_1/2) \Gamma(\nu_2/2) } 
    \nu_1^{\nu_1/2} \nu_2^{\nu_2/2} 
    x^{\nu_1/2 - 1} (\nu_2 + \nu_1 x)^{-\nu_1/2 -\nu_2/2}
-   @math{x >= 0}.")
+   for x >= 0.")
 
 (defun-gsl fdist-pdf (x nu1 nu2)
   "gsl_ran_fdist_pdf" ((x :double) (nu1 :double) (nu2 :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for an F-distribution with @var{nu1} and @var{nu2} degrees of freedom,
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for an F-distribution with nu1 and nu2 degrees of freedom,
    using the formula given #'fdist.")
 
 (defun-gsl fdist-P (x nu1 nu2)
   "gsl_cdf_fdist_P" ((x :double) (nu1 :double) (nu2 :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the fdist distribution with
-  @var{nu1} and @var{nu2} degrees of freedom.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the fdist distribution with
+  nu1 and nu2 degrees of freedom.")
 
 (defun-gsl fdist-Q (x nu1 nu2)
   "gsl_cdf_fdist_Q" ((x :double) (nu1 :double) (nu2 :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the fdist distribution with
-  @var{nu1} and @var{nu2} degrees of freedom.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the fdist distribution with
+  nu1 and nu2 degrees of freedom.")
 
 (defun-gsl fdist-Pinv (P nu1 nu2)
   "gsl_cdf_fdist_Pinv" ((P :double) (nu1 :double) (nu2 :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the fdist distribution with
-  @var{nu1} and @var{nu2} degrees of freedom.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the fdist distribution with
+  nu1 and nu2 degrees of freedom.")
 
 (defun-gsl fdist-Qinv (Q nu1 nu2)
   "gsl_cdf_fdist_Qinv" ((Q :double) (nu1 :double) (nu2 :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the fdist distribution with
-   @var{nu1} and @var{nu2} degrees of freedom.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the fdist distribution with
+   nu1 and nu2 degrees of freedom.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test fdist
@@ -67,11 +67,10 @@
      "0.635729092565d-01" "0.477966365217d+00" "0.347211676079d-01"
      "0.486974823041d+00" "0.253179451696d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (fdist *rng-mt19937* 1.0d0 2.0d0)))))
+	    (fdist rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.159471988462d+00"
    (fdist-pdf 1.2d0 1.0d0 2.0d0))
diff --git a/random/flat.lisp b/random/flat.lisp
index d19686ad..4cc16d5b 100644
--- a/random/flat.lisp
+++ b/random/flat.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        flat.lisp                          
-; description: Flat distribution                  
-; date:        Oct  7 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Oct  7 2006 - 16:13
-;********************************************************
-;;; $Id: $
+;; Flat distribution
+;; Liam Healy, Oct  7 2006
+;; Time-stamp: <2008-02-03 11:18:30EST flat.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,43 +9,47 @@
   "gsl_ran_flat"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the flat (uniform)
-   distribution from @var{a} to @var{b}. The distribution is,
+   distribution from a to b.  The distribution is
    p(x) dx = {1 \over (b-a)} dx
-   if @c{$a \le x < b$} @math{a <= x < b} and 0 otherwise.")
+   if a <= x < b, and 0 otherwise.")
 
 (defun-gsl flat-pdf (x a b)
   "gsl_ran_flat_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a uniform distribution from @var{a} to @var{b}, using the formula
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a uniform distribution from a to b, using the formula
    given for #'flat.")
 
 (defun-gsl flat-P (x a b)
   "gsl_cdf_flat_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for a uniform distribution from @var{a} to @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+   P(x) for a uniform distribution from a to b.")
 
 (defun-gsl flat-Q (x a b)
   "gsl_cdf_flat_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for a uniform distribution from @var{a} to @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+   Q(x) for a uniform distribution from a to b.")
 
 (defun-gsl flat-Pinv (P a b)
   "gsl_cdf_flat_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for a uniform distribution from @var{a} to @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   P(x) for a uniform distribution from a to b.")
 
 (defun-gsl flat-Qinv (Q a b)
   "gsl_cdf_flat_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for a uniform distribution from @var{a} to @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for a uniform distribution from a to b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test flat
@@ -59,11 +59,10 @@
      "0.195747695654d+01" "0.174430534313d+01" "0.154004365834d+01"
      "0.173995298147d+01" "0.175994379818d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (flat *rng-mt19937* 1.0d0 2.0d0)))))
+	    (flat rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.100000000000d+01"
    (flat-pdf 1.2d0 1.0d0 2.0d0))
diff --git a/random/gamma.lisp b/random/gamma.lisp
index a439c6a6..0df209fe 100644
--- a/random/gamma.lisp
+++ b/random/gamma.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gamma.lisp                          
-; description: Gamma distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Oct  7 2006 - 11:53
-;********************************************************
-;;; $Id: $
+;; Gamma distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 10:43:54EST gamma.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -14,11 +10,11 @@
   "gsl_ran_gamma"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the gamma distribution.
    The distribution function is
    p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx
-   for @math{x > 0}. The gamma distribution with an integer parameter @var{a}
+   for x > 0. The gamma distribution with an integer parameter a
    is known as the Erlang distribution.  The variates are computed using
    the algorithms from Knuth (vol 2).")
 
@@ -26,40 +22,44 @@
   "gsl_ran_gamma_mt"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A gamma variate using the Marsaglia-Tsang fast gamma method.")
 
 (defun-gsl gamma-pdf (x a b)
   "gsl_ran_gamma_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a gamma distribution with parameters @var{a} and @var{b}, using the
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a gamma distribution with parameters a and b, using the
    formula given in #'gamma.")
 
 (defun-gsl gamma-P (x a b)
   "gsl_cdf_gamma_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the Gamma distribution with parameters a and b.")
 
 (defun-gsl gamma-Q (x a b)
   "gsl_cdf_gamma_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the Gamma distribution with parameters a and b.")
 
 (defun-gsl gamma-Pinv (P a b)
   "gsl_cdf_gamma_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the Gamma distribution with parameters a and b.")
 
 (defun-gsl gamma-Qinv (Q a b)
   "gsl_cdf_gamma_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the Gamma distribution with parameters a and b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test gamma-randist
@@ -69,22 +69,20 @@
      "0.869072489937d-01" "0.590607841809d+00" "0.123221058781d+01"
      "0.602337266708d+00" "0.549021596387d+00")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (gamma-rd *rng-mt19937* 1.0d0 2.0d0)))))
+	    (gamma-rd rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-equal
    '("0.260013787613d+01" "0.252226669542d+01" "0.773142209213d+01"
      "0.422727992649d+01" "0.951930434749d-01" "0.571092010687d+00"
      "0.891063771946d+00" "0.826322120255d+00" "0.318306657206d+01"
      "0.380840036132d-02" "0.103201173341d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (gamma-mt *rng-mt19937* 1.0d0 2.0d0)))))
+	    (gamma-mt rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.475614712250d+00"
    (gamma-pdf 0.1d0 1.0d0 2.0d0))
diff --git a/random/gaussian-bivariate.lisp b/random/gaussian-bivariate.lisp
index e4f1e1ff..432ccb5e 100644
--- a/random/gaussian-bivariate.lisp
+++ b/random/gaussian-bivariate.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gaussian-bivariate.lisp                   
-; description: Gaussian bivariate distribution           
-; date:        Sat Sep  2 2006 - 16:32                   
-; author:      Liam M. Healy                             
-; modified:    Sat Sep  2 2006 - 21:34
-;********************************************************
-;;; $Id: $
+;; Gaussian bivariate distribution
+;; Liam Healy, Sat Sep  2 2006 - 16:32
+;; Time-stamp: <2008-02-03 10:30:48EST gaussian-bivariate.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -14,27 +10,27 @@
   (((generator generator) :pointer) (sigma-x :double) (sigma-y :double) (rho :double)
    (x :double) (y :double))
   :c-return :void
-  :documentation
+  :documentation			; FDL
   "Generate a pair of correlated Gaussian variates, with
-   mean zero, correlation coefficient @var{rho} and standard deviations
-   @var{sigma_x} and @var{sigma_y} in the @math{x} and @math{y} directions.
+   mean zero, correlation coefficient rho and standard deviations
+   sigma_x and sigma_y in the x and y directions.
    The probability distribution for bivariate Gaussian random variates is,
    p(x,y) dx dy
    = {1 \over 2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}}
    \exp \left(-{(x^2/\sigma_x^2 + y^2/\sigma_y^2 - 2 \rho x y/(\sigma_x\sigma_y))
    \over 2(1-\rho^2)}\right) dx dy
-   for @math{x,y} in the range @math{-\infty} to @math{+\infty}.  The
-   correlation coefficient @var{rho} should lie between @math{1} and @math{-1}.")
+   for x,y in the range -\infty to +\infty.  The
+   correlation coefficient rho should lie between 1 and -1.")
 
 (defun-gsl bivariate-gaussian-pdf (x y sigma-x sigma-y rho)
   "gsl_ran_bivariate_gaussian_pdf"
   ((x :double) (y :double) (sigma-x :double) (sigma-y :double) (rho :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x,y)} at
-   (@var{x},@var{y}) for a bivariate Gaussian distribution with standard
-   deviations @var{sigma_x}, @var{sigma_y} and correlation coefficient
-   @var{rho}, using the formula given for bivariate-gaussian.")
+  :documentation			; FDL
+  "The probability density p(x,y) at
+   (x,y) for a bivariate Gaussian distribution with standard
+   deviations sigma_x, sigma_y and correlation coefficient
+   rho, using the formula given for bivariate-gaussian.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test gaussian-bivariate
@@ -44,11 +40,10 @@
      "-0.544622941217d+00" "-0.659202684161d+00" "-0.110295166108d+00"
      "0.179318404121d+00" "0.210251049803d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (bivariate-gaussian *rng-mt19937* 1.0d0 0.75d0 0.25d0)))))
+	    (bivariate-gaussian rng 1.0d0 0.75d0 0.25d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.554826555797d+00"
    (bivariate-gaussian-pdf 0.25d0 0.5d0 0.25d0
diff --git a/random/gaussian-tail.lisp b/random/gaussian-tail.lisp
index 583d1ff6..3e1715ae 100644
--- a/random/gaussian-tail.lisp
+++ b/random/gaussian-tail.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gaussian-tail.lisp                        
-; description: Gaussian tail distribution                
-; date:        Mon Aug 21 2006 - 21:52                   
-; author:      Liam M. Healy                             
-; modified:    Mon Aug 21 2006 - 22:20
-;********************************************************
-;;; $Id: $
+;; Gaussian tail distribution
+;; Liam Healy, Mon Aug 21 2006 - 21:52
+;; Time-stamp: <2008-02-03 10:15:01EST gaussian-tail.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,38 +9,38 @@
   "gsl_ran_gaussian_tail"
   (((generator generator) :pointer) (a :double) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Random variates from the upper tail of a Gaussian
-   distribution with standard deviation @var{sigma}.  The values returned
-   are larger than the lower limit @var{a}, which must be positive.  The
+   distribution with standard deviation sigma.  The values returned
+   are larger than the lower limit a, which must be positive.  The
    method is based on Marsaglia's famous rectangle-wedge-tail algorithm (Ann. 
    Math. Stat. 32, 894--899 (1961)), with this aspect explained in Knuth, v2,
    3rd ed, p139,586 (exercise 11).
    The probability distribution for Gaussian tail random variates is,
    p(x) dx = {1 \over N(a;\sigma) \sqrt{2 \pi \sigma^2}}
                   \exp (- x^2 / 2\sigma^2) dx
-   for @math{x > a} where @math{N(a;\sigma)} is the normalization constant,
+   for x > a where N(a;\sigma) is the normalization constant,
    N(a;\sigma) = (1/2) erfc(a / sqrt(2 sigma^2)).")
 
 (defun-gsl gaussian-tail-pdf (x a sigma)
   "gsl_ran_gaussian_tail_pdf" ((x :double) (a :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-  for a Gaussian tail distribution with standard deviation @var{sigma} and
-  lower limit @var{a}, using the formula given for gaussian-tail.")
+  :documentation			; FDL
+  "The probability density p(x) at x
+  for a Gaussian tail distribution with standard deviation sigma and
+  lower limit a, using the formula given for gaussian-tail.")
 
 (defun-gsl ugaussian-tail (generator a)
   "gsl_ran_ugaussian_tail"
   (((generator generator) :pointer) (a :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Equivalent to gaussian-tail with sigma=1.")
 
 (defun-gsl ugaussian-tail-pdf (x a)
   "gsl_ran_ugaussian_tail_pdf" ((x :double) (a :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Equivalent to gaussian-tail-pdf with sigma=1.")
 
 ;;; Examples and unit test
@@ -55,11 +51,10 @@
      "0.534422862873d+02" "0.518376171418d+02" "0.530010742143d+02"
      "0.521497741699d+02" "0.501157244350d+02")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (gaussian-tail *rng-mt19937* 50.0d0 10.0d0)))))
+	    (gaussian-tail rng 50.0d0 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.187022708773d+00"
    (gaussian-tail-pdf 52.0d0 50.0d0 10.0d0))
@@ -69,11 +64,10 @@
      "0.534422862873d+01" "0.518376171418d+01" "0.530010742143d+01"
      "0.521497741699d+01" "0.501157244350d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (ugaussian-tail *rng-mt19937* 5.0d0)))))
+	    (ugaussian-tail rng 5.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.187022708773d+01"
    (ugaussian-tail-pdf 5.2d0 5.0d0)))
diff --git a/random/gaussian.lisp b/random/gaussian.lisp
index c43bcbf2..6f0bc6e5 100644
--- a/random/gaussian.lisp
+++ b/random/gaussian.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gaussian.lisp                             
-; description: Gaussian distribution                     
-; date:        Sun Jul 16 2006 - 22:09                   
-; author:      Liam M. Healy                             
-; modified:    Sun Jul 16 2006 - 22:42
-;********************************************************
-;;; $Id: $
+;; Gaussian distribution
+;; Liam Healy, Sun Jul 16 2006 - 22:09
+;; Time-stamp: <2008-02-03 10:21:12EST gaussian.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,29 +9,29 @@
   "gsl_ran_gaussian"
   (((generator generator) :pointer) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A Gaussian random variate, with mean zero and
-  standard deviation @var{sigma}.  The probability distribution for
+  standard deviation sigma.  The probability distribution for
   Gaussian random variates is
   p(x) dx = {1 \over \sqrt{2 \pi \sigma^2}} \exp (-x^2 / 2\sigma^2) dx
-  for @math{x} in the range @math{-\infty} to @math{+\infty}.  Use the
-  transformation @math{z = \mu + x} on the numbers returned by
-  @code{gsl_ran_gaussian} to obtain a Gaussian distribution with mean
-  @math{\mu}.  This function uses the Box-Mueller algorithm which requires two
-  calls to the random number generator @var{r}.")
+  for x in the range -\infty to +\infty.  Use the
+  transformation z = \mu + x on the numbers returned by
+  #'gaussian to obtain a Gaussian distribution with mean
+  mu.  This function uses the Box-Mueller algorithm which requires two
+  calls to the random number generator r.")
 
 (defun-gsl gaussian-pdf (x sigma)
   "gsl_ran_gaussian_pdf" ((x :double) (sigma :double))
   :c-return :double
-  :documentation
-  "Compute the probability density @math{p(x)} at @var{x}
-   for a Gaussian distribution with standard deviation @var{sigma}.")
+  :documentation			; FDL
+  "Compute the probability density p(x) at x
+   for a Gaussian distribution with standard deviation sigma.")
 
 (defun-gsl gaussian-ziggurat (generator sigma)
   "gsl_ran_gaussian_ziggurat"
   (((generator generator) :pointer) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Compute a Gaussian random variate using the alternative
    Marsaglia-Tsang ziggurat method. The Ziggurat algorithm
    is the fastest available algorithm in most cases.")
@@ -44,89 +40,89 @@
   "gsl_ran_gaussian_ratio_method"
   (((generator generator) :pointer) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Compute a Gaussian random variate using the Kinderman-Monahan-Leva
    ratio method.")
 
 (defun-gsl ugaussian (generator)
   "gsl_ran_ugaussian" (((generator generator) :pointer))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Compute results for the unit Gaussian distribution,
    equivalent to the #'gaussian with a standard deviation of one,
-   @var{sigma} = 1.")
+   sigma = 1.")
 
 (defun-gsl ugaussian-pdf (x)
   "gsl_ran_ugaussian_pdf" ((x :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Compute results for the unit Gaussian distribution,
    equivalent to the #'gaussian-pdf with a standard deviation of one,
-   @var{sigma} = 1.")
+   sigma = 1.")
 
 (defun-gsl ugaussian-ratio-method (generator)
   "gsl_ran_ugaussian_ratio_method"
   (((generator generator) :pointer))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "Compute results for the unit Gaussian distribution,
    equivalent to the #'gaussian-ration-method with a
-   standard deviation of one, @var{sigma} = 1.")
+   standard deviation of one, sigma = 1.")
 
 (defun-gsl gaussian-P (x sigma)
   "gsl_cdf_gaussian_P" ((x :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The cumulative distribution function @math{P(x)} for the Gaussian
-   distribution with standard deviation @var{sigma}.")
+  :documentation			; FDL
+  "The cumulative distribution function P(x) for the Gaussian
+   distribution with standard deviation sigma.")
 
 (defun-gsl gaussian-Q (x sigma)
   "gsl_cdf_gaussian_Q" ((x :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The cumulative distribution function @math{Q(x)} for the Gaussian
-   distribution with standard deviation @var{sigma}.")
+  :documentation			; FDL
+  "The cumulative distribution function Q(x) for the Gaussian
+   distribution with standard deviation sigma.")
 
 (defun-gsl gaussian-Pinv (P sigma)
   "gsl_cdf_gaussian_Pinv" ((P :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The inverse cumulative distribution function @math{P(x)} for the Gaussian
-   distribution with standard deviation @var{sigma}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function P(x) for the Gaussian
+   distribution with standard deviation sigma.")
 
 (defun-gsl gaussian-Qinv (Q sigma)
   "gsl_cdf_gaussian_Qinv" ((Q :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The inverse cumulative distribution function @math{Q(x)} for the Gaussian
-   distribution with standard deviation @var{sigma}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function Q(x) for the Gaussian
+   distribution with standard deviation sigma.")
 
 (defun-gsl ugaussian-P (x)
   "gsl_cdf_ugaussian_P" ((x :double))
   :c-return :double
-  :documentation
-  "The cumulative distribution function @math{P(x)} for the Gaussian
+  :documentation			; FDL
+  "The cumulative distribution function P(x) for the Gaussian
    distribution with unit standard deviation.")
 
 (defun-gsl ugaussian-Q (x)
   "gsl_cdf_ugaussian_Q" ((x :double))
   :c-return :double
-  :documentation
-  "The cumulative distribution function @math{Q(x)} for the Gaussian
+  :documentation			; FDL
+  "The cumulative distribution function Q(x) for the Gaussian
    distribution with unit standard deviation.")
 
 (defun-gsl ugaussian-Pinv (P)
   "gsl_cdf_ugaussian_Pinv" ((P :double))
   :c-return :double
-  :documentation
-  "The inverse cumulative distribution function @math{P(x)} for the Gaussian
+  :documentation			; FDL
+  "The inverse cumulative distribution function P(x) for the Gaussian
    distribution with unit standard deviation.")
 
 (defun-gsl ugaussian-Qinv (Q)
   "gsl_cdf_ugaussian_Qinv" ((Q :double))
   :c-return :double
-  :documentation
-  "The inverse cumulative distribution function @math{Q(x)} for the Gaussian
+  :documentation			; FDL
+  "The inverse cumulative distribution function Q(x) for the Gaussian
    distribution with unit standard deviation.")
 
 ;;; Examples and unit test
@@ -137,11 +133,10 @@
      "-0.239671528273d+02" "-0.679280164729d+01" "-0.390913184336d+00"
      "0.893555545521d+01" "-0.176477945898d+00")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
-	 collect
-	 (gaussian *rng-mt19937* 10.0d0)))))
+	    collect
+	    (gaussian rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.398942280401d-01"
    (gaussian-pdf 0.0d0 10.0d0))
@@ -151,11 +146,10 @@
      "-0.143332423262d+02" "0.945508871923d+00" "-0.111203439859d+02"
      "-0.307919574280d+01" "-0.209068077307d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
-	 collect
-	 (gaussian-ziggurat *rng-mt19937* 10.0d0)))))
+	    collect
+	    (gaussian-ziggurat rng 10.0d0)))))
   ;; Given in examples in GSL documentation
   (lisp-unit:assert-first-fp-equal "0.977249868052d+00" (ugaussian-p 2.0d0))
   (lisp-unit:assert-first-fp-equal "0.227501319482d-01" (ugaussian-q 2.0d0))
diff --git a/random/generators.lisp b/random/generators.lisp
index 7b46f30b..52f412db 100644
--- a/random/generators.lisp
+++ b/random/generators.lisp
@@ -1,6 +1,6 @@
 ;; Generators of random numbers.
 ;; Liam Healy, Sat Jul 15 2006 - 14:43
-;; Time-stamp: <2008-02-02 23:01:02EST generators.lisp>
+;; Time-stamp: <2008-02-03 11:36:39EST generators.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -244,10 +244,6 @@
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-;;; (defparameter *rng-mt19937* (make-random-number-generator *mt19937*))
-;;; (defparameter *rng-cmrg* (make-random-number-generator *cmrg*))
-;;; (defparameter *rng-default* (make-random-number-generator *default-type*))
-
 (lisp-unit:define-test random-number-generators
   (lisp-unit:assert-equal
    '(999 162 282 947 231 484 957 744 540 739 759)
diff --git a/random/geometric.lisp b/random/geometric.lisp
index b94dbac3..fde5be59 100644
--- a/random/geometric.lisp
+++ b/random/geometric.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        geometric.lisp                          
-; description: Geometric distribution
-; date:        Sat Nov 25 2006 - 16:00
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 21:56
-;********************************************************
-;;; $Id: $
+;; Geometric distribution
+;; Liam Healy, Sat Nov 25 2006 - 16:00
+;; Time-stamp: <2008-02-03 11:04:00EST geometric.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,44 +9,45 @@
   "gsl_ran_geometric"
   (((generator generator) :pointer) (p :double))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "A random integer from the geometric distribution,
-   the number of independent trials with probability @var{p} until the
+   the number of independent trials with probability p until the
    first success.  The probability distribution for geometric variates
-   is p(k) =  p (1-p)^{k-1} for @math{k >= 1}.
-   Note that the distribution begins with @math{k=1} with this
-   definition.  There is another convention in which the exponent @math{k-1} 
-   is replaced by @math{k}.")
+   is p(k) =  p (1-p)^{k-1} for k >= 1.
+   Note that the distribution begins with k=1 with this
+   definition.  There is another convention in which the exponent k-1
+   is replaced by k.")
 
 (defun-gsl geometric-pdf (k p)
   "gsl_ran_geometric_pdf" ((k :uint) (p :double))
   :c-return :double
-  :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a geometric distribution with probability parameter @var{p}, using
+  :documentation			; FDL
+  "The probability p(k) of obtaining k
+   from a geometric distribution with probability parameter p, using
    the formula given in #'geometric.")
 
 (defun-gsl geometric-P (k p)
   "gsl_cdf_geometric_P" ((k :uint) (p :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(k)} for the Geometric distribution with parameter @var{p}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(k) for the geometric distribution with parameter p.")
 
 (defun-gsl geometric-Q (k p)
   "gsl_cdf_geometric_Q" ((k :uint) (p :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(k)} for the Geometric distribution with parameters @var{p}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(k) for the geometric distribution with parameters p.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test geometric
   (lisp-unit:assert-equal
    '(1 4 3 1 3 2 1 1 2 1 1)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(geometric *rng-mt19937* 0.4d0))))
+	   collect
+	   (geometric rng 0.4d0))))
   (lisp-unit:assert-first-fp-equal
    "0.240000000000d+00"
    (geometric-pdf 2 0.4d0))
diff --git a/random/gumbel1.lisp b/random/gumbel1.lisp
index d259d6f0..61f207ac 100644
--- a/random/gumbel1.lisp
+++ b/random/gumbel1.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gumbel1.lisp
-; description: Beta distribution                  
-; date:        Sun Oct 29 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 11 2006 - 20:58
-;********************************************************
-;;; $Id: $
+;; The Gumbel type 1 random number distribution
+;; Liam Healy, Sun Oct 29 2006
+;; Time-stamp: <2008-02-03 09:46:20EST gumbel1.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,47 +9,51 @@
   "gsl_ran_gumbel1"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Type-1 Gumbel
    distribution,
    p(x) dx = a b \exp(-(b \exp(-ax) + ax)) dx
-   for @math{-\infty < x < \infty}.")
+   for -\infty < x < \infty.")
 
 (defun-gsl gumbel1-pdf (x a b)
   "gsl_ran_gumbel1_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-  for a Type-1 Gumbel distribution with parameters @var{a} and @var{b},
+  :documentation			; FDL
+  "The probability density p(x) at x
+  for a Type-1 Gumbel distribution with parameters a and b,
   using the formula given for #'gumbel1.")
 
 (defun-gsl gumbel1-P (x a b)
   "gsl_cdf_gumbel1_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the Type-1 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the Type-1 Gumbel distribution with
+  parameters a and b.")
 
 (defun-gsl gumbel1-Q (x a b)
   "gsl_cdf_gumbel1_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the Type-1 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the Type-1 Gumbel distribution with
+  parameters a and b.")
 
 (defun-gsl gumbel1-Pinv (P a b)
   "gsl_cdf_gumbel1_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the Type-1 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the Type-1 Gumbel distribution with
+  parameters a and b.")
 
 (defun-gsl gumbel1-Qinv (Q a b)
   "gsl_cdf_gumbel1_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{Q(x)} for the Type-1 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  Q(x) for the Type-1 Gumbel distribution with
+  parameters a and b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test gumbel1
diff --git a/random/gumbel2.lisp b/random/gumbel2.lisp
index 7a4be3cd..23109738 100644
--- a/random/gumbel2.lisp
+++ b/random/gumbel2.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gumbel2.lisp
-; description: Beta distribution                  
-; date:        Sun Oct 29 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 11 2006 - 20:58
-;********************************************************
-;;; $Id: $
+;; The Gumbel type 2 random number distribution
+;; Liam Healy, Sun Oct 29 2006
+;; Time-stamp: <2008-02-03 10:10:26EST gumbel2.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,47 +9,49 @@
   "gsl_ran_gumbel2"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Type-2 Gumbel
-   distribution,
-   p(x) dx = a b x^@{-a-1@} \exp(-b x^@{-a@}) dx
-   for @math{0 < x < \infty}.")
+   distribution, p(x) dx = a b x^{-a-1} \exp(-b x^{-a}) dx
+   for 0 < x < \infty.")
 
 (defun-gsl gumbel2-pdf (x a b)
   "gsl_ran_gumbel2_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Type-2 Gumbel distribution with parameters @var{a} and @var{b},
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Type-2 Gumbel distribution with parameters a and b,
    using the formula given in #'gumbel2.")
 
 (defun-gsl gumbel2-P (x a b)
   "gsl_cdf_gumbel2_P" ((x :double) (a :double) (b :double))
-  :c-return :double
+  :c-return :double			; FDL
   :documentation "The cumulative distribution functions
-  @math{P(x)} for the Type-2 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  P(x) for the Type-2 Gumbel distribution with
+  parameters a and b.")
 
 (defun-gsl gumbel2-Q (x a b)
   "gsl_cdf_gumbel2_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the Type-2 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the Type-2 Gumbel distribution with
+  parameters a and b.")
 
 (defun-gsl gumbel2-Pinv (P a b)
   "gsl_cdf_gumbel2_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the Type-2 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the Type-2 Gumbel distribution with
+  parameters a and b.")
 
 (defun-gsl gumbel2-Qinv (Q a b)
   "gsl_cdf_gumbel2_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{Q(x)} for the Type-2 Gumbel distribution with
-  parameters @var{a} and @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  Q(x) for the Type-2 Gumbel distribution with
+  parameters a and b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test gumbel2
@@ -63,11 +61,10 @@
      "0.460260800603d+02" "0.677268352507d+01" "0.324619836866d+01"
      "0.664079780729d+01" "0.728568789702d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
-	 collect
-	 (gumbel2 *rng-mt19937* 1.0d0 2.0d0)))))
+	    collect
+	    (gumbel2 rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.536256036829d-01"
    (gumbel2-pdf 5.0d0 1.0d0 2.0d0))
diff --git a/random/hypergeometric.lisp b/random/hypergeometric.lisp
index fc1aa1b0..242c48b7 100644
--- a/random/hypergeometric.lisp
+++ b/random/hypergeometric.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        hypergeometric.lisp                          
-; description: Hypergeometric distribution
-; date:        Sat Nov 25 2006 - 16:00
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 22:40
-;********************************************************
-;;; $Id: $
+;; Hypergeometric distribution
+;; Liam Healy, Sat Nov 25 2006 - 16:00
+;; Time-stamp: <2008-02-03 10:12:46EST hypergeometric.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,53 +9,50 @@
   "gsl_ran_hypergeometric"
   (((generator generator) :pointer) (n1 :uint) (n2 :uint)(tt :uint))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "A random integer from the hypergeometric
    distribution.  The probability distribution for hypergeometric
    random variates is
    p(k) =  C(n_1, k) C(n_2, t - k) / C(n_1 + n_2, t)
-   where @math{C(a,b) = a!/(b!(a-b)!)} and 
-   @math{t <= n_1 + n_2}.  The domain of @math{k} is 
-   @math{max(0,t-n_2), ..., min(t,n_1)}.
-   If a population contains @math{n_1} elements of ``type 1'' and
-   @math{n_2} elements of ``type 2'' then the hypergeometric
-   distribution gives the probability of obtaining @math{k} elements of
-   ``type 1'' in @math{t} samples from the population without
+   where C(a,b) = a!/(b!(a-b)!) and 
+   t <= n_1 + n_2.  The domain of k is 
+   max(0,t-n_2), ..., min(t,n_1).
+   If a population contains n_1 elements of ``type 1'' and
+   n_2 elements of ``type 2'' then the hypergeometric
+   distribution gives the probability of obtaining k elements of
+   ``type 1'' in t samples from the population without
    replacement.")
 
 (defun-gsl hypergeometric-pdf (k n1 n2 tt)
   "gsl_ran_hypergeometric_pdf" ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
   :c-return :double
-  :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a hypergeometric distribution with parameters @var{n1}, @var{n2},
-   @var{tt}, using the formula given in #'hypergeometric.")
+  :documentation			; FDL
+  "The probability p(k) of obtaining k
+   from a hypergeometric distribution with parameters n1, n2,
+   tt, using the formula given in #'hypergeometric.")
 
 (defun-gsl hypergeometric-P (k n1 n2 tt)
   "gsl_cdf_hypergeometric_P" ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
   :c-return :double
-  :documentation
-  "The cumulative distribution functions @math{P(k)} for the
-   hypergeometric distribution with parameters @var{n1}, @var{n2}
-   and @var{tt}.")
+  :documentation			; FDL
+  "The cumulative distribution functions P(k) for the
+   hypergeometric distribution with parameters n1, n2 and tt.")
 
 (defun-gsl hypergeometric-Q  (k n1 n2 tt)
   "gsl_cdf_hypergeometric_Q"  ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
   :c-return :double
-  :documentation
-  "The cumulative distribution functions @math{Q(k)} for the
-   hypergeometric distribution with parameters @var{n1}, @var{n2}
-   and @var{tt}.")
+  :documentation			; FDL
+  "The cumulative distribution functions Q(k) for the
+   hypergeometric distribution with parameters n1, n2, and tt.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test hypergeometric-randist
   (lisp-unit:assert-equal
    '(2 1 0 0 1 1 3 1 0 1 3)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(hypergeometric *rng-mt19937* 3 6 3))))
+	   collect
+	   (hypergeometric rng 3 6 3))))
   (lisp-unit:assert-first-fp-equal
    "0.357142857143d+00"
    (hypergeometric-pdf 0 2 6 3))
diff --git a/random/landau.lisp b/random/landau.lisp
index 9d63ffa6..d15104b1 100644
--- a/random/landau.lisp
+++ b/random/landau.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        landau.lisp                          
-; description: Landau distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Sep 30 2006 - 21:52
-;********************************************************
-;;; $Id: $
+;; Landau distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 11:30:46EST landau.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,7 +9,7 @@
   "gsl_ran_landau"
   (((generator generator) :pointer))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Landau distribution.  The
    probability distribution for Landau random variates is defined
    analytically by the complex integral,
@@ -25,8 +21,8 @@
 (defun-gsl landau-pdf (x)
   "gsl_ran_landau_pdf" ((x :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
+  :documentation			; FDL
+  "The probability density p(x) at x
    for the Landau distribution using an approximation to the formula given
    in #'landau.")
 
@@ -38,11 +34,9 @@
      "0.261461684796d+02" "0.433721764097d+01" "0.167995462811d+01"
      "0.424757192183d+01" "0.468150620898d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
-	    collect
-	    (landau *rng-mt19937*)))))
+	    collect (landau rng)))))
   (lisp-unit:assert-first-fp-equal
    "0.173319689959d+00"
    (landau-pdf 0.25d0)))
diff --git a/random/laplace.lisp b/random/laplace.lisp
index 778b9893..343ddf32 100644
--- a/random/laplace.lisp
+++ b/random/laplace.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        exponential.lisp                          
-; description: Exponential distribution                  
-; date:        Sun Sep 17 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Sep 30 2006 - 19:17
-;********************************************************
-;;; $Id: $
+;; Exponential distribution
+;; Liam Healy, Sun Sep 17 2006
+;; Time-stamp: <2008-02-03 11:06:31EST laplace.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,43 +9,47 @@
   "gsl_ran_laplace"
   (((generator generator) :pointer) (a :double))
   :c-return :double
-  :documentation
-  "A random variate from the Laplace distribution with width @var{a}.
+  :documentation			; FDL
+  "A random variate from the Laplace distribution with width a.
    The distribution is
    p(x) dx = {1 \over 2 a}  \exp(-|x/a|) dx
-   for @math{-\infty < x < \infty}.")
+   for -\infty < x < \infty.")
 
 (defun-gsl laplace-pdf (x a)
   "gsl_ran_laplace_pdf" ((x :double) (a :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Laplace distribution with width @var{a}, using the formula
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Laplace distribution with width a, using the formula
    given for #'laplace.")
 
 (defun-gsl laplace-P (x a)
   "gsl_cdf_laplace_P" ((x :double) (a :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-   @math{P(x)} for the laplace distribution with width @var{a}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+   P(x) for the laplace distribution with width a.")
 
 (defun-gsl laplace-Q (x a)
   "gsl_cdf_laplace_Q" ((x :double) (a :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-   @math{Q(x)} for the laplace distribution with width @var{a}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+   Q(x) for the laplace distribution with width a.")
 
 (defun-gsl laplace-Pinv (P a)
   "gsl_cdf_laplace_Pinv" ((P :double) (a :double))
   :c-return :double
-  :documentation "The inverse cumulative distribution function
-   @math{P(x)} for the laplace distribution with width @var{a}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function
+   P(x) for the laplace distribution with width a.")
 
 (defun-gsl laplace-Qinv (Q a)
   "gsl_cdf_laplace_Qinv" ((Q :double) (a :double))
   :c-return :double
-  :documentation "The inverse cumulative distribution function
-   @math{Q(x)} for the laplace distribution with width @var{a}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function
+   Q(x) for the laplace distribution with width a.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test laplace
@@ -59,11 +59,10 @@
      "0.888815832003d+00" "0.716189249197d+01" "0.252463778091d+02"
      "0.734165104806d+01" "0.654142651602d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (laplace *rng-mt19937* 10.0d0)))))
+	    (laplace rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.500000000000d-01"
    (laplace-pdf 0.0d0 10.0d0))
diff --git a/random/levy.lisp b/random/levy.lisp
index 8d5eb22f..969887b1 100644
--- a/random/levy.lisp
+++ b/random/levy.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        levy.lisp                          
-; description: Levy distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Sep 30 2006 - 22:45
-;********************************************************
-;;; $Id: $
+;; Levy distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 10:44:49EST levy.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,39 +9,38 @@
   "gsl_ran_levy"
   (((generator generator) :pointer) (c :double) (alpha :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Levy symmetric stable
-   distribution with scale @var{c} and exponent @var{alpha}.  The symmetric
-   stable probability distribution is defined by a fourier transform
+   distribution with scale c and exponent alpha.  The symmetric
+   stable probability distribution is defined by a fourier transform,
    p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^\alpha)
-   There is no explicit solution for the form of @math{p(x)} and the
-   library does not define a corresponding @code{pdf} function.  For
-   @math{\alpha = 1} the distribution reduces to the Cauchy distribution.  For
-   @math{\alpha = 2} it is a Gaussian distribution with @c{$\sigma = \sqrt{2} c$} 
-   @math{\sigma = \sqrt@{2@} c}.  For @math{\alpha < 1} the tails of the
-   distribution become extremely wide.
-   The algorithm only works for @c{$0 < \alpha \le 2$} @math{0 < alpha <= 2}.")
+   There is no explicit solution for the form of p(x) and the
+   library does not define a corresponding pdf function.  For
+   \alpha = 1 the distribution reduces to the Cauchy distribution.  For
+   \alpha = 2 it is a Gaussian distribution with \sigma = \sqrt{2} c
+   For \alpha < 1 the tails of the distribution become extremely wide.
+   The algorithm only works for 0 < alpha <= 2.")
 
 (defun-gsl levy-skew (generator c alpha beta)
   "gsl_ran_levy_skew"
   (((generator generator) :pointer) (c :double) (alpha :double) (beta :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Levy skew stable
-   distribution with scale @var{c}, exponent @var{alpha} and skewness
-   parameter @var{beta}.  The skewness parameter must lie in the range
-   @math{[-1,1]}.  The Levy skew stable probability distribution is defined
+   distribution with scale c exponent alpha and skewness
+   parameter beta.  The skewness parameter must lie in the range
+   [-1,1].  The Levy skew stable probability distribution is defined
    by a fourier transform,
-   p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt \exp(-it x - |c t|^\alpha (1-i \beta \sign(t) \tan(\pi\alpha/2)))
-   When @math{\alpha = 1} the term @math{\tan(\pi \alpha/2)} is replaced by
-   @math{-(2/\pi)\log|t|}.  There is no explicit solution for the form of
-   @math{p(x)} and the library does not define a corresponding @code{pdf}
-   function.  For @math{\alpha = 2} the distribution reduces to a Gaussian
-   distribution with @c{$\sigma = \sqrt{2} c$} 
-   @math{\sigma = \sqrt@{2@} c} and the skewness parameter has no effect.  
-   For @math{\alpha < 1} the tails of the distribution become extremely
-   wide.  The symmetric distribution corresponds to @math{\beta = 0}.
-   The algorithm only works for @c{$0 < \alpha \le 2$} @math{0 < alpha <= 2}.")
+   p(x) = {1 \over 2 \pi} \int_{-\infty}^{+\infty} dt
+        \exp(-it x - |c t|^\alpha (1-i \beta \sign(t) \tan(\pi\alpha/2)))
+   When \alpha = 1 the term \tan(\pi \alpha/2) is replaced by
+   -(2/\pi)\log|t|.  There is no explicit solution for the form of
+   p(x)} and the library does not define a corresponding pdf
+   function.  For \alpha = 2 the distribution reduces to a Gaussian
+   distribution with \sigma = \sqrt{2} c and the skewness parameter
+   has no effect.   For \alpha < 1 the tails of the distribution
+   become extremely wide.  The symmetric distribution corresponds to \beta = 0.
+   The algorithm only works for 0 < \alpha \le 2.")
 
 
 ;;; Examples and unit test
@@ -56,19 +51,17 @@
      "-0.510713467479d+00" "0.164820785369d+00" "-0.148578990410d+00"
      "-0.190748857444d+01" "-0.208619521400d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (levy *rng-mt19937* 1.0d0 2.0d0)))))
+	    (levy rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-equal
    '("0.269410983324d+01" "-0.293954386447d+00" "-0.127034013523d+01"
      "0.107715386401d+01" "0.137712184069d+00" "0.941972843811d+00"
      "-0.510713467479d+00" "0.164820785369d+00" "-0.148578990410d+00"
      "-0.190748857444d+01" "-0.208619521400d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (levy-skew *rng-mt19937* 1.0d0 2.0d0 1.0d0))))))
+	    (levy-skew rng 1.0d0 2.0d0 1.0d0))))))
diff --git a/random/logarithmic.lisp b/random/logarithmic.lisp
index 8a722802..e034879a 100644
--- a/random/logarithmic.lisp
+++ b/random/logarithmic.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        logarithmic.lisp                          
-; description: Logarithmic distribution
-; date:        Sat Nov 25 2006 - 16:00
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 22:57
-;********************************************************
-;;; $Id: $
+;; Logarithmic distribution
+;; Liam Healy, Sat Nov 25 2006 - 16:00
+;; Time-stamp: <2008-02-03 11:10:27EST logarithmic.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,29 +9,28 @@
   "gsl_ran_logarithmic"
   (((generator generator) :pointer) (p :double))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "A random integer from the logarithmic distribution.
    The probability distribution for logarithmic random variates
    is p(k) = {-1 \over \log(1-p)} {\left( p^k \over k \right)}
-   for @math{k >= 1}.")
+   for k >= 1.")
 
 (defun-gsl logarithmic-pdf (k p)
   "gsl_ran_logarithmic_pdf" ((k :uint) (p :double))
   :c-return :double
   :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a logarithmic distribution with probability parameter @var{p},
+  "The probability p(k) of obtaining k
+   from a logarithmic distribution with probability parameter p,
    using the formula given in #'logarithmic.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test logarithmic
   (lisp-unit:assert-equal
    '(1 3 1 4 1 1 2 1 1 5 2)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(logarithmic *rng-mt19937* 0.9d0))))
+	   collect
+	   (logarithmic rng 0.9d0))))
   (lisp-unit:assert-first-fp-equal
    "0.156609215118d+00"
    (logarithmic-pdf 2 0.4d0)))
diff --git a/random/logistic.lisp b/random/logistic.lisp
index 53b31266..a61b4ddb 100644
--- a/random/logistic.lisp
+++ b/random/logistic.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        logistic.lisp                          
-; description: logistic distribution                  
-; date:        Sat Oct  7 2006 - 16:13
-; author:      Liam M. Healy                             
-; modified:    Sun Oct  8 2006 - 17:32
-;********************************************************
-;;; $Id: $
+;; Logistic distribution
+;; Liam Healy, Sat Oct  7 2006 - 16:13
+;; Time-stamp: <2008-02-03 10:35:08EST logistic.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,42 +9,46 @@
   "gsl_ran_logistic"
   (((generator generator) :pointer) (a :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the logistic distribution.  The distribution function is
    p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx
-   for @math{-\infty < x < +\infty}.")
+   for -\infty < x < +\infty.")
 
 (defun-gsl logistic-pdf (x a)
   "gsl_ran_logistic_pdf" ((x :double) (a :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a logistic distribution with scale parameter @var{a}, using the
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a logistic distribution with scale parameter a, using the
    formula given in #'logistic.")
 
 (defun-gsl logistic-P (x a)
   "gsl_cdf_logistic_P" ((x :double) (a :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the logistic distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the logistic distribution with scale parameter a.")
 
 (defun-gsl logistic-Q (x a)
   "gsl_cdf_logistic_Q" ((x :double) (a :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the logistic distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the logistic distribution with scale parameter a.")
 
 (defun-gsl logistic-Pinv (P a)
   "gsl_cdf_logistic_Pinv" ((P :double) (a :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the logistic distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the logistic distribution with scale parameter a.")
 
 (defun-gsl logistic-Qinv (Q a)
   "gsl_cdf_logistic_Qinv" ((Q :double) (a :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the logistic distribution with scale parameter @var{a}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the logistic distribution with scale parameter a.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test logistic
@@ -58,11 +58,10 @@
      "0.311425552636d+02" "0.106846737210d+02" "0.160518409540d+01"
      "0.104572419047d+02" "0.115237141063d+02")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (logistic *rng-mt19937* 10.0d0)))))
+	    (logistic rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.235003712202d+00"
    (logistic-pdf 0.5d0 1.0d0))
diff --git a/random/lognormal.lisp b/random/lognormal.lisp
index b23c9796..a0dd1170 100644
--- a/random/lognormal.lisp
+++ b/random/lognormal.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        lognormal.lisp                          
-; description: Lognormal distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Oct  7 2006 - 15:18
-;********************************************************
-;;; $Id: $
+;; Lognormal distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 11:08:34EST lognormal.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,47 +9,50 @@
   "gsl_ran_lognormal"
   (((generator generator) :pointer) (zeta :double) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the lognormal distribution.
    The distribution function is
    p(x) dx = {1 \over x \sqrt{2 \pi \sigma^2}} \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx
-   for @math{x > 0}.")
+   for x > 0.")
 
 (defun-gsl lognormal-pdf (x zeta sigma)
   "gsl_ran_lognormal_pdf" ((x :double) (zeta :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a lognormal distribution with parameters @var{zeta} and @var{sigma},
+  :documentation			; FDL
+  "The probability density p(x) at X
+   for a lognormal distribution with parameters zeta and sigma,
    using the formula given in #'lognormal.")
 
 (defun-gsl lognormal-P (x zeta sigma)
   "gsl_cdf_lognormal_P" ((x :double) (zeta :double) (sigma :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the lognormal distribution with parameters
-  @var{zeta} and @var{sigma}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the lognormal distribution with parameters zeta and sigma.")
 
 (defun-gsl lognormal-Q (x zeta sigma)
   "gsl_cdf_lognormal_Q" ((x :double) (zeta :double) (sigma :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the lognormal distribution with parameters
-  @var{zeta} and @var{sigma}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the lognormal distribution with parameters
+  zeta and sigma.")
 
 (defun-gsl lognormal-Pinv (P zeta sigma)
   "gsl_cdf_lognormal_Pinv" ((P :double) (zeta :double) (sigma :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the lognormal distribution with parameters
-  @var{zeta} and @var{sigma}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the lognormal distribution with parameters
+  zeta and sigma.")
 
 (defun-gsl lognormal-Qinv (Q zeta sigma)
   "gsl_cdf_lognormal_Qinv" ((Q :double) (zeta :double) (sigma :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the lognormal distribution with parameters
-   @var{zeta} and @var{sigma}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the lognormal distribution with parameters
+   zeta and sigma.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test lognormal
@@ -63,11 +62,10 @@
      "0.914620656772d+00" "0.727307901066d+00" "0.218018485218d+01"
      "0.389088566169d+01" "0.182184697889d+03")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (lognormal *rng-mt19937* 1.0d0 2.0d0)))))
+	    (lognormal rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.152898339657d+00"
    (lognormal-pdf 1.2d0 1.0d0 2.0d0))
diff --git a/random/multinomial.lisp b/random/multinomial.lisp
index 37ce552e..68c5c344 100644
--- a/random/multinomial.lisp
+++ b/random/multinomial.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        multinomial.lisp                          
-; description: Multinomial distribution
-; date:        Sat Nov 25 2006 - 16:00
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 21:16
-;********************************************************
-;;; $Id: $
+;; Multinomial distribution
+;; Liam Healy, Sat Nov 25 2006 - 16:00
+;; Time-stamp: <2008-02-03 10:28:38EST multinomial.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -18,63 +14,57 @@
    ;; technically, n should be a uint array, but integers work
    ((gsl-array n) :pointer))
   :c-return :void
-  :documentation
-  "Returns an array n of @var{K} random variates from a 
+  :documentation			; FDL
+  "Returns an array n of K random variates from a 
    multinomial distribution.  The sum of the array n is specified
    by sum=N.  The distribution function is
    P(n_1, n_2, ..., n_K) = 
    (N!/(n_1! n_2! ... n_K!)) p_1^n_1 p_2^n_2 ... p_K^n_K
-   where @math{(n_1, n_2, ..., n_K)} 
-   are nonnegative integers with 
-   @math{sum_@{k=1@}^K n_k = N},
-   and @math{(p_1, p_2, ..., p_K)}
-   is a probability distribution with @math{\sum p_i = 1}.  
-   If the array @var{p}[@var{K}] is not normalized then its entries will be
+   where (n_1, n_2, ..., n_K) are nonnegative integers with 
+   sum_@{k=1@}^K n_k = N, and (p_1, p_2, ..., p_K)
+   is a probability distribution with \sum p_i = 1.  
+   If the array p[K] is not normalized then its entries will be
    treated as weights and normalized appropriately.
    Random variates are generated using the conditional binomial method (see
-   C.S. David, @cite{The computer generation of multinomial random
-   variates}, Comp. Stat. Data Anal. 16 (1993) 205--217 for details).")
+   C.S. David, \"The computer generation of multinomial random
+   variates,\" Comp. Stat. Data Anal. 16 (1993) 205--217 for details).")
 
 (defun-gsl multinomial-pdf (p n)
   "gsl_ran_multinomial_pdf"
   (((dim0 p) :uint) ((gsl-array p) :pointer) ((gsl-array n) :pointer))
   :c-return :double
-  :documentation
-  "Compute the probability @math{P(n_1, n_2, ..., n_K)}
-   of sampling @var{n}[@var{K}] from a multinomial distribution 
-   with parameters @var{p}[@var{K}], using the formula given for #'multinomial.")
+  :documentation			; FDL
+  "Compute the probability P(n_1, n_2, ..., n_K)
+   of sampling n[K] from a multinomial distribution 
+   with parameters p[K], using the formula given for #'multinomial.")
 
 (defun-gsl multinomial-log-pdf (p n)
   "gsl_ran_multinomial_lnpdf"
   (((dim0 p) :uint) ((gsl-array p) :pointer) ((gsl-array n) :pointer))
   :c-return :double
-  :documentation
-  "Compute the natural logarithm of the probability @math{P(n_1, n_2, ..., n_K)}
-   of sampling @var{n}[@var{K}] from a multinomial distribution 
-   with parameters @var{p}[@var{K}], using the formula given for #'multinomial.")
+  :documentation			; FDL
+  "Compute the natural logarithm of the probability P(n_1, n_2, ..., n_K)
+   of sampling n[K] from a multinomial distribution 
+   with parameters p[K], using the formula given for #'multinomial.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test multinomial
   (lisp-unit:assert-equalp
    #(5 0 1 2)
-   (progn
-     (rng-set *rng-mt19937* 0)
-     (with-data (p vector-double 4)
-       (setf (data p) #(0.1d0 0.2d0 0.3d0 0.4d0))
-       (with-data (n vector-fixnum 4)
-	 (multinomial *rng-mt19937* 8 p n)
-	 (data n)))))
+   (letm ((rng (random-number-generator *mt19937* 0))
+	  (p (vector-double #(0.1d0 0.2d0 0.3d0 0.4d0)))
+	  (n (vector-fixnum 4)))
+     (multinomial rng 8 p n)
+     (data n)))
   (lisp-unit:assert-first-fp-equal
    "0.806400000000d-04"
-   (with-data (p vector-double 4)
-     (setf (data p) #(0.1d0 0.2d0 0.3d0 0.4d0))
-     (with-data (n vector-fixnum 4)
-       (setf (data n) #(5 0 1 2))
-       (multinomial-pdf p N))))
+   (letm ((p (vector-double #(0.1d0 0.2d0 0.3d0 0.4d0)))
+	  (n (vector-fixnum 4)))
+     (setf (data n) #(5 0 1 2))
+     (multinomial-pdf p N)))
   (lisp-unit:assert-first-fp-equal
    "-0.942551575364d+01"
-   (with-data (p vector-double 4)
-     (setf (data p) #(0.1d0 0.2d0 0.3d0 0.4d0))
-     (with-data (n vector-fixnum 4)
-       (setf (data n) #(5 0 1 2))
-       (multinomial-log-pdf p n)))))
+   (letm ((p (vector-double #(0.1d0 0.2d0 0.3d0 0.4d0)))
+	  (n (vector-fixnum 4)))
+     (setf (data n) #(5 0 1 2))
+     (multinomial-log-pdf p n))))
diff --git a/random/negative-binomial.lisp b/random/negative-binomial.lisp
index 3fed45e1..1a5b5fb0 100644
--- a/random/negative-binomial.lisp
+++ b/random/negative-binomial.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        negative-binomial.lisp                          
-; description: Negative binomial and Pascal distributions
-; date:        Sat Nov 25 2006 - 16:00
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 25 2006 - 21:42
-;********************************************************
-;;; $Id: $
+;; Negative binomial and Pascal distributions
+;; Liam Healy, Sat Nov 25 2006 - 16:00
+;; Time-stamp: <2008-02-03 11:25:48EST negative-binomial.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -17,35 +13,37 @@
   "gsl_ran_negative_binomial"
   (((generator generator) :pointer) (p :double) (n :double))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "A random integer from the negative binomial
-   distribution, the number of failures occurring before @var{n} successes
-   in independent trials with probability @var{p} of success.  The
+   distribution, the number of failures occurring before n successes
+   in independent trials with probability p of success.  The
    probability distribution for negative binomial variates is,
    p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k
-   Note that @math{n} is not required to be an integer.")
+   Note that n is not required to be an integer.")
 
 (defun-gsl negative-binomial-pdf (k p n)
   "gsl_ran_negative_binomial_pdf" ((k :uint) (p :double) (n :double))
   :c-return :double
-  :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a negative binomial distribution with parameters @var{p} and
-   @var{n}, using the formula given in #'negative-binomial.")
+  :documentation			; FDL
+  "The probability p(k) of obtaining k
+   from a negative binomial distribution with parameters p and
+   n, using the formula given in #'negative-binomial.")
 
 (defun-gsl negative-binomial-P (k p n)
   "gsl_cdf_negative_binomial_P" ((k :uint) (p :double) (n :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(k)} for the negative binomial distribution
-  with parameters @var{p} and @var{n}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(k) for the negative binomial distribution
+  with parameters p and n.")
 
 (defun-gsl negative-binomial-Q (k p n)
   "gsl_cdf_negative_binomial_Q" ((k :uint) (p :double) (n :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(k)} for the negative binomial distribution
-   with parameters @var{p} and @var{n}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(k) for the negative binomial distribution
+   with parameters p and n.")
 
 ;;;;****************************************************************************
 ;;;; Pascal
@@ -55,34 +53,36 @@
   "gsl_ran_pascal"
   (((generator generator) :pointer) (p :double) (n :uint))
   :c-return :uint
-  :documentation
+  :documentation			; FDL
   "A random integer from the Pascal distribution.  The
    Pascal distribution is simply a negative binomial distribution with an
-   integer value of @math{n}.
+   integer value of n.
    p(k) = {(n + k - 1)! \over k! (n - 1)! } p^n (1-p)^k
-   @math{k >= 0}.")
+   k >= 0.")
 
 (defun-gsl pascal-pdf (k p n)
   "gsl_ran_pascal_pdf" ((k :uint) (p :double) (n :uint))
   :c-return :double
-  :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a Pascal distribution with parameters @var{p} and
-   @var{n}, using the formula given in #'pascal.")
+  :documentation			; FDL
+  "The probability p(k) of obtaining k
+   from a Pascal distribution with parameters p and
+   n, using the formula given in #'pascal.")
 
 (defun-gsl pascal-P (k p n)
   "gsl_cdf_pascal_P" ((k :uint) (p :double) (n :uint))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(k)} for the Pascal distribution
-  with parameters @var{p} and @var{n}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(k) for the Pascal distribution
+  with parameters p and n.")
 
 (defun-gsl pascal-Q (k p n)
   "gsl_cdf_pascal_Q" ((k :uint) (p :double) (n :uint))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(k)} for the Pascal distribution
-   with parameters @var{p} and @var{n}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+   Q(k) for the Pascal distribution
+   with parameters p and n.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
@@ -91,11 +91,10 @@
 (lisp-unit:define-test negative-binomial
   (lisp-unit:assert-equal
    '(10 7 12 23 20 24 18 12 4 22 15)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(negative-binomial *rng-mt19937* 0.4d0 12.0d0))))
+	   collect
+	   (negative-binomial rng 0.4d0 12.0d0))))
   (lisp-unit:assert-first-fp-equal
    "0.569847670899d-02"
    (negative-binomial-pdf 5 0.4d0 12.0d0))
@@ -107,11 +106,10 @@
    (negative-binomial-Q 5 0.4d0 12.0d0))
   (lisp-unit:assert-equal
    '(10 7 12 23 20 24 18 12 4 22 15)
-   (progn
-     (rng-set *rng-mt19937* 0)
+   (letm ((rng (random-number-generator *mt19937* 0)))
      (loop for i from 0 to 10
-	collect
-	(pascal *rng-mt19937* 0.4d0 12))))
+	   collect
+	   (pascal rng 0.4d0 12))))
   (lisp-unit:assert-first-fp-equal
    "0.569847670899d-02"
    (pascal-pdf 5 0.4d0 12))
diff --git a/random/pareto.lisp b/random/pareto.lisp
index 465b126c..f269b185 100644
--- a/random/pareto.lisp
+++ b/random/pareto.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        pareto.lisp                          
-; description: Pareto distribution                  
-; date:        Sat Oct  8 2006 - 21:23
-; author:      Liam M. Healy                             
-; modified:    Sun Oct  8 2006 - 21:45
-;********************************************************
-;;; $Id: $
+;; Pareto distribution
+;; Liam Healy, Sat Oct  8 2006 - 21:23
+;; Time-stamp: <2008-02-03 10:33:14EST pareto.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,43 +9,47 @@
   "gsl_ran_pareto"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "A random variate from the Pareto distribution of order @var{a}.
+  :documentation			; FDL
+  "A random variate from the Pareto distribution of order a.
    The distribution function is
    p(x) dx = (a/b) / (x/b)^{a+1} dx
-   @math{x >= b}.")
+   x >= b.")
 
 (defun-gsl pareto-pdf (x a b)
   "gsl_ran_pareto_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Pareto distribution with exponent @var{a} and scale @var{b}, using
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Pareto distribution with exponent a and scale b, using
    the formula given in #'pareto.")
 
 (defun-gsl pareto-P (x a b)
   "gsl_cdf_pareto_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the Pareto distribution with exponent a and scale b.")
 
 (defun-gsl pareto-Q (x a b)
   "gsl_cdf_pareto_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the Pareto distribution with exponent a and scale b.")
 
 (defun-gsl pareto-Pinv (P a b)
   "gsl_cdf_pareto_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the Pareto distribution with exponent a and scale b.")
 
 (defun-gsl pareto-Qinv (Q a b)
   "gsl_cdf_pareto_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the Pareto distribution with exponent a and scale b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test pareto
@@ -59,11 +59,10 @@
      "0.208882311615d+01" "0.268706924980d+01" "0.370340428797d+01"
      "0.270287443943d+01" "0.263177356639d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (pareto *rng-mt19937* 1.0d0 2.0d0)))))
+	    (pareto rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.511603440571d+00"
    (pareto-pdf 1.5d0 1.3d0 1.0d0))
diff --git a/random/poisson.lisp b/random/poisson.lisp
index 70d5bf26..e12c9fa7 100644
--- a/random/poisson.lisp
+++ b/random/poisson.lisp
@@ -1,6 +1,6 @@
 ;; Poisson distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2008-02-02 23:02:43EST poisson.lisp>
+;; Time-stamp: <2008-02-03 09:50:16EST poisson.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -9,31 +9,33 @@
   "gsl_ran_poisson"
   (((generator generator) :pointer) (mu :double))
   :c-return :uint
-  :documentation
-  "A random integer from the Poisson distribution with mean @var{mu}.
+  :documentation			; FDL
+  "A random integer from the Poisson distribution with mean mu.
    The probability distribution for Poisson variates is
    p(k) = {\mu^k \over k!} \exp(-\mu)
-   @math{k >= 0}.")
+   k >= 0.")
 
 (defun-gsl poisson-pdf (k mu)
   "gsl_ran_poisson_pdf" ((k :uint) (mu :double))
   :c-return :double
-  :documentation
-  "The probability @math{p(k)} of obtaining @var{k}
-   from a Poisson distribution with mean @var{mu} using the formula
+  :documentation			; FDL
+  "The probability p(k) of obtaining k
+   from a Poisson distribution with mean mu using the formula
    given in #'poisson.")
 
 (defun-gsl poisson-P (k mu)
   "gsl_cdf_poisson_P" ((k :uint) (mu :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(k)} for the Poisson distribution with parameter @var{mu}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(k) for the Poisson distribution with parameter mu.")
 
 (defun-gsl poisson-Q (k mu)
   "gsl_cdf_poisson_Q" ((k :uint) (mu :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(k)} for the Poisson distribution with parameter @var{mu}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(k) for the Poisson distribution with parameter mu.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test poisson
diff --git a/random/rayleigh-tail.lisp b/random/rayleigh-tail.lisp
index a4671180..2d6f00ed 100644
--- a/random/rayleigh-tail.lisp
+++ b/random/rayleigh-tail.lisp
@@ -1,6 +1,6 @@
 ;; Rayleigh tail distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2008-02-02 23:00:39EST rayleigh-tail.lisp>
+;; Time-stamp: <2008-02-03 09:47:12EST rayleigh-tail.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -9,20 +9,20 @@
   "gsl_ran_rayleigh_tail"
   (((generator generator) :pointer) (a :double) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the tail of the Rayleigh
-  distribution with scale parameter @var{sigma} and a lower limit of
-  @var{a}.  The distribution is
+  distribution with scale parameter sigma and a lower limit of
+  a.  The distribution is
   p(x) dx = {x \over \sigma^2} \exp ((a^2 - x^2) /(2 \sigma^2)) dx
-  for @math{x > a}.")
+  for x > a.")
 
 (defun-gsl rayleigh-tail-pdf (x a sigma)
   "gsl_ran_rayleigh_tail_pdf" ((x :double) (a :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Rayleigh tail distribution with scale parameter @var{sigma} and
-   lower limit @var{a}, using the formula given in #'rayleigh-tail.")
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Rayleigh tail distribution with scale parameter sigma and
+   lower limit a, using the formula given in #'rayleigh-tail.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test rayleigh-tail
diff --git a/random/rayleigh.lisp b/random/rayleigh.lisp
index 693b7e83..7f66f1e3 100644
--- a/random/rayleigh.lisp
+++ b/random/rayleigh.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        rayleigh.lisp                          
-; description: Rayleigh distribution                  
-; date:        Sat Sep 30 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Sep 30 2006 - 19:41
-;********************************************************
-;;; $Id: $
+;; Rayleigh distribution
+;; Liam Healy, Sat Sep 30 2006
+;; Time-stamp: <2008-02-03 11:29:48EST rayleigh.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,47 +9,51 @@
   "gsl_ran_rayleigh"
   (((generator generator) :pointer) (sigma :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Rayleigh distribution with
-   scale parameter @var{sigma}.  The distribution is
+   scale parameter sigma.  The distribution is
    p(x) dx = {x \over \sigma^2} \exp(- x^2/(2 \sigma^2)) dx
-   for @math{x > 0}.")
+   for x > 0.")
 
 (defun-gsl rayleigh-pdf (x sigma)
   "gsl_ran_rayleigh_pdf" ((x :double) (sigma :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Rayleigh distribution with scale parameter @var{sigma}, using the
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Rayleigh distribution with scale parameter sigma, using the
    formula given for #'rayleigh.")
 
 (defun-gsl rayleigh-P (x sigma)
   "gsl_cdf_rayleigh_P" ((x :double) (sigma :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-  @math{P(x)} for the Rayleigh distribution with scale
-  parameter @var{sigma}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+  P(x) for the Rayleigh distribution with scale
+  parameter sigma.")
 
 (defun-gsl rayleigh-Q (x sigma)
   "gsl_cdf_rayleigh_Q" ((x :double) (sigma :double))
   :c-return :double
-  :documentation "The cumulative distribution function
-  @math{Q(x)} for the Rayleigh distribution with scale
-  parameter @var{sigma}.")
+  :documentation			; FDL
+  "The cumulative distribution function
+  Q(x) for the Rayleigh distribution with scale
+  parameter sigma.")
 
 (defun-gsl rayleigh-Pinv (P sigma)
   "gsl_cdf_rayleigh_Pinv" ((P :double) (sigma :double))
   :c-return :double
-  :documentation "The inverse cumulative distribution function
-  @math{P(x)} for the Rayleigh distribution with scale
-  parameter @var{sigma}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function
+  P(x)} for the Rayleigh distribution with scale
+  parameter sigma.")
 
 (defun-gsl rayleigh-Qinv (Q sigma)
   "gsl_cdf_rayleigh_Qinv" ((Q :double) (sigma :double))
   :c-return :double
-  :documentation "The inverse cumulative distribution function
-  @math{Q(x)} for the Rayleigh distribution with scale
-  parameter @var{sigma}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution function
+  Q(x) for the Rayleigh distribution with scale
+  parameter sigma.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test rayleigh
@@ -63,11 +63,10 @@
      "0.294800354467d+01" "0.768510144246d+01" "0.111004981321d+02"
      "0.776103902005d+01" "0.740959915506d+01")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (rayleigh *rng-mt19937* 10.0d0)))))
+	    (rayleigh rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.441248451292d+00"
    (rayleigh-pdf 0.5d0 1.0d0))
diff --git a/random/spherical-vector.lisp b/random/spherical-vector.lisp
index 3b8147eb..12a54f7e 100644
--- a/random/spherical-vector.lisp
+++ b/random/spherical-vector.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        spherical-vector.lisp                          
-; description: Spherical Vector distribution                  
-; date:        Sun Oct  22 2006
-; author:      Liam M. Healy                             
-; modified:    Sun Oct 22 2006 - 21:58
-;********************************************************
-;;; $Id: $
+;; Spherical Vector distribution
+;; Liam Healy, Sun Oct  22 2006
+;; Time-stamp: <2008-02-03 10:38:08EST spherical-vector.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -15,28 +11,28 @@
   "gsl_ran_dir_2d"
   (((generator generator) :pointer) (x :double) (y :double))
   :c-return :void
-  :documentation
-  "A random direction vector @math{v} = (@var{x},@var{y}) in
+  :documentation			; FDL
+  "A random direction vector v = (x,y) in
    two dimensions.  The vector is normalized such that
-   @math{|v|^2 = x^2 + y^2 = 1}.")
+   |v|^2 = x^2 + y^2 = 1.")
 
 (defun-gsl direction-2d-trig-method (generator)
   "gsl_ran_dir_2d_trig_method"
   (((generator generator) :pointer) (x :double) (y :double))
   :c-return :void
-  :documentation
-  "A random direction vector @math{v} = (@var{x},@var{y}) in
+  :documentation			; FDL
+  "A random direction vector v = (x,y) in
    two dimensions.  The vector is normalized such that
-   @math{|v|^2 = x^2 + y^2 = 1}.  Uses trigonometric functions.")
+   |v|^2 = x^2 + y^2 = 1.  Uses trigonometric functions.")
 
 (defun-gsl direction-3d (generator)
   "gsl_ran_dir_3d"
   (((generator generator) :pointer) (x :double) (y :double) (z :double))
   :c-return :void
-  :documentation
-  "A random direction vector @math{v} =
-  (@var{x},@var{y},@var{z}) in three dimensions.  The vector is normalized
-  such that @math{|v|^2 = x^2 + y^2 + z^2 = 1}.  The method employed is
+  :documentation			; FDL
+  "A random direction vector v =
+  (x,y,z) in three dimensions.  The vector is normalized
+  such that |v|^2 = x^2 + y^2 + z^2 = 1.  The method employed is
   due to Robert E. Knop (CACM 13, 326 (1970)), and explained in Knuth, v2,
   3rd ed, p136.  It uses the surprising fact that the distribution
   projected along any axis is actually uniform (this is only true for 3
@@ -47,10 +43,10 @@
   (((generator generator) :pointer) ((dim0 x) :size) ((gsl-array x) :pointer))
   :c-return :void
   :return (x)
-  :documentation
-  "A random direction vector @math{v = (x_1,x_2,...,x_n)} in @var{n} dimensions,
+  :documentation			; FDL
+  "A random direction vector v = (x_1,x_2,...,x_n) in n dimensions,
    where n is the length of the vector x passed in. The vector is normalized such that 
-   @math{|v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1}.  The method
+   |v|^2 = x_1^2 + x_2^2 + ... + x_n^2 = 1.  The method
    uses the fact that a multivariate gaussian distribution is spherically
    symmetric.  Each component is generated to have a gaussian distribution,
    and then the components are normalized.  The method is described by
@@ -65,29 +61,26 @@
      "0.457266229462d+00" "0.889329857473d+00" "-0.463256161598d+00"
      "-0.886224423462d+00")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 4
 	    append
-	    (multiple-value-list (direction-2d *rng-mt19937*))))))
+	    (multiple-value-list (direction-2d rng))))))
   (lisp-unit:assert-equal
    '("0.999998683521d+00" "-0.162263876311d-02" "0.520301010608d+00"
      "0.853982937980d+00" "-0.203512053104d+00" "0.979072440753d+00"
      "0.945475322749d+00" "-0.325693742761d+00" "0.115000339166d+00"
      "0.993365452385d+00")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 4
 	    append
-	    (multiple-value-list (direction-2d-trig-method *rng-mt19937*))))))
+	    (multiple-value-list (direction-2d-trig-method rng))))))
   (lisp-unit:assert-equal
    '("-0.912992575045d-01" "0.187821853572d+00" "0.977950610665d+00"
      "-0.905118296156d+00" "-0.506837644858d-01" "-0.422127973465d+00"
      "0.139937665360d+00" "0.838546262052d+00" "-0.526552576873d+00")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 2
 	    append
-	    (multiple-value-list (direction-3d *rng-mt19937*)))))))
+	    (multiple-value-list (direction-3d rng)))))))
diff --git a/random/tdist.lisp b/random/tdist.lisp
index 667f8f7c..9a88f88c 100644
--- a/random/tdist.lisp
+++ b/random/tdist.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        tdist.lisp                          
-; description: tdist distribution                  
-; date:        Sat Oct  7 2006 - 16:13
-; author:      Liam M. Healy                             
-; modified:    Sun Oct  8 2006 - 16:44
-;********************************************************
-;;; $Id: $
+;; Tdist distribution
+;; Liam Healy, Sat Oct  7 2006 - 16:13
+;; Time-stamp: <2008-02-03 11:27:51EST tdist.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,44 +9,48 @@
   "gsl_ran_tdist"
   (((generator generator) :pointer) (nu :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Student t-distribution.  The
    distribution function is,
    p(x) dx = {\Gamma((\nu + 1)/2) \over \sqrt{\pi \nu} \Gamma(\nu/2)}
    (1 + x^2/\nu)^{-(\nu + 1)/2} dx
-  for @math{-\infty < x < +\infty}.")
+   for -\infty < x < +\infty.")
 
 (defun-gsl tdist-pdf (x nu)
   "gsl_ran_tdist_pdf" ((x :double) (nu :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a t-distribution with @var{nu} degrees of freedom, using the formula
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a t-distribution with nu degrees of freedom, using the formula
    given in #'tdist.")
 
 (defun-gsl tdist-P (x nu)
   "gsl_cdf_tdist_P" ((x :double) (nu :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  P(x) for the tdist distribution with nu degrees of freedom.")
 
 (defun-gsl tdist-Q (x nu)
   "gsl_cdf_tdist_Q" ((x :double) (nu :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+   Q(x) for the tdist distribution with nu degrees of freedom.")
 
 (defun-gsl tdist-Pinv (P nu)
   "gsl_cdf_tdist_Pinv" ((P :double) (nu :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   P(x) for the tdist distribution with nu degrees of freedom.")
 
 (defun-gsl tdist-Qinv (Q nu)
   "gsl_cdf_tdist_Qinv" ((Q :double) (nu :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the tdist distribution with nu degrees of freedom.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test tdist
@@ -59,12 +59,11 @@
      "-0.160088628258d+01" "-0.170109355058d+01" "-0.437095974981d-01"
      "0.127611592766d+00" "-0.197312182555d-01" "-0.653466611720d+00"
      "0.203577132452d+00" "0.177650300478d+01")
-   (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
-      (loop for i from 0 to 10
-	    collect
-	    (tdist *rng-mt19937* 10.0d0)))))
+   (letm ((rng (random-number-generator *mt19937* 0)))
+     (rng-set rng 0)
+     (loop for i from 0 to 10
+	   collect
+	   (tdist rng 10.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.254647908947d+00"
    (tdist-pdf 0.5d0 1.0d0))
diff --git a/random/weibull.lisp b/random/weibull.lisp
index 44da1a63..c857e819 100644
--- a/random/weibull.lisp
+++ b/random/weibull.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        weibull.lisp                          
-; description: Weibull distribution                  
-; date:        Sun Oct 22 2006
-; author:      Liam M. Healy                             
-; modified:    Sat Nov 11 2006 - 20:59
-;********************************************************
-;;; $Id: $
+;; Weibull distribution
+;; Liam Healy, Sun Oct 22 2006
+;; Time-stamp: <2008-02-03 11:15:10EST weibull.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,42 +9,46 @@
   "gsl_ran_weibull"
   (((generator generator) :pointer) (a :double) (b :double))
   :c-return :double
-  :documentation
+  :documentation			; FDL
   "A random variate from the Weibull distribution.  The distribution function is
    p(x) dx = {b \over a^b} x^{b-1}  \exp(-(x/a)^b) dx
-   for @math{x >= 0}.")
+   for x >= 0.")
 
 (defun-gsl weibull-pdf (x a b)
   "gsl_ran_weibull_pdf" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation
-  "The probability density @math{p(x)} at @var{x}
-   for a Weibull distribution with scale @var{a} and exponent @var{b},
+  :documentation			; FDL
+  "The probability density p(x) at x
+   for a Weibull distribution with scale a and exponent b,
    using the formula given in #'weibull.")
 
 (defun-gsl weibull-P (x a b)
   "gsl_cdf_weibull_P" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{P(x)} for the Weibull distribution with scale @var{a} and exponent @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+   P(x) for the Weibull distribution with scale a and exponent b.")
 
 (defun-gsl weibull-Q (x a b)
   "gsl_cdf_weibull_Q" ((x :double) (a :double) (b :double))
   :c-return :double
-  :documentation "The cumulative distribution functions
-  @math{Q(x)} for the Weibull distribution with scale @var{a} and exponent @var{b}.")
+  :documentation			; FDL
+  "The cumulative distribution functions
+  Q(x) for the Weibull distribution with scale a and exponent b.")
 
 (defun-gsl weibull-Pinv (P a b)
   "gsl_cdf_weibull_Pinv" ((P :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-  @math{P(x)} for the Weibull distribution scale @var{a} and exponent @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+  P(x) for the Weibull distribution scale a and exponent b.")
 
 (defun-gsl weibull-Qinv (Q a b)
   "gsl_cdf_weibull_Qinv" ((Q :double) (a :double) (b :double))
   :c-return :double
-  :documentation  "The inverse cumulative distribution functions
-   @math{Q(x)} for the Weibull distribution exponent @var{a} and scale @var{b}.")
+  :documentation			; FDL
+  "The inverse cumulative distribution functions
+   Q(x) for the Weibull distribution exponent a and scale b.")
 
 ;;; Examples and unit test
 (lisp-unit:define-test weibull
@@ -58,11 +58,10 @@
      "0.208455329740d+00" "0.543418734407d+00" "0.784923750377d+00"
      "0.548788332013d+00" "0.523937780842d+00")
    (lisp-unit:fp-sequence
-    (progn
-      (rng-set *rng-mt19937* 0)
+    (letm ((rng (random-number-generator *mt19937* 0)))
       (loop for i from 0 to 10
 	    collect
-	    (weibull *rng-mt19937* 1.0d0 2.0d0)))))
+	    (weibull rng 1.0d0 2.0d0)))))
   (lisp-unit:assert-first-fp-equal
    "0.242631749722d+00"
    (weibull-pdf 1.5d0 1.3d0 1.0d0))
-- 
GitLab