diff --git a/gsll.asd b/gsll.asd
index 4fb9cf221da5d470d79f391e02b15fe476bcecd7..f94c90e959f2050fc717c988ca70fe8a8bf0450e 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Mon Jun 12 2006 - 23:43
+; modified:    Tue Jun 13 2006 - 22:58
 ;********************************************************
 ;;; $Id: $
 
@@ -38,10 +38,8 @@
 	     (:file "matrix" :depends-on (data vector))
 	     (:file "permutation" :depends-on (data vector))
 	     (:file "combination" :depends-on (data))))
-   #+future
-   (:file "cffi-array")
-   #+future
-   (:file "polynomial" :depends-on (init cffi-array)) ; see file
+   #+unnecessary (:file "cffi-array")
+   (:file "polynomial" :depends-on (init data))
    (:module special-functions
 	    :depends-on (init)
 	    :components
@@ -53,15 +51,15 @@
 	     (:file "dawson")
 	     (:file "debye")
 	     (:file "dilogarithm")
-	     ;;(:file "elementary")
-	     ;;(:file "elliptic-integrals")
-	     ;;(:file "elliptic-functions")
-	     ;;(:file "error-functions")
-	     ;;(:file "exponential-functions")
-	     ;;(:file "exponential-integrals")
-	     ;;(:file "fermi-dirac")
-	     ;;(:file "gamma")
-	     ;;(:file "gegenbauer")
+	     (:file "elementary")
+	     (:file "elliptic-integrals")
+	     (:file "elliptic-functions")
+	     (:file "error-functions")
+	     (:file "exponential-functions")
+	     (:file "exponential-integrals")
+	     (:file "fermi-dirac")
+	     (:file "gamma")
+	     (:file "gegenbauer")
 	     ;;(:file "hypergeometric")
 	     ;;(:file "laguerre")
 	     ;;(:file "lambert")
diff --git a/init/structures.lisp b/init/structures.lisp
index 4710cfc2bc788daffd8dcbc7699cc3741b58f0e7..0e0028aad377b2e136586df676a1ba876513c43a 100644
--- a/init/structures.lisp
+++ b/init/structures.lisp
@@ -3,7 +3,7 @@
 ; description: Common structures defined by GSL.         
 ; date:        Sun May 28 2006 - 22:04                   
 ; author:      Liam M. Healy                             
-; modified:    Thu Jun  1 2006 - 14:32
+; modified:    Tue Jun 13 2006 - 21:45
 ;********************************************************
 ;;; $Id: $
 
@@ -41,7 +41,7 @@ and a scaling exponent e10, such that the value is val*10^e10."
   (cffi:foreign-slot-value sf-result type 'err))
 
 (defun e10 (sf-result)
-  (cffi:foreign-slot-value sf-result 'sf-result-e10 'err))
+  (cffi:foreign-slot-value sf-result 'sf-result-e10 'e10))
 
 ;;;;****************************************************************************
 ;;;; Complex numbers
diff --git a/polynomial.lisp b/polynomial.lisp
index a0d5cd7761f41a26bfa5fde67adc3e2371c34ade..9136a7a7e22107c0a2bda5d76732d67ab35785d4 100644
--- a/polynomial.lisp
+++ b/polynomial.lisp
@@ -3,15 +3,14 @@
 ; description: Polynomials                               
 ; date:        Tue Mar 21 2006 - 18:33                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Jun 10 2006 - 23:37
+; modified:    Tue Jun 13 2006 - 22:57
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-;;; To do: finish divided differences, which requires figuring out how
-;;; to handle raw C arrays, and deciding if/how to provide
-;;; autotranslation from CL pure arrays.
+;;; Provide autotranslation from CL pure arrays?
+;;; Divided differences not complete/tested.
 
 ;;;;****************************************************************************
 ;;;; Polynomial Evaluation
@@ -28,10 +27,6 @@
 ;;;; Divided Difference Representation of Polynomials
 ;;;;****************************************************************************
 
-;;; Use with-divided-difference to compute the divided difference,
-;;; which may be passed to eval-divided-difference or
-;;; taylor-divided-difference in the body.
-
 (defun-gsl divided-difference-int (dd xa ya)
   "gsl_poly_dd_init"
   (((gsl-array dd) :pointer)
@@ -67,19 +62,22 @@
   "Evaluate the polynomial stored in divided-difference form
    in the arrays @var{dd} and @var{xa} at the point @var{x}.")
 
-#+development
-(defun-gsl taylor-divided-difference (coefs dd xp workspace)
+(defun-gsl taylor-divided-difference (coefs xp dd xa workspace)
   "gsl_poly_dd_taylor"
   (((gsl-array coefs) :pointer)
-   ((gsl-array xp) :pointer)
+   (xp :double)
    ((gsl-array dd) :pointer)
-   ((gsl-array x) :pointer)
+   ((gsl-array xa) :pointer)
    ((dim0 xa) :size)
    ((gsl-array workspace) :pointer))
+  :invalidate (coefs)
   :documentation
-  "Convert the divided-difference representation of a polynomial
-   to a Taylor expansion about the point xp.  Call only within a
-   with-divided-difference form.")
+  "Convert the divided-difference representation of a
+  polynomial to a Taylor expansion.  The divided-difference representation
+  is supplied in the arrays @var{dd} and @var{xa} of the same length.
+  On output the Taylor coefficients of the polynomial expanded about the
+  point @var{xp} are stored in the array coefs which has the same length
+  as xa and dd.  A workspace of length @var{size} must be provided.")
 
 ;;;;****************************************************************************
 ;;;; Quadratic Equations
@@ -219,5 +217,5 @@
      ("0.309016994375d+00" "-0.951056516295d+00")
      ("0.100000000000d+01" "0.000000000000d+01"))
    ;; Example from GSL manual
-   (lisp-unit:fp-values (polynomial-solve #(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
-   ))
+   (lisp-unit:fp-values
+    (polynomial-solve #(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))))
diff --git a/special-functions/elementary.lisp b/special-functions/elementary.lisp
index 417cf2bb93f4d176103863d7b23c0400e3c89c66..17cdd46b30e13deb704722b732fb0f7ad3b8d8e6 100644
--- a/special-functions/elementary.lisp
+++ b/special-functions/elementary.lisp
@@ -3,25 +3,32 @@
 ; description: Elementary functions                      
 ; date:        Mon Mar 20 2006 - 21:43                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Mar 25 2006 - 22:11
+; modified:    Tue Jun 13 2006 - 21:04
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl multiply ((x :double) (y :double))
+(defun-gsl multiply (x y)
   "gsl_sf_multiply_e"
+  ((x :double) (y :double) (ret sf-result))
   :documentation
   "Multiplies @var{x} and @var{y} returning the product and
-associated error."
-  :return (sf-result))
+  associated error.")
 
-(defun-gsl multiply-err ((x :double) (dx :double) (y :double)  (dy :double))
-  "gsl_sf_multiply_err_e"
+(defun-gsl multiply-err (x dx y dy)
+    "gsl_sf_multiply_err_e"
+  ((x :double) (dx :double) (y :double)  (dy :double) (ret sf-result))
   :documentation
   "Multiplies @var{x} and @var{y} with associated absolute
-errors @var{dx} and @var{dy}.  The product 
-@c{$xy \pm xy \sqrt{(dx/x)^2 +(dy/y)^2}$} 
-@math{xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2)} 
-is returned."
-  :return (sf-result))
+   errors @var{dx} and @var{dy}.  The product 
+   @math{xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2)} 
+   is returned.")
+
+(lisp-unit:define-test elementary
+  (lisp-unit:assert-first-fp-equal
+   "0.600000000000d+01"
+   (multiply 3.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.600000000000d+01"
+   (multiply-err 3.0d0 0.1d0 2.0d0 0.1d0)))
diff --git a/special-functions/elliptic-functions.lisp b/special-functions/elliptic-functions.lisp
index 6c0ed1386455acc77c85b857ed148da652e7d5a7..af3a0f057c80a13088e18cbc7f5b9f438ccfd0bf 100644
--- a/special-functions/elliptic-functions.lisp
+++ b/special-functions/elliptic-functions.lisp
@@ -3,19 +3,20 @@
 ; description: Jacobian elliptic functions                        
 ; date:        Mon Mar 20 2006 - 22:21                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Mar 25 2006 - 22:11
+; modified:    Tue Jun 13 2006 - 21:16
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl jacobian-elliptic-functions ((u :double) (m :double))
+(defun-gsl jacobian-elliptic-functions (u m)
   "gsl_sf_elljac_e"
+  ((u :double) (m :double) (sn sf-result) (cn sf-result) (dn sf-result))
   :documentation
   "The Jacobian elliptic functions @math{sn(u|m)},
-@math{cn(u|m)}, @math{dn(u|m)} computed by descending Landen
-transformations."
-  :return (:double :double :double))
+  @math{cn(u|m)}, @math{dn(u|m)} computed by descending Landen
+  transformations."
+  :return ((val sn) (val cn) (val dn) (err sn) (err cn) (err dn)))
 
 ;;; > (jacobian-elliptic-functions 0.61802d0 0.5d0)
 ;;; 0.564575752943391
@@ -27,3 +28,14 @@ transformations."
 ;;; 0.9840560289645665
 ;;; > (jacobian-elliptic-functions 0.61802d0 1.5d0)
 ;;; ;;;error
+
+
+(lisp-unit:define-test elliptic-functions
+  (lisp-unit:assert-equal
+   '("0.197620823672d+00" "0.980278536974d+00" "0.984056028965d+00")
+   (subseq
+    (lisp-unit:fp-values (jacobian-elliptic-functions 0.2d0 0.81d0))
+    0 3))
+  (lisp-unit:assert-error
+   'gsl-error
+   (jacobian-elliptic-functions 0.61802d0 1.5d0)))
diff --git a/special-functions/elliptic-integrals.lisp b/special-functions/elliptic-integrals.lisp
index 4f4415453b4d180069a860ced2ee28713685459b..6efb5c3b0362a15ca60ce42c53205f779668083f 100644
--- a/special-functions/elliptic-integrals.lisp
+++ b/special-functions/elliptic-integrals.lisp
@@ -3,7 +3,7 @@
 ; description: Elliptic integrals                        
 ; date:        Mon Mar 20 2006 - 21:50                   
 ; author:      Liam M. Healy                             
-; modified:    Sun May 21 2006 - 19:06
+; modified:    Tue Jun 13 2006 - 21:09
 ;********************************************************
 ;;; $Id: $
 
@@ -13,86 +13,70 @@
 ;;;; Legendre form of complete elliptic integrals
 ;;;;****************************************************************************
 
-(defun-gsl elliptic-integral-K-complete ((k :double))
-  "gsl_sf_ellint_Kcomp_e"
+(defun-gsl elliptic-integral-K-complete (k)
+  "gsl_sf_ellint_Kcomp_e" ((k :double) :mode (ret sf-result))
   :documentation
-  "The complete elliptic integral @math{K(k)}."
-  :mode t
-  :return (sf-result))
+  "The complete elliptic integral @math{K(k)}.")
 
-(defun-gsl elliptic-integral-E-complete ((k :double))
-  "gsl_sf_ellint_Ecomp_e"
-  :documentation
-  "The complete elliptic integral @math{E(k)}."
-  :mode t
-  :return (sf-result))
+(defun-gsl elliptic-integral-E-complete (k)
+  "gsl_sf_ellint_Ecomp_e" ((k :double) :mode (ret sf-result))
+  :documentation "The complete elliptic integral @math{E(k)}.")
 
 ;;;;****************************************************************************
 ;;;; Legendre form of incomplete elliptic integrals
 ;;;;****************************************************************************
 
-(defun-gsl elliptic-integral-F ((phi :double) (k :double))
-  "gsl_sf_ellint_F_e"
+(defun-gsl elliptic-integral-F (phi k)
+  "gsl_sf_ellint_F_e" ((phi :double) (k :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{F(\phi,k)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{F(\phi,k)}.")
 
-(defun-gsl elliptic-integral-E ((phi :double) (k :double))
-  "gsl_sf_ellint_E_e"
+(defun-gsl elliptic-integral-E (phi k)
+  "gsl_sf_ellint_E_e" ((phi :double) (k :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{E(\phi,k)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{E(\phi,k)}.")
 
-(defun-gsl elliptic-integral-P ((phi :double) (k :double) (n :double))
+(defun-gsl elliptic-integral-P (phi k n)
   "gsl_sf_ellint_P_e"
+  ((phi :double) (k :double) (n :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{P(\phi,k,n)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{P(\phi,k,n)}.")
 
-(defun-gsl elliptic-integral-D ((phi :double) (k :double) (n :double))
+(defun-gsl elliptic-integral-D (phi k n)
   "gsl_sf_ellint_D_e"
+  ((phi :double) (k :double) (n :double) :mode (ret sf-result))
   :documentation
   "The incomplete elliptic integral @math{D(\phi,k,n)} which is
    defined through the Carlson form @math{RD(x,y,z)}
    by the following relation, 
-   D(\phi,k,n) = RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1)."
-  :mode t
-  :return (sf-result))
+   D(\phi,k,n) = RD (1-\sin^2(\phi), 1-k^2 \sin^2(\phi), 1).")
 
 ;;;;****************************************************************************
 ;;;; Carlson forms
 ;;;;****************************************************************************
 
-(defun-gsl elliptic-integral-RC ((x :double) (y :double))
-  "gsl_sf_ellint_RC_e"
+(defun-gsl elliptic-integral-RC (x y)
+  "gsl_sf_ellint_RC_e" ((x :double) (y :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{RC(x,y)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{RC(x,y)}.")
 
-(defun-gsl elliptic-integral-RD ((x :double) (y :double) (z :double))
+(defun-gsl elliptic-integral-RD (x y z)
   "gsl_sf_ellint_RD_e"
+  ((x :double) (y :double) (z :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{RD(x,y,z)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{RD(x,y,z)}.")
 
-(defun-gsl elliptic-integral-RF ((x :double) (y :double) (z :double))
+(defun-gsl elliptic-integral-RF (x y z)
   "gsl_sf_ellint_RF_e"
+  ((x :double) (y :double) (z :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{RF(x,y,z)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{RF(x,y,z)}.")
 
-(defun-gsl elliptic-integral-RJ ((x :double) (y :double) (z :double) (p :double))
+(defun-gsl elliptic-integral-RJ (x y z p)
   "gsl_sf_ellint_RJ_e"
+  ((x :double) (y :double) (z :double) (p :double) :mode (ret sf-result))
   :documentation
-  "The incomplete elliptic integral @math{RJ(x,y,z,p)}."
-  :mode t
-  :return (sf-result))
+  "The incomplete elliptic integral @math{RJ(x,y,z,p)}.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
diff --git a/special-functions/error-functions.lisp b/special-functions/error-functions.lisp
index 053dd68288bf2ea5d3e5c7a28638c7bd1573c922..480b01de5482ea76f64f4f4df77036ef19719fad 100644
--- a/special-functions/error-functions.lisp
+++ b/special-functions/error-functions.lisp
@@ -3,62 +3,50 @@
 ; description: Error functions                           
 ; date:        Mon Mar 20 2006 - 22:31                   
 ; author:      Liam M. Healy                             
-; modified:    Sun May 21 2006 - 19:08
+; modified:    Tue Jun 13 2006 - 21:20
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl erf ((x :double))
-  "gsl_sf_erf_e"
+(defun-gsl erf (x)
+  "gsl_sf_erf_e" ((x :double) (ret sf-result))
   :documentation
-  "The error function @c{$\erf(x)$}
-@math{erf(x)}, where
-@c{$\erf(x) = (2/\sqrt{\pi}) \int_0^x dt \exp(-t^2)$}
-@math{erf(x) = (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2)}."
-  :return (sf-result))
+  "The error function @math{erf(x)}, where
+  @math{erf(x) = (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2)}.")
 
-(defun-gsl erfc ((x :double))
-  "gsl_sf_erfc_e"
+(defun-gsl erfc (x)
+  "gsl_sf_erfc_e" ((x :double) (ret sf-result))
   :documentation
   "The complementary error function 
-@c{$\erfc(x) = 1 - \erf(x) = (2/\sqrt{\pi}) \int_x^\infty \exp(-t^2)$}
-@math{erfc(x) = 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2)}."
-  :return (sf-result))
+  @math{erfc(x) = 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2)}.")
 
-(defun-gsl log-erfc ((x :double))
-  "gsl_sf_log_erfc_e"
+(defun-gsl log-erfc (x)
+  "gsl_sf_log_erfc_e" ((x :double) (ret sf-result))
   :documentation
-  "The logarithm of the complementary error function @math{\log(\erfc(x))}."
-  :return (sf-result))
+  "The logarithm of the complementary error function @math{\log(\erfc(x))}.")
 
-(defun-gsl erf-Z ((x :double))
-  "gsl_sf_erf_Z_e"
+(defun-gsl erf-Z (x)
+  "gsl_sf_erf_Z_e" ((x :double) (ret sf-result))
   :documentation
   "The Gaussian probability density function 
-@c{$Z(x) = (1/\sqrt{2\pi}) \exp(-x^2/2)$} 
-@math{Z(x) = (1/\sqrt@{2\pi@}) \exp(-x^2/2)}."
-  :return (sf-result))
+  @math{Z(x) = (1/\sqrt@{2\pi@}) \exp(-x^2/2)}.")
 
-(defun-gsl erf-Q ((x :double))
-  "gsl_sf_erf_Q_e"
+(defun-gsl erf-Q (x)
+  "gsl_sf_erf_Q_e" ((x :double) (ret sf-result))
   :documentation
-"The upper tail of the Gaussian probability
-function 
-@c{$Q(x) = (1/\sqrt{2\pi}) \int_x^\infty dt \exp(-t^2/2)$}
-@math{Q(x) = (1/\sqrt@{2\pi@}) \int_x^\infty dt \exp(-t^2/2)}."
-  :return (sf-result))
+  "The upper tail of the Gaussian probability function 
+  @math{Q(x) = (1/\sqrt@{2\pi@}) \int_x^\infty dt \exp(-t^2/2)}.")
 
-(defun-gsl hazard ((x :double))
-  "gsl_sf_hazard_e"
+(defun-gsl hazard (x)
+  "gsl_sf_hazard_e" ((x :double) (ret sf-result))
   :documentation
-  "The hazard function for the normal distribution."
-  :return (sf-result))
+  "The hazard function for the normal distribution.")
 
 (lisp-unit:define-test error-functions
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.842700792950d+00" (ERF 1.0d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.157299207050d+00" (ERFC 1.0d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "-0.184960550993d+01" (LOG-ERFC 1.0d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.241970724519d+00" (ERF-Z 1.0d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.158655253931d+00" (ERF-Q 1.0d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL "0.152513527616d+01" (HAZARD 1.0d0)))
+  (lisp-unit:assert-first-fp-equal "0.842700792950d+00" (erf 1.0d0))
+  (lisp-unit:assert-first-fp-equal "0.157299207050d+00" (erfc 1.0d0))
+  (lisp-unit:assert-first-fp-equal "-0.184960550993d+01" (log-erfc 1.0d0))
+  (lisp-unit:assert-first-fp-equal "0.241970724519d+00" (erf-z 1.0d0))
+  (lisp-unit:assert-first-fp-equal "0.158655253931d+00" (erf-q 1.0d0))
+  (lisp-unit:assert-first-fp-equal "0.152513527616d+01" (hazard 1.0d0)))
diff --git a/special-functions/exponential-functions.lisp b/special-functions/exponential-functions.lisp
index ec6d589f921994350f3382825fe94340373817f4..81ea52ae116a98f07ee565f85a20ba9109628e78 100644
--- a/special-functions/exponential-functions.lisp
+++ b/special-functions/exponential-functions.lisp
@@ -3,7 +3,7 @@
 ; description: Exponential functions                     
 ; date:        Tue Mar 21 2006 - 17:05                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Apr 29 2006 - 19:07
+; modified:    Tue Jun 13 2006 - 21:36
 ;********************************************************
 ;;; $Id: $
 
@@ -13,86 +13,84 @@
 ;;;; Exponential Functions
 ;;;;****************************************************************************
 
-(defun-gsl gsl-exp ((x :double))
-  "gsl_sf_exp_e"
-  :documentation
-  "The exponential function."
-  :return (sf-result))
+(defun-gsl gsl-exp (x)
+  "gsl_sf_exp_e" ((x :double) (ret sf-result))
+  :documentation "The exponential function.")
 
-(defun-gsl exp-scaled ((x :double))
-  "gsl_sf_exp_e10_e"
+(defun-gsl exp-scaled (x)
+  "gsl_sf_exp_e10_e" ((x :double) (ret sf-result-e10))
   :documentation
   "The exponential function scaled. This function may be useful if the value
-of @math{\exp(x)} would overflow the  numeric range of @code{double}."
-  :return (sf-result-e10))
+   of @math{\exp(x)} would overflow the  numeric range of @code{double}.")
 
-(defun-gsl exp-mult ((x :double) (y :double))
-  "gsl_sf_exp_mult_e"
-  :documentation
-  "Exponentiate @var{x} and multiply by the factor @var{y} to return the product @math{y \exp(x)}."
-  :return (sf-result))
+(defun-gsl exp-mult (x y)
+  "gsl_sf_exp_mult_e" ((x :double) (y :double) (ret sf-result))
+  :documentation "Exponentiate @var{x} and multiply by the
+  factor @var{y} to return the product @math{y \exp(x)}.")
 
-(defun-gsl exp-mult-scaled ((x :double) (y :double))
-  "gsl_sf_exp_mult_e10_e"
+(defun-gsl exp-mult-scaled (x y)
+  "gsl_sf_exp_mult_e10_e" ((x :double) (y :double) (ret sf-result-e10))
   :documentation
-  "The product @math{y \exp(x)} with extended numeric range."
-  :return (sf-result-e10))
+  "The product @math{y \exp(x)} with extended numeric range.")
 
 ;;;;****************************************************************************
 ;;;; Relative Exponential Functions
 ;;;;****************************************************************************
 
-(defun-gsl expm1 ((x :double))
-  "gsl_sf_expm1_e"
+(defun-gsl expm1 (x)
+  "gsl_sf_expm1_e" ((x :double) (ret sf-result))
   :documentation
-  "@math{\exp(x)-1} using an algorithm that is accurate for small @math{x}."
-  :return (sf-result))
+  "@math{\exp(x)-1} using an algorithm that is accurate for small @math{x}.")
 
-(defun-gsl exprel ((x :double))
-  "gsl_sf_exprel_e"
+(defun-gsl exprel (x)
+  "gsl_sf_exprel_e" ((x :double) (ret sf-result))
   :documentation
-  "@math{(\exp(x)-1)/x} using an algorithm that is accurate for small @math{x}.  For small @math{x} the algorithm is based on the expansion @math{(\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + \dots}."
-  :return (sf-result))
+  "@math{(\exp(x)-1)/x} using an algorithm that is accurate for small @math{x}.
+  For small @math{x} the algorithm is based on the expansion
+  @math{(\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + \dots}.")
 
-(defun-gsl exprel-2 ((x :double))
-  "gsl_sf_exprel_2_e"
+(defun-gsl exprel-2 (x)
+  "gsl_sf_exprel_2_e" ((x :double) (ret sf-result))
   :documentation
-  "@math{2(\exp(x)-1-x)/x^2} using an algorithm that is accurate for small @math{x}.  For small @math{x} the algorithm is based on the expansion @math{2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + \dots}."
-  :return (sf-result))
+  "@math{2(\exp(x)-1-x)/x^2} using an algorithm that is accurate for small
+   @math{x}.  For small @math{x} the algorithm is based on the expansion
+   @math{2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + \dots}.")
 
-(defun-gsl exprel-n ((n :int) (x :double))
-  "gsl_sf_exprel_n_e"
+(defun-gsl exprel-n (n x)
+  "gsl_sf_exprel_n_e" ((n :int) (x :double) (ret sf-result))
   :documentation
-  "@math{N}-relative exponential, which is the @var{n}-th generalization of the functions @code{gsl_sf_exprel} and @code{gsl_sf_exprel2}."
-  :return (sf-result))
+  "@math{N}-relative exponential, which is the @var{n}-th generalization
+   of the functions @code{gsl_sf_exprel} and @code{gsl_sf_exprel2}.")
 
 ;;;;****************************************************************************
 ;;;; Exponentiation With Error Estimate
 ;;;;****************************************************************************
 
-(defun-gsl exp-err ((x :double) (dx :double))
-  "gsl_sf_exp_err_e"
+(defun-gsl exp-err (x dx)
+  "gsl_sf_exp_err_e" ((x :double) (dx :double) (ret sf-result))
   :documentation
-  "Exponentiate @var{x} with an associated absolute error @var{dx}."
-  :return (sf-result))
+  "Exponentiate @var{x} with an associated absolute error @var{dx}.")
 
-(defun-gsl exp-err-scaled ((x :double) (dx :double))
+(defun-gsl exp-err-scaled (x dx)
   "gsl_sf_exp_err_e10_e"
+  ((x :double) (dx :double) (ret sf-result))
   :documentation
-  "Exponentiate @var{x} with an associated absolute error @var{dx} and with extended numeric range."
-  :return (sf-result))
+  "Exponentiate @var{x} with an associated absolute error @var{dx}
+  and with extended numeric range.")
 
-(defun-gsl exp-mult-err ((x :double) (dx :double) (y :double) (dy :double))
+(defun-gsl exp-mult-err (x dx y dy)
   "gsl_sf_exp_mult_err_e"
+  ((x :double) (dx :double) (y :double) (dy :double) (ret sf-result))
   :documentation
-  "The product @math{y \exp(x)} for the quantities @var{x}, @var{y} with associated absolute errors @var{dx}, @var{dy}."
-  :return (sf-result))
+  "The product @math{y \exp(x)} for the quantities @var{x},
+   @var{y} with associated absolute errors @var{dx}, @var{dy}.")
 
-(defun-gsl exp-mult-err-scaled ((x :double) (y :double))
-  "gsl_sf_exp_mult_err_e10_e"
+(defun-gsl exp-mult-err-scaled (x y)
+  "gsl_sf_exp_mult_err_e10_e" ((x :double) (y :double) (ret sf-result-e10))
   :documentation
-  "The product @math{y \exp(x)} for the quantities @var{x}, @var{y} with associated absolute errors @var{dx}, @var{dy} and with extended numeric range."
-  :return (sf-result-e10))
+  "The product @math{y \exp(x)} for the quantities @var{x}, @var{y}
+   with associated absolute errors @var{dx}, @var{dy} and with
+   extended numeric range.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
@@ -102,9 +100,30 @@ of @math{\exp(x)} would overflow the  numeric range of @code{double}."
   (lisp-unit:assert-first-fp-equal
    "0.200855369232d+02"
    (gsl-exp 3.0d0))
+  (lisp-unit:assert-equal
+   '("0.108003407162d+01" "0.241000000000e+03")
+   (subseq (lisp-unit:fp-values (exp-scaled 555.0d0)) 0 2))
+  (lisp-unit:assert-first-fp-equal
+   "0.365352998968d+45"
+   (exp-mult 101.0d0 5.0d0))
+  (lisp-unit:assert-equal
+   '("0.109083441234d+01" "0.243000000000e+03")
+   (subseq (lisp-unit:fp-values (exp-mult-scaled 555.0d0 101.0d0)) 0 2))
   (lisp-unit:assert-first-fp-equal
    "0.100005000167d-03"
    (expm1 0.0001d0))
   (lisp-unit:assert-first-fp-equal
    "0.100005000167d+01"
-   (exprel 0.0001d0)))
+   (exprel 0.0001d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.100033341668d+01"
+   (exprel-2 0.001d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.100025005001d+01"
+   (exprel-n 3 0.001d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.200855369232d+02"
+   (exp-err 3.0d0 0.001d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.461967349233d+03"
+   (exp-mult-err 3.0d0 0.001d0 23.0d0 0.001d0)))
diff --git a/special-functions/exponential-integrals.lisp b/special-functions/exponential-integrals.lisp
index 0367e59bdc67ac7f71df51adbb7e34fed7d614ec..ed430e60e604138048d10446690271d637159040 100644
--- a/special-functions/exponential-integrals.lisp
+++ b/special-functions/exponential-integrals.lisp
@@ -3,7 +3,7 @@
 ; description: Exponential integrals                     
 ; date:        Tue Mar 21 2006 - 17:37                   
 ; author:      Liam M. Healy                             
-; modified:    Wed Apr 26 2006 - 10:13
+; modified:    Tue Jun 13 2006 - 21:40
 ;********************************************************
 ;;; $Id: $
 
@@ -13,79 +13,78 @@
 ;;;; Exponential Integral
 ;;;;****************************************************************************
 
-(defun-gsl expint-E1 ((x :double))
-  "gsl_sf_expint_E1_e"
+(defun-gsl expint-E1 (x)
+  "gsl_sf_expint_E1_e" ((x :double) (ret sf-result))
   :documentation
-  "The exponential integral @math{E_1(x)}, E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t.."
-  :return (sf-result))
+  "The exponential integral
+   @math{E_1(x)}, E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t..")
 
-(defun-gsl expint-E2 ((x :double))
-  "gsl_sf_expint_E2_e"
+(defun-gsl expint-E2 (x)
+  "gsl_sf_expint_E2_e" ((x :double) (ret sf-result))
   :documentation
-  "The second-order exponential integral @math{E_2(x)}, E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2."
-  :return (sf-result))
+  "The second-order exponential integral
+   @math{E_2(x)}, E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.")
 
 ;;;;****************************************************************************
 ;;;; Ei
 ;;;;****************************************************************************
 
-(defun-gsl expint-Ei ((x :double))
-  "gsl_sf_expint_Ei_e"
+(defun-gsl expint-Ei (x)
+    "gsl_sf_expint_Ei_e" ((x :double) (ret sf-result))
   :documentation
-  "The exponential integral @math{Ei(x)}, Ei(x) := - PV\left(\int_{-x}^\infty dt \exp(-t)/t\right)."
-  :return (sf-result))
+  "The exponential integral @math{Ei(x)},
+   Ei(x) := - PV\left(\int_{-x}^\infty dt \exp(-t)/t\right).")
 
 ;;;;****************************************************************************
 ;;;; Hyperbolic Integrals
 ;;;;****************************************************************************
 
-(defun-gsl Shi ((x :double))
-  "gsl_sf_Shi_e"
+(defun-gsl Shi (x)
+  "gsl_sf_Shi_e" ((x :double) (ret sf-result))
   :documentation
-  "The integral @math{Shi(x) = \int_0^x dt \sinh(t)/t}."
-  :return (sf-result))
+  "The integral @math{Shi(x) = \int_0^x dt \sinh(t)/t}.")
 
-(defun-gsl Chi ((x :double))
-  "gsl_sf_Chi_e"
+(defun-gsl Chi (x)
+  "gsl_sf_Chi_e" ((x :double) (ret sf-result))
   :documentation
-  "The integral @math{ Chi(x) := \Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh[t]-1)/t] }, where @math{\gamma_E} is the Euler constant."
-  :return (sf-result))
+  "The integral
+   @math{ Chi(x) := \Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh[t]-1)/t] },
+   where @math{\gamma_E} is the Euler constant.")
 
 ;;;;****************************************************************************
 ;;;; Ei-3
 ;;;;****************************************************************************
 
-(defun-gsl expint-3 ((x :double))
-  "gsl_sf_expint_3_e"
+(defun-gsl expint-3 (x)
+  "gsl_sf_expint_3_e" ((x :double) (ret sf-result))
   :documentation
-  "The third-order exponential integral @math{Ei_3(x) = \int_0^xdt \exp(-t^3)} for @c{$x \ge 0$} @math{x >= 0}."
-  :return (sf-result))
+  "The third-order exponential integral @math{Ei_3(x) = \int_0^xdt \exp(-t^3)}
+  for @math{x >= 0}.")
 
 ;;;;****************************************************************************
 ;;;; Trigonometric Integrals
 ;;;;****************************************************************************
 
-(defun-gsl Si ((x :double))
-  "gsl_sf_Si_e"
+(defun-gsl Si (x)
+  "gsl_sf_Si_e" ((x :double) (ret sf-result))
   :documentation
-  "The Sine integral @math{Si(x) = \int_0^x dt \sin(t)/t}."
-  :return (sf-result))
+  "The Sine integral @math{Si(x) = \int_0^x dt \sin(t)/t}.")
 
-(defun-gsl Ci ((x :double))
-  "gsl_sf_Ci_e"
+(defun-gsl Ci (x)
+  "gsl_sf_Ci_e" ((x :double) (ret sf-result))
   :documentation
-  "The Cosine integral @math{Ci(x) = -\int_x^\infty dt \cos(t)/t} for @math{x > 0}."
-  :return (sf-result))
+  "The Cosine integral @math{Ci(x) = -\int_x^\infty dt \cos(t)/t}
+   for @math{x > 0}.")
 
 ;;;;****************************************************************************
 ;;;; Trigonometric Integrals
 ;;;;****************************************************************************
 
-(defun-gsl atanint ((x :double))
-  "gsl_sf_atanint_e"
+(defun-gsl atanint (x)
+  "gsl_sf_atanint_e" ((x :double) (ret sf-result))
   :documentation
-  "The Arctangent integral, which is defined as @math{AtanInt(x) = \int_0^x dt \arctan(t)/t}."
-  :return (sf-result))
+  "The Arctangent integral, which is defined as
+   @math{AtanInt(x) = \int_0^x dt \arctan(t)/t}.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
diff --git a/special-functions/fermi-dirac.lisp b/special-functions/fermi-dirac.lisp
index 19037dc64685951c6e19b35ef802a1deabe5b2a5..2f793d6f7e809122e2b0dd7aeb84b274e8ce7938 100644
--- a/special-functions/fermi-dirac.lisp
+++ b/special-functions/fermi-dirac.lisp
@@ -3,7 +3,7 @@
 ; description: Fermi-Dirac function.                     
 ; date:        Sat Apr 22 2006 - 16:12                   
 ; author:      Liam M. Healy                             
-; modified:    Wed Apr 26 2006 - 09:51
+; modified:    Tue Jun 13 2006 - 21:45
 ;********************************************************
 ;;; $Id: $
 
@@ -13,70 +13,60 @@
 ;;;; Complete Fermi-Dirac Integrals
 ;;;;****************************************************************************
 
-(defun-gsl fermi-dirac-m1 ((x :double))
-  "gsl_sf_fermi_dirac_m1_e"
+(defun-gsl fermi-dirac-m1 (x)
+  "gsl_sf_fermi_dirac_m1_e" ((x :double) (ret sf-result))
   :documentation
   "The complete Fermi-Dirac integral with an index of @math{-1}. 
-  This integral is given by @c{$F_{-1}(x) = e^x / (1 + e^x)$}
-  @math{F_@{-1@}(x) = e^x / (1 + e^x)}."
-  :return (sf-result))
+  This integral is given by @math{F_@{-1@}(x) = e^x / (1 + e^x)}.")
 
-(defun-gsl fermi-dirac-0 ((x :double))
-  "gsl_sf_fermi_dirac_0_e"
+(defun-gsl fermi-dirac-0 (x)
+  "gsl_sf_fermi_dirac_0_e" ((x :double) (ret sf-result))
   :documentation
   "The complete Fermi-Dirac integral with an index of @math{0}. 
-   This integral is given by @math{F_0(x) = \ln(1 + e^x)}."
-  :return (sf-result))
+   This integral is given by @math{F_0(x) = \ln(1 + e^x)}.")
 
-(defun-gsl fermi-dirac-1 ((x :double))
-  "gsl_sf_fermi_dirac_1_e"
+(defun-gsl fermi-dirac-1 (x)
+  "gsl_sf_fermi_dirac_1_e" ((x :double) (ret sf-result))
   :documentation
   "The complete Fermi-Dirac integral with an index of @math{1},
-   @math{F_1(x) = \int_0^\infty dt (t /(\exp(t-x)+1))}."
-  :return (sf-result))
+   @math{F_1(x) = \int_0^\infty dt (t /(\exp(t-x)+1))}.")
 
-(defun-gsl fermi-dirac-2 ((x :double))
-  "gsl_sf_fermi_dirac_2_e"
+(defun-gsl fermi-dirac-2 (x)
+  "gsl_sf_fermi_dirac_2_e" ((x :double) (ret sf-result))
   :documentation
   "The complete Fermi-Dirac integral with an index of @math{2},
-  @math{F_2(x) = (1/2) \int_0^\infty dt (t^2 /(\exp(t-x)+1))}."
-  :return (sf-result))
+  @math{F_2(x) = (1/2) \int_0^\infty dt (t^2 /(\exp(t-x)+1))}.")
 
-(defun-gsl fermi-dirac-integral ((j :int) (x :double))
-  "gsl_sf_fermi_dirac_int_e"
+(defun-gsl fermi-dirac-integral (j x)
+  "gsl_sf_fermi_dirac_int_e" ((j :int) (x :double) (ret sf-result))
   :documentation
   "The complete Fermi-Dirac integral with an integer index of @math{j},
-  @math{F_j(x) = (1/\Gamma(j+1)) \int_0^\infty dt (t^j /(\exp(t-x)+1))}."
-  :return (sf-result))
+  @math{F_j(x) = (1/\Gamma(j+1)) \int_0^\infty dt (t^j /(\exp(t-x)+1))}.")
 
-(defun-gsl fermi-dirac-m1/2 ((x :double))
-  "gsl_sf_fermi_dirac_mhalf_e"
+(defun-gsl fermi-dirac-m1/2 (x)
+  "gsl_sf_fermi_dirac_mhalf_e" ((x :double) (ret sf-result))
   :documentation
-  "The complete Fermi-Dirac integral @c{$F_{-1/2}(x)$}"
-  :return (sf-result))
+  "The complete Fermi-Dirac integral @c{$F_{-1/2}(x)$}")
 
-(defun-gsl fermi-dirac-1/2 ((x :double))
-  "gsl_sf_fermi_dirac_half_e"
+(defun-gsl fermi-dirac-1/2 (x)
+  "gsl_sf_fermi_dirac_half_e" ((x :double) (ret sf-result))
   :documentation
-  "The complete Fermi-Dirac integral @c{$F_{1/2}(x)$}."
-  :return (sf-result))
+  "The complete Fermi-Dirac integral @c{$F_{1/2}(x)$}.")
 
-(defun-gsl fermi-dirac-3/2 ((x :double))
-  "gsl_sf_fermi_dirac_3half_e"
+(defun-gsl fermi-dirac-3/2 (x)
+  "gsl_sf_fermi_dirac_3half_e" ((x :double) (ret sf-result))
   :documentation
-  "The complete Fermi-Dirac integral @c{$F_{3/2}(x)$}."
-  :return (sf-result))
+  "The complete Fermi-Dirac integral @c{$F_{3/2}(x)$}.")
 
 ;;;;****************************************************************************
 ;;;; Incomplete Fermi-Dirac Integrals
 ;;;;****************************************************************************
 
-(defun-gsl fermi-dirac-inc-0 ((x :double) (b :double))
-  "gsl_sf_fermi_dirac_inc_0_e"
+(defun-gsl fermi-dirac-inc-0 (x b)
+  "gsl_sf_fermi_dirac_inc_0_e" ((x :double) (b :double) (ret sf-result))
   :documentation
   "The incomplete Fermi-Dirac integral with an index
-  of zero, @c{$F_0(x,b) = \ln(1 + e^{b-x}) - (b-x)$}."
-  :return (sf-result))
+  of zero, @c{$F_0(x,b) = \ln(1 + e^{b-x}) - (b-x)$}.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
@@ -98,9 +88,15 @@
   (lisp-unit:assert-first-fp-equal
    "0.666882708765d+04"
    (fermi-dirac-integral 5 12.35d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.146429458909d+01"
+   (fermi-dirac-m1/2 2.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.282372127740d+01"
    (fermi-dirac-1/2 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.416541445987d+01"
+   (fermi-dirac-3/2 2.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.170141327798d+01"
    (fermi-dirac-inc-0 2.0d0 0.5d0)))
diff --git a/special-functions/gamma.lisp b/special-functions/gamma.lisp
index 697280c114a8cef22f8af74d642d95674b6c70c3..b96533eeccf5025ec9035c0cf02205f7a9831e1b 100644
--- a/special-functions/gamma.lisp
+++ b/special-functions/gamma.lisp
@@ -3,7 +3,7 @@
 ; description: Gamma functions                           
 ; date:        Thu Apr 27 2006 - 22:06                   
 ; author:      Liam M. Healy                             
-; modified:    Fri Apr 28 2006 - 00:09
+; modified:    Tue Jun 13 2006 - 22:16
 ;********************************************************
 ;;; $Id: $
 
@@ -19,51 +19,47 @@
 
 (defconstant +gamma-xmax+ 171.0d0)
 
-(defun-gsl gamma ((x :double))
-  "gsl_sf_gamma_e"
-  :return (sf-result)
+(defun-gsl gamma (x)
+  "gsl_sf_gamma_e" ((x :double) (ret sf-result))
   :documentation "The Gamma function @math{\Gamma(x)}, subject to x
    not being a negative integer.  The function is computed using the real
    Lanczos method. The maximum value of @math{x} such that
    @math{\Gamma(x)} is not considered an overflow is given by +gamma-xmax+.")
 
-(defun-gsl log-gamma ((x :double))
-  "gsl_sf_lngamma_e"
-  :return (sf-result)
+(defun-gsl log-gamma (x)
+  "gsl_sf_lngamma_e" ((x :double) (ret sf-result))
   :documentation "The logarithm of the Gamma function,
    @math{\log(\Gamma(x))}, subject to @math{x} not a being negative
    integer.  For @math{x<0} the real part of @math{\log(\Gamma(x))} is
    returned, which is equivalent to @math{\log(|\Gamma(x)|)}.  The function
    is computed using the real Lanczos method.")
 
-(defun-gsl log-gamma-sign ((x :double))
-  "gsl_sf_lngamma_sgn_e"
-  :return (sf-result :double)
+(defun-gsl log-gamma-sign (x)
+  "gsl_sf_lngamma_sgn_e" ((x :double) (ret sf-result) (sign :double))
   :documentation "Compute the sign of the gamma function and the logarithm of
   its magnitude, subject to @math{x} not being a negative integer.  The
   function is computed using the real Lanczos method.  The value of the
   gamma function can be reconstructed using the relation @math{\Gamma(x) =
-  sgn * \exp(resultlg)}.")
+  sgn * \exp(resultlg)}."
+  :return ((val ret) (double-to-cl sign) (err ret)))
 
-(defun-gsl gamma* ((x :double))
-  "gsl_sf_gammastar_e"
-  :return (sf-result)
+(defun-gsl gamma* (x)
+  "gsl_sf_gammastar_e" ((x :double) (ret sf-result))
   :Documentation "The regulated Gamma Function @math{\Gamma^*(x)}
   for @math{x > 0}, given by
   \Gamma^*(x) &= \Gamma(x)/(\sqrt{2\pi} x^{(x-1/2)} \exp(-x))\cr
             &= \left(1 + {1 \over 12x} + ...\right)
   \quad\hbox{for~} x\to \infty\cr.")
 
-(defun-gsl 1/gamma ((x :double))
-  "gsl_sf_gammainv_e"
-  :return (sf-result)
+(defun-gsl 1/gamma (x)
+  "gsl_sf_gammainv_e" ((x :double) (ret sf-result))
   :documentation "The reciprocal of the gamma function,
   @math{1/\Gamma(x)} using the real Lanczos method.")
 
-#|
-(defun-gsl log-gamma-complex ((z gsl-complex))
+(defun-gsl log-gamma-complex (z)
   "gsl_sf_lngamma_complex_e"
-  :return (sf-result)
+  (((realpart z) :double) ((imagpart z) :double)
+   (lnr sf-result) (arg sf-result))
   :documentation "Compute @math{\log(\Gamma(z))} for complex @math{z=z_r+i
   z_i} and @math{z} not a negative integer, using the complex Lanczos
   method.  The returned parameters are @math{lnr = \log|\Gamma(z)|} and
@@ -71,121 +67,108 @@
   part (@var{arg}) is not well-determined when @math{|z|} is very large,
   due to inevitable roundoff in restricting to @math{(-\pi,\pi]}.  This
   will result in a @code{GSL_ELOSS} error when it occurs.  The absolute
-  value part (@var{lnr}), however, never suffers from loss of precision.")
-|#
+  value part (@var{lnr}), however, never suffers from loss of precision."
+  :return
+  ((val lnr) (val arg) (err lnr) (err arg)))
 
-(defun-gsl taylor-coefficient ((n :int) (x :double))
-  "gsl_sf_taylorcoeff_e"
-  :return (sf-result)
+(defun-gsl taylor-coefficient (n x)
+  "gsl_sf_taylorcoeff_e" ((n :int) (x :double) (ret sf-result))
   :documentatiOn "Compute the Taylor coefficient @math{x^n / n!} for 
   @math{x >= 0}, @math{n >= 0}.")
 
-(defun-gsl factorial ((n :size))
-  "gsl_sf_fact_e"
-  :return (sf-result)
+(defun-gsl factorial (n)
+  "gsl_sf_fact_e" ((n :size) (ret sf-result))
   :documentation "The factorial @math{n!},
   related to the Gamma function by @math{n! = \Gamma(n+1)}.")
 
-(defun-gsl double-factorial ((n :size))
-  "gsl_sf_doublefact_e"
-  :return (sf-result)
+(defun-gsl double-factorial (n)
+  "gsl_sf_doublefact_e" ((n :size) (ret sf-result))
   :documentation "The double factorial @math{n!! = n(n-2)(n-4) \dots}.")
 
-(defun-gsl log-factorial ((n :size))
-  "gsl_sf_lnfact_e"
-  :return (sf-result)
+(defun-gsl log-factorial (n)
+  "gsl_sf_lnfact_e" ((n :size) (ret sf-result))
   :documentation "The logarithm of the factorial of @var{n},
   @math{\log(n!)}.  The algorithm is faster than computing
   @math{\ln(\Gamma(n+1))} via @code{gsl_sf_lngamma} for @math{n < 170},
   but defers for larger @var{n}.")
 
-(defun-gsl log-double-factorial ((n :size))
-  "gsl_sf_lndoublefact_e"
-  :return (sf-result)
+(defun-gsl log-double-factorial (n)
+  "gsl_sf_lndoublefact_e" ((n :size) (ret sf-result))
   :documentation "These routines compute the logarithm of
   the double factorial of @var{n}, @math{\log(n!!)}.")
 
-(defun-gsl choose ((n :size) (m :size))
-  "gsl_sf_choose_e"
-  :return (sf-result)
+(defun-gsl choose (n m)
+  "gsl_sf_choose_e" ((n :size) (m :size) (ret sf-result))
   :documentation "The combinatorial factor @code{n choose m}
   @math{= n!/(m!(n-m)!)}")
 
-(defun-gsl log-choose ((n :size) (m :size))
-  "gsl_sf_lnchoose_e"
-  :return (sf-result)
+(defun-gsl log-choose (n m)
+  "gsl_sf_lnchoose_e" ((n :size) (m :size) (ret sf-result))
   :documentation "The logarithm of @code{n choose m}.  This is
   equivalent to the sum @math{\log(n!) - \log(m!) - \log((n-m)!)}.")
 
-(defun-gsl pochammer ((a :double) (x :double))
-  "gsl_sf_poch_e"
-  :return (sf-result)
+(defun-gsl pochammer (a x)
+  "gsl_sf_poch_e"  ((a :double) (x :double) (ret sf-result))
   :documentation "The Pochhammer symbol @math{(a)_x := \Gamma(a +
    x)/\Gamma(a)}, subject to @math{a} and @math{a+x} not being negative
    integers. The Pochhammer symbol is also known as the Apell symbol and
    sometimes written as @math{(a,x)}.")
 
-(defun-gsl log-pochammer ((a :double) (x :double))
-  "gsl_sf_lnpoch_e"
-  :return (sf-result)
+(defun-gsl log-pochammer (a x)
+  "gsl_sf_lnpoch_e" ((a :double) (x :double) (ret sf-result))
   :documentation "The logarithm of the Pochhammer symbol,
   @math{\log((a)_x) = \log(\Gamma(a + x)/\Gamma(a))} for @math{a > 0},
   @math{a+x > 0}.")
 
-(defun-gsl log-pochammer-sign ((a :double) (x :double))
-  "gsl_sf_lnpoch_e"
-  :return (sf-result :double)
+(defun-gsl log-pochammer-sign (a x)
+  "gsl_sf_lnpoch_sgn_e"
+  ((a :double) (x :double) (ret sf-result) (sign :double))
   :documentation "The logarithm of the Pochhammer symbol and its sign.
   The computed parameters are @math{result =
   \log(|(a)_x|)} and @math{sgn = \sgn((a)_x)} where @math{(a)_x :=
   \Gamma(a + x)/\Gamma(a)}, subject to @math{a}, @math{a+x} not being
-  negative integers.")
+  negative integers."
+  :return ((val ret) (double-to-cl sign) (err ret)))
 
-(defun-gsl relative-pochammer ((a :double) (x :double))
-  "gsl_sf_pochrel_e"
-  :return (sf-result)
+(defun-gsl relative-pochammer (a x)
+  "gsl_sf_pochrel_e" ((a :double) (x :double) (ret sf-result))
   :documentation "The relative Pochhammer symbol @math{((a)_x -
   1)/x} where @math{(a)_x := \Gamma(a + x)/\Gamma(a)}.")
 
-(defun-gsl incomplete-gamma ((a :double) (x :double))
-  "gsl_sf_gamma_inc_Q_e"
-  :return (sf-result)
+(defun-gsl incomplete-gamma (a x)
+  "gsl_sf_gamma_inc_Q_e" ((a :double) (x :double) (ret sf-result))
   :documentation "The normalized incomplete Gamma Function
   @math{Q(a,x) = 1/\Gamma(a) \int_x^\infty dt t^@{a-1@} \exp(-t)}
   for @math{a > 0}, @math{x >= 0}.")
 
-(defun-gsl complementary-incomplete-gamma ((a :double) (x :double))
-  "gsl_sf_gamma_inc_P_e"
-  :return (sf-result)
+(defun-gsl complementary-incomplete-gamma (a x)
+  "gsl_sf_gamma_inc_P_e" ((a :double) (x :double) (ret sf-result))
   :documentation "The complementary normalized incomplete Gamma Function
   @math{P(a,x) = 1/\Gamma(a) \int_0^x dt t^@{a-1@} \exp(-t)}
   for @math{a > 0}, @math{x >= 0}.  Note that Abramowitz & Stegun
   call @math{P(a,x)} the incomplete gamma function (section 6.5).")
 
-(defun-gsl nonnormalized-incomplete-gamma ((a :double) (x :double))
-  "gsl_sf_gamma_inc_e"
-  :return (sf-result)
+(defun-gsl nonnormalized-incomplete-gamma (a x)
+  "gsl_sf_gamma_inc_e" ((a :double) (x :double) (ret sf-result))
   :documentation "The incomplete Gamma Function
    @math{\Gamma(a,x)}, without the normalization factor
    included in the previously defined functions:
    @math{\Gamma(a,x) = \int_x^\infty dt t^@{a-1@} \exp(-t)}
    for @math{a} real and @math{x >= 0}.")
 
-(defun-gsl beta ((a :double) (b :double))
-  "gsl_sf_beta_e"
-  :return (sf-result)
+(defun-gsl beta (a b)
+  "gsl_sf_beta_e" ((a :double) (b :double) (ret sf-result))
   :documentation "The Beta Function, @math{B(a,b) =
   \Gamma(a)\Gamma(b)/\Gamma(a+b)} for @math{a > 0}, @math{b > 0}.")
 
-(defun-gsl log-beta ((a :double) (b :double))
-  "gsl_sf_lnbeta_e"
-  :return (sf-result)
+(defun-gsl log-beta (a b)
+  "gsl_sf_lnbeta_e" ((a :double) (b :double) (ret sf-result))
   :documentation "The logarithm of the Beta Function,
   @math{\log(B(a,b))} for @math{a > 0}, @math{b > 0}.")
 
-(defun-gsl incomplete-beta ((a :double) (b :double) (x :double))
+(defun-gsl incomplete-beta (a b x)
   "gsl_sf_beta_inc_e"
-  :return (sf-result)
+  ((a :double) (b :double) (x :double) (ret sf-result))
   :documentation "The normalized incomplete Beta function
    @math{B_x(a,b)/B(a,b)} where
    @math{B_x(a,b) = \int_0^x t^@{a-1@} (1-t)^@{b-1@} dt}
@@ -204,12 +187,20 @@
   (lisp-unit:assert-first-fp-equal
    "0.359134205370d+03"
    (log-gamma 100.0d0))
+  (lisp-unit:assert-equal
+   '("0.359134205370d+03" "0.100000000000d+01")
+   (subseq (lisp-unit:fp-values (log-gamma-sign 100.0d0)) 0 2))
   (lisp-unit:assert-first-fp-equal
    "0.100347805583d+01"
    (gamma* 24.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.198412698413d-03"
    (1/gamma 8.0d0))
+  (lisp-unit:assert-equal
+   '("0.823613175045d+01" "-0.118403781494d+01")
+   (subseq
+    (lisp-unit:fp-values (log-gamma-complex #C(10.0d0 10.0d0)))
+    0 2))
   (lisp-unit:assert-first-fp-equal
    "0.110947646104d-02"
    (taylor-coefficient 12 3.0d0))
@@ -222,6 +213,9 @@
   (lisp-unit:assert-first-fp-equal
    "0.857933669826d+03"
    (log-factorial 199))
+  (lisp-unit:assert-first-fp-equal
+   "0.430177893581d+03"
+   (log-double-factorial 199))
   (lisp-unit:assert-first-fp-equal
    "0.560000000000d+02"
    (choose 8 3))
@@ -232,9 +226,32 @@
   (lisp-unit:assert-first-fp-equal
    "0.120000000000d+02"
    (pochammer 3.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.863231987192d+03"
+   (log-pochammer 2.0d0 199.0d0))
+  (lisp-unit:assert-equal
+   '("0.863231987192d+03" "0.100000000000d+01")
+   (subseq (lisp-unit:fp-values
+	    (log-pochammer-sign 2.0d0 199.0d0))
+	   0 2))
+  (lisp-unit:assert-first-fp-equal
+   "0.403199888889d+06"
+   (relative-pochammer 2.0d0 9.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.406005849710d+00"
+   (incomplete-gamma 2.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.593994150290d+00"
+   (complementary-incomplete-gamma 2.0d0 2.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.406005849710d+00"
+   (nonnormalized-incomplete-gamma 2.0d0 2.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.181818181818d+00"
    (beta 5.50d0 1.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "-0.170474809224d+01"
+   (log-beta 5.5d0 1.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.646446609407d+00"
    (incomplete-beta 1.0d0 1.50d0 0.50d0)))
diff --git a/special-functions/gegenbauer.lisp b/special-functions/gegenbauer.lisp
index 7344e60031d9a98cc61e88c9bcde8a28fe58218d..8f97259763c0266278bc48d5e8dd5608b64fb5c8 100644
--- a/special-functions/gegenbauer.lisp
+++ b/special-functions/gegenbauer.lisp
@@ -3,45 +3,41 @@
 ; description: Gegenbauer polynomials                    
 ; date:        Fri Apr 28 2006 - 20:40                   
 ; author:      Liam M. Healy                             
-; modified:    Fri Apr 28 2006 - 22:44
+; modified:    Tue Jun 13 2006 - 22:56
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl gegenbauer-1 ((lambda :double) (x :double))
-  "gsl_sf_gegenpoly_1_e"
-  :return (sf-result)
+(defun-gsl gegenbauer-1 (lambda x)
+  "gsl_sf_gegenpoly_1_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_1(x)}.")
 
-(defun-gsl gegenbauer-2 ((lambda :double) (x :double))
-  "gsl_sf_gegenpoly_2_e"
-  :return (sf-result)
+(defun-gsl gegenbauer-2 (lambda x)
+    "gsl_sf_gegenpoly_2_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_2(x)}.")
 
-(defun-gsl gegenbauer-3 ((lambda :double) (x :double))
-  "gsl_sf_gegenpoly_3_e"
-  :return (sf-result)
+(defun-gsl gegenbauer-3 (lambda x)
+  "gsl_sf_gegenpoly_3_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The Gegenbauer polynomial @math{C^@{(\lambda)@}_3(x)}.")
 
-(defun-gsl gegenbauer ((n :int) (lambda :double) (x :double))
+(defun-gsl gegenbauer (n lambda x)
   "gsl_sf_gegenpoly_n_e"
-  :return (sf-result)
+  ((n :int) (lambda :double) (x :double) (ret sf-result))
   :documentation "The Gegenbauer polynomial
   @math{C^@{(\lambda)@}_n(x)} for a specific value of @var{n},
   @var{lambda}, @var{x} subject to @math{\lambda > -1/2},
   @math{n >= 0}.")
 
-(defun-gsl gegenbauer-array
-    (((dim0 result) :int)
-     (lambda :double) (x :double) ((gsl-array result) :pointer))
+(defun-gsl gegenbauer-array (lambda x result)
   "gsl_sf_gegenpoly_array"
-  :function (lambda x result)
+  (((1- (dim0 result)) :int)
+   (lambda :double) (x :double) ((gsl-array result) :pointer))
   :documentation "Compute an array of Gegenbauer polynomials
-  @math{C^@{(\lambda)@}_n(x)} for @math{n = 0, 1, 2, \dots, nmax}, subject
-  to @math{\lambda > -1/2}, @math{nmax >= 0}."
-  :after ((cl-invalidate result))
-  :return-input (result))
+  @math{C^@{(\lambda)@}_n(X)} for
+  @math{n = 0, 1, 2, \dots, length(result)-1}, subject
+  to @math{\lambda > -1/2}."
+  :invalidate (result))
 
 ;;; (defparameter vec (make-data 'vector nil 3))
 ;;; (gegenbauer-array 1.0d0 3.0d0 vec)
@@ -59,4 +55,13 @@
    (gegenbauer-2 1.0d0 3.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.204000000000d+03"
-   (gegenbauer-3 1.0d0 3.0d0)))
+   (gegenbauer-3 1.0d0 3.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.118900000000d+04"
+   (gegenbauer 4 1.0d0 3.0d0))
+  (lisp-unit:assert-equal
+   '("0.100000000000d+01" "0.600000000000d+01" "0.350000000000d+02"
+     "0.204000000000d+03")
+   (lisp-unit:fp-sequence
+    (with-data (arr vector-double 4)
+      (gegenbauer-array 1.0d0 3.0d0 arr) (data arr)))))