/[cl-gsl]/cl-gsl/poly.lisp
ViewVC logotype

Contents of /cl-gsl/poly.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.3 - (show annotations)
Mon Apr 4 00:46:42 2005 UTC (9 years ago) by edenny
Branch: MAIN
CVS Tags: HEAD
Changes since 1.2: +93 -104 lines
Replace some functions with macros that clean up after
themselves. Plug a leak.
1 ;;;; -*- Mode: Lisp; Synatx: ANSI-Common-Lisp; Base: 10 -*-
2 ;;;;
3 ;;;; Copyright (C) 2005 Edgar Denny <edgardenny@comcast.net>
4 ;;;; This file is part of CL-GSL.
5 ;;;;
6 ;;;; This program is free software; you can redistribute it and/or modify
7 ;;;; it under the terms of the GNU General Public License as published by
8 ;;;; the Free Software Foundation; either version 2 of the License, or
9 ;;;; (at your option) any later version.
10 ;;;;
11 ;;;; This program is distributed in the hope that it will be useful,
12 ;;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
13 ;;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 ;;;; GNU General Public License for more details.
15 ;;;;
16 ;;;; You should have received a copy of the GNU General Public License
17 ;;;; along with this program; if not, write to the Free Software
18 ;;;; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19
20 (in-package #:cl-gsl-poly)
21
22 (defun-foreign "gsl_poly_eval"
23 ((c (* :double))
24 (len :int)
25 (x :double))
26 :double)
27
28 (defun poly-eval (coefficients x)
29 "Returns the value of the polynomial
30 c[0] + c[1] X + c[2] X^2 + ... + c[n-1] X^{n-1}
31 where COEFFICIENTS is a vector of the coefficients of length n."
32 (with-lisp-vec->c-array (c-ptr coefficients)
33 (gsl-poly-eval c-ptr (length coefficients) x)))
34
35 ;; ----------------------------------------------------------------------
36
37 (defun-foreign "gsl_poly_solve_quadratic"
38 ((a :double)
39 (b :double)
40 (c :double)
41 (x0 double-ptr)
42 (x1 double-ptr))
43 :int)
44
45 (defun solve-quadratic (a b c)
46 "Computes the real roots of the quadratic equation A x^2 + B x + C = 0.
47 Returns three values. The first two values are the real roots of the equation.
48 The third value is the number of roots (either 2 or 0).
49 If there are 0 real roots, the first two values are 0.0d0. When there are
50 2 real roots, their values are returned in ascending order."
51 (declare (double-float a) (double-float b) (double-float c))
52 (uffi:with-foreign-object (x0-ptr :double)
53 (uffi:with-foreign-object (x1-ptr :double)
54 (setf (uffi:deref-pointer x0-ptr :double) 0.0d0)
55 (setf (uffi:deref-pointer x1-ptr :double) 0.0d0)
56 (let ((num-roots (gsl-poly-solve-quadratic a b c x0-ptr x1-ptr)))
57 (values (uffi:deref-pointer x0-ptr :double)
58 (uffi:deref-pointer x1-ptr :double)
59 num-roots)))))
60
61 ;; ----------------------------------------------------------------------
62
63 (defun-foreign "gsl_poly_solve_cubic"
64 ((a :double)
65 (b :double)
66 (c :double)
67 (x0 double-ptr)
68 (x1 double-ptr)
69 (x2 double-ptr))
70 :int)
71
72 (defun solve-cubic (a b c)
73 "Computes the real roots of the cubic equation, x^3 + A x^2 + B x + C = 0
74 with a leading coefficient of unity.
75 Returns four values. The first 3 values are the real roots of the equation.
76 The fourth value is the number of real roots (either 1 or 3).
77 If 1 real root is found, the other two roots are 0.0d0. When 3 real
78 roots are found, they are returned in ascending order."
79 (declare (double-float a) (double-float b) (double-float c))
80 (uffi:with-foreign-object (x0-ptr :double)
81 (uffi:with-foreign-object (x1-ptr :double)
82 (uffi:with-foreign-object (x2-ptr :double)
83 (setf (uffi:deref-pointer x0-ptr :double) 0.0d0)
84 (setf (uffi:deref-pointer x1-ptr :double) 0.0d0)
85 (setf (uffi:deref-pointer x2-ptr :double) 0.0d0)
86 (let ((num-roots (gsl-poly-solve-cubic a b c x0-ptr x1-ptr x2-ptr)))
87 (values (uffi:deref-pointer x0-ptr :double)
88 (uffi:deref-pointer x1-ptr :double)
89 (uffi:deref-pointer x2-ptr :double)
90 num-roots))))))
91
92
93 ;; ----------------------------------------------------------------------
94
95 (defun-foreign "gsl_poly_complex_solve_quadratic"
96 ((a :double)
97 (b :double)
98 (c :double)
99 (z0 gsl-complex-ptr)
100 (z1 gsl-complex-ptr))
101 :int)
102
103 (defun complex-solve-quadratic (a b c)
104 (declare (double-float a) (double-float b) (double-float c))
105 (uffi:with-foreign-object (z0-ptr 'gsl-complex)
106 (uffi:with-foreign-object (z1-ptr 'gsl-complex)
107 (let ((num-roots (gsl-poly-complex-solve-quadratic a b c z0-ptr z1-ptr)))
108 (values (gsl-complex->complex (uffi:deref-pointer z0-ptr 'gsl-complex))
109 (gsl-complex->complex (uffi:deref-pointer z1-ptr 'gsl-complex))
110 num-roots)))))
111
112 ;; ----------------------------------------------------------------------
113
114 (defun-foreign "gsl_poly_complex_solve_cubic"
115 ((a :double)
116 (b :double)
117 (c :double)
118 (z0 gsl-complex-ptr)
119 (z1 gsl-complex-ptr)
120 (z2 gsl-complex-ptr))
121 :int)
122
123 (defun complex-solve-cubic (a b c)
124 (declare (double-float a) (double-float b) (double-float c))
125 (uffi:with-foreign-object (z0-ptr 'gsl-complex)
126 (uffi:with-foreign-object (z1-ptr 'gsl-complex)
127 (uffi:with-foreign-object (z2-ptr 'gsl-complex)
128 (let ((num-roots (gsl-poly-complex-solve-cubic a b c
129 z0-ptr z1-ptr z2-ptr)))
130 (values (gsl-complex->complex (uffi:deref-pointer z0-ptr 'gsl-complex))
131 (gsl-complex->complex (uffi:deref-pointer z1-ptr 'gsl-complex))
132 (gsl-complex->complex (uffi:deref-pointer z2-ptr 'gsl-complex))
133 num-roots))))))
134
135 ;; ----------------------------------------------------------------------
136
137 (defun-foreign "gsl_poly_complex_workspace_alloc"
138 ((size :unsigned-long))
139 gsl-poly-complex-workspace-ptr)
140
141 (defun-foreign "gsl_poly_complex_workspace_free"
142 ((w gsl-poly-complex-workspace-ptr))
143 :void)
144
145 (defun-foreign "gsl_poly_complex_solve"
146 ((a double-ptr)
147 (size :unsigned-long)
148 (w gsl-poly-complex-workspace-ptr)
149 (z gsl-complex-packed-ptr))
150 :int)
151
152 (defun complex-solve (a)
153 (with-lisp-vec->c-array (a-ptr a)
154 (let* ((len (length a))
155 (w (gsl-poly-complex-workspace-alloc len))
156 (z-ptr (uffi:allocate-foreign-object :double (* 2 (1- len))))
157 (ret-val (gsl-poly-complex-solve a-ptr len w z-ptr)))
158 (gsl-poly-complex-workspace-free w)
159 (multiple-value-prog1
160 (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val)
161 (uffi:free-foreign-object z-ptr)))))
162
163 ;; ----------------------------------------------------------------------
164
165 (defun-foreign "gsl_poly_dd_init"
166 ((dd double-ptr)
167 (xa double-ptr)
168 (ya double-ptr)
169 (size :unsigned-long))
170 :int)
171
172
173 (defun dd-init (xa ya)
174 "Computes a divided-difference representation of the interpolating polynomial
175 for the points (xa, ya) stored in the vectors XA and YA of equal length.
176 Returns two values: the divided differences as a vector of length equal to XA,
177 and the status, indicating the success of the computation."
178 (with-lisp-vec->c-array (xa-ptr xa)
179 (with-lisp-vec->c-array (ya-ptr ya)
180 (let* ((len (length xa))
181 (dd-ptr (uffi:allocate-foreign-object :double len))
182 (ret-val (gsl-poly-dd-init dd-ptr xa-ptr ya-ptr len)))
183 (multiple-value-prog1
184 (values (c-array->lisp-vec dd-ptr len) ret-val)
185 (uffi:free-foreign-object dd-ptr))))))
186
187 ;; ----------------------------------------------------------------------
188
189 (defun-foreign "gsl_poly_dd_eval"
190 ((dd double-ptr)
191 (xa double-ptr)
192 (size :unsigned-long)
193 (x :double))
194 :double)
195
196
197 (defun dd-eval (dd xa x)
198 "Returns the value of the polynomial at point X. The vectors DD and XA,
199 of equal length, store the divided difference representation of the polynomial."
200 (with-lisp-vec->c-array (dd-ptr dd)
201 (with-lisp-vec->c-array (xa-ptr xa)
202 (gsl-poly-dd-eval dd-ptr xa-ptr (length dd) x))))
203
204 ;; ----------------------------------------------------------------------
205
206 (defun-foreign "gsl_poly_dd_taylor"
207 ((c double-ptr)
208 (xp :double)
209 (dd double-ptr)
210 (xa double-ptr)
211 (size :unsigned-long)
212 (w double-ptr))
213 :int)
214
215
216 (defun dd-taylor (xp dd xa)
217 "Converts the divided-difference representation of a polynomial to
218 a Taylor expansion. The divided-difference representation is supplied in the
219 vectors DD and XA of equal length. Returns a vector of the Taylor coefficients
220 of the polynomial expanded about the point XP."
221 (with-lisp-vec->c-array (dd-ptr dd)
222 (with-lisp-vec->c-array (xa-ptr xa)
223 (let* ((len (length dd))
224 (w-ptr (uffi:allocate-foreign-object :double len))
225 (c-ptr (uffi:allocate-foreign-object :double len))
226 (ret-val (gsl-poly-dd-taylor c-ptr xp dd-ptr xa-ptr len w-ptr)))
227 (multiple-value-prog1
228 (values (c-array->lisp-vec c-ptr len) ret-val)
229 (uffi:free-foreign-object w-ptr)
230 (uffi:free-foreign-object c-ptr))))))

  ViewVC Help
Powered by ViewVC 1.1.5