Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv30704
Modified Files: poly.lisp Log Message: Replace some functions with macros that clean up after themselves. Plug a leak.
Date: Mon Apr 4 02:46:44 2005 Author: edenny
Index: cl-gsl/poly.lisp diff -u cl-gsl/poly.lisp:1.2 cl-gsl/poly.lisp:1.3 --- cl-gsl/poly.lisp:1.2 Sun Mar 13 01:48:25 2005 +++ cl-gsl/poly.lisp Mon Apr 4 02:46:42 2005 @@ -26,10 +26,11 @@ :double)
(defun poly-eval (coefficients x) - (let ((c-ptr (lisp-vec->c-array coefficients))) - (prog1 - (gsl-poly-eval c-ptr (length coefficients) x) - (uffi:free-foreign-object c-ptr)))) + "Returns the value of the polynomial +c[0] + c[1] X + c[2] X^2 + ... + c[n-1] X^{n-1} +where COEFFICIENTS is a vector of the coefficients of length n." + (with-lisp-vec->c-array (c-ptr coefficients) + (gsl-poly-eval c-ptr (length coefficients) x)))
;; ----------------------------------------------------------------------
@@ -42,19 +43,20 @@ :int)
(defun solve-quadratic (a b c) + "Computes the real roots of the quadratic equation A x^2 + B x + C = 0. +Returns three values. The first two values are the real roots of the equation. +The third value is the number of roots (either 2 or 0). +If there are 0 real roots, the first two values are 0.0d0. When there are +2 real roots, their values are returned in ascending order." (declare (double-float a) (double-float b) (double-float c)) - (let ((x0-ptr (uffi:allocate-foreign-object :double)) - (x1-ptr (uffi:allocate-foreign-object :double)) - (num-roots) - (x0) - (x1)) - (declare (double-ptr-def x0-ptr) (double-ptr-def x1-ptr)) - (setq num-roots (gsl-poly-solve-quadratic a b c x0-ptr x1-ptr) - x0 (uffi:deref-pointer x0-ptr :double) - x1 (uffi:deref-pointer x1-ptr :double)) - (uffi:free-foreign-object x0-ptr) - (uffi:free-foreign-object x1-ptr) - (values x0 x1 num-roots))) + (uffi:with-foreign-object (x0-ptr :double) + (uffi:with-foreign-object (x1-ptr :double) + (setf (uffi:deref-pointer x0-ptr :double) 0.0d0) + (setf (uffi:deref-pointer x1-ptr :double) 0.0d0) + (let ((num-roots (gsl-poly-solve-quadratic a b c x0-ptr x1-ptr))) + (values (uffi:deref-pointer x0-ptr :double) + (uffi:deref-pointer x1-ptr :double) + num-roots)))))
;; ----------------------------------------------------------------------
@@ -68,24 +70,24 @@ :int)
(defun solve-cubic (a b c) + "Computes the real roots of the cubic equation, x^3 + A x^2 + B x + C = 0 +with a leading coefficient of unity. +Returns four values. The first 3 values are the real roots of the equation. +The fourth value is the number of real roots (either 1 or 3). +If 1 real root is found, the other two roots are 0.0d0. When 3 real +roots are found, they are returned in ascending order." (declare (double-float a) (double-float b) (double-float c)) - (let ((x0-ptr (uffi:allocate-foreign-object :double)) - (x1-ptr (uffi:allocate-foreign-object :double)) - (x2-ptr (uffi:allocate-foreign-object :double)) - (num-roots) - (x0) - (x1) - (x2)) - (declare (double-ptr-def x0-ptr) (double-ptr-def x1-ptr) - (double-ptr-def x2-ptr)) - (setq num-roots (gsl-poly-solve-cubic a b c x0-ptr x1-ptr x2-ptr) - x0 (uffi:deref-pointer x0-ptr :double) - x1 (uffi:deref-pointer x1-ptr :double) - x2 (uffi:deref-pointer x2-ptr :double)) - (uffi:free-foreign-object x0-ptr) - (uffi:free-foreign-object x1-ptr) - (uffi:free-foreign-object x2-ptr) - (values x0 x1 x2 num-roots))) + (uffi:with-foreign-object (x0-ptr :double) + (uffi:with-foreign-object (x1-ptr :double) + (uffi:with-foreign-object (x2-ptr :double) + (setf (uffi:deref-pointer x0-ptr :double) 0.0d0) + (setf (uffi:deref-pointer x1-ptr :double) 0.0d0) + (setf (uffi:deref-pointer x2-ptr :double) 0.0d0) + (let ((num-roots (gsl-poly-solve-cubic a b c x0-ptr x1-ptr x2-ptr))) + (values (uffi:deref-pointer x0-ptr :double) + (uffi:deref-pointer x1-ptr :double) + (uffi:deref-pointer x2-ptr :double) + num-roots))))))
;; ---------------------------------------------------------------------- @@ -100,18 +102,12 @@
(defun complex-solve-quadratic (a b c) (declare (double-float a) (double-float b) (double-float c)) - (let ((z0-ptr (uffi:allocate-foreign-object 'gsl-complex)) - (z1-ptr (uffi:allocate-foreign-object 'gsl-complex)) - (num-roots) - (z0) - (z1)) - (declare (gsl-complex-ptr-def z0-ptr) (gsl-complex-ptr-def z1-ptr)) - (setq num-roots (gsl-poly-complex-solve-quadratic a b c z0-ptr z1-ptr) - z0 (uffi:deref-pointer z0-ptr 'gsl-complex) - z1 (uffi:deref-pointer z1-ptr 'gsl-complex)) - (uffi:free-foreign-object z0-ptr) - (uffi:free-foreign-object z1-ptr) - (values (gsl-complex->complex z0) (gsl-complex->complex z1) num-roots))) + (uffi:with-foreign-object (z0-ptr 'gsl-complex) + (uffi:with-foreign-object (z1-ptr 'gsl-complex) + (let ((num-roots (gsl-poly-complex-solve-quadratic a b c z0-ptr z1-ptr))) + (values (gsl-complex->complex (uffi:deref-pointer z0-ptr 'gsl-complex)) + (gsl-complex->complex (uffi:deref-pointer z1-ptr 'gsl-complex)) + num-roots)))))
;; ----------------------------------------------------------------------
@@ -126,24 +122,15 @@
(defun complex-solve-cubic (a b c) (declare (double-float a) (double-float b) (double-float c)) - (let ((z0-ptr (uffi:allocate-foreign-object 'gsl-complex)) - (z1-ptr (uffi:allocate-foreign-object 'gsl-complex)) - (z2-ptr (uffi:allocate-foreign-object 'gsl-complex)) - (num-roots) - (z0) - (z1) - (z2)) - (declare (gsl-complex-ptr-def z0-ptr) (gsl-complex-ptr-def z1-ptr) - (gsl-complex-ptr-def z2-ptr)) - (setq num-roots (gsl-poly-complex-solve-cubic a b c z0-ptr z1-ptr z2-ptr) - z0 (uffi:deref-pointer z0-ptr 'gsl-complex) - z1 (uffi:deref-pointer z1-ptr 'gsl-complex) - z2 (uffi:deref-pointer z2-ptr 'gsl-complex)) - (uffi:free-foreign-object z0-ptr) - (uffi:free-foreign-object z1-ptr) - (uffi:free-foreign-object z2-ptr) - (values (gsl-complex->complex z0) (gsl-complex->complex z1) - (gsl-complex->complex z2) num-roots))) + (uffi:with-foreign-object (z0-ptr 'gsl-complex) + (uffi:with-foreign-object (z1-ptr 'gsl-complex) + (uffi:with-foreign-object (z2-ptr 'gsl-complex) + (let ((num-roots (gsl-poly-complex-solve-cubic a b c + z0-ptr z1-ptr z2-ptr))) + (values (gsl-complex->complex (uffi:deref-pointer z0-ptr 'gsl-complex)) + (gsl-complex->complex (uffi:deref-pointer z1-ptr 'gsl-complex)) + (gsl-complex->complex (uffi:deref-pointer z2-ptr 'gsl-complex)) + num-roots))))))
;; ----------------------------------------------------------------------
@@ -163,16 +150,15 @@ :int)
(defun complex-solve (a) - (let* ((a-ptr (lisp-vec->c-array a)) - (len (length a)) - (w (gsl-poly-complex-workspace-alloc len)) - (z-ptr (uffi:allocate-foreign-object :double (* 2 (1- len)))) - (ret-val)) - (setq ret-val (gsl-poly-complex-solve a-ptr len w z-ptr)) - (gsl-poly-complex-workspace-free w) - (multiple-value-prog1 - (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val) - (uffi:free-foreign-object z-ptr)))) + (with-lisp-vec->c-array (a-ptr a) + (let* ((len (length a)) + (w (gsl-poly-complex-workspace-alloc len)) + (z-ptr (uffi:allocate-foreign-object :double (* 2 (1- len)))) + (ret-val (gsl-poly-complex-solve a-ptr len w z-ptr))) + (gsl-poly-complex-workspace-free w) + (multiple-value-prog1 + (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val) + (uffi:free-foreign-object z-ptr)))))
;; ----------------------------------------------------------------------
@@ -183,18 +169,20 @@ (size :unsigned-long)) :int)
+ (defun dd-init (xa ya) - (let* ((xa-ptr (lisp-vec->c-array xa)) - (ya-ptr (lisp-vec->c-array ya)) - (len (length xa)) - (dd-ptr (uffi:allocate-foreign-object :double len)) - (ret-val)) - (setq ret-val (gsl-poly-dd-init dd-ptr xa-ptr ya-ptr len)) - (multiple-value-prog1 - (values (c-array->lisp-vec dd-ptr len) ret-val) - (uffi:free-foreign-object xa-ptr) - (uffi:free-foreign-object ya-ptr) - (uffi:free-foreign-object dd-ptr)))) + "Computes a divided-difference representation of the interpolating polynomial +for the points (xa, ya) stored in the vectors XA and YA of equal length. +Returns two values: the divided differences as a vector of length equal to XA, +and the status, indicating the success of the computation." + (with-lisp-vec->c-array (xa-ptr xa) + (with-lisp-vec->c-array (ya-ptr ya) + (let* ((len (length xa)) + (dd-ptr (uffi:allocate-foreign-object :double len)) + (ret-val (gsl-poly-dd-init dd-ptr xa-ptr ya-ptr len))) + (multiple-value-prog1 + (values (c-array->lisp-vec dd-ptr len) ret-val) + (uffi:free-foreign-object dd-ptr))))))
;; ----------------------------------------------------------------------
@@ -205,14 +193,13 @@ (x :double)) :double)
+ (defun dd-eval (dd xa x) - (let ((dd-ptr (lisp-vec->c-array dd)) - (xa-ptr (lisp-vec->c-array xa)) - (len (length dd))) - (prog1 - (gsl-poly-dd-eval dd-ptr xa-ptr len x) - (uffi:free-foreign-object xa-ptr) - (uffi:free-foreign-object dd-ptr)))) + "Returns the value of the polynomial at point X. The vectors DD and XA, +of equal length, store the divided difference representation of the polynomial." + (with-lisp-vec->c-array (dd-ptr dd) + (with-lisp-vec->c-array (xa-ptr xa) + (gsl-poly-dd-eval dd-ptr xa-ptr (length dd) x))))
;; ----------------------------------------------------------------------
@@ -225,17 +212,19 @@ (w double-ptr)) :int)
+ (defun dd-taylor (xp dd xa) - (let* ((dd-ptr (lisp-vec->c-array dd)) - (xa-ptr (lisp-vec->c-array xa)) - (len (length dd)) - (w-ptr (uffi:allocate-foreign-object :double len)) - (c-ptr (uffi:allocate-foreign-object :double len)) - (ret-val)) - (setq ret-val (gsl-poly-dd-taylor c-ptr xp dd-ptr xa-ptr len w-ptr)) - (multiple-value-prog1 - (values (c-array->lisp-vec c-ptr len) ret-val) - (uffi:free-foreign-object w-ptr) - (uffi:free-foreign-object xa-ptr) - (uffi:free-foreign-object dd-ptr) - (uffi:free-foreign-object c-ptr)))) + "Converts the divided-difference representation of a polynomial to +a Taylor expansion. The divided-difference representation is supplied in the +vectors DD and XA of equal length. Returns a vector of the Taylor coefficients +of the polynomial expanded about the point XP." + (with-lisp-vec->c-array (dd-ptr dd) + (with-lisp-vec->c-array (xa-ptr xa) + (let* ((len (length dd)) + (w-ptr (uffi:allocate-foreign-object :double len)) + (c-ptr (uffi:allocate-foreign-object :double len)) + (ret-val (gsl-poly-dd-taylor c-ptr xp dd-ptr xa-ptr len w-ptr))) + (multiple-value-prog1 + (values (c-array->lisp-vec c-ptr len) ret-val) + (uffi:free-foreign-object w-ptr) + (uffi:free-foreign-object c-ptr))))))