diff --git a/README b/README
index a5d76e549ffbc4747947e28c5a1ab4d876b9e131..c5ff361d35f0a50a7b8604fa74ed9a04a4782a4a 100644
--- a/README
+++ b/README
@@ -1,5 +1,4 @@
 TODO
-remove cffi-array
 Document letm
 Document chapters
 organize presentation
@@ -9,6 +8,8 @@ Other regression errors
 Do tests in ia32.
 normalize terminology for solve-minimize-fit; make subdir?
 Add bspline.
+Include required parts per FDL.
+
 
 Wishlist: 
 Sun Feb 10 2008 GSLL creation of objects within letm as from e.g. #'covariance.
@@ -16,6 +17,8 @@ Sun Feb 10 2008 GSLL creation of objects within letm as from e.g. #'covariance.
 Bulk changes needed:
 Rename defun-gsl, defvariable.
 Clean up markup and header files.
+:size to size
+String comparison of floats?
 
 
 Features:
diff --git a/init/init.lisp b/init/init.lisp
index 61e4df103401b53fdeefc01c06c59c3071e54633..3a14bb2eff2b4f3bf8f477974dddfb7425395b38 100644
--- a/init/init.lisp
+++ b/init/init.lisp
@@ -1,6 +1,6 @@
 ;; Load GSL
 ;; Liam Healy Sat Mar  4 2006 - 18:53
-;; Time-stamp: <2008-01-21 12:33:11EST init.lisp>
+;; Time-stamp: <2008-02-16 17:56:56EST init.lisp>
 ;; $Id: $
 
 (defpackage gsll
@@ -25,6 +25,17 @@
 ;;; be changed to something outside the keyword package to suppress
 ;;; warnings from CFFI.
 
+(cffi:defctype uint32 :unsigned-int)
+(cffi:defctype uint64 :unsigned-long)
+(cffi:defctype size
+   #-cffi-features:no-long-long uint64
+   #+cffi-features:no-long-long
+   #.(progn (cerror "Use uint32 instead."
+   "This platform does not support long long types.")
+   uint32))
+
+;;; obsolete
+
 (cffi:defctype :uint32 :unsigned-int)
 (cffi:defctype :uint64 :unsigned-long)
 (cffi:defctype :size
diff --git a/init/interface.lisp b/init/interface.lisp
index 78a88fd56dc21707cc9641cb5c3d6e73d31b8357..8cfaa6f1daa322904f054c530ee9aa2e862cf350 100644
--- a/init/interface.lisp
+++ b/init/interface.lisp
@@ -1,6 +1,6 @@
 ;; Macros to interface GSL functions.
 ;; Liam Healy 
-;; Time-stamp: <2008-02-16 10:44:10EST interface.lisp>
+;; Time-stamp: <2008-02-16 17:47:30EST interface.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -152,6 +152,107 @@
 ;;; invalidate   Invalidate the CL array/matrix cache.
 ;;; after        After method.
 ;;; enumeration  The name of the enumeration return.
+(defmacro defmfun
+    (name arglist gsl-name c-arguments
+     &key (c-return :error-code)
+     (return nil return-supplied-p)
+     (type :function) (index t) (export (not (eq type :method)))
+     null-pointer-info documentation invalidate after enumeration)
+  (let* ((cargs (substitute '(mode sf-mode) :mode c-arguments))
+	 (carg-symbs
+	  (remove-if-not #'symbolp
+			 (mapcar #'st-symbol (remove :mode c-arguments))))
+	 (clargs
+	  (or arglist carg-symbs))
+	 (arglist-symbs
+	  (when arglist			; can be method, so get symbol
+	    (mapcar (lambda (x) (if (listp x) (first x) x)) arglist)))
+	 (cret-type (if (member c-return *special-c-return*)
+			:int
+			(if (listp c-return) (st-type c-return) c-return)))
+	 (cret-name
+	  (if (listp c-return) (st-symbol c-return) (make-symbol "CRETURN")))
+	 (allocated		     ; Foreign objects to be allocated
+	  (remove-if (lambda (s) (member s arglist-symbs)) carg-symbs))
+	 (allocated-decl
+	  (mapcar
+	   (lambda (s) (find s cargs :key #'st-symbol))
+	   allocated))
+	 (clret (or (substitute cret-name :c-return return) ; better as a symbol macro
+		    (mapcan #'cl-convert-form allocated-decl)
+		    invalidate
+		    (unless (eq c-return :void)
+		      (list cret-name))))
+	 (clargs-types (cl-argument-types clargs cargs)))
+    `(progn
+      (,(if (eq type :function) 'defun 'defmethod)
+       ,name
+       ,(if (member :mode c-arguments)
+	    `(,@clargs &optional (mode :double-prec))
+	    `(,@clargs))
+       ,(declaration-form clargs-types)
+       ,@(when documentation (list documentation))
+       (,@(if allocated
+	      `(cffi:with-foreign-objects
+		,(mapcar #'wfo-declare allocated-decl))
+	      '(let ()))
+	(let ((,cret-name
+	       (cffi:foreign-funcall
+		,gsl-name
+		,@(mapcan
+		   (lambda (arg)
+		     (list (if (member (st-symbol arg) allocated)
+			       :pointer
+			       (st-type arg))
+			   (st-symbol arg)))
+		   cargs)
+		,cret-type)))
+	  ,@(case c-return
+		  (:void `((declare (ignore ,cret-name))))
+		  (:error-code		; fill in arguments
+		   `((check-gsl-status ,cret-name ',name))))
+	  ,@(when invalidate `((cl-invalidate ,@invalidate)))
+	  ,@(when (or null-pointer-info (eq c-return :pointer))
+		  `((check-null-pointer ,cret-name
+		     ,@(or null-pointer-info
+			   '(:ENOMEM "No memory allocated")))))
+	  ,@after
+	  (values
+	   ,@(case c-return
+		   (:number-of-answers
+		    (mapcan
+		     (lambda (vbl seq)
+		       `((when (> ,cret-name ,seq) ,vbl)))
+		     clret
+		     (loop for i below (length clret) collect i)))
+		   (:success-failure
+		    (if (equal clret invalidate)
+			;; success-failure more important than passed-in
+			`((success-failure ,cret-name))
+			(remove cret-name ; don't return c-return itself
+				`(,@clret (success-failure ,cret-name)))))
+		   (:success-continue
+		    (if (equal clret invalidate)
+			;; success-failure more important than passed-in
+			`((success-continue ,cret-name))
+			(remove cret-name ; don't return c-return itself
+				`(,@clret (success-continue ,cret-name)))))
+		   (:true-false
+		    `((not (zerop ,cret-name))))
+		   (:enumerate
+		    `((cffi:foreign-enum-keyword ',enumeration ,cret-name)))
+		   (t (unless
+			  (or
+			   (and (eq c-return :error-code)
+				(not allocated)
+				(not return-supplied-p))
+			   (and (null return) return-supplied-p))
+			clret)))))))
+      ,@(when index `((map-name ',(if (eql index t) name index) ,gsl-name)))
+      ,@(when export `((export ',name))))))
+
+
+;;; to be removed
 (defmacro defun-gsl
     (name arglist gsl-name c-arguments
      &key (c-return :error-code)
@@ -251,6 +352,8 @@
       ,@(when index `((map-name ',(if (eql index t) name index) ,gsl-name)))
       ,@(when export `((export ',name))))))
 
+
+
 (defmacro defun-optionals
     (name arglist no-optional optionals &optional documentation)
   "Define a function with and without optional arguments."
@@ -272,6 +375,15 @@
 ;;;; Variables in library
 ;;;;****************************************************************************
 
+(defmacro defmpar (cl-symbol gsl-symbol documentation)
+  "Define a library variable pointer."
+  `(progn
+    (cffi:defcvar (,gsl-symbol ,cl-symbol :read-only t) :pointer
+      ,documentation)
+    (map-name ',cl-symbol ,gsl-symbol)
+    (export ',cl-symbol)))
+
+;;; to be removed
 (defmacro defvariable (cl-symbol gsl-symbol documentation)
   "Define a library variable pointer."
   `(progn
diff --git a/init/tests.lisp b/init/tests.lisp
index f373cf916217c61ccb73c93d4e4f45362160522c..e5cde39135af2fbf31b06138e5eddd0979b1b043 100644
--- a/init/tests.lisp
+++ b/init/tests.lisp
@@ -11,15 +11,83 @@
 ;;; point numbers and a form for generating floating point tests.
 
 (in-package :lisp-unit)
+
+(defparameter *zero-threshold* 1.0d-12
+  "Threshold below which a number is to be considered zero.")
+
+(defparameter *acceptable-fraction-error* 1.0d-12
+  "Fractional error which is considered acceptable when
+  comparing floating point numbers.")
+
+(defun fp-equal (fp1 fp2)
+  "The floats fp1 and fp2 are to be considered equal."
+  ;; For now, fp1 and fp2 must be actual floats and not nans/infs.
+  (or (and (<= (abs fp1) *zero-threshold*)
+	   (<= (abs fp2) *zero-threshold*))
+      (<= (abs (/ (- fp1 fp2) fp1)) *acceptable-fraction-error*)))
+
+(defun numerical-equal (result1 result2)
+  (and
+   (or
+     (typep result1 'vector) (typep result2 'vector)
+     (typep result1 'complex) (typep result2 'complex)
+     (eql (type-of result1) (type-of result2)))
+   (typecase result1
+     (integer (= result1 result2))
+     (float (fp-equal result1 result2))
+     (complex (and (fp-equal (realpart result1) (realpart result2))
+		   (fp-equal (imagpart result1) (imagpart result2))))
+     (sequence
+      (and (eql (length result1) (length result2))
+	   (every #'numerical-equal result1 result2)))
+     (array
+      (and (= (array-rank result1) (array-rank result2) 2)
+	   (loop for i below (array-dimension result1 0)
+		 always
+		 (loop for j below (array-dimension result1 1)
+		       always (numerical-equal (aref result1 i j)
+					       (aref result2 i j)))))))))
+
+(defun numerical-serialize (form)
+  (if (typep form 'list)
+      (cons 'list (mapcar #'numerical-serialize form))
+      form))
+
+(defmacro assert-numerical-equal (expected form &rest extras)
+  (lisp-unit::expand-assert
+    :equal form form expected extras
+    :test #'numerical-equal))
+
+
+;;; (make-test '(legendre-conicalP-half 3.5d0 10.0d0))
+(defun gsl::make-test (form)
+  "Make a test for lisp-unit."
+  (let ((vals (multiple-value-list (ignore-errors (eval form)))))
+    (if (typep (second vals) 'error)
+	`(lisp-unit::assert-error
+	  ',(type-of (second vals))
+	  ,form)
+	`(lisp-unit::assert-numerical-equal
+	  ,(numerical-serialize vals)
+	  (multiple-value-list ,form)))))
+
+(defmacro gsl::make-tests (name &rest forms)
+  (append
+   `(lisp-unit:define-test ,name)
+   (mapcar #'gsl::make-test forms)))
+
+
+;;;;;;;;;;;;;; OBSOLETE
+
 (export '(assert-first-fp-equal fp-values fp-sequence))
 
+(defmacro fp-values (results)
+  `(mapcar #'fp-record (multiple-value-list ,results)))
+
 (defparameter *test-fp-decimal-digits* 12
   "The number of decimal digits on which floating point
    number must agree.")
 
-(defparameter *zero-threshold* 1.0d-15
-  "Threshold below which a number is tpo be considered zero.")
-
 (defun fp-string (fp &optional (decimal-digits *test-fp-decimal-digits*))
   "Format the floating point number as a string for comparison."
   (if (typep fp 'complex)
@@ -42,14 +110,50 @@
   (map 'list #'fp-string results))
 
 
+
+
+
+#|
+
+
+
 ;;(gsl:double-float-unequal x y double-float-epsilon)
 
-;;; (make-fp-test '(legendre-conicalP-half 3.5d0 10.0d0))
-(defun gsl::make-fp-test (form)
-  "Make a test."
-  `(lisp-unit::assert-first-fp-equal
-    ,(lisp-unit::fp-string (eval form))
-    ,form))
+
+
+;;; Tests:
+;;; assert-float-equal: First value only
+;;; turn multiple values into a list
+;;; compare sequences including structure.
+
+
+
+(defmacro assert-fp-equal (expected form &rest extras)
+  (lisp-unit::expand-assert
+   :equal form `(fp-string (nth-value 0 ,form)) expected extras
+   :test #'fp-equal))
+
+(defun fp-record (fp)
+  "Record the floating point number for comparison."
+  (if (typep fp 'complex)
+      (list (fp-record (realpart fp)) (fp-record (imagpart fp)))
+      ;; relies on ~g printing full precision per spec
+      (format nil "~g" fp)))
+
+(defun fp-record (fp)
+  "Record the floating point number for comparison."
+  (etypecase fp
+    (complex
+     (list (fp-record (realpart fp)) (fp-record (imagpart fp))))
+    (sequence (map (type-of fp) #'fp-record))
+    ;; relies on ~g printing full precision per spec
+    (float (format nil "~g" fp))))
+
+
+(defun fp-sequence (results)
+  (map 'list #'fp-string results))
+
+
 
 ;;;;;;;;;;;;;;;;;
 ;;; Thinking about how to do floating point comparisons
@@ -61,7 +165,6 @@
 	       double-float-epsilon
 	       single-float-epsilon))))
 
-;;; Tests:
-;;; assert-float-equal: First value only
-;;; turn multiple values into a list
-;;; compare sequences including structure.
+
+    
+|#
diff --git a/special-functions/airy.lisp b/special-functions/airy.lisp
index fcf745b5f4288ba022a5adfb9c49288803fcb6ae..1c6b8cfef6ac549a48b289e411f5d884b4e5fe13 100644
--- a/special-functions/airy.lisp
+++ b/special-functions/airy.lisp
@@ -1,10 +1,7 @@
-;********************************************************
-; file:        airy.lisp                                 
-; description: Airy functions                            
-; date:        Fri Mar 17 2006 - 18:41                   
-; author:      Liam M. Healy
-; modified:    Sun Feb 18 2007 - 12:16
-;********************************************************
+;; Airy functions
+;; Liam Healy, Fri Mar 17 2006 - 18:41
+;; Time-stamp: <2008-02-16 17:57:19EST airy.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -12,86 +9,130 @@
 ;;;; Airy functions
 ;;;;****************************************************************************
 
-(defun-gsl airy-Ai (x)
+(defmfun airy-Ai (x)
     "gsl_sf_airy_Ai_e" ((x :double) :mode (ret sf-result))
-  :documentation "The Airy function Ai(x).")
+  :documentation			; FDL
+  "The Airy function Ai(x).")
 
-(defun-gsl airy-Bi (x)
+(defmfun airy-Bi (x)
   "gsl_sf_airy_Bi_e" ((x :double) :mode (ret sf-result))
-  :documentation "The Airy function Bi(x).")
+  :documentation			; FDL
+  "The Airy function Bi(x).")
 
-(defun-gsl airy-Ai-scaled (x)
+(defmfun airy-Ai-scaled (x)
   "gsl_sf_airy_Ai_scaled_e" ((x :double) :mode (ret sf-result))
-  :documentation
-  "The scaled Airy function @math{S_A(x) Ai(x)}.
-   For @math{x>0} the scaling factor @math{S_A(x)} is @math{\exp(+(2/3) x^(3/2))},
-   and is 1 for @math{x<0}.")
+  :documentation			; FDL
+  "The scaled Airy function S_A(x) Ai(x).
+   For x>0 the scaling factor S_A(x) is \exp(+(2/3) x^(3/2)),
+   and is 1 for x<0.")
 
-(defun-gsl airy-Bi-scaled (x)
+(defmfun airy-Bi-scaled (x)
   "gsl_sf_airy_Bi_scaled_e" ((x :double) :mode (ret sf-result))
-  :documentation "The scaled Airy function @math{S_B(x) Bi(x)}.
-   For @math{x>0} the scaling factor @math{S_B(x)} is @math{exp(-(2/3) x^(3/2))},
-   and is 1 for @math{x<0}.")
+  :documentation			; FDL
+  "The scaled Airy function S_B(x) Bi(x).
+   For x>0 the scaling factor S_B(x) is exp(-(2/3) x^(3/2)),
+   and is 1 for x<0.")
 
-(defun-gsl airy-Ai-deriv (x)
+(defmfun airy-Ai-deriv (x)
   "gsl_sf_airy_Ai_deriv_e" ((x :double) :mode (ret sf-result))
-  :documentation "The Airy function derivative Ai'(x).")
+  :documentation			; FDL
+  "The Airy function derivative Ai'(x).")
 
-(defun-gsl airy-Bi-deriv (x)
+(defmfun airy-Bi-deriv (x)
   "gsl_sf_airy_Bi_deriv_e" ((x :double) :mode (ret sf-result))
   :documentation "The Airy function derivative Bi'(x).")
 
-(defun-gsl airy-Ai-deriv-scaled (x)
+(defmfun airy-Ai-deriv-scaled (x)
   "gsl_sf_airy_Ai_deriv_scaled_e" ((x :double) :mode (ret sf-result))
-  :documentation "The scaled Airy function derivative S_A(x) Ai'(x).
-  For @math{x>0} the scaling factor @math{S_A(x)} is @math{\exp(+(2/3) x^(3/2))},
-  and is 1 for @math{x<0}.")
+  :documentation			; FDL
+  "The scaled Airy function derivative S_A(x) Ai'(x).
+  For x>0 the scaling factor S_A(x) is exp(+(2/3) x^(3/2)),
+  and is 1 for x<0.")
 
-(defun-gsl airy-Bi-deriv-scaled (x)
+(defmfun airy-Bi-deriv-scaled (x)
   "gsl_sf_airy_Bi_deriv_scaled_e" ((x :double) :mode (ret sf-result))
-  :documentation "The scaled Airy function derivative S_B(x) Bi'(x).
-   For @math{x>0} the scaling factor @math{S_B(x)} is
-   @math{exp(-(2/3) x^(3/2))}, and is 1 for @math{x<0}.")
-
-(defun-gsl airy-zero-Ai (s)
-  "gsl_sf_airy_zero_Ai_e" ((s :size) (ret sf-result))
-  :documentation "The location of the @var{s}-th zero of the Airy
-  function @math{Ai(x)}.")
-
-(defun-gsl airy-zero-Bi (s)
-  "gsl_sf_airy_zero_Bi_e" ((s :size) (ret sf-result))
-  :documentation "The location of the @var{s}-th zero of the Airy
-   function @math{Bi(x)}.")
-
-(defun-gsl airy-zero-Ai-deriv (s)
-  "gsl_sf_airy_zero_Ai_deriv_e" ((s :size) (ret sf-result))
-  :documentation "The location of the @var{s}-th zero of the Airy
-   function derivative @math{Ai'(x)}.")
-
-(defun-gsl airy-zero-Bi-deriv (s)
-  "gsl_sf_airy_zero_Bi_deriv_e" ((s :size) (ret sf-result))
-  :documentation "The location of the @var{s}-th zero of the Airy
-   function derivative @math{Bi'(x)}.")
+  :documentation			; FDL
+  "The scaled Airy function derivative S_B(x) Bi'(x).
+   For x>0 the scaling factor S_B(x) is
+   exp(-(2/3) x^(3/2)), and is 1 for x<0.")
+
+(defmfun airy-zero-Ai (s)
+  "gsl_sf_airy_zero_Ai_e" ((s size) (ret sf-result))
+  :documentation			; FDL
+  "The location of the s-th zero of the Airy function Ai(x).")
+
+(defmfun airy-zero-Bi (s)
+  "gsl_sf_airy_zero_Bi_e" ((s size) (ret sf-result))
+  :documentation			; FDL
+  "The location of the s-th zero of the Airy function Bi(x).")
+
+(defmfun airy-zero-Ai-deriv (s)
+  "gsl_sf_airy_zero_Ai_deriv_e" ((s size) (ret sf-result))
+  :documentation			; FDL
+  "The location of the s-th zero of the Airy function derivative Ai'(x).")
+
+(defmfun airy-zero-Bi-deriv (s)
+  "gsl_sf_airy_zero_Bi_deriv_e" ((s size) (ret sf-result))
+  :documentation			; FDL
+  "The location of the s-th zero of the Airy function derivative Bi'(x).")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test airy
-  (lisp-unit:assert-first-fp-equal "0.157259233805d-01" (airy-ai 2.5d0))
-  (lisp-unit:assert-first-fp-equal "0.648166073846d+01" (airy-bi 2.5d0))
-  (lisp-unit:assert-first-fp-equal "0.219322205129d+00" (airy-ai-scaled 2.5d0))
-  (lisp-unit:assert-first-fp-equal "0.464750480196d+00" (airy-bi-scaled 2.5d0))
-  (lisp-unit:assert-first-fp-equal "-0.262508810359d-01" (airy-ai-deriv 2.5d0))
-  (lisp-unit:assert-first-fp-equal "0.942142331733d+01" (airy-bi-deriv 2.5d0))
-  (lisp-unit:assert-first-fp-equal "-0.366108938475d+00"
-				   (airy-ai-deriv-scaled 2.5d0))
-  (lisp-unit:assert-first-fp-equal "0.675538444164d+00"
-				   (airy-bi-deriv-scaled 2.5d0))
-  (lisp-unit:assert-first-fp-equal "-0.233810741046d+01" (airy-zero-ai 1))
-  (lisp-unit:assert-first-fp-equal "-0.117371322271d+01" (airy-zero-bi 1))
-  (lisp-unit:assert-first-fp-equal "-0.101879297165d+01" (airy-zero-ai-deriv 1))
-  (lisp-unit:assert-first-fp-equal "-0.229443968261d+01" (airy-zero-bi-deriv 1)))
+#|
+(make-tests airy
+  (airy-ai 2.5d0)
+  (airy-bi 2.5d0)
+  (airy-ai-scaled 2.5d0)
+  (airy-bi-scaled 2.5d0)
+  (airy-ai-deriv 2.5d0)
+  (airy-bi-deriv 2.5d0)
+  (airy-ai-deriv-scaled 2.5d0)
+  (airy-bi-deriv-scaled 2.5d0)
+  (airy-zero-ai 1)
+  (airy-zero-bi 1)
+  (airy-zero-ai-deriv 1)
+  (airy-zero-bi-deriv 1))
+|#
+
+(LISP-UNIT:DEFINE-TEST AIRY
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(0.01572592338047048d0 2.1573014423586447d-17)
+   (MULTIPLE-VALUE-LIST (AIRY-AI 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(6.48166073846058d0 8.77836191730236d-15)
+   (MULTIPLE-VALUE-LIST (AIRY-BI 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(0.21932220512871203d0 5.966864198455728d-17)
+   (MULTIPLE-VALUE-LIST (AIRY-AI-SCALED 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(0.4647504801960925d0 1.1831869152362144d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-BI-SCALED 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(-0.02625088103590322d0 4.306971270221159d-17)
+   (MULTIPLE-VALUE-LIST (AIRY-AI-DERIV 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(9.421423317334305d0 1.5213125884867257d-14)
+   (MULTIPLE-VALUE-LIST (AIRY-BI-DERIV 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(-0.36610893847516224d0 1.167515239400716d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-AI-DERIV-SCALED 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(0.6755384441644995d0 1.978922049880242d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-BI-DERIV-SCALED 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(-2.338107410459767d0 5.19164136227827d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-ZERO-AI 1)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(-1.173713222709128d0 2.606166888317336d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-ZERO-BI 1)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(-1.018792971647471d0 2.2621748288986134d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-ZERO-AI-DERIV 1)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   '(-2.294439682614123d0 5.094679528503672d-16)
+   (MULTIPLE-VALUE-LIST (AIRY-ZERO-BI-DERIV 1))))
 
 
 #|
diff --git a/special-functions/bessel.lisp b/special-functions/bessel.lisp
index 1adcd96ae94dcb8c6402148ca5bfd41c9febb7ce..c1db0afa9ef806a9eacf3b31911036945118251a 100644
--- a/special-functions/bessel.lisp
+++ b/special-functions/bessel.lisp
@@ -1,6 +1,6 @@
 ;; Bessel functions
 ;; Liam Healy, Fri Mar 17 2006 - 18:42
-;; Time-stamp: <2008-02-03 21:04:32EST bessel.lisp>
+;; Time-stamp: <2008-02-16 19:16:11EST bessel.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -9,23 +9,23 @@
 ;;;; Regular Cylindrical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl cylindrical-bessel-J0 (x)
+(defmfun cylindrical-bessel-J0 (x)
     "gsl_sf_bessel_J0_e"
   ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular cylindrical Bessel function of zeroth order, J_0(x).")
 
-(defun-gsl cylindrical-bessel-J1 (x)
+(defmfun cylindrical-bessel-J1 (x)
   "gsl_sf_bessel_J1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular cylindrical Bessel function of first order, J_1(x).")
 
-(defun-gsl cylindrical-bessel-Jn (n x)
+(defmfun cylindrical-bessel-Jn (n x)
   "gsl_sf_bessel_Jn_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The regular cylindrical Bessel function of order n, J_n(x).")
 
-(defun-gsl cylindrical-bessel-Jn-array (x array &optional (nmin 0))
+(defmfun cylindrical-bessel-Jn-array (x array &optional (nmin 0))
   "gsl_sf_bessel_Jn_array"
   ((nmin :int) ((+ nmin (dim0 array) -1) :int) (x :double)
    ((gsl-array array) :pointer))
@@ -40,23 +40,23 @@
 ;;;; Irregular Cylindrical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl cylindrical-bessel-Y0 (x)
+(defmfun cylindrical-bessel-Y0 (x)
   "gsl_sf_bessel_Y0_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular cylindrical Bessel function of zeroth order, Y_0(x).")
 
-(defun-gsl cylindrical-bessel-Y1 (x)
+(defmfun cylindrical-bessel-Y1 (x)
   "gsl_sf_bessel_Y1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular cylindrical Bessel function of first order, Y_1(x).")
 
-(defun-gsl cylindrical-bessel-Yn (n x)
+(defmfun cylindrical-bessel-Yn (n x)
   "gsl_sf_bessel_Yn_e"
   ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular cylindrical Bessel function of order n, Y_n(x).")
 
-(defun-gsl cylindrical-bessel-Yn-array (x array &optional (nmin 0))
+(defmfun cylindrical-bessel-Yn-array (x array &optional (nmin 0))
   "gsl_sf_bessel_Yn_array"
   ((nmin :int) ((+ nmin (dim0 array) -1) :int) (x :double)
    ((gsl-array array) :pointer))
@@ -72,22 +72,22 @@
 ;;;; Regular Modified Cylindrical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl cylindrical-bessel-I0 (x)
+(defmfun cylindrical-bessel-I0 (x)
   "gsl_sf_bessel_I0_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular modified cylindrical Bessel function of zeroth order, I_0(x).")
 
-(defun-gsl cylindrical-bessel-I1 (x)
+(defmfun cylindrical-bessel-I1 (x)
   "gsl_sf_bessel_I1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular modified cylindrical Bessel function of first order, I_1(x).")
 
-(defun-gsl cylindrical-bessel-In (n x)
+(defmfun cylindrical-bessel-In (n x)
   "gsl_sf_bessel_In_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The regular modified cylindrical Bessel function of order n, I_n(x).")
 
-(defun-gsl cylindrical-bessel-In-array (x array &optional (nmin 0))
+(defmfun cylindrical-bessel-In-array (x array &optional (nmin 0))
   "gsl_sf_bessel_In_array"
   ((nmin :int) ((+ nmin (dim0 array) -1) :int) (x :double)
    ((gsl-array array) :pointer))
@@ -98,25 +98,25 @@
    The values are computed using recurrence relations for efficiency, and
    therefore may differ slightly from the exact values.")
 
-(defun-gsl cylindrical-bessel-I0-scaled (x)
+(defmfun cylindrical-bessel-I0-scaled (x)
   "gsl_sf_bessel_I0_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified cylindrical Bessel function of zeroth order,
   \exp(-|x|) I_0(x).")
 
-(defun-gsl cylindrical-bessel-I1-scaled (x)
+(defmfun cylindrical-bessel-I1-scaled (x)
   "gsl_sf_bessel_I1_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified cylindrical Bessel function of first order,
   \exp(-|x|) I_1(x).")
 
-(defun-gsl cylindrical-bessel-In-scaled (n x)
+(defmfun cylindrical-bessel-In-scaled (n x)
   "gsl_sf_bessel_In_scaled_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified cylindrical Bessel function of order n,
   \exp(-|x|) I_n(x)}.")
 
-(defun-gsl cylindrical-bessel-In-scaled-array (x array &optional (nmin 0))
+(defmfun cylindrical-bessel-In-scaled-array (x array &optional (nmin 0))
   "gsl_sf_bessel_In_scaled_array"
   ((nmin :int) ((+ nmin (dim0 array) -1) :int) (x :double)
    ((gsl-array array) :pointer))
@@ -133,23 +133,23 @@
 ;;;; Irregular Modified Cylindrical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl cylindrical-bessel-K0 (x)
+(defmfun cylindrical-bessel-K0 (x)
   "gsl_sf_bessel_K0_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular modified cylindrical Bessel function of zeroth order,
   K_0(x).")
 
-(defun-gsl cylindrical-bessel-K1 (x)
+(defmfun cylindrical-bessel-K1 (x)
   "gsl_sf_bessel_K1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular modified cylindrical Bessel function of first order, K_1(x).")
 
-(defun-gsl cylindrical-bessel-Kn (n x)
+(defmfun cylindrical-bessel-Kn (n x)
   "gsl_sf_bessel_Kn_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular modified cylindrical Bessel function of order n, K_n(x).")
 
-(defun-gsl cylindrical-bessel-Kn-array (x array &optional (nmin 0))
+(defmfun cylindrical-bessel-Kn-array (x array &optional (nmin 0))
   "gsl_sf_bessel_Kn_array"
   ((nmin :int) ((+ nmin (dim0 array) -1) :int) (x :double)
    ((gsl-array array) :pointer))
@@ -160,25 +160,25 @@
    The values are computed using recurrence relations for efficiency, and
    therefore may differ slightly from the exact values.")
 
-(defun-gsl cylindrical-bessel-K0-scaled (x)
+(defmfun cylindrical-bessel-K0-scaled (x)
   "gsl_sf_bessel_K0_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified cylindrical Bessel function of zeroth order,
   \exp(-|x|) K_0(x).")
 
-(defun-gsl cylindrical-bessel-K1-scaled (x)
+(defmfun cylindrical-bessel-K1-scaled (x)
   "gsl_sf_bessel_K1_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified cylindrical Bessel function of first order,
    \exp(-|x|) K_1(x).")
 
-(defun-gsl cylindrical-bessel-Kn-scaled (n x)
+(defmfun cylindrical-bessel-Kn-scaled (n x)
   "gsl_sf_bessel_Kn_scaled_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified cylindrical Bessel function of order n,
   \exp(-|x|) K_n(x).")
 
-(defun-gsl cylindrical-bessel-Kn-scaled-array (x array &optional (nmin 0))
+(defmfun cylindrical-bessel-Kn-scaled-array (x array &optional (nmin 0))
   "gsl_sf_bessel_Kn_scaled_array"
   ((nmin :int) ((+ nmin (dim0 array) -1) :int) (x :double)
    ((gsl-array array) :pointer))
@@ -196,30 +196,30 @@
 ;;;; Regular Spherical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl spherical-bessel-j0 (x)
+(defmfun spherical-bessel-j0 (x)
   "gsl_sf_bessel_j0_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular spherical Bessel function of zeroth order, j_0(x) = \sin(x)/x.")
 
-(defun-gsl spherical-bessel-j1 (x)
+(defmfun spherical-bessel-j1 (x)
   "gsl_sf_bessel_j1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular spherical Bessel function of first order, j_1(x)
    = (\sin(x)/x - \cos(x))/x.")
 
-(defun-gsl spherical-bessel-j2 (x)
+(defmfun spherical-bessel-j2 (x)
   "gsl_sf_bessel_j2_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The regular spherical Bessel function of second order, j_2(x)
    = ((3/x^2 - 1)\sin(x) - 3\cos(x)/x)/x.")
 
-(defun-gsl spherical-bessel-jl (l x)
+(defmfun spherical-bessel-jl (l x)
   "gsl_sf_bessel_jl_e" ((l :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The regular spherical Bessel function of order
    l, j_l(x), for l >= 0 and x >= 0.")
 
-(defun-gsl spherical-bessel-jl-array (x array)
+(defmfun spherical-bessel-jl-array (x array)
   "gsl_sf_bessel_jl_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :invalidate (array)
@@ -229,7 +229,7 @@
   The values are computed using recurrence relations for
   efficiency, and therefore may differ slightly from the exact values.")
 
-(defun-gsl spherical-bessel-jl-steed-array (x array)
+(defmfun spherical-bessel-jl-steed-array (x array)
   "gsl_sf_bessel_jl_steed_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :invalidate (array)
@@ -245,30 +245,30 @@
 ;;;; Irregular Spherical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl spherical-bessel-y0 (x)
+(defmfun spherical-bessel-y0 (x)
   "gsl_sf_bessel_y0_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular spherical Bessel function of zeroth order,
   y_0(x) = -\cos(x)/x.")
 
-(defun-gsl spherical-bessel-y1 (x)
+(defmfun spherical-bessel-y1 (x)
   "gsl_sf_bessel_y1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular spherical Bessel function of first order,
   y_1(x) = -(\cos(x)/x + \sin(x))/x.")
 
-(defun-gsl spherical-bessel-y2 (x)
+(defmfun spherical-bessel-y2 (x)
   "gsl_sf_bessel_y2_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular spherical Bessel function of second order,
   y_2(x) = (-3/x^3 + 1/x)\cos(x) - (3/x^2)\sin(x).")
 
-(defun-gsl spherical-bessel-yl (l x)
+(defmfun spherical-bessel-yl (l x)
   "gsl_sf_bessel_yl_e" ((l :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular spherical Bessel function of order l, y_l(x), for l >= 0.")
 
-(defun-gsl spherical-bessel-yl-array (x array)
+(defmfun spherical-bessel-yl-array (x array)
   "gsl_sf_bessel_yl_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :invalidate (array)
@@ -282,31 +282,31 @@
 ;;;; Regular Modified Spherical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl spherical-bessel-i0-scaled (x)
+(defmfun spherical-bessel-i0-scaled (x)
   "gsl_sf_bessel_i0_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified spherical Bessel function of zeroth
   order, \exp(-|x|) i_0(x).")
 
-(defun-gsl spherical-bessel-i1-scaled (x)
+(defmfun spherical-bessel-i1-scaled (x)
   "gsl_sf_bessel_i1_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified spherical Bessel function of first order,
   \exp(-|x|) i_1(x).")
 
-(defun-gsl spherical-bessel-i2-scaled (x)
+(defmfun spherical-bessel-i2-scaled (x)
   "gsl_sf_bessel_i2_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified spherical Bessel function of second order,
    \exp(-|x|) i_2(x).")
 
-(defun-gsl spherical-bessel-il-scaled (n x)
+(defmfun spherical-bessel-il-scaled (n x)
   "gsl_sf_bessel_il_scaled_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified spherical Bessel function of order l,
    \exp(-|x|) i_l(x).")
 
-(defun-gsl spherical-bessel-il-scaled-array (x array)
+(defmfun spherical-bessel-il-scaled-array (x array)
   "gsl_sf_bessel_il_scaled_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :invalidate (array)
@@ -320,31 +320,31 @@
 ;;;; Irregular Modified Spherical Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl spherical-bessel-k0-scaled (x)
+(defmfun spherical-bessel-k0-scaled (x)
   "gsl_sf_bessel_k0_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified spherical Bessel function of zeroth
   order, \exp(x) k_0(x), for x>0.")
 
-(defun-gsl spherical-bessel-k1-scaled (x)
+(defmfun spherical-bessel-k1-scaled (x)
   "gsl_sf_bessel_k1_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified spherical Bessel function of first order,
    \exp(x) k_1(x), for x>0.")
 
-(defun-gsl spherical-bessel-k2-scaled (x)
+(defmfun spherical-bessel-k2-scaled (x)
   "gsl_sf_bessel_k2_scaled_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified spherical Bessel function of second order,
   \exp(x) k_2(x), for x>0.")
 
-(defun-gsl spherical-bessel-kl-scaled (n x)
+(defmfun spherical-bessel-kl-scaled (n x)
   "gsl_sf_bessel_il_scaled_e" ((n :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified spherical Bessel function of order l,
    \exp(x) k_l(x), for x>0.")
 
-(defun-gsl spherical-bessel-kl-scaled-array (x array)
+(defmfun spherical-bessel-kl-scaled-array (x array)
   "gsl_sf_bessel_kl_scaled_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :invalidate (array)
@@ -359,13 +359,13 @@
 ;;;; Regular Bessel Function - Fractional Order
 ;;;;****************************************************************************
 
-(defun-gsl bessel-Jnu (nu x)
+(defmfun bessel-Jnu (nu x)
   "gsl_sf_bessel_Jnu_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The regular cylindrical Bessel function of fractional order
   \nu, J_\nu(x).")
 
-(defun-gsl spherical-Jnu-array (nu v)
+(defmfun spherical-Jnu-array (nu v)
   "gsl_sf_bessel_sequence_Jnu_e"
   ((nu :double) :mode ((dim0 v) :int) ((gsl-array v) :pointer))
   :invalidate (v)
@@ -380,7 +380,7 @@
 ;;;; Irregular Bessel Function - Fractional Order
 ;;;;****************************************************************************
 
-(defun-gsl bessel-Ynu (nu x)
+(defmfun bessel-Ynu (nu x)
   "gsl_sf_bessel_Ynu_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular cylindrical Bessel function of fractional order
@@ -390,13 +390,13 @@
 ;;;; Regular Modified Bessel Functions - Fractional Order
 ;;;;****************************************************************************
 
-(defun-gsl bessel-Inu (nu x)
+(defmfun bessel-Inu (nu x)
   "gsl_sf_bessel_Inu_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The regular modified Bessel function of fractional order
   \nu, I_\nu(x) for x>0, \nu>0.")
 
-(defun-gsl bessel-Inu-scaled (nu x)
+(defmfun bessel-Inu-scaled (nu x)
   "gsl_sf_bessel_Inu_scaled_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled regular modified Bessel function of fractional order
@@ -406,19 +406,19 @@
 ;;;; Irregular Modified Bessel Functions - Fractional Order
 ;;;;****************************************************************************
 
-(defun-gsl bessel-Knu (nu x)
+(defmfun bessel-Knu (nu x)
   "gsl_sf_bessel_Knu_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular modified Bessel function of fractional order \nu,
    K_\nu(x) for x>0, \nu>0.")
 
-(defun-gsl bessel-lnKnu (nu x)
+(defmfun bessel-lnKnu (nu x)
   "gsl_sf_bessel_lnKnu_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The logarithm of the irregular modified Bessel function of fractional
    order \nu, \ln(K_\nu(x)) for x>0, \nu>0.")
 
-(defun-gsl bessel-Knu-scaled (nu x)
+(defmfun bessel-Knu-scaled (nu x)
   "gsl_sf_bessel_Knu_scaled_e" ((nu :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The scaled irregular modified Bessel function of fractional order
@@ -428,19 +428,19 @@
 ;;;; Zeros of Regular Bessel Functions
 ;;;;****************************************************************************
 
-(defun-gsl bessel-zero-J0 (s)
+(defmfun bessel-zero-J0 (s)
   "gsl_sf_bessel_zero_J0_e" ((s :int) (ret sf-result))
   :documentation			; FDL
   "The location of the s-th positive zero of the Bessel function
   J_0(x).")
 
-(defun-gsl bessel-zero-J1 (s)
+(defmfun bessel-zero-J1 (s)
   "gsl_sf_bessel_zero_J1_e" ((s :int) (ret sf-result))
   :documentation			; FDL
   "The location of the s-th positive zero of the Bessel function
   J_1(x).")
 
-(defun-gsl bessel-zero-Jnu (nu s)
+(defmfun bessel-zero-Jnu (nu s)
   "gsl_sf_bessel_zero_Jnu_e" ((nu :double) (s :int) (ret sf-result))
   :documentation			; FDL
   "These routines compute the location of the s-th positive zero
@@ -451,219 +451,341 @@
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test bessel
-  (lisp-unit:assert-first-fp-equal
-   "-0.397149809864d+00"
-   (cylindrical-bessel-J0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.660433280235d-01"
-   (cylindrical-bessel-J1 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.364128145852d+00"
-   (cylindrical-bessel-Jn 2 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.352834028616d+00" "0.128943249474d+00" "0.339957198076d-01"
-     "0.703962975587d-02")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (cylindrical-bessel-Jn-array 2.0d0 besarr 2)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "-0.169407393251d-01"
-   (cylindrical-bessel-Y0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.397925710557d+00"
-   (cylindrical-bessel-Y1 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.182022115953d+00"
-   (cylindrical-bessel-Yn 3 4.0d0))
-  (lisp-unit:assert-equal
-   '("-0.617408104191d+00" "-0.112778377684d+01" "-0.276594322633d+01"
-     "-0.993598912848d+01")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (cylindrical-bessel-Yn-array 2.0d0 besarr 2)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.113019219521d+02"
-   (cylindrical-bessel-I0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.975946515370d+01"
-   (cylindrical-bessel-I1 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.333727577842d+01"
-   (cylindrical-bessel-In 3 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.688948447699d+00" "0.212739959240d+00" "0.507285699792d-01"
-     "0.982567932313d-02")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (cylindrical-bessel-In-array 2.0d0 besarr 2)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.207001921224d+00"
-   (cylindrical-bessel-I0-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.178750839502d+00"
-   (cylindrical-bessel-I1-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.611243380297d-01"
-   (cylindrical-bessel-In-scaled 3 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.932390333047d-01" "0.287912226395d-01" "0.686536538632d-02"
-     "0.132976109419d-02")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (cylindrical-bessel-In-scaled-array 2.0d0 besarr 2)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.111596760859d-01"
-   (cylindrical-bessel-K0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.124834988873d-01"
-   (cylindrical-bessel-K1 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.174014255295d-01"
-   (cylindrical-bessel-Kn 2 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.253759754566d+00" "0.647385390949d+00" "0.219591592741d+01"
-     "0.943104910060d+01")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (cylindrical-bessel-Kn-array 2.0d0 besarr 2)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.609297669257d+00"
-   (cylindrical-bessel-K0-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.681575945186d+00"
-   (cylindrical-bessel-K1-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.950085641850d+00"
-   (cylindrical-bessel-Kn-scaled 2 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.253759754566d+00" "0.647385390949d+00" "0.219591592741d+01"
-     "0.943104910060d+01")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (cylindrical-bessel-Kn-array 2.0d0 besarr 2)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "-0.189200623827d+00"
-   (spherical-bessel-j0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.116110749259d+00"
-   (spherical-bessel-j1 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.276283685771d+00"
-   (spherical-bessel-j2 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.229243857955d+00"
-   (spherical-bessel-jl 3 4.0d0))
-  (lisp-unit:assert-equal
-   '("-0.189200623827d+00" "0.116110749259d+00" "0.276283685771d+00"
-     "0.229243857955d+00")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (spherical-bessel-jl-array 4.0d0 besarr)
-      (data besarr))))
-  (lisp-unit:assert-equal
-   '("-0.189200623827d+00" "0.116110749259d+00" "0.276283685771d+00"
-     "0.229243857955d+00")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (spherical-bessel-jl-steed-array 4.0d0 besarr)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.163410905216d+00"
-   (spherical-bessel-y0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.230053350131d+00"
-   (spherical-bessel-y1 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.912910738232d-02"
-   (spherical-bessel-y2 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.912910738232d-02"
-   (spherical-bessel-yl 2 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.163410905216d+00" "0.230053350131d+00" "0.912910738232d-02"
-     "-0.218641965903d+00")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (spherical-bessel-yl-array 4.0d0 besarr)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.124958067172d+00"
-   (spherical-bessel-i0-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.938024160356d-01"
-   (spherical-bessel-i1-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.546062551448d-01"
-   (spherical-bessel-i2-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.255445971046d-01"
-   (spherical-bessel-il-scaled 3 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.124958067172d+00" "0.938024160356d-01" "0.546062551448d-01"
-     "0.255445971046d-01")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (spherical-bessel-il-scaled-array 4.0d0 besarr)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.392699081699d+00"
-   (spherical-bessel-k0-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.490873852123d+00"
-   (spherical-bessel-k1-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.760854470791d+00"
-   (spherical-bessel-k2-scaled 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.255445971046d-01"
-   (spherical-bessel-kl-scaled 3 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.392699081699d+00" "0.490873852123d+00" "0.760854470791d+00"
-     "0.144194194061d+01")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double 4)))
-      (spherical-bessel-kl-scaled-array 4.0d0 besarr)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.430171473876d+00"
-   (bessel-jnu 3.0d0 4.0d0))
-  (lisp-unit:assert-equal
-   '("0.671396707142d+00" "0.513016136562d+00" "0.650081828774d-01")
-   (lisp-unit:fp-sequence
-    (letm ((besarr (vector-double #(1.0d0 2.0d0 3.0d0))))
-      (spherical-Jnu-array 0.5d0 besarr)
-      (data besarr))))
-  (lisp-unit:assert-first-fp-equal
-   "-0.182022115953d+00"
-   (bessel-Ynu 3.0d0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.333727577842d+01"
-   (bessel-Inu 3.0d0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.611243380297d-01"
-   (bessel-Inu-scaled 3.0d0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.298849244168d-01"
-   (bessel-Knu 3.0d0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.351040112585d+01"
-   (bessel-lnKnu 3.0d0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.163166158704d+01"
-   (bessel-Knu-scaled 3.0d0 4.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.149309177085d+02"
-   (bessel-zero-J0 5))
-  (lisp-unit:assert-first-fp-equal
-   "0.164706300509d+02"
-   (bessel-zero-J1 5))
-  (lisp-unit:assert-first-fp-equal
-   "0.179598194950d+02"
-   (bessel-zero-Jnu 2.0d0 5)))
+#|
+
+;; This bizarrely causes a division-by-zero error 
+;;	    (bessel-jnu 3.0d0 4.0d0)
+;; in macroexpansion the form from within SLIME using SBCL (expands
+;; fine from shell)
+
+(make-tests bessel
+	    (cylindrical-bessel-J0 4.0d0)
+	    (cylindrical-bessel-J1 4.0d0)
+	    (cylindrical-bessel-Jn 2 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (cylindrical-bessel-Jn-array 2.0d0 besarr 2)
+	      (data besarr))
+	    (cylindrical-bessel-Y0 4.0d0)
+	    (cylindrical-bessel-Y1 4.0d0)
+	    (cylindrical-bessel-Yn 3 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (cylindrical-bessel-Yn-array 2.0d0 besarr 2)
+	      (data besarr))
+	    (cylindrical-bessel-I0 4.0d0)
+	    (cylindrical-bessel-I1 4.0d0)
+	    (cylindrical-bessel-In 3 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (cylindrical-bessel-In-array 2.0d0 besarr 2)
+	      (data besarr))
+	    (cylindrical-bessel-I0-scaled 4.0d0)
+	    (cylindrical-bessel-I1-scaled 4.0d0)
+	    (cylindrical-bessel-In-scaled 3 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (cylindrical-bessel-In-scaled-array 2.0d0 besarr 2)
+	      (data besarr))
+	    (cylindrical-bessel-K0 4.0d0)
+	    (cylindrical-bessel-K1 4.0d0)
+	    (cylindrical-bessel-Kn 2 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (cylindrical-bessel-Kn-array 2.0d0 besarr 2)
+	      (data besarr))
+	    (cylindrical-bessel-K0-scaled 4.0d0)
+	    (cylindrical-bessel-K1-scaled 4.0d0)
+	    (cylindrical-bessel-Kn-scaled 2 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (cylindrical-bessel-Kn-array 2.0d0 besarr 2)
+	      (data besarr))
+	    (spherical-bessel-j0 4.0d0)
+	    (spherical-bessel-j1 4.0d0)
+	    (spherical-bessel-j2 4.0d0)
+	    (spherical-bessel-jl 3 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (spherical-bessel-jl-array 4.0d0 besarr)
+	      (data besarr))
+	    (letm ((besarr (vector-double 4)))
+	      (spherical-bessel-jl-steed-array 4.0d0 besarr)
+	      (data besarr))
+	    (spherical-bessel-y0 4.0d0)
+	    (spherical-bessel-y1 4.0d0)
+	    (spherical-bessel-y2 4.0d0)
+	    (spherical-bessel-yl 2 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (spherical-bessel-yl-array 4.0d0 besarr)
+	      (data besarr))
+	    (spherical-bessel-i0-scaled 4.0d0)
+	    (spherical-bessel-i1-scaled 4.0d0)
+	    (spherical-bessel-i2-scaled 4.0d0)
+	    (spherical-bessel-il-scaled 3 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (spherical-bessel-il-scaled-array 4.0d0 besarr)
+	      (data besarr))
+	    (spherical-bessel-k0-scaled 4.0d0)
+	    (spherical-bessel-k1-scaled 4.0d0)
+
+	    (spherical-bessel-k2-scaled 4.0d0)
+	    (spherical-bessel-kl-scaled 3 4.0d0)
+	    (letm ((besarr (vector-double 4)))
+	      (spherical-bessel-kl-scaled-array 4.0d0 besarr)
+	      (data besarr))
+	    (bessel-jnu 3.0d0 4.0d0)
+	    (letm ((besarr (vector-double #(1.0d0 2.0d0 3.0d0))))
+	      (spherical-Jnu-array 0.5d0 besarr)
+	      (data besarr))
+	    (bessel-Ynu 3.0d0 4.0d0)
+	    (bessel-Inu 3.0d0 4.0d0)
+	    (bessel-Inu-scaled 3.0d0 4.0d0)
+	    (bessel-Knu 3.0d0 4.0d0)
+	    (bessel-lnKnu 3.0d0 4.0d0)
+	    (bessel-Knu-scaled 3.0d0 4.0d0)
+	    (bessel-zero-J0 5)
+	    (bessel-zero-J1 5)
+	    (bessel-zero-Jnu 2.0d0 5))
+
+|#
+
+(LISP-UNIT:DEFINE-TEST BESSEL
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.3971498098638474d0 4.334456411751256d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-J0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.0660433280235491d0 2.1409770694795335d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-J1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.3641281458520729d0 3.974061014982464d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-JN 2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.35283402861563773d0 0.12894324947440206d0
+      0.033995719807568436d0 0.007039629755871686d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (CYLINDRICAL-BESSEL-JN-ARRAY 2.0d0 BESARR 2)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.016940739325064968d0 1.8993556609468549d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-Y0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.3979257105570999d0 3.1396236150465943d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-Y1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.18202211595348539d0 3.355735727760045d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-YN 3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(-0.6174081041906827d0 -1.127783776840428d0
+      -2.7659432263306014d0 -9.935989128481978d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (CYLINDRICAL-BESSEL-YN-ARRAY 2.0d0 BESARR 2)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 11.301921952136329d0 2.7297681442535893d-14)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-I0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 9.759465153704449d0 1.9210136786427457d-14)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-I1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.3372757784203446d0 8.06056628872663d-15)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-IN 3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.6889484476987382d0 0.21273995923985267d0
+      0.05072856997918024d0 0.009825679323131702d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (CYLINDRICAL-BESSEL-IN-ARRAY 2.0d0 BESARR 2)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.2070019212239867d0 2.241925168997723d-16)
+   (MULTIPLE-VALUE-LIST
+    (CYLINDRICAL-BESSEL-I0-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1787508395024353d0 1.1370197115937822d-16)
+   (MULTIPLE-VALUE-LIST
+    (CYLINDRICAL-BESSEL-I1-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0611243380296663d0 9.334510342661594d-17)
+   (MULTIPLE-VALUE-LIST
+    (CYLINDRICAL-BESSEL-IN-SCALED 3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.09323903330473338d0 0.028791222639470898d0
+      0.006865365386320685d0 0.0013297610941881578d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (CYLINDRICAL-BESSEL-IN-SCALED-ARRAY 2.0d0 BESARR
+					  2)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.011159676085853023d0 2.0424662435034432d-17)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-K0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.01248349888726843d0 1.767412161819488d-17)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-K1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.017401425529487143d0 2.257461693414273d-16)
+   (MULTIPLE-VALUE-LIST (CYLINDRICAL-BESSEL-KN 2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.2537597545660558d0 0.6473853909486341d0
+      2.1959159274119586d0 9.431049100596468d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (CYLINDRICAL-BESSEL-KN-ARRAY 2.0d0 BESARR 2)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6092976692566953d0 3.0340122249326356d-16)
+   (MULTIPLE-VALUE-LIST
+    (CYLINDRICAL-BESSEL-K0-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.681575945185671d0 3.596132979136138d-16)
+   (MULTIPLE-VALUE-LIST
+    (CYLINDRICAL-BESSEL-K1-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9500856418495256d0 1.1481477659153143d-14)
+   (MULTIPLE-VALUE-LIST
+    (CYLINDRICAL-BESSEL-KN-SCALED 2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.2537597545660558d0 0.6473853909486341d0
+      2.1959159274119586d0 9.431049100596468d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (CYLINDRICAL-BESSEL-KN-ARRAY 2.0d0 BESARR 2)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.18920062382698205d0 1.6804391107692678d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-J0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.11611074925915747d0 2.387125482192573d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-J1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.27628368577135015d0 3.680838111259856d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-J2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.22924385795503022d0 4.581212568400321d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-JL 3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(-0.18920062382698202d0 0.11611074925915743d0
+      0.2762836857713501d0 0.22924385795503022d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (SPHERICAL-BESSEL-JL-ARRAY 4.0d0 BESARR)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(-0.18920062382698208d0 0.11611074925915742d0
+      0.27628368577135015d0 0.22924385795503024d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (SPHERICAL-BESSEL-JL-STEED-ARRAY 4.0d0 BESARR)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.16341090521590299d0 1.4513803955642766d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-Y0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.2300533501309578d0 1.5324631572452525d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-Y1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.009129107382315343d0 1.6876604955506113d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-Y2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.009129107382315343d0 1.6876604955506113d-16)
+   (MULTIPLE-VALUE-LIST (SPHERICAL-BESSEL-YL 2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.16341090521590299d0 0.2300533501309578d0
+      0.009129107382315343d0 -0.21864196590306362d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (SPHERICAL-BESSEL-YL-ARRAY 4.0d0 BESARR)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.12495806717151219d0 5.5492529314587895d-17)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-I0-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.09380241603560975d0 4.165664081928078d-17)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-I1-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.05460625514480487d0 2.425004870012731d-17)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-I2-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.02554459710460367d0 5.842201171222646d-16)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-IL-SCALED 3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.12495806717151212d0 0.09380241603560971d0
+      0.05460625514480483d0 0.02554459710460367d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (SPHERICAL-BESSEL-IL-SCALED-ARRAY 4.0d0 BESARR)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.39269908169872414d0 1.743934249004316d-16)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-K0-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4908738521234052d0 2.1799178112553949d-16)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-K1-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.760854470791278d0 3.378872607445862d-16)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-K2-SCALED 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.02554459710460367d0 5.842201171222646d-16)
+   (MULTIPLE-VALUE-LIST
+    (SPHERICAL-BESSEL-KL-SCALED 3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.39269908169872414d0 0.4908738521234052d0
+      0.760854470791278d0 1.4419419406125027d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE 4)))
+      (SPHERICAL-BESSEL-KL-SCALED-ARRAY 4.0d0 BESARR)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.43017147387562193d0 7.641380397338472d-16)
+   (MULTIPLE-VALUE-LIST (BESSEL-JNU 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.6713967071418024d0 0.5130161365618323d0
+      0.06500818287738516d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((BESARR (VECTOR-DOUBLE #(1.0d0 2.0d0 3.0d0))))
+      (SPHERICAL-JNU-ARRAY 0.5d0 BESARR)
+      (DATA BESARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.1820221159534852d0 2.020851441225493d-15)
+   (MULTIPLE-VALUE-LIST (BESSEL-YNU 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.3372757784203437d0 1.1856385307923545d-14)
+   (MULTIPLE-VALUE-LIST (BESSEL-INU 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.061124338029666284d0 1.3572329489101316d-16)
+   (MULTIPLE-VALUE-LIST (BESSEL-INU-SCALED 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.029884924416755682d0 1.0617257976532701d-16)
+   (MULTIPLE-VALUE-LIST (BESSEL-KNU 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -3.5104011258456183d0 4.776268519767339d-15)
+   (MULTIPLE-VALUE-LIST (BESSEL-LNKNU 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.6316615870352025d0 5.072223134504136d-15)
+   (MULTIPLE-VALUE-LIST (BESSEL-KNU-SCALED 3.0d0 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 14.930917708487813d0 4.4792753125463437d-14)
+   (MULTIPLE-VALUE-LIST (BESSEL-ZERO-J0 5)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 16.470630050877624d0 3.2941260101755246d-13)
+   (MULTIPLE-VALUE-LIST (BESSEL-ZERO-J1 5)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 17.95981949498783d0 3.591963898997566d-14)
+   (MULTIPLE-VALUE-LIST (BESSEL-ZERO-JNU 2.0d0 5))))
+
+
diff --git a/special-functions/clausen.lisp b/special-functions/clausen.lisp
index 267fd4737ce50e2c017dc5302d4ef7753f57b767..8a432325719c7f08a7dd932f2269846b4054dfeb 100644
--- a/special-functions/clausen.lisp
+++ b/special-functions/clausen.lisp
@@ -1,18 +1,21 @@
-;********************************************************
-; file:        clausen.lisp                              
-; description: Clausen function                          
-; date:        Sat Mar 18 2006 - 23:18                   
-; author:      Liam M. Healy                             
-; modified:    Sun Jun 11 2006 - 21:03
-;********************************************************
-;;; $Id:$
+;; Clausen function
+;; Liam Healy, Sat Mar 18 2006 - 23:18
+;; Time-stamp: <2008-02-16 19:20:29EST clausen.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl clausen (x)
+(defmfun clausen (x)
   "gsl_sf_clausen_e" ((x :double) (ret sf-result))
-  :documentation
-  "The Clausen integral @math{Cl_2(x)}.")
+  :documentation			; FDL
+  "The Clausen integral Cl_2(x).")
 
-(lisp-unit:define-test clausen
-  (lisp-unit:assert-first-fp-equal "0.433598203236d+00" (clausen 2.5d0)))
+#|
+(make-tests clausen
+	    (clausen 2.5d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST CLAUSEN
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4335982032355329d0 1.2032502825912828d-15)
+   (MULTIPLE-VALUE-LIST (CLAUSEN 2.5d0))))
diff --git a/special-functions/coulomb.lisp b/special-functions/coulomb.lisp
index 50cce726e31f751d3912f555380cedd8b4cd2475..10f0f69f227639d22a8d34fd2f65f7a2eed0cde9 100644
--- a/special-functions/coulomb.lisp
+++ b/special-functions/coulomb.lisp
@@ -1,6 +1,6 @@
 ;; Coulumb functions
 ;; Liam Healy, Sat Mar 18 2006 - 23:23
-;; Time-stamp: <2008-02-03 19:18:56EST coulomb.lisp>
+;; Time-stamp: <2008-02-16 19:57:49EST coulomb.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -42,8 +42,8 @@
    (err F) (err Fp) (err G) (err Gp))
   :documentation			; FDL
   "The Coulomb wave functions F_L(\eta,x),
-  G_@{L-k@}(\eta,x) and their derivatives F'_L(\eta,x), 
-  G'_@{L-k@}(\eta,x) with respect to x.  The parameters are restricted to
+  G_{L-k}(\eta,x) and their derivatives F'_L(\eta,x), 
+  G'_{L-k}(\eta,x) with respect to x.  The parameters are restricted to
   L, L-k > -1/2}, x > 0 and integer k.  Note that L
   itself is not restricted to being an integer. The results are stored in
   the parameters F, G for the function values and Fp,
@@ -127,44 +127,83 @@
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test coulomb
-  (lisp-unit:assert-first-fp-equal
-   "0.164169997248d+00"
-   (hydrogenicr-1 1.0d0 2.5d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.766646808415d-01"
-   (hydrogenicr 3 1 1.0d0 2.5d0))
-  (lisp-unit:assert-equal
-   '("0.620350520114d-01" "0.177098574917d+00" "0.360501756616d+01"
-     "-0.582826184164d+01" "0.000000000000d+01" "0.000000000000d+01")
-   (lisp-unit:fp-sequence
-    (subseq (multiple-value-list (coulomb-wave-FG 0.0d0 1.0d0 2.0d0 0)) 0 6)))
-  (lisp-unit:assert-equal
-   '("0.661781613833d+00" "0.361412857745d+00" "0.132677576099d+00")
-   (lisp-unit:fp-sequence
-    (letm ((arr (vector-double 3)))
-      (coulomb-wave-F-array  0.0d0 1.0d0 2.0d0 arr)
-      (data arr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.716177996747d-01"
-   (coulomb-wave-fg 1.0d0 2.0d0 2.5d0 1))
-  (lisp-unit:assert-equal
-   '("0.350215846039d-01" "0.575250061420d-02" "0.711695560198d-03"
-     "0.647172649613d+01" "0.275745747216d+02" "0.170560372931d+03")
-   (lisp-unit:fp-sequence 
-    (letm ((Farr (vector-double 3)) (Garr (vector-double 3)))
-      (coulomb-wave-FG-array 1.5d0 1.0d0 1.0d0 Farr Garr)
-      (append (coerce (data Farr) 'list) (coerce (data Garr) 'list)))))
-  (lisp-unit:assert-equal
-   '("0.330890806916d+00" "0.180706428873d+00" "0.663387880496d-01")
-   (lisp-unit:fp-sequence
-    (letm ((arr (vector-double 3)))
-      (coulomb-wave-sphF-array  0.0d0 1.0d0 2.0d0 arr) (data arr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.138091464419d-02"
-   (coulomb-cl 1.0d0 2.5d0))
-  (lisp-unit:assert-equal
-   '("0.108422513102d+00" "0.511108628318d-01" "0.114287363681d-01")
-   (lisp-unit:fp-sequence
-    (letm ((cl (vector-double 3)))
-      (coulomb-CL-array 0.0d0 1.0d0 cl) (data cl)))))
+#|
+(make-tests coulomb
+  (hydrogenicr-1 1.0d0 2.5d0)
+  (hydrogenicr 3 1 1.0d0 2.5d0)
+  (coulomb-wave-FG 0.0d0 1.0d0 2.0d0 0)
+  (letm ((arr (vector-double 3)))
+    (coulomb-wave-F-array  0.0d0 1.0d0 2.0d0 arr)
+    (data arr))
+  (coulomb-wave-fg 1.0d0 2.0d0 2.5d0 1)
+  (letm ((Farr (vector-double 3)) (Garr (vector-double 3)))
+    (coulomb-wave-FG-array 1.5d0 1.0d0 1.0d0 Farr Garr)
+    (append (coerce (data Farr) 'list) (coerce (data Garr) 'list)))
+  (letm ((arr (vector-double 3)))
+    (coulomb-wave-sphF-array  0.0d0 1.0d0 2.0d0 arr) (data arr))
+  (coulomb-cl 1.0d0 2.5d0)
+  (letm ((cl (vector-double 3)))
+    (coulomb-CL-array 0.0d0 1.0d0 cl) (data cl)))
+|#
+
+(LISP-UNIT:DEFINE-TEST COULOMB
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1641699972477976d0 1.8226531089715347d-16)
+   (MULTIPLE-VALUE-LIST (HYDROGENICR-1 1.0d0 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07666468084146327d0 4.203054519905891d-16)
+   (MULTIPLE-VALUE-LIST (HYDROGENICR 3 1 1.0d0 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.06203505201137388d0 0.17709857491700906d0
+	 3.6050175661599693d0 -5.8282618416439025d0 0.0d0
+	 0.0d0 4.3793707124032857d-16
+	 1.2502291640824433d-15 3.5328525523867457d-14
+	 5.711592064490999d-14)
+   (MULTIPLE-VALUE-LIST
+    (COULOMB-WAVE-FG 0.0d0 1.0d0 2.0d0 0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.6617816138326813d0 0.3614128577450535d0
+      0.13267757609917497d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 3)))
+      (COULOMB-WAVE-F-ARRAY 0.0d0 1.0d0 2.0d0 ARR)
+      (DATA ARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07161779967468254d0 0.1278101031568499d0
+	 2.2510703464871114d0 -1.4245543587641651d0 0.0d0
+	 0.0d0 8.269219937869759d-15
+	 1.4757362807662922d-14 2.5991577338697564d-13
+	 1.644835970887306d-13)
+   (MULTIPLE-VALUE-LIST
+    (COULOMB-WAVE-FG 1.0d0 2.0d0 2.5d0 1)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    (LIST 0.035021584603913254d0 0.005752500614202263d0
+	  7.116955601984411d-4 6.471726496134135d0
+	  27.57457472159366d0 170.56037293106908d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM
+	((FARR (VECTOR-DOUBLE 3)) (GARR (VECTOR-DOUBLE 3)))
+      (COULOMB-WAVE-FG-ARRAY 1.5d0 1.0d0 1.0d0 FARR GARR)
+      (APPEND (COERCE (DATA FARR) 'LIST)
+	      (COERCE (DATA GARR) 'LIST)))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.33089080691634065d0 0.18070642887252675d0
+      0.06633878804958748d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 3)))
+      (COULOMB-WAVE-SPHF-ARRAY 0.0d0 1.0d0 2.0d0 ARR)
+      (DATA ARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0013809146441856027d0 2.759621819430441d-17)
+   (MULTIPLE-VALUE-LIST (COULOMB-CL 1.0d0 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.10842251310207264d0 0.05111086283184191d0
+      0.011428736368066591d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((CL (VECTOR-DOUBLE 3)))
+      (COULOMB-CL-ARRAY 0.0d0 1.0d0 CL) (DATA CL)))))
+
diff --git a/special-functions/coupling.lisp b/special-functions/coupling.lisp
index 68b66c657ffa0104cbee2a6716ab4a22c96218b4..baa4b62dd563f4a8ed300a1ccd940d5efd25b326 100644
--- a/special-functions/coupling.lisp
+++ b/special-functions/coupling.lisp
@@ -1,71 +1,75 @@
-;********************************************************
-; file:        coupling.lisp                             
-; description: Coupling coefficients                     
-; date:        Sun Mar 19 2006 - 13:30                   
-; author:      Liam M. Healy                             
-; modified:    Mon Jun 12 2006 - 23:01
-;********************************************************
-;;; $Id: $
+;; Coupling coefficients
+;; Liam Healy, Sun Mar 19 2006 - 13:30
+;; Time-stamp: <2008-02-16 20:02:09EST coupling.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
 #|
+;;; FDL
 The Wigner 3-j, 6-j and 9-j symbols give the coupling coefficients for
 combined angular momentum vectors.  Since the arguments of the standard
 coupling coefficient functions are integer or half-integer, the
 arguments of the following functions are, by convention, integers equal
 to twice the actual spin value.  For information on the 3-j coefficients
 see Abramowitz & Stegun, Section 27.9.  The functions described in this
-section are declared in the header file @file{gsl_sf_coupling.h}.
+section are declared in the header file gsl_sf_coupling.h.
 |#
 
-(defun-gsl coupling-3j (two-ja two-jb two-jc two-ma two-mb two-mc)
+(defmfun coupling-3j (two-ja two-jb two-jc two-ma two-mb two-mc)
   "gsl_sf_coupling_3j_e" 
   ((two-ja :int) (two-jb :int) (two-jc :int)
    (two-ma :int) (two-mb :int) (two-mc :int)
    (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The Wigner 3-j coefficient, 
   \pmatrix{ja & jb & jc\cr
          ma & mb & mc\cr}
-  where the arguments are given in half-integer units, @math{ja} =
-  @var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc.")
+  where the arguments are given in half-integer units,
+  ja = two_ja/2, ma = two_ma/2, etc.")
 
-(defun-gsl coupling-6j (two-ja two-jb two-jc two-jd two-je two-jf)
+(defmfun coupling-6j (two-ja two-jb two-jc two-jd two-je two-jf)
   "gsl_sf_coupling_6j_e" 
   ((two-ja :int) (two-jb :int) (two-jc :int)
    (two-jd :int) (two-je :int) (two-jf :int)
    (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The Wigner 6-j coefficient, 
-  \left\{\matrix{ja & jb & jc\cr
-               jd & je & jf\cr}\right\}
-  where the arguments are given in half-integer units, @math{ja} =
-  @var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc.")
+   ja & jb & jc
+   jd & je & jf
+  where the arguments are given in half-integer units, ja =
+  two_ja/2, ma = two_ma/2, etc.")
 
-(defun-gsl coupling-9j
+(defmfun coupling-9j
     (two-ja two-jb two-jc two-jd two-je two-jf two-jg two-jh two-ji)
   "gsl_sf_coupling_9j_e" 
   ((two-ja :int) (two-jb :int) (two-jc :int)
    (two-jd :int) (two-je :int) (two-jf :int)
    (two-jg :int) (two-jh :int) (two-ji :int)
    (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The Wigner 9-j coefficient, 
-  \left\{\matrix{ja & jb & jc\cr
-               jd & je & jf\cr
-               jg & jh & ji\cr}\right\}
-  where the arguments are given in half-integer units, @math{ja} =
-  @var{two_ja}/2, @math{ma} = @var{two_ma}/2, etc.")
+  ja & jb & jc
+  jd & je & jf
+  jg & jh & ji
+  where the arguments are given in half-integer units,
+  ja = two_ja/2, ma = two_ma/2, etc.")
 
 ;; Check with online calculator http://www-stone.ch.cam.ac.uk/wigner.html
-(lisp-unit:define-test coupling
-  (lisp-unit:assert-first-fp-equal
-   "0.707106781187d+00"
-   (coupling-3j 0 1 1 0 1 -1))
-  (lisp-unit:assert-first-fp-equal
-   "0.408248290464d+00"
-   (coupling-6j 1 1 2 0 2 1))
-  (lisp-unit:assert-first-fp-equal
-   "0.138888888889d+00"
-   (coupling-9j 1 1 2 1 2 1 2 1 1)))
+#|
+(make-tests coupling
+  (coupling-3j 0 1 1 0 1 -1)
+  (coupling-6j 1 1 2 0 2 1)
+  (coupling-9j 1 1 2 1 2 1 2 1 1))
+|#
+
+(LISP-UNIT:DEFINE-TEST COUPLING
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7071067811865475d0 3.14018491736755d-16)
+   (MULTIPLE-VALUE-LIST (COUPLING-3J 0 1 1 0 1 -1)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.408248290463863d0 5.438959822042073d-16)
+   (MULTIPLE-VALUE-LIST (COUPLING-6J 1 1 2 0 2 1)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1388888888888889d0 6.638400825147663d-16)
+   (MULTIPLE-VALUE-LIST (COUPLING-9J 1 1 2 1 2 1 2 1 1))))
diff --git a/special-functions/dawson.lisp b/special-functions/dawson.lisp
index 49348f274e2dcd3d360159ba8db2cafc6448e8c3..0a15c13328ed64d265f9e5cacfe0546e22e8b327 100644
--- a/special-functions/dawson.lisp
+++ b/special-functions/dawson.lisp
@@ -1,25 +1,30 @@
-;********************************************************
-; file:        dawson.lisp                               
-; description: Dawson function                           
-; date:        Sun Mar 19 2006 - 14:31                   
-; author:      Liam M. Healy                             
-; modified:    Mon Jun 12 2006 - 23:17
-;********************************************************
-;;; $Id: $
+;; Dawson function
+;; Liam Healy, Sun Mar 19 2006 - 14:31
+;; Time-stamp: <2008-02-16 20:04:51EST dawson.lisp>
+;; $Id: $
 
 #|
-The Dawson integral is defined by @math{\exp(-x^2) \int_0^x dt
-\exp(t^2)}.  A table of Dawson's integral can be found in Abramowitz &
+;;; FDL
+The Dawson integral is defined by \exp(-x^2) \int_0^x dt
+\exp(t^2).  A table of Dawson's integral can be found in Abramowitz &
 Stegun, Table 7.5.  The Dawson functions are declared in the header file
-@file{gsl_sf_dawson.h}.
+gsl_sf_dawson.h.
 |#
 
 (in-package :gsl)
 
-(defun-gsl dawson (x)
+(defmfun dawson (x)
   "gsl_sf_dawson_e" ((x :double) (ret sf-result))
-  :documentation
-  "Dawson's integral for @var{x}.")
+  :documentation			; FDL
+  "Dawson's integral for x.")
+
+#|
+(make-tests dawson
+	    (dawson 1.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST DAWSON
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5380795069127684d0 1.424354102650492d-15)
+   (MULTIPLE-VALUE-LIST (DAWSON 1.0d0))))
 
-(lisp-unit:define-test dawson
-  (lisp-unit:assert-first-fp-equal "0.538079506913d+00" (dawson 1.0d0)))
diff --git a/special-functions/debye.lisp b/special-functions/debye.lisp
index 159e7cccf31ce68b129228255b14758638ef2006..18354a787448a23f3f06bc7fe4d869d7221b50bb 100644
--- a/special-functions/debye.lisp
+++ b/special-functions/debye.lisp
@@ -1,14 +1,11 @@
-;********************************************************
-; file:        debye.lisp                                
-; description: Deybe functions                           
-; date:        Sun Mar 19 2006 - 14:34                   
-; author:      Liam M. Healy                             
-; modified:    Mon Jun 12 2006 - 23:22
-;********************************************************
-;;; $Id: $
+;; Deybe functions
+;; Liam Healy, Sun Mar 19 2006 - 14:34
+;; Time-stamp: <2008-02-16 20:07:26EST debye.lisp>
+;; $Id: $
 
 #|
-The Debye functions @math{D_n(x)} are defined by the following integral,
+;;; FDL
+The Debye functions D_n(x) are defined by the following integral,
 D_n(x) = n/x^n \int_0^x dt (t^n/(e^t - 1))
 For further information see Abramowitz &
 Stegun, Section 27.1.
@@ -16,31 +13,47 @@ Stegun, Section 27.1.
 
 (in-package :gsl)
 
-(defun-gsl debye-1 (x)
+(defmfun debye-1 (x)
   "gsl_sf_debye_1_e" ((x :double) (ret sf-result))
-  :documentation
-  "The first-order Debye function @math{D_1(x) = (1/x) \int_0^x dt (t/(e^t - 1))}.")
+  :documentation			; FDL
+  "The first-order Debye function D_1(x) = (1/x) \int_0^x dt (t/(e^t - 1)).")
 
-(defun-gsl debye-2 (x)
+(defmfun debye-2 (x)
   "gsl_sf_debye_2_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The second-order Debye function
-   @math{D_2(x) = (2/x^2) \int_0^x dt (t^2/(e^t - 1))}.")
+   D_2(x) = (2/x^2) \int_0^x dt (t^2/(e^t - 1)).")
 
-(defun-gsl debye-3 (x)
+(defmfun debye-3 (x)
   "gsl_sf_debye_3_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The third-order Debye function
-   @math{D_3(x) = (3/x^3) \int_0^x dt (t^3/(e^t - 1))}.")
+   D_3(x) = (3/x^3) \int_0^x dt (t^3/(e^t - 1)).")
 
-(defun-gsl debye-4 (x)
+(defmfun debye-4 (x)
   "gsl_sf_debye_4_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The fourth-order Debye function
-  @math{D_4(x) = (4/x^4) \int_0^x dt (t^4/(e^t - 1))}.")
+  D_4(x) = (4/x^4) \int_0^x dt (t^4/(e^t - 1)).")
 
-(lisp-unit:define-test debye
-  (lisp-unit:assert-first-fp-equal "0.777504634112d+00" (debye-1 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.707878475628d+00" (debye-2 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.674415564078d+00" (debye-3 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.654874068887d+00" (debye-4 1.0d0)))
+#|
+(make-tests debye
+  (debye-1 1.0d0)
+  (debye-2 1.0d0)
+  (debye-3 1.0d0)
+  (debye-4 1.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST DEBYE
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7775046341122482d0 4.117962028082377d-16)
+   (MULTIPLE-VALUE-LIST (DEBYE-1 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7078784756278294d0 4.948781517277596d-16)
+   (MULTIPLE-VALUE-LIST (DEBYE-2 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6744155640778147d0 5.450931753448871d-16)
+   (MULTIPLE-VALUE-LIST (DEBYE-3 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.654874068886737d0 5.769372522202218d-16)
+   (MULTIPLE-VALUE-LIST (DEBYE-4 1.0d0))))
diff --git a/special-functions/dilogarithm.lisp b/special-functions/dilogarithm.lisp
index 6db149da5b9a91c9a6e601ef13040acce6983b11..9d5b7cb133096b397bcb24693f755ee44616568a 100644
--- a/special-functions/dilogarithm.lisp
+++ b/special-functions/dilogarithm.lisp
@@ -1,31 +1,37 @@
-;********************************************************
-; file:        dilogarithm.lisp                          
-; description: Dilogarithm                               
-; date:        Fri Mar 17 2006 - 18:44                   
-; author:      Liam M. Healy
-; modified:    Mon Oct  8 2007 - 11:27
-;********************************************************
+;; Dilogarithm
+;; Liam Healy, Fri Mar 17 2006 - 18:44
+;; Time-stamp: <2008-02-16 20:19:55EST dilogarithm.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
 ;;; dilog merge complex and real
 (defgeneric dilogarithm (x)
-  (:documentation "The dilogarithm."))
+  (:documentation			; FDL
+   "The dilogarithm."))
 
-(defun-gsl dilogarithm ((x float))
+(defmfun dilogarithm ((x float))
   "gsl_sf_dilog_e" ((x :double) (ret sf-result))
   :type :method)
 
-(defun-gsl dilogarithm ((x complex))
+(defmfun dilogarithm ((x complex))
   "gsl_sf_complex_dilog_e"
   (((abs x) :double) ((phase x) :double) (re sf-result) (im sf-result))
   :type :method
   :return ((complex (val re) (val im)) (complex (err re) (err im))))
 
-(lisp-unit:define-test dilogarithm
-  (lisp-unit:assert-first-fp-equal
-   "0.164493406685d+01" (dilogarithm 1.0d0))
-  (lisp-unit:assert-equal
-   '("-0.205616758356d+00" "0.915965594177d+00")
-   (let ((c (dilogarithm #c(0.0d0 1.0d0))))
-     (lisp-unit:fp-sequence (list (realpart c) (imagpart c))))))
+#|
+(make-tests dilogarithm
+	    (dilogarithm 1.0d0)
+	    (dilogarithm #c(0.0d0 1.0d0)))
+|#
+
+(LISP-UNIT:DEFINE-TEST DILOGARITHM
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.6449340668482264d0 7.304974700020789d-16)
+   (MULTIPLE-VALUE-LIST (DILOGARITHM 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(-0.20561675835602822d0 0.915965594177219d0)
+	 #C(2.100180226255977d-15 7.618282373747058d-16))
+   (MULTIPLE-VALUE-LIST (DILOGARITHM #C(0.0d0 1.0d0)))))
+
diff --git a/special-functions/elementary.lisp b/special-functions/elementary.lisp
index 17cdd46b30e13deb704722b732fb0f7ad3b8d8e6..c6bbbc3f9d707d672d979e7b931c8a5146563225 100644
--- a/special-functions/elementary.lisp
+++ b/special-functions/elementary.lisp
@@ -1,34 +1,35 @@
-;********************************************************
-; file:        elementary.lisp                           
-; description: Elementary functions                      
-; date:        Mon Mar 20 2006 - 21:43                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:04
-;********************************************************
-;;; $Id: $
+;; Elementary functions
+;; Liam Healy, Mon Mar 20 2006 - 21:43
+;; Time-stamp: <2008-02-16 20:22:17EST elementary.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl multiply (x y)
+(defmfun 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.")
+  :documentation			; FDL
+  "Multiplies x and y returning the product and associated error.")
 
-(defun-gsl multiply-err (x dx y dy)
+(defmfun 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 
-   @math{xy +/- xy \sqrt((dx/x)^2 +(dy/y)^2)} 
+  :documentation			; FDL
+  "Multiplies x and y with associated absolute
+   errors dx and dy.  The product 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)))
+#|
+(make-tests elementary
+	    (multiply 3.0d0 2.0d0)
+	    (multiply-err 3.0d0 0.1d0 2.0d0 0.1d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST ELEMENTARY
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 6.0d0 2.6645352591003757d-15)
+   (MULTIPLE-VALUE-LIST (MULTIPLY 3.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 6.0d0 0.5000000000000027d0)
+   (MULTIPLE-VALUE-LIST
+    (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 af3a0f057c80a13088e18cbc7f5b9f438ccfd0bf..38f0044ff55fd5bcd7dee8af55d33a473e7f27ec 100644
--- a/special-functions/elliptic-functions.lisp
+++ b/special-functions/elliptic-functions.lisp
@@ -1,22 +1,17 @@
-;********************************************************
-; file:        elliptic-functions.lisp
-; description: Jacobian elliptic functions                        
-; date:        Mon Mar 20 2006 - 22:21                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:16
-;********************************************************
-;;; $Id: $
+;; Jacobian elliptic functions
+;; Liam Healy, Mon Mar 20 2006 - 22:21
+;; Time-stamp: <2008-02-16 20:42:55EST elliptic-functions.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl jacobian-elliptic-functions (u m)
+(defmfun 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 ((val sn) (val cn) (val dn) (err sn) (err cn) (err dn)))
+  :return ((val sn) (val cn) (val dn) (err sn) (err cn) (err dn))
+  :documentation			; FDL
+  "The Jacobian elliptic functions sn(u|m),
+  cn(u|m), dn(u|m) computed by descending Landen transformations.")
 
 ;;; > (jacobian-elliptic-functions 0.61802d0 0.5d0)
 ;;; 0.564575752943391
@@ -30,12 +25,19 @@
 ;;; ;;;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)))
+#|
+(make-tests elliptic-functions
+	    (jacobian-elliptic-functions 0.2d0 0.81d0)
+	    (jacobian-elliptic-functions 0.61802d0 1.5d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST ELLIPTIC-FUNCTIONS
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19762082367187703d0 0.9802785369736752d0
+	 0.9840560289645665d0 0.0d0 0.0d0 0.0d0)
+   (MULTIPLE-VALUE-LIST
+    (JACOBIAN-ELLIPTIC-FUNCTIONS 0.2d0 0.81d0)))
+  (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 6efb5c3b0362a15ca60ce42c53205f779668083f..01b97b731925b5f0e2cf439ac0f5b06b2887cf51 100644
--- a/special-functions/elliptic-integrals.lisp
+++ b/special-functions/elliptic-integrals.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        elliptic-integrals.lisp                   
-; description: Elliptic integrals                        
-; date:        Mon Mar 20 2006 - 21:50                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:09
-;********************************************************
-;;; $Id: $
+;; Elliptic integrals
+;; Liam Healy, Mon Mar 20 2006 - 21:50
+;; Time-stamp: <2008-02-16 20:47:06EST elliptic-integrals.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,93 +9,128 @@
 ;;;; Legendre form of complete elliptic integrals
 ;;;;****************************************************************************
 
-(defun-gsl elliptic-integral-K-complete (k)
+(defmfun elliptic-integral-K-complete (k)
   "gsl_sf_ellint_Kcomp_e" ((k :double) :mode (ret sf-result))
-  :documentation
-  "The complete elliptic integral @math{K(k)}.")
+  :documentation			; FDL
+  "The complete elliptic integral K(k).")
 
-(defun-gsl elliptic-integral-E-complete (k)
+(defmfun elliptic-integral-E-complete (k)
   "gsl_sf_ellint_Ecomp_e" ((k :double) :mode (ret sf-result))
-  :documentation "The complete elliptic integral @math{E(k)}.")
+  :documentation			; FDL
+  "The complete elliptic integral E(k).")
 
 ;;;;****************************************************************************
 ;;;; Legendre form of incomplete elliptic integrals
 ;;;;****************************************************************************
 
-(defun-gsl elliptic-integral-F (phi k)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral F(\phi,k).")
 
-(defun-gsl elliptic-integral-E (phi k)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral E(\phi,k).")
 
-(defun-gsl elliptic-integral-P (phi k n)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral P(\phi,k,n).")
 
-(defun-gsl elliptic-integral-D (phi k n)
+(defmfun 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).")
+  :documentation			; FDL
+  "The incomplete elliptic integral D(phi,k,n) which is
+   defined through the Carlson form RD(x,y,z)
+   by the following relation:
+   D(phi,k,n) = RD (1-sin^2(phi), 1-k^2 sin^2(phi), 1).")
 
 ;;;;****************************************************************************
 ;;;; Carlson forms
 ;;;;****************************************************************************
 
-(defun-gsl elliptic-integral-RC (x y)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral RC(x,y).")
 
-(defun-gsl elliptic-integral-RD (x y z)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral RD(x,y,z).")
 
-(defun-gsl elliptic-integral-RF (x y z)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral RF(x,y,z).")
 
-(defun-gsl elliptic-integral-RJ (x y z p)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The incomplete elliptic integral RJ(x,y,z,p).")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test elliptic-integrals
-  (lisp-unit:assert-first-fp-equal "0.157079632679d+01"
-				   (elliptic-integral-K-complete 0.0d0))
-  (lisp-unit:assert-first-fp-equal "0.157079632679d+01"
-				   (elliptic-integral-E-complete 0.0d0))
-  (lisp-unit:assert-first-fp-equal "-0.677417538204d+00"
-				   (elliptic-integral-F -0.5d0 2.0d0))
-  (lisp-unit:assert-first-fp-equal "-0.401819480553d+00"
-				   (elliptic-integral-E -0.5d0 2.0d0))
-  (lisp-unit:assert-first-fp-equal "-0.617791316339d+00"
-				   (elliptic-integral-P -0.5d0 2.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal "-0.688995144126d-01"
-				   (elliptic-integral-D -0.5d0 2.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.881373587020d+00"
-				   (elliptic-integral-RC 2.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.799259963030d+00"
-				   (elliptic-integral-RD 2.0d0 1.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.881373587020d+00"
-				   (elliptic-integral-RF 2.0d0 1.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.522800417499d+00"
-				   (elliptic-integral-RJ 2.0d0 1.0d0 1.0d0 2.0d0)))
+#|
+(make-tests elliptic-integrals
+  (elliptic-integral-K-complete 0.0d0)
+  (elliptic-integral-E-complete 0.0d0)
+  (elliptic-integral-F -0.5d0 2.0d0)
+  (elliptic-integral-E -0.5d0 2.0d0)
+  (elliptic-integral-P -0.5d0 2.0d0 1.0d0)
+  (elliptic-integral-D -0.5d0 2.0d0 1.0d0)
+  (elliptic-integral-RC 2.0d0 1.0d0)
+  (elliptic-integral-RD 2.0d0 1.0d0 1.0d0)
+  (elliptic-integral-RF 2.0d0 1.0d0 1.0d0)
+  (elliptic-integral-RJ 2.0d0 1.0d0 1.0d0 2.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST ELLIPTIC-INTEGRALS
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.5707963267948966d0 4.598091522633788d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-K-COMPLETE 0.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.5707963267948966d0 3.487868498008632d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-E-COMPLETE 0.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.6774175382039307d0 3.008338192795582d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-F -0.5d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.4018194805534952d0 4.232239429377521d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-E -0.5d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.617791316339182d0 1.6365659051690947d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-P -0.5d0 2.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.06889951441260889d0 3.059753091454848d-17)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-D -0.5d0 2.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8813735870195432d0 1.9570424992111216d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-RC 2.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7992599630303281d0 1.7747136272346433d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-RD 2.0d0 1.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8813735870195432d0 1.9570424992111216d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-RF 2.0d0 1.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5228004174989865d0 1.160850121582039d-16)
+   (MULTIPLE-VALUE-LIST
+    (ELLIPTIC-INTEGRAL-RJ 2.0d0 1.0d0 1.0d0 2.0d0))))
diff --git a/special-functions/error-functions.lisp b/special-functions/error-functions.lisp
index 480b01de5482ea76f64f4f4df77036ef19719fad..c1543bbe8836275cb6219525d34240c4fe606a6a 100644
--- a/special-functions/error-functions.lisp
+++ b/special-functions/error-functions.lisp
@@ -1,52 +1,71 @@
-;********************************************************
-; file:        error-functions.lisp                      
-; description: Error functions                           
-; date:        Mon Mar 20 2006 - 22:31                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:20
-;********************************************************
-;;; $Id: $
+;; Error functions
+;; Liam Healy, Mon Mar 20 2006 - 22:31
+;; Time-stamp: <2008-02-16 20:54:49EST error-functions.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl erf (x)
+(defmfun erf (x)
   "gsl_sf_erf_e" ((x :double) (ret sf-result))
-  :documentation
-  "The error function @math{erf(x)}, where
-  @math{erf(x) = (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2)}.")
+  :documentation			; FDL
+  "The error function erf(x), where
+  erf(x) = (2/\sqrt(\pi)) \int_0^x dt \exp(-t^2).")
 
-(defun-gsl erfc (x)
+(defmfun erfc (x)
   "gsl_sf_erfc_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The complementary error function 
-  @math{erfc(x) = 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2)}.")
+  erfc(x) = 1 - erf(x) = (2/\sqrt(\pi)) \int_x^\infty \exp(-t^2).")
 
-(defun-gsl log-erfc (x)
+(defmfun 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))}.")
+  :documentation			; FDL
+  "The logarithm of the complementary error function \log(\erfc(x)).")
 
-(defun-gsl erf-Z (x)
+(defmfun erf-Z (x)
   "gsl_sf_erf_Z_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The Gaussian probability density function 
-  @math{Z(x) = (1/\sqrt@{2\pi@}) \exp(-x^2/2)}.")
+  Z(x) = (1/sqrt{2\pi}) \exp(-x^2/2)}.")
 
-(defun-gsl erf-Q (x)
+(defmfun erf-Q (x)
   "gsl_sf_erf_Q_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The upper tail of the Gaussian probability function 
-  @math{Q(x) = (1/\sqrt@{2\pi@}) \int_x^\infty dt \exp(-t^2/2)}.")
+  Q(x) = (1/\sqrt{2\pi}) \int_x^\infty dt \exp(-t^2/2)}.")
 
-(defun-gsl hazard (x)
+(defmfun hazard (x)
   "gsl_sf_hazard_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "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)))
+#|
+(make-tests error-functions
+  (erf 1.0d0)
+  (erfc 1.0d0)
+  (log-erfc 1.0d0)
+  (erf-z 1.0d0)
+  (erf-q 1.0d0)
+  (hazard 1.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST ERROR-FUNCTIONS
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8427007929497149d0 7.789237746491556d-16)
+   (MULTIPLE-VALUE-LIST (ERF 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1572992070502851d0 4.0468944536809554d-16)
+   (MULTIPLE-VALUE-LIST (ERFC 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -1.8496055099332485d0 3.394126565390616d-15)
+   (MULTIPLE-VALUE-LIST (LOG-ERFC 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.24197072451914334d0 1.611848817878303d-16)
+   (MULTIPLE-VALUE-LIST (ERF-Z 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.15865525393145707d0 2.832400331480832d-16)
+   (MULTIPLE-VALUE-LIST (ERF-Q 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.5251352761609807d0 5.532094155354489d-15)
+   (MULTIPLE-VALUE-LIST (HAZARD 1.0d0))))
+
diff --git a/special-functions/exponential-functions.lisp b/special-functions/exponential-functions.lisp
index 81ea52ae116a98f07ee565f85a20ba9109628e78..3e30cf46fa3572fca2e4dba1a320a72422d56e67 100644
--- a/special-functions/exponential-functions.lisp
+++ b/special-functions/exponential-functions.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        exponential-functions.lisp                
-; description: Exponential functions                     
-; date:        Tue Mar 21 2006 - 17:05                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:36
-;********************************************************
-;;; $Id: $
+;; Exponential functions
+;; Liam Healy, Tue Mar 21 2006 - 17:05
+;; Time-stamp: <2008-02-16 21:01:05EST exponential-functions.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,117 +9,136 @@
 ;;;; Exponential Functions
 ;;;;****************************************************************************
 
-(defun-gsl gsl-exp (x)
+(defmfun gsl-exp (x)
   "gsl_sf_exp_e" ((x :double) (ret sf-result))
-  :documentation "The exponential function.")
+  :documentation			; FDL
+  "The exponential function.")
 
-(defun-gsl exp-scaled (x)
+(defmfun exp-scaled (x)
   "gsl_sf_exp_e10_e" ((x :double) (ret sf-result-e10))
-  :documentation
+  :documentation			; FDL
   "The exponential function scaled. This function may be useful if the value
-   of @math{\exp(x)} would overflow the  numeric range of @code{double}.")
+   of exp(x) would overflow the numeric range of double.")
 
-(defun-gsl exp-mult (x y)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "Exponentiate x and multiply by the factor y to
+   return the product y \exp(x).")
 
-(defun-gsl exp-mult-scaled (x y)
+(defmfun 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.")
+  :documentation			; FDL
+  "The product y \exp(x) with extended numeric range.")
 
 ;;;;****************************************************************************
 ;;;; Relative Exponential Functions
 ;;;;****************************************************************************
 
-(defun-gsl expm1 (x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "\exp(x)-1 using an algorithm that is accurate for small x.")
 
-(defun-gsl exprel (x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "(\exp(x)-1)/x using an algorithm that is accurate for small x.
+  For small x the algorithm is based on the expansion
+  (\exp(x)-1)/x = 1 + x/2 + x^2/(2*3) + x^3/(2*3*4) + ...")
 
-(defun-gsl exprel-2 (x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "2(\exp(x)-1-x)/x^2 using an algorithm that is accurate for small
+   x.  For small x the algorithm is based on the expansion
+   2(\exp(x)-1-x)/x^2 = 1 + x/3 + x^2/(3*4) + x^3/(3*4*5) + ...")
 
-(defun-gsl exprel-n (n x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "N-relative exponential, which is the n-th generalization
+   of the functions #'exprel and #'exprel-2.")
 
 ;;;;****************************************************************************
 ;;;; Exponentiation With Error Estimate
 ;;;;****************************************************************************
 
-(defun-gsl exp-err (x dx)
+(defmfun 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}.")
+  :documentation			; FDL
+  "Exponentiate x with an associated absolute error dx.")
 
-(defun-gsl exp-err-scaled (x dx)
+(defmfun 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}
+  :documentation			; FDL
+  "Exponentiate x with an associated absolute error dx
   and with extended numeric range.")
 
-(defun-gsl exp-mult-err (x dx y dy)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The product y \exp(x) for the quantities x, y
+   with associated absolute errors dx, dy.")
 
-(defun-gsl exp-mult-err-scaled (x y)
+(defmfun 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
+  :documentation			; FDL
+  "The product y \exp(x) for the quantities x, y
+   with associated absolute errors dx, dy and with
    extended numeric range.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test exponential-functions
-  (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))
-  (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)))
+#|
+(make-tests exponential-functions
+  (gsl-exp 3.0d0)
+  (exp-scaled 555.0d0)
+  (exp-mult 101.0d0 5.0d0)
+  (exp-mult-scaled 555.0d0 101.0d0)
+  (expm1 0.0001d0)
+  (exprel 0.0001d0)
+  (exprel-2 0.001d0)
+  (exprel-n 3 0.001d0)
+  (exp-err 3.0d0 0.001d0)
+  (exp-mult-err 3.0d0 0.001d0 23.0d0 0.001d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST EXPONENTIAL-FUNCTIONS
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 20.085536923187668d0 8.91977022163267d-15)
+   (MULTIPLE-VALUE-LIST (GSL-EXP 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0800340716201098d0 241 2.666751014771678d-13)
+   (MULTIPLE-VALUE-LIST (EXP-SCALED 555.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.6535299896840335d44 8.355840218353793d30)
+   (MULTIPLE-VALUE-LIST (EXP-MULT 101.0d0 5.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0908344123363103d0 243 2.7201204352005766d-15)
+   (MULTIPLE-VALUE-LIST
+    (EXP-MULT-SCALED 555.0d0 101.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0000500016667085d-4 4.441114150507224d-20)
+   (MULTIPLE-VALUE-LIST (EXPM1 1.d-4)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0000500016667084d0 4.4411141505072235d-16)
+   (MULTIPLE-VALUE-LIST (EXPREL 1.d-4)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0003334166833362d0 4.442372766015162d-16)
+   (MULTIPLE-VALUE-LIST (EXPREL-2 0.001d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0002500500083344d0 2.665201526164121d-15)
+   (MULTIPLE-VALUE-LIST (EXPREL-N 3 0.001d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 20.085536923187668d0 0.04017108054156605d0)
+   (MULTIPLE-VALUE-LIST (EXP-ERR 3.0d0 0.001d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 461.9673492333164d0 0.4820528861567092d0)
+   (MULTIPLE-VALUE-LIST
+    (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 ed430e60e604138048d10446690271d637159040..2fa7364c78fd214824f49ad1f1db725e455acfcd 100644
--- a/special-functions/exponential-integrals.lisp
+++ b/special-functions/exponential-integrals.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        exponential-integrals.lisp                
-; description: Exponential integrals                     
-; date:        Tue Mar 21 2006 - 17:37                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:40
-;********************************************************
-;;; $Id: $
+;; Exponential integrals
+;; Liam Healy, Tue Mar 21 2006 - 17:37
+;; Time-stamp: <2008-02-16 22:16:29EST exponential-integrals.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,106 +9,124 @@
 ;;;; Exponential Integral
 ;;;;****************************************************************************
 
-(defun-gsl expint-E1 (x)
+(defmfun expint-E1 (x)
   "gsl_sf_expint_E1_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The exponential integral
-   @math{E_1(x)}, E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t..")
+   E_1(x)}, E_1(x) := \Re \int_1^\infty dt \exp(-xt)/t..")
 
-(defun-gsl expint-E2 (x)
+(defmfun expint-E2 (x)
   "gsl_sf_expint_E2_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The second-order exponential integral
-   @math{E_2(x)}, E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.")
+   E_2(x)}, E_2(x) := \Re \int_1^\infty dt \exp(-xt)/t^2.")
 
 ;;;;****************************************************************************
 ;;;; Ei
 ;;;;****************************************************************************
 
-(defun-gsl expint-Ei (x)
+(defmfun expint-Ei (x)
     "gsl_sf_expint_Ei_e" ((x :double) (ret sf-result))
-  :documentation
-  "The exponential integral @math{Ei(x)},
+  :documentation			; FDL
+  "The exponential integral Ei(x),
    Ei(x) := - PV\left(\int_{-x}^\infty dt \exp(-t)/t\right).")
 
 ;;;;****************************************************************************
 ;;;; Hyperbolic Integrals
 ;;;;****************************************************************************
 
-(defun-gsl Shi (x)
+(defmfun Shi (x)
   "gsl_sf_Shi_e" ((x :double) (ret sf-result))
-  :documentation
-  "The integral @math{Shi(x) = \int_0^x dt \sinh(t)/t}.")
+  :documentation			; FDL
+  "The integral Shi(x) = \int_0^x dt \sinh(t)/t.")
 
-(defun-gsl Chi (x)
+(defmfun Chi (x)
   "gsl_sf_Chi_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "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.")
+   Chi(x) := \Re[ \gamma_E + \log(x) + \int_0^x dt (\cosh[t]-1)/t],
+   where \gamma_E} is the Euler constant.")
 
 ;;;;****************************************************************************
 ;;;; Ei-3
 ;;;;****************************************************************************
 
-(defun-gsl expint-3 (x)
+(defmfun 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 @math{x >= 0}.")
+  :documentation			; FDL
+  "The third-order exponential integral Ei_3(x) = \int_0^xdt \exp(-t^3)
+  for x >= 0.")
 
 ;;;;****************************************************************************
 ;;;; Trigonometric Integrals
 ;;;;****************************************************************************
 
-(defun-gsl Si (x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The Sine integral Si(x) = \int_0^x dt \sin(t)/t.")
 
-(defun-gsl Ci (x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The Cosine integral Ci(x) = -\int_x^\infty dt \cos(t)/t
+   for x > 0.")
 
 ;;;;****************************************************************************
 ;;;; Trigonometric Integrals
 ;;;;****************************************************************************
 
-(defun-gsl atanint (x)
+(defmfun atanint (x)
   "gsl_sf_atanint_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The Arctangent integral, which is defined as
-   @math{AtanInt(x) = \int_0^x dt \arctan(t)/t}.")
+   AtanInt(x) = \int_0^x dt \arctan(t)/t.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test exponential-integrals
-  (lisp-unit:assert-error 'gsl-error (expint-E1 0.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.219383934396d+00"
-   (expint-E1 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.495423435600d+01"
-   (expint-Ei 2.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.136373067344d+01"
-   (Shi 1.25d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.121731730091d+01"
-   (Chi 1.25d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.868892654126d+00"
-   (expint-3 1.25d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.114644641567d+01"
-   (si 1.25d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.434300724034d+00"
-   (ci 1.25d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.110361916168d+01"
-   (atanint 1.25d0)))
+#|
+;;; Macroexpanding in SLIME with slime-macroexpand-1 will produce the
+;;; wrong error type for the first test.  Instead, evaluate in
+;;; listener with (macroexpand-1 '(make-tests ... )).
+
+(make-tests exponential-integrals
+  (expint-E1 0.0d0)
+  (expint-E1 1.0d0)
+  (expint-Ei 2.0d0)
+  (Shi 1.25d0)
+  (Chi 1.25d0)
+  (expint-3 1.25d0)
+  (si 1.25d0)
+  (ci 1.25d0)
+  (atanint 1.25d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST EXPONENTIAL-INTEGRALS
+  (LISP-UNIT:ASSERT-ERROR 'GSL-ERROR (EXPINT-E1 0.0d0))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.21938393439552029d0 2.6541220085226265d-16)
+   (MULTIPLE-VALUE-LIST (EXPINT-E1 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 4.954234356001891d0 7.64289440273947d-15)
+   (MULTIPLE-VALUE-LIST (EXPINT-EI 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.363730673440693d0 4.275641706089105d-15)
+   (MULTIPLE-VALUE-LIST (SHI 1.25d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.217317300914783d0 4.21062110717259d-15)
+   (MULTIPLE-VALUE-LIST (CHI 1.25d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8688926541262327d0 4.083149724141479d-16)
+   (MULTIPLE-VALUE-LIST (EXPINT-3 1.25d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.1464464156732344d0 7.36847290397885d-16)
+   (MULTIPLE-VALUE-LIST (SI 1.25d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.4343007240335524d0 1.2207688418479174d-15)
+   (MULTIPLE-VALUE-LIST (CI 1.25d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.103619161676528d0 6.711588415833395d-16)
+   (MULTIPLE-VALUE-LIST (ATANINT 1.25d0))))
+
diff --git a/special-functions/fermi-dirac.lisp b/special-functions/fermi-dirac.lisp
index 2f793d6f7e809122e2b0dd7aeb84b274e8ce7938..95eb5255d35186a931c97c933c10a93e04d190f2 100644
--- a/special-functions/fermi-dirac.lisp
+++ b/special-functions/fermi-dirac.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        fermi-dirac.lisp                          
-; description: Fermi-Dirac function.                     
-; date:        Sat Apr 22 2006 - 16:12                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 21:45
-;********************************************************
-;;; $Id: $
+;; Fermi-Dirac function.
+;; Liam Healy, Sat Apr 22 2006 - 16:12
+;; Time-stamp: <2008-02-16 21:37:25EST fermi-dirac.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,90 +9,103 @@
 ;;;; Complete Fermi-Dirac Integrals
 ;;;;****************************************************************************
 
-(defun-gsl fermi-dirac-m1 (x)
+(defmfun 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 @math{F_@{-1@}(x) = e^x / (1 + e^x)}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral with an index of -1. 
+  This integral is given by F_{-1}(x) = e^x / (1 + e^x).")
 
-(defun-gsl fermi-dirac-0 (x)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral with an index of 0. 
+   This integral is given by F_0(x) = \ln(1 + e^x).")
 
-(defun-gsl fermi-dirac-1 (x)
+(defmfun 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))}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral with an index of 1,
+   F_1(x) = \int_0^\infty dt (t /(\exp(t-x)+1)).")
 
-(defun-gsl fermi-dirac-2 (x)
+(defmfun 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))}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral with an index of 2,
+   F_2(x) = (1/2) \int_0^\infty dt (t^2 /(\exp(t-x)+1)).")
 
-(defun-gsl fermi-dirac-integral (j x)
+(defmfun 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))}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral with an integer index of j,
+   F_j(x) = (1/\Gamma(j+1)) \int_0^\infty dt (t^j /(\exp(t-x)+1)).")
 
-(defun-gsl fermi-dirac-m1/2 (x)
+(defmfun 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)$}")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral F_{-1/2}(x).")
 
-(defun-gsl fermi-dirac-1/2 (x)
+(defmfun 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)$}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral F_{1/2}(x).")
 
-(defun-gsl fermi-dirac-3/2 (x)
+(defmfun 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)$}.")
+  :documentation			; FDL
+  "The complete Fermi-Dirac integral F_{3/2}(x).")
 
 ;;;;****************************************************************************
 ;;;; Incomplete Fermi-Dirac Integrals
 ;;;;****************************************************************************
 
-(defun-gsl fermi-dirac-inc-0 (x b)
+(defmfun fermi-dirac-inc-0 (x b)
   "gsl_sf_fermi_dirac_inc_0_e" ((x :double) (b :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The incomplete Fermi-Dirac integral with an index
-  of zero, @c{$F_0(x,b) = \ln(1 + e^{b-x}) - (b-x)$}.")
+  of zero, F_0(x,b) = \ln(1 + e^{b-x}) - (b-x).")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test fermi-dirac
-  (lisp-unit:assert-first-fp-equal
-   "0.622459331202d+00"
-   (fermi-dirac-m1 0.5d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.113687100611d+01"
-   (fermi-dirac-0 0.75d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.425894306124d+00"   
-   (fermi-dirac-1 -0.75d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.113016301626d+01"
-   (fermi-dirac-2 0.25d0))
-  (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)))
+#|
+(make-tests fermi-dirac
+  (fermi-dirac-m1 0.5d0)
+  (fermi-dirac-0 0.75d0)
+  (fermi-dirac-1 -0.75d0)
+  (fermi-dirac-2 0.25d0)
+  (fermi-dirac-integral 5 12.35d0)
+  (fermi-dirac-m1/2 2.0d0)
+  (fermi-dirac-1/2 2.0d0)
+  (fermi-dirac-3/2 2.0d0)
+  (fermi-dirac-inc-0 2.0d0 0.5d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST FERMI-DIRAC
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6224593312018546d0 4.040305821324309d-16)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-M1 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.1368710061148999d0 1.6653345369377348d-16)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-0 0.75d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.42589430612383183d0 5.594591187346375d-16)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-1 -0.75d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.1301630162645293d0 6.04246936785309d-16)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-2 0.25d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 6668.827087650077d0 2.1696654306572242d-10)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-INTEGRAL 5 12.35d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.4642945890876293d0 4.053118505064653d-15)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-M1/2 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.8237212774015843d0 2.909989668869937d-15)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-1/2 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 4.165414459868321d0 3.666976569654393d-15)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-3/2 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.7014132779827524d0 8.881784197001252d-16)
+   (MULTIPLE-VALUE-LIST (FERMI-DIRAC-INC-0 2.0d0 0.5d0))))
diff --git a/special-functions/gamma.lisp b/special-functions/gamma.lisp
index b39e4fd2c969361ae5b7606f44b32e09943912d5..18e1fcf895c5fd564755743381b0c6fa8ad83f66 100644
--- a/special-functions/gamma.lisp
+++ b/special-functions/gamma.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        gamma.lisp                                
-; description: Gamma functions                           
-; date:        Thu Apr 27 2006 - 22:06                   
-; author:      Liam M. Healy                             
-; modified:    Tue Jun 13 2006 - 22:16
-;********************************************************
-;;; $Id: $
+;; Gamma functions
+;; Liam Healy, Thu Apr 27 2006 - 22:06
+;; Time-stamp: <2008-02-16 21:58:18EST gamma.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -19,239 +15,285 @@
 
 (defconstant +gamma-xmax+ 171.0d0)
 
-(defun-gsl gamma (x)
+(defmfun gamma (x)
   "gsl_sf_gamma_e" ((x :double) (ret sf-result))
-  :documentation "The Gamma function @math{\Gamma(x)}, subject to x
+  :documentation			; FDL
+  "The Gamma function 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+.")
+   Lanczos method. The maximum value of x such that
+   Gamma(x) is not considered an overflow is given by +gamma-xmax+.")
 
-(defun-gsl log-gamma (x)
+(defmfun 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
+  :documentation			; FDL
+  "The logarithm of the Gamma function,
+   log(Gamma(x)), subject to x not a being negative
+   integer.  For x<0 the real part of log(Gamma(x)) is
+   returned, which is equivalent to log(|Gamma(x)|).  The function
    is computed using the real Lanczos method.")
 
-(defun-gsl log-gamma-sign (x)
+(defmfun 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
+  :return ((val ret) (dcref sign) (err ret))
+  :documentation			; FDL
+  "Compute the sign of the gamma function and the logarithm of
+  its magnitude, subject to 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)}."
-  :return ((val ret) (dcref sign) (err ret)))
+  gamma function can be reconstructed using the relation Gamma(x) =
+  sgn * exp(resultlg)}.")
 
-(defun-gsl gamma* (x)
+(defmfun 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.")
+  :documentation			; FDL
+  "The regulated Gamma Function Gamma^*(x)
+  for x > 0, given by
+  Gamma^*(x) = Gamma(x)/(sqrt{2\pi} x^{(x-1/2)} \exp(-x))
+            = (1 + {1 \over 12x} + ...)
+  for x to infinity.")
 
-(defun-gsl 1/gamma (x)
+(defmfun 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.")
+  :documentation			; FDL
+  "The reciprocal of the gamma function,
+  1/\Gamma(x) using the real Lanczos method.")
 
-(defun-gsl log-gamma-complex (z)
+(defmfun log-gamma-complex (z)
   "gsl_sf_lngamma_complex_e"
   (((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
-  @math{arg = \arg(\Gamma(z))} in @math{(-\pi,\pi]}.  Note that the phase
-  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."
+  :documentation			; FDL
+  "Compute log(Gamma(z)) for complex z=z_r+i z_i
+  and z not a negative integer, using the complex Lanczos
+  method.  The returned parameters are lnr = log|\Gamma(z)| and
+  arg = arg(Gamma(z)) in (-pi,pi].  Note that the phase
+  part (arg) is not well-determined when |z| is very large,
+  due to inevitable roundoff in restricting to (-\pi,\pi].  This
+  will result in a :ELOSS error when it occurs.  The absolute
+  value part (lnr), however, never suffers from loss of precision."
   :return
   ((val lnr) (val arg) (err lnr) (err arg)))
 
-(defun-gsl taylor-coefficient (n x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "Compute the Taylor coefficient x^n / n! for x >= 0, n >= 0.")
 
-(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)}.")
+(defmfun factorial (n)
+  "gsl_sf_fact_e" ((n size) (ret sf-result))
+  :documentation			; FDL
+  "The factorial n!, related to the Gamma function by n! = \Gamma(n+1).")
 
-(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}.")
+(defmfun double-factorial (n)
+  "gsl_sf_doublefact_e" ((n size) (ret sf-result))
+  :documentation			; FDL
+  "The double factorial n!! = n(n-2)(n-4) \dots.")
 
-(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}.")
+(defmfun log-factorial (n)
+  "gsl_sf_lnfact_e" ((n size) (ret sf-result))
+  :documentation			; FDL
+  "The logarithm of the factorial of n, log(n!).
+  The algorithm is faster than computing
+  ln(Gamma(n+1)) via #'log-gamma for n < 170, but defers for larger n.")
 
-(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!!)}.")
+(defmfun log-double-factorial (n)
+  "gsl_sf_lndoublefact_e" ((n size) (ret sf-result))
+  :documentation			; FDL
+  "Compute the logarithm of the double factorial of n, log(n!!).")
 
-(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)!)}")
+(defmfun choose (n m)
+  "gsl_sf_choose_e" ((n size) (m size) (ret sf-result))
+  :documentation			; FDL
+  "The combinatorial factor (n choose m) = n!/(m!(n-m)!).")
 
-(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)!)}.")
+(defmfun log-choose (n m)
+  "gsl_sf_lnchoose_e" ((n size) (m size) (ret sf-result))
+  :documentation			; FDL
+  "The logarithm of (n choose m).  This is
+  equivalent to the sum log(n!) - log(m!) - log((n-m)!).")
 
-(defun-gsl pochammer (a x)
+(defmfun 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
+  :documentation			; FDL
+  "The Pochhammer symbol (a)_x := Gamma(a +
+   x)/Gamma(a), subject to a and a+x not being negative
    integers. The Pochhammer symbol is also known as the Apell symbol and
-   sometimes written as @math{(a,x)}.")
+   sometimes written as (a,x).")
 
-(defun-gsl log-pochammer (a x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The logarithm of the Pochhammer symbol,
+  log((a)_x) = log(Gamma(a + x)/Gamma(a)) for a > 0, a+x > 0.")
 
-(defun-gsl log-pochammer-sign (a x)
+(defmfun 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."
+  :documentation			; FDL
+  "The logarithm of the Pochhammer symbol and its sign.
+  The computed parameters are result =
+  log(|(a)_x|) and sgn = sgn((a)_x) where (a)_x :=
+  Gamma(a + x)/Gamma(a), subject to a, a+x not being negative integers."
   :return ((val ret) (dcref sign) (err ret)))
 
-(defun-gsl relative-pochammer (a x)
+(defmfun 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)}.")
+  :documentation			; FDL
+  "The relative Pochhammer symbol ((a)_x -
+  1)/x where (a)_x := Gamma(a + x)/Gamma(a)}.")
 
-(defun-gsl incomplete-gamma (a x)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The normalized incomplete Gamma Function
+  Q(a,x) = 1/Gamma(a) \int_x^\infty dt t^{a-1} \exp(-t) for a > 0, x >= 0.")
 
-(defun-gsl complementary-incomplete-gamma (a x)
+(defmfun 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).")
+  :documentation			; FDL
+  "The complementary normalized incomplete Gamma Function
+  P(a,x) = 1/\Gamma(a) \int_0^x dt t^{a-1} \exp(-t)}
+  for a > 0, x >= 0.  Note that Abramowitz & Stegun
+  call P(a,x) the incomplete gamma function (section 6.5).")
 
-(defun-gsl nonnormalized-incomplete-gamma (a x)
+(defmfun 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
+  :documentation			; FDL
+  "The incomplete Gamma Function
+   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}.")
+   Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t)
+   for a real and x >= 0.")
 
-(defun-gsl beta (a b)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The Beta Function, B(a,b) = Gamma(a)Gamma(b)/Gamma(a+b)} for a > 0, b > 0.")
 
-(defun-gsl log-beta (a b)
+(defmfun 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}.")
+  :documentation			; FDL
+  "The logarithm of the Beta Function, log(B(a,b)) for a > 0, b > 0.")
 
-(defun-gsl incomplete-beta (a b x)
+(defmfun incomplete-beta (a b x)
   "gsl_sf_beta_inc_e"
   ((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}
-   for @math{a > 0}, @math{b > 0}, and @math{0 <= x <= 1}.")
+  :documentation			; FDL
+  "The normalized incomplete Beta function
+   B_x(a,b)/B(a,b) where
+   B_x(a,b) = \int_0^x t^{a-1} (1-t)^{b-1} dt
+   for a > 0, b > 0, and 0 <= x <= 1.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test gamma
-  (lisp-unit:assert-error 'gsl-error (gamma -1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.120000000000d+03"
-   (gamma 6.0d0))
-  (lisp-unit:assert-error 'gsl-error (log-gamma -100.0d0))
-  (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))
-  (lisp-unit:assert-first-fp-equal
-   "0.479001600000d+09"
-   (factorial 12))
-  (lisp-unit:assert-first-fp-equal
-   "0.460800000000d+05"
-   (double-factorial 12))
-  (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))
-  (lisp-unit:assert-error 'gsl-error (choose 3 8))
-  (lisp-unit:assert-first-fp-equal
-   "0.294222741699d+02"
-   (log-choose 67 12))
-  (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)))
+#|
+(make-tests gamma
+  (gamma -1.0d0)
+  (gamma 6.0d0)
+  (log-gamma -100.0d0)
+  (log-gamma 100.0d0)
+  (log-gamma-sign 100.0d0)
+  (gamma* 24.0d0)
+  (1/gamma 8.0d0)
+  (log-gamma-complex #C(10.0d0 10.0d0))
+  (taylor-coefficient 12 3.0d0)
+  (factorial 12)
+  (double-factorial 12)
+  (log-factorial 199)
+  (log-double-factorial 199)
+  (choose 8 3)
+  (choose 3 8)
+  (log-choose 67 12)
+  (pochammer 3.0d0 2.0d0)
+  (log-pochammer 2.0d0 199.0d0)
+  (log-pochammer-sign 2.0d0 199.0d0)
+  (relative-pochammer 2.0d0 9.0d0)
+  (incomplete-gamma 2.0d0 2.0d0)
+  (complementary-incomplete-gamma 2.0d0 2.0d0)
+  (nonnormalized-incomplete-gamma 2.0d0 2.0d0)
+  (beta 5.50d0 1.0d0)
+  (log-beta 5.5d0 1.0d0)
+  (incomplete-beta 1.0d0 1.50d0 0.50d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST GAMMA
+  (LISP-UNIT:ASSERT-ERROR 'GSL-ERROR (GAMMA -1.0d0))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 120.0d0 2.6645352591003757d-14)
+   (MULTIPLE-VALUE-LIST (GAMMA 6.0d0)))
+  (LISP-UNIT:ASSERT-ERROR 'GSL-ERROR (LOG-GAMMA -100.0d0))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 359.13420536957534d0 2.4544868717695813d-13)
+   (MULTIPLE-VALUE-LIST (LOG-GAMMA 100.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 359.13420536957534d0 1.0d0
+	 2.4544868717695813d-13)
+   (MULTIPLE-VALUE-LIST (LOG-GAMMA-SIGN 100.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0034780558311105d0 4.456337769159149d-16)
+   (MULTIPLE-VALUE-LIST (GAMMA* 24.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.984126984126984d-4 1.3216940769347103d-19)
+   (MULTIPLE-VALUE-LIST (1/GAMMA 8.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 8.236131750448724d0 -1.1840378149363078d0
+	 7.315154482555574d-15 2.1796539969587586d-14)
+   (MULTIPLE-VALUE-LIST
+    (LOG-GAMMA-COMPLEX #C(10.0d0 10.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0011094764610389606d0 2.956239149580215d-18)
+   (MULTIPLE-VALUE-LIST (TAYLOR-COEFFICIENT 12 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 4.790016d8 0.0d0)
+   (MULTIPLE-VALUE-LIST (FACTORIAL 12)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 46080.0d0 0.0d0)
+   (MULTIPLE-VALUE-LIST (DOUBLE-FACTORIAL 12)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 857.9336698258575d0 5.777158772429952d-13)
+   (MULTIPLE-VALUE-LIST (LOG-FACTORIAL 199)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 430.17789358084747d0 1.9103736085528287d-13)
+   (MULTIPLE-VALUE-LIST (LOG-DOUBLE-FACTORIAL 199)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 56.0d0 7.460698725481052d-14)
+   (MULTIPLE-VALUE-LIST (CHOOSE 8 3)))
+  (LISP-UNIT:ASSERT-ERROR 'GSL-ERROR (CHOOSE 3 8))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 29.422274169864693d0 1.9338924605168215d-13)
+   (MULTIPLE-VALUE-LIST (LOG-CHOOSE 67 12)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 12.0d0 6.927791673660977d-14)
+   (MULTIPLE-VALUE-LIST (POCHAMMER 3.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 863.2319871924054d0 9.645972767118375d-13)
+   (MULTIPLE-VALUE-LIST (LOG-POCHAMMER 2.0d0 199.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 863.2319871924054d0 1.0d0 9.645972767118375d-13)
+   (MULTIPLE-VALUE-LIST
+    (LOG-POCHAMMER-SIGN 2.0d0 199.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 403199.88888888917d0 3.5998302729212807d-9)
+   (MULTIPLE-VALUE-LIST (RELATIVE-POCHAMMER 2.0d0 9.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.40600584970983766d0 9.568405127077496d-16)
+   (MULTIPLE-VALUE-LIST (INCOMPLETE-GAMMA 2.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5939941502901611d0 3.510166705531295d-15)
+   (MULTIPLE-VALUE-LIST
+    (COMPLEMENTARY-INCOMPLETE-GAMMA 2.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.40600584970983766d0 1.2272947381959672d-15)
+   (MULTIPLE-VALUE-LIST
+    (NONNORMALIZED-INCOMPLETE-GAMMA 2.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.1818181818181818d0 1.0189782228738177d-15)
+   (MULTIPLE-VALUE-LIST (BETA 5.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -1.7047480922384253d0 3.81791013808375d-15)
+   (MULTIPLE-VALUE-LIST (LOG-BETA 5.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6464466094067263d0 1.0178662453689468d-14)
+   (MULTIPLE-VALUE-LIST
+    (INCOMPLETE-BETA 1.0d0 1.5d0 0.5d0))))
diff --git a/special-functions/gegenbauer.lisp b/special-functions/gegenbauer.lisp
index 8556857fb80e12ff5dd33f317cad36871b720af0..a7bfa1b1f19150e27f20fe15da6e0b553026e176 100644
--- a/special-functions/gegenbauer.lisp
+++ b/special-functions/gegenbauer.lisp
@@ -1,33 +1,33 @@
 ;; Gegenbauer polynomials
 ;; Liam Healy, Fri Apr 28 2006 - 20:40
-;; Time-stamp: <2008-02-03 19:23:36EST gegenbauer.lisp>
+;; Time-stamp: <2008-02-16 22:04:38EST gegenbauer.lisp>
 ;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl gegenbauer-1 (lambda x)
+(defmfun gegenbauer-1 (lambda x)
   "gsl_sf_gegenpoly_1_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Gegenbauer polynomial C^{(\lambda)}_1(x)}.")
 
-(defun-gsl gegenbauer-2 (lambda x)
+(defmfun gegenbauer-2 (lambda x)
   "gsl_sf_gegenpoly_2_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Gegenbauer polynomial C^{(\lambda)}_2(x)}.")
 
-(defun-gsl gegenbauer-3 (lambda x)
+(defmfun gegenbauer-3 (lambda x)
   "gsl_sf_gegenpoly_3_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Gegenbauer polynomial C^{(\lambda)}_3(x)}.")
 
-(defun-gsl gegenbauer (n lambda x)
+(defmfun gegenbauer (n lambda x)
   "gsl_sf_gegenpoly_n_e"
   ((n :int) (lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Gegenbauer polynomial C^{(\lambda)}_n(x)} for a specific value of n,
   lambda, x subject to \lambda > -1/2, n >= 0.")
 
-(defun-gsl gegenbauer-array (lambda x result)
+(defmfun gegenbauer-array (lambda x result)
   "gsl_sf_gegenpoly_array"
   (((1- (dim0 result)) :int)
    (lambda :double) (x :double) ((gsl-array result) :pointer))
@@ -43,22 +43,31 @@
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test gegenbauer
-  (lisp-unit:assert-first-fp-equal
-   "0.600000000000d+01"
-   (gegenbauer-1 1.0d0 3.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.350000000000d+02"
-   (gegenbauer-2 1.0d0 3.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.204000000000d+03"
-   (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
-    (letm ((arr (vector-double 4)))
-      (gegenbauer-array 1.0d0 3.0d0 arr) (data arr)))))
+#|
+(make-tests gegenbauer
+  (gegenbauer-1 1.0d0 3.0d0)
+  (gegenbauer-2 1.0d0 3.0d0)
+  (gegenbauer-3 1.0d0 3.0d0)
+  (gegenbauer 4 1.0d0 3.0d0)
+  (letm ((arr (vector-double 4)))
+      (gegenbauer-array 1.0d0 3.0d0 arr) (data arr)))
+|#
+
+(LISP-UNIT:DEFINE-TEST GEGENBAUER
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 6.0d0 5.329070518200751d-15)
+   (MULTIPLE-VALUE-LIST (GEGENBAUER-1 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 35.0d0 1.5765166949677223d-14)
+   (MULTIPLE-VALUE-LIST (GEGENBAUER-2 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 204.0d0 9.126033262418787d-14)
+   (MULTIPLE-VALUE-LIST (GEGENBAUER-3 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1189.0d0 1.056044141023449d-12)
+   (MULTIPLE-VALUE-LIST (GEGENBAUER 4 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(1.0d0 6.0d0 35.0d0 204.0d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 4)))
+      (GEGENBAUER-ARRAY 1.0d0 3.0d0 ARR) (DATA ARR)))))
diff --git a/special-functions/hypergeometric.lisp b/special-functions/hypergeometric.lisp
index a50fbf3840234a76938520019e241d0831fd53e8..554625bed229f8071e523eb6157309c38c4548c7 100644
--- a/special-functions/hypergeometric.lisp
+++ b/special-functions/hypergeometric.lisp
@@ -1,152 +1,180 @@
-;********************************************************
-; file:        hypergeometric.lisp                       
-; description: Hypergeometric function                   
-; date:        Fri Apr 28 2006 - 23:00                   
-; author:      Liam M. Healy                             
-; modified:    Mon Oct  8 2007 - 11:30
-;********************************************************
-;;; $Id: $
+;; Hypergeometric function
+;; Liam Healy, Fri Apr 28 2006 - 23:00
+;; Time-stamp: <2008-02-16 22:11:13EST hypergeometric.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl hypergeometric-0F1 (c x)
+(defmfun hypergeometric-0F1 (c x)
   "gsl_sf_hyperg_0F1_e" ((c :double) (x :double) (ret sf-result))
-  :documentation "The hypergeometric function @math{0F1(c,x)}.")
+  :documentation			; FDL
+  "The hypergeometric function 0F1(c,x).")
 
 (defgeneric hypergeometric-1F1 (m n x)
-  (:documentation "The confluent hypergeometric function
-   @math{1F1(m,n,x) = M(m,n,x)}."))
+  (:documentation			; FDL
+   "The confluent hypergeometric function 1F1(m,n,x) = M(m,n,x)."))
 
-(defun-gsl hypergeometric-1F1 ((m integer) (n integer) x)
+(defmfun hypergeometric-1F1 ((m integer) (n integer) x)
   "gsl_sf_hyperg_1F1_int_e" ((m :int) (n :int) (x :double) (ret sf-result))
   :type :method
   :export t
-  :documentation "The confluent hypergeometric function
-   @math{1F1(m,n,x) = M(m,n,x)} for integer parameters @var{m}, @var{n}.")
+  :documentation			; FDL
+  "The confluent hypergeometric function 1F1(m,n,x) = M(m,n,x)
+   for integer parameters m, n.")
 
-(defun-gsl hypergeometric-1F1 ((a float) (b float) x)
+(defmfun hypergeometric-1F1 ((a float) (b float) x)
   "gsl_sf_hyperg_1F1_e" ((a :double) (b :double) (x :double) (ret sf-result))
   :type :method
-  :documentation "The confluent hypergeometric function
-  @math{1F1(a,b,x) = M(a,b,x)} for general parameters @var{a}, @var{b}.")
+  :documentation			; FDL
+  "The confluent hypergeometric function
+  1F1(a,b,x) = M(a,b,x) for general parameters a, b.")
 
 (defgeneric hypergeometric-U (m n x)
-  (:documentation "The confluent hypergeometric function
-   @math{U(m,n,x)}."))
+  (:documentation			; FDL
+   "The confluent hypergeometric function U(m,n,x)."))
 
-(defun-gsl hypergeometric-U ((m integer) (n integer) x)
+(defmfun hypergeometric-U ((m integer) (n integer) x)
   "gsl_sf_hyperg_U_int_e" ((m :int) (n :int) (x :double) (ret sf-result))
   :type :method
   :export t
-  :documentation "The confluent hypergeometric function
-   @math{U(m,n,x)} for integer parameters @var{m}, @var{n}.")
+  :documentation			; FDL
+  "The confluent hypergeometric function U(m,n,x) for integer parameters m, n.")
 
-(defun-gsl hypergeometric-U ((a float) (b float) x)
+(defmfun hypergeometric-U ((a float) (b float) x)
   "gsl_sf_hyperg_U_e" ((a :double) (b :double) (x :double) (ret sf-result))
   :type :method
   :documentation "The confluent hypergeometric function @math{U(a,b,x)}.")
 
 (defgeneric hypergeometric-U-e10 (m n x)
-  (:documentation "The confluent hypergeometric function
-   @math{U(m,n,x)} using the @code{gsl_sf_result_e10} type
-   to return a result with extended range."))
+  (:documentation			; FDL
+   "The confluent hypergeometric function
+   U(m,n,x) that returns a result with extended range."))
 
-(defun-gsl hypergeometric-U-e10 ((m integer) (n integer) x)
+(defmfun hypergeometric-U-e10 ((m integer) (n integer) x)
   "gsl_sf_hyperg_U_int_e10_e"
   ((m :int) (n :int) (x :double) (ret sf-result-e10))
   :type :method
   :export t
-  :documentation "The confluent hypergeometric function
-  @math{U(m,n,x)} for integer parameters @var{m}, @var{n} using the
-  @code{gsl_sf_result_e10} type to return a result with extended range.")
+  :documentation			; FDL
+  "The confluent hypergeometric function
+   U(m,n,x) for integer parameters m, n that returns a
+   result with extended range.")
 
-(defun-gsl hypergeometric-U-e10 ((a float) (b float) x)
+(defmfun hypergeometric-U-e10 ((a float) (b float) x)
   "gsl_sf_hyperg_U_e10_e"
   ((a :double) (b :double) (x :double) (ret sf-result-e10))
   :type :method
-  :documentation "The confluent hypergeometric function
-  @math{U(a,b,x)} using the @code{gsl_sf_result_e10} type to return a
-  result with extended range.")
+  :documentation			; FDL
+  "The confluent hypergeometric function
+  U(a,b,x) using that returns a result with extended range.")
 
-(defun-gsl hypergeometric-2F1 (a b c x)
+(defmfun hypergeometric-2F1 (a b c x)
   "gsl_sf_hyperg_2F1_e"
   ((a :double) (b :double) (c :double) (x :double) (ret sf-result))
-  :documentation "The Gauss hypergeometric function
-  @math{2F1(a,b,c,x)} for @math{|x| < 1}. If the arguments
-  @math{(a,b,c,x)} are too close to a singularity then the function can
+  :documentation			; FDL
+  "The Gauss hypergeometric function
+  2F1(a,b,c,x) for |x| < 1. If the arguments
+  (a,b,c,x) are too close to a singularity then the function can
   return the error code :EMAXITER when the series
   approximation converges too slowly.  This occurs in the region of
-  @math{x=1}, @math{c - a - b = m} for integer m.")
+  x=1, c - a - b = m for integer m.")
 
-(defun-gsl hypergeometric-2F1-conj (a c x)
+(defmfun hypergeometric-2F1-conj (a c x)
   "gsl_sf_hyperg_2F1_conj_e"
   (((realpart a) :double) ((imagpart a) :double)
    (c :double) (x :double) (ret sf-result))
-  :documentation "The Gauss hypergeometric function
-  @math{2F1(a, a*, c, x)} with complex parameters 
-  for @math{|x| < 1}.")
+  :documentation			; FDL
+  "The Gauss hypergeometric function 2F1(a, a*, c, x) with complex parameters 
+  for |x| < 1.")
 
-(defun-gsl hypergeometric-2F1-renorm (a b c x)
+(defmfun hypergeometric-2F1-renorm (a b c x)
   "gsl_sf_hyperg_2F1_renorm_e"
   ((a :double) (b :double) (c :double) (x :double) (ret sf-result))
-  :documentation "The renormalized Gauss hypergeometric function
-  @math{2F1(a,b,c,x) / \Gamma(c)} for @math{|x| < 1}.")
+  :documentation			; FDL
+  "The renormalized Gauss hypergeometric function
+  2F1(a,b,c,x) / Gamma(c) for |x| < 1.")
 
-(defun-gsl hypergeometric-2F1-conj-renorm (a c x)
+(defmfun hypergeometric-2F1-conj-renorm (a c x)
   "gsl_sf_hyperg_2F1_conj_renorm_e"
   (((realpart a) :double) ((imagpart a) :double)
    (c :double) (x :double) (ret sf-result))
-  :documentation "The renormalized Gauss hypergeometric function
-  @math{2F1(a, a*, c, x) / \Gamma(c)} for @math{|x| < 1}.")
+  :documentation			; FDL
+  "The renormalized Gauss hypergeometric function
+  2F1(a, a*, c, x) / Gamma(c) for |x| < 1.")
 
-(defun-gsl hypergeometric-2F0 (a b x)
+(defmfun hypergeometric-2F0 (a b x)
   "gsl_sf_hyperg_2F0_e"
   ((a :double) (b :double) (x :double) (ret sf-result))
-  :documentation "The hypergeometric function 
-  @math{2F0(a,b,x)}.  The series representation
-  is a divergent hypergeometric series.  However, for @math{x < 0} we
-  have @math{2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)}")
+  :documentation 			; FDL
+  "The hypergeometric function 
+  2F0(a,b,x).  The series representation
+  is a divergent hypergeometric series.  However, for x < 0 we
+  have 2F0(a,b,x) = (-1/x)^a U(a,1+a-b,-1/x)")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test hypergeometric
-  (lisp-unit:assert-first-fp-equal
-   "0.376219569108d+01"
-   (hypergeometric-0f1 0.5d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.543656365692d+01"
-   (hypergeometric-1F1 2 1 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.543656365692d+01"
-   (hypergeometric-1F1 2.0d0 1.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.192694724646d+00"
-   (hypergeometric-U 2.0d0 1.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.192694724646d+00"
-   (hypergeometric-U 2 1 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.192694724646d+00"
-   (hypergeometric-U-e10 2.0d0 1.0d0 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.192694724646d+00"
-   (hypergeometric-U-e10 2 1 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.229739670999d+01"
-   (hypergeometric-2F1 1.0d0 1.2d0 1.0d0 0.5d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.662959493455d+01"
-   (hypergeometric-2F1-conj #c(1.0d0 0.5d0) 0.5d0 0.6d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.374034840521d+01"
-   (hypergeometric-2F1-conj-renorm
-    #C(1.0d0 0.5d0) 0.5d0 0.6d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.229739670999d+01"
-   (hypergeometric-2F1-renorm
-    1.0d0 1.2d0 1.0d0 0.5d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.435139241256d-01"
-   (hypergeometric-2F0 1.0d0 2.0d0 -20.0d0)))
+#|
+(make-tests hypergeometric
+  (hypergeometric-0f1 0.5d0 1.0d0)
+  (hypergeometric-1F1 2 1 1.0d0)
+  (hypergeometric-1F1 2.0d0 1.0d0 1.0d0)
+  (hypergeometric-U 2.0d0 1.0d0 1.0d0)
+  (hypergeometric-U 2 1 1.0d0)
+  (hypergeometric-U-e10 2.0d0 1.0d0 1.0d0)
+  (hypergeometric-U-e10 2 1 1.0d0)
+  (hypergeometric-2F1 1.0d0 1.2d0 1.0d0 0.5d0)
+  (hypergeometric-2F1-conj #C(1.0d0 0.5d0) 0.5d0 0.6d0)
+  (hypergeometric-2F1-conj-renorm #C(1.0d0 0.5d0) 0.5d0 0.6d0)
+  (hypergeometric-2F1-renorm 1.0d0 1.2d0 1.0d0 0.5d0)
+  (hypergeometric-2F0 1.0d0 2.0d0 -20.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST HYPERGEOMETRIC
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.7621956910836287d0 9.207469176703522d-14)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-0F1 0.5d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 5.43656365691809d0 6.0357981467508045d-15)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-1F1 2 1 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 5.43656365691809d0 6.0357981467508045d-15)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-1F1 2.0d0 1.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19269472464638984d0 2.5468520714552053d-12)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-U 2.0d0 1.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19269472464638984d0 2.5468520714552053d-12)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-U 2 1 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19269472464638984d0 0 2.5468520714552053d-12)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-U-E10 2.0d0 1.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.19269472464638984d0 0 2.5468520714552053d-12)
+   (MULTIPLE-VALUE-LIST (HYPERGEOMETRIC-U-E10 2 1 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.29739670999407d0 2.0404981793068d-15)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-2F1 1.0d0 1.2d0 1.0d0 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 6.629594934547447d0 5.761143589129193d-14)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-2F1-CONJ #C(1.0d0 0.5d0) 0.5d0
+			     0.6d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 3.7403484052126372d0 5.884558632199498d-14)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-2F1-CONJ-RENORM #C(1.0d0 0.5d0) 0.5d0
+				    0.6d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.29739670999407d0 3.0607472689602005d-15)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-2F1-RENORM 1.0d0 1.2d0 1.0d0 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.043513924125598485d0 3.81571654717325d-15)
+   (MULTIPLE-VALUE-LIST
+    (HYPERGEOMETRIC-2F0 1.0d0 2.0d0 -20.0d0))))
diff --git a/special-functions/laguerre.lisp b/special-functions/laguerre.lisp
index 9557cbcb3f2606b79ad22bd80b42400f39b7c5c9..9248eea9c5e104780c8147efd3a898529683d9c9 100644
--- a/special-functions/laguerre.lisp
+++ b/special-functions/laguerre.lisp
@@ -1,48 +1,58 @@
-;********************************************************
-; file:        laguerre.lisp                           
-; description: Laguerre polynomials                    
-; date:        Fri Apr 28 2006 - 20:40                   
-; author:      Liam M. Healy                             
-; modified:    Fri Jun 16 2006 - 21:12
-;********************************************************
-;;; $Id: $
+;; Laguerre polynomials
+;; Liam Healy, Fri Apr 28 2006 - 20:40
+;; Time-stamp: <2008-02-16 22:16:49EST laguerre.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl laguerre-1 (a x)
+(defmfun laguerre-1 (a x)
   "gsl_sf_laguerre_1_e" ((a :double) (x :double) (ret sf-result))
-  :documentation "The generalized Laguerre polynomial
-   @math{L^a_1(x)} using explicit representations.")
+  :documentation			; FDL
+  "The generalized Laguerre polynomial L^a_1(x) using
+   explicit representations.")
 
-(defun-gsl laguerre-2 (a x)
+(defmfun laguerre-2 (a x)
   "gsl_sf_laguerre_2_e" ((a :double) (x :double) (ret sf-result))
-  :documentation  "The generalized Laguerre polynomial
-   @math{L^a_2(x)} using explicit representations.")
+  :documentation			; FDL
+  "The generalized Laguerre polynomial
+   L^a_2(x) using explicit representations.")
 
-(defun-gsl laguerre-3 (a x)
+(defmfun laguerre-3 (a x)
   "gsl_sf_laguerre_3_e" ((a :double) (x :double) (ret sf-result))
-  :documentation "The generalized Laguerre polynomial
-   @math{L^a_3(x)} using explicit representations.")
+  :documentation			; FDL
+  "The generalized Laguerre polynomial
+   L^a_3(x) using explicit representations.")
 
-(defun-gsl laguerre (n a x)
+(defmfun laguerre (n a x)
   "gsl_sf_laguerre_n_e" ((n :int) (a :double) (x :double) (ret sf-result))
-  :documentatiOn "The generalized Laguerre polynomials
-  @math{L^a_n(x)} for @math{a > -1}, @math{n >= 0}.")
+  :documentation			; FDL
+  "The generalized Laguerre polynomials L^a_n(x) for a > -1, n >= 0.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test laguerre
-  (lisp-unit:assert-first-fp-equal
-   "-0.100000000000d+01"
-   (laguerre-1 1.0d0 3.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.150000000000d+01"
-   (laguerre-2 1.0d0 3.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.500000000000d+00"
-   (laguerre-3 1.0d0 3.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.875000000000d+00"
-   (laguerre 4 1.0d0 3.0d0)))
+#|
+;;; Don't slime-macroexpand-1, the last one will produce an error that
+;;; shouldn't be there.
+
+(make-tests laguerre
+  (laguerre-1 1.0d0 3.0d0)
+  (laguerre-2 1.0d0 3.0d0)
+  (laguerre-3 1.0d0 3.0d0)
+  (laguerre 4 1.0d0 3.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST LAGUERRE
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -1.0d0 2.220446049250313d-15)
+   (MULTIPLE-VALUE-LIST (LAGUERRE-1 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -1.5d0 1.7985612998927536d-14)
+   (MULTIPLE-VALUE-LIST (LAGUERRE-2 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.5d0 6.59472476627343d-14)
+   (MULTIPLE-VALUE-LIST (LAGUERRE-3 1.0d0 3.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.875d0 6.793349802129357d-14)
+   (MULTIPLE-VALUE-LIST (LAGUERRE 4 1.0d0 3.0d0))))
diff --git a/special-functions/lambert.lisp b/special-functions/lambert.lisp
index b5bbaed7e130db5e896cc4f0de9c902bf5372e71..f495f03ba7f819252d8a1efd4bb44dd22a781a1b 100644
--- a/special-functions/lambert.lisp
+++ b/special-functions/lambert.lisp
@@ -1,42 +1,45 @@
-;********************************************************
-; file:        lambert.lisp
-; description: Lambert's W functions
-; date:        Fri Apr 28 2006 - 20:40                   
-; author:      Liam M. Healy                             
-; modified:    Fri Jun 16 2006 - 21:13
-;********************************************************
-;;; $Id: $
+;; Lambert's W functions
+;; Liam Healy, Fri Apr 28 2006 - 20:40
+;; Time-stamp: <2008-02-16 22:19:40EST lambert.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-;;; Lambert's W functions, @math{W(x)}, are defined to be solutions
-;;; of the equation @math{W(x) \exp(W(x)) = x}. This function has
-;;; multiple branches for @math{x < 0}; however, it has only
-;;; two real-valued branches. We define @math{W_0(x)} to be the
-;;; principal branch, where @math{W > -1} for @math{x < 0}, and 
-;;; @c{$W_{-1}(x)$}
-;;; @math{W_@{-1@}(x)} to be the other real branch, where
-;;; @math{W < -1} for @math{x < 0}.  
+;;; FDL
+;;; Lambert's W functions, W(x), are defined to be solutions
+;;; of the equation W(x) exp(W(x)) = x. This function has
+;;; multiple branches for x < 0; however, it has only
+;;; two real-valued branches. We define W_0(x) to be the
+;;; principal branch, where W > -1 for x < 0, and 
+;;; W_{-1}(x)
+;;; W_{-1}(x) to be the other real branch, where
+;;; W < -1 for x < 0.  
 
-(defun-gsl lambert-W0 (x)
+(defmfun lambert-W0 (x)
   "gsl_sf_lambert_W0_e" ((x :double) (ret sf-result))
-  :documentation
-  "The principal branch of the Lambert W function, @math{W_0(x)}.")
+  :documentation			; FDL
+  "The principal branch of the Lambert W function, W_0(x).")
 
-(defun-gsl lambert-Wm1 (x)
+(defmfun lambert-Wm1 (x)
   "gsl_sf_lambert_Wm1_e" ((x :double) (ret sf-result))
-  :documentation
+  :documentation			; FDL
   "The secondary real-valued branch of the Lambert W function, 
-   @math{W_@{-1@}(x)}.")
+   W_{-1}(x).")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test lambert
-  (lisp-unit:assert-first-fp-equal
-   "0.567143290410d+00"
-   (lambert-W0 1.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.567143290410d+00"
-   (lambert-Wm1 1.0d0)))
+#|
+(make-tests lambert
+  (lambert-W0 1.0d0)
+  (lambert-Wm1 1.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST LAMBERT
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5671432904097838d0 2.518622157098455d-16)
+   (MULTIPLE-VALUE-LIST (LAMBERT-W0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5671432904097838d0 2.518622157098455d-16)
+   (MULTIPLE-VALUE-LIST (LAMBERT-WM1 1.0d0))))
diff --git a/special-functions/legendre.lisp b/special-functions/legendre.lisp
index 54d065fdc787e4d819fe04377b66e5feef371295..f560da3b562a4e019192f6287988600a35a8ad90 100644
--- a/special-functions/legendre.lisp
+++ b/special-functions/legendre.lisp
@@ -1,6 +1,6 @@
 ;; Legendre functions
 ;; Liam Healy, Sat Apr 29 2006 - 19:16
-;; Time-stamp: <2008-02-03 16:58:53EST legendre.lisp>
+;; Time-stamp: <2008-02-16 22:30:23EST legendre.lisp>
 ;; $Id: $
 
 (in-package :gsl)
@@ -11,31 +11,31 @@
 ;;;; Legendre polynomials
 ;;;;****************************************************************************
 
-(defun-gsl legendre-P1 (x)
+(defmfun legendre-P1 (x)
   "gsl_sf_legendre_P1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre polynomials P_1(x) using an explicit
    representation.")
 
-(defun-gsl legendre-P2 (x)
+(defmfun legendre-P2 (x)
   "gsl_sf_legendre_P2_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre polynomials P_2(x) using an explicit
    representation.")
 
-(defun-gsl legendre-P3 (x)
+(defmfun legendre-P3 (x)
   "gsl_sf_legendre_P3_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre polynomials P_3(x) using an explicit
    representation.")
 
-(defun-gsl legendre-Pl (l x)
+(defmfun legendre-Pl (l x)
   "gsl_sf_legendre_Pl_e" ((l :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre polynomial P_l(x) for a specific value of l,
    x subject to l >= 0, |x| <= 1.")
 
-(defun-gsl legendre-Pl-array (x array)
+(defmfun legendre-Pl-array (x array)
   "gsl_sf_legendre_Pl_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :documentation			; FDL
@@ -43,7 +43,7 @@
   P_l(x) for l = 0, ..., length(array), |x| <= 1."
   :invalidate (array))
 
-(defun-gsl legendre-Pl-deriv-array (x array)
+(defmfun legendre-Pl-deriv-array (x array)
   "gsl_sf_legendre_Pl_deriv_array"
   (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :documentation			; FDL
@@ -51,19 +51,19 @@
   dP_l(x)/dx, for l = 0, ...,  length(array), |x| <= 1."
   :invalidate (array))
 
-(defun-gsl legendre-Q0 (x)
+(defmfun legendre-Q0 (x)
   "gsl_sf_legendre_Q0_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre function Q_0(x) for x > -1,
    x /= 1.")
 
-(defun-gsl legendre-Q1 (x)
+(defmfun legendre-Q1 (x)
   "gsl_sf_legendre_Q1_e" ((x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre function Q_1(x) for x > -1,
    x /= 1.")
 
-(defun-gsl legendre-Ql (l x)
+(defmfun legendre-Ql (l x)
   "gsl_sf_legendre_Ql_e" ((l :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Legendre function Q_l(x) for x > -1, x /= 1, l >= 0.")
@@ -85,13 +85,13 @@
 ;;; these functions.  Instead use legendre-sphPlm below,
 ;;; which uses a similar recursion, but with the normalized functions.
 
-(defun-gsl legendre-Plm (l m x)
+(defmfun legendre-Plm (l m x)
   "gsl_sf_legendre_Plm_e" ((l :int) (m :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The associated Legendre polynomial
    P_l^m(x) for m >= 0, l >= m, |x| <= 1.")
 
-(defun-gsl legendre-Plm-array (m x array)
+(defmfun legendre-Plm-array (m x array)
   "gsl_sf_legendre_Plm_array"
   (((+ (dim0 array) m -1) :int) (m :int) (x :double)
    ((gsl-array array) :pointer))
@@ -101,7 +101,7 @@
     l = |m|, ..., |m|+length(array)-1} and |x| <= 1."
   :invalidate (array))
 
-(defun-gsl legendre-Plm-deriv-array (m x values derivatives)
+(defmfun legendre-Plm-deriv-array (m x values derivatives)
   "gsl_sf_legendre_Plm_deriv_array"
   (((+ (dim0 values) m -1) :int) (m :int) (x :double)
    ((gsl-array values) :pointer) ((gsl-array derivatives) :pointer))
@@ -111,7 +111,7 @@
     l = |m|, ..., length(values) and |x| <= 1."
   :invalidate (values derivatives))
 
-(defun-gsl legendre-sphPlm (l m x)
+(defmfun legendre-sphPlm (l m x)
   "gsl_sf_legendre_sphPlm_e" ((l :int) (m :int) (x :double) (ret sf-result))
   :documentation			; FDL
   "The normalized associated Legendre polynomial
@@ -120,7 +120,7 @@
    m >= 0, l >= m, |x| <= 1.  These routines avoid the overflows
    that occur for the standard normalization of P_l^m(x).")
 
-(defun-gsl legendre-sphPlm-array (m x array)
+(defmfun legendre-sphPlm-array (m x array)
   "gsl_sf_legendre_sphPlm_array"
   (((+ (dim0 array) m -1) :int) (m :int) (x :double)
    ((gsl-array array) :pointer))
@@ -130,7 +130,7 @@
    for m >= 0, l = |m|, ..., length(array)}, |x| <= 1.0."
   :invalidate (array))
 
-(defun-gsl legendre-sphPlm-deriv-array (m x values derivatives)
+(defmfun legendre-sphPlm-deriv-array (m x values derivatives)
   "gsl_sf_legendre_sphPlm_deriv_array"
   (((+ (dim0 values) m -1) :int) (m :int) (x :double)
    ((gsl-array values) :pointer) ((gsl-array derivatives) :pointer))
@@ -140,7 +140,7 @@
    l = |m|, ..., length(array)}, |x| <= 1.0."
   :invalidate (values derivatives))
 
-(defun-gsl legendre-array-size (lmax m)
+(defmfun legendre-array-size (lmax m)
   "gsl_sf_legendre_array_size" ((lmax :int) (m :int))
   :documentation			; FDL
   "The size of result array needed for the array
@@ -153,40 +153,40 @@
 
 ;;; FDL
 ;;; The Conical Functions P^\mu_{-(1/2)+i\lambda}(x)} and 
-;;; Q^\mu_@{-(1/2)+i\lambda}
+;;; Q^\mu_{-(1/2)+i\lambda}
 ;;; are described in Abramowitz & Stegun, Section 8.12.
 
-(defun-gsl legendre-conicalP-half (lambda x)
+(defmfun legendre-conicalP-half (lambda x)
   "gsl_sf_conicalP_half_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The irregular Spherical Conical Function
    P^{1/2}_{-1/2 + i \lambda}(x) for x > -1.")
 
-(defun-gsl legendre-conicalP-mhalf (lambda x)
+(defmfun legendre-conicalP-mhalf (lambda x)
   "gsl_sf_conicalP_mhalf_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The regular Spherical Conical Function
    P^{-1/2}_{-1/2 + i \lambda}(x) for x > -1.")
 
-(defun-gsl legendre-conicalP-0 (lambda x)
+(defmfun legendre-conicalP-0 (lambda x)
   "gsl_sf_conicalP_0_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The conical function P^0_{-1/2 + i \lambda(x)} for x > -1.")
 
-(defun-gsl legendre-conicalP-1 (lambda x)
+(defmfun legendre-conicalP-1 (lambda x)
   "gsl_sf_conicalP_1_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The conical function 
    P^1_{-1/2 + i \lambda}(x)} for x > -1.")
 
-(defun-gsl legendre-regular-spherical-conical (l lambda x)
+(defmfun legendre-regular-spherical-conical (l lambda x)
   "gsl_sf_conicalP_sph_reg_e"
   ((l :int) (lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
   "The Regular Spherical Conical Function
    P^{-1/2-l}_{-1/2 + i \lambda}(x) for x > -1, l >= -1.")
 
-(defun-gsl legendre-regular-cylindrical-conical (l lambda x)
+(defmfun legendre-regular-cylindrical-conical (l lambda x)
   "gsl_sf_conicalP_cyl_reg_e"
   ((l :int) (lambda :double) (x :double) (ret sf-result))
   :documentation			; FDL
@@ -204,7 +204,7 @@
 ;;; is the flat limit, \lambda \to \infty, \eta \to 0, \lambda\eta
 ;;; fixed.
 
-(defun-gsl legendre-H3d-0 (lambda eta)
+(defmfun legendre-H3d-0 (lambda eta)
   "gsl_sf_legendre_H3d_0_e"
   ((lambda :double) (eta :double) (ret sf-result))
   :documentation			; FDL
@@ -212,9 +212,9 @@
    3-dimensional hyperbolic space,
    L^{H3d}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta))
    for \eta >= 0. In the flat limit this takes the form
-   L^{H3d@}_0(\lambda,\eta) = j_0(\lambda\eta).")
+   L^{H3d}_0(\lambda,\eta) = j_0(\lambda\eta).")
 
-(defun-gsl legendre-H3d-1 (lambda eta)
+(defmfun legendre-H3d-1 (lambda eta)
   "gsl_sf_legendre_H3d_1_e"
   ((lambda :double) (eta :double) (ret sf-result))
   :documentation			; FDL
@@ -225,7 +225,7 @@
    for \eta >= 0.  In the flat limit this takes the form 
    L^{H3d}_1(\lambda,\eta) = j_1(\lambda\eta)}.")
 
-(defun-gsl legendre-H3d (l lambda eta)
+(defmfun legendre-H3d (l lambda eta)
   "gsl_sf_legendre_H3d_e"
   ((l :int) (lambda :double) (eta :double) (ret sf-result))
   :documentation			; FDL
@@ -234,7 +234,7 @@
    \eta >= 0, l >= 0.  In the flat limit this takes the form
    L^{H3d}_l(\lambda,\eta) = j_l(\lambda\eta).")
 
-(defun-gsl legendre-H3d-array (lambda eta array)
+(defmfun legendre-H3d-array (lambda eta array)
   "gsl_sf_legendre_H3d_array"
   (((1- (dim0 array)) :int) (lambda :double) (eta :double)
    ((gsl-array array) :pointer))
@@ -251,105 +251,157 @@
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test legendre
-  (lisp-unit:assert-first-fp-equal
-   "0.300000000000d+00"
-   (legendre-P1 0.3d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.365000000000d+00"
-   (legendre-P2 0.3d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.382500000000d+00"
-   (legendre-P3 0.3d0))
-  (lisp-unit:assert-error 'gsl-error (legendre-Pl -4 0.3d0))
-  (lisp-unit:assert-error 'gsl-error (legendre-Pl 4 3.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.729375000000d-01"
-   (legendre-Pl 4 0.3d0))
-  (lisp-unit:assert-equal
-   '("0.100000000000d+01" "0.500000000000d+00"
-     "-0.125000000000d+00" "-0.437500000000d+00")
-   (lisp-unit:fp-sequence
-    (letm ((arr (vector-double 4)))
+#|
+(make-tests legendre
+  (legendre-P1 0.3d0)
+  (legendre-P2 0.3d0)
+  (legendre-P3 0.3d0)
+  (legendre-Pl -4 0.3d0)
+  (legendre-Pl 4 3.0d0)
+  (legendre-Pl 4 0.3d0)
+  (letm ((arr (vector-double 4)))
       (legendre-Pl-array 0.5d0 arr)
-      (data arr))))
-  (lisp-unit:assert-first-fp-equal
-   "0.312852949882d+00"
-   (legendre-Q0 3.3d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.324147346113d-01"
-   (legendre-Q1 3.3d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.402646138474d-02"
-   (legendre-Ql 2 3.3d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.340997502740d+02"
-   (legendre-Plm 4 3 0.5d0))
-  (lisp-unit:assert-equal
-   '("0.225000000000d+01" "0.562500000000d+01"
-     "0.421875000000d+01" "-0.492187500000d+01")
-   (lisp-unit:fp-sequence
-    (letm ((arr (vector-double 4)))
+      (data arr))
+  (legendre-Q0 3.3d0)
+  (legendre-Q1 3.3d0)
+  (legendre-Ql 2 3.3d0)
+  (legendre-Plm 4 3 0.5d0)
+  (letm ((arr (vector-double 4)))
       (legendre-Plm-array 2 0.5d0 arr)
-      (data arr))))
-  (lisp-unit:assert-equal
-   '("-0.300000000000d+01" "0.375000000000d+01"
-     "0.337500000000d+02" "0.557812500000d+02")
-   (lisp-unit:fp-sequence
-    (letm ((val (vector-double 4))
+      (data arr))
+  (letm ((val (vector-double 4))
 	   (deriv (vector-double 4)))
       (legendre-Plm-deriv-array 2 0.5d0 val deriv)
-      (data deriv))))
-  (lisp-unit:assert-first-fp-equal
-   "0.398506257222d-13"
-   (legendre-sphplm 1200 1100 0.5d0))
-  (lisp-unit:assert-equal
-   '("0.248924639500d+00" "0.412794815148d+00"
-     "0.351206555622d+00" "0.515993518936d-01")
-   (lisp-unit:fp-sequence
-    (letm ((arr (vector-double 4)))
+      (data deriv))
+  (legendre-sphplm 1200 1100 0.3d0)
+  (letm ((arr (vector-double 4)))
       (legendre-sphPlm-array 4 0.5d0 arr)
-      (data arr))))
-  ;; suspicious? same answer as legendre-sphPlm-array?
-  (lisp-unit:assert-equal
-   '("-0.663799038667d+00" "-0.275196543432d+00"
-     "0.127103324892d+01" "0.264876673054d+01")
-   (lisp-unit:fp-sequence
-    (letm ((val (vector-double 4))
+      (data arr))
+  (letm ((val (vector-double 4))
 	   (deriv (vector-double 4)))
 	(legendre-sphPlm-deriv-array 4 0.5d0 val deriv)
-	(data deriv))))
-  (lisp-unit:assert-first-fp-equal
-   "-0.125529904888d+00"
-   (legendre-conicalp-half 3.5d0 10.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.627433629279d-01"
-   (legendre-conicalp-mhalf 3.5d0 10.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "-0.131636618937d+00"
-   (legendre-conicalp-0 3.5d0 10.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.174071955601d+00"
-   (legendre-conicalp-1 3.5d0 10.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.898079795297d-03"
-   (legendre-regular-spherical-conical 3 3.5d0 10.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.230602506199d-02"
-   (legendre-regular-cylindrical-conical 3 3.5d0 10.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.920034269259d+00"
-   (legendre-h3d-0 1.0d0 0.5d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.216940264504d+00"
-   (legendre-h3d-1 1.0d0 0.5d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.240061623900d-02"
-   (legendre-h3d 4 1.0d0 0.5d0))
-  (lisp-unit:assert-equal
-   '("0.920034269259d+00" "0.216940264504d+00"
-     "0.479506604883d-01" "0.106637690961d-01")
-   (lisp-unit:fp-sequence
-    (letm ((arr (vector-double 4)))
+	(data deriv))
+  (legendre-conicalp-half 3.5d0 10.0d0)
+  (legendre-conicalp-mhalf 3.5d0 10.0d0)
+  (legendre-conicalp-0 3.5d0 10.0d0)
+  (legendre-conicalp-1 3.5d0 10.0d0)
+  (legendre-regular-spherical-conical 3 3.5d0 10.0d0)
+  (legendre-regular-cylindrical-conical 3 3.5d0 10.0d0)
+  (legendre-h3d-0 1.0d0 0.5d0)
+  (legendre-h3d-1 1.0d0 0.5d0)
+  (legendre-h3d 4 1.0d0 0.5d0)
+  (letm ((arr (vector-double 4)))
       (legendre-h3d-array 1.0d0 0.5d0 arr)
-      (data arr)))))
+      (data arr)))
+|#
+
+(LISP-UNIT:DEFINE-TEST LEGENDRE
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL (LIST 0.3d0 0.0d0)
+				     (MULTIPLE-VALUE-LIST
+				      (LEGENDRE-P1
+				       0.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.365d0 2.8199664825478977d-16)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-P2 0.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.38249999999999995d0 1.9984014443252816d-16)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-P3 0.3d0)))
+  (LISP-UNIT:ASSERT-ERROR 'GSL-ERROR
+			  (LEGENDRE-PL -4 0.3d0))
+  (LISP-UNIT:ASSERT-ERROR 'GSL-ERROR
+			  (LEGENDRE-PL 4 3.0d0))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.07293749999999999d0 5.668382430101814d-17)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-PL 4 0.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(1.0d0 0.5d0 -0.125d0 -0.4375d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 4)))
+      (LEGENDRE-PL-ARRAY 0.5d0 ARR) (DATA ARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.3128529498822064d0 1.3893461931245028d-16)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-Q0 3.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.03241473461128108d0 1.4395033881023292d-17)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-Q1 3.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.004026461384737812d0 1.788108054840004d-18)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-QL 2 3.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -34.099750274012266d0 3.0286662310541114d-14)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-PLM 4 3 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(2.25d0 5.625d0 4.21875d0 -4.921875d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 4)))
+      (LEGENDRE-PLM-ARRAY 2 0.5d0 ARR) (DATA ARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #(-3.0d0 3.75d0 33.75d0 55.78125d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM
+	((VAL (VECTOR-DOUBLE 4)) (DERIV (VECTOR-DOUBLE 4)))
+      (LEGENDRE-PLM-DERIV-ARRAY 2 0.5d0 VAL DERIV)
+      (DATA DERIV))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.30366280894310793d0 3.438761110552081d-15)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-SPHPLM 1200 1100 0.3d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.24892463950030283d0 0.4127948151484927d0
+      0.35120655562190445d0 0.051599351893561574d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 4)))
+      (LEGENDRE-SPHPLM-ARRAY 4 0.5d0 ARR)
+      (DATA ARR))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(-0.6637990386674741d0 -0.2751965434323283d0
+      1.2710332489173686d0 2.648766730536161d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM
+	((VAL (VECTOR-DOUBLE 4)) (DERIV (VECTOR-DOUBLE 4)))
+      (LEGENDRE-SPHPLM-DERIV-ARRAY 4 0.5d0 VAL DERIV)
+      (DATA DERIV))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.1255299048878925d0 1.3395992025650077d-15)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-CONICALP-HALF 3.5d0 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.06274336292793128d0 2.504525777730328d-16)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-CONICALP-MHALF 3.5d0 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.1316366189368757d0 3.0865532615549887d-15)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-CONICALP-0 3.5d0 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.17407195560076358d0 4.024312727070929d-15)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-CONICALP-1 3.5d0 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 8.980797952969897d-4 7.157154480497032d-17)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-REGULAR-SPHERICAL-CONICAL 3 3.5d0 10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.0023060250619859907d0 2.7345630541297237d-16)
+   (MULTIPLE-VALUE-LIST
+    (LEGENDRE-REGULAR-CYLINDRICAL-CONICAL 3 3.5d0
+					  10.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.9200342692589382d0 1.7018240558144874d-15)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-H3D-0 1.0d0 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.2169402645039212d0 1.418962379417407d-15)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-H3D-1 1.0d0 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.002400616238997978d0 4.3381136468045d-17)
+   (MULTIPLE-VALUE-LIST (LEGENDRE-H3D 4 1.0d0 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST
+    #(0.9200342692589379d0 0.21694026450392115d0
+      0.04795066048830776d0 0.010663769096144337d0))
+   (MULTIPLE-VALUE-LIST
+    (LETM ((ARR (VECTOR-DOUBLE 4)))
+      (LEGENDRE-H3D-ARRAY 1.0d0 0.5d0 ARR)
+      (DATA ARR)))))
diff --git a/special-functions/logarithm.lisp b/special-functions/logarithm.lisp
index a9960523dec42c576e0cd02b6874718ce14daec4..f0843c6bc566bf7e480e4e079992231a600dd2e7 100644
--- a/special-functions/logarithm.lisp
+++ b/special-functions/logarithm.lisp
@@ -1,65 +1,71 @@
-;********************************************************
-; file:        logarithm.lisp                            
-; description: Logarithm                                 
-; date:        Sun Apr 30 2006 - 22:08                   
-; author:      Liam M. Healy                             
-; modified:    Mon Oct  8 2007 - 11:27
-;********************************************************
-;;; $Id: $
+;; Logarithm
+;; Liam Healy, Sun Apr 30 2006 - 22:08
+;; Time-stamp: <2008-02-16 22:34:59EST logarithm.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
 (defgeneric gsl-log (x)
-  (:documentation
-   "The natural logarithm of @var{x}, @math{\log(x)}, for @math{x > 0}."))
+  (:documentation			; FDL
+   "The natural logarithm of x, log(x), for x > 0."))
 
-(defun-gsl gsl-log ((x float))
-  "gsl_sf_log_e" ((x :double) (ret sf-result))
+(defmfun gsl-log ((x float))
+  "gsl_sf_log_e"
+  ((x :double) (ret sf-result))
   :type :method
   :export t)
 
-(defun-gsl gsl-log ((x complex))
+(defmfun gsl-log ((x complex))
   "gsl_sf_complex_log_e"
   (((realpart x) :double) ((imagpart x) :double)
    (re-ret sf-result) (im-ret sf-result))
   :type :method
-  :documentation "Results are returned as @var{lnr}, @var{theta} such that
-  @math{\exp(lnr + i \theta) = z_r + i z_i}, where @math{\theta} lies in
-  the range @math{[-\pi,\pi]}."
   :return
-  ((complex (val re-ret) (val im-ret)) (complex (err re-ret) (err im-ret))))
+  ((complex (val re-ret) (val im-ret)) (complex (err re-ret) (err im-ret)))
+  :documentation			; FDL
+  "Results are returned as lnr, theta such that
+  exp(lnr + i \theta) = z_r + i z_i, where theta lies in the range [-\pi,\pi].")
 
-(defun-gsl log-abs (x)
+(defmfun log-abs (x)
   "gsl_sf_log_abs_e" ((x :double) (ret sf-result))
-  :documentation
-  "The natural logarithm of the magnitude of @var{x},
-  @math{\log(|x|)}, for @math{x \ne 0}.")
+  :documentation			; FDL
+  "The natural logarithm of the magnitude of x, log(|x|), for x ne 0.")
 
-(defun-gsl log-1+x (x)
+(defmfun log-1+x (x)
   "gsl_sf_log_1plusx_e" ((x :double) (ret sf-result))
-  :documentation
-  "@math{\log(1 + x)} for @math{x > -1} using an
-   algorithm that is accurate for small @math{x}.")
+  :documentation			; FDL
+  "log(1 + x) for x > -1 using an algorithm that is accurate for small x.")
 
-(defun-gsl log-1+x-m1 (x)
+(defmfun log-1+x-m1 (x)
   "gsl_sf_log_1plusx_mx_e" ((x :double) (ret sf-result))
-  :documentation
-  "@math{\log(1 + x) - x} for @math{x > -1} using an
-  algorithm that is accurate for small @math{x}.")
+  :documentation			; FDL
+  "log(1 + x) - x for x > -1 using an algorithm that is accurate for small x.")
 
-(lisp-unit:define-test logarithm
-  (lisp-unit:assert-first-fp-equal
-   "0.693147180560d+00"
-   (gsl-log 2.0d0))
-  (lisp-unit:assert-equal
-   '("0.346573590280d+00" "0.785398163397d+00")
-   (lisp-unit::fp-string (gsl-log #C(1.0d0 1.0d0))))
-  (lisp-unit:assert-first-fp-equal
-   "0.693147180560d+00"
-   (log-abs -2.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.999950003333d-04"
-   (log-1+x 1.d-4))
-  (lisp-unit:assert-first-fp-equal
-   "-0.499966669166d-08"
-   (log-1+x-m1 1.d-4)))
+;;; Examples and unit test
+
+#|
+(make-tests logarithm
+  (gsl-log 2.0d0)
+  (gsl-log #C(1.0d0 1.0d0))
+  (log-abs -2.0d0)
+  (log-1+x 1.d-4)
+  (log-1+x-m1 1.d-4))
+|#
+
+(LISP-UNIT:DEFINE-TEST LOGARITHM
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6931471805599453d0 3.078191837246648d-16)
+   (MULTIPLE-VALUE-LIST (GSL-LOG 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(0.34657359027997264d0 0.7853981633974483d0)
+	 #C(1.539095918623324d-16 7.69547959311662d-17))
+   (MULTIPLE-VALUE-LIST (GSL-LOG #C(1.0d0 1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6931471805599453d0 3.078191837246648d-16)
+   (MULTIPLE-VALUE-LIST (LOG-ABS -2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 9.999500033330834d-5 2.2203350343487824d-20)
+   (MULTIPLE-VALUE-LIST (LOG-1+X 1.d-4)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -4.999666691664667d-9 1.1101490153075193d-24)
+   (MULTIPLE-VALUE-LIST (LOG-1+X-M1 1.d-4))))
diff --git a/special-functions/power.lisp b/special-functions/power.lisp
index 932a58f51fe3a6b94ff733b7125e145f9d57f6e5..1655460f096669830dfd916a0e6e50f9f0882812 100644
--- a/special-functions/power.lisp
+++ b/special-functions/power.lisp
@@ -1,23 +1,25 @@
-;********************************************************
-; file:        power.lisp                                
-; description: Power                                     
-; date:        Sun Apr 30 2006 - 22:46                   
-; author:      Liam M. Healy                             
-; modified:    Sat Jun 17 2006 - 22:26
-;********************************************************
-;;; $Id: $
+;; Integer powers
+;; Liam Healy, Sun Apr 30 2006 - 22:46
+;; Time-stamp: <2008-02-16 22:39:17EST power.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl pow (x n)
+(defmfun pow (x n)
   "gsl_sf_pow_int_e" ((x :double) (n :int) (ret sf-result))
-  :documentation "The power @math{x^n} for integer @var{n}.  The
+  :documentation			; FDL
+  "The power x^n for integer n.  The
   power is computed using the minimum number of multiplications. For
-  example, @math{x^8} is computed as @math{((x^2)^2)^2}, requiring only 3
+  example, x^8 is computed as ((x^2)^2)^2, requiring only 3
   multiplications.  For reasons of efficiency, these functions do not
   check for overflow or underflow conditions.")
 
-(lisp-unit:define-test power
-  (lisp-unit:assert-first-fp-equal
-   "0.525218750000d+03"
-   (pow 3.5d0 5)))
+#|
+(make-tests power
+  (pow 3.5d0 5))
+|#
+
+(LISP-UNIT:DEFINE-TEST POWER
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 525.21875d0 9.329759187437503d-13)
+   (MULTIPLE-VALUE-LIST (POW 3.5d0 5))))
diff --git a/special-functions/psi.lisp b/special-functions/psi.lisp
index e5ae44b8dfbd31a1a453dfb7a0d4bac0e2cf6691..01ff70a447baa35bccfd7cdaa2087d0eca0a0b46 100644
--- a/special-functions/psi.lisp
+++ b/special-functions/psi.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        psi.lisp                                  
-; description: Psi (digamma) functions                   
-; date:        Mon May  1 2006 - 22:11                   
-; author:      Liam M. Healy                             
-; modified:    Mon Oct  8 2007 - 11:30
-;********************************************************
-;;; $Id: $
+;; Psi (digamma) functions
+;; Liam Healy, Mon May  1 2006 - 22:11
+;; Time-stamp: <2008-02-16 22:43:28EST psi.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -14,59 +10,87 @@
 ;;;;****************************************************************************
 
 (defgeneric psi (x)
+  ;; FDL
   (:documentation "The psi, or digamma, function."))
 
-(defun-gsl psi  ((n integer))
+(defmfun psi  ((n integer))
   "gsl_sf_psi_int_e" ((n :int) (ret sf-result))
   :type :method
   :export t
-  :documentation "Domain: n integer, n > 0.")
+  :documentation			; FDL
+  "Domain: n integer, n > 0.")
 
-(defun-gsl psi ((x float))
+(defmfun psi ((x float))
   "gsl_sf_psi_e" ((x :double) (ret sf-result))
   :type :method 
-  :documentation "Domain: x /= 0.0, -1.0, -2.0, ...")
+  :documentation			; FDL
+  "Domain: x /= 0.0, -1.0, -2.0, ...")
 
-(defun-gsl psi-1+iy (x)
+(defmfun psi-1+iy (x)
   "gsl_sf_psi_1piy_e" ((x :double) (ret sf-result))
-  :documentation "The real part of the digamma function
-  on the line @math{1+i y}, @math{\Re[\psi(1 + i y)]}.")
+  :documentation			; FDL
+  "The real part of the digamma function
+  on the line 1+i y, Re[psi(1 + i y)].")
 
 ;;;;****************************************************************************
 ;;;; Trigamma Function
 ;;;;****************************************************************************
 
 (defgeneric psi-1 (x)
+  ;; FDL
   (:documentation "The Trigamma function."))
 
-(defun-gsl psi-1 ((n integer))
+(defmfun psi-1 ((n integer))
   "gsl_sf_psi_1_int_e" ((n :int) (ret sf-result))
   :type :method 
-  :documentation "Domain: n integer, n > 0.")
+  :documentation			; FDL
+  "Domain: n integer, n > 0.")
 
-(defun-gsl psi-1 ((x float))
+(defmfun psi-1 ((x float))
   "gsl_sf_psi_1_e" ((x :double) (ret sf-result))
   :type :method
-  :documentation "Domain: x /= 0.0, -1.0, -2.0, ...")
+  :documentation			; FDL
+  "Domain: x /= 0.0, -1.0, -2.0, ...")
 
 ;;;;****************************************************************************
 ;;;; Polygamma
 ;;;;****************************************************************************
 
-(defun-gsl psi-n (m x)
+(defmfun psi-n (m x)
   "gsl_sf_psi_n_e" ((m :int) (x :double) (ret sf-result))
-  :documentation "The polygamma function @math{\psi^@{(m)@}(x)} for
-   @math{m >= 0}, @math{x > 0}.")
+  :documentation			; FDL
+  "The polygamma function psi^{(m)}(x)} for m >= 0, x > 0.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test psi
-  (lisp-unit:assert-first-fp-equal "0.125611766843d+01" (psi 4))
-  (lisp-unit:assert-first-fp-equal "0.125611766843d+01" (psi 4.0d0))
-  (lisp-unit:assert-first-fp-equal "0.714591515374d+00" (psi-1+iy 2.0d0))
-  (lisp-unit:assert-first-fp-equal "0.283822955737d+00" (psi-1 4))
-  (lisp-unit:assert-first-fp-equal "0.283822955737d+00" (psi-1 4.0d0))
-  (lisp-unit:assert-first-fp-equal "-0.800397322451d-01" (psi-n 2 4.0d0)))
+#|
+(make-tests psi
+  (psi 4)
+  (psi 4.0d0)
+  (psi-1+iy 2.0d0)
+  (psi-1 4)
+  (psi-1 4.0d0)
+  (psi-n 2 4.0d0))
+|#
 
+(LISP-UNIT:DEFINE-TEST PSI
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.2561176684318005d0 2.789141514262906d-16)
+   (MULTIPLE-VALUE-LIST (PSI 4)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.2561176684318005d0 2.846229783626858d-16)
+   (MULTIPLE-VALUE-LIST (PSI 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7145915153739806d0 1.925418559790263d-14)
+   (MULTIPLE-VALUE-LIST (PSI-1+IY 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.2838229557371153d0 6.302135607530242d-17)
+   (MULTIPLE-VALUE-LIST (PSI-1 4)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.28382295573711525d0 1.764597970108467d-15)
+   (MULTIPLE-VALUE-LIST (PSI-1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.0800397322451145d0 5.222647053360405d-16)
+   (MULTIPLE-VALUE-LIST (PSI-N 2 4.0d0))))
diff --git a/special-functions/return-structures.lisp b/special-functions/return-structures.lisp
index 60cbc08397ee53b532dce1f58961e4e69f01b26f..453cf7d0f0a5dcc5d947f9d595c12cb9c1ca25c4 100644
--- a/special-functions/return-structures.lisp
+++ b/special-functions/return-structures.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        return-structures.lisp                    
-; description: Structures returned by special functions. 
-; date:        Mon Jan  1 2007 - 11:35                   
-; author:      Liam M. Healy                             
-; modified:    Mon Jan  1 2007 - 11:35
-;********************************************************
-;;; $Id: $
+;; Structures returned by special functions.
+;; Liam Healy, Mon Jan  1 2007 - 11:35
+;; Time-stamp: <2008-02-16 22:45:31EST return-structures.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -21,7 +17,7 @@
 
 (cffi:defcstruct sf-result-e10
   "Results from special functions with value, error estimate
-and a scaling exponent e10, such that the value is val*10^e10."
+   and a scaling exponent e10, such that the value is val*10^e10."
   ;; file:///usr/share/doc/gsl-ref-html/gsl-ref_7.html#SEC61
   (val :double)
   (err :double)
@@ -43,13 +39,3 @@ and a scaling exponent e10, such that the value is val*10^e10."
 (defun e10 (sf-result)
   (cffi:foreign-slot-value sf-result 'sf-result-e10 'e10))
 
-;;; obsolete; just do this manually for the few cases of multiple sf-result returns?
-(defun rearrange-sf-result-err (return-list)
-  "Put the 'err values from the sf-results at the end of the return values."
-  (flet ((sf-err (x)
-	   (and (eq (first x) 'cffi:foreign-slot-value)
-		(member (third x) '('sf-result 'sf-result-e10) :test #'equal)
-		(equal (fourth x) ''err))))
-    (append
-     (remove-if #'sf-err return-list)
-     (remove-if-not #'sf-err return-list))))
diff --git a/special-functions/synchrotron.lisp b/special-functions/synchrotron.lisp
index ae5c3b15e16cedf914a2e6f6af2cfe5aa5c2a5d7..35fa515bf625a1f740aa75137ff82b2b3402eefb 100644
--- a/special-functions/synchrotron.lisp
+++ b/special-functions/synchrotron.lisp
@@ -1,28 +1,34 @@
-;********************************************************
-; file:        synchrotron.lisp                                
-; description: Synchrotron functions
-; date:        Mon May  1 2006 - 22:29
-; author:      Liam M. Healy                             
-; modified:    Sat Jun 17 2006 - 22:32
-;********************************************************
-;;; $Id: $
+;; Synchrotron functions
+;; Liam Healy, Mon May  1 2006 - 22:29
+;; Time-stamp: <2008-02-16 22:48:24EST synchrotron.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl synchrotron-1 (x)
+(defmfun synchrotron-1 (x)
   "gsl_sf_synchrotron_1_e" ((x :double) (ret sf-result))
-  :documentation "The first synchrotron function 
-  @math{x \int_x^\infty dt K_@{5/3@}(t)} for @math{x >= 0}.")
+  :documentation			; FDL
+  "The first synchrotron function x \int_x^\infty dt K_{5/3}(t)} for x >= 0.")
 
-(defun-gsl synchrotron-2 (x)
+(defmfun synchrotron-2 (x)
   "gsl_sf_synchrotron_2_e" ((x :double) (ret sf-result))
-  :documentation "The second synchrotron function 
-  @math{x K_@{2/3@}(x)} for @math{x >= 0}.")
+  :documentation			; FDL
+  "The second synchrotron function x K_{2/3}(x)} for x >= 0.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test synchrotron
-  (lisp-unit:assert-first-fp-equal "0.528273966979d-01" (synchrotron-1 4.0d0))
-  (lisp-unit:assert-first-fp-equal "0.469232058261d-01" (synchrotron-2 4.0d0)))
+#|
+(make-tests synchrotron
+  (synchrotron-1 4.0d0)
+  (synchrotron-2 4.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST SYNCHROTRON
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.052827396697912476d0 4.825849878208132d-14)
+   (MULTIPLE-VALUE-LIST (SYNCHROTRON-1 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.04692320582614684d0 6.854449168174663d-14)
+   (MULTIPLE-VALUE-LIST (SYNCHROTRON-2 4.0d0))))
diff --git a/special-functions/transport.lisp b/special-functions/transport.lisp
index 6be0ba459e40bb2c490ec5a3fe15cfc39e3b2591..0f30950d84597797522c9ccc415e23fba7e29d7e 100644
--- a/special-functions/transport.lisp
+++ b/special-functions/transport.lisp
@@ -1,41 +1,57 @@
-;********************************************************
-; file:        transport.lisp                                
-; description: Transport functions
-; date:        Mon May  1 2006 - 22:29
-; author:      Liam M. Healy                             
-; modified:    Sat Jun 17 2006 - 22:33
-;********************************************************
-;;; $Id: $
+;; Transport functions
+;; Liam Healy, Mon May  1 2006 - 22:29
+;; Time-stamp: <2008-02-16 22:50:39EST transport.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
-;;; The transport functions @math{J(n,x)} are defined by the integral 
+;;;  FDL
+;;; The transport functions J(n,x) are defined by the integral 
 ;;; representations
-;;; @c{$J(n,x) := \int_0^x dt \, t^n e^t /(e^t - 1)^2$}
-;;; @math{J(n,x) := \int_0^x dt t^n e^t /(e^t - 1)^2}.
+;;; J(n,x) := \int_0^x dt \, t^n e^t /(e^t - 1)^2.
 
-(defun-gsl transport-2 (x)
+(defmfun transport-2 (x)
   "gsl_sf_transport_2_e" ((x :double) (ret sf-result))
-  :documentation "The transport function @math{J(2,x)}.")
+  :documentation			; FDL
+  "The transport function J(2,x).")
 
-(defun-gsl transport-3 (x)
+(defmfun transport-3 (x)
   "gsl_sf_transport_3_e" ((x :double) (ret sf-result))
-  :documentation "The transport function @math{J(3,x)}.")
+  :documentation			; FDL
+  "The transport function J(3,x).")
 
-(defun-gsl transport-4 (x)
+(defmfun transport-4 (x)
   "gsl_sf_transport_4_e" ((x :double) (ret sf-result))
-  :documentation "The transport function @math{J(4,x)}.")
+  :documentation			; FDL
+  "The transport function J(4,x).")
 
-(defun-gsl transport-5 (x)
+(defmfun transport-5 (x)
   "gsl_sf_transport_5_e" ((x :double) (ret sf-result))
-  :documentation "The transport function @math{J(5,x)}.")
+  :documentation			; FDL
+  "The transport function J(5,x).")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test transport
-  (lisp-unit:assert-first-fp-equal "0.280666640456d+01" (transport-2 4.0d0))
-  (lisp-unit:assert-first-fp-equal "0.457921743723d+01" (transport-3 4.0d0))
-  (lisp-unit:assert-first-fp-equal "0.107319323930d+02" (transport-4 4.0d0))
-  (lisp-unit:assert-first-fp-equal "0.294883390152d+02" (transport-5 4.0d0)))
+#|
+(make-tests transport
+  (transport-2 4.0d0)
+  (transport-3 4.0d0)
+  (transport-4 4.0d0)
+  (transport-5 4.0d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST TRANSPORT
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.806666404563118d0 2.2867923780257255d-15)
+   (MULTIPLE-VALUE-LIST (TRANSPORT-2 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 4.579217437229157d0 3.242324689309112d-15)
+   (MULTIPLE-VALUE-LIST (TRANSPORT-3 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 10.731932392998623d0 1.0925209116254758d-14)
+   (MULTIPLE-VALUE-LIST (TRANSPORT-4 4.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 29.488339015245842d0 3.204450601879883d-14)
+   (MULTIPLE-VALUE-LIST (TRANSPORT-5 4.0d0))))
diff --git a/special-functions/trigonometry.lisp b/special-functions/trigonometry.lisp
index 15e10b81e7577769552a8de7690a24faab2b8543..4c15d31992d4241f73e7a1e3c86fe1249a3c88d2 100644
--- a/special-functions/trigonometry.lisp
+++ b/special-functions/trigonometry.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        trigonometry.lisp                         
-; description: Trigonometry                              
-; date:        Thu May  4 2006 - 22:58                   
-; author:      Liam M. Healy                             
-; modified:    Mon Oct  8 2007 - 11:27
-;********************************************************
-;;; $Id: $
+;; Trigonometry
+;; Liam Healy, Thu May  4 2006 - 22:58
+;; Time-stamp: <2008-02-16 22:59:12EST trigonometry.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -19,12 +15,12 @@
 (defgeneric gsl-cos (x)
   (:documentation "The cosine function @math{\sin(x)}."))
 
-(defun-gsl gsl-sin ((x float))
+(defmfun gsl-sin ((x float))
   "gsl_sf_sin_e" ((x :double) (ret sf-result))
   :type :method
   :export t)
 
-(defun-gsl gsl-sin ((x complex))
+(defmfun gsl-sin ((x complex))
   "gsl_sf_complex_sin_e"
   (((realpart x) :double) ((imagpart x) :double)
    (re-ret sf-result) (im-ret sf-result))
@@ -32,12 +28,12 @@
   :return ((complex (val re-ret) (val im-ret))
 	   (complex (err re-ret) (err im-ret))))
 
-(defun-gsl gsl-cos ((x float))
+(defmfun gsl-cos ((x float))
   "gsl_sf_cos_e" ((x :double) (ret sf-result))
   :type :method
   :export t)
 
-(defun-gsl gsl-cos ((x complex))
+(defmfun gsl-cos ((x complex))
   "gsl_sf_complex_cos_e"
   (((realpart x) :double) ((imagpart x) :double)
    (re-ret sf-result) (im-ret sf-result))
@@ -45,21 +41,24 @@
   :return ((complex (val re-ret) (val im-ret))
 	   (complex (err re-ret) (err im-ret))))
 
-(defun-gsl hypotenuse (x y)
+(defmfun hypotenuse (x y)
   "gsl_sf_hypot_e" ((x :double) (y :double) (ret sf-result))
-  :documentation "The hypotenuse function @math{\sqrt@{x^2 + y^2@}}.")
+  :documentation			; FDL
+  "The hypotenuse function sqrt{x^2 + y^2}.")
 
-(defun-gsl sinc (x)
+(defmfun sinc (x)
   "gsl_sf_sinc_e" ((x :double) (ret sf-result))
-  :documentation "@math{\sinc(x) = \sin(\pi x) / (\pi x)}")
+  :documentation			; FDL
+  "sinc(x) = sin(pi x) / (pi x)}")
 
-(defun-gsl log-sin (x)
+(defmfun log-sin (x)
   "gsl_sf_complex_logsin_e"
   (((realpart x) :double) ((imagpart x) :double)
    (re-ret sf-result) (im-ret sf-result))
-  :documentation "This function computes the logarithm of the complex sine,
-  @math{\log(\sin(z_r + i z_i))} storing the real and imaginary parts in
-  @var{szr}, @var{szi}."
+  :documentation			; FDL
+  "This function computes the logarithm of the complex sine,
+  \log(\sin(z_r + i z_i)) storing the real and imaginary parts in
+  szr, szi."
   :return ((complex (val re-ret) (val im-ret))
 	   (complex (err re-ret) (err im-ret))))
 
@@ -67,104 +66,145 @@
 ;;;; Hyperbolic Trigonometric Functions
 ;;;;****************************************************************************
 
-(defun-gsl log-sinh (x)
+(defmfun log-sinh (x)
   "gsl_sf_lnsinh_e" ((x :double) (ret sf-result))
-  :documentation "Logarithm of sinh function, special functions
-  These routines compute @math{\log(\sinh(x))} for @math{x > 0}.")
+  :documentation			; FDL
+  "Logarithm of sinh function, special functions
+  These routines compute log(\sinh(x)) for x > 0.")
 
-(defun-gsl log-cosh (x)
+(defmfun log-cosh (x)
   "gsl_sf_lncosh_e" ((x :double) (ret sf-result))
-  :documentation "Logarithm of cosh function, special functions
-  These routines compute @math{\log(\cosh(x))} for any @var{x}.")
+  :documentation			; FDL
+  "Logarithm of cosh function, special functions
+   These routines compute log(cosh(x)) for any x.")
 
 ;;;;****************************************************************************
 ;;;; Conversion Functions
 ;;;;****************************************************************************
 
-(defun-gsl polar-to-rectangular (r theta)
+(defmfun polar-to-rectangular (r theta)
   "gsl_sf_polar_to_rect"
   ((r :double) (theta :double) (x sf-result) (y sf-result))
-  :documentation "Convert the polar coordinates (@var{r},@var{theta}) to
-  rectilinear coordinates (@var{x},@var{y}), @math{x = r\cos(\theta)},
-  @math{y = r\sin(\theta)}."
-  :return ((val x) (val y) (err x) (err y)))
+  :return ((val x) (val y) (err x) (err y))
+  :documentation			; FDL
+  "Convert the polar coordinates (r, theta) to
+  rectilinear coordinates (x, y), x = r\cos(\theta), y = r\sin(\theta).")
 
-(defun-gsl rectangular-to-polar (x y)
+(defmfun rectangular-to-polar (x y)
   "gsl_sf_rect_to_polar"
   ((x :double) (y :double) (r sf-result) (theta sf-result))
-  :documentation "Convert the rectilinear coordinates (@var{x},@var{y}) to
-  polar coordinates (@var{r},@var{theta}), such that @math{x =
-  r\cos(\theta)}, @math{y = r\sin(\theta)}.  The argument @var{theta}
-  lies in the range @math{[-\pi, \pi]}."
-  :return ((val r) (val theta) (err r) (err theta)))
+  :return ((val r) (val theta) (err r) (err theta))
+  :documentation			; FDL
+  "Convert the rectilinear coordinates (x, y) to
+  polar coordinates (r, theta), such that x =
+  r cos(theta)}, y = r sin(theta).  The argument theta
+  lies in the range [-\pi, \pi].")
 
 ;;;;****************************************************************************
 ;;;; Restriction Functions
 ;;;;****************************************************************************
 
-(defun-gsl restrict-symmetric (theta)
+(defmfun restrict-symmetric (theta)
   "gsl_sf_angle_restrict_symm" ((theta :double))
   :c-return :double
-  :documentation "Force the angle @var{theta} to lie in the range
-  @math{(-\pi,\pi]}.")
+  :documentation			; FDL
+  "Force the angle theta to lie in the range (-\pi,\pi].")
 
-(defun-gsl restrict-positive (theta)
+(defmfun restrict-positive (theta)
   "gsl_sf_angle_restrict_pos" ((theta :double))
   :c-return :double
-  :documentation "Force the angle @var{theta} to lie in
-  the range @math{[0,2\pi)}.")
+  :documentation			; FDL
+  "Force the angle theta to lie in the range [0,2\pi).")
 
 ;;;;****************************************************************************
 ;;;; Trigonometric Functions With Error Estimates
 ;;;;****************************************************************************
 
-(defun-gsl sin-err (x dx)
+(defmfun sin-err (x dx)
   "gsl_sf_sin_err_e" ((x :double) (dx :double) (ret sf-result))
-  :documentation "Compute the sine of an angle @var{x} with
-   an associated absolute error @var{dx}, @math{\sin(x \pm dx)}.")
+  :documentation			; FDL
+  "Compute the sine of an angle x with
+   an associated absolute error dx, sin(x \pm dx).")
 
-(defun-gsl cos-err (x dx)
+(defmfun cos-err (x dx)
   "gsl_sf_cos_err_e" ((x :double) (dx :double) (ret sf-result))
-  :documentation "The cosine of an angle @var{x} with an associated
-  absolute error @var{dx}, @math{\cos(x \pm dx)}.")
+  :documentation			; FDL
+  "The cosine of an angle x with an associated
+  absolute error dx, cos(x \pm dx).")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test trigonometry
-  (lisp-unit:assert-first-fp-equal "0.841470984808d+00" (gsl-sin 1.0d0))
-  (lisp-unit:assert-equal
-   '("0.129845758142d+01" "0.634963914785d+00")
-   (lisp-unit::fp-string (gsl-sin #C(1.0d0 1.0d0))))
-  (lisp-unit:assert-first-fp-equal "0.540302305868d+00" (gsl-cos 1.0d0))
-  (lisp-unit:assert-equal
-   '("0.833730025131d+00" "-0.988897705763d+00")
-   (lisp-unit::fp-string (gsl-cos #C(1.0d0 1.0d0))))
-  (lisp-unit:assert-first-fp-equal
-   "0.223606797750d+01"
-   (hypotenuse 1.0d0 2.0d0))
-  (lisp-unit:assert-first-fp-equal "0.636619772368d+00" (sinc 0.5d0))
-  (lisp-unit:assert-equal
-   '("0.368383731425d+00" "0.454820233310d+00")
-   (lisp-unit::fp-string (log-sin #C(1.0d0 1.0d0))))
-  (lisp-unit:assert-first-fp-equal "-0.651822325947d+00" (log-sinh 0.5d0))
-  (lisp-unit:assert-first-fp-equal "0.120114506958d+00" (log-cosh 0.5d0))
-  (lisp-unit:assert-equal
-   '("0.108060461174d+01" "0.168294196962d+01")
-   (lisp-unit::fp-sequence
-    (subseq (multiple-value-list
-	     (polar-to-rectangular 2.0d0 1.0d0))
-	    0 2)))
-  (lisp-unit:assert-equal
-   '("0.223606797750d+01" "0.463647609001d+00")
-   (lisp-unit::fp-sequence
-    (subseq (multiple-value-list (rectangular-to-polar 2.0d0 1.0d0)) 0 2)))
-  (lisp-unit:assert-first-fp-equal
-   "-0.128318530718d+01"
-   (restrict-symmetric 5.0d0))
-  (lisp-unit:assert-first-fp-equal
-   "0.528318530718d+01"
-   (restrict-positive -1.0d0))
-  (lisp-unit:assert-first-fp-equal "0.479425538604d+00" (sin-err 0.5d0 0.01d0))
-  (lisp-unit:assert-first-fp-equal "0.877582561890d+00" (cos-err 0.5d0 0.01d0)))
+#|
+(make-tests trigonometry
+  (gsl-sin 1.0d0)
+  (gsl-sin #C(1.0d0 1.0d0))
+  (gsl-cos 1.0d0)
+  (gsl-cos #C(1.0d0 1.0d0))
+  (hypotenuse 1.0d0 2.0d0)
+  (sinc 0.5d0)
+  (log-sin #C(1.0d0 1.0d0))
+  (log-sinh 0.5d0)
+  (log-cosh 0.5d0)
+  (polar-to-rectangular 2.0d0 1.0d0)
+  (rectangular-to-polar 2.0d0 1.0d0)
+  (restrict-symmetric 5.0d0)
+  (restrict-positive -1.0d0)
+  (sin-err 0.5d0 0.01d0)
+  (cos-err 0.5d0 0.01d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST TRIGONOMETRY
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8414709848078965d0 3.736881847550928d-16)
+   (MULTIPLE-VALUE-LIST (GSL-SIN 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(1.2984575814159773d0 0.6349639147847361d0)
+	 #C(5.766310013548447d-16 2.8198062320005596d-16))
+   (MULTIPLE-VALUE-LIST (GSL-SIN #C(1.0d0 1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.5403023058681398d0 2.3994242409314904d-16)
+   (MULTIPLE-VALUE-LIST (GSL-COS 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(0.8337300251311491d0 -0.9888977057628651d0)
+	 #C(3.7025050808876487d-16 4.3915880077477046d-16))
+   (MULTIPLE-VALUE-LIST (GSL-COS #C(1.0d0 1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.23606797749979d0 9.930136612989092d-16)
+   (MULTIPLE-VALUE-LIST (HYPOTENUSE 1.0d0 2.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6366197723675814d0 3.0072729231305663d-16)
+   (MULTIPLE-VALUE-LIST (SINC 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST #C(0.3683837314249251d0 0.4548202333099499d0)
+	 #C(1.6359524021011266d-16 8.179762010505633d-17))
+   (MULTIPLE-VALUE-LIST (LOG-SIN #C(1.0d0 1.0d0))))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -0.6518223259470272d0 2.8946726169244526d-16)
+   (MULTIPLE-VALUE-LIST (LOG-SINH 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.12011450695827751d0 4.401938023864669d-17)
+   (MULTIPLE-VALUE-LIST (LOG-COSH 0.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.0806046117362795d0 1.682941969615793d0
+	 8.535730329413909d-16 9.873187936033346d-16)
+   (MULTIPLE-VALUE-LIST
+    (POLAR-TO-RECTANGULAR 2.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.23606797749979d0 0.4636476090008061d0
+	 9.930136612989092d-16 2.0590090033003876d-16)
+   (MULTIPLE-VALUE-LIST
+    (RECTANGULAR-TO-POLAR 2.0d0 1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST -1.2831853071795862d0)
+   (MULTIPLE-VALUE-LIST (RESTRICT-SYMMETRIC 5.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 5.283185307179586d0)
+   (MULTIPLE-VALUE-LIST (RESTRICT-POSITIVE -1.0d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.479425538604203d0 0.008775825618904047d0)
+   (MULTIPLE-VALUE-LIST (SIN-ERR 0.5d0 0.01d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.8775825618903728d0 0.004794255386042614d0)
+   (MULTIPLE-VALUE-LIST (COS-ERR 0.5d0 0.01d0))))
diff --git a/special-functions/zeta.lisp b/special-functions/zeta.lisp
index 6cbb365bb5bcc802d075ddce9a646e1eacd84287..86a840cd5da363e4302a0eaef8a1d2c6e86e32e4 100644
--- a/special-functions/zeta.lisp
+++ b/special-functions/zeta.lisp
@@ -1,11 +1,7 @@
-;********************************************************
-; file:        zeta.lisp                         
-; description: Zeta functions                              
-; date:        Sat May 13 2006 - 23:27
-; author:      Liam M. Healy                             
-; modified:    Mon Oct  8 2007 - 11:30
-;********************************************************
-;;; $Id: $
+;; Zeta functions
+;; Liam Healy, Sat May 13 2006 - 23:27
+;; Time-stamp: <2008-02-16 23:04:22EST zeta.lisp>
+;; $Id: $
 
 (in-package :gsl)
 
@@ -13,24 +9,26 @@
 ;;;; Riemann Zeta Function
 ;;;;****************************************************************************
 
+;;; FDL
 ;;; The Riemann zeta function is defined by the infinite sum 
-;;; @math{\zeta(s) = \sum_@{k=1@}^\infty k^@{-s@}}.  
+;;; zeta(s) = \sum_{k=1}^\infty k^{-s}.  
 
 (defgeneric zeta (x)
-  (:documentation "The Riemann zeta function @math{\zeta(n)}."))
+  (:documentation			; FDL
+   "The Riemann zeta function zeta(n)."))
 
-(defun-gsl zeta ((n integer))
+(defmfun zeta ((n integer))
   "gsl_sf_zeta_int_e" ((n :int) (ret sf-result))
   :type :method
   :export t
-  :documentation "The Riemann zeta function @math{\zeta(n)} 
-   for integer @var{n}, @math{n \ne 1}.")
+  :documentation			; FDL
+  "The Riemann zeta function zeta(n) for integer n, n \ne 1.")
 
-(defun-gsl zeta ((s float))
+(defmfun zeta ((s float))
   "gsl_sf_zeta_e" ((s :double) (ret sf-result))
-  :documentation "The Riemann zeta function @math{\zeta(s)}
-   for arbitrary @var{s}, @math{s \ne 1}."
-  :type :method)
+  :type :method
+  :documentation			; FDL
+  "The Riemann zeta function zeta(s) for arbitrary s, s \ne 1.")
 
 ;;;;****************************************************************************
 ;;;; Riemann Zeta Function Minus One
@@ -39,45 +37,67 @@
 (defgeneric zeta-1 (x)
   (:documentation "zeta - 1."))
 
-(defun-gsl zeta-1 ((n integer))
+(defmfun zeta-1 ((n integer))
   "gsl_sf_zetam1_int_e" ((n :int) (ret sf-result))
-  :documentation "The Riemann zeta function @math{\zeta(n)} 
-   for integer @var{n}, @math{n \ne 1}."
   :type :method
-  :export t)
+  :export t
+  :documentation			; FDL
+  "The Riemann zeta function zeta(n) for integer n, n \ne 1.")
 
-(defun-gsl zeta-1 ((s float))
+(defmfun zeta-1 ((s float))
   "gsl_sf_zetam1_e" ((s :double) (ret sf-result))
-  :documentation "The Riemann zeta function @math{\zeta(s)}
-   for arbitrary @var{s}, @math{s \ne 1}."
-  :type :method)
+  :type :method
+  :documentation			; FDL
+  "The Riemann zeta function zeta(s) for arbitrary s, s \ne 1.")
 
 ;;;;****************************************************************************
 ;;;; Hurwitz Zeta Function
 ;;;;****************************************************************************
 
-(defun-gsl hurwitz-zeta (s q)
+(defmfun hurwitz-zeta (s q)
   "gsl_sf_hzeta_e" ((s :double) (q :double) (ret sf-result))
-  :documentation "The Hurwitz zeta function @math{\zeta(s,q)} for
-  @math{s > 1}, @math{q > 0}.")
+  :documentation			; FDL
+  "The Hurwitz zeta function zeta(s,q) for s > 1, q > 0.")
 
 ;;;;****************************************************************************
 ;;;; Eta Function
 ;;;;****************************************************************************
 
-(defun-gsl eta (s)
+(defmfun eta (s)
   "gsl_sf_eta_e" ((s :double) (ret sf-result))
-  :documentation "The eta function @math{\eta(s)} for arbitrary @var{s}.")
+  :documentation			; FDL
+  "The eta function eta(s) for arbitrary s.")
 
 ;;;;****************************************************************************
 ;;;; Examples and unit test
 ;;;;****************************************************************************
 
-(lisp-unit:define-test zeta
-  (lisp-unit:assert-first-fp-equal "0.164493406685d+01" (zeta 2))
-  (lisp-unit:assert-first-fp-equal "0.261237534869d+01" (zeta 1.5d0))
-  (lisp-unit:assert-first-fp-equal "0.644934066848d+00" (zeta-1 2))
-  (lisp-unit:assert-first-fp-equal "0.161237534869d+01" (zeta-1 1.5d0))
-  (lisp-unit:assert-first-fp-equal "0.140377976886d+01"
-				   (hurwitz-zeta 1.5d0 2.5d0))
-  (lisp-unit:assert-first-fp-equal "0.765147024625d+00" (eta 1.5d0)))
+#|
+(make-tests zeta
+  (zeta 2)
+  (zeta 1.5d0)
+  (zeta-1 2)
+  (zeta-1 1.5d0)
+  (hurwitz-zeta 1.5d0 2.5d0)
+  (eta 1.5d0))
+|#
+
+(LISP-UNIT:DEFINE-TEST ZETA
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.6449340668482264d0 7.304974700020789d-16)
+   (MULTIPLE-VALUE-LIST (ZETA 2)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 2.612375348685493d0 1.419471177334903d-14)
+   (MULTIPLE-VALUE-LIST (ZETA 1.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.6449340668482264d0 2.8640826015201633d-16)
+   (MULTIPLE-VALUE-LIST (ZETA-1 2)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.612375348685493d0 1.419471177334903d-14)
+   (MULTIPLE-VALUE-LIST (ZETA-1 1.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 1.4037797688568256d0 8.104244828616706d-15)
+   (MULTIPLE-VALUE-LIST (HURWITZ-ZETA 1.5d0 2.5d0)))
+  (LISP-UNIT::ASSERT-NUMERICAL-EQUAL
+   (LIST 0.7651470246254092d0 1.0661276646941275d-14)
+   (MULTIPLE-VALUE-LIST (ETA 1.5d0))))