April 2005 - cl-gsl-cvs@common-lisp.net - mailman3.common-lisp.net

[cl-gsl-cvs] CVS update: cl-gsl/vector.lisp
by cl-gsl-cvs＠common-lisp.net 07 Apr '05

07 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv2175 Modified Files: vector.lisp Log Message: Fixes as a result of unit tests. Date: Thu Apr 7 04:37:14 2005 Author: edenny Index: cl-gsl/vector.lisp diff -u cl-gsl/vector.lisp:1.5 cl-gsl/vector.lisp:1.6 --- cl-gsl/vector.lisp:1.5 Mon Apr 4 02:47:39 2005 +++ cl-gsl/vector.lisp Thu Apr 7 04:37:13 2005 @@ -332,6 +332,7 @@ (error "No matching type")))) +;; TODO: make a (setf (get-element v i) x) version. (defun set-element (v i x) (assert (eq 'gsl-vec (type-of v))) (assert (typep x (gsl-vec-element-type v))) @@ -351,7 +352,8 @@ (with-complex-double-float->gsl-complex-ptr (c-ptr x) (wrap-gsl-vector-complex-set (gsl-vec-ptr v) i c-ptr))) (t - (error "No matching type")))) + (error "No matching type"))) + v) (defun set-all (v x) @@ -371,7 +373,8 @@ (with-complex-double-float->gsl-complex-ptr (c-ptr x) (wrap-gsl-vector-complex-set-all (gsl-vec-ptr v) c-ptr))) (t - (error "No matching type")))) + (error "No matching type"))) + v) (defun set-zero (v) @@ -388,7 +391,8 @@ ((equal (gsl-vec-element-type v) '(complex (double-float))) (gsl-vector-complex-set-zero (gsl-vec-ptr v))) (t - (error "No matching type")))) + (error "No matching type"))) + v) (defun set-basis (v i) @@ -407,11 +411,57 @@ ((equal (gsl-vec-element-type v) '(complex (double-float))) (gsl-vector-complex-set-basis (gsl-vec-ptr v) i)) (t - (error "No matching type")))) + (error "No matching type"))) + v) + + +(defun read-from-binary-file (v file-name size) + (assert (eq 'gsl-vec (type-of v))) + (assert (<= size (gsl-vec-size v))) + (let ((status)) + (uffi:with-cstring (c-file-name file-name) + (setq status + (cond + ((eq (gsl-vec-element-type v) 'integer) + (wrap-gsl-vector-int-fread c-file-name (gsl-vec-ptr v))) + ((eq (gsl-vec-element-type v) 'single-float) + (wrap-gsl-vector-float-fread c-file-name (gsl-vec-ptr v))) + ((eq (gsl-vec-element-type v) 'double-float) + (wrap-gsl-vector-fread c-file-name (gsl-vec-ptr v))) + ((equal (gsl-vec-element-type v) '(complex (single-float))) + (wrap-gsl-vector-complex-float-fread c-file-name (gsl-vec-ptr v))) + ((equal (gsl-vec-element-type v) '(complex (double-float))) + (wrap-gsl-vector-complex-fread c-file-name (gsl-vec-ptr v))) + (t + (error "No matching type"))))) + (values v status))) + + +(defun read-from-file (v file-name size) + (assert (eq 'gsl-vec (type-of v))) + (assert (<= size (gsl-vec-size v))) + (let ((status)) + (uffi:with-cstring (c-file-name file-name) + (setq status + (cond + ((eq (gsl-vec-element-type v) 'integer) + (wrap-gsl-vector-int-fscanf c-file-name (gsl-vec-ptr v))) + ((eq (gsl-vec-element-type v) 'single-float) + (wrap-gsl-vector-float-fscanf c-file-name (gsl-vec-ptr v))) + ((eq (gsl-vec-element-type v) 'double-float) + (wrap-gsl-vector-fscanf c-file-name (gsl-vec-ptr v))) + ((equal (gsl-vec-element-type v) '(complex (single-float))) + (wrap-gsl-vector-complex-float-fscanf c-file-name + (gsl-vec-ptr v))) + ((equal (gsl-vec-element-type v) '(complex (double-float))) + (wrap-gsl-vector-complex-fscanf c-file-name (gsl-vec-ptr v))) + (t + (error "No matching type"))))) + (values v status))) (defun make-vector (size &key (element-type 'double-float) initial-element - initial-contents) + initial-contents from-file from-binary-file) (assert (typep size 'integer)) (assert (find element-type '(integer single-float double-float (complex (single-float)) @@ -419,8 +469,8 @@ (let ((v (make-gsl-vec :size size :element-type element-type))) (setf (gsl-vec-ptr v) (alloc v)) (cond - ((and initial-element initial-contents) - (error "cannot define both initial-element and initial-contents keys")) + ((and initial-element initial-contents from-file from-binary-file) + (error "can only define one of the keys: initial-element, initial-contents, from-file, from-binary-file.")) (initial-element (gsl-vector:set-all v initial-element)) (initial-contents @@ -435,18 +485,25 @@ ((= i size)) (gsl-vector:set-element v i (aref initial-contents i)))) (t - (error "initial-contents must be either a list or a vector."))))) + (error "initial-contents must be either a list or a vector.")))) + (from-file + (read-from-file v from-file size)) + (from-binary-file + (read-from-binary-file v from-binary-file size))) v)) -(defmacro with-vector ((vec size &key element-type initial-element - initial-contents) &body body) +(defmacro with-vector + ((vec size &key element-type initial-element initial-contents from-file + from-binary-file) &body body) `(let ((,vec (make-vector ,size :element-type (or ,element-type 'double-float) :initial-element ,initial-element - :initial-contents ,initial-contents))) + :initial-contents ,initial-contents + :from-file ,from-file + :from-binary-file ,from-binary-file))) (unwind-protect - ,@body + (progn ,@body) (free ,vec)))) @@ -495,49 +552,6 @@ status)) -(defun read-from-binary-file (file-name size element-type) - (let ((v (make-vector size :element-type element-type)) - (status)) - (uffi:with-cstring (c-file-name file-name) - (setq status - (cond - ((eq (gsl-vec-element-type v) 'integer) - (wrap-gsl-vector-int-fread c-file-name (gsl-vec-ptr v))) - ((eq (gsl-vec-element-type v) 'single-float) - (wrap-gsl-vector-float-fread c-file-name (gsl-vec-ptr v))) - ((eq (gsl-vec-element-type v) 'double-float) - (wrap-gsl-vector-fread c-file-name (gsl-vec-ptr v))) - ((equal (gsl-vec-element-type v) '(complex (single-float))) - (wrap-gsl-vector-complex-float-fread c-file-name (gsl-vec-ptr v))) - ((equal (gsl-vec-element-type v) '(complex (double-float))) - (wrap-gsl-vector-complex-fread c-file-name (gsl-vec-ptr v))) - (t - (error "No matching type"))))) - (values v status))) - - -(defun read-from-file (file-name size element-type) - (let ((v (make-vector size :element-type element-type)) - (status)) - (uffi:with-cstring (c-file-name file-name) - (setq status - (cond - ((eq (gsl-vec-element-type v) 'integer) - (wrap-gsl-vector-int-fscanf c-file-name (gsl-vec-ptr v))) - ((eq (gsl-vec-element-type v) 'single-float) - (wrap-gsl-vector-float-fscanf c-file-name (gsl-vec-ptr v))) - ((eq (gsl-vec-element-type v) 'double-float) - (wrap-gsl-vector-fscanf c-file-name (gsl-vec-ptr v))) - ((equal (gsl-vec-element-type v) '(complex (single-float))) - (wrap-gsl-vector-complex-float-fscanf c-file-name - (gsl-vec-ptr v))) - ((equal (gsl-vec-element-type v) '(complex (double-float))) - (wrap-gsl-vector-complex-fscanf c-file-name (gsl-vec-ptr v))) - (t - (error "No matching type"))))) - (values v status))) - - (defun subvector (v offset n) (assert (eq 'gsl-vec (type-of v))) (assert (typep offset 'integer)) @@ -891,6 +905,13 @@ (gsl-vector-complex-isnull (gsl-vec-ptr v))) (t (error "No matching type"))))) + + +(defun gsl-vector->lisp-vector (v) + (assert (eq 'gsl-vec (type-of v))) + (let ((a (make-array (gsl-vec-size v) :element-type (gsl-vec-element-type v)))) + (dotimes (i (gsl-vec-size v) a) + (setf (aref a i) (get-element v i))))) ;; Function: gsl_vector_view gsl_vector_complex_real (gsl_vector_complex *v) ;; Function: gsl_vector_view gsl_vector_complex_imag (gsl_vector_complex *v)

1 0

[cl-gsl-cvs] CVS update: cl-gsl/package.lisp
by cl-gsl-cvs＠common-lisp.net 07 Apr '05

07 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv2156 Modified Files: package.lisp Log Message: Add additional symbols to vector package. Date: Thu Apr 7 04:36:20 2005 Author: edenny Index: cl-gsl/package.lisp diff -u cl-gsl/package.lisp:1.4 cl-gsl/package.lisp:1.5 --- cl-gsl/package.lisp:1.4 Mon Apr 4 02:45:25 2005 +++ cl-gsl/package.lisp Thu Apr 7 04:36:20 2005 @@ -133,4 +133,5 @@ #:min-max-indicies #:min-max-values #:isnull + #:gsl-vector->lisp-vector ))

1 0

[cl-gsl-cvs] CVS update: cl-gsl/ChangeLog
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv2630 Modified Files: ChangeLog Log Message: *** empty log message *** Date: Mon Apr 4 04:02:35 2005 Author: edenny Index: cl-gsl/ChangeLog diff -u cl-gsl/ChangeLog:1.11 cl-gsl/ChangeLog:1.12 --- cl-gsl/ChangeLog:1.11 Mon Apr 4 02:53:28 2005 +++ cl-gsl/ChangeLog Mon Apr 4 04:02:35 2005 @@ -3,7 +3,7 @@ * doc/math.tex: Use \mathalt macro for equations. * doc/cl-gsl-ref.tex: - Add \mathalt macro. It create a html image alt that preserves the + Add \mathalt macro. It creates a html image alt that preserves the original latex. * doc/Makefile: Add another file type to clean dependency. @@ -18,7 +18,7 @@ * ffi.lisp: Added macros which allocate, bind, and free foreign structures. -2005-03-29 Edgar Denny <edenny(a)parenchyma.comcast.net> +2005-03-29 Edgar Denny <edenny(a)comcast.net> * doc/index.html, doc/introduction.html, doc/poly.html, doc/sf.html: Remove html files.

1 0

[cl-gsl-cvs] CVS update: cl-gsl/ChangeLog
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv30822 Modified Files: ChangeLog Log Message: *** empty log message *** Date: Mon Apr 4 02:53:28 2005 Author: edenny Index: cl-gsl/ChangeLog diff -u cl-gsl/ChangeLog:1.10 cl-gsl/ChangeLog:1.11 --- cl-gsl/ChangeLog:1.10 Sun Mar 27 05:49:46 2005 +++ cl-gsl/ChangeLog Mon Apr 4 02:53:28 2005 @@ -1,3 +1,28 @@ +2005-04-04 Edgar Denny <edgardenny(a)comcast.net> + + * doc/math.tex: Use \mathalt macro for equations. + + * doc/cl-gsl-ref.tex: + Add \mathalt macro. It create a html image alt that preserves the + original latex. + + * doc/Makefile: Add another file type to clean dependency. + + * vector.lisp: Add macros that automatically free foreign objects. + + * poly.lisp: Replace some functions with macros that clean up after + themselves. Plug a leak. + + * package.lisp: Add additional symbols to package cl-gsl-vector. + + * ffi.lisp: + Added macros which allocate, bind, and free foreign structures. + +2005-03-29 Edgar Denny <edenny(a)parenchyma.comcast.net> + + * doc/index.html, doc/introduction.html, doc/poly.html, doc/sf.html: + Remove html files. + 2005-03-27 Edgar Denny <edgardenny(a)comcast.net> * sf.tex, poly.tex, math.tex, constants.tex: tex4ht fixes.

1 0

[cl-gsl-cvs] CVS update: cl-gsl/doc/math.tex
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl/doc In directory common-lisp.net:/tmp/cvs-serv30800 Modified Files: math.tex Log Message: Use \mathalt macro for equations. Date: Mon Apr 4 02:52:04 2005 Author: edenny Index: cl-gsl/doc/math.tex diff -u cl-gsl/doc/math.tex:1.2 cl-gsl/doc/math.tex:1.3 --- cl-gsl/doc/math.tex:1.2 Sun Mar 27 05:47:52 2005 +++ cl-gsl/doc/math.tex Mon Apr 4 02:52:03 2005 @@ -6,38 +6,39 @@ \begin{description} -\item [\const{+e+}] The base of exponentials, $e$ +\item [\const{+e+}] The base of exponentials, \mathalt{$e$} -\item [\const{+log2e+}] The base-2 logarithm of e, $\log_2(e)$ +\item [\const{+log2e+}] The base-2 logarithm of e, \mathalt{$\log_2(e)$} -\item [\const{+log10e+}] The base-10 logarithm of e, $\log_{10}(e)$ +\item [\const{+log10e+}] The base-10 logarithm of e, \mathalt{$\log_{10}(e)$} -\item [\const{+sqrt2+}] The square root of two, $\sqrt 2$ +\item [\const{+sqrt2+}] The square root of two, \mathalt{$\sqrt 2$} -\item [\const{+sqrt1/2+}] The square root of one-half, $\sqrt{1/2}$ +\item [\const{+sqrt1/2+}] The square root of one-half,\mathalt{$\sqrt{1/2}$} -\item [\const{+sqrt3+}] The square root of three, $\sqrt 3$ +\item [\const{+sqrt3+}] The square root of three, \mathalt{$\sqrt 3$} -\item [\const{+pi+}] The constant pi, $\pi$ +\item [\const{+pi+}] The constant pi, \mathalt{$\pi$} -\item [\const{+pi/2+}] Pi divided by two, $\pi/2$ +\item [\const{+pi/2+}] Pi divided by two, \mathalt{$\pi/2$} -\item [\const{+pi/4+}] Pi divided by four, $\pi/4$ +\item [\const{+pi/4+}] Pi divided by four, \mathalt{$\pi/4$} -\item [\const{+sqrtpi+}] The square root of pi, $\sqrt\pi$ +\item [\const{+sqrtpi+}] The square root of pi, \mathalt{$\sqrt\pi$} -\item [\const{+2/sqrtpi+}] Two divided by the square root of pi, $2/\sqrt\pi$ +\item [\const{+2/sqrtpi+}] Two divided by the square root of pi, +\mathalt{$2/\sqrt\pi$} -\item [\const{+1/pi+}] The reciprocal of pi, $1/\pi$ +\item [\const{+1/pi+}] The reciprocal of pi, \mathalt{$1/\pi$} -\item [\const{+2/pi+}] Twice the reciprocal of pi, $2/\pi$ +\item [\const{+2/pi+}] Twice the reciprocal of pi, \mathalt{$2/\pi$} -\item [\const{+ln10+}] The natural logarithm of ten, $\ln(10)$ +\item [\const{+ln10+}] The natural logarithm of ten, \mathalt{$\ln(10)$} -\item [\const{+ln2+}] The natural logarithm of two, $\ln(2)$ +\item [\const{+ln2+}] The natural logarithm of two, \mathalt{$\ln(2)$} -\item [\const{+lnpi+}] The natural logarithm of pi, $\ln(\pi)$ +\item [\const{+lnpi+}] The natural logarithm of pi, \mathalt{$\ln(\pi)$} -\item [\const{+euler+}] Euler's constant, $\gamma$ +\item [\const{+euler+}] Euler's constant, \mathalt{$\gamma$} \end{description}

1 0

[cl-gsl-cvs] CVS update: cl-gsl/doc/cl-gsl-ref.tex
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl/doc In directory common-lisp.net:/tmp/cvs-serv30781 Added Files: cl-gsl-ref.tex Log Message: Add math alt macro. It create a html image alt that preserves the original latex. Date: Mon Apr 4 02:51:18 2005 Author: edenny

1 0

[cl-gsl-cvs] CVS update: cl-gsl/doc/Makefile
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl/doc In directory common-lisp.net:/tmp/cvs-serv30753 Modified Files: Makefile Log Message: Add another file type to clean dependency. Date: Mon Apr 4 02:48:27 2005 Author: edenny Index: cl-gsl/doc/Makefile diff -u cl-gsl/doc/Makefile:1.2 cl-gsl/doc/Makefile:1.3 --- cl-gsl/doc/Makefile:1.2 Sun Mar 27 05:46:32 2005 +++ cl-gsl/doc/Makefile Mon Apr 4 02:48:26 2005 @@ -3,7 +3,8 @@ #### GNU Makefile clean: - $(RM) *.png *.4ct *.4tc *.aux *.css *.dvi *.html *.idv *.lg *.log *.tmp *.xref *.toc + $(RM) *.png *.4ct *.4tc *.aux *.css *.dvi *.html *.idv *.lg *.log *.tmp *.xref *.toc *.idx + html: htlatex cl-gsl-ref.tex

1 0

[cl-gsl-cvs] CVS update: cl-gsl/vector.lisp
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv30732 Modified Files: vector.lisp Log Message: Add macros that automatically free foreign objects. Date: Mon Apr 4 02:47:40 2005 Author: edenny Index: cl-gsl/vector.lisp diff -u cl-gsl/vector.lisp:1.4 cl-gsl/vector.lisp:1.5 --- cl-gsl/vector.lisp:1.4 Tue Mar 15 04:17:29 2005 +++ cl-gsl/vector.lisp Mon Apr 4 02:47:39 2005 @@ -345,11 +345,11 @@ ((eq (gsl-vec-element-type v) 'double-float) (gsl-vector-set (gsl-vec-ptr v) i x)) ((equal (gsl-vec-element-type v) '(complex (single-float))) - (wrap-gsl-vector-complex-float-set (gsl-vec-ptr v) i - (complex->gsl-complex-float-ptr x))) + (with-complex-single-float->gsl-complex-float-ptr (c-ptr x) + (wrap-gsl-vector-complex-float-set (gsl-vec-ptr v) i c-ptr))) ((equal (gsl-vec-element-type v) '(complex (double-float))) - (wrap-gsl-vector-complex-set (gsl-vec-ptr v) i - (complex->gsl-complex-ptr x))) + (with-complex-double-float->gsl-complex-ptr (c-ptr x) + (wrap-gsl-vector-complex-set (gsl-vec-ptr v) i c-ptr))) (t (error "No matching type")))) @@ -365,11 +365,11 @@ ((eq (gsl-vec-element-type v) 'double-float) (gsl-vector-set-all (gsl-vec-ptr v) x)) ((equal (gsl-vec-element-type v) '(complex (single-float))) - (wrap-gsl-vector-complex-float-set-all (gsl-vec-ptr v) - (complex->gsl-complex-float-ptr x))) + (with-complex-single-float->gsl-complex-float-ptr (c-ptr x) + (wrap-gsl-vector-complex-float-set-all (gsl-vec-ptr v) c-ptr))) ((equal (gsl-vec-element-type v) '(complex (double-float))) - (wrap-gsl-vector-complex-set-all (gsl-vec-ptr v) - (complex->gsl-complex-ptr x))) + (with-complex-double-float->gsl-complex-ptr (c-ptr x) + (wrap-gsl-vector-complex-set-all (gsl-vec-ptr v) c-ptr))) (t (error "No matching type")))) @@ -439,6 +439,17 @@ v)) +(defmacro with-vector ((vec size &key element-type initial-element + initial-contents) &body body) + `(let ((,vec (make-vector ,size + :element-type (or ,element-type 'double-float) + :initial-element ,initial-element + :initial-contents ,initial-contents))) + (unwind-protect + ,@body + (free ,vec)))) + + (defun write-to-binary-file (file-name v) (assert (eq 'gsl-vec (type-of v))) (let ((status)) @@ -608,6 +619,13 @@ (t (error "No matching type"))))) (values v-dest status))) + + +(defmacro with-vector-copy ((vec-dest vec-src) &body body) + `(let ((,vec-dest (copy ,vec-src))) + (unwind-protect + ,@body + (free ,vec-dest)))) (defun swap (va vb)

1 0

[cl-gsl-cvs] CVS update: cl-gsl/poly.lisp
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

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))))))

1 0

[cl-gsl-cvs] CVS update: cl-gsl/package.lisp
by cl-gsl-cvs＠common-lisp.net 04 Apr '05

04 Apr '05

Update of /project/cl-gsl/cvsroot/cl-gsl In directory common-lisp.net:/tmp/cvs-serv30183 Modified Files: package.lisp Log Message: Add additional symbols to package cl-gsl-vector. Date: Mon Apr 4 02:45:25 2005 Author: edenny Index: cl-gsl/package.lisp diff -u cl-gsl/package.lisp:1.3 cl-gsl/package.lisp:1.4 --- cl-gsl/package.lisp:1.3 Tue Mar 15 04:15:57 2005 +++ cl-gsl/package.lisp Mon Apr 4 02:45:25 2005 @@ -38,10 +38,10 @@ #:def-foreign-struct #:gsl-complex->complex #:gsl-complex-float->complex - #:lisp-vec->c-array + #:with-lisp-vec->c-array #:complex-packed-array->lisp-vec - #:complex->gsl-complex-ptr - #:complex->gsl-complex-float-ptr + #:with-complex-double-float->gsl-complex-ptr + #:with-complex-single-float->gsl-complex-float-ptr #:c-array->lisp-vec #:defconstant-export #:register-constants @@ -69,6 +69,9 @@ #:complex-solve-quadratic #:complex-solve-cubic #:complex-solve + #:dd-init + #:dd-eval + #:dd-taylor )) (defpackage #:cl-gsl-sf @@ -100,6 +103,7 @@ #:free #:make-vector + #:with-vector #:get-element #:set-element #:set-all @@ -112,6 +116,7 @@ #:subvector #:subvector-with-stride #:copy + #:with-vector-copy #:swap #:swap-elements #:reverse-vector

1 0