diff --git a/README b/README new file mode 100644 index 0000000000000000000000000000000000000000..f121543fd4bb550dea713867f895c684c53acb5d --- /dev/null +++ b/README @@ -0,0 +1,19 @@ +Features: + +- Designed so that usage of GSL is as if GSL were written in CL: + information flow correctly through functions, full condition + handling, all information returned is available. + +- Conditions are handled in a Lisp-friendly way; full use of the CL + condition system with full explanations. + +- Uses CFFI. Should work with any Common Lisp that supports CFFI with + callbacks and foreign-funcall. + +- This library is pure common lisp. No C is used; the only requirement +is that the GSL library must be loaded and accessible. + +- ??? No actual functions are defined, just the macros for gaining +access to the functions. Defining the functions you need yourself is +very easy. + diff --git a/conditions.lisp b/conditions.lisp new file mode 100644 index 0000000000000000000000000000000000000000..11c95bd30a52a4cb8455d3c81592d4ece8cf2cb8 --- /dev/null +++ b/conditions.lisp @@ -0,0 +1,96 @@ +;******************************************************** +; file: conditions.lisp +; description: GSL errors +; date: Sat Mar 4 2006 - 18:33 +; author: Liam M. Healy +; modified: Wed Mar 8 2006 - 23:13 +;******************************************************** + +(in-package :gsl) + +(cffi:defcenum gsl-errorno + "Error codes for GSL, from /usr/include/gsl/gsl_errno.h." + (:CONTINUE -2) + :FAILURE :SUCCESS :EDOM :ERANGE :EFAULT :EINVAL :EFAILED :EFACTOR + :ESANITY :ENOMEM :EBADFUNC :ERUNAWAY :EMAXITER :EZERODIV :EBADTOL + :ETOL :EUNDRFLW :EOVRFLW :ELOSS :EROUND :EBADLEN :ENOTSQR :ESING + :EDIVERGE :EUNSUP :EUNIMPL :ECACHE :ETABLE :ENOPROG :ENOPROGJ + :ETOLF :ETOLX :ETOLG :EOF) + +(defparameter *gsl-error-alist* + '((-2 . "Iteration has not converged") + (-1 . "Failure") + (0 . "Success") + (1 . "Input domain error") + (2 . "Output range error") + (3 . "Invalid pointer") + (4 . "Invalid argument supplied by user") + (5 . "Generic failure") + (6 . "Factorization failed") + (7 . "Sanity check failed - shouldn't happen") + (8 . "Malloc failed") + (9 . "Problem with user-supplied function") + (10 . "Iterative process is out of control") + (11 . "Exceeded max number of iterations") + (12 . "Tried to divide by zero") + (13 . "User specified an invalid tolerance") + (14 . "Failed to reach the specified tolerance") + (15 . "Underflow") + (16 . "Overflow ") + (17 . "Loss of accuracy") + (18 . "Failed because of roundoff error") + (19 . "Matrix, vector lengths are not conformant") + (20 . "Matrix not square") + (21 . "Apparent singularity detected") + (22 . "Integral or series is divergent") + (23 . "Requested feature is not supported by the hardware") + (24 . "Requested feature not (yet) implemented") + (25 . "Cache limit exceeded") + (26 . "Table limit exceeded") + (27 . "Iteration is not making progress towards solution") + (28 . "Jacobian evaluations are not improving the solution") + (29 . "Cannot reach the specified tolerance in F") + (30 . "Cannot reach the specified tolerance in X") + (31 . "Cannot reach the specified tolerance in gradient") + (32 . "End of file"))) + +(define-condition gsl-error (arithmetic-error) + ((gsl-errno :initarg :gsl-errno :reader gsl-errno) + (gsl-reason :initarg :gsl-reason :reader gsl-reason) + (gsl-source-file :initarg :gsl-source-file :reader gsl-source-file) + (gsl-line-number :initarg :gsl-line-number :reader gsl-line-number)) + (:report + (lambda (condition stream) + (format stream "~a (~a), ~a in ~a at line ~d" + (rest (assoc (gsl-errno condition) *gsl-error-alist*)) + (cffi:foreign-enum-keyword 'gsl-errorno (gsl-errno condition)) + (gsl-reason condition) + (gsl-source-file condition) + (gsl-line-number condition)))) + (:documentation + "An error that has been signalled by the GNU Scientific Library.")) + +(define-condition gsl-warning (warning) + ((gsl-errno :initarg :gsl-errno :reader gsl-errno) + (gsl-context :initarg :gsl-context :reader gsl-context)) + (:report + (lambda (condition stream) + (format stream "GSL condition ~a (~d), ~a in ~a" + (cffi:foreign-enum-keyword 'gsl-errorno (gsl-errno condition)) + (gsl-errno condition) + (rest (assoc (gsl-errno condition) *gsl-error-alist*)) + (gsl-context condition)))) + (:documentation + "An warning that has been signalled by the GNU Scientific Library.")) + +(cffi:defcallback gsl-error :void + ((reason :string) (file :string) (line :int) (error-number :int)) + (error 'gsl-error + :gsl-errno error-number + :gsl-reason reason + :gsl-source-file file + :gsl-line-number line)) + +(cffi:foreign-funcall + "gsl_set_error_handler" + :pointer (cffi:callback gsl-error)) diff --git a/gsll.asd b/gsll.asd new file mode 100644 index 0000000000000000000000000000000000000000..0e261373627dff981cb26fa4ad203edaeaecd7b0 --- /dev/null +++ b/gsll.asd @@ -0,0 +1,24 @@ +;******************************************************** +; file: numerica.asd +; description: Definition of GSLL system +; date: +; author: Liam Healy +; modified: Thu Mar 9 2006 - 23:14 +;******************************************************** +;;; $Id: $ + +(asdf:defsystem "gsll" + :name "gsll" + :description "GNU Scientific Library for Lisp." + :version "0" + :author "Liam M. Healy" + :licence "GPL" + :depends-on (cffi) + :components + ((:file "init") + (:file "conditions" :depends-on (init)) + (:file "interface" :depends-on (init)) + (:file "mathematical" :depends-on (init)) + ;; complex numbers not necessary? Just make a struct. + (:file "special-functions" :depends-on ("interface")) + )) diff --git a/init.lisp b/init.lisp new file mode 100644 index 0000000000000000000000000000000000000000..62f67f61e926944a207b0daac481264f2d030120 --- /dev/null +++ b/init.lisp @@ -0,0 +1,23 @@ +;******************************************************** +; file: library.lisp +; description: Load GSL +; date: Sat Mar 4 2006 - 18:53 +; author: Liam M. Healy +; modified: Wed Mar 8 2006 - 22:27 +;******************************************************** + +(defpackage gsll + (:nicknames :gsl) + (:use :common-lisp :cffi)) + +(cffi:define-foreign-library libgslcblas + (:unix (:or "/usr/lib/libgslcblas.so.0" "/usr/lib/libgslcblas.so")) + (t (:default "libgslcblas"))) + +(cffi:use-foreign-library libgslcblas) + +(cffi:define-foreign-library libgsl + (:unix (:or "/usr/lib/libgsl.so.0" "/usr/lib/libgsl.so")) + (t (:default "libgsl"))) + +(cffi:use-foreign-library libgsl) diff --git a/interface.lisp b/interface.lisp new file mode 100644 index 0000000000000000000000000000000000000000..9a33ada5f291c0b748779afa9468aa071dc2dc2c --- /dev/null +++ b/interface.lisp @@ -0,0 +1,81 @@ +;******************************************************** +; file: interface.lisp +; description: Macros to interface GSL functions. +; date: Mon Mar 6 2006 - 22:35 +; author: Liam M. Healy +; modified: Sat Mar 11 2006 - 22:15 +;******************************************************** + +(in-package :gsl) + +;;; http://common-lisp.net/project/cffi/ + + +(eval-when (:compile-toplevel :load-toplevel :execute) + (export '())) + +(cffi:defcstruct sf-result + "Results from special functions with value and error estimate." + ;; file:///usr/share/doc/gsl-ref-html/gsl-ref_7.html#SEC61 + (val :double) + (err :double)) + +(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." + ;; file:///usr/share/doc/gsl-ref-html/gsl-ref_7.html#SEC61 + (val :double) + (err :double) + (e10 :int)) + +(cffi:defcenum sf-mode + "Numerical precision modes with which to calculate special functions." + ;; file:///usr/share/doc/gsl-ref-html/gsl-ref_7.html#SEC62 + :double-prec + :single-prec + :approx-prec) + +(defun pick-result (decl) + (case (second decl) + (sf-result + `((cffi:foreign-slot-value ,(first decl) 'sf-result 'val) + (cffi:foreign-slot-value ,(first decl) 'sf-result 'err))) + (:double `((cffi:mem-ref ,(first decl) :double))))) + +;;; Warning isn't quite right for lambdas. +;;; New name? +(defmacro defun-sf + (cl-name arguments gsl-name &key documentation return mode) + "Define a mathematical special function from GSL using the _e form + GSL function definition. If cl-name is :lambda, make a lambda." + (let ((args (mapcar #'first arguments)) + (return-symb-type + (mapcar (lambda (typ) (list (gensym "RET") typ)) return))) + `(,@(if (eq cl-name :lambda) + '(lambda) + `(defun ,cl-name)) + ,(if mode + `(,@args &optional (mode :double-prec)) + `(,@args)) + ,@(when documentation (list documentation)) + (cffi:with-foreign-objects + (,@(mapcar (lambda (d) `(,(first d) ',(second d))) + return-symb-type)) + (let ((status + (cffi:foreign-funcall + ,gsl-name + ,@(mapcan (lambda (ar) (list (second ar) (first ar))) + arguments) + ,@(when mode '(sf-mode mode)) + ,@(mapcan (lambda (r) `(:pointer ,(first r))) return-symb-type) + :int))) + (unless (eql :success (cffi:foreign-enum-keyword 'gsl-errorno status)) + (warn 'gsl-warning + :gsl-errno status :gsl-context `(,',cl-name ,,@args))) + (values + ,@(mapcan #'pick-result return-symb-type))))))) + + +;;; arguments: list like ((x :double) (y :double)) +;;; mode: t or nil + diff --git a/mathematical.lisp b/mathematical.lisp new file mode 100644 index 0000000000000000000000000000000000000000..46d7d19b972569748c938ab40527d8d6a7346708 --- /dev/null +++ b/mathematical.lisp @@ -0,0 +1,156 @@ +;******************************************************** +; file: mathematical.lisp +; description: Mathematical functions +; date: Wed Mar 8 2006 - 22:09 +; author: Liam M. Healy +; modified: Thu Mar 9 2006 - 23:10 +;******************************************************** + +(in-package :gsl) + +(eval-when (:compile-toplevel :load-toplevel :execute) + (export '(nanp infinityp finitep log+1 exp-1 hypotenuse approximately=))) + +;;;;**************************************************************************** +;;; Infinities and Not-a-number +;;;;**************************************************************************** + +(defmacro pmnil (x) + "+1, -1, or nil" + `(let ((v ,x)) + (when (or (= 1 v) (= -1 v)) + v))) + +(defun nanp (x) + "Return T if x is a double-float NaN." + (= 1 + (cffi:foreign-funcall + "gsl_isnan" + :double x + :int))) + +(defun infinityp (x) + "Return +1 if x is positive infinity, -1 if negative infinity + nil if finite." + (pmnil + (cffi:foreign-funcall + "gsl_isinf" + :double x + :int))) + +(defun finitep (x) + "Return T if finite." + (= 1 + (cffi:foreign-funcall + "gsl_finite" + :double x + :int))) + +;;;;**************************************************************************** +;;; Elementary functions +;;;;**************************************************************************** + +(defun log+1 (x) + "log(1+x), computed in a way that is accurate for small x." + (cffi:foreign-funcall + "gsl_log1p" + :double x + :double)) + +(defun exp-1 (x) + "exp(x)-1, computed in a way that is accurate for small x." + (cffi:foreign-funcall + "gsl_expm1" + :double x + :double)) + +(defun hypotenuse (x y) + "sqrt{x^2 + y^2} computed in a way that avoids overflow." + (cffi:foreign-funcall + "gsl_hypot" + :double x + :double y + :double)) + +;; Not clear why this function exists +(defun gsl-asinh (x) + "arcsinh" + (cffi:foreign-funcall + "gsl_asinh" + :double x + :double)) + +;;; Not clear why this function exists +(defun gsl-atanh (x) + "arcsinh" + (cffi:foreign-funcall + "gsl_atanh" + :double x + :double)) + +;;; gsl_ldexp +;;; gsl_frexp +;;; not mapped because CL has equivalents. + +;;;;**************************************************************************** +;;; Small integer powers +;;;;**************************************************************************** + +;;; Does CL need these? + +#| +A common complaint about the standard C library is its lack of a +function for calculating (small) integer powers. GSL provides a +simple functions to fill this gap. For reasons of efficiency, +these functions do not check for overflow or underflow +conditions. + +Function: double gsl_pow_int (double x, int n) + This routine computes the power x^n for integer n. The power is computed efficiently--for example, x^8 is computed as ((x^2)^2)^2, requiring only 3 multiplications. A version of this function which also computes the numerical error in the result is available as gsl_sf_pow_int_e. + +Function: double gsl_pow_2 (const double x) +Function: double gsl_pow_3 (const double x) +Function: double gsl_pow_4 (const double x) +Function: double gsl_pow_5 (const double x) +Function: double gsl_pow_6 (const double x) +Function: double gsl_pow_7 (const double x) +Function: double gsl_pow_8 (const double x) +Function: double gsl_pow_9 (const double x) + These functions can be used to compute small integer powers x^2, x^3, etc. efficiently. The functions will be inlined when possible so that use of these functions should be as efficient as explicitly writing the corresponding product expression. + +|# + +;;;;**************************************************************************** +;;; Approximate Comparison of Floating Point Numbers +;;;;**************************************************************************** + +;;; It is sometimes useful to be able to compare two floating point +;;; numbers approximately, to allow for rounding and truncation +;;; errors. This function implements the approximate +;;; floating-point comparison algorithm proposed by D.E. Knuth in +;;; Section 4.2.2 of Seminumerical Algorithms (3rd edition). + +(defun approximately= (x y epsilon) + "This function determines whether x and y are approximately equal + to a relative accuracy epsilon. + + The relative accuracy is measured using an interval of size 2 + \delta, where \delta = 2^k \epsilon and k is the maximum + base-2 exponent of x and y as computed by the function + frexp(). + + If x and y lie within this interval, they are considered + approximately equal and the function returns 0. Otherwise if + x < y, the function returns -1, or if x > y, the function + returns +1. + + The implementation is based on the package fcmp by + T.C. Belding." + (pmnil + (cffi:foreign-funcall + "gsl_fcmp" + :double x + :double y + :double epsilon + :double))) + diff --git a/special-functions.lisp b/special-functions.lisp new file mode 100644 index 0000000000000000000000000000000000000000..2b16eeb2efae7cae4744a44e44abe46991309621 --- /dev/null +++ b/special-functions.lisp @@ -0,0 +1,206 @@ +;******************************************************** +; file: special-functions.lisp +; description: GSL types +; date: Sat Mar 4 2006 - 21:07 +; author: Liam M. Healy +; modified: Sun Mar 12 2006 - 00:40 +;******************************************************** + +(in-package :gsl) + +(eval-when (:compile-toplevel :load-toplevel :execute) + (export + '(airy-ai airy-bi airy-ai-scaled airy-bi-scaled + airy-ai-deriv airy-ai-deriv-scaled airy-bi-deriv airy-bi-deriv-scaled + airy-zero-ai airy-zero-bi airy-zero-ai-deriv airy-zero-bi-deriv))) + +;;;;**************************************************************************** +;;;; Airy functions +;;;;**************************************************************************** + +(defun-sf airy-ai ((x :double)) + "gsl_sf_airy_Ai_e" + :return (sf-result) + :documentation "The Airy function Ai(x)." + :mode t) + +;;; (airy-ai 0.45d0) +;;; 0.24308135437540998d0 +;;; 6.288670879282072d-17 + +(defun-sf airy-bi ((x :double)) + "gsl_sf_airy_Bi_e" + :return (sf-result) + :documentation "The Airy function Bi(x)." + :mode t) + +(defun-sf airy-ai-scaled ((x :double)) + "gsl_sf_airy_Ai_scaled_e" + :return (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 @c{$\exp(+(2/3) x^{3/2})$} @math{\exp(+(2/3) x^(3/2))}, and is 1 for @math{x<0}." + :mode t) + +(defun-sf airy-bi-scaled ((x :double)) + "gsl_sf_airy_Bi_scaled_e" + :return (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 @c{$\exp(-(2/3) x^{3/2})$} @math{exp(-(2/3) x^(3/2))}, and is 1 for @math{x<0}." + :mode t) + +(defun-sf airy-ai-deriv ((x :double)) + "gsl_sf_airy_Ai_deriv_e" + :return (sf-result) + :documentation "The Airy function derivative Ai'(x)." + :mode t) + +(defun-sf airy-bi-deriv ((x :double)) + "gsl_sf_airy_Bi_deriv_e" + :return (sf-result) + :documentation "The Airy function derivative Bi'(x)." + :mode t) + +(defun-sf airy-ai-deriv-scaled ((x :double)) + "gsl_sf_airy_Ai_deriv_scaled_e" + :return (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 @c{$\exp(+(2/3) x^{3/2})$} @math{\exp(+(2/3) x^(3/2))}, and is 1 for @math{x<0}." + :mode t) + +(defun-sf airy-bi-deriv-scaled ((x :double)) + "gsl_sf_airy_Bi_deriv_scaled_e" + :return (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 @c{$\exp(-(2/3) x^{3/2})$} @math{exp(-(2/3) x^(3/2))}, and is 1 for @math{x<0}." + :mode t) + +(defun-sf airy-zero-ai ((s :unsigned-int)) + "gsl_sf_airy_zero_Ai_e" + :return (sf-result) + :documentation "The location of the @var{s}-th zero of the Airy function @math{Ai(x)}." + :mode t) + +(defun-sf airy-zero-bi ((s :unsigned-int)) + "gsl_sf_airy_zero_Bi_e" + :return (sf-result) + :documentation "The location of the @var{s}-th zero of the Airy function @math{Bi(x)}." + :mode t) + +(defun-sf airy-zero-ai-deriv ((s :unsigned-int)) + "gsl_sf_airy_zero_Ai_deriv_e" + :return (sf-result) + :documentation "The location of the @var{s}-th zero of the Airy +function derivative @math{Ai'(x)}." + :mode t) + +(defun-sf airy-zero-bi-deriv ((s :unsigned-int)) + "gsl_sf_airy_zero_Bi_deriv_e" + :return (sf-result) + :documentation "The location of the @var{s}-th zero of the Airy function derivative @math{Bi'(x)}." + :mode t) + +;;;;**************************************************************************** +;;;; Bessel functions +;;;;**************************************************************************** + +(eval-when (:compile-toplevel :load-toplevel :execute) + (export + '(bessel-J0 bessel-J1 bessel-Jn))) + +(defun-sf bessel-J0 ((x :double)) + "gsl_sf_bessel_J0_e" + :documentation + "The regular cylindrical Bessel function of zeroth order, @math{J_0(x)}." + :return (sf-result)) + +(defun-sf bessel-J1 ((x :double)) + "gsl_sf_bessel_J1_e" + :documentation + "The regular cylindrical Bessel function of first order, @math{J_1(x)}." + :return (sf-result)) + +(defun-sf bessel-Jn ((n :int) (x :double)) + "gsl_sf_bessel_Jn_e" + :documentation + "The regular cylindrical Bessel function of order @var{n}, @math{J_n(x)}." + :return (sf-result)) + +;;; Doesn't work; how do I get an array of doubles back? +(defun bessel-Jn-array (nmin nmax x) + "The values of the regular cylindrical Bessel functions @math{J_n(x)} for @math{n} from @var{nmin} to @var{nmax} inclusive. The values are computed using recurrence relations for efficiency, and therefore may differ slightly from the exact values." + (with-foreign-objects ((ptr :double)) + (let ((status + (foreign-funcall "gsl_sf_bessel_Jn_array" + :int + nmin + :int + nmax + :double + x + :pointer + ptr + :int))) + (unless + (eql :success + (foreign-enum-keyword 'gsl-errorno status)) + (warn 'gsl-warning + :gsl-errno + status + :gsl-context + `(foobar ,nmin ,nmax ,x))) + (loop for i from 0 below (- nmax nmin) + collect (cffi:mem-aref ptr :double i))))) + + +(defun-sf regular-spherical-bessel ((order :int) (x :double)) + "gsl_sf_bessel_jl_e" + :documentation "The regular spherical Bessel function j_l." + :return (sf-result)) + +;;; (regular-spherical-bessel 3 12.4d0) +;;; 0.07813007015855176d0 +;;; 1.9083196617232666d-16 + +;;;;**************************************************************************** +;;;; Dilogarithm +;;;;**************************************************************************** + +;;; dilog merge complex and real +(defun dilogarithm (x) + "The dilogarithm." + (etypecase x + (double-float + (funcall + (defun-sf :lambda ((x :double)) + "gsl_sf_dilog_e" + :return (sf-result)) + x)) + (complex + (multiple-value-bind (re re-err im im-err) + (funcall + ;; returns two gsl_sf_result + (defun-sf :lambda + ((radius :double) (angle :double)) + "gsl_sf_complex_dilog_e" + :return (sf-result sf-result)) + (abs x) + (phase x)) + (values + (complex re im) + (complex re-err im-err)))))) + +(defun-sf jacobian-elliptic-functions ; CL name + ((u :double) (m :double)) ; inputs + "gsl_sf_elljac_e" ; GSL name + :documentation + "The Jacobian elliptic functions sn(u|m), cn(u|m), dn(u|m) + computed by descending Landen transformations. + See Abramowitz & Stegun, Chapter 16" + :return (:double :double :double)) + +;;; > (jacobian-elliptic-functions 0.61802d0 0.5d0) +;;; 0.564575752943391 +;;; 0.8253812568676386 +;;; 0.916857191493965 +;;; > (jacobian-elliptic-functions 0.2d0 0.81d0) +;;; 0.19762082367187697 +;;; 0.9802785369736752 +;;; 0.9840560289645665 +;;; > (jacobian-elliptic-functions 0.61802d0 1.5d0) +;;; ;;;error