diff --git a/basis-splines.lisp b/basis-splines.lisp
index e93a2fb8bdc21edc7c601b11613c4f9da8c26050..ed1154d03ccbf1bd16ef41ddc3fd4dc8f76f6a3f 100644
--- a/basis-splines.lisp
+++ b/basis-splines.lisp
@@ -1,6 +1,6 @@
 ;; Basis splines.
 ;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
-;; Time-stamp: <2009-02-21 16:46:31EST basis-splines.lisp>
+;; Time-stamp: <2009-05-24 22:58:41EDT basis-splines.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -109,7 +109,7 @@
     (dotimes (i ndata)
       (let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float))
 	     (yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
-		    (gaussian rng sigma))))
+		    (sample rng 'gaussian :sigma sigma))))
 	(setf (maref x i) xi
 	      (maref y i) yi
 	      (maref w i) (/ (expt sigma 2)))))
diff --git a/gsll.asd b/gsll.asd
index 3858069425ca18ef88d333b2d9eb541418b73df1..f9e59502a6216170adfc590596aa2c7924efec98 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -1,6 +1,6 @@
 ;; Definition of GSLL system 
 ;; Liam Healy
-;; Time-stamp: <2009-05-19 21:06:12EDT gsll.asd>
+;; Time-stamp: <2009-05-24 16:07:08EDT gsll.asd>
 ;; $Id$
 
 (when (asdf:find-system :fsbv nil)
@@ -172,7 +172,7 @@
 	     ;; minimum and maximum values provided in vector.lisp
 	     (:file "median-percentile")))
    (:module histogram
-	    :depends-on (init linear-algebra)
+	    :depends-on (init linear-algebra random)
 	    :components
 	    ((:file "histogram")
 	     (:file "updating-accessing" :depends-on (histogram))
diff --git a/histogram/probability-distribution.lisp b/histogram/probability-distribution.lisp
index cf5eba9e0024fc554dfc97e7ca57ff1ba69978b3..0498d658fc2f7127f3104bc7f1d5671e7da7cfdb 100644
--- a/histogram/probability-distribution.lisp
+++ b/histogram/probability-distribution.lisp
@@ -1,6 +1,6 @@
 ;; Histogram probability distribution.
 ;; Liam Healy, Mon Jan  1 2007 - 17:51
-;; Time-stamp: <2009-02-21 17:12:56EST probability-distribution.lisp>
+;; Time-stamp: <2009-05-24 16:28:15EDT probability-distribution.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -32,25 +32,28 @@
   :initialize-suffix "init"
   :initialize-args (((mpointer histogram) :pointer)))
 
-(export 'sample)
-(defgeneric sample (pdf value)
-  (:documentation ;; FDL
-   "Given a uniform random number between zero and one,
+(defmfun sample ((value number) (pdf histogram-pdf) &key)
+  "gsl_histogram_pdf_sample"
+  (((mpointer pdf) :pointer) (value :double))
+  :definition :method
+  :c-return :double
+  :documentation ;; FDL
+  "Given a uniform random number (source) between zero and one,
    compute a single random sample from the probability distribution
    'pdf.  The algorithm used to compute the sample s is given by
    s = range[i] + delta * (range[i+1] - range[i])
    where i is the index which satisfies 
    sum[i] <=  value < sum[i+1] and delta is 
-   (value - sum[i])/(sum[i+1] - sum[i])."))
-
-(defmfun sample ((pdf histogram-pdf) value)
-  "gsl_histogram_pdf_sample"
-  (((mpointer pdf) :pointer) (value :double))
-  :definition :method
-  :c-return :double)
+   (value - sum[i])/(sum[i+1] - sum[i]).")
 
-(defmfun sample ((pdf histogram2d-pdf) value)
+(defmfun sample ((value number) (pdf histogram2d-pdf) &key)
   "gsl_histogram2d_pdf_sample"
   (((mpointer pdf) :pointer) (value :double))
   :definition :method
   :c-return :double)
+
+(defmethod sample ((source random-number-generator) (pdf histogram-pdf) &key)
+  (sample (sample source 'uniform) pdf))
+
+(defmethod sample ((source random-number-generator) (pdf histogram2d-pdf) &key)
+  (sample (sample source 'uniform) pdf))
diff --git a/init/defmfun.lisp b/init/defmfun.lisp
index b044d6e87f036355483d3d03390461b82db4c115..d6f76dc572cf9a990254ebe9acde184dcde7ed60 100644
--- a/init/defmfun.lisp
+++ b/init/defmfun.lisp
@@ -1,6 +1,6 @@
 ;; Macro for defining GSL functions.
 ;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
-;; Time-stamp: <2009-05-03 15:38:45EDT defmfun.lisp>
+;; Time-stamp: <2009-05-24 16:54:33EDT defmfun.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -171,7 +171,7 @@
     (if (listp gsl-name)
 	(mapc (lambda (n) (push n indexed-functions)) gsl-name)
 	(push gsl-name indexed-functions))
-    (remf key-args :documentation)
+    ;;(remf key-args :documentation)
       (complete-definition
        'cl:defmethod
        name
diff --git a/random/bernoulli.lisp b/random/bernoulli.lisp
index 2aeaad1a9bc2b92364c93987632fcc9a156482bb..34d51fbbc881cf5f6811a312ab6d708c66fdd85d 100644
--- a/random/bernoulli.lisp
+++ b/random/bernoulli.lisp
@@ -1,13 +1,15 @@
 ;; Bernoulli distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:59
-;; Time-stamp: <2009-02-16 10:08:17EST bernoulli.lisp>
+;; Time-stamp: <2009-05-24 22:16:20EDT bernoulli.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
-(defmfun bernoulli (generator p)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'bernoulli)) &key p)
   "gsl_ran_bernoulli"
   (((mpointer generator) :pointer) (p :double))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "Returns either 0 or 1, the result of a Bernoulli trial
@@ -29,5 +31,5 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (bernoulli rng 0.5d0)))
+	   (sample rng 'bernoulli :p 0.5d0)))
   (bernoulli-pdf 0 0.5d0))
diff --git a/random/beta.lisp b/random/beta.lisp
index a58e5356c6a24ea4c925e5a1b229148fce12afd3..0a4ef57e78e64aae9735040b78d8a984273e40a2 100644
--- a/random/beta.lisp
+++ b/random/beta.lisp
@@ -1,14 +1,15 @@
 ;; Beta distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:10:32EST beta.lisp>
+;; Time-stamp: <2009-05-24 21:34:36EDT beta.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
-(defmfun beta-rd (generator a b)
-  ;; Named #'beta-rd to avoid confusion with the special function #'beta.
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'beta)) &key a b)
   "gsl_ran_beta"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the beta distribution.  The distribution function is
@@ -56,7 +57,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (beta-rd rng 1.0d0 2.0d0)))
+	    (sample rng 'beta :a 1.0d0 :b 2.0d0)))
   (beta-pdf 0.1d0 1.0d0 2.0d0)
   (beta-P 0.1d0 1.0d0 2.0d0)
   (beta-Q 0.1d0 1.0d0 2.0d0)
diff --git a/random/cauchy.lisp b/random/cauchy.lisp
index 8d82c56f3ed24d94a7af9ee13c1d0402ee087536..47c23eb9bb488f0b3950597e1ab399904e476099 100644
--- a/random/cauchy.lisp
+++ b/random/cauchy.lisp
@@ -1,6 +1,6 @@
 ;; Cauchy distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:08:14EST cauchy.lisp>
+;; Time-stamp: <2009-05-24 20:04:40EDT cauchy.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun cauchy (generator a)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'cauchy)) &key a)
   "gsl_ran_cauchy"
   (((mpointer generator) :pointer) (a :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Cauchy distribution with
@@ -61,7 +63,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (cauchy rng 10.0d0)))
+	    (sample rng 'cauchy :a 10.0d0)))
   (cauchy-pdf 0.0d0 10.0d0)
   (cauchy-P 1.0d0 2.0d0)
   (cauchy-Q 1.0d0 2.0d0)
diff --git a/random/chi-squared.lisp b/random/chi-squared.lisp
index bcfacb68119cecf1ed69a41d19d09c659f57a782..a8ee5d3f906181fa6b829add101243d0c5263379 100644
--- a/random/chi-squared.lisp
+++ b/random/chi-squared.lisp
@@ -1,6 +1,6 @@
 ;; Chi-squared distribution
 ;; Liam Healy, Sat Oct  7 2006 - 16:13
-;; Time-stamp: <2009-02-16 10:10:33EST chi-squared.lisp>
+;; Time-stamp: <2009-05-24 20:23:33EDT chi-squared.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun chi-squared (generator nu)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'chi-squared)) &key nu)
   "gsl_ran_chisq"
   (((mpointer generator) :pointer) (nu :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the chi-squared distribution
@@ -59,7 +61,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (chi-squared rng 10.0d0)))
+	    (sample rng 'chi-squared :nu 10.0d0)))
   (chi-squared-pdf 0.5d0 1.0d0)
   (chi-squared-P 0.5d0 1.0d0)
   (chi-squared-Q 0.5d0 1.0d0)
diff --git a/random/dirichlet.lisp b/random/dirichlet.lisp
index d25b512ea734681f14e1f4c61e1e10d8e087893a..b94e25b592bb022c339d67497edbcdb40d7b2b05 100644
--- a/random/dirichlet.lisp
+++ b/random/dirichlet.lisp
@@ -1,19 +1,22 @@
 ;; Dirichlet distribution
 ;; Liam Healy, Sun Oct 29 2006
-;; Time-stamp: <2009-02-16 10:10:31EST dirichlet.lisp>
+;; Time-stamp: <2009-05-24 22:55:07EDT dirichlet.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun dirichlet (generator alpha theta)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'dirichlet))
+     &key alpha (theta (vdf (dim0 alpha))))
   "gsl_ran_dirichlet"
   (((mpointer generator) :pointer)
    ((dim0 alpha) sizet)
    ((c-pointer alpha) :pointer)
    ;; theta had better be at least as long as alpha, or they'll be trouble
    ((c-pointer theta) :pointer))
+  :definition :method
   :inputs (alpha)
   :outputs (theta)
   :c-return :void
@@ -64,10 +67,8 @@
 ;;; Examples and unit test
 (save-test dirichlet
  (let ((rng (make-random-number-generator +mt19937+ 0))
-       (alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
-       (theta (make-marray 'double-float :dimensions 4)))
-   (dirichlet rng alpha theta)
-   (cl-array theta))
+       (alpha #m(1.0d0 2.0d0 3.0d0 4.0d0)))
+   (cl-array (sample rng 'dirichlet :alpha alpha)))
  (let ((alpha #m(1.0d0 2.0d0 3.0d0 4.0d0))
        (theta #m(1.0d0 2.0d0 3.0d0 4.0d0)))
    (dirichlet-pdf alpha theta))
diff --git a/random/discrete.lisp b/random/discrete.lisp
index e65d0ff884e3aa6086e0b6e6c2cddf8730aa60e3..7e9a5f31a78c9886af0365847c1e2ca57ea39ee2 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: <2009-02-16 10:08:14EST discrete.lisp>
+;; Time-stamp: <2009-05-24 21:56:42EDT discrete.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -20,9 +20,12 @@
   generally as ``weights'')---the preprocessor will normalize appropriately.
   This return value is used as an argument to #'discrete.")
 
-(defmfun discrete (generator table)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'discrete))
+     &key table)
   "gsl_ran_discrete"
   (((mpointer generator) :pointer) ((mpointer table) :pointer))
+  :definition :method
   :c-return sizet
   :documentation
   "Generate discrete random numbers.")
@@ -46,7 +49,7 @@
 	(rng (make-random-number-generator +mt19937+ 0)))
    (loop for i from 0 to 10
       collect
-      (discrete rng table)))
+      (sample rng 'discrete :table table)))
  (let* ((probabilities #m(0.25d0 0.5d0 0.25d0))
 	(table (make-discrete-random probabilities)))
    (discrete-pdf 1 table)))
diff --git a/random/exponential-power.lisp b/random/exponential-power.lisp
index 7f93fcb58b02326ff781af4ee1a2ab35298e2b9e..e46b5c27a112d388ae4949c85a43ac47918288e1 100644
--- a/random/exponential-power.lisp
+++ b/random/exponential-power.lisp
@@ -1,15 +1,18 @@
 ;; Exponential power distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:08:16EST exponential-power.lisp>
+;; Time-stamp: <2009-05-24 20:02:56EDT exponential-power.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun exponential-power (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'exponential-power))
+     &key a b)
   "gsl_ran_exppow"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the exponential power distribution
@@ -49,8 +52,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (exponential-power rng 1.0d0 2.0d0)))
+	    (sample rng 'exponential-power :a 1.0d0 :b 2.0d0)))
   (exponential-power-pdf 0.0d0 1.0d0 2.0d0)
   (exponential-power-P 1.0d0 1.0d0 2.0d0)
   (exponential-power-Q 1.0d0 1.0d0 2.0d0))
-
diff --git a/random/exponential.lisp b/random/exponential.lisp
index 5e8451dd86f5b158c45e16a8ff4065fe7acb249a..a1459fc9f2b45d3c4bfc7155313fc32c33e8bc76 100644
--- a/random/exponential.lisp
+++ b/random/exponential.lisp
@@ -1,6 +1,6 @@
 ;; Exponential distribution
 ;; Liam Healy, Sat Sep  2 2006 - 19:04
-;; Time-stamp: <2009-02-16 10:08:11EST exponential.lisp>
+;; Time-stamp: <2009-05-24 19:31:30EDT exponential.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,12 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun exponential (generator mu)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'exponential))
+     &key mu)
   "gsl_ran_exponential"
   (((mpointer generator) :pointer) (mu :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the exponential distribution
@@ -59,7 +62,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (exponential rng 10.0d0)))
+	    (sample rng 'exponential :mu 10.0d0)))
   (exponential-pdf 0.0d0 10.0d0)
   (exponential-p 1.0d0 2.0d0)
   (exponential-q 1.0d0 2.0d0)
diff --git a/random/fdist.lisp b/random/fdist.lisp
index 38bbc88b830bf3171eed672b683954bac2b68345..665715320af958510f2b11d8565d6c648f794492 100644
--- a/random/fdist.lisp
+++ b/random/fdist.lisp
@@ -1,6 +1,6 @@
 ;; Fdist distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:08:13EST fdist.lisp>
+;; Time-stamp: <2009-05-24 20:24:40EDT fdist.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun fdist (generator nu1 nu2)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'fdist)) &key nu1 nu2)
   "gsl_ran_fdist"
   (((mpointer generator) :pointer) (nu1 :double) (nu2 :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the F-distribution with degrees of freedom nu1
@@ -67,7 +69,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (fdist rng 1.0d0 2.0d0)))
+	    (sample rng 'fdist :nu1 1.0d0 :nu2 2.0d0)))
   (fdist-pdf 1.2d0 1.0d0 2.0d0)
   (fdist-P 1.2d0 1.0d0 2.0d0)
   (fdist-Q 1.2d0 1.0d0 2.0d0)
diff --git a/random/flat.lisp b/random/flat.lisp
index 37c3d9acfee4fb8555f81557dc0449623bfb35d7..181715c5313504a669763a6dfd2a78a3d16a987e 100644
--- a/random/flat.lisp
+++ b/random/flat.lisp
@@ -1,15 +1,17 @@
 ;; Flat distribution
 ;; Liam Healy, Oct  7 2006
-;; Time-stamp: <2009-02-16 10:08:17EST flat.lisp>
+;; Time-stamp: <2009-05-24 22:58:39EDT flat.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun flat (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'flat)) &key a b)
   "gsl_ran_flat"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the flat (uniform)
@@ -58,7 +60,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (flat rng 1.0d0 2.0d0)))
+	    (sample rng 'flat :a 1.0d0 :b 2.0d0)))
   (flat-pdf 1.2d0 1.0d0 2.0d0)
   (flat-P 1.2d0 1.0d0 2.0d0)
   (flat-Q 1.2d0 1.0d0 2.0d0)
diff --git a/random/gamma.lisp b/random/gamma.lisp
index ed76b49a302fe71617f26e3c90d0ee35a4c2d1d5..d4b77992d90f5e9314c2f7457d6d96764ad92f3b 100644
--- a/random/gamma.lisp
+++ b/random/gamma.lisp
@@ -1,6 +1,6 @@
 ;; Gamma distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:08:11EST gamma.lisp>
+;; Time-stamp: <2009-05-24 20:16:57EDT gamma.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,10 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun gamma-rd (generator a b)
-  ;; Named #'gamma-rd to avoid confusion with the special function #'gamma.
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gamma)) &key a b)
   "gsl_ran_gamma"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the gamma distribution.
@@ -21,9 +22,11 @@
    is known as the Erlang distribution.  The variates are computed using
    the algorithms from Knuth (vol 2).")
 
-(defmfun gamma-mt (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gamma-mt)) &key a b)
   "gsl_ran_gamma_mt"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A gamma variate using the Marsaglia-Tsang fast gamma method.")
@@ -69,11 +72,11 @@
  (let ((rng (make-random-number-generator +mt19937+ 0)))
    (loop for i from 0 to 10
 	 collect
-	 (gamma-rd rng 1.0d0 2.0d0)))
+	 (sample rng 'gamma :a 1.0d0 :b 2.0d0)))
  (let ((rng (make-random-number-generator +mt19937+ 0)))
    (loop for i from 0 to 10
 	 collect
-	 (gamma-mt rng 1.0d0 2.0d0)))
+	 (sample rng 'gamma-mt :a 1.0d0 :b 2.0d0)))
  (gamma-pdf 0.1d0 1.0d0 2.0d0)
  (gamma-P 0.1d0 1.0d0 2.0d0)
  (gamma-Q 0.1d0 1.0d0 2.0d0)
diff --git a/random/gaussian-bivariate.lisp b/random/gaussian-bivariate.lisp
index c1000cd68ffd7baf5f70d84fde4e85ba1fecd91a..5aa075583ab04a0020c9da254472ca4ba747d747 100644
--- a/random/gaussian-bivariate.lisp
+++ b/random/gaussian-bivariate.lisp
@@ -1,16 +1,19 @@
 ;; Gaussian bivariate distribution
 ;; Liam Healy, Sat Sep  2 2006 - 16:32
-;; Time-stamp: <2009-04-26 23:06:34EDT gaussian-bivariate.lisp>
+;; Time-stamp: <2009-05-24 22:54:07EDT gaussian-bivariate.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun bivariate-gaussian (generator sigma-x sigma-y rho)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gaussian-bivariate))
+     &key sigma-x sigma-y rho)
   "gsl_ran_bivariate_gaussian"
   (((mpointer generator) :pointer) (sigma-x :double) (sigma-y :double) (rho :double)
    (x (:pointer :double)) (y (:pointer :double)))
+  :definition :method
   :c-return :void
   :documentation			; FDL
   "Generate a pair of correlated Gaussian variates, with
@@ -39,6 +42,8 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (bivariate-gaussian rng 1.0d0 0.75d0 0.25d0)))
+	    (sample
+	     rng 'gaussian-bivariate
+	     :sigma-x 1.0d0 :sigma-y 0.75d0 :rho 0.25d0)))
   (bivariate-gaussian-pdf 0.25d0 0.5d0 0.25d0
 			   0.4d0 0.2d0))
diff --git a/random/gaussian-tail.lisp b/random/gaussian-tail.lisp
index 596852b8404306213fcb121977bd8bebc0625eb9..e932b26c0c3ec79264520bab41419cfea28d3722 100644
--- a/random/gaussian-tail.lisp
+++ b/random/gaussian-tail.lisp
@@ -1,15 +1,18 @@
 ;; Gaussian tail distribution
 ;; Liam Healy, Mon Aug 21 2006 - 21:52
-;; Time-stamp: <2009-02-16 10:08:12EST gaussian-tail.lisp>
+;; Time-stamp: <2009-05-24 16:53:42EDT gaussian-tail.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun gaussian-tail (generator a sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gaussian-tail))
+     &key a sigma)
   "gsl_ran_gaussian_tail"
   (((mpointer generator) :pointer) (a :double) (sigma :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "Random variates from the upper tail of a Gaussian
@@ -32,9 +35,12 @@
   for a Gaussian tail distribution with standard deviation sigma and
   lower limit a, using the formula given for gaussian-tail.")
 
-(defmfun ugaussian-tail (generator a)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'ugaussian-tail))
+     &key a)
   "gsl_ran_ugaussian_tail"
   (((mpointer generator) :pointer) (a :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "Equivalent to gaussian-tail with sigma=1.")
@@ -51,10 +57,10 @@
  (let ((rng (make-random-number-generator +mt19937+ 0)))
    (loop for i from 0 to 10
 	 collect
-	 (gaussian-tail rng 50.0d0 10.0d0)))
+	 (sample rng 'gaussian-tail :a 50.0d0 :sigma 10.0d0)))
  (gaussian-tail-pdf 52.0d0 50.0d0 10.0d0)
  (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (ugaussian-tail rng 5.0d0)))
+	   (sample rng 'ugaussian-tail :a 5.0d0)))
  (ugaussian-tail-pdf 5.2d0 5.0d0))
diff --git a/random/gaussian.lisp b/random/gaussian.lisp
index 314f78f84470056b63a1d05c2b3258d891fa11c4..4e55d90ab245be779d89cbd478906ffa0650d5e2 100644
--- a/random/gaussian.lisp
+++ b/random/gaussian.lisp
@@ -1,6 +1,6 @@
 ;; Gaussian distribution
 ;; Liam Healy, Sun Jul 16 2006 - 22:09
-;; Time-stamp: <2009-02-16 10:08:13EST gaussian.lisp>
+;; Time-stamp: <2009-05-24 16:42:13EDT gaussian.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,12 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun gaussian (generator sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gaussian))
+     &key sigma)
   "gsl_ran_gaussian"
   (((mpointer generator) :pointer) (sigma :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A Gaussian random variate, with mean zero and
@@ -30,25 +33,33 @@
   "Compute the probability density p(x) at x
    for a Gaussian distribution with standard deviation sigma.")
 
-(defmfun gaussian-ziggurat (generator sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gaussian-ziggurat))
+     &key sigma)
   "gsl_ran_gaussian_ziggurat"
   (((mpointer generator) :pointer) (sigma :double))
+  :definition :method
   :c-return :double
   :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.")
 
-(defmfun gaussian-ratio-method (generator sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gaussian-ratio-method))
+     &key sigma)
   "gsl_ran_gaussian_ratio_method"
   (((mpointer generator) :pointer) (sigma :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "Compute a Gaussian random variate using the Kinderman-Monahan-Leva
    ratio method.")
 
-(defmfun ugaussian (generator)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'ugaussian)) &key)
   "gsl_ran_ugaussian" (((mpointer generator) :pointer))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "Compute results for the unit Gaussian distribution,
@@ -63,9 +74,12 @@
    equivalent to the #'gaussian-pdf with a standard deviation of one,
    sigma = 1.")
 
-(defmfun ugaussian-ratio-method (generator)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'ugaussian-ratio-method))
+     &key)
   "gsl_ran_ugaussian_ratio_method"
   (((mpointer generator) :pointer))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "Compute results for the unit Gaussian distribution,
@@ -134,12 +148,12 @@
  (let ((rng (make-random-number-generator +mt19937+ 0)))
    (loop for i from 0 to 10
 	 collect
-	 (gaussian rng 10.0d0)))
+	 (sample rng 'gaussian :sigma 10.0d0)))
  (gaussian-pdf 0.0d0 10.0d0)
  (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (gaussian-ziggurat rng 10.0d0)))
+	   (sample rng 'gaussian-ziggurat :sigma 10.0d0)))
  ;; Given in examples in GSL documentation
  (ugaussian-p 2.0d0)
  (ugaussian-q 2.0d0)
diff --git a/random/generators.lisp b/random/generators.lisp
index 74cd95ac407c2afecaddae76857f2f1cb400fc88..f531f39122b43594dd22b30fe8fe64ce40b918a4 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: <2009-05-07 22:42:30EDT generators.lisp>
+;; Time-stamp: <2009-05-24 22:50:40EDT generators.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -56,8 +56,14 @@
    values of min and max can determined using the auxiliary
    functions #'rng-max and #'rng-min.")
 
-(defmfun uniform (generator)
-  "gsl_rng_uniform" (((mpointer generator) :pointer))
+(export 'sample)
+(defgeneric sample (source distribution &key #+sbcl &allow-other-keys)
+  (:documentation
+   "Sample from the probability distribution."))
+
+(defmfun sample ((source random-number-generator) (type (eql 'uniform)) &key)
+  "gsl_rng_uniform" (((mpointer source) :pointer))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A double precision floating point number uniformly
@@ -69,19 +75,23 @@
    (the maximum number of bits that can be portably represented in a single
    :ulong.")
 
-(defmfun uniform>0 (generator)
-  "gsl_rng_uniform_pos" (((mpointer generator) :pointer))
+(defmfun sample ((source random-number-generator) (type (eql 'uniform>0)) &key)
+  "gsl_rng_uniform_pos" (((mpointer source) :pointer))
+  :definition :method
   :c-return :double
   :documentation			; FDL
-  "Return a positive double precision floating point number
-   uniformly distributed in the range (0,1), excluding both 0.0 and 1.0.
-   The number is obtained by sampling the generator with the algorithm of
-   #'uniform until a non-zero value is obtained.  You can use
-   this function if you need to avoid a singularity at 0.0.")
-
-(defmfun uniform-fixnum (generator upperbound)
+  "Return a positive double precision floating point number uniformly
+   distributed in the range (0,1), excluding both 0.0 and 1.0.  The
+   number is obtained by sampling the generator with the algorithm for
+   type 'uniform until a non-zero value is obtained.  You can use this
+   function if you need to avoid a singularity at 0.0.")
+
+(defmfun sample
+    ((source random-number-generator) (type (eql 'uniform-fixnum))
+     &key upperbound)
   "gsl_rng_uniform_int"
-  (((mpointer generator) :pointer) (upperbound :ulong))
+  (((mpointer source) :pointer) (upperbound :ulong))
+  :definition :method
   :c-return :ulong
   :documentation			; FDL
   "Generate a random integer from 0 to upperbound-1 inclusive.
@@ -181,6 +191,6 @@
  (let ((rng (make-random-number-generator +mt19937+ 0)))
    (loop for i from 0 to 10
 	 collect
-	 (uniform-fixnum rng 1000)))
+	 (sample rng 'uniform-fixnum :upperbound 1000)))
  (let ((rng (make-random-number-generator *cmrg* 0)))
-   (loop for i from 0 to 10 collect (uniform rng))))
+   (loop for i from 0 to 10 collect (sample rng 'uniform))))
diff --git a/random/geometric.lisp b/random/geometric.lisp
index 0cf3e9f380326c38fdaa4180d627ae5004b9d120..adeb277d92981c043adeef48b22867330cbcf033 100644
--- a/random/geometric.lisp
+++ b/random/geometric.lisp
@@ -1,15 +1,18 @@
 ;; Geometric distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2009-02-16 10:08:14EST geometric.lisp>
+;; Time-stamp: <2009-05-24 22:28:31EDT geometric.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun geometric (generator p)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'geometric))
+     &key probability)
   "gsl_ran_geometric"
-  (((mpointer generator) :pointer) (p :double))
+  (((mpointer generator) :pointer) (probability :double))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "A random integer from the geometric distribution,
@@ -47,7 +50,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (geometric rng 0.4d0)))
+	   (sample rng 'geometric :probability 0.4d0)))
   (geometric-pdf 2 0.4d0)
   (geometric-P 2 0.4d0)
   (geometric-Q 2 0.4d0))
diff --git a/random/gumbel1.lisp b/random/gumbel1.lisp
index 54ca833888825c0dfbd6426c5d92cb2ddf696f96..6fd442e308fd4a3e7f6741ed83747bcd469b0d76 100644
--- a/random/gumbel1.lisp
+++ b/random/gumbel1.lisp
@@ -1,6 +1,6 @@
 ;; The Gumbel type 1 random number distribution
 ;; Liam Healy, Sun Oct 29 2006
-;; Time-stamp: <2009-02-16 10:10:30EST gumbel1.lisp>
+;; Time-stamp: <2009-05-24 21:45:10EDT gumbel1.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun gumbel1 (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gumbel1)) &key a b)
   "gsl_ran_gumbel1"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Type-1 Gumbel
@@ -62,7 +64,7 @@
 (save-test gumbel1
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
-	    collect (gumbel1 rng 1.0d0 2.0d0)))
+	    collect (sample rng 'gumbel1 :a 1.0d0 :b 2.0d0)))
   (gumbel1-pdf 0.1d0 1.0d0 2.0d0)
   (gumbel1-P 0.1d0 1.0d0 2.0d0)
   (gumbel1-Q 0.1d0 1.0d0 2.0d0)
diff --git a/random/gumbel2.lisp b/random/gumbel2.lisp
index 3b10191eb60f2dfe828e04702220261318230ce1..49de78544ac8966d5f08b362ee1f1866c55117fb 100644
--- a/random/gumbel2.lisp
+++ b/random/gumbel2.lisp
@@ -1,6 +1,6 @@
 ;; The Gumbel type 2 random number distribution
 ;; Liam Healy, Sun Oct 29 2006
-;; Time-stamp: <2009-02-16 10:08:12EST gumbel2.lisp>
+;; Time-stamp: <2009-05-24 21:47:03EDT gumbel2.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun gumbel2 (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'gumbel2)) &key a b)
   "gsl_ran_gumbel2"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Type-2 Gumbel
@@ -61,7 +63,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (gumbel2 rng 1.0d0 2.0d0)))
+	    (sample rng 'gumbel2 :a 1.0d0 :b 2.0d0)))
   (gumbel2-pdf 5.0d0 1.0d0 2.0d0)
   (gumbel2-P 10.0d0 1.0d0 2.0d0)
   (gumbel2-Q 10.0d0 1.0d0 2.0d0)
diff --git a/random/hypergeometric.lisp b/random/hypergeometric.lisp
index 4d04edc529af6c0ba480dbc9ec24f8982eb9f1ea..95b20f577b78a0391cf77668ed1c0a6e8fe58257 100644
--- a/random/hypergeometric.lisp
+++ b/random/hypergeometric.lisp
@@ -1,6 +1,6 @@
 ;; Hypergeometric distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2009-02-16 10:08:12EST hypergeometric.lisp>
+;; Time-stamp: <2009-05-24 22:31:19EDT hypergeometric.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,12 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun hypergeometric (generator n1 n2 tt)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'hypergeometric))
+     &key n1 n2 tt)
   "gsl_ran_hypergeometric"
   (((mpointer generator) :pointer) (n1 :uint) (n2 :uint)(tt :uint))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "A random integer from the hypergeometric
@@ -41,7 +44,7 @@
   "The cumulative distribution functions P(k) for the
    hypergeometric distribution with parameters n1, n2 and tt.")
 
-(defmfun hypergeometric-Q  (k n1 n2 tt)
+(defmfun hypergeometric-Q (k n1 n2 tt)
   "gsl_cdf_hypergeometric_Q"  ((k :uint) (n1 :uint) (n2 :uint)(tt :uint))
   :c-return :double
   :documentation			; FDL
@@ -53,9 +56,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (hypergeometric rng 3 6 3)))
+	   (sample rng 'hypergeometric :n1 3 :n2 6 :tt 3)))
   (hypergeometric-pdf 0 2 6 3)
   (hypergeometric-P 1 2 6 3)
   (hypergeometric-Q 1 2 6 3))
-
-
diff --git a/random/landau.lisp b/random/landau.lisp
index 40959af7d589b915a8a2ee72f907285d702d3dff..39ed973ec159059376e461b5d125d6ef746377c6 100644
--- a/random/landau.lisp
+++ b/random/landau.lisp
@@ -1,15 +1,17 @@
 ;; Landau distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:05:57EST landau.lisp>
+;; Time-stamp: <2009-05-24 20:08:47EDT landau.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun landau (generator)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'landau)) &key)
   "gsl_ran_landau"
   (((mpointer generator) :pointer))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Landau distribution.  The
@@ -32,5 +34,5 @@
 (save-test landau
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
-	    collect (landau rng)))
+	    collect (sample rng 'landau)))
   (landau-pdf 0.25d0))
diff --git a/random/laplace.lisp b/random/laplace.lisp
index 31628e8e4ad3b656584e20a4227a4ca60827e739..167321fc82124b7287bf269f7e74eae48480a214 100644
--- a/random/laplace.lisp
+++ b/random/laplace.lisp
@@ -1,6 +1,6 @@
 ;; Exponential distribution
 ;; Liam Healy, Sun Sep 17 2006
-;; Time-stamp: <2009-02-16 10:08:15EST laplace.lisp>
+;; Time-stamp: <2009-05-24 19:32:56EDT laplace.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun laplace (generator a)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'laplace)) &key a)
   "gsl_ran_laplace"
   (((mpointer generator) :pointer) (a :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Laplace distribution with width a.
@@ -59,7 +61,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (laplace rng 10.0d0)))
+	    (sample rng 'laplace :a 10.0d0)))
   (laplace-pdf 0.0d0 10.0d0)
   (laplace-p 1.0d0 2.0d0)
   (laplace-q 1.0d0 2.0d0)
diff --git a/random/levy.lisp b/random/levy.lisp
index d680fd6f13df2015eb65bad34189ab2d541abfdb..adacff521319ad3b3267a1bdb82014579ac23863 100644
--- a/random/levy.lisp
+++ b/random/levy.lisp
@@ -1,15 +1,17 @@
 ;; Levy distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:08:11EST levy.lisp>
+;; Time-stamp: <2009-05-24 20:10:32EDT levy.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun levy (generator c alpha)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'levy)) &key c alpha)
   "gsl_ran_levy"
   (((mpointer generator) :pointer) (c :double) (alpha :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Levy symmetric stable
@@ -23,9 +25,12 @@
    For \alpha < 1 the tails of the distribution become extremely wide.
    The algorithm only works for 0 < alpha <= 2.")
 
-(defmfun levy-skew (generator c alpha beta)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'levy-skew))
+     &key c alpha beta)
   "gsl_ran_levy_skew"
   (((mpointer generator) :pointer) (c :double) (alpha :double) (beta :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Levy skew stable
@@ -49,8 +54,8 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (levy rng 1.0d0 2.0d0)))
+	    (sample rng 'levy :c 1.0d0 :alpha 2.0d0)))
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (levy-skew rng 1.0d0 2.0d0 1.0d0))))
+	    (sample rng 'levy-skew :c 1.0d0 :alpha 2.0d0 :beta 1.0d0))))
diff --git a/random/logarithmic.lisp b/random/logarithmic.lisp
index b3e7f6571ed5638c33d77f1f8bf5f8d3bfb3cca3..07928b07afa9e6ef6988f5c141d3c29d143075f7 100644
--- a/random/logarithmic.lisp
+++ b/random/logarithmic.lisp
@@ -1,15 +1,18 @@
 ;; Logarithmic distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2009-02-16 10:08:16EST logarithmic.lisp>
+;; Time-stamp: <2009-05-24 22:49:00EDT logarithmic.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun logarithmic (generator p)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'logarithmic))
+     &key probability)
   "gsl_ran_logarithmic"
-  (((mpointer generator) :pointer) (p :double))
+  (((mpointer generator) :pointer) (probability :double))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "A random integer from the logarithmic distribution.
@@ -30,5 +33,5 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (logarithmic rng 0.9d0)))
+	   (sample rng 'logarithmic :probability 0.9d0)))
   (logarithmic-pdf 2 0.4d0))
diff --git a/random/logistic.lisp b/random/logistic.lisp
index 6ef3bf92c1375322060ce7497ffd85642056406f..5d7fe42e64a233f451e06ea4e3131e74415d4b4c 100644
--- a/random/logistic.lisp
+++ b/random/logistic.lisp
@@ -1,6 +1,6 @@
 ;; Logistic distribution
 ;; Liam Healy, Sat Oct  7 2006 - 16:13
-;; Time-stamp: <2009-02-16 10:08:10EST logistic.lisp>
+;; Time-stamp: <2009-05-24 21:36:39EDT logistic.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun logistic (generator a)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'logistic)) &key a)
   "gsl_ran_logistic"
   (((mpointer generator) :pointer) (a :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the logistic distribution.  The distribution function is
@@ -58,7 +60,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (logistic rng 10.0d0)))
+	    (sample rng 'logistic :a 10.0d0)))
   (logistic-pdf 0.5d0 1.0d0)
   (logistic-P 0.5d0 1.0d0)
   (logistic-Q 0.5d0 1.0d0)
diff --git a/random/lognormal.lisp b/random/lognormal.lisp
index 23e052946b32cb2e8c187b5562c86483c49c76b8..18581360663e1e3cb648823984d6af542926e388 100644
--- a/random/lognormal.lisp
+++ b/random/lognormal.lisp
@@ -1,13 +1,16 @@
 ;; Lognormal distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:08:15EST lognormal.lisp>
+;; Time-stamp: <2009-05-24 20:22:32EDT lognormal.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
-(defmfun lognormal (generator zeta sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'lognormal))
+     &key zeta sigma)
   "gsl_ran_lognormal"
   (((mpointer generator) :pointer) (zeta :double) (sigma :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the lognormal distribution.
@@ -59,7 +62,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (lognormal rng 1.0d0 2.0d0)))
+	    (sample rng 'lognormal :zeta 1.0d0 :sigma 2.0d0)))
   (lognormal-pdf 1.2d0 1.0d0 2.0d0)
   (lognormal-P 1.2d0 1.0d0 2.0d0)
   (lognormal-Q 1.2d0 1.0d0 2.0d0)
diff --git a/random/multinomial.lisp b/random/multinomial.lisp
index b5edc8f0dd7a39dc58184e45e1b42684ff432818..5645d1fe76ecd6945958af408ea05ad23f79200e 100644
--- a/random/multinomial.lisp
+++ b/random/multinomial.lisp
@@ -1,28 +1,33 @@
 ;; Multinomial distribution
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2009-02-16 10:10:33EST multinomial.lisp>
+;; Time-stamp: <2009-05-24 23:00:01EDT multinomial.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun multinomial (generator sum p n)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'multinomial))
+     &key sum probabilities
+     (n (make-marray
+	 '(signed-byte 32) :dimensions (dim0 probabilities))))
   "gsl_ran_multinomial"
   (((mpointer generator) :pointer) 
-   ((dim0 p) sizet)
+   ((dim0 probabilities) sizet)
    (sum sizet)
-   ((c-pointer p) :pointer)
+   ((c-pointer probabilities) :pointer)
    ;; technically, n should be a uint array, but integers work
    ((c-pointer n) :pointer))
+  :definition :method
   :inputs (p)
   :outputs (n)
   :return (n)
   :c-return :void
   :documentation			; FDL
-  "Returns an array n of K random variates from a 
+  "Returns an array n of (dim0 probabilities) random variates from a 
    multinomial distribution.  The sum of the array n is specified
-   by sum=N.  The distribution function is
+   by sum.  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 (n_1, n_2, ..., n_K) are nonnegative integers with 
@@ -57,13 +62,11 @@
 ;;; Examples and unit test
 (save-test multinomial
  (let ((rng (make-random-number-generator +mt19937+ 0))
-	(p #m(0.1d0 0.2d0 0.3d0 0.4d0))
-	(n (make-marray '(signed-byte 32) :dimensions 4)))
-   (multinomial rng 8 p n)
-   (cl-array n))
+       (p #m(0.1d0 0.2d0 0.3d0 0.4d0)))
+   (cl-array (sample rng 'multinomial :sum 8 :probabilities p)))
  (let ((p #m(0.1d0 0.2d0 0.3d0 0.4d0))
-	(n #31m(5 0 1 2)))
+       (n #31m(5 0 1 2)))
    (multinomial-pdf p N))
  (let ((p #m(0.1d0 0.2d0 0.3d0 0.4d0))
-	(n #31m(5 0 1 2)))
+       (n #31m(5 0 1 2)))
    (multinomial-log-pdf p n)))
diff --git a/random/negative-binomial.lisp b/random/negative-binomial.lisp
index 8b4a709f553f87faf5fc92bda875ab86a11822bf..273e90ce482a3a37659cff46f74d7fd29b197c7f 100644
--- a/random/negative-binomial.lisp
+++ b/random/negative-binomial.lisp
@@ -1,6 +1,6 @@
 ;; Negative binomial and Pascal distributions
 ;; Liam Healy, Sat Nov 25 2006 - 16:00
-;; Time-stamp: <2009-02-16 10:05:56EST negative-binomial.lisp>
+;; Time-stamp: <2009-05-24 22:27:29EDT negative-binomial.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -11,15 +11,19 @@
 ;;;; Negative binomial
 ;;;;****************************************************************************
 
-(defmfun negative-binomial (generator p n)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'negative-binomial))
+     &key probability n)
   "gsl_ran_negative_binomial"
-  (((mpointer generator) :pointer) (p :double) (n :double))
+  (((mpointer generator) :pointer) (probability :double) (n :double))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "A random integer from the negative binomial
    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,
+   in independent trials with probability of success.  The
+   probability distribution for negative binomial variates is
+   given by probability (p):
    p(k) = {\Gamma(n + k) \over \Gamma(k+1) \Gamma(n) } p^n (1-p)^k
    Note that n is not required to be an integer.")
 
@@ -51,9 +55,12 @@
 ;;;; Pascal
 ;;;;****************************************************************************
 
-(defmfun pascal (generator p n)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'pascal))
+     &key probability n)
   "gsl_ran_pascal"
-  (((mpointer generator) :pointer) (p :double) (n :uint))
+  (((mpointer generator) :pointer) (probability :double) (n :uint))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "A random integer from the Pascal distribution.  The
@@ -94,14 +101,14 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (negative-binomial rng 0.4d0 12.0d0)))
+	   (sample rng 'negative-binomial :probability 0.4d0 :n 12.0d0)))
   (negative-binomial-pdf 5 0.4d0 12.0d0)
   (negative-binomial-P 5 0.4d0 12.0d0)
   (negative-binomial-Q 5 0.4d0 12.0d0)
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (pascal rng 0.4d0 12)))
+	   (sample rng 'pascal :probability 0.4d0 :n 12)))
   (pascal-pdf 5 0.4d0 12)
   (pascal-P 5 0.4d0 12)
   (pascal-Q 5 0.4d0 12))
diff --git a/random/pareto.lisp b/random/pareto.lisp
index 505e768a10a2290be3cfb02f2ee15b4086aa7c7d..1a35cc7a32a1e7da42593ea03b9260a150c194c1 100644
--- a/random/pareto.lisp
+++ b/random/pareto.lisp
@@ -1,6 +1,6 @@
 ;; Pareto distribution
 ;; Liam Healy, Sat Oct  8 2006 - 21:23
-;; Time-stamp: <2009-02-16 10:10:34EST pareto.lisp>
+;; Time-stamp: <2009-05-24 21:38:03EDT pareto.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun pareto (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'pareto)) &key a b)
   "gsl_ran_pareto"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Pareto distribution of order a.
@@ -59,7 +61,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (pareto rng 1.0d0 2.0d0)))
+	    (sample rng 'pareto :a 1.0d0 :b 2.0d0)))
   (pareto-pdf 1.5d0 1.3d0 1.0d0)
   (pareto-P 3.5d0 1.3d0 2.0d0)
   (pareto-Q 3.5d0 1.3d0 2.0d0)
diff --git a/random/poisson.lisp b/random/poisson.lisp
index e4e8da91ffd0f858f1ca588ae0e3dc0c3e34f828..7a3a66e40158ec0030901fd93a6fc8927874cc0d 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: <2009-02-16 10:10:31EST poisson.lisp>
+;; Time-stamp: <2009-05-24 21:58:22EDT poisson.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun poisson (generator mu)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'poisson)) &key mu)
   "gsl_ran_poisson"
   (((mpointer generator) :pointer) (mu :double))
+  :definition :method
   :c-return :uint
   :documentation			; FDL
   "A random integer from the Poisson distribution with mean mu.
@@ -45,7 +47,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
      (loop for i from 0 to 10
 	   collect
-	   (poisson rng 10.0d0)))
+	   (sample rng 'poisson :mu 10.0d0)))
   (poisson-pdf 8 10.0d0)
   (poisson-P 8 10.0d0)
   (poisson-Q 8 10.0d0))
diff --git a/random/rayleigh-tail.lisp b/random/rayleigh-tail.lisp
index e8218dfbca50870050f856ec13c5cb489a8c28c0..548b10a18a19b26c9d0c4e3f41846d0cc1aabff6 100644
--- a/random/rayleigh-tail.lisp
+++ b/random/rayleigh-tail.lisp
@@ -1,15 +1,18 @@
 ;; Rayleigh tail distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:10:31EST rayleigh-tail.lisp>
+;; Time-stamp: <2009-05-24 20:07:27EDT rayleigh-tail.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun rayleigh-tail (generator a sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'rayleigh-tail))
+     &key a sigma)
   "gsl_ran_rayleigh_tail"
   (((mpointer generator) :pointer) (a :double) (sigma :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the tail of the Rayleigh
@@ -30,5 +33,5 @@
 (save-test rayleigh-tail
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
-	    collect (rayleigh-tail rng 1.0d0 10.0d0)))
+	    collect (sample rng 'rayleigh-tail :a 1.0d0 :sigma 10.0d0)))
   (rayleigh-tail-pdf 0.25d0 -2.0d0 2.0d0))
diff --git a/random/rayleigh.lisp b/random/rayleigh.lisp
index 233a1156e8f4574afbad9977e0b9fedbf11e0523..1e460c2bc82a0d298d82af102f5257a00efa2e51 100644
--- a/random/rayleigh.lisp
+++ b/random/rayleigh.lisp
@@ -1,15 +1,18 @@
 ;; Rayleigh distribution
 ;; Liam Healy, Sat Sep 30 2006
-;; Time-stamp: <2009-02-16 10:05:57EST rayleigh.lisp>
+;; Time-stamp: <2009-05-24 20:05:46EDT rayleigh.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun rayleigh (generator sigma)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'rayleigh))
+     &key sigma)
   "gsl_ran_rayleigh"
   (((mpointer generator) :pointer) (sigma :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Rayleigh distribution with
@@ -62,7 +65,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (rayleigh rng 10.0d0)))
+	    (sample rng 'rayleigh :sigma 10.0d0)))
   (rayleigh-pdf 0.5d0 1.0d0)
   (rayleigh-P 1.0d0 2.0d0)
   (rayleigh-Q 1.0d0 2.0d0)
diff --git a/random/shuffling-sampling.lisp b/random/shuffling-sampling.lisp
index 852ece5c426120796ee1020c48537c92149b6924..d9be085f0da6d5e14381682604beab36e232d8e9 100644
--- a/random/shuffling-sampling.lisp
+++ b/random/shuffling-sampling.lisp
@@ -1,16 +1,19 @@
 ;; Shuffling and sampling
 ;; Liam Healy, Sat Dec  2 2006 - 18:40
-;; Time-stamp: <2009-02-16 10:10:30EST shuffling-sampling.lisp>
+;; Time-stamp: <2009-05-24 23:13:40EDT shuffling-sampling.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; These are currently defined only for vectors.
 
-(defmfun shuffle (generator base)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'shuffle))
+     &key base)
   "gsl_ran_shuffle"
   (((mpointer generator) :pointer)
    ((c-pointer base) :pointer) ((dim0 base) sizet) ((element-size base) sizet))
+  :definition :method
   :inputs (base)
   :outputs (base)
   :c-return :void
@@ -22,34 +25,49 @@
    permutations with equal probability, assuming a perfect source of random
    numbers.")
 
-(defmfun choose-random (generator dest src)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'choose-random))
+     &key src (dest (dim0 src))
+     &aux
+     (destarr
+      (if (integerp dest)
+	  (make-marray (element-type src) :dimensions dest)
+	  dest)))
   "gsl_ran_choose"
   (((mpointer generator) :pointer)
-   ((c-pointer dest) :pointer) ((dim0 dest) sizet)
+   ((c-pointer destarr) :pointer) ((dim0 destarr) sizet)
    ((c-pointer src) :pointer) ((dim0 src) sizet) ((element-size src) sizet))
-  :inputs (dest src)
-  :outputs (dest)
+  :definition :method
+  :inputs (destarr src)
+  :outputs (destarr)
   :documentation			; FDL
-  "Fill the array dest[k] with k objects taken
-   randomly from the n elements of the array
-   src[0...n-1].  The objects are each of size size.  The
-   output of the random number generator r is used to make the
-   selection.  The algorithm ensures all possible samples are equally
-   likely, assuming a perfect source of randomness.
+  "Fill the array dest[k] with k objects taken randomly from the n
+   elements of the array src[0...n-1].  The output of the random
+   number generator r is used to make the selection.  The algorithm
+   ensures all possible samples are equally likely, assuming a perfect
+   source of randomness.
 
    The objects are sampled without replacement, thus each object can
    only appear once in dest[k].  It is required that k be less
    than or equal to n.  The objects in dest will be in the
    same relative order as those in src.  You will need to call
-   #'shuffle if you want to randomize the order.")
+   with 'shuffle if you want to randomize the order.")
 
-(defmfun random-sample (generator dest src)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'random-sample))
+     &key src (dest (dim0 src))
+     &aux
+     (destarr
+      (if (integerp dest)
+	  (make-marray (element-type src) :dimensions dest)
+	  dest)))
   "gsl_ran_sample"
   (((mpointer generator) :pointer)
-   ((c-pointer dest) :pointer) ((dim0 dest) sizet)
+   ((c-pointer destarr) :pointer) ((dim0 destarr) sizet)
    ((c-pointer src) :pointer) ((dim0 src) sizet) ((element-size src) sizet))
-  :inputs (dest src)
-  :outputs (dest)
+  :definition :method
+  :inputs (destarr src)
+  :outputs (destarr)
   :c-return :void
   :documentation
   "Like #'choose-random, but samples k items
@@ -61,16 +79,11 @@
 ;;; Examples and unit test
 (save-test shuffling-sampling
  (let ((rng (make-random-number-generator +mt19937+ 0))
-	(v1 #31m(1 2 3 4 5 6 7 8)))
-   (shuffle rng v1)
-   (cl-array v1))
+       (v1 #31m(1 2 3 4 5 6 7 8)))
+   (cl-array (sample rng 'shuffle :base v1)))
  (let ((rng (make-random-number-generator +mt19937+ 0))
-	(v1 #31m(1 2 3 4 5 6 7 8))
-	(v2 (make-marray '(signed-byte 32) :dimensions 4)))
-   (choose-random rng v2 v1)
-   (cl-array v2))
+       (v1 #31m(1 2 3 4 5 6 7 8)))
+   (cl-array (sample rng 'choose-random :src v1 :dest 4)))
  (let ((rng (make-random-number-generator +mt19937+ 0))
-	(v1 #31m(1 2 3 4 5 6 7 8))
-	(v2 (make-marray '(signed-byte 32) :dimensions 10)))
-   (random-sample rng v2 v1)
-   (cl-array v2)))
+       (v1 #31m(1 2 3 4 5 6 7 8)))
+   (cl-array (sample rng 'random-sample :src v1 :dest 10))))
diff --git a/random/spherical-vector.lisp b/random/spherical-vector.lisp
index d563ba6b19172d889ec4655891dc1514dd4cc3e0..5d242c87fd5e72ffc8321af643550282489bf616 100644
--- a/random/spherical-vector.lisp
+++ b/random/spherical-vector.lisp
@@ -1,6 +1,6 @@
 ;; Spherical Vector distribution
 ;; Liam Healy, Sun Oct  22 2006
-;; Time-stamp: <2009-04-26 23:08:37EDT spherical-vector.lisp>
+;; Time-stamp: <2009-05-24 21:41:40EDT spherical-vector.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -9,30 +9,36 @@
 
 ;;; No test for #'direction-Nd yet.
 
-(defmfun direction-2d (generator)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'direction-2d)) &key)
   "gsl_ran_dir_2d"
   (((mpointer generator) :pointer)
    (x (:pointer :double)) (y (:pointer :double)))
+  :definition :method
   :c-return :void
   :documentation			; FDL
   "A random direction vector v = (x,y) in
    two dimensions.  The vector is normalized such that
    |v|^2 = x^2 + y^2 = 1.")
 
-(defmfun direction-2d-trig-method (generator)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'direction-2d-trig-method)) &key)
   "gsl_ran_dir_2d_trig_method"
   (((mpointer generator) :pointer)
    (x (:pointer :double)) (y (:pointer :double)))
+  :definition :method
   :c-return :void
   :documentation			; FDL
   "A random direction vector v = (x,y) in
    two dimensions.  The vector is normalized such that
    |v|^2 = x^2 + y^2 = 1.  Uses trigonometric functions.")
 
-(defmfun direction-3d (generator)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'direction-3d)) &key)
   "gsl_ran_dir_3d"
   (((mpointer generator) :pointer)
    (x (:pointer :double)) (y (:pointer :double)) (z (:pointer :double)))
+  :definition :method
   :c-return :void
   :documentation			; FDL
   "A random direction vector v =
@@ -43,10 +49,12 @@
   projected along any axis is actually uniform (this is only true for 3
   dimensions).")
 
-(defmfun direction-Nd (generator vector)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'direction-Nd)) &key vector)
   "gsl_ran_dir_nd"
   (((mpointer generator) :pointer) ((dim0 vector) sizet)
    ((c-pointer vector) :pointer))
+  :definition :method
   :c-return :void
   :outputs (vector)
   :return (vector)
@@ -65,12 +73,12 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 4
 	    append
-	    (multiple-value-list (direction-2d rng))))
+	    (multiple-value-list (sample rng 'direction-2d))))
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 4
 	    append
-	    (multiple-value-list (direction-2d-trig-method rng))))
+	    (multiple-value-list (sample rng 'direction-2d-trig-method))))
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 2
 	    append
-	    (multiple-value-list (direction-3d rng)))))
+	    (multiple-value-list (sample rng 'direction-3d)))))
diff --git a/random/tdist.lisp b/random/tdist.lisp
index 4e30266c442ff46d4d31a0f2593ac375f5b8f6f1..91f49a761b2076a0c9ed5c9959a442163afb197c 100644
--- a/random/tdist.lisp
+++ b/random/tdist.lisp
@@ -1,15 +1,17 @@
 ;; Tdist distribution
 ;; Liam Healy, Sat Oct  7 2006 - 16:13
-;; Time-stamp: <2009-02-16 10:05:56EST tdist.lisp>
+;; Time-stamp: <2009-05-24 21:33:28EDT tdist.lisp>
 ;; $Id$
 
 (in-package :gsl)
 
 ;;; /usr/include/gsl/gsl_randist.h
 
-(defmfun tdist (generator nu)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'tdist)) &key nu)
   "gsl_ran_tdist"
   (((mpointer generator) :pointer) (nu :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Student t-distribution.  The
@@ -59,7 +61,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (tdist rng 10.0d0)))
+	    (sample rng 'tdist :nu 10.0d0)))
   (tdist-pdf 0.5d0 1.0d0)
   (tdist-P 0.5d0 1.0d0)
   (tdist-Q 0.5d0 1.0d0)
diff --git a/random/weibull.lisp b/random/weibull.lisp
index 7108e61c10a0fde6f4528537876512b246f792e8..549dd937520af50ebf6582ca8a50614a445a1130 100644
--- a/random/weibull.lisp
+++ b/random/weibull.lisp
@@ -1,6 +1,6 @@
 ;; Weibull distribution
 ;; Liam Healy, Sun Oct 22 2006
-;; Time-stamp: <2009-02-16 10:08:16EST weibull.lisp>
+;; Time-stamp: <2009-05-24 21:43:44EDT weibull.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -8,9 +8,11 @@
 ;;; /usr/include/gsl/gsl_randist.h
 ;;; /usr/include/gsl/gsl_cdf.h
 
-(defmfun weibull (generator a b)
+(defmfun sample
+    ((generator random-number-generator) (type (eql 'weibull)) &key a b)
   "gsl_ran_weibull"
   (((mpointer generator) :pointer) (a :double) (b :double))
+  :definition :method
   :c-return :double
   :documentation			; FDL
   "A random variate from the Weibull distribution.  The distribution function is
@@ -58,7 +60,7 @@
   (let ((rng (make-random-number-generator +mt19937+ 0)))
       (loop for i from 0 to 10
 	    collect
-	    (weibull rng 1.0d0 2.0d0)))
+	    (sample rng 'weibull :a 1.0d0 :b 2.0d0)))
   (weibull-pdf 1.5d0 1.3d0 1.0d0)
   (weibull-P 3.5d0 1.3d0 2.0d0)
   (weibull-Q 3.5d0 1.3d0 2.0d0)
diff --git a/solve-minimize-fit/linear-least-squares.lisp b/solve-minimize-fit/linear-least-squares.lisp
index 6d1c4d6fff0b370029ee7ba13d892ba008fdc50c..91435af23e76a01ba52ac1d16dc3cbe97f642dc5 100644
--- a/solve-minimize-fit/linear-least-squares.lisp
+++ b/solve-minimize-fit/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: <2009-04-26 23:36:23EDT linear-least-squares.lisp>
+;; Time-stamp: <2009-05-24 22:56:00EDT linear-least-squares.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -326,7 +326,7 @@
 	  for y0 = (exp xd)
 	  for sigma = (* 0.1d0 y0)
 	  collect
-	  (list xd (+ y0 (gaussian rng sigma)) sigma))))
+	  (list xd (+ y0 (sample rng 'gaussian :sigma sigma)) sigma))))
 
 (defun linear-least-squares-multivariate-example (data &optional (print-details t))
   "Second example in Section 36.5 of the GSL manual.  Returns the
diff --git a/solve-minimize-fit/nonlinear-least-squares.lisp b/solve-minimize-fit/nonlinear-least-squares.lisp
index bec7440cf5fdd5e03b87ba547a379848fdad8aee..6e41da22dd44bd503192c2185af9938ae18c8309 100644
--- a/solve-minimize-fit/nonlinear-least-squares.lisp
+++ b/solve-minimize-fit/nonlinear-least-squares.lisp
@@ -1,6 +1,6 @@
 ;; Nonlinear least squares fitting.
 ;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
-;; Time-stamp: <2009-04-04 22:30:04EDT nonlinear-least-squares.lisp>
+;; Time-stamp: <2009-05-24 22:56:35EDT nonlinear-least-squares.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -294,7 +294,8 @@
 	 (rng (make-random-number-generator +mt19937+ 0)))
      (dotimes (i number-of-observations arr)
        (setf (maref arr i)
-	     (+ 1 (* 5 (exp (* -1/10 i))) (gaussian rng 0.1d0)))))
+	     (+ 1 (* 5 (exp (* -1/10 i)))
+		(sample rng 'gaussian :sigma 0.1d0)))))
    :sigma
    (make-marray
     'double-float :dimensions number-of-observations :initial-element 0.1d0)))
diff --git a/solve-minimize-fit/simulated-annealing.lisp b/solve-minimize-fit/simulated-annealing.lisp
index 1b9d724503831efd53c2be35cf91a8f8f878f86c..8528831c24a8dd6378b0107d1f3799b94848603d 100644
--- a/solve-minimize-fit/simulated-annealing.lisp
+++ b/solve-minimize-fit/simulated-annealing.lisp
@@ -1,6 +1,6 @@
 ;; Simulated Annealing
 ;; Liam Healy Sun Feb 11 2007 - 17:23
-;; Time-stamp: <2009-05-23 21:49:10EDT simulated-annealing.lisp>
+;; Time-stamp: <2009-05-24 23:17:07EDT simulated-annealing.lisp>
 ;; $Id$
 
 (in-package :gsl)
@@ -75,7 +75,7 @@
 
 ;;; The user-step-function should take a rng pointer, a state pointer,
 ;;; and a double-float, and return nothing.  The rng pointer should be
-;;; used in call to one of the distribution functions like #'uniform.
+;;; used in call to one of the distribution functions like 'uniform.
 ;;; typedef void (*gsl_siman_step_t) (const gsl_rng *r, void *xp, double step_size);
 (cffi:defcallback sa-step-function :void
     ((rng :pointer) (state :pointer) (step-size :double))
@@ -182,7 +182,7 @@
 (defun trivial-example-step (generator state step-size)
   (declare (type double-float step-size) (optimize (speed 3) (safety 1)))
   (symbol-macrolet ((x (maref state 0)))
-    (let ((rand (uniform generator)))
+    (let ((rand (sample generator 'uniform)))
       (declare (type double-float rand))
       (setf x (+  (the double-float x) (- (* 2.0d0 rand step-size) step-size))))))
 
@@ -218,7 +218,7 @@
 (defun simulated-annealing-test (initial-value)
   (simulated-annealing
    (list initial-value)
-   (list 200 1000 1.0d0 1.0d0 0.008d0 1.003d0 2.0d-6) ; parameters
+   200 1000 1.0d0 1.0d0 0.008d0 1.003d0 2.0d-6 ; parameters
    (make-random-number-generator +mt19937+ 0)
    (lambda (&optional initial)
      (make-marray 'double-float :dimensions 1 :initial-contents initial))
diff --git a/tests/bernoulli.lisp b/tests/bernoulli.lisp
index 85d5972735ff2fc0bff9d9e5b6663159c224844f..3be67a655e44b3192f84ee522abab1a7373a32a4 100644
--- a/tests/bernoulli.lisp
+++ b/tests/bernoulli.lisp
@@ -3,14 +3,13 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST BERNOULLI
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 0 1 1 0 1 1 0 0 0 0 0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (BERNOULLI RNG 0.5d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.5d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (BERNOULLI-PDF 0
-                                                                          0.5d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 0 1 1 0 1 1 0 0 0 0 0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'bernoulli :p 0.5d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5d0)
+   (MULTIPLE-VALUE-LIST (BERNOULLI-PDF 0 0.5d0))))
 
diff --git a/tests/beta.lisp b/tests/beta.lisp
index 70e95789e21c3f02dc9a924295a4c95ea0b6b6df..fc983ee8bce795856c2170b4aac1d1f9aa28822b 100644
--- a/tests/beta.lisp
+++ b/tests/beta.lisp
@@ -3,31 +3,31 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST BETA
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 0.3935474359990073d0 0.7063621551518341d0
-                               0.044515648447265056d0 0.09286083229785232d0
-                               0.210544366728104d0 0.010114317425185686d0
-                               0.4595767375719009d0 0.1515157002550483d0
-                               0.1731331145031117d0 0.4270743075655188d0
-                               0.3353314142479658d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (BETA-RD RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 1.8000000000000016d0)
-                        (MULTIPLE-VALUE-LIST (BETA-PDF 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.19000000000000017d0)
-                        (MULTIPLE-VALUE-LIST (BETA-P 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.8099999999999998d0)
-                        (MULTIPLE-VALUE-LIST (BETA-Q 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.09999999999999992d0)
-                        (MULTIPLE-VALUE-LIST (BETA-PINV 0.19d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.09999999999999988d0)
-                        (MULTIPLE-VALUE-LIST (BETA-QINV 0.81d0 1.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.3935474359990073d0 0.7063621551518341d0
+	  0.044515648447265056d0 0.09286083229785232d0
+	  0.210544366728104d0 0.010114317425185686d0
+	  0.4595767375719009d0 0.1515157002550483d0
+	  0.1731331145031117d0 0.4270743075655188d0
+	  0.3353314142479658d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'beta :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.8000000000000016d0)
+   (MULTIPLE-VALUE-LIST (BETA-PDF 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19000000000000017d0)
+   (MULTIPLE-VALUE-LIST (BETA-P 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8099999999999998d0)
+   (MULTIPLE-VALUE-LIST (BETA-Q 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.09999999999999992d0)
+   (MULTIPLE-VALUE-LIST (BETA-PINV 0.19d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.09999999999999988d0)
+   (MULTIPLE-VALUE-LIST (BETA-QINV 0.81d0 1.0d0 2.0d0))))
 
diff --git a/tests/cauchy.lisp b/tests/cauchy.lisp
index f9c778da500858d5c502033cd5716b1caefe04e4..2558f563ac5d26a70d8bac60a82cc893c0dbe2d3 100644
--- a/tests/cauchy.lisp
+++ b/tests/cauchy.lisp
@@ -3,33 +3,33 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST CAUCHY
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST -0.00811319915595434d0 5.617196410586812d0
-                               12.292369828923075d0 -1.6741088357812182d0
-                               8.909104486260928d0 211.6765861544609d0
-                               -1.3439049184367153d0 -10.364363282910663d0
-                               -79.0709314248171d0 -10.652071087998578d0
-                               -9.393948243493877d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (CAUCHY RNG 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.03183098861837907d0)
-                        (MULTIPLE-VALUE-LIST (CAUCHY-PDF 0.0d0 10.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.6475836176504333d0)
-                        (MULTIPLE-VALUE-LIST (CAUCHY-P 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.35241638234956674d0)
-                        (MULTIPLE-VALUE-LIST (CAUCHY-Q 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.9999999999999998d0)
-                        (MULTIPLE-VALUE-LIST
-                         (CAUCHY-PINV 0.6475836176504333d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 1.0000000000000002d0)
-                        (MULTIPLE-VALUE-LIST
-                         (CAUCHY-QINV 0.35241638234956674d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST -0.00811319915595434d0 5.617196410586812d0
+	  12.292369828923075d0 -1.6741088357812182d0
+	  8.909104486260928d0 211.6765861544609d0
+	  -1.3439049184367153d0 -10.364363282910663d0
+	  -79.0709314248171d0 -10.652071087998578d0
+	  -9.393948243493877d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'cauchy :a 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.03183098861837907d0)
+   (MULTIPLE-VALUE-LIST (CAUCHY-PDF 0.0d0 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6475836176504333d0)
+   (MULTIPLE-VALUE-LIST (CAUCHY-P 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.35241638234956674d0)
+   (MULTIPLE-VALUE-LIST (CAUCHY-Q 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9999999999999998d0)
+   (MULTIPLE-VALUE-LIST
+    (CAUCHY-PINV 0.6475836176504333d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0000000000000002d0)
+   (MULTIPLE-VALUE-LIST
+    (CAUCHY-QINV 0.35241638234956674d0 2.0d0))))
 
diff --git a/tests/chi-squared.lisp b/tests/chi-squared.lisp
index ed0448c31f3eb0a1eca0c31b731f531f736aa6f1..e86fab7d612962f85701a46b9d320bcabf116cac 100644
--- a/tests/chi-squared.lisp
+++ b/tests/chi-squared.lisp
@@ -3,33 +3,33 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST CHI-SQUARED
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 13.043328884186328d0 11.427227829236712d0
-                               16.55811815484942d0 7.128795406995407d0
-                               5.120266499239882d0 10.464572605669142d0
-                               5.8126929867006405d0 8.784940866479005d0
-                               7.559275305609187d0 8.35181083950897d0
-                               4.140798004825149d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (CHI-SQUARED RNG 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.43939128946772227d0)
-                        (MULTIPLE-VALUE-LIST (CHI-SQUARED-PDF 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5204998778130462d0)
-                        (MULTIPLE-VALUE-LIST (CHI-SQUARED-P 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.4795001221869538d0)
-                        (MULTIPLE-VALUE-LIST (CHI-SQUARED-Q 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5000000000000003d0)
-                        (MULTIPLE-VALUE-LIST
-                         (CHI-SQUARED-PINV 0.5204998778130463d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5000000000000003d0)
-                        (MULTIPLE-VALUE-LIST
-                         (CHI-SQUARED-QINV 0.4795001221869537d0 1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 13.043328884186328d0 11.427227829236712d0
+	  16.55811815484942d0 7.128795406995407d0
+	  5.120266499239882d0 10.464572605669142d0
+	  5.8126929867006405d0 8.784940866479005d0
+	  7.559275305609187d0 8.35181083950897d0
+	  4.140798004825149d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'chi-squared :nu 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.43939128946772227d0)
+   (MULTIPLE-VALUE-LIST (CHI-SQUARED-PDF 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5204998778130462d0)
+   (MULTIPLE-VALUE-LIST (CHI-SQUARED-P 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4795001221869538d0)
+   (MULTIPLE-VALUE-LIST (CHI-SQUARED-Q 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5000000000000003d0)
+   (MULTIPLE-VALUE-LIST
+    (CHI-SQUARED-PINV 0.5204998778130463d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5000000000000003d0)
+   (MULTIPLE-VALUE-LIST
+    (CHI-SQUARED-QINV 0.4795001221869537d0 1.0d0))))
 
diff --git a/tests/dirichlet.lisp b/tests/dirichlet.lisp
index b26a0a80fe1fd733e446daf7808d74c3bfc3f5b8..b96d9cb76bec987779ae1e0f2265fdd26ae7d032 100644
--- a/tests/dirichlet.lisp
+++ b/tests/dirichlet.lisp
@@ -3,37 +3,32 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST DIRICHLET
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         #(0.12448073544131681d0 0.19182353706734917d0
-                           0.46054388544826397d0 0.22315184204307006d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0))
-                               (ALPHA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(1.0d0 2.0d0 3.0d0 4.0d0)))
-                               (THETA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :DIMENSIONS 4)))
-                           (DIRICHLET RNG ALPHA THETA)
-                           (CL-ARRAY THETA))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 1080.0000000000025d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((ALPHA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(1.0d0 2.0d0 3.0d0 4.0d0)))
-                               (THETA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(1.0d0 2.0d0 3.0d0 4.0d0))))
-                           (DIRICHLET-PDF ALPHA THETA))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 6.984716320118268d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((ALPHA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(1.0d0 2.0d0 3.0d0 4.0d0)))
-                               (THETA
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(1.0d0 2.0d0 3.0d0 4.0d0))))
-                           (DIRICHLET-LOG-PDF ALPHA THETA)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.12448073544131681d0 0.19182353706734917d0
+      0.46054388544826397d0 0.22315184204307006d0))
+   (MULTIPLE-VALUE-LIST
+    (let ((rng (make-random-number-generator +mt19937+ 0))
+	  (alpha #m(1.0d0 2.0d0 3.0d0 4.0d0)))
+      (cl-array (sample rng 'dirichlet :alpha alpha)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1080.0000000000025d0)
+   (MULTIPLE-VALUE-LIST
+    (LET ((ALPHA
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(1.0d0 2.0d0 3.0d0 4.0d0)))
+	  (THETA
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(1.0d0 2.0d0 3.0d0 4.0d0))))
+      (DIRICHLET-PDF ALPHA THETA))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 6.984716320118268d0)
+   (MULTIPLE-VALUE-LIST
+    (LET ((ALPHA
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(1.0d0 2.0d0 3.0d0 4.0d0)))
+	  (THETA
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(1.0d0 2.0d0 3.0d0 4.0d0))))
+      (DIRICHLET-LOG-PDF ALPHA THETA)))))
 
diff --git a/tests/discrete.lisp b/tests/discrete.lisp
index 4de841123326c68c46943782009e5769d4965c9c..b2357d482c1b735263358c6f93aa5b1c9ced98e5 100644
--- a/tests/discrete.lisp
+++ b/tests/discrete.lisp
@@ -3,29 +3,22 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST DISCRETE
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 1 0 1 1 0 1 1 2 1 2 2))
-                        (MULTIPLE-VALUE-LIST
-                         (LET* ((PROBABILITIES
-                                 (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                              '(0.25d0 0.5d0 0.25d0)))
-                                (TABLE (MAKE-DISCRETE-RANDOM PROBABILITIES))
-                                (RNG
-                                 (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (DISCRETE RNG TABLE)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.5d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (LET* ((PROBABILITIES
-                                                                   (MAKE-MARRAY
-                                                                    'DOUBLE-FLOAT
-                                                                    :INITIAL-CONTENTS
-                                                                    '(0.25d0
-                                                                      0.5d0
-                                                                      0.25d0)))
-                                                                  (TABLE
-                                                                   (MAKE-DISCRETE-RANDOM
-                                                                    PROBABILITIES)))
-                                                             (DISCRETE-PDF 1
-                                                                           TABLE)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 1 0 1 1 0 1 1 2 1 2 2))
+   (MULTIPLE-VALUE-LIST
+    (LET* ((PROBABILITIES
+	    (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			 '(0.25d0 0.5d0 0.25d0)))
+	   (TABLE (MAKE-DISCRETE-RANDOM PROBABILITIES))
+	   (RNG
+	    (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'discrete :table table)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5d0)
+   (MULTIPLE-VALUE-LIST
+    (LET* ((PROBABILITIES
+	    (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS '(0.25d0 0.5d0 0.25d0)))
+	   (TABLE (MAKE-DISCRETE-RANDOM PROBABILITIES)))
+      (DISCRETE-PDF 1 TABLE)))))
 
diff --git a/tests/exponential-power.lisp b/tests/exponential-power.lisp
index 80fc3172cad2eeee10efeef0b932572badfbf94d..d50698cca9b27dff82a5468ae15f897da872d1f3 100644
--- a/tests/exponential-power.lisp
+++ b/tests/exponential-power.lisp
@@ -3,28 +3,28 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST EXPONENTIAL-POWER
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 0.09469475592777954d0 -0.06229680875327071d0
-                               1.183985538537803d0 0.5187626019237904d0
-                               0.7053564314063956d0 -0.9033303844569821d0
-                               -1.6947336289940842d0 -0.4803236108055401d0
-                               -0.027641736349912214d0 0.6318391856046153d0
-                               -0.012478875227423025d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (EXPONENTIAL-POWER RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5641895835477557d0)
-                        (MULTIPLE-VALUE-LIST
-                         (EXPONENTIAL-POWER-PDF 0.0d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.9213503964748571d0)
-                        (MULTIPLE-VALUE-LIST
-                         (EXPONENTIAL-POWER-P 1.0d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.07864960352514248d0)
-                        (MULTIPLE-VALUE-LIST
-                         (EXPONENTIAL-POWER-Q 1.0d0 1.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.09469475592777954d0 -0.06229680875327071d0
+	  1.183985538537803d0 0.5187626019237904d0
+	  0.7053564314063956d0 -0.9033303844569821d0
+	  -1.6947336289940842d0 -0.4803236108055401d0
+	  -0.027641736349912214d0 0.6318391856046153d0
+	  -0.012478875227423025d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'exponential-power :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5641895835477557d0)
+   (MULTIPLE-VALUE-LIST
+    (EXPONENTIAL-POWER-PDF 0.0d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9213503964748571d0)
+   (MULTIPLE-VALUE-LIST
+    (EXPONENTIAL-POWER-P 1.0d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07864960352514248d0)
+   (MULTIPLE-VALUE-LIST
+    (EXPONENTIAL-POWER-Q 1.0d0 1.0d0 2.0d0))))
 
diff --git a/tests/exponential.lisp b/tests/exponential.lisp
index b417234c19e8cbaf4f01f64dde645e919bd0c224..83afeed010b17119270a2bd2e60730832710b7d2 100644
--- a/tests/exponential.lisp
+++ b/tests/exponential.lisp
@@ -14,7 +14,7 @@
                         (MULTIPLE-VALUE-LIST
                          (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
                            (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (EXPONENTIAL RNG 10.0d0)))))
+                                 (sample rng 'exponential :mu 10.0d0)))))
                        (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.1d0)
                                                           (MULTIPLE-VALUE-LIST
                                                            (EXPONENTIAL-PDF
diff --git a/tests/fdist.lisp b/tests/fdist.lisp
index d9a63e5b992638134b5ae848945c736179d7f684..2050dc7ce7d621c7de0ab5e08360f0ccc6d517bb 100644
--- a/tests/fdist.lisp
+++ b/tests/fdist.lisp
@@ -3,34 +3,34 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST FDIST
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 18.411308904958695d0 5.120881381831058d0
-                               18.104535265974707d0 0.15934280606960227d0
-                               0.06272171507636773d0 2.2809441726456456d0
-                               0.5259458753395939d0 1.8256109001076744d0
-                               0.845346894458977d0 2.5212086970057763d0
-                               0.5212415547032052d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (FDIST RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.1594719884624466d0)
-                        (MULTIPLE-VALUE-LIST (FDIST-PDF 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.6123724356957948d0)
-                        (MULTIPLE-VALUE-LIST (FDIST-P 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.3876275643042052d0)
-                        (MULTIPLE-VALUE-LIST (FDIST-Q 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.2d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (FDIST-PINV
-                                                            0.612372435695795d0
-                                                            1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 1.1999999999999995d0)
-                        (MULTIPLE-VALUE-LIST
-                         (FDIST-QINV 0.38762756430420503d0 1.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 18.411308904958695d0 5.120881381831058d0
+	  18.104535265974707d0 0.15934280606960227d0
+	  0.06272171507636773d0 2.2809441726456456d0
+	  0.5259458753395939d0 1.8256109001076744d0
+	  0.845346894458977d0 2.5212086970057763d0
+	  0.5212415547032052d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'fdist :nu1 1.0d0 :nu2 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1594719884624466d0)
+   (MULTIPLE-VALUE-LIST (FDIST-PDF 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6123724356957948d0)
+   (MULTIPLE-VALUE-LIST (FDIST-P 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.3876275643042052d0)
+   (MULTIPLE-VALUE-LIST (FDIST-Q 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.2d0)
+				     (MULTIPLE-VALUE-LIST
+				      (FDIST-PINV
+				       0.612372435695795d0
+				       1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.1999999999999995d0)
+   (MULTIPLE-VALUE-LIST
+    (FDIST-QINV 0.38762756430420503d0 1.0d0 2.0d0))))
 
diff --git a/tests/flat.lisp b/tests/flat.lisp
index 49a2b530b3110dfe9a37e9987bff7cf319788474..9f6a9b638b9ac7ec95620a8e1f37724d2c5ba66d 100644
--- a/tests/flat.lisp
+++ b/tests/flat.lisp
@@ -3,38 +3,32 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST FLAT
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 1.999741748906672d0 1.1629098753910512d0
-                               1.2826178052928299d0 1.9472010820172727d0
-                               1.2316565427463502d0 1.4849736143369228d0
-                               1.9574769565369934d0 1.7443053431343287d0
-                               1.540043658344075d0 1.7399529814720154d0
-                               1.7599437981843948d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (FLAT RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.0d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (FLAT-PDF 1.2d0
-                                                                     1.0d0
-                                                                     2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.19999999999999996d0)
-                        (MULTIPLE-VALUE-LIST (FLAT-P 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.8d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (FLAT-Q 1.2d0 1.0d0
-                                                                   2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.2d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (FLAT-PINV
-                                                            0.19999999999999996d0
-                                                            1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.2d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (FLAT-QINV 0.8d0
-                                                                      1.0d0
-                                                                      2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 1.999741748906672d0 1.1629098753910512d0
+	  1.2826178052928299d0 1.9472010820172727d0
+	  1.2316565427463502d0 1.4849736143369228d0
+	  1.9574769565369934d0 1.7443053431343287d0
+	  1.540043658344075d0 1.7399529814720154d0
+	  1.7599437981843948d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'flat :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0d0)
+   (MULTIPLE-VALUE-LIST (FLAT-PDF 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19999999999999996d0)
+   (MULTIPLE-VALUE-LIST (FLAT-P 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8d0)
+   (MULTIPLE-VALUE-LIST (FLAT-Q 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.2d0)
+   (MULTIPLE-VALUE-LIST
+    (FLAT-PINV 0.19999999999999996d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.2d0)
+   (MULTIPLE-VALUE-LIST (FLAT-QINV 0.8d0 1.0d0 2.0d0))))
 
diff --git a/tests/gamma-randist.lisp b/tests/gamma-randist.lisp
index 3f5c4a42f29a465d9216ddd65bd7107fee6bdde9..51142b490605c227c1fd23e2858c0b5712ed8e10 100644
--- a/tests/gamma-randist.lisp
+++ b/tests/gamma-randist.lisp
@@ -3,45 +3,44 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST GAMMA-RANDIST
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 3.012983063768798d0 2.216796987787054d0
-                               5.033971231985026d0 0.6152002566487763d0
-                               0.1876159751197978d0 1.7884608326846099d0
-                               0.30812625873110316d0 1.1328459017528132d0
-                               0.7363931539298727d0 0.9843618987581162d0
-                               0.06871686155296197d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GAMMA-RD RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 3.012983063768798d0 2.216796987787054d0
-                               5.033971231985026d0 0.6152002566487763d0
-                               0.1876159751197978d0 1.7884608326846099d0
-                               0.30812625873110316d0 1.1328459017528132d0
-                               0.7363931539298727d0 0.9843618987581162d0
-                               0.06871686155296197d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GAMMA-MT RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.475614712250357d0)
-                        (MULTIPLE-VALUE-LIST (GAMMA-PDF 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.04877057549928599d0)
-                        (MULTIPLE-VALUE-LIST (GAMMA-P 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.951229424500714d0)
-                        (MULTIPLE-VALUE-LIST (GAMMA-Q 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.10000000000000006d0)
-                        (MULTIPLE-VALUE-LIST
-                         (GAMMA-PINV 0.048770575499286005d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.10000000000000003d0)
-                        (MULTIPLE-VALUE-LIST
-                         (GAMMA-QINV 0.951229424500714d0 1.0d0 2.0d0))))
-
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 3.012983063768798d0 2.216796987787054d0
+	  5.033971231985026d0 0.6152002566487763d0
+	  0.1876159751197978d0 1.7884608326846099d0
+	  0.30812625873110316d0 1.1328459017528132d0
+	  0.7363931539298727d0 0.9843618987581162d0
+	  0.06871686155296197d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'gamma :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 3.012983063768798d0 2.216796987787054d0
+	  5.033971231985026d0 0.6152002566487763d0
+	  0.1876159751197978d0 1.7884608326846099d0
+	  0.30812625873110316d0 1.1328459017528132d0
+	  0.7363931539298727d0 0.9843618987581162d0
+	  0.06871686155296197d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'gamma-mt :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.475614712250357d0)
+   (MULTIPLE-VALUE-LIST (GAMMA-PDF 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.04877057549928599d0)
+   (MULTIPLE-VALUE-LIST (GAMMA-P 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.951229424500714d0)
+   (MULTIPLE-VALUE-LIST (GAMMA-Q 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.10000000000000006d0)
+   (MULTIPLE-VALUE-LIST
+    (GAMMA-PINV 0.048770575499286005d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.10000000000000003d0)
+   (MULTIPLE-VALUE-LIST
+    (GAMMA-QINV 0.951229424500714d0 1.0d0 2.0d0))))
diff --git a/tests/gaussian-bivariate.lisp b/tests/gaussian-bivariate.lisp
index ab12f4b96cfcd14324d2a4ce938d8cb82a862e8f..4f6fa13a01728d24a8be1f654ca11ee84f3a6e32 100644
--- a/tests/gaussian-bivariate.lisp
+++ b/tests/gaussian-bivariate.lisp
@@ -3,22 +3,23 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST GAUSSIAN-BIVARIATE
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST -0.06509716124488897d0 -1.5733207749096374d0
-                               0.27942740172325414d0 1.2021528358889673d0
-                               -0.6041530626907894d0 0.07582702719413444d0
-                               -0.5446229412165632d0 -0.6592026841613081d0
-                               -0.11029516610819164d0 0.17931840412143885d0
-                               2.1025104980291696d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (BIVARIATE-GAUSSIAN RNG 1.0d0 0.75d0
-                                                     0.25d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5548265557970462d0)
-                        (MULTIPLE-VALUE-LIST
-                         (BIVARIATE-GAUSSIAN-PDF 0.25d0 0.5d0 0.25d0 0.4d0
-                                                 0.2d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST -0.06509716124488897d0 -1.5733207749096374d0
+	  0.27942740172325414d0 1.2021528358889673d0
+	  -0.6041530626907894d0 0.07582702719413444d0
+	  -0.5446229412165632d0 -0.6592026841613081d0
+	  -0.11029516610819164d0 0.17931840412143885d0
+	  2.1025104980291696d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample
+	    rng 'gaussian-bivariate
+	    :sigma-x 1.0d0 :sigma-y 0.75d0 :rho 0.25d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5548265557970462d0)
+   (MULTIPLE-VALUE-LIST
+    (BIVARIATE-GAUSSIAN-PDF 0.25d0 0.5d0 0.25d0 0.4d0
+			    0.2d0))))
 
diff --git a/tests/gaussian-tail.lisp b/tests/gaussian-tail.lisp
index 7d7ff5e3b603da83b84760aa5b296c10c3513983..7851a2f0c9d660cb6c83be3844a684aad722072a 100644
--- a/tests/gaussian-tail.lisp
+++ b/tests/gaussian-tail.lisp
@@ -14,7 +14,7 @@
                         (MULTIPLE-VALUE-LIST
                          (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
                            (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GAUSSIAN-TAIL RNG 50.0d0 10.0d0)))))
+                                 (sample rng 'gaussian-tail :a 50.0d0 :sigma 10.0d0)))))
                        (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
                         (LIST 0.18702270877331703d0)
                         (MULTIPLE-VALUE-LIST
@@ -30,7 +30,7 @@
                         (MULTIPLE-VALUE-LIST
                          (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
                            (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (UGAUSSIAN-TAIL RNG 5.0d0)))))
+                                 (sample rng 'ugaussian-tail :a 5.0d0)))))
                        (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
                         (LIST 1.8702270877331704d0)
                         (MULTIPLE-VALUE-LIST (UGAUSSIAN-TAIL-PDF 5.2d0 5.0d0))))
diff --git a/tests/gaussian.lisp b/tests/gaussian.lisp
index 67ef24f3e1b2755ba906ba36f68636783bcbada4..29e6b6920f5f225cf0a47476a75a176339af457b 100644
--- a/tests/gaussian.lisp
+++ b/tests/gaussian.lisp
@@ -14,7 +14,7 @@
                         (MULTIPLE-VALUE-LIST
                          (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
                            (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GAUSSIAN RNG 10.0d0)))))
+                                 (sample rng 'gaussian :sigma 10.0d0)))))
                        (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
                         (LIST 0.039894228040143274d0)
                         (MULTIPLE-VALUE-LIST (GAUSSIAN-PDF 0.0d0 10.0d0)))
@@ -29,7 +29,7 @@
                         (MULTIPLE-VALUE-LIST
                          (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
                            (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GAUSSIAN-ZIGGURAT RNG 10.0d0)))))
+                                 (sample rng 'gaussian-ziggurat :sigma 10.0d0)))))
                        (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
                         (LIST 0.9772498680518208d0)
                         (MULTIPLE-VALUE-LIST (UGAUSSIAN-P 2.0d0)))
diff --git a/tests/geometric.lisp b/tests/geometric.lisp
index 99325708502e31e5de38830f2ced49f6caba4655..9f9774410279068bd06168eedb3199b1a71be87a 100644
--- a/tests/geometric.lisp
+++ b/tests/geometric.lisp
@@ -3,22 +3,18 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST GEOMETRIC
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 1 4 3 1 3 2 1 1 2 1 1))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GEOMETRIC RNG 0.4d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.24d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (GEOMETRIC-PDF 2
-                                                                          0.4d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.64d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (GEOMETRIC-P 2
-                                                                        0.4d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.36d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (GEOMETRIC-Q 2
-                                                                        0.4d0))))
-
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 1 4 3 1 3 2 1 1 2 1 1))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'geometric :probability 0.4d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.24d0)
+   (MULTIPLE-VALUE-LIST (GEOMETRIC-PDF 2 0.4d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.64d0)
+   (MULTIPLE-VALUE-LIST (GEOMETRIC-P 2 0.4d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.36d0)
+   (MULTIPLE-VALUE-LIST (GEOMETRIC-Q 2 0.4d0))))
diff --git a/tests/gumbel1.lisp b/tests/gumbel1.lisp
index 076e757e21adff72ad6425669c52b545469add96..8291ce9a707427cbb40d9c4298e144f1fc564d48 100644
--- a/tests/gumbel1.lisp
+++ b/tests/gumbel1.lisp
@@ -3,33 +3,33 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST GUMBEL1
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 8.954596257487015d0 0.0973051899750762d0
-                               0.45913506233088003d0 3.6074124224293223d0
-                               0.31300027468174807d0 1.0165796949651174d0
-                               3.8292081936610396d0 1.912897393181305d0
-                               1.17748457894919d0 1.893232107970416d0
-                               1.9859118616847695d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GUMBEL1 RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.29625708964974956d0)
-                        (MULTIPLE-VALUE-LIST (GUMBEL1-PDF 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.1637073598773166d0)
-                        (MULTIPLE-VALUE-LIST (GUMBEL1-P 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.8362926401226833d0)
-                        (MULTIPLE-VALUE-LIST (GUMBEL1-Q 0.1d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.10000000000000007d0)
-                        (MULTIPLE-VALUE-LIST
-                         (GUMBEL1-PINV 0.1637073598773166d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.10000000000000028d0)
-                        (MULTIPLE-VALUE-LIST
-                         (GUMBEL1-QINV 0.8362926401226833d0 1.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 8.954596257487015d0 0.0973051899750762d0
+	  0.45913506233088003d0 3.6074124224293223d0
+	  0.31300027468174807d0 1.0165796949651174d0
+	  3.8292081936610396d0 1.912897393181305d0
+	  1.17748457894919d0 1.893232107970416d0
+	  1.9859118616847695d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'gumbel1 :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.29625708964974956d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL1-PDF 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1637073598773166d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL1-P 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8362926401226833d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL1-Q 0.1d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.10000000000000007d0)
+   (MULTIPLE-VALUE-LIST
+    (GUMBEL1-PINV 0.1637073598773166d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.10000000000000028d0)
+   (MULTIPLE-VALUE-LIST
+    (GUMBEL1-QINV 0.8362926401226833d0 1.0d0 2.0d0))))
 
diff --git a/tests/gumbel2.lisp b/tests/gumbel2.lisp
index 93ee230256805c3c2587e36f285bec72c9a580f6..90ac66762ea1b7c96a2e6222b4528dee06e8223a 100644
--- a/tests/gumbel2.lisp
+++ b/tests/gumbel2.lisp
@@ -3,34 +3,32 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST GUMBEL2
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 7743.400858519516d0 1.102196701680339d0
-                               1.5827044520998628d0 36.87052393317972d0
-                               1.3675219066608615d0 2.7637257945633085d0
-                               46.026080060263446d0 6.772683525074477d0
-                               3.2461983686562204d0 6.640797807286079d0
-                               7.285687897019733d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (GUMBEL2 RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.053625603682851145d0)
-                        (MULTIPLE-VALUE-LIST (GUMBEL2-PDF 5.0d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.8187307530779818d0)
-                        (MULTIPLE-VALUE-LIST (GUMBEL2-P 10.0d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.18126924692201815d0)
-                        (MULTIPLE-VALUE-LIST (GUMBEL2-Q 10.0d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 9.999999999999998d0)
-                        (MULTIPLE-VALUE-LIST
-                         (GUMBEL2-PINV 0.8187307530779818d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 10.0d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (GUMBEL2-QINV
-                                                            0.18126924692201815d0
-                                                            1.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 7743.400858519516d0 1.102196701680339d0
+	  1.5827044520998628d0 36.87052393317972d0
+	  1.3675219066608615d0 2.7637257945633085d0
+	  46.026080060263446d0 6.772683525074477d0
+	  3.2461983686562204d0 6.640797807286079d0
+	  7.285687897019733d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'gumbel2 :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.053625603682851145d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL2-PDF 5.0d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8187307530779818d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL2-P 10.0d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.18126924692201815d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL2-Q 10.0d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 9.999999999999998d0)
+   (MULTIPLE-VALUE-LIST
+    (GUMBEL2-PINV 0.8187307530779818d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 10.0d0)
+   (MULTIPLE-VALUE-LIST (GUMBEL2-QINV 0.18126924692201815d0 1.0d0 2.0d0))))
 
diff --git a/tests/hypergeometric-randist.lisp b/tests/hypergeometric-randist.lisp
index 3802af8befa29ed23cd50406054a425c8056eb35..1b7fdc7aa1ec6d98db2d52fc4f003a0676f8929c 100644
--- a/tests/hypergeometric-randist.lisp
+++ b/tests/hypergeometric-randist.lisp
@@ -3,19 +3,19 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST HYPERGEOMETRIC-RANDIST
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 2 1 0 0 1 1 3 1 0 1 3))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (HYPERGEOMETRIC RNG 3 6 3)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.35714285714285693d0)
-                        (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-PDF 0 2 6 3)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.892857142857143d0)
-                        (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-P 1 2 6 3)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.10714285714285704d0)
-                        (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-Q 1 2 6 3))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 2 1 0 0 1 1 3 1 0 1 3))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'hypergeometric :n1 3 :n2 6 :tt 3)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.35714285714285693d0)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-PDF 0 2 6 3)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.892857142857143d0)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-P 1 2 6 3)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.10714285714285704d0)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-Q 1 2 6 3))))
 
diff --git a/tests/landau.lisp b/tests/landau.lisp
index 7b5d27565dcd320881380172670c505bbced3404..8d9ee91819b58b5c7a846e4382dcd04681781444 100644
--- a/tests/landau.lisp
+++ b/tests/landau.lisp
@@ -3,18 +3,18 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST LANDAU
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 3880.037426254597d0 -0.6953200314545297d0
-                               -0.02354364646600932d0 21.329209630030316d0
-                               -0.3062224704714883d0 1.2424186669362394d0
-                               26.146168479649152d0 4.337217640968217d0
-                               1.6799546281085946d0 4.2475719218268395d0
-                               4.681506208977819d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT (LANDAU RNG)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.17331968995860203d0)
-                        (MULTIPLE-VALUE-LIST (LANDAU-PDF 0.25d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 3880.037426254597d0 -0.6953200314545297d0
+	  -0.02354364646600932d0 21.329209630030316d0
+	  -0.3062224704714883d0 1.2424186669362394d0
+	  26.146168479649152d0 4.337217640968217d0
+	  1.6799546281085946d0 4.2475719218268395d0
+	  4.681506208977819d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT (sample rng 'landau)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.17331968995860203d0)
+   (MULTIPLE-VALUE-LIST (LANDAU-PDF 0.25d0))))
 
diff --git a/tests/laplace.lisp b/tests/laplace.lisp
index e542c4636a39e42294be304cb9046cf951a83028..960325d1860e821a006215206a74571f3d6deb74 100644
--- a/tests/laplace.lisp
+++ b/tests/laplace.lisp
@@ -14,7 +14,7 @@
                         (MULTIPLE-VALUE-LIST
                          (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
                            (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (LAPLACE RNG 10.0d0)))))
+                                 (sample rng 'laplace :a 10.0d0)))))
                        (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.05d0)
                                                           (MULTIPLE-VALUE-LIST
                                                            (LAPLACE-PDF 0.0d0
diff --git a/tests/levy.lisp b/tests/levy.lisp
index b97aa985715a22d151f042a1042c18a4b9172dc1..a74dfa2f5cb9319cc3e4f5c87440a12e2e2e1beb 100644
--- a/tests/levy.lisp
+++ b/tests/levy.lisp
@@ -3,28 +3,28 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST LEVY
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 2.6941098332360465d0 -0.29395438644676647d0
-                               -1.2703401352272083d0 1.0771538640113538d0
-                               0.13771218406916103d0 0.9419728438107844d0
-                               -0.5107134674789159d0 0.1648207853689268d0
-                               -0.14857899041035147d0 -1.9074885744364487d0
-                               -2.086195213997167d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (LEVY RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 2.6941098332360465d0 -0.2939543864467665d0
-                               -1.2703401352272083d0 1.0771538640113538d0
-                               0.13771218406916097d0 0.9419728438107844d0
-                               -0.510713467478916d0 0.1648207853689266d0
-                               -0.14857899041035158d0 -1.907488574436449d0
-                               -2.086195213997167d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (LEVY-SKEW RNG 1.0d0 2.0d0 1.0d0))))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 2.6941098332360465d0 -0.29395438644676647d0
+	  -1.2703401352272083d0 1.0771538640113538d0
+	  0.13771218406916103d0 0.9419728438107844d0
+	  -0.5107134674789159d0 0.1648207853689268d0
+	  -0.14857899041035147d0 -1.9074885744364487d0
+	  -2.086195213997167d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'levy :c 1.0d0 :alpha 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 2.6941098332360465d0 -0.2939543864467665d0
+	  -1.2703401352272083d0 1.0771538640113538d0
+	  0.13771218406916097d0 0.9419728438107844d0
+	  -0.510713467478916d0 0.1648207853689266d0
+	  -0.14857899041035158d0 -1.907488574436449d0
+	  -2.086195213997167d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'levy-skew :c 1.0d0 :alpha 2.0d0 :beta 1.0d0))))))
 
diff --git a/tests/logarithmic.lisp b/tests/logarithmic.lisp
index 6dce09cc36d321d589d7f1e196c2fd7589a4b9fb..6e408ec6e0d118c875b3400f9046e349f9179fcf 100644
--- a/tests/logarithmic.lisp
+++ b/tests/logarithmic.lisp
@@ -3,13 +3,13 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST LOGARITHMIC
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 1 3 1 4 1 1 2 1 1 5 2))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (LOGARITHMIC RNG 0.9d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.15660921511769743d0)
-                        (MULTIPLE-VALUE-LIST (LOGARITHMIC-PDF 2 0.4d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 1 3 1 4 1 1 2 1 1 5 2))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'logarithmic :probability 0.9d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.15660921511769743d0)
+   (MULTIPLE-VALUE-LIST (LOGARITHMIC-PDF 2 0.4d0))))
 
diff --git a/tests/logistic.lisp b/tests/logistic.lisp
index 665c89f5636235a49bbc9b53b3a7950098072d9f..f3f17e8b85a5e82fa9f091988e88d6890dc6c57b 100644
--- a/tests/logistic.lisp
+++ b/tests/logistic.lisp
@@ -3,33 +3,31 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST LOGISTIC
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 82.6131993192451d0 -16.367346042668906d0
-                               -9.31513272043762d0 28.87020708710654d0
-                               -11.989809875784625d0 -0.6012364762000397d0
-                               31.142555263622d0 10.684673721048895d0
-                               1.6051840954024446d0 10.457241904701199d0
-                               11.523714106294097d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (LOGISTIC RNG 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.2350037122015945d0)
-                        (MULTIPLE-VALUE-LIST (LOGISTIC-PDF 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.6224593312018546d0)
-                        (MULTIPLE-VALUE-LIST (LOGISTIC-P 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.37754066879814546d0)
-                        (MULTIPLE-VALUE-LIST (LOGISTIC-Q 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5000000000000001d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LOGISTIC-PINV 0.6224593312018546d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.4999999999999998d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LOGISTIC-QINV 0.37754066879814546d0 1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 82.6131993192451d0 -16.367346042668906d0
+	  -9.31513272043762d0 28.87020708710654d0
+	  -11.989809875784625d0 -0.6012364762000397d0
+	  31.142555263622d0 10.684673721048895d0
+	  1.6051840954024446d0 10.457241904701199d0
+	  11.523714106294097d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'logistic :a 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.2350037122015945d0)
+   (MULTIPLE-VALUE-LIST (LOGISTIC-PDF 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6224593312018546d0)
+   (MULTIPLE-VALUE-LIST (LOGISTIC-P 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.37754066879814546d0)
+   (MULTIPLE-VALUE-LIST (LOGISTIC-Q 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5000000000000001d0)
+   (MULTIPLE-VALUE-LIST (LOGISTIC-PINV 0.6224593312018546d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4999999999999998d0)
+   (MULTIPLE-VALUE-LIST (LOGISTIC-QINV 0.37754066879814546d0 1.0d0))))
 
diff --git a/tests/lognormal.lisp b/tests/lognormal.lisp
index b9ac8e19f74c166bb266a33ccfa67984214911f0..446394ca328c060f1631b32c1d93a22564807b38 100644
--- a/tests/lognormal.lisp
+++ b/tests/lognormal.lisp
@@ -3,35 +3,35 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST LOGNORMAL
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 2.386447068127768d0 0.11687602167359055d0
-                               4.753374578796263d0 30.093393775755004d0
-                               0.8119584375760986d0 3.163421055157545d0
-                               0.9146206567715651d0 0.727307901065758d0
-                               2.180184852178898d0 3.8908856616896017d0
-                               182.18469788916977d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (LOGNORMAL RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.15289833965691607d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LOGNORMAL-PDF 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.3413288272347351d0)
-                        (MULTIPLE-VALUE-LIST (LOGNORMAL-P 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.6586711727652649d0)
-                        (MULTIPLE-VALUE-LIST (LOGNORMAL-Q 1.2d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 1.2000000000000004d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LOGNORMAL-PINV 0.3413288272347352d0 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.2d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (LOGNORMAL-QINV
-                                                            0.6586711727652649d0
-                                                            1.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 2.386447068127768d0 0.11687602167359055d0
+	  4.753374578796263d0 30.093393775755004d0
+	  0.8119584375760986d0 3.163421055157545d0
+	  0.9146206567715651d0 0.727307901065758d0
+	  2.180184852178898d0 3.8908856616896017d0
+	  182.18469788916977d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'lognormal :zeta 1.0d0 :sigma 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.15289833965691607d0)
+   (MULTIPLE-VALUE-LIST
+    (LOGNORMAL-PDF 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.3413288272347351d0)
+   (MULTIPLE-VALUE-LIST (LOGNORMAL-P 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6586711727652649d0)
+   (MULTIPLE-VALUE-LIST (LOGNORMAL-Q 1.2d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.2000000000000004d0)
+   (MULTIPLE-VALUE-LIST
+    (LOGNORMAL-PINV 0.3413288272347352d0 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 1.2d0)
+				     (MULTIPLE-VALUE-LIST
+				      (LOGNORMAL-QINV
+				       0.6586711727652649d0
+				       1.0d0 2.0d0))))
 
diff --git a/tests/multinomial.lisp b/tests/multinomial.lisp
index adf8786e5606a32b7d97c191954fc1a4a883c77e..5a5163a301de0a31eeac944979d330da45811124 100644
--- a/tests/multinomial.lisp
+++ b/tests/multinomial.lisp
@@ -3,47 +3,30 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST MULTINOMIAL
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST #(5 0 1 2))
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (LET ((RNG
-                                                                  (MAKE-RANDOM-NUMBER-GENERATOR
-                                                                   +MT19937+
-                                                                   0))
-                                                                 (P
-                                                                  (MAKE-MARRAY
-                                                                   'DOUBLE-FLOAT
-                                                                   :INITIAL-CONTENTS
-                                                                   '(0.1d0
-                                                                     0.2d0
-                                                                     0.3d0
-                                                                     0.4d0)))
-                                                                 (N
-                                                                  (MAKE-MARRAY
-                                                                   '(SIGNED-BYTE
-                                                                     32)
-                                                                   :DIMENSIONS
-                                                                   4)))
-                                                             (MULTINOMIAL RNG 8
-                                                                          P N)
-                                                             (CL-ARRAY N))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 8.064000000000026d-5)
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((P
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(0.1d0 0.2d0 0.3d0 0.4d0)))
-                               (N
-                                (MAKE-MARRAY '(SIGNED-BYTE 32)
-                                             :INITIAL-CONTENTS '(5 0 1 2))))
-                           (MULTINOMIAL-PDF P N))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST -9.425515753641212d0)
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((P
-                                (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
-                                             '(0.1d0 0.2d0 0.3d0 0.4d0)))
-                               (N
-                                (MAKE-MARRAY '(SIGNED-BYTE 32)
-                                             :INITIAL-CONTENTS '(5 0 1 2))))
-                           (MULTINOMIAL-LOG-PDF P N)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(5 0 1 2))
+   (MULTIPLE-VALUE-LIST
+    (let ((rng (make-random-number-generator +mt19937+ 0))
+	  (p #m(0.1d0 0.2d0 0.3d0 0.4d0)))
+      (cl-array (sample rng 'multinomial :sum 8 :probabilities p)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 8.064000000000026d-5)
+   (MULTIPLE-VALUE-LIST
+    (LET ((P
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(0.1d0 0.2d0 0.3d0 0.4d0)))
+	  (N
+	   (MAKE-MARRAY '(SIGNED-BYTE 32)
+			:INITIAL-CONTENTS '(5 0 1 2))))
+      (MULTINOMIAL-PDF P N))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -9.425515753641212d0)
+   (MULTIPLE-VALUE-LIST
+    (LET ((P
+	   (MAKE-MARRAY 'DOUBLE-FLOAT :INITIAL-CONTENTS
+			'(0.1d0 0.2d0 0.3d0 0.4d0)))
+	  (N
+	   (MAKE-MARRAY '(SIGNED-BYTE 32)
+			:INITIAL-CONTENTS '(5 0 1 2))))
+      (MULTINOMIAL-LOG-PDF P N)))))
 
diff --git a/tests/negative-binomial.lisp b/tests/negative-binomial.lisp
index 99174922ecdd6794af4abb498fc622b4a98090f7..72a18a49392242a5fe5dbf7fa8ad6deacc7c48db 100644
--- a/tests/negative-binomial.lisp
+++ b/tests/negative-binomial.lisp
@@ -3,37 +3,37 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST NEGATIVE-BINOMIAL
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 15 21 19 15 8 18 23 18 33 16 10))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (NEGATIVE-BINOMIAL RNG 0.4d0 12.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.0056984767089869d0)
-                        (MULTIPLE-VALUE-LIST
-                         (NEGATIVE-BINOMIAL-PDF 5 0.4d0 12.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.010594202555514881d0)
-                        (MULTIPLE-VALUE-LIST
-                         (NEGATIVE-BINOMIAL-P 5 0.4d0 12.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.9894057974444851d0)
-                        (MULTIPLE-VALUE-LIST
-                         (NEGATIVE-BINOMIAL-Q 5 0.4d0 12.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 15 21 19 15 8 18 23 18 33 16 10))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (PASCAL RNG 0.4d0 12)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.0056984767089869d0)
-                        (MULTIPLE-VALUE-LIST (PASCAL-PDF 5 0.4d0 12)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.010594202555514881d0)
-                        (MULTIPLE-VALUE-LIST (PASCAL-P 5 0.4d0 12)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.9894057974444851d0)
-                        (MULTIPLE-VALUE-LIST (PASCAL-Q 5 0.4d0 12))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 15 21 19 15 8 18 23 18 33 16 10))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'negative-binomial :probability 0.4d0 :n 12.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0056984767089869d0)
+   (MULTIPLE-VALUE-LIST
+    (NEGATIVE-BINOMIAL-PDF 5 0.4d0 12.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.010594202555514881d0)
+   (MULTIPLE-VALUE-LIST
+    (NEGATIVE-BINOMIAL-P 5 0.4d0 12.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9894057974444851d0)
+   (MULTIPLE-VALUE-LIST
+    (NEGATIVE-BINOMIAL-Q 5 0.4d0 12.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 15 21 19 15 8 18 23 18 33 16 10))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'pascal :probability 0.4d0 :n 12)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0056984767089869d0)
+   (MULTIPLE-VALUE-LIST (PASCAL-PDF 5 0.4d0 12)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.010594202555514881d0)
+   (MULTIPLE-VALUE-LIST (PASCAL-P 5 0.4d0 12)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9894057974444851d0)
+   (MULTIPLE-VALUE-LIST (PASCAL-Q 5 0.4d0 12))))
 
diff --git a/tests/pareto.lisp b/tests/pareto.lisp
index 661a543b12b6673fa18995961d3143f7a7d5e753..595abe1c558814b19b19e484eda463f523f09f9b 100644
--- a/tests/pareto.lisp
+++ b/tests/pareto.lisp
@@ -3,33 +3,33 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST PARETO
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 2.0005166356083666d0 12.276726596218747d0
-                               7.076694965937239d0 2.111484074469764d0
-                               8.633470811095883d0 4.123935696449152d0
-                               2.0888231161547828d0 2.6870692498025632d0
-                               3.703404287965457d0 2.7028744394290123d0
-                               2.631773566385122d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (PARETO RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5116034405707658d0)
-                        (MULTIPLE-VALUE-LIST (PARETO-PDF 1.5d0 1.3d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.5168849835182453d0)
-                        (MULTIPLE-VALUE-LIST (PARETO-P 3.5d0 1.3d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.4831150164817547d0)
-                        (MULTIPLE-VALUE-LIST (PARETO-Q 3.5d0 1.3d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 3.5000000000000004d0)
-                        (MULTIPLE-VALUE-LIST
-                         (PARETO-PINV 0.5168849835182453d0 1.3d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 3.5000000000000004d0)
-                        (MULTIPLE-VALUE-LIST
-                         (PARETO-QINV 0.4831150164817547d0 1.3d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 2.0005166356083666d0 12.276726596218747d0
+	  7.076694965937239d0 2.111484074469764d0
+	  8.633470811095883d0 4.123935696449152d0
+	  2.0888231161547828d0 2.6870692498025632d0
+	  3.703404287965457d0 2.7028744394290123d0
+	  2.631773566385122d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'pareto :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5116034405707658d0)
+   (MULTIPLE-VALUE-LIST (PARETO-PDF 1.5d0 1.3d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5168849835182453d0)
+   (MULTIPLE-VALUE-LIST (PARETO-P 3.5d0 1.3d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4831150164817547d0)
+   (MULTIPLE-VALUE-LIST (PARETO-Q 3.5d0 1.3d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.5000000000000004d0)
+   (MULTIPLE-VALUE-LIST
+    (PARETO-PINV 0.5168849835182453d0 1.3d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.5000000000000004d0)
+   (MULTIPLE-VALUE-LIST
+    (PARETO-QINV 0.4831150164817547d0 1.3d0 2.0d0))))
 
diff --git a/tests/poisson.lisp b/tests/poisson.lisp
index 47ba650012db8e76728144dfdb64f33eda11123a..6874a07abaf77d9a67696687430fbeeef3e88580 100644
--- a/tests/poisson.lisp
+++ b/tests/poisson.lisp
@@ -3,19 +3,19 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST POISSON
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST (LIST 15 6 9 9 5 8 11 9 11 5 10))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (POISSON RNG 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.11259903214902009d0)
-                        (MULTIPLE-VALUE-LIST (POISSON-PDF 8 10.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.3328196787507177d0)
-                        (MULTIPLE-VALUE-LIST (POISSON-P 8 10.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.6671803212492823d0)
-                        (MULTIPLE-VALUE-LIST (POISSON-Q 8 10.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST (LIST 15 6 9 9 5 8 11 9 11 5 10))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'poisson :mu 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.11259903214902009d0)
+   (MULTIPLE-VALUE-LIST (POISSON-PDF 8 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.3328196787507177d0)
+   (MULTIPLE-VALUE-LIST (POISSON-P 8 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6671803212492823d0)
+   (MULTIPLE-VALUE-LIST (POISSON-Q 8 10.0d0))))
 
diff --git a/tests/random-number-generators.lisp b/tests/random-number-generators.lisp
index 4327715499237791eace603b3f869f4e79764be4..eacb324b3edd120b339b4c63901acd2cb51892bb 100644
--- a/tests/random-number-generators.lisp
+++ b/tests/random-number-generators.lisp
@@ -9,7 +9,7 @@
    (MULTIPLE-VALUE-LIST
     (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
       (LOOP FOR I FROM 0 TO 10 COLLECT
-	   (UNIFORM-FIXNUM RNG 1000)))))
+	   (sample rng 'uniform-fixnum :upperbound 1000)))))
   (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
    (LIST
     (LIST 0.11177622997750353d0 0.9591667949963206d0
@@ -20,5 +20,4 @@
 	  0.9311084621265104d0))
    (MULTIPLE-VALUE-LIST
     (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +CMRG+ 0)))
-      (LOOP FOR I FROM 0 TO 10 COLLECT (UNIFORM RNG))))))
-
+      (LOOP FOR I FROM 0 TO 10 COLLECT (sample rng 'uniform))))))
diff --git a/tests/rayleigh-tail.lisp b/tests/rayleigh-tail.lisp
index b7ced5a26fc6126aa9750668f1b3622357019e32..4b486f814673922234c50b5aa3c49a3bc8fb06dd 100644
--- a/tests/rayleigh-tail.lisp
+++ b/tests/rayleigh-tail.lisp
@@ -3,20 +3,20 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST RAYLEIGH-TAIL
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 1.0255032370386696d0 19.0764679267351d0
-                               15.928966102255199d0 3.4422048899106383d0
-                               17.131838333441106d0 12.071957529361999d0
-                               3.112992916690818d0 7.749889301203328d0
-                               11.145450138119857d0 7.825198187316554d0
-                               7.476774681552917d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (RAYLEIGH-TAIL RNG 1.0d0 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.10224317624874313d0)
-                        (MULTIPLE-VALUE-LIST
-                         (RAYLEIGH-TAIL-PDF 0.25d0 -2.0d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 1.0255032370386696d0 19.0764679267351d0
+	  15.928966102255199d0 3.4422048899106383d0
+	  17.131838333441106d0 12.071957529361999d0
+	  3.112992916690818d0 7.749889301203328d0
+	  11.145450138119857d0 7.825198187316554d0
+	  7.476774681552917d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'rayleigh-tail :a 1.0d0 :sigma 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.10224317624874313d0)
+   (MULTIPLE-VALUE-LIST
+    (RAYLEIGH-TAIL-PDF 0.25d0 -2.0d0 2.0d0))))
 
diff --git a/tests/rayleigh.lisp b/tests/rayleigh.lisp
index 7f5e9719982891d1bab926bbb4d4808f9645eef6..57c78a96dabaaf1b9f98f6f2906caafa5160a565 100644
--- a/tests/rayleigh.lisp
+++ b/tests/rayleigh.lisp
@@ -3,33 +3,33 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST RAYLEIGH
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 0.22728151965522753d0 19.05023959323748d0
-                               15.897545756713367d0 3.2937477899992147d0
-                               17.102628005168157d0 12.030467929000928d0
-                               2.9480035446666624d0 7.6851014424603274d0
-                               11.100498132125239d0 7.76103902005281d0
-                               7.409599155063027d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (RAYLEIGH RNG 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.4412484512922977d0)
-                        (MULTIPLE-VALUE-LIST (RAYLEIGH-PDF 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.1175030974154046d0)
-                        (MULTIPLE-VALUE-LIST (RAYLEIGH-P 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.8824969025845955d0)
-                        (MULTIPLE-VALUE-LIST (RAYLEIGH-Q 1.0d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 1.0000000000000002d0)
-                        (MULTIPLE-VALUE-LIST
-                         (RAYLEIGH-PINV 0.1175030974154046d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.9999999999999998d0)
-                        (MULTIPLE-VALUE-LIST
-                         (RAYLEIGH-QINV 0.8824969025845955d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.22728151965522753d0 19.05023959323748d0
+	  15.897545756713367d0 3.2937477899992147d0
+	  17.102628005168157d0 12.030467929000928d0
+	  2.9480035446666624d0 7.6851014424603274d0
+	  11.100498132125239d0 7.76103902005281d0
+	  7.409599155063027d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'rayleigh :sigma 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4412484512922977d0)
+   (MULTIPLE-VALUE-LIST (RAYLEIGH-PDF 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1175030974154046d0)
+   (MULTIPLE-VALUE-LIST (RAYLEIGH-P 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8824969025845955d0)
+   (MULTIPLE-VALUE-LIST (RAYLEIGH-Q 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0000000000000002d0)
+   (MULTIPLE-VALUE-LIST
+    (RAYLEIGH-PINV 0.1175030974154046d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9999999999999998d0)
+   (MULTIPLE-VALUE-LIST
+    (RAYLEIGH-QINV 0.8824969025845955d0 2.0d0))))
 
diff --git a/tests/shuffling-sampling.lisp b/tests/shuffling-sampling.lisp
index dc24734751427d70de1d2d7c2f38c44d2f7d4ecf..e15ef2fb3180727e69af27c212b6b4c5d2d3b0c4 100644
--- a/tests/shuffling-sampling.lisp
+++ b/tests/shuffling-sampling.lisp
@@ -3,49 +3,21 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST SHUFFLING-SAMPLING
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST #(4 3 6 1 5 7 2 8))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0))
-                               (V1
-                                (MAKE-MARRAY '(SIGNED-BYTE 32)
-                                             :INITIAL-CONTENTS
-                                             '(1 2 3 4 5 6 7 8))))
-                           (SHUFFLE RNG V1)
-                           (CL-ARRAY V1))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST #(2 3 5 8))
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (LET ((RNG
-                                                                  (MAKE-RANDOM-NUMBER-GENERATOR
-                                                                   +MT19937+
-                                                                   0))
-                                                                 (V1
-                                                                  (MAKE-MARRAY
-                                                                   '(SIGNED-BYTE
-                                                                     32)
-                                                                   :INITIAL-CONTENTS
-                                                                   '(1 2 3 4 5
-                                                                     6 7 8)))
-                                                                 (V2
-                                                                  (MAKE-MARRAY
-                                                                   '(SIGNED-BYTE
-                                                                     32)
-                                                                   :DIMENSIONS
-                                                                   4)))
-                                                             (CHOOSE-RANDOM RNG
-                                                                            V2
-                                                                            V1)
-                                                             (CL-ARRAY V2))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST #(8 2 3 8 2 4 8 6 5 6))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0))
-                               (V1
-                                (MAKE-MARRAY '(SIGNED-BYTE 32)
-                                             :INITIAL-CONTENTS
-                                             '(1 2 3 4 5 6 7 8)))
-                               (V2
-                                (MAKE-MARRAY '(SIGNED-BYTE 32) :DIMENSIONS 10)))
-                           (RANDOM-SAMPLE RNG V2 V1)
-                           (CL-ARRAY V2)))))
-
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(4 3 6 1 5 7 2 8))
+   (MULTIPLE-VALUE-LIST
+    (let ((rng (make-random-number-generator +mt19937+ 0))
+	(v1 #31m(1 2 3 4 5 6 7 8)))
+   (cl-array (sample rng 'shuffle :base v1)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(2 3 5 8))
+   (MULTIPLE-VALUE-LIST
+    (let ((rng (make-random-number-generator +mt19937+ 0))
+	(v1 #31m(1 2 3 4 5 6 7 8)))
+   (cl-array (sample rng 'choose-random :src v1 :dest 4)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(8 2 3 8 2 4 8 6 5 6))
+   (MULTIPLE-VALUE-LIST
+    (let ((rng (make-random-number-generator +mt19937+ 0))
+	(v1 #31m(1 2 3 4 5 6 7 8)))
+   (cl-array (sample rng 'random-sample :src v1 :dest 10))))))
diff --git a/tests/spherical-vector.lisp b/tests/spherical-vector.lisp
index b99939ec1fd17b7b610df178a1f95cad0d90b7ad..f15c8500cd7a252521067f516c247cb15d2c7baa 100644
--- a/tests/spherical-vector.lisp
+++ b/tests/spherical-vector.lisp
@@ -3,38 +3,37 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST SPHERICAL-VECTOR
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST -0.617745613497854d0 -0.7863779988047479d0
-                               0.993748310886084d0 0.1116436053298841d0
-                               -0.9458104280982743d0 0.3247193158722761d0
-                               0.45726622946182216d0 0.8893298574734622d0
-                               -0.46325616159849964d0 -0.8862244234622655d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 4 APPEND
-                                 (MULTIPLE-VALUE-LIST (DIRECTION-2D RNG))))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 0.9999986835208556d0 -0.0016226387631051197d0
-                               0.5203010106077766d0 0.8539829379797504d0
-                               -0.2035120531038584d0 0.9790724407527016d0
-                               0.9454753227485545d0 -0.3256937427607672d0
-                               0.11500033916619544d0 0.9933654523848008d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 4 APPEND
-                                 (MULTIPLE-VALUE-LIST
-                                  (DIRECTION-2D-TRIG-METHOD RNG))))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST -0.09129925750445994d0 0.18782185357162273d0
-                               0.977950610665004d0 -0.9051182961559773d0
-                               -0.050683764485791594d0 -0.4221279734645046d0
-                               0.13993766535985133d0 0.8385462620524484d0
-                               -0.526552576872909d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 2 APPEND
-                                 (MULTIPLE-VALUE-LIST (DIRECTION-3D RNG)))))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST -0.617745613497854d0 -0.7863779988047479d0
+	  0.993748310886084d0 0.1116436053298841d0
+	  -0.9458104280982743d0 0.3247193158722761d0
+	  0.45726622946182216d0 0.8893298574734622d0
+	  -0.46325616159849964d0 -0.8862244234622655d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 4 APPEND
+	   (MULTIPLE-VALUE-LIST (sample rng 'direction-2d))))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.9999986835208556d0 -0.0016226387631051197d0
+	  0.5203010106077766d0 0.8539829379797504d0
+	  -0.2035120531038584d0 0.9790724407527016d0
+	  0.9454753227485545d0 -0.3256937427607672d0
+	  0.11500033916619544d0 0.9933654523848008d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 4 APPEND
+	   (MULTIPLE-VALUE-LIST (sample rng 'direction-2d-trig-method))))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST -0.09129925750445994d0 0.18782185357162273d0
+	  0.977950610665004d0 -0.9051182961559773d0
+	  -0.050683764485791594d0 -0.4221279734645046d0
+	  0.13993766535985133d0 0.8385462620524484d0
+	  -0.526552576872909d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 2 APPEND
+	   (MULTIPLE-VALUE-LIST (sample rng 'direction-3d)))))))
 
diff --git a/tests/tdist.lisp b/tests/tdist.lisp
index 41c44b038e4e71bbf8943ca93925f1ac10e838a8..972d9bfacffd99431dc1cf3133ddaa8b5bc4cd4f 100644
--- a/tests/tdist.lisp
+++ b/tests/tdist.lisp
@@ -3,34 +3,32 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST TDIST
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 0.14989366374481017d0 0.6794142879291215d0
-                               -1.615833951108472d0 -1.6008862825783456d0
-                               -1.7010935505767397d0 -0.04370959749808691d0
-                               0.12761159276595174d0 -0.019731218255494867d0
-                               -0.6534666117199732d0 0.2035771324523077d0
-                               1.77650300477611d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (TDIST RNG 10.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.2546479089470325d0)
-                        (MULTIPLE-VALUE-LIST (TDIST-PDF 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.6475836176504332d0)
-                        (MULTIPLE-VALUE-LIST (TDIST-P 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.35241638234956685d0)
-                        (MULTIPLE-VALUE-LIST (TDIST-Q 0.5d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.500000000000001d0)
-                        (MULTIPLE-VALUE-LIST
-                         (TDIST-PINV 0.6475836176504334d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.5d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (TDIST-QINV
-                                                            0.3524163823495667d0
-                                                            1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.14989366374481017d0 0.6794142879291215d0
+	  -1.615833951108472d0 -1.6008862825783456d0
+	  -1.7010935505767397d0 -0.04370959749808691d0
+	  0.12761159276595174d0 -0.019731218255494867d0
+	  -0.6534666117199732d0 0.2035771324523077d0
+	  1.77650300477611d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'tdist :nu 10.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.2546479089470325d0)
+   (MULTIPLE-VALUE-LIST (TDIST-PDF 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6475836176504332d0)
+   (MULTIPLE-VALUE-LIST (TDIST-P 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.35241638234956685d0)
+   (MULTIPLE-VALUE-LIST (TDIST-Q 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.500000000000001d0)
+   (MULTIPLE-VALUE-LIST
+    (TDIST-PINV 0.6475836176504334d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5d0)
+   (MULTIPLE-VALUE-LIST (TDIST-QINV 0.3524163823495667d0 1.0d0))))
 
diff --git a/tests/weibull.lisp b/tests/weibull.lisp
index e23e307a8efdc409fa4a4cf56a49282faf3559a2..264d4305ef042d3ef55c3bdfe30585f58c9de310 100644
--- a/tests/weibull.lisp
+++ b/tests/weibull.lisp
@@ -3,34 +3,32 @@
 (in-package :gsl)
 
 (LISP-UNIT:DEFINE-TEST WEIBULL
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST
-                         (LIST 0.0160712303786595d0 1.347055359960668d0
-                               1.1241262408795445d0 0.2329031397826649d0
-                               1.209338423856536d0 0.8506825453443836d0
-                               0.20845532973957762d0 0.5434187344070215d0
-                               0.7849237503774361d0 0.5487883320132739d0
-                               0.5239377808419179d0))
-                        (MULTIPLE-VALUE-LIST
-                         (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
-                           (LOOP FOR I FROM 0 TO 10 COLLECT
-                                 (WEIBULL RNG 1.0d0 2.0d0)))))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.24263174972226745d0)
-                        (MULTIPLE-VALUE-LIST (WEIBULL-PDF 1.5d0 1.3d0 1.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 0.9992887742799077d0)
-                        (MULTIPLE-VALUE-LIST (WEIBULL-P 3.5d0 1.3d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 7.112257200923508d-4)
-                        (MULTIPLE-VALUE-LIST (WEIBULL-Q 3.5d0 1.3d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
-                        (LIST 3.5000000000000164d0)
-                        (MULTIPLE-VALUE-LIST
-                         (WEIBULL-PINV 0.9992887742799077d0 1.3d0 2.0d0)))
-                       (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 3.5d0)
-                                                          (MULTIPLE-VALUE-LIST
-                                                           (WEIBULL-QINV
-                                                            7.112257200923508d-4
-                                                            1.3d0 2.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.0160712303786595d0 1.347055359960668d0
+	  1.1241262408795445d0 0.2329031397826649d0
+	  1.209338423856536d0 0.8506825453443836d0
+	  0.20845532973957762d0 0.5434187344070215d0
+	  0.7849237503774361d0 0.5487883320132739d0
+	  0.5239377808419179d0))
+   (MULTIPLE-VALUE-LIST
+    (LET ((RNG (MAKE-RANDOM-NUMBER-GENERATOR +MT19937+ 0)))
+      (LOOP FOR I FROM 0 TO 10 COLLECT
+	   (sample rng 'weibull :a 1.0d0 :b 2.0d0)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.24263174972226745d0)
+   (MULTIPLE-VALUE-LIST (WEIBULL-PDF 1.5d0 1.3d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9992887742799077d0)
+   (MULTIPLE-VALUE-LIST (WEIBULL-P 3.5d0 1.3d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 7.112257200923508d-4)
+   (MULTIPLE-VALUE-LIST (WEIBULL-Q 3.5d0 1.3d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.5000000000000164d0)
+   (MULTIPLE-VALUE-LIST
+    (WEIBULL-PINV 0.9992887742799077d0 1.3d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.5d0)
+   (MULTIPLE-VALUE-LIST (WEIBULL-QINV 7.112257200923508d-4 1.3d0 2.0d0))))