# Diff of /cl-gsl/poly.lisp

revision 1.1.1.1 by edenny, Wed Mar 2 01:04:53 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      (prog1            (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val)
161          ;; FIXME: where is ret-val in result?          (uffi:free-foreign-object z-ptr)))))
(values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val)
(uffi:free-foreign-object z-ptr))))

;; ----------------------------------------------------------------------

;; (defun-foreign "gsl_poly_dd_int"
;;     ((dd double-ptr)
;;      (xa double-ptr)
;;      (ya double-ptr)
;;      (size :unsigned-long))
;;   :int)

;; (defun-foreign "gsl_poly_dd_eval"
;;     ((dd double-ptr)
;;      (xa double-ptr)
;;      (size :unsigned-long)
;;      (x :double))
;;   :double)

;; (defun-foreign "gsl_poly_dd_taylor"
;;     ((c double-ptr)
;;      (xp :double)
;;      (dd double-ptr)
;;      (xa double-ptr)
;;      (size :unsigned-long)
;;      (w double-ptr))
;;   :int)

;; (defun dd-int (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))

;; Function: int gsl_poly_dd_init (double dd[], const double xa[], const double ya[], size_t size)
;;     This function computes a divided-difference representation of the interpolating polynomial for the points (xa, ya) stored in the arrays xa and ya of length size.
;;     On output the divided-differences of (xa,ya) are stored in the array dd, also of length size.

;; Function: double gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x)
;;     This function evaluates the polynomial stored in divided-difference form in the arrays dd and xa of length size at the point x.

;; Function: int gsl_poly_dd_taylor (double c[], double xp, const double dd[], const double xa[], size_t size, double w[])
;;     This function converts the divided-difference representation of a polynomial to a Taylor expansion. The divided-difference representation is supplied in the
;;     arrays dd and xa of length size. On output the Taylor coefficients of the polynomial expanded about the point xp are stored in the array c also of length size. A
;;     workspace of length size must be provided in the array w.
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
;; #define GSL_REAL(z)     ((z).dat[0])
;; #define GSL_IMAG(z)     ((z).dat[1])
;; #define GSL_COMPLEX_P(zp) ((zp)->dat)
;; #define GSL_COMPLEX_P_REAL(zp)  ((zp)->dat[0])
;; #define GSL_COMPLEX_P_IMAG(zp)  ((zp)->dat[1])
;; #define GSL_COMPLEX_EQ(z1,z2) (((z1).dat[0] == (z2).dat[0]) && ((z1).dat[1] == (z2).dat[1]))

;; #define GSL_SET_COMPLEX(zp,x,y) do {(zp)->dat[0]=(x); (zp)->dat[1]=(y);} while(0)
;; #define GSL_SET_REAL(zp,x) do {(zp)->dat[0]=(x);} while(0)
;; #define GSL_SET_IMAG(zp,y) do {(zp)->dat[1]=(y);} while(0)
215
216  ;; #define GSL_SET_COMPLEX_PACKED(zp,n,x,y) do {*((zp)+2*(n))=(x); *((zp)+(2*(n)+1))=(y);} while(0)  (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))))))

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