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

Contents of /cl-gsl/poly.lisp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations)
Wed Mar 2 01:04:53 2005 UTC (9 years, 1 month ago) by edenny
Branch: MAIN
Branch point for: cl-gsl
Initial revision
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 (let ((c-ptr (lisp-vec->c-array coefficients)))
30 (prog1
31 (gsl-poly-eval c-ptr (length coefficients) x)
32 (uffi:free-foreign-object c-ptr))))
33
34 ;; ----------------------------------------------------------------------
35
36 (defun-foreign "gsl_poly_solve_quadratic"
37 ((a :double)
38 (b :double)
39 (c :double)
40 (x0 double-ptr)
41 (x1 double-ptr))
42 :int)
43
44 (defun solve-quadratic (a b c)
45 (declare (double-float a) (double-float b) (double-float c))
46 (let ((x0-ptr (uffi:allocate-foreign-object :double))
47 (x1-ptr (uffi:allocate-foreign-object :double))
48 (num-roots)
49 (x0)
50 (x1))
51 (declare (double-ptr-def x0-ptr) (double-ptr-def x1-ptr))
52 (setq num-roots (gsl-poly-solve-quadratic a b c x0-ptr x1-ptr)
53 x0 (uffi:deref-pointer x0-ptr :double)
54 x1 (uffi:deref-pointer x1-ptr :double))
55 (uffi:free-foreign-object x0-ptr)
56 (uffi:free-foreign-object x1-ptr)
57 (values x0 x1 num-roots)))
58
59 ;; ----------------------------------------------------------------------
60
61 (defun-foreign "gsl_poly_solve_cubic"
62 ((a :double)
63 (b :double)
64 (c :double)
65 (x0 double-ptr)
66 (x1 double-ptr)
67 (x2 double-ptr))
68 :int)
69
70 (defun solve-cubic (a b c)
71 (declare (double-float a) (double-float b) (double-float c))
72 (let ((x0-ptr (uffi:allocate-foreign-object :double))
73 (x1-ptr (uffi:allocate-foreign-object :double))
74 (x2-ptr (uffi:allocate-foreign-object :double))
75 (num-roots)
76 (x0)
77 (x1)
78 (x2))
79 (declare (double-ptr-def x0-ptr) (double-ptr-def x1-ptr)
80 (double-ptr-def x2-ptr))
81 (setq num-roots (gsl-poly-solve-cubic a b c x0-ptr x1-ptr x2-ptr)
82 x0 (uffi:deref-pointer x0-ptr :double)
83 x1 (uffi:deref-pointer x1-ptr :double)
84 x2 (uffi:deref-pointer x2-ptr :double))
85 (uffi:free-foreign-object x0-ptr)
86 (uffi:free-foreign-object x1-ptr)
87 (uffi:free-foreign-object x2-ptr)
88 (values x0 x1 x2 num-roots)))
89
90
91 ;; ----------------------------------------------------------------------
92
93 (defun-foreign "gsl_poly_complex_solve_quadratic"
94 ((a :double)
95 (b :double)
96 (c :double)
97 (z0 gsl-complex-ptr)
98 (z1 gsl-complex-ptr))
99 :int)
100
101 (defun complex-solve-quadratic (a b c)
102 (declare (double-float a) (double-float b) (double-float c))
103 (let ((z0-ptr (uffi:allocate-foreign-object 'gsl-complex))
104 (z1-ptr (uffi:allocate-foreign-object 'gsl-complex))
105 (num-roots)
106 (z0)
107 (z1))
108 (declare (gsl-complex-ptr-def z0-ptr) (gsl-complex-ptr-def z1-ptr))
109 (setq num-roots (gsl-poly-complex-solve-quadratic a b c z0-ptr z1-ptr)
110 z0 (uffi:deref-pointer z0-ptr 'gsl-complex)
111 z1 (uffi:deref-pointer z1-ptr 'gsl-complex))
112 (uffi:free-foreign-object z0-ptr)
113 (uffi:free-foreign-object z1-ptr)
114 (values (gsl-complex->complex z0) (gsl-complex->complex z1) num-roots)))
115
116 ;; ----------------------------------------------------------------------
117
118 (defun-foreign "gsl_poly_complex_solve_cubic"
119 ((a :double)
120 (b :double)
121 (c :double)
122 (z0 gsl-complex-ptr)
123 (z1 gsl-complex-ptr)
124 (z2 gsl-complex-ptr))
125 :int)
126
127 (defun complex-solve-cubic (a b c)
128 (declare (double-float a) (double-float b) (double-float c))
129 (let ((z0-ptr (uffi:allocate-foreign-object 'gsl-complex))
130 (z1-ptr (uffi:allocate-foreign-object 'gsl-complex))
131 (z2-ptr (uffi:allocate-foreign-object 'gsl-complex))
132 (num-roots)
133 (z0)
134 (z1)
135 (z2))
136 (declare (gsl-complex-ptr-def z0-ptr) (gsl-complex-ptr-def z1-ptr)
137 (gsl-complex-ptr-def z2-ptr))
138 (setq num-roots (gsl-poly-complex-solve-cubic a b c z0-ptr z1-ptr z2-ptr)
139 z0 (uffi:deref-pointer z0-ptr 'gsl-complex)
140 z1 (uffi:deref-pointer z1-ptr 'gsl-complex)
141 z2 (uffi:deref-pointer z2-ptr 'gsl-complex))
142 (uffi:free-foreign-object z0-ptr)
143 (uffi:free-foreign-object z1-ptr)
144 (uffi:free-foreign-object z2-ptr)
145 (values (gsl-complex->complex z0) (gsl-complex->complex z1)
146 (gsl-complex->complex z2) num-roots)))
147
148 ;; ----------------------------------------------------------------------
149
150 (defun-foreign "gsl_poly_complex_workspace_alloc"
151 ((size :unsigned-long))
152 gsl-poly-complex-workspace-ptr)
153
154 (defun-foreign "gsl_poly_complex_workspace_free"
155 ((w gsl-poly-complex-workspace-ptr))
156 :void)
157
158 (defun-foreign "gsl_poly_complex_solve"
159 ((a double-ptr)
160 (size :unsigned-long)
161 (w gsl-poly-complex-workspace-ptr)
162 (z gsl-complex-packed-ptr))
163 :int)
164
165 (defun complex-solve (a)
166 (let* ((a-ptr (lisp-vec->c-array a))
167 (len (length a))
168 (w (gsl-poly-complex-workspace-alloc len))
169 (z-ptr (uffi:allocate-foreign-object :double (* 2 (1- len))))
170 (ret-val))
171 (setq ret-val (gsl-poly-complex-solve a-ptr len w z-ptr))
172 (gsl-poly-complex-workspace-free w)
173 (prog1
174 ;; FIXME: where is ret-val in result?
175 (values (complex-packed-array->lisp-vec z-ptr (* 2 (1- len))) ret-val)
176 (uffi:free-foreign-object z-ptr))))
177
178 ;; ----------------------------------------------------------------------
179
180 ;; (defun-foreign "gsl_poly_dd_int"
181 ;; ((dd double-ptr)
182 ;; (xa double-ptr)
183 ;; (ya double-ptr)
184 ;; (size :unsigned-long))
185 ;; :int)
186
187 ;; (defun-foreign "gsl_poly_dd_eval"
188 ;; ((dd double-ptr)
189 ;; (xa double-ptr)
190 ;; (size :unsigned-long)
191 ;; (x :double))
192 ;; :double)
193
194 ;; (defun-foreign "gsl_poly_dd_taylor"
195 ;; ((c double-ptr)
196 ;; (xp :double)
197 ;; (dd double-ptr)
198 ;; (xa double-ptr)
199 ;; (size :unsigned-long)
200 ;; (w double-ptr))
201 ;; :int)
202
203
204 ;; (defun dd-int (xa ya)
205 ;; (let* ((xa-ptr (lisp-vec->c-array xa))
206 ;; (ya-ptr (lisp-vec->c-array ya))
207 ;; (len (length xa))
208 ;; (dd-ptr (uffi:allocate-foreign-object :double len))
209 ;; (ret-val))
210 ;; (setq ret-val (gsl-poly-dd-init dd-ptr xa-ptr ya-ptr len))
211
212
213
214 ;; Function: int gsl_poly_dd_init (double dd[], const double xa[], const double ya[], size_t size)
215 ;; 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.
216 ;; On output the divided-differences of (xa,ya) are stored in the array dd, also of length size.
217
218 ;; Function: double gsl_poly_dd_eval (const double dd[], const double xa[], const size_t size, const double x)
219 ;; This function evaluates the polynomial stored in divided-difference form in the arrays dd and xa of length size at the point x.
220
221 ;; Function: int gsl_poly_dd_taylor (double c[], double xp, const double dd[], const double xa[], size_t size, double w[])
222 ;; This function converts the divided-difference representation of a polynomial to a Taylor expansion. The divided-difference representation is supplied in the
223 ;; 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
224 ;; workspace of length size must be provided in the array w.
225
226
227
228
229
230 ;; #define GSL_REAL(z) ((z).dat[0])
231 ;; #define GSL_IMAG(z) ((z).dat[1])
232 ;; #define GSL_COMPLEX_P(zp) ((zp)->dat)
233 ;; #define GSL_COMPLEX_P_REAL(zp) ((zp)->dat[0])
234 ;; #define GSL_COMPLEX_P_IMAG(zp) ((zp)->dat[1])
235 ;; #define GSL_COMPLEX_EQ(z1,z2) (((z1).dat[0] == (z2).dat[0]) && ((z1).dat[1] == (z2).dat[1]))
236
237 ;; #define GSL_SET_COMPLEX(zp,x,y) do {(zp)->dat[0]=(x); (zp)->dat[1]=(y);} while(0)
238 ;; #define GSL_SET_REAL(zp,x) do {(zp)->dat[0]=(x);} while(0)
239 ;; #define GSL_SET_IMAG(zp,y) do {(zp)->dat[1]=(y);} while(0)
240
241 ;; #define GSL_SET_COMPLEX_PACKED(zp,n,x,y) do {*((zp)+2*(n))=(x); *((zp)+(2*(n)+1))=(y);} while(0)

  ViewVC Help
Powered by ViewVC 1.1.5