From 95e9f944af1c807524d33c245ac1e679d0b1d34b Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Mon, 9 Oct 2006 02:22:49 +0000
Subject: [PATCH] New random distributions.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3148 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
---
 gsll.asd                       |  13 ++++-
 random/beta.lisp               |  81 ++++++++++++++++++++++++++
 random/chi-squared.lisp        |  81 ++++++++++++++++++++++++++
 random/fdist.lisp              |  89 ++++++++++++++++++++++++++++
 random/flat.lisp               |  81 ++++++++++++++++++++++++++
 random/gamma.lisp              | 102 +++++++++++++++++++++++++++++++++
 random/logistic.lisp           |  80 ++++++++++++++++++++++++++
 random/lognormal.lisp          |  85 +++++++++++++++++++++++++++
 random/pareto.lisp             |  81 ++++++++++++++++++++++++++
 random/tdist.lisp              |  82 ++++++++++++++++++++++++++
 special-functions/coulomb.lisp |   4 +-
 11 files changed, 775 insertions(+), 4 deletions(-)
 create mode 100644 random/beta.lisp
 create mode 100644 random/chi-squared.lisp
 create mode 100644 random/fdist.lisp
 create mode 100644 random/flat.lisp
 create mode 100644 random/gamma.lisp
 create mode 100644 random/logistic.lisp
 create mode 100644 random/lognormal.lisp
 create mode 100644 random/pareto.lisp
 create mode 100644 random/tdist.lisp

diff --git a/gsll.asd b/gsll.asd
index 98743f16..8ecaefff 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Sat Sep 30 2006 - 22:46
+; modified:    Sun Oct  8 2006 - 21:46
 ;********************************************************
 ;;; $Id: $
 
@@ -104,4 +104,13 @@
 	     (:file "rayleigh" :depends-on (rng-types))
 	     (:file "rayleigh-tail" :depends-on (rng-types))
 	     (:file "landau" :depends-on (rng-types))
-	     (:file "levy" :depends-on (rng-types))))))
+	     (:file "levy" :depends-on (rng-types))
+	     (:file "gamma" :depends-on (rng-types))
+	     (:file "flat" :depends-on (rng-types))
+	     (:file "lognormal" :depends-on (rng-types))
+	     (:file "chi-squared" :depends-on (rng-types))
+	     (:file "fdist" :depends-on (rng-types))
+	     (:file "tdist" :depends-on (rng-types))
+	     (:file "beta" :depends-on (rng-types))
+	     (:file "logistic" :depends-on (rng-types))
+	     (:file "pareto" :depends-on (rng-types))))))
diff --git a/random/beta.lisp b/random/beta.lisp
new file mode 100644
index 00000000..8d3defaf
--- /dev/null
+++ b/random/beta.lisp
@@ -0,0 +1,81 @@
+;********************************************************
+; file:        beta.lisp                          
+; description: Beta distribution                  
+; date:        Sat Sep 30 2006
+; author:      Liam M. Healy                             
+; modified:    Sun Oct  8 2006 - 17:05
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl beta-rd (generator a b)
+  ;; Named #'beta-rd to avoid confusion with the special function #'beta.
+  "gsl_ran_beta"
+  (((generator generator) :pointer) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "A random variate from the beta distribution.  The distribution function is
+   p(x) dx = {\Gamma(a+b) \over \Gamma(a) \Gamma(b)} x^{a-1} (1-x)^{b-1} dx
+   @math{0 <= x <= 1}.")
+
+(defun-gsl beta-pdf (x a b)
+  "gsl_ran_beta_pdf" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a beta distribution with parameters @var{a} and @var{b}, using the
+   formula given in #'beta.")
+
+(defun-gsl beta-P (x a b)
+  "gsl_cdf_beta_P" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+
+(defun-gsl beta-Q (x a b)
+  "gsl_cdf_beta_Q" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+
+(defun-gsl beta-Pinv (P a b)
+  "gsl_cdf_beta_Pinv" ((P :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+
+(defun-gsl beta-Qinv (Q a b)
+  "gsl_cdf_beta_Qinv" ((Q :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the beta distribution with parameters @var{a} and @var{b}.")
+
+;;; Examples and unit test
+(lisp-unit:define-test beta
+  (lisp-unit:assert-equal
+   '("0.839000941902d-04" "0.242116468139d-01" "0.455077134726d-01"
+     "0.303211445340d+00" "0.569357215111d+00" "0.514651520667d+00"
+     "0.230096194773d+00" "0.392834882565d+00" "0.514387412254d+00"
+     "0.233783685805d+00" "0.198512886686d+00")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (beta-rd *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.180000000000d+01"
+   (beta-pdf 0.1d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.190000000000d+00"
+   (beta-P 0.1d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.810000000000d+00"
+   (beta-Q 0.1d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "1.000000000000d-01"
+   (beta-Pinv 0.19d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "1.000000000000d-01"
+   (beta-Qinv 0.81d0 1.0d0 2.0d0)))
diff --git a/random/chi-squared.lisp b/random/chi-squared.lisp
new file mode 100644
index 00000000..7bc51cb9
--- /dev/null
+++ b/random/chi-squared.lisp
@@ -0,0 +1,81 @@
+;********************************************************
+; file:        chi-squared.lisp                          
+; description: Chi-squared distribution                  
+; date:        Sat Oct  7 2006 - 16:13
+; author:      Liam M. Healy                             
+; modified:    Sun Oct  8 2006 - 16:41
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl chi-squared (generator nu)
+  "gsl_ran_chisq"
+  (((generator generator) :pointer) (nu :double))
+  :c-return :double
+  :documentation
+  "A random variate from the chi-squared distribution
+  with @var{nu} degrees of freedom. The distribution function is
+  p(x) dx = {1 \over 2 \Gamma(\nu/2) } (x/2)^{\nu/2 - 1} \exp(-x/2) dx
+  @math{x >= 0}. ")
+
+(defun-gsl chi-squared-pdf (x nu)
+  "gsl_ran_chisq_pdf" ((x :double) (nu :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a chi-squared distribution with @var{nu} degrees of freedom, using
+   the formula given in #'chi-squared.")
+
+(defun-gsl chi-squared-P (x nu)
+  "gsl_cdf_chisq_P" ((x :double) (nu :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+
+(defun-gsl chi-squared-Q (x nu)
+  "gsl_cdf_chisq_Q" ((x :double) (nu :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+
+(defun-gsl chi-squared-Pinv (P nu)
+  "gsl_cdf_chisq_Pinv" ((P :double) (nu :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+
+(defun-gsl chi-squared-Qinv (Q nu)
+  "gsl_cdf_chisq_Qinv" ((Q :double) (nu :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the chi-squared distribution with @var{nu} degrees of freedom.")
+
+;;; Examples and unit test
+(lisp-unit:define-test chi-squared
+  (lisp-unit:assert-equal
+   '("0.919043905723d+01" "0.395938453123d+01" "0.543493338508d+01"
+     "0.130231074199d+02" "0.136706908417d+02" "0.698235971535d+01"
+     "0.137566397216d+02" "0.113594580454d+02" "0.160682825906d+02"
+     "0.129995142525d+02" "0.111473474244d+02")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (chi-squared *rng-mt19937* 10.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.439391289468d+00"
+   (chi-squared-pdf 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.520499877813d+00"
+   (chi-squared-P 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.479500122187d+00"
+   (chi-squared-Q 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.500000000000d+00"
+   (chi-squared-Pinv 0.5204998778130463d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.500000000000d+00"
+   (chi-squared-Qinv 0.4795001221869537d0 1.0d0)))
diff --git a/random/fdist.lisp b/random/fdist.lisp
new file mode 100644
index 00000000..979b6988
--- /dev/null
+++ b/random/fdist.lisp
@@ -0,0 +1,89 @@
+;********************************************************
+; file:        fdist.lisp                          
+; description: Fdist distribution                  
+; date:        Sat Sep 30 2006
+; author:      Liam M. Healy                             
+; modified:    Sun Oct  8 2006 - 16:31
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl fdist (generator nu1 nu2)
+  "gsl_ran_fdist"
+  (((generator generator) :pointer) (nu1 :double) (nu2 :double))
+  :c-return :double
+  :documentation
+  "A random variate from the F-distribution with degrees of freedom @var{nu1}
+   and @var{nu2}. The distribution function is
+   p(x) dx = 
+   { \Gamma((\nu_1 + \nu_2)/2)
+        \over \Gamma(\nu_1/2) \Gamma(\nu_2/2) } 
+   \nu_1^{\nu_1/2} \nu_2^{\nu_2/2} 
+   x^{\nu_1/2 - 1} (\nu_2 + \nu_1 x)^{-\nu_1/2 -\nu_2/2}
+   @math{x >= 0}.")
+
+(defun-gsl fdist-pdf (x nu1 nu2)
+  "gsl_ran_fdist_pdf" ((x :double) (nu1 :double) (nu2 :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for an F-distribution with @var{nu1} and @var{nu2} degrees of freedom,
+   using the formula given #'fdist.")
+
+(defun-gsl fdist-P (x nu1 nu2)
+  "gsl_cdf_fdist_P" ((x :double) (nu1 :double) (nu2 :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the fdist distribution with
+  @var{nu1} and @var{nu2} degrees of freedom.")
+
+(defun-gsl fdist-Q (x nu1 nu2)
+  "gsl_cdf_fdist_Q" ((x :double) (nu1 :double) (nu2 :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the fdist distribution with
+  @var{nu1} and @var{nu2} degrees of freedom.")
+
+(defun-gsl fdist-Pinv (P nu1 nu2)
+  "gsl_cdf_fdist_Pinv" ((P :double) (nu1 :double) (nu2 :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the fdist distribution with
+  @var{nu1} and @var{nu2} degrees of freedom.")
+
+(defun-gsl fdist-Qinv (Q nu1 nu2)
+  "gsl_cdf_fdist_Qinv" ((Q :double) (nu1 :double) (nu2 :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the fdist distribution with
+   @var{nu1} and @var{nu2} degrees of freedom.")
+
+;;; Examples and unit test
+(lisp-unit:define-test fdist
+  (lisp-unit:assert-equal
+   '("0.103774233365d+03" "0.212485013221d+01" "0.304410694709d+00"
+     "0.300188687388d+00" "0.112282068448d-02" "0.292109400785d+00"
+     "0.635729092565d-01" "0.477966365217d+00" "0.347211676079d-01"
+     "0.486974823041d+00" "0.253179451696d+01")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (fdist *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.159471988462d+00"
+   (fdist-pdf 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.612372435696d+00"
+   (fdist-P 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.387627564304d+00"
+   (fdist-Q 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.120000000000d+01"
+   (fdist-Pinv 0.612372435695795d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.120000000000d+01"
+   (fdist-Qinv 0.38762756430420503d0 1.0d0 2.0d0)))
diff --git a/random/flat.lisp b/random/flat.lisp
new file mode 100644
index 00000000..d19686ad
--- /dev/null
+++ b/random/flat.lisp
@@ -0,0 +1,81 @@
+;********************************************************
+; file:        flat.lisp                          
+; description: Flat distribution                  
+; date:        Oct  7 2006
+; author:      Liam M. Healy                             
+; modified:    Sat Oct  7 2006 - 16:13
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl flat (generator a b)
+  "gsl_ran_flat"
+  (((generator generator) :pointer) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "A random variate from the flat (uniform)
+   distribution from @var{a} to @var{b}. The distribution is,
+   p(x) dx = {1 \over (b-a)} dx
+   if @c{$a \le x < b$} @math{a <= x < b} and 0 otherwise.")
+
+(defun-gsl flat-pdf (x a b)
+  "gsl_ran_flat_pdf" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a uniform distribution from @var{a} to @var{b}, using the formula
+   given for #'flat.")
+
+(defun-gsl flat-P (x a b)
+  "gsl_cdf_flat_P" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for a uniform distribution from @var{a} to @var{b}.")
+
+(defun-gsl flat-Q (x a b)
+  "gsl_cdf_flat_Q" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for a uniform distribution from @var{a} to @var{b}.")
+
+(defun-gsl flat-Pinv (P a b)
+  "gsl_cdf_flat_Pinv" ((P :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for a uniform distribution from @var{a} to @var{b}.")
+
+(defun-gsl flat-Qinv (Q a b)
+  "gsl_cdf_flat_Qinv" ((Q :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for a uniform distribution from @var{a} to @var{b}.")
+
+;;; Examples and unit test
+(lisp-unit:define-test flat
+  (lisp-unit:assert-equal
+   '("0.199974174891d+01" "0.116290987539d+01" "0.128261780529d+01"
+     "0.194720108202d+01" "0.123165654275d+01" "0.148497361434d+01"
+     "0.195747695654d+01" "0.174430534313d+01" "0.154004365834d+01"
+     "0.173995298147d+01" "0.175994379818d+01")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (flat *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.100000000000d+01"
+   (flat-pdf 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.200000000000d+00"
+   (flat-P 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.800000000000d+00"
+   (flat-Q 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.120000000000d+01"
+   (flat-Pinv 0.19999999999999996d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.120000000000d+01"
+   (flat-Qinv 0.8d0 1.0d0 2.0d0)))
diff --git a/random/gamma.lisp b/random/gamma.lisp
new file mode 100644
index 00000000..a439c6a6
--- /dev/null
+++ b/random/gamma.lisp
@@ -0,0 +1,102 @@
+;********************************************************
+; file:        gamma.lisp                          
+; description: Gamma distribution                  
+; date:        Sat Sep 30 2006
+; author:      Liam M. Healy                             
+; modified:    Sat Oct  7 2006 - 11:53
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl gamma-rd (generator a b)
+  ;; Named #'gamma-rd to avoid confusion with the special function #'gamma.
+  "gsl_ran_gamma"
+  (((generator generator) :pointer) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "A random variate from the gamma distribution.
+   The distribution function is
+   p(x) dx = {1 \over \Gamma(a) b^a} x^{a-1} e^{-x/b} dx
+   for @math{x > 0}. The gamma distribution with an integer parameter @var{a}
+   is known as the Erlang distribution.  The variates are computed using
+   the algorithms from Knuth (vol 2).")
+
+(defun-gsl gamma-mt (generator a b)
+  "gsl_ran_gamma_mt"
+  (((generator generator) :pointer) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "A gamma variate using the Marsaglia-Tsang fast gamma method.")
+
+(defun-gsl gamma-pdf (x a b)
+  "gsl_ran_gamma_pdf" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a gamma distribution with parameters @var{a} and @var{b}, using the
+   formula given in #'gamma.")
+
+(defun-gsl gamma-P (x a b)
+  "gsl_cdf_gamma_P" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+
+(defun-gsl gamma-Q (x a b)
+  "gsl_cdf_gamma_Q" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+
+(defun-gsl gamma-Pinv (P a b)
+  "gsl_cdf_gamma_Pinv" ((P :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+
+(defun-gsl gamma-Qinv (Q a b)
+  "gsl_cdf_gamma_Qinv" ((Q :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the Gamma distribution with parameters @var{a} and @var{b}.")
+
+;;; Examples and unit test
+(lisp-unit:define-test gamma-randist
+  (lisp-unit:assert-equal
+   '("0.516568891768d-03" "0.362911628560d+01" "0.252731961087d+01"
+     "0.108487745041d+00" "0.292499884683d+01" "0.144732158591d+01"
+     "0.869072489937d-01" "0.590607841809d+00" "0.123221058781d+01"
+     "0.602337266708d+00" "0.549021596387d+00")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (gamma-rd *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-equal
+   '("0.260013787613d+01" "0.252226669542d+01" "0.773142209213d+01"
+     "0.422727992649d+01" "0.951930434749d-01" "0.571092010687d+00"
+     "0.891063771946d+00" "0.826322120255d+00" "0.318306657206d+01"
+     "0.380840036132d-02" "0.103201173341d+01")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (gamma-mt *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.475614712250d+00"
+   (gamma-pdf 0.1d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.487705754993d-01"
+   (gamma-P 0.1d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.951229424501d+00"
+   (gamma-Q 0.1d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.100000000000d+00"
+   (gamma-Pinv 0.048770575499286005d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.100000000000d+00"
+   (gamma-Qinv 0.951229424500714d0 1.0d0 2.0d0)))
diff --git a/random/logistic.lisp b/random/logistic.lisp
new file mode 100644
index 00000000..53b31266
--- /dev/null
+++ b/random/logistic.lisp
@@ -0,0 +1,80 @@
+;********************************************************
+; file:        logistic.lisp                          
+; description: logistic distribution                  
+; date:        Sat Oct  7 2006 - 16:13
+; author:      Liam M. Healy                             
+; modified:    Sun Oct  8 2006 - 17:32
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl logistic (generator a)
+  "gsl_ran_logistic"
+  (((generator generator) :pointer) (a :double))
+  :c-return :double
+  :documentation
+  "A random variate from the logistic distribution.  The distribution function is
+   p(x) dx = { \exp(-x/a) \over a (1 + \exp(-x/a))^2 } dx
+   for @math{-\infty < x < +\infty}.")
+
+(defun-gsl logistic-pdf (x a)
+  "gsl_ran_logistic_pdf" ((x :double) (a :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a logistic distribution with scale parameter @var{a}, using the
+   formula given in #'logistic.")
+
+(defun-gsl logistic-P (x a)
+  "gsl_cdf_logistic_P" ((x :double) (a :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the logistic distribution with scale parameter @var{a}.")
+
+(defun-gsl logistic-Q (x a)
+  "gsl_cdf_logistic_Q" ((x :double) (a :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the logistic distribution with scale parameter @var{a}.")
+
+(defun-gsl logistic-Pinv (P a)
+  "gsl_cdf_logistic_Pinv" ((P :double) (a :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the logistic distribution with scale parameter @var{a}.")
+
+(defun-gsl logistic-Qinv (Q a)
+  "gsl_cdf_logistic_Qinv" ((Q :double) (a :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the logistic distribution with scale parameter @var{a}.")
+
+;;; Examples and unit test
+(lisp-unit:define-test logistic
+  (lisp-unit:assert-equal
+   '("0.826131993192d+02" "-0.163673460427d+02" "-0.931513272044d+01"
+     "0.288702070871d+02" "-0.119898098758d+02" "-0.601236476200d+00"
+     "0.311425552636d+02" "0.106846737210d+02" "0.160518409540d+01"
+     "0.104572419047d+02" "0.115237141063d+02")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (logistic *rng-mt19937* 10.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.235003712202d+00"
+   (logistic-pdf 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.622459331202d+00"
+   (logistic-P 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.377540668798d+00"
+   (logistic-Q 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.500000000000d+00"
+   (logistic-Pinv 0.6224593312018546d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.500000000000d+00"
+   (logistic-Qinv 0.37754066879814546d0 1.0d0)))
diff --git a/random/lognormal.lisp b/random/lognormal.lisp
new file mode 100644
index 00000000..b23c9796
--- /dev/null
+++ b/random/lognormal.lisp
@@ -0,0 +1,85 @@
+;********************************************************
+; file:        lognormal.lisp                          
+; description: Lognormal distribution                  
+; date:        Sat Sep 30 2006
+; author:      Liam M. Healy                             
+; modified:    Sat Oct  7 2006 - 15:18
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl lognormal (generator zeta sigma)
+  "gsl_ran_lognormal"
+  (((generator generator) :pointer) (zeta :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "A random variate from the lognormal distribution.
+   The distribution function is
+   p(x) dx = {1 \over x \sqrt{2 \pi \sigma^2}} \exp(-(\ln(x) - \zeta)^2/2 \sigma^2) dx
+   for @math{x > 0}.")
+
+(defun-gsl lognormal-pdf (x zeta sigma)
+  "gsl_ran_lognormal_pdf" ((x :double) (zeta :double) (sigma :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a lognormal distribution with parameters @var{zeta} and @var{sigma},
+   using the formula given in #'lognormal.")
+
+(defun-gsl lognormal-P (x zeta sigma)
+  "gsl_cdf_lognormal_P" ((x :double) (zeta :double) (sigma :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the lognormal distribution with parameters
+  @var{zeta} and @var{sigma}.")
+
+(defun-gsl lognormal-Q (x zeta sigma)
+  "gsl_cdf_lognormal_Q" ((x :double) (zeta :double) (sigma :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the lognormal distribution with parameters
+  @var{zeta} and @var{sigma}.")
+
+(defun-gsl lognormal-Pinv (P zeta sigma)
+  "gsl_cdf_lognormal_Pinv" ((P :double) (zeta :double) (sigma :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the lognormal distribution with parameters
+  @var{zeta} and @var{sigma}.")
+
+(defun-gsl lognormal-Qinv (Q zeta sigma)
+  "gsl_cdf_lognormal_Qinv" ((Q :double) (zeta :double) (sigma :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the lognormal distribution with parameters
+   @var{zeta} and @var{sigma}.")
+
+;;; Examples and unit test
+(lisp-unit:define-test lognormal
+  (lisp-unit:assert-equal
+   '("0.238644706813d+01" "0.116876021674d+00" "0.475337457880d+01"
+     "0.300933937758d+02" "0.811958437576d+00" "0.316342105516d+01"
+     "0.914620656772d+00" "0.727307901066d+00" "0.218018485218d+01"
+     "0.389088566169d+01" "0.182184697889d+03")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (lognormal *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.152898339657d+00"
+   (lognormal-pdf 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.341328827235d+00"
+   (lognormal-P 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.658671172765d+00"
+   (lognormal-Q 1.2d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.120000000000d+01"
+   (lognormal-Pinv 0.3413288272347352d0 1.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.120000000000d+01"
+   (lognormal-Qinv 0.6586711727652649d0 1.0d0 2.0d0)))
diff --git a/random/pareto.lisp b/random/pareto.lisp
new file mode 100644
index 00000000..465b126c
--- /dev/null
+++ b/random/pareto.lisp
@@ -0,0 +1,81 @@
+;********************************************************
+; file:        pareto.lisp                          
+; description: Pareto distribution                  
+; date:        Sat Oct  8 2006 - 21:23
+; author:      Liam M. Healy                             
+; modified:    Sun Oct  8 2006 - 21:45
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl pareto (generator a b)
+  "gsl_ran_pareto"
+  (((generator generator) :pointer) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "A random variate from the Pareto distribution of order @var{a}.
+   The distribution function is
+   p(x) dx = (a/b) / (x/b)^{a+1} dx
+   @math{x >= b}.")
+
+(defun-gsl pareto-pdf (x a b)
+  "gsl_ran_pareto_pdf" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a Pareto distribution with exponent @var{a} and scale @var{b}, using
+   the formula given in #'pareto.")
+
+(defun-gsl pareto-P (x a b)
+  "gsl_cdf_pareto_P" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+
+(defun-gsl pareto-Q (x a b)
+  "gsl_cdf_pareto_Q" ((x :double) (a :double) (b :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+
+(defun-gsl pareto-Pinv (P a b)
+  "gsl_cdf_pareto_Pinv" ((P :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+
+(defun-gsl pareto-Qinv (Q a b)
+  "gsl_cdf_pareto_Qinv" ((Q :double) (a :double) (b :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the Pareto distribution with exponent @var{a} and scale @var{b}.")
+
+;;; Examples and unit test
+(lisp-unit:define-test pareto
+  (lisp-unit:assert-equal
+   '("0.200051663561d+01" "0.122767265962d+02" "0.707669496594d+01"
+     "0.211148407447d+01" "0.863347081110d+01" "0.412393569645d+01"
+     "0.208882311615d+01" "0.268706924980d+01" "0.370340428797d+01"
+     "0.270287443943d+01" "0.263177356639d+01")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (pareto *rng-mt19937* 1.0d0 2.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.511603440571d+00"
+   (pareto-pdf 1.5d0 1.3d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.516884983518d+00"
+   (pareto-P 3.5d0 1.3d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.483115016482d+00"
+   (pareto-Q 3.5d0 1.3d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.350000000000d+01"
+   (pareto-Pinv 0.5168849835182453d0 1.3d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.350000000000d+01"
+   (pareto-Qinv 0.4831150164817547d0 1.3d0 2.0d0)))
diff --git a/random/tdist.lisp b/random/tdist.lisp
new file mode 100644
index 00000000..667f8f7c
--- /dev/null
+++ b/random/tdist.lisp
@@ -0,0 +1,82 @@
+;********************************************************
+; file:        tdist.lisp                          
+; description: tdist distribution                  
+; date:        Sat Oct  7 2006 - 16:13
+; author:      Liam M. Healy                             
+; modified:    Sun Oct  8 2006 - 16:44
+;********************************************************
+;;; $Id: $
+
+(in-package :gsl)
+
+(defun-gsl tdist (generator nu)
+  "gsl_ran_tdist"
+  (((generator generator) :pointer) (nu :double))
+  :c-return :double
+  :documentation
+  "A random variate from the Student t-distribution.  The
+   distribution function is,
+   p(x) dx = {\Gamma((\nu + 1)/2) \over \sqrt{\pi \nu} \Gamma(\nu/2)}
+   (1 + x^2/\nu)^{-(\nu + 1)/2} dx
+  for @math{-\infty < x < +\infty}.")
+
+(defun-gsl tdist-pdf (x nu)
+  "gsl_ran_tdist_pdf" ((x :double) (nu :double))
+  :c-return :double
+  :documentation
+  "The probability density @math{p(x)} at @var{x}
+   for a t-distribution with @var{nu} degrees of freedom, using the formula
+   given in #'tdist.")
+
+(defun-gsl tdist-P (x nu)
+  "gsl_cdf_tdist_P" ((x :double) (nu :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{P(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+
+(defun-gsl tdist-Q (x nu)
+  "gsl_cdf_tdist_Q" ((x :double) (nu :double))
+  :c-return :double
+  :documentation "The cumulative distribution functions
+  @math{Q(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+
+(defun-gsl tdist-Pinv (P nu)
+  "gsl_cdf_tdist_Pinv" ((P :double) (nu :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+  @math{P(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+
+(defun-gsl tdist-Qinv (Q nu)
+  "gsl_cdf_tdist_Qinv" ((Q :double) (nu :double))
+  :c-return :double
+  :documentation  "The inverse cumulative distribution functions
+   @math{Q(x)} for the tdist distribution with @var{nu} degrees of freedom.")
+
+;;; Examples and unit test
+(lisp-unit:define-test tdist
+  (lisp-unit:assert-equal
+   '("0.149893663745d+00" "0.679414287929d+00" "-0.161583395111d+01"
+     "-0.160088628258d+01" "-0.170109355058d+01" "-0.437095974981d-01"
+     "0.127611592766d+00" "-0.197312182555d-01" "-0.653466611720d+00"
+     "0.203577132452d+00" "0.177650300478d+01")
+   (lisp-unit:fp-sequence
+    (progn
+      (rng-set *rng-mt19937* 0)
+      (loop for i from 0 to 10
+	    collect
+	    (tdist *rng-mt19937* 10.0d0)))))
+  (lisp-unit:assert-first-fp-equal
+   "0.254647908947d+00"
+   (tdist-pdf 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.647583617650d+00"
+   (tdist-P 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.352416382350d+00"
+   (tdist-Q 0.5d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.500000000000d+00"
+   (tdist-Pinv 0.6475836176504334d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.500000000000d+00"
+   (tdist-Qinv 0.3524163823495667d0 1.0d0)))
diff --git a/special-functions/coulomb.lisp b/special-functions/coulomb.lisp
index 5b9c7617..17b25c7a 100644
--- a/special-functions/coulomb.lisp
+++ b/special-functions/coulomb.lisp
@@ -3,7 +3,7 @@
 ; description: Coulumb functions                         
 ; date:        Sat Mar 18 2006 - 23:23                   
 ; author:      Liam M. Healy                             
-; modified:    Mon Jun 12 2006 - 08:55
+; modified:    Sat Oct  7 2006 - 15:21
 ;********************************************************
 ;;; $Id:$
 
@@ -156,7 +156,7 @@
    (coulomb-wave-fg 1.0d0 2.0d0 2.5d0 1))
   (lisp-unit:assert-equal
    '("0.350215846039d-01" "0.575250061420d-02" "0.711695560198d-03"
-     "0.000000000000d+01" "0.000000000000d+01" "0.000000000000d+01")
+     "0.000000000000d+01" "0.800000000000d+01" "0.110000000000d+02")
    (lisp-unit:fp-sequence 
     (with-data (Farr vector-double 3)
       (with-data (Garr vector-double 3)
-- 
GitLab