# Diff of /cl-gsl/poly.lisp

revision 1.2 by edenny, Sun Mar 13 00:48:25 2005 UTC revision 1.3 by edenny, Mon Apr 4 00:46:42 2005 UTC
# Line 26  Line 26
26    :double)    :double)
27
28  (defun poly-eval (coefficients x)  (defun poly-eval (coefficients x)
29    (let ((c-ptr (lisp-vec->c-array coefficients)))    "Returns the value of the polynomial
30      (prog1  c[0] + c[1] X + c[2] X^2 + ... + c[n-1] X^{n-1}
31          (gsl-poly-eval c-ptr (length coefficients) x)  where COEFFICIENTS is a vector of the coefficients of length n."
32        (uffi:free-foreign-object c-ptr))))    (with-lisp-vec->c-array (c-ptr coefficients)
33        (gsl-poly-eval c-ptr (length coefficients) x)))
34
35  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
36
# Line 42  Line 43
43    :int)    :int)
44
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))    (declare (double-float a) (double-float b) (double-float c))
52    (let ((x0-ptr (uffi:allocate-foreign-object :double))    (uffi:with-foreign-object (x0-ptr :double)
53          (x1-ptr (uffi:allocate-foreign-object :double))      (uffi:with-foreign-object (x1-ptr :double)
54          (num-roots)        (setf (uffi:deref-pointer x0-ptr :double) 0.0d0)
55          (x0)        (setf (uffi:deref-pointer x1-ptr :double) 0.0d0)
56          (x1))        (let ((num-roots (gsl-poly-solve-quadratic a b c x0-ptr x1-ptr)))
57      (declare (double-ptr-def x0-ptr) (double-ptr-def x1-ptr))          (values (uffi:deref-pointer x0-ptr :double)
58      (setq num-roots (gsl-poly-solve-quadratic a b c x0-ptr x1-ptr)                  (uffi:deref-pointer x1-ptr :double)
59            x0 (uffi:deref-pointer x0-ptr :double)                  num-roots)))))
x1 (uffi:deref-pointer x1-ptr :double))
(uffi:free-foreign-object x0-ptr)
(uffi:free-foreign-object x1-ptr)
(values x0 x1 num-roots)))
60
61  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
62
# Line 68  Line 70
70    :int)    :int)
71
72  (defun solve-cubic (a b c)  (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))    (declare (double-float a) (double-float b) (double-float c))
80    (let ((x0-ptr (uffi:allocate-foreign-object :double))    (uffi:with-foreign-object (x0-ptr :double)
81          (x1-ptr (uffi:allocate-foreign-object :double))      (uffi:with-foreign-object (x1-ptr :double)
82          (x2-ptr (uffi:allocate-foreign-object :double))        (uffi:with-foreign-object (x2-ptr :double)
83          (num-roots)          (setf (uffi:deref-pointer x0-ptr :double) 0.0d0)
84          (x0)          (setf (uffi:deref-pointer x1-ptr :double) 0.0d0)
85          (x1)          (setf (uffi:deref-pointer x2-ptr :double) 0.0d0)
86          (x2))          (let ((num-roots (gsl-poly-solve-cubic a b c x0-ptr x1-ptr x2-ptr)))
87      (declare (double-ptr-def x0-ptr) (double-ptr-def x1-ptr)            (values (uffi:deref-pointer x0-ptr :double)
88               (double-ptr-def x2-ptr))                    (uffi:deref-pointer x1-ptr :double)
89      (setq num-roots (gsl-poly-solve-cubic a b c x0-ptr x1-ptr x2-ptr)                    (uffi:deref-pointer x2-ptr :double)
90            x0 (uffi:deref-pointer x0-ptr :double)                    num-roots))))))
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)))
91
92
93  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
# Line 100  Line 102
102
104    (declare (double-float a) (double-float b) (double-float c))    (declare (double-float a) (double-float b) (double-float c))
105    (let ((z0-ptr (uffi:allocate-foreign-object 'gsl-complex))    (uffi:with-foreign-object (z0-ptr 'gsl-complex)
106          (z1-ptr (uffi:allocate-foreign-object 'gsl-complex))      (uffi:with-foreign-object (z1-ptr 'gsl-complex)
107          (num-roots)        (let ((num-roots (gsl-poly-complex-solve-quadratic a b c z0-ptr z1-ptr)))
108          (z0)          (values (gsl-complex->complex (uffi:deref-pointer z0-ptr 'gsl-complex))
109          (z1))                  (gsl-complex->complex (uffi:deref-pointer z1-ptr 'gsl-complex))
110      (declare (gsl-complex-ptr-def z0-ptr) (gsl-complex-ptr-def z1-ptr))                  num-roots)))))
(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)))
111
112  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
113
# Line 126  Line 122
122
123  (defun complex-solve-cubic (a b c)  (defun complex-solve-cubic (a b c)
124    (declare (double-float a) (double-float b) (double-float c))    (declare (double-float a) (double-float b) (double-float c))
125    (let ((z0-ptr (uffi:allocate-foreign-object 'gsl-complex))    (uffi:with-foreign-object (z0-ptr 'gsl-complex)
126          (z1-ptr (uffi:allocate-foreign-object 'gsl-complex))      (uffi:with-foreign-object (z1-ptr 'gsl-complex)
127          (z2-ptr (uffi:allocate-foreign-object 'gsl-complex))        (uffi:with-foreign-object (z2-ptr 'gsl-complex)
128          (num-roots)          (let ((num-roots (gsl-poly-complex-solve-cubic a b c
129          (z0)                                                         z0-ptr z1-ptr z2-ptr)))
130          (z1)            (values (gsl-complex->complex (uffi:deref-pointer z0-ptr 'gsl-complex))
131          (z2))                    (gsl-complex->complex (uffi:deref-pointer z1-ptr 'gsl-complex))
132      (declare (gsl-complex-ptr-def z0-ptr) (gsl-complex-ptr-def z1-ptr)                    (gsl-complex->complex (uffi:deref-pointer z2-ptr 'gsl-complex))
133               (gsl-complex-ptr-def z2-ptr))                    num-roots))))))
(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)))
134
135  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
136
# Line 163  Line 150
150    :int)    :int)
151
152  (defun complex-solve (a)  (defun complex-solve (a)
153    (let* ((a-ptr (lisp-vec->c-array a))    (with-lisp-vec->c-array (a-ptr a)
154           (len (length a))      (let* ((len (length a))
155           (w (gsl-poly-complex-workspace-alloc len))             (w (gsl-poly-complex-workspace-alloc len))
156           (z-ptr (uffi:allocate-foreign-object :double (* 2 (1- len))))             (z-ptr (uffi:allocate-foreign-object :double (* 2 (1- len))))
157           (ret-val))             (ret-val (gsl-poly-complex-solve a-ptr len w z-ptr)))
158      (setq ret-val (gsl-poly-complex-solve a-ptr len w z-ptr))        (gsl-poly-complex-workspace-free w)
159      (gsl-poly-complex-workspace-free w)        (multiple-value-prog1
160      (multiple-value-prog1            (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val)
161          (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val)          (uffi:free-foreign-object z-ptr)))))
(uffi:free-foreign-object z-ptr))))
162
163  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
164
# Line 183  Line 169
169       (size :unsigned-long))       (size :unsigned-long))
170    :int)    :int)
171
172
173  (defun dd-init (xa ya)  (defun dd-init (xa ya)
174    (let* ((xa-ptr (lisp-vec->c-array xa))    "Computes a divided-difference representation of the interpolating polynomial
175           (ya-ptr (lisp-vec->c-array ya))  for the points (xa, ya) stored in the vectors XA and YA of equal length.
176           (len (length xa))  Returns two values: the divided differences as a vector of length equal to XA,
177           (dd-ptr (uffi:allocate-foreign-object :double len))  and the status, indicating the success of the computation."
178           (ret-val))    (with-lisp-vec->c-array (xa-ptr xa)
179      (setq ret-val (gsl-poly-dd-init dd-ptr xa-ptr ya-ptr len))      (with-lisp-vec->c-array (ya-ptr ya)
180      (multiple-value-prog1        (let* ((len (length xa))
181          (values (c-array->lisp-vec dd-ptr len) ret-val)               (dd-ptr (uffi:allocate-foreign-object :double len))
182        (uffi:free-foreign-object xa-ptr)               (ret-val (gsl-poly-dd-init dd-ptr xa-ptr ya-ptr len)))
183        (uffi:free-foreign-object ya-ptr)          (multiple-value-prog1
184        (uffi:free-foreign-object dd-ptr))))              (values (c-array->lisp-vec dd-ptr len) ret-val)
185              (uffi:free-foreign-object dd-ptr))))))
186
187  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
188
# Line 205  Line 193
193       (x :double))       (x :double))
194    :double)    :double)
195
196
197  (defun dd-eval (dd xa x)  (defun dd-eval (dd xa x)
198    (let ((dd-ptr (lisp-vec->c-array dd))    "Returns the value of the polynomial at point X. The vectors DD and XA,
199          (xa-ptr (lisp-vec->c-array xa))  of equal length, store the divided difference representation of the polynomial."
200          (len (length dd)))    (with-lisp-vec->c-array (dd-ptr dd)
201      (prog1      (with-lisp-vec->c-array (xa-ptr xa)
202          (gsl-poly-dd-eval dd-ptr xa-ptr len x)        (gsl-poly-dd-eval dd-ptr xa-ptr (length dd) x))))
(uffi:free-foreign-object xa-ptr)
(uffi:free-foreign-object dd-ptr))))
203
204  ;; ----------------------------------------------------------------------  ;; ----------------------------------------------------------------------
205
# Line 225  Line 212
212       (w double-ptr))       (w double-ptr))
213    :int)    :int)
214
215
216  (defun dd-taylor (xp dd xa)  (defun dd-taylor (xp dd xa)
217    (let* ((dd-ptr (lisp-vec->c-array dd))    "Converts the divided-difference representation of a polynomial to
218           (xa-ptr (lisp-vec->c-array xa))  a Taylor expansion. The divided-difference representation is supplied in the
219           (len (length dd))  vectors DD and XA of equal length. Returns a vector of the Taylor coefficients
220           (w-ptr (uffi:allocate-foreign-object :double len))  of the polynomial expanded about the point XP."
221           (c-ptr (uffi:allocate-foreign-object :double len))    (with-lisp-vec->c-array (dd-ptr dd)
222           (ret-val))      (with-lisp-vec->c-array (xa-ptr xa)
223      (setq ret-val (gsl-poly-dd-taylor c-ptr xp dd-ptr xa-ptr len w-ptr))        (let* ((len (length dd))
224      (multiple-value-prog1               (w-ptr (uffi:allocate-foreign-object :double len))
225          (values (c-array->lisp-vec c-ptr len) ret-val)               (c-ptr (uffi:allocate-foreign-object :double len))
226        (uffi:free-foreign-object w-ptr)               (ret-val (gsl-poly-dd-taylor c-ptr xp dd-ptr xa-ptr len w-ptr)))
227        (uffi:free-foreign-object xa-ptr)          (multiple-value-prog1
228        (uffi:free-foreign-object dd-ptr)              (values (c-array->lisp-vec c-ptr len) ret-val)
229        (uffi:free-foreign-object c-ptr))))            (uffi:free-foreign-object w-ptr)
230              (uffi:free-foreign-object c-ptr))))))

Legend:
 Removed from v.1.2 changed lines Added in v.1.3