[cl-gsl-cvs] CVS update: cl-gsl/poly.lisp
cl-gsl-cvs at common-lisp.net
cl-gsl-cvs at common-lisp.net
Mon Apr 4 00:46:44 UTC 2005
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))))))
More information about the Cl-gsl-cvs
mailing list