From 313da4b9ec7b90a6c71984c311da0a9aa177f668 Mon Sep 17 00:00:00 2001
From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30>
Date: Sat, 17 Jun 2006 02:57:50 +0000
Subject: [PATCH] Port special-functions/hypergeometric.lisp, laguerre.lisp,
 lambert.lisp, legendre.lisp, logarithm.lisp to new defun-gsl, add tests.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3095 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
---
 gsll.asd                              |  12 +-
 special-functions/hypergeometric.lisp | 114 +++++----
 special-functions/laguerre.lisp       |  27 +--
 special-functions/lambert.lisp        |  12 +-
 special-functions/legendre.lisp       | 333 ++++++++++++++------------
 special-functions/logarithm.lisp      |  57 +++--
 6 files changed, 309 insertions(+), 246 deletions(-)

diff --git a/gsll.asd b/gsll.asd
index f94c90e9..47a6cb1e 100644
--- a/gsll.asd
+++ b/gsll.asd
@@ -3,7 +3,7 @@
 ; description: Definition of GSLL system 
 ; date:        
 ; author:      Liam Healy
-; modified:    Tue Jun 13 2006 - 22:58
+; modified:    Fri Jun 16 2006 - 22:21
 ;********************************************************
 ;;; $Id: $
 
@@ -60,11 +60,11 @@
 	     (:file "fermi-dirac")
 	     (:file "gamma")
 	     (:file "gegenbauer")
-	     ;;(:file "hypergeometric")
-	     ;;(:file "laguerre")
-	     ;;(:file "lambert")
-	     ;;(:file "legendre")
-	     ;;(:file "logarithm")
+	     (:file "hypergeometric")
+	     (:file "laguerre")
+	     (:file "lambert")
+	     (:file "legendre")
+	     (:file "logarithm")
 	     ;;(:file "power")
 	     ;;(:file "psi")
 	     ;;(:file "synchrotron")
diff --git a/special-functions/hypergeometric.lisp b/special-functions/hypergeometric.lisp
index 69bd7fc0..95380bf5 100644
--- a/special-functions/hypergeometric.lisp
+++ b/special-functions/hypergeometric.lisp
@@ -3,58 +3,74 @@
 ; description: Hypergeometric function                   
 ; date:        Fri Apr 28 2006 - 23:00                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Apr 29 2006 - 18:35
+; modified:    Fri Jun 16 2006 - 21:05
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl hypergeometric-0F1 ((c :double) (x :double))
-  "gsl_sf_hyperg_0F1_e"
-  :return (sf-result)
+(defun-gsl hypergeometric-0F1 (c x)
+  "gsl_sf_hyperg_0F1_e" ((c :double) (x :double) (ret sf-result))
   :documentation "The hypergeometric function @math{0F1(c,x)}.")
 
-(defun-gsl hypergeometric-1F1-int ((m :int) (n :int) (x :double))
-  "gsl_sf_hyperg_1F1_int_e"
-  :return (sf-result)
+(defgeneric hypergeometric-1F1 (m n x)
+  (:documentation "The confluent hypergeometric function
+   @math{1F1(m,n,x) = M(m,n,x)}."))
+
+(defun-gsl hypergeometric-1F1 ((m fixnum) (n fixnum) 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}.")
 
-(defun-gsl hypergeometric-1F1 ((a :double) (b :double) (x :double))
-  "gsl_sf_hyperg_1F1_e"
-  :return (sf-result)
-  :documentation "Compute the confluent hypergeometric function
+(defun-gsl hypergeometric-1F1 ((a double-float) (b double-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}.")
 
-(defun-gsl hypergeometric-U-int ((m :int) (n :int) (x :double))
-  "gsl_sf_hyperg_U_int_e"
-  :return (sf-result)
+(defgeneric hypergeometric-U (m n x)
+  (:documentation "The confluent hypergeometric function
+   @math{U(m,n,x)}."))
+
+(defun-gsl hypergeometric-U ((m fixnum) (n fixnum) 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}.")
 
-(defun-gsl hypergeometric-U-int-e10 ((m :int) (n :int) (x :double))
+(defun-gsl hypergeometric-U ((a double-float) (b double-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."))
+
+(defun-gsl hypergeometric-U-e10 ((m fixnum) (n fixnum) x)
   "gsl_sf_hyperg_U_int_e10_e"
-  :return (sf-result-e10)
+  ((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.")
 
-(defun-gsl hypergeometric-U ((a :double) (b :double) (x :double))
-  "gsl_sf_hyperg_U_e"
-  :return (sf-result)
-  :documentation "The confluent hypergeometric function @math{U(a,b,x)}.")
-
-(defun-gsl hypergeometric-U-e10 ((a :double) (b :double) (x :double))
+(defun-gsl hypergeometric-U-e10 ((a double-float) (b double-float) x)
   "gsl_sf_hyperg_U_e10_e"
-  :return (sf-result-e10)
+  ((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.")
 
-(defun-gsl hypergeometric-2F1
-    ((a :double) (b :double) (c :double) (x :double))
+(defun-gsl hypergeometric-2F1 (a b c x)
   "gsl_sf_hyperg_2F1_e"
-  :return (sf-result)
+  ((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
@@ -62,34 +78,30 @@
   approximation converges too slowly.  This occurs in the region of
   @math{x=1}, @math{c - a - b = m} for integer m.")
 
-(defun-gsl hypergeometric-2F1-conj
-    (((realpart a) :double) ((imagpart a) :double) (c :double) (x :double))
+(defun-gsl hypergeometric-2F1-conj (a c x)
   "gsl_sf_hyperg_2F1_conj_e"
-  :function (a c x)
-  :return (sf-result)
+  (((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}.")
 
-(defun-gsl hypergeometric-renorm
-    ((a :double) (b :double) (c :double) (x :double))
+(defun-gsl hypergeometric-2F1-renorm (a b c x)
   "gsl_sf_hyperg_2F1_renorm_e"
-  :return (sf-result)
+  ((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}.")
 
-(defun-gsl hypergeometric-conj-renorm
-    (((realpart a) :double) ((imagpart a) :double) (c :double) (x :double))
+(defun-gsl hypergeometric-2F1-conj-renorm (a c x)
   "gsl_sf_hyperg_2F1_conj_renorm_e"
-  :function (a b c x)
-  :return (sf-result)
+  (((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}.")
 
-(defun-gsl hypergeometric-2F0
-    ((a :double) (b :double) (x :double))
+(defun-gsl hypergeometric-2F0 (a b x)
   "gsl_sf_hyperg_2F0_e"
-  :return (sf-result)
+  ((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
@@ -105,16 +117,36 @@
    (hypergeometric-0f1 0.5d0 1.0d0))
   (lisp-unit:assert-first-fp-equal
    "0.543656365692d+01"
-   (hypergeometric-1F1-int 2 1 1.0d0))
+   (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)))
diff --git a/special-functions/laguerre.lisp b/special-functions/laguerre.lisp
index 2261eff5..9557cbcb 100644
--- a/special-functions/laguerre.lisp
+++ b/special-functions/laguerre.lisp
@@ -3,33 +3,29 @@
 ; description: Laguerre polynomials                    
 ; date:        Fri Apr 28 2006 - 20:40                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Apr 29 2006 - 18:32
+; modified:    Fri Jun 16 2006 - 21:12
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
-(defun-gsl laguerre-1 ((a :double) (x :double))
-  "gsl_sf_laguerre_1_e"
-  :return (sf-result)
+(defun-gsl 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.")
 
-(defun-gsl laguerre-2 ((a :double) (x :double))
-  "gsl_sf_laguerre_2_e"
-  :return (sf-result)
+(defun-gsl 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.")
 
-(defun-gsl laguerre-3 ((a :double) (x :double))
-  "gsl_sf_laguerre_3_e"
-  :return (sf-result)
+(defun-gsl 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.")
 
-(defun-gsl laguerre ((n :int) (a :double) (x :double))
-  "gsl_sf_laguerre_n_e"
-  :return (sf-result)
+(defun-gsl 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}.")
 
@@ -46,4 +42,7 @@
    (laguerre-2 1.0d0 3.0d0))
   (lisp-unit:assert-first-fp-equal
    "-0.500000000000d+00"
-   (laguerre-3 1.0d0 3.0d0)))
+   (laguerre-3 1.0d0 3.0d0))
+  (lisp-unit:assert-first-fp-equal
+   "0.875000000000d+00"
+   (laguerre 4 1.0d0 3.0d0)))
diff --git a/special-functions/lambert.lisp b/special-functions/lambert.lisp
index 20498745..b5bbaed7 100644
--- a/special-functions/lambert.lisp
+++ b/special-functions/lambert.lisp
@@ -3,7 +3,7 @@
 ; description: Lambert's W functions
 ; date:        Fri Apr 28 2006 - 20:40                   
 ; author:      Liam M. Healy                             
-; modified:    Sat Apr 29 2006 - 19:00
+; modified:    Fri Jun 16 2006 - 21:13
 ;********************************************************
 ;;; $Id: $
 
@@ -18,15 +18,13 @@
 ;;; @math{W_@{-1@}(x)} to be the other real branch, where
 ;;; @math{W < -1} for @math{x < 0}.  
 
-(defun-gsl lambert-W0 ((x :double))
-  "gsl_sf_lambert_W0_e"
-  :return (sf-result)
+(defun-gsl 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)}.")
 
-(defun-gsl lambert-Wm1 ((x :double))
-  "gsl_sf_lambert_Wm1_e"
-  :return (sf-result)
+(defun-gsl lambert-Wm1 (x)
+  "gsl_sf_lambert_Wm1_e" ((x :double) (ret sf-result))
   :documentation
   "The secondary real-valued branch of the Lambert W function, 
    @math{W_@{-1@}(x)}.")
diff --git a/special-functions/legendre.lisp b/special-functions/legendre.lisp
index d29f655f..5a891d59 100644
--- a/special-functions/legendre.lisp
+++ b/special-functions/legendre.lisp
@@ -3,87 +3,73 @@
 ; description: Legendre functions                        
 ; date:        Sat Apr 29 2006 - 19:16                   
 ; author:      Liam M. Healy                             
-; modified:    Sun Apr 30 2006 - 12:25
+; modified:    Fri Jun 16 2006 - 22:56
 ;********************************************************
 ;;; $Id: $
 
 (in-package :gsl)
 
+;;; legendre-Plm-deriv-array same answer as legendre-Plm-array?
+
 ;;;;****************************************************************************
 ;;;; Legendre polynomials
 ;;;;****************************************************************************
 
-(defun-gsl legendre-P1 ((x :double))
-  "gsl_sf_legendre_P1_e"
+(defun-gsl legendre-P1 (x)
+  "gsl_sf_legendre_P1_e" ((x :double) (ret sf-result))
   :documentation
   "The Legendre polynomials @math{P_1(x)} using an explicit
-   representation."
-  :return (sf-result))
+   representation.")
 
-(defun-gsl legendre-P2 ((x :double))
-  "gsl_sf_legendre_P2_e"
+(defun-gsl legendre-P2 (x)
+  "gsl_sf_legendre_P2_e" ((x :double) (ret sf-result))
   :documentation
   "The Legendre polynomials @math{P_2(x)} using an explicit
-   representation."
-  :return (sf-result))
+   representation.")
 
-(defun-gsl legendre-P3 ((x :double))
-  "gsl_sf_legendre_P3_e"
+(defun-gsl legendre-P3 (x)
+  "gsl_sf_legendre_P3_e" ((x :double) (ret sf-result))
   :documentation
   "The Legendre polynomials @math{P_3(x)} using an explicit
-   representation."
-  :return (sf-result))
+   representation.")
 
-(defun-gsl legendre-Pl ((l :int) (x :double))
-  "gsl_sf_legendre_Pl_e"
+(defun-gsl legendre-Pl (l x)
+  "gsl_sf_legendre_Pl_e" ((l :int) (x :double) (ret sf-result))
   :documentation
   "The Legendre polynomial @math{P_l(x)} for a specific value of @var{l},
-   @var{x} subject to @math{l >= 0}, @math{|x| <= 1}."
-  :return (sf-result))
+   @var{x} subject to @math{l >= 0}, @math{|x| <= 1}.")
 
-(defun-gsl legendre-Pl-array
-    (((dim0 array) :int) (x :double) ((gsl-array array) :pointer))
+(defun-gsl legendre-Pl-array (x array)
   "gsl_sf_legendre_Pl_array"
+  (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :documentation "Compute an array of Legendre polynomials
-  @math{P_l(x)} for @math{l = 0, \dots, lmax}, @math{|x| <= 1}."
-  :function (x array)
-  :invalidate (array)
-  :return-input (array))
-
-;;; (defparameter leg (make-data 'vector nil 8))
-;;; (legendre-Pl-array 0.5d0 leg)
-;;; #<GSL-VECTOR #(1.0d0 0.5d0 -0.125d0 -0.4375d0 -0.2890625d0 0.08984375d0
-;;;                0.3232421875d0 0.22314453125d0) {BCC2319}>
-
-(defun-gsl legendre-Pl-deriv-array
-    (((dim0 array) :int) (x :double) ((gsl-array array) :pointer))
+  @math{P_l(x)} for @math{l = 0, \dots, length(array)}, @math{|x| <= 1}."
+  :invalidate (array))
+
+(defun-gsl legendre-Pl-deriv-array (x array)
   "gsl_sf_legendre_Pl_deriv_array"
+  (((1- (dim0 array)) :int) (x :double) ((gsl-array array) :pointer))
   :documentation "Compute an array of Legendre polynomials derivatives
-  @math{dP_l(x)/dx}, for @math{l = 0, \dots, lmax}, @math{|x| <= 1}."
-  :function (x array)
-  :invalidate (array)
-  :return-input (array))
+  @math{dP_l(x)/dx}, for @math{l = 0, \dots,  length(array)}, @math{|x| <= 1}."
+  :invalidate (array))
 
-(defun-gsl legendre-Q0 ((x :double))
-  "gsl_sf_legendre_Q0_e"
+(defun-gsl legendre-Q0 (x)
+  "gsl_sf_legendre_Q0_e" ((x :double) (ret sf-result))
   :documentation
   "The Legendre function @math{Q_0(x)} for @math{x > -1},
-   @math{x /= 1}."
-  :return (sf-result))
+   @math{x /= 1}.")
 
-(defun-gsl legendre-Q1 ((x :double))
-  "gsl_sf_legendre_Q1_e"
+(defun-gsl legendre-Q1 (x)
+  "gsl_sf_legendre_Q1_e" ((x :double) (ret sf-result))
   :documentation
   "The Legendre function @math{Q_1(x)} for @math{x > -1},
-   @math{x /= 1}."
-  :return (sf-result))
+   @math{x /= 1}.")
 
-(defun-gsl legendre-Ql ((l :int) (x :double))
-  "gsl_sf_legendre_Ql_e"
+(defun-gsl legendre-Ql (l x)
+  "gsl_sf_legendre_Ql_e" ((l :int) (x :double) (ret sf-result))
   :documentation
   "The Legendre function @math{Q_l(x)} for @math{x > -1},
-   @math{x /= 1}, @math{l >= 0}."
-  :return (sf-result))
+   @math{x /= 1}, @math{l >= 0}.")
 
 ;;;;****************************************************************************
 ;;;; Associated Legendre Polynomials and Spherical Harmonics
@@ -94,87 +80,70 @@
 ;;; @math{l} and can overflow for @math{l} larger than about 150.  There is
 ;;; no trouble for small @math{m}, but overflow occurs when @math{m} and
 ;;; @math{l} are both large.  Rather than allow overflows, these functions
-;;; refuse to calculate @math{P_l^m(x)} and return @code{GSL_EOVRFLW} when
+;;; refuse to calculate @math{P_l^m(x)} and return :EOVRFLW when
 ;;; they can sense that @math{l} and @math{m} are too big.
 
 ;;; If you want to calculate a spherical harmonic, then @emph{do not} use
-;;; these functions.  Instead use @code{gsl_sf_legendre_sphPlm()} below,
+;;; these functions.  Instead use legendre-sphPlm below,
 ;;; which uses a similar recursion, but with the normalized functions.
 
-(defun-gsl legendre-Plm ((l :int) (m :int) (x :double))
-  "gsl_sf_legendre_Plm_e"
+(defun-gsl legendre-Plm (l m x)
+  "gsl_sf_legendre_Plm_e" ((l :int) (m :int) (x :double) (ret sf-result))
   :documentation
   "The associated Legendre polynomial
-   @math{P_l^m(x)} for @math{m >= 0}, @math{l >= m}, @math{|x| <= 1}."
-  :return (sf-result))
-
+   @math{P_l^m(x)} for @math{m >= 0}, @math{l >= m}, @math{|x| <= 1}.")
 
-(defun-gsl legendre-Plm-array
-    (((+ (dim0 array) m -1) :int) (m :int) (x :double)
-     ((gsl-array array) :pointer))
+(defun-gsl legendre-Plm-array (m x array)
   "gsl_sf_legendre_Plm_array"
+  (((+ (dim0 array) m -1) :int) (m :int) (x :double)
+   ((gsl-array array) :pointer))
   :documentation "An array of Legendre polynomials
     @math{P_l^m(x)}, for @math{m >= 0}, 
-    @math{l = |m|, ..., lmax}, where lmax is the length
-    of the vector array, and @math{|x| <= 1}."
-  :function (m x array)
-  :invalidate (array)
-  :return-input (array))
-
-;;; (defparameter aleg (make-data 'vector nil 3))
-;;; (legendre-plm-array 3 0.5d0 aleg)
-;;; #<GSL-VECTOR #(-9.742785792574935d0 -34.09975027401227d0 -42.62468784251534d0) {BD445E9}>
-
-(defun-gsl legendre-Plm-deriv-array
-    (((+ (dim0 array) m -1) :int) (m :int) (x :double)
-     ((gsl-array array) :pointer))
+    @math{l = |m|, ..., |m|+length(array)-1} and @math{|x| <= 1}."
+  :invalidate (array))
+
+(defun-gsl legendre-Plm-deriv-array (m x array)
   "gsl_sf_legendre_Plm_deriv_array"
+  (((+ (dim0 array) m -1) :int) (m :int) (x :double)
+   ((gsl-array array) :pointer))
   :documentation "An array of Legendre polynomials derivatives
     @math{dP_l^m(x)/dx} for @math{m >= 0}, 
-    @math{l = |m|, ..., lmax}, where lmax is the length
-    of the vector array, and @math{|x| <= 1}."
-  :function (m x array)
-  :invalidate (array)
-  :return-input (array))
-
-(defun-gsl legendre-sphPlm ((l :int) (m :int) (x :double))
-  "gsl_sf_legendre_sphPlm_e"
+    @math{l = |m|, ..., length(array)} and @math{|x| <= 1}."
+  :invalidate (array))
+
+(defun-gsl legendre-sphPlm (l m x)
+  "gsl_sf_legendre_sphPlm_e" ((l :int) (m :int) (x :double) (ret sf-result))
   :documentation
   "The normalized associated Legendre polynomial
    @math{$\sqrt@{(2l+1)/(4\pi)@} \sqrt@{(l-m)!/(l+m)!@} P_l^m(x)$} suitable
    for use in spherical harmonics.  The parameters must satisfy
    @math{m >= 0}, @math{l >= m}, @math{|x| <= 1}.
    These routines avoid the overflows that occur for the standard
-   normalization of @math{P_l^m(x)}."
-  :return (sf-result))
+   normalization of @math{P_l^m(x)}.")
 
-(defun-gsl legendre-sphPlm-array
-    (((+ (dim0 array) m -1) :int) (m :int) (x :double)
-     ((gsl-array array) :pointer))
+(defun-gsl legendre-sphPlm-array (m x array)
   "gsl_sf_legendre_sphPlm_array"
+  (((+ (dim0 array) m -1) :int) (m :int) (x :double)
+   ((gsl-array array) :pointer))
   :documentation "An array of normalized associated Legendre functions
    @math{$\sqrt@{(2l+1)/(4\pi)@} \sqrt@{(l-m)!/(l+m)!@} P_l^m(x)$},
-   for @math{m >= 0}, @math{l = |m|, ..., lmax}, @math{|x| <= 1.0}."
-  :function (m x array)
-  :invalidate (array)
-  :return-input (array))
-
-(defun-gsl legendre-sphPlm-deriv-array
-    (((+ (dim0 array) m -1) :int) (m :int) (x :double)
-     ((gsl-array array) :pointer))
+   for @math{m >= 0}, @math{l = |m|, ..., length(array)}, @math{|x| <= 1.0}."
+  :invalidate (array))
+
+(defun-gsl legendre-sphPlm-deriv-array (m x array)
   "gsl_sf_legendre_sphPlm_deriv_array"
+  (((+ (dim0 array) m -1) :int) (m :int) (x :double)
+   ((gsl-array array) :pointer))
   :documentation "An array of normalized associated Legendre functions
-   derivatives for @math{m >= 0}, @math{l = |m|, ..., lmax}, @math{|x| <= 1.0}."
-  :function (m x array)
-  :invalidate (array)
-  :return-input (array))
+   derivatives for @math{m >= 0}, @math{l = |m|, ..., length(array)},
+   @math{|x| <= 1.0}."
+  :invalidate (array))
 
-(defun-gsl legendre-array-size ((lmax :int) (m :int))
-  "gsl_sf_legendre_array_size"
+(defun-gsl legendre-array-size (lmax m)
+  "gsl_sf_legendre_array_size" ((lmax :int) (m :int))
   :documentation "The size of @var{result_array}[] needed for the array
     versions of @math{P_l^m(x)}, @math{@var{lmax} - @var{m} + 1}."
-  :c-return-value :return
-  :return (:int))
+  :c-return :int)
 
 ;;;;****************************************************************************
 ;;;; Conical Functions
@@ -184,45 +153,39 @@
 ;;; @math{Q^\mu_@{-(1/2)+i\lambda@}} 
 ;;; are described in Abramowitz & Stegun, Section 8.12.
 
-(defun-gsl legendre-conicalP-half ((lambda :double) (x :double))
-  "gsl_sf_conicalP_half_e"
+(defun-gsl legendre-conicalP-half (lambda x)
+  "gsl_sf_conicalP_half_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The irregular Spherical Conical Function
-   @math{P^@{1/2@}_@{-1/2 + i \lambda@}(x)} for @math{x > -1}."
-  :return (sf-result))
+   @math{P^@{1/2@}_@{-1/2 + i \lambda@}(x)} for @math{x > -1}.")
 
-(defun-gsl legendre-conicalP-mhalf ((lambda :double) (x :double))
-  "gsl_sf_conicalP_mhalf_e"
+(defun-gsl legendre-conicalP-mhalf (lambda x)
+  "gsl_sf_conicalP_mhalf_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The regular Spherical Conical Function
-  @math{P^@{-1/2@}_@{-1/2 + i \lambda@}(x)} for @math{x > -1}."
-  :return (sf-result))
+  @math{P^@{-1/2@}_@{-1/2 + i \lambda@}(x)} for @math{x > -1}.")
 
-(defun-gsl legendre-conicalP-0 ((lambda :double) (x :double))
-  "gsl_sf_conicalP_0_e"
+(defun-gsl legendre-conicalP-0 (lambda x)
+  "gsl_sf_conicalP_0_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The conical function @math{P^0_@{-1/2 + i \lambda@}(x)}
-  for @math{x > -1}."
-  :return (sf-result))
+  for @math{x > -1}.")
 
-(defun-gsl legendre-conicalP-1 ((lambda :double) (x :double))
-  "gsl_sf_conicalP_1_e"
+(defun-gsl legendre-conicalP-1 (lambda x)
+  "gsl_sf_conicalP_1_e" ((lambda :double) (x :double) (ret sf-result))
   :documentation "The conical function 
-   @math{P^1_@{-1/2 + i \lambda@}(x)} for @math{x > -1}."
-  :return (sf-result))
+   @math{P^1_@{-1/2 + i \lambda@}(x)} for @math{x > -1}.")
 
-(defun-gsl legendre-regular-spherical-conical
-    ((l :int) (lambda :double) (x :double))
+(defun-gsl legendre-regular-spherical-conical (l lambda x)
   "gsl_sf_conicalP_sph_reg_e"
+  ((l :int) (lambda :double) (x :double) (ret sf-result))
   :documentation "The Regular Spherical Conical Function
    @math{P^@{-1/2-l@}_@{-1/2 + i \lambda@}(x)} for @math{x > -1},
-   @math{l >= -1}."
-  :return (sf-result))
+   @math{l >= -1}.")
 
-(defun-gsl legendre-regular-cylindrical-conical
-    ((l :int) (lambda :double) (x :double))
+(defun-gsl legendre-regular-cylindrical-conical (l lambda x)
   "gsl_sf_conicalP_cyl_reg_e"
+  ((l :int) (lambda :double) (x :double) (ret sf-result))
   :documentation "The Regular Cylindrical Conical Function
    @math{P^@{-m@}_@{-1/2 + i \lambda@}(x)} for @math{x > -1},
-   @math{m >= -1}."
-  :return (sf-result))
+   @math{m >= -1}.")
 
 ;;;;****************************************************************************
 ;;;; Radial Functions for Hyperbolic Space
@@ -234,43 +197,41 @@
 ;;; the flat limit, @math{\lambda \to \infty}, @math{\eta \to 0},
 ;;; @math{\lambda\eta} fixed.
 
-(defun-gsl legendre-H3d-0 ((lambda :double) (eta :double))
+(defun-gsl legendre-H3d-0 (lambda eta)
   "gsl_sf_legendre_H3d_0_e"
+  ((lambda :double) (eta :double) (ret sf-result))
   :documentation "The zeroth radial eigenfunction of the Laplacian on the
    3-dimensional hyperbolic space,
    @math{L^@{H3d@}_0(\lambda,\eta) := \sin(\lambda\eta)/(\lambda\sinh(\eta))}
    for @math{\eta >= 0}. In the flat limit this takes the form
-   @math{L^@{H3d@}_0(\lambda,\eta) = j_0(\lambda\eta)}."
-  :return (sf-result))
+   @math{L^@{H3d@}_0(\lambda,\eta) = j_0(\lambda\eta)}.")
 
-(defun-gsl legendre-H3d-1 ((lambda :double) (eta :double))
+(defun-gsl legendre-H3d-1 (lambda eta)
   "gsl_sf_legendre_H3d_1_e"
+  ((lambda :double) (eta :double) (ret sf-result))
   :documentation "The first radial eigenfunction of the Laplacian on
    the 3-dimensional hyperbolic space,
    @math{L^@{H3d@}_1(\lambda,\eta) := 1/\sqrt@{\lambda^2 + 1@}
    \sin(\lambda \eta)/(\lambda \sinh(\eta)) (\coth(\eta) - \lambda \cot(\lambda\eta))}
    for @math{\eta >= 0}.  In the flat limit this takes the form 
-   @math{L^@{H3d@}_1(\lambda,\eta) = j_1(\lambda\eta)}."
-  :return (sf-result))
+   @math{L^@{H3d@}_1(\lambda,\eta) = j_1(\lambda\eta)}.")
 
-(defun-gsl legendre-H3d ((l :int) (lambda :double) (eta :double))
+(defun-gsl legendre-H3d (l lambda eta)
   "gsl_sf_legendre_H3d_e"
+  ((l :int) (lambda :double) (eta :double) (ret sf-result))
   :documentation "The @var{l}-th radial eigenfunction of the
    Laplacian on the 3-dimensional hyperbolic space
    @math{\eta >= 0}, @c{$l \ge 0$} @math{l >= 0}.
    In the flat limit this takes the form
-   @math{L^@{H3d@}_l(\lambda,\eta) = j_l(\lambda\eta)}."
-  :return (sf-result))
+   @math{L^@{H3d@}_l(\lambda,\eta) = j_l(\lambda\eta)}.")
 
-(defun-gsl legendre-H3d-array
-    (((dim0 array) :int) (lambda :double) (eta :double)
-     ((gsl-array array) :pointer))
+(defun-gsl legendre-H3d-array (lambda eta array)
   "gsl_sf_legendre_H3d_array"
+  (((1- (dim0 array)) :int) (lambda :double) (eta :double)
+   ((gsl-array array) :pointer))
   :documentation "An array of radial eigenfunctions
-   @math{L^@{H3d@}_l(\lambda, \eta)} for @math{0 <= l <= lmax}."
-  :function (lambda eta array)
-  :invalidate (array)
-  :return-input (array))
+   @math{L^@{H3d@}_l(\lambda, \eta)} for @math{0 <= l <= length(array)}."
+  :invalidate (array))
 
 ;;; (defparameter hleg (make-data 'vector nil 3))
 ;;; (legendre-H3d-array 1.0d0 0.5d0 hleg)
@@ -281,47 +242,103 @@
 ;;;;****************************************************************************
 
 (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
+    (with-data (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))
+   (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
+    (with-data (arr vector-double 4)
+      (legendre-Plm-array 2 0.5d0 arr)
+      (data arr))))
+  (lisp-unit:assert-equal
+   '("0.225000000000d+01" "0.562500000000d+01"
+     "0.421875000000d+01" "-0.492187500000d+01")
+   ;; suspicious? same answer as legendre-Plm-array?
+   (lisp-unit:fp-sequence
+    (with-data (arr vector-double 4)
+      (legendre-Plm-deriv-array 2 0.5d0 arr)
+      (data arr))))
   (lisp-unit:assert-first-fp-equal
    "0.398506257222d-13"
    (legendre-sphplm 1200 1100 0.5d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+  (lisp-unit:assert-equal
+   '("0.248924639500d+00" "0.412794815148d+00"
+     "0.351206555622d+00" "0.515993518936d-01")
+   (lisp-unit:fp-sequence
+    (with-data (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.248924639500d+00" "0.412794815148d+00"
+     "0.351206555622d+00" "0.515993518936d-01")
+   (lisp-unit:fp-sequence
+    (with-data (arr vector-double 4)
+      (legendre-sphPlm-deriv-array 4 0.5d0 arr)
+      (data arr))))
+  (lisp-unit:assert-first-fp-equal
    "-0.125529904888d+00"
-   (LEGENDRE-CONICALP-HALF 3.5d0 10.0d0))
-  (LISP-UNIT:ASSERT-FIRST-FP-EQUAL
+   (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
+   (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
+   (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
+   (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
+   (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
+   (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
+   (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
+   (legendre-h3d-1 1.0d0 0.5d0))
+  (lisp-unit:assert-first-fp-equal
    "0.240061623900d-02"
-   (LEGENDRE-H3D 4 1.0d0 0.5d0)))
+   (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
+    (with-data (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 aaeeb637..a2eb2918 100644
--- a/special-functions/logarithm.lisp
+++ b/special-functions/logarithm.lisp
@@ -3,7 +3,7 @@
 ; description: Logarithm                                 
 ; date:        Sun Apr 30 2006 - 22:08                   
 ; author:      Liam M. Healy                             
-; modified:    Sun Apr 30 2006 - 22:43
+; modified:    Fri Jun 16 2006 - 22:20
 ;********************************************************
 ;;; $Id: $
 
@@ -13,36 +13,53 @@
   (:documentation
    "The natural logarithm of @var{x}, @math{\log(x)}, for @math{x > 0}."))
 
-(defun-gsl gsl-log ((x :double))
-  "gsl_sf_log_e"
-  :method ((x double-float))
-  :return (sf-result))
+(defun-gsl gsl-log ((x double-float))
+  "gsl_sf_log_e" ((x :double) (ret sf-result))
+  :type :method
+  :export t)
 
-(defun-gsl gsl-log (((realpart x) :double) ((imagpart x) :double))
+(defun-gsl gsl-log ((x complex))
   "gsl_sf_complex_log_e"
-  :method ((x complex))
+  (((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 (sf-result sf-result))
+  :return
+  ((complex (val re-ret) (val im-ret)) (complex (err re-ret) (err im-ret))))
 
-(defun-gsl log-abs ((x :double))
-  "gsl_sf_log_abs_e"
+(defun-gsl 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}."
-  :return (sf-result))
+  @math{\log(|x|)}, for @math{x \ne 0}.")
 
-(defun-gsl log-1+x ((x :double))
-  "gsl_sf_log_1plusx_e"
+(defun-gsl 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}."
-  :return (sf-result))
+   algorithm that is accurate for small @math{x}.")
 
-(defun-gsl log-1+x-m1 ((x :double))
-  "gsl_sf_log_1plusx_mx_e"
+(defun-gsl 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}."
-  :return (sf-result))
+  algorithm that is accurate for small @math{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)))
-- 
GitLab