Skip to content
Snippets Groups Projects
Commit ac31280d authored by liam's avatar liam
Browse files

Port more special functions to new defun-gsl. Make polynomial

workable if not completely done (divided differences compiles but not
tested).


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