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

Add polynomial; needs an array type in CFFI. GSL complex type.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@2990 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent c39c124b
No related branches found
No related tags found
No related merge requests found
;********************************************************
; 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)))
......@@ -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
......
......@@ -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)
......
;********************************************************
; 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))
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