diff --git a/cffi-array.lisp b/cffi-array.lisp new file mode 100644 index 0000000000000000000000000000000000000000..7bc96fd616a602aa5a732529eec1c7eed247ac7d --- /dev/null +++ b/cffi-array.lisp @@ -0,0 +1,81 @@ +;******************************************************** +; file: cffi-array.lisp +; description: Addition to CFFI +; date: Tue Mar 21 2006 - 18:57 +; author: Liam M. Healy +; modified: Tue Mar 21 2006 - 19:02 +;******************************************************** +;;; $Id: $ + +;;; http://article.gmane.org/gmane.lisp.cffi.devel/283 + +(in-package :cffi) + +(defun indexes-to-row-major-index (dimensions &rest subscripts) + (apply #'+ (maplist #'(lambda (x y) + (* (car x) (apply #'* (cdr y)))) + subscripts + dimensions))) + +(defun row-major-index-to-indexes (index dimensions) + (loop with idx = index + with rank = (length dimensions) + with indexes = (make-list rank) + for dim-index from (- rank 1) downto 0 do + (setf (values idx (nth dim-index indexes)) (floor idx (nth dim-index +dimensions))) + finally (return indexes))) + +(defun lisp-array-to-foreign (array ptr cffi-type dimensions) + "Copy elements from a Lisp array to ptr." + (loop with foreign-type-size = (foreign-type-size cffi-type) + with size = (reduce #'* dimensions) + for i from 0 below size + for offset = (* i foreign-type-size) + for element = (apply #'aref array + (row-major-index-to-indexes i dimensions)) + do (setf (mem-ref ptr cffi-type offset) element))) + +(defun foreign-array-to-lisp (ptr cffi-type dimensions) + "Copy elements from ptr into a Lisp array. +If ptr is a null pointer, returns nil." + (unless (null-pointer-p ptr) + (let ((array (make-array dimensions))) + (loop with foreign-type-size = (foreign-type-size cffi-type) + with size = (reduce #'* dimensions) + for i from 0 below size + for offset = (* i foreign-type-size) + for element = (mem-ref ptr cffi-type offset) + do (setf (apply #'aref array (row-major-index-to-indexes i +dimensions)) element)) + array))) + +(defun foreign-array-alloc (array cffi-type dimensions) + "Allocate a foreign array containing the elements of lisp array. +The foreign array must be freed with foreign-array-free." + (check-type array array) + (let ((ptr (foreign-alloc cffi-type :count (reduce #'* dimensions)))) + (lisp-array-to-foreign array ptr cffi-type dimensions) + ptr)) + +(defun foreign-array-free (ptr) + "Free a foreign array allocated by foreign-array-alloc." + (foreign-free ptr)) + +(defmacro with-foreign-array + ((var lisp-array cffi-type dimensions) &body body) + "Bind var to a foreign array containing lisp-array elements in body." + `(with-foreign-pointer (,var (* (reduce #'* ,dimensions) + (foreign-type-size ,cffi-type))) + (lisp-array-to-foreign ,lisp-array ,var ,cffi-type ,dimensions) + ,@body)) + +(defun foreign-aref (ptr cffi-type dimensions &rest indexes) + (let ((offset (* (foreign-type-size cffi-type) + (apply #'indexes-to-row-major-index dimensions indexes)))) + (mem-ref ptr cffi-type offset))) + +(defun (setf foreign-aref) (value ptr cffi-type dimensions &rest indexes) + (let ((offset (* (foreign-type-size cffi-type) + (apply #'indexes-to-row-major-index dimensions indexes)))) + (setf (mem-ref ptr cffi-type offset) value))) diff --git a/gsll.asd b/gsll.asd index f42a2ebdf652d9b1d972e8b3ae16a4e9b2c8e4a2..e423e258684d411f1abe1ea2b65ffa2dba4e29f8 100644 --- a/gsll.asd +++ b/gsll.asd @@ -3,7 +3,7 @@ ; description: Definition of GSLL system ; date: ; author: Liam Healy -; modified: Tue Mar 21 2006 - 17:49 +; modified: Tue Mar 21 2006 - 23:10 ;******************************************************** ;;; $Id: $ @@ -23,6 +23,8 @@ ((:file "conditions") (:file "mathematical"))) ;; complex numbers not necessary? Just make a struct. + (:file "cffi-array") + (:file "polynomial" :depends-on (init interface cffi-array)) ; see file (:module special-functions :depends-on (init interface) :components diff --git a/interface.lisp b/interface.lisp index ceac491e873a9fd7fe842c561c01bd8402999e52..80c2e6805340d3ce69fff03e6c51d4104c95315d 100644 --- a/interface.lisp +++ b/interface.lisp @@ -3,7 +3,7 @@ ; description: Macros to interface GSL functions. ; date: Mon Mar 6 2006 - 22:35 ; author: Liam M. Healy -; modified: Tue Mar 21 2006 - 17:11 +; modified: Tue Mar 21 2006 - 23:10 ;******************************************************** (in-package :gsl) @@ -28,6 +28,10 @@ and a scaling exponent e10, such that the value is val*10^e10." (err :double) (e10 :int)) +(cffi:defcstruct gsl-complex + "A complex number in GSL." + (dat :double :count 2)) + (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 @@ -74,7 +78,7 @@ and a scaling exponent e10, such that the value is val*10^e10." (defun pick-result (decl) (if (listp (second decl)) `((loop for i from 0 below ,(second (second decl)) - collect (cffi:mem-aref ,(first decl) ,(first (second decl)) i))) + collect (cffi:mem-aref ,(first decl) ,(first (second decl)) i))) (case (second decl) (sf-result `((cffi:foreign-slot-value ,(first decl) 'sf-result 'val) @@ -83,6 +87,12 @@ and a scaling exponent e10, such that the value is val*10^e10." `((cffi:foreign-slot-value ,(first decl) 'sf-result-e10 'val) (cffi:foreign-slot-value ,(first decl) 'sf-result-e10 'e10) (cffi:foreign-slot-value ,(first decl) 'sf-result-e10 'err))) + (gsl-complex + `((let ((carr + (cffi:foreign-slot-value ,(first decl) 'gsl-complex 'dat))) + (complex + (cffi:mem-aref carr :double 0) + (cffi:mem-aref carr :double 1))))) (:double `((cffi:mem-ref ,(first decl) :double)))))) (defun wfo-declare (d) diff --git a/polynomial.lisp b/polynomial.lisp new file mode 100644 index 0000000000000000000000000000000000000000..91c2f907cdebb737d2a604d4ec11453e60d7430d --- /dev/null +++ b/polynomial.lisp @@ -0,0 +1,176 @@ +;******************************************************** +; file: polynomial.lisp +; description: Polynomials +; date: Tue Mar 21 2006 - 18:33 +; author: Liam M. Healy +; modified: Tue Mar 21 2006 - 23:07 +;******************************************************** +;;; $Id: $ + +(in-package :gsl) + +;;; To do: +;;; Awaits incorporation of arrays into CFFI. +;;; Divided difference needs to be checked. +;;; Quadratic solver always gives output range warning, except when it needs to. +;;; Cubic solver always gives output range warning. + +;;; finish + +;;;;**************************************************************************** +;;;; Polynomial Evaluation +;;;;**************************************************************************** + +(export 'polynomial-eval) +(defun polynomial-eval (coefficients x) + "Evaluate the polyonomial with coefficients at the point x." + (let ((len (length coefficients))) + (cffi::with-foreign-array (coef coefficients :double (list len)) + (cffi:foreign-funcall + "gsl_poly_eval" + :pointer coef + :int len + :double x + :double)))) + +;;;;**************************************************************************** +;;;; Divided Difference Representation of Polynomials +;;;;**************************************************************************** + +;;; Use with-divided-difference to compute the divided difference, an +;;; opaque object that is then passed to +;;; get-divided-difference, eval-divided-difference +;;; in the body. + +(export + '(with-divided-difference get-divided-difference + eval-divided-difference taylor-divided-difference)) + +(defmacro with-divided-difference ((dd xa ya) &body body) + "Compute the divided difference and bind dd to it. + This variable may then be passed to divided-difference + functions in the body." + `(let ((,dd (find-divided-difference ,xa ,ya))) + (unwind-protect + (progn ,@body) + (cffi::foreign-array-free (first ,dd)) + (cffi::foreign-array-free (second ,dd))))) + +;;; Do not call this function directly; use with-divided-difference. +(defun find-divided-difference (xa ya) + "Find the divided difference representation of the interpolating polynomial + for the points stored in the arrays @var{xa} and @var{ya}." + (let* ((size (length xa)) + (dd (foreign-alloc :double :count size)) + (xac (foreign-alloc :double :count size))) + (cffi::with-foreign-array (yac ya :double (list size)) + ;; return code not checked + (cffi:foreign-funcall + "gsl_poly_dd_init" + :pointer dd + :pointer xac + :pointer yac + :uint size ; "size_t is almost always an unsigned int" + :int)) + (list dd xac size))) + +(defun get-divided-difference (dd) + "Convert the divided difference into an array." + (cffi::foreign-array-to-lisp + (first dd) :double (list (third dd)))) + +(defun eval-divided-difference (dd x) + "Evaluate the polynomial stored in divided-difference form + at the point @var{x}." + (cffi:foreign-funcall + "gsl_poly_dd_eval" + :pointer (first dd) + :pointer (second dd) + :uint (third dd) ; "size_t is almost always an unsigned int" + :double x + :double)) + +(defun taylor-divided-difference (dd xp) + "Convert the divided-difference representation of a polynomial + to a Taylor expansion about the point xp." + (let ((cc (foreign-alloc :double :count (third dd))) + (workspace (foreign-alloc :double :count (third dd)))) + (unwind-protect + (progn + ;; Return value not checked. + (cffi:foreign-funcall + "gsl_poly_dd_taylor" + :pointer cc + :double xp + :pointer (first dd) + :pointer (second dd) + :uint (third dd) ; "size_t is almost always an unsigned int" + :pointer workspace + :int) + (cffi::foreign-array-to-lisp + cc :double (list (third dd)))) + (cffi::foreign-array-free workspace) + (cffi::foreign-array-free cc)))) + +;;;;**************************************************************************** +;;;; Quadratic Equations +;;;;**************************************************************************** + +;;; GSL seems to always give a output range error warning: +;;; WARNING: +;;; GSL condition ERANGE (2), Output range error in +;;; (SOLVE-QUADRATIC 1.0d0 2.0d0 1.0d0) +;;; Oddly, there is NO warning if a complex answer is generated +;;; from the real root solver. + +(defun-sf solve-quadratic ((a :double) (b :double) (c :double)) + "gsl_poly_solve_quadratic" + :documentation + "The real roots of the quadratic equation a x^2 + b x + c = 0. + Warning: if roots are not real, nonsense is produced and + no warning is given." + :return (:double :double)) + +(defun-sf solve-quadratic-complex ((a :double) (b :double) (c :double)) + "gsl_poly_complex_solve_quadratic" + :documentation + "The complex roots of the quadratic equation a x^2 + b x + c = 0." + :return (gsl-complex gsl-complex)) + +;;;;**************************************************************************** +;;;; Cubic Equations +;;;;**************************************************************************** + +;;; GSL seems to always give warning: +;;; WARNING: +;;; GSL condition EFAULT (3), Invalid pointer in +;;; (SOLVE-CUBIC -3.0d0 3.0d0 -1.0d0) + +;;; GSL appears to give a different return code (EDOM) if there is +;;; only one real root. + +;;; (solve-cubic -6.0d0 -13.0d0 42.0d0) +;;; -3.0d0 +;;; 1.9999999999999996d0 +;;; 7.0d0 + +;;; (solve-cubic -1.0d0 1.0d0 -1.0d0) +;;; 1.0d0 +;;; 1.9999999999999996d0 +;;; 7.0d0 + +(defun-sf solve-cubic ((a :double) (b :double) (c :double)) + "gsl_poly_solve_cubic" + :documentation + "Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0 + with a leading coefficient of unity. The real roots are returned. + No special indication is given if there is only one real root." + :return (:double :double :double)) + +(defun-sf solve-cubic-complex ((a :double) (b :double) (c :double)) + "gsl_poly_complex_solve_cubic" + :documentation + "Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0 + with a leading coefficient of unity. The real roots are returned. + No special indication is given if there is only one real root." + :return (gsl-complex gsl-complex gsl-complex))