Commit 3a7fc273 authored by liam's avatar liam

GSLL initial macro, conditions, and some definitions.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@2975 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 646c623c
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.
;********************************************************
; 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))
;********************************************************
; 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"))
))
;********************************************************
; 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)
;********************************************************
; 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
;********************************************************
; 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)))
;********************************************************
; 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
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment