Newer
Older
;********************************************************
; file: polynomial.lisp
; description: Polynomials
; date: Tue Mar 21 2006 - 18:33
; author: Liam M. Healy
; modified: Sat Jun 10 2006 - 23:37
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; To do: finish divided differences, which requires figuring out how
;;; to handle raw C arrays, and deciding if/how to provide
;;; autotranslation from CL pure arrays.
;;;;****************************************************************************
;;;; Polynomial Evaluation
;;;;****************************************************************************
(defun-gsl polynomial-eval (coefficients x)
"gsl_poly_eval"
(((gsl-array coefficients) :pointer) ((dim0 coefficients) :size) (x :double))
:documentation
"Evaluate the polyonomial with coefficients at the point x."
:c-return :double)
;;;;****************************************************************************
;;;; Divided Difference Representation of Polynomials
;;;;****************************************************************************
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
;;; Use with-divided-difference to compute the divided difference,
;;; which may be passed to eval-divided-difference or
;;; taylor-divided-difference in the body.
(defun-gsl divided-difference-int (dd xa ya)
"gsl_poly_dd_init"
(((gsl-array dd) :pointer)
((gsl-array xa) :pointer) ((gsl-array ya) :pointer)
((dim0 xa) :size))
:return (dd)
:export nil
:index divided-difference)
(export '(divided-difference))
(defun divided-difference (xa ya)
"Compute a divided-difference representation of the
interpolating polynomial for the points (@var{xa}, @var{ya}) stored in
the arrays @var{xa} and @var{ya}. The output is the
divided-differences of (@var{xa},@var{ya}) stored in an gsl-vector
of the same length as xa and ya."
(let ((len (length xa)))
(with-data (xad vector-double len)
(with-data (yad vector-double len)
(setf (data xad) xa (data yad) ya)
(divided-difference-int
(make-data 'vector-double nil len)
xad yad)))))
(defun-gsl polynomial-eval-divided-difference (dd xa x)
"gsl_poly_dd_eval"
(((gsl-array dd) :pointer)
((gsl-array xa) :pointer)
((dim0 xa) :size)
(x :double))
:c-return :double
:documentation
"Evaluate the polynomial stored in divided-difference form
in the arrays @var{dd} and @var{xa} at the point @var{x}.")
#+development
(defun-gsl taylor-divided-difference (coefs dd xp workspace)
"gsl_poly_dd_taylor"
(((gsl-array coefs) :pointer)
((gsl-array xp) :pointer)
((gsl-array dd) :pointer)
((gsl-array x) :pointer)
((dim0 xa) :size)
((gsl-array workspace) :pointer))
:documentation
"Convert the divided-difference representation of a polynomial
to a Taylor expansion about the point xp. Call only within a
with-divided-difference form.")
;;;;****************************************************************************
;;;; Quadratic Equations
;;;;****************************************************************************
(defun-gsl solve-quadratic (a b c)
"gsl_poly_solve_quadratic"
((a :double) (b :double) (c :double) (root1 :double) (root2 :double))
:documentation
"The real roots of the quadratic equation a x^2 + b x + c = 0.
Two values are always returned; if the roots are not real, these
values are NIL."
:c-return :number-of-answers)
(defun-gsl solve-quadratic-complex (a b c)
"gsl_poly_complex_solve_quadratic"
((a :double) (b :double) (c :double) (root1 gsl-complex) (root2 gsl-complex))
"The complex roots of the quadratic equation a x^2 + b x + c = 0.
Two values are always returned; if a root does not exist, the
value returned will be NIL."
:c-return :number-of-answers)
;;;;****************************************************************************
;;;; Cubic Equations
;;;;****************************************************************************
(defun-gsl solve-cubic (a b c)
"gsl_poly_solve_cubic"
((a :double) (b :double) (c :double)
(root1 :double) (root2 :double) (root3 :double))
:documentation
"Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0
with a leading coefficient of unity. The roots are given
in ascending order. Three values are always returned;
if a root is not real, the value returned for it will be NIL."
:c-return :number-of-answers)
(defun-gsl solve-cubic-complex (a b c)
"gsl_poly_complex_solve_cubic"
((a :double) (b :double) (c :double)
(root1 gsl-complex) (root2 gsl-complex) (root3 gsl-complex))
:documentation
"Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0
with a leading coefficient of unity. Three values are always returned;
if a root does not exist, the value returned for it will be NIL."
:c-return :number-of-answers)
;;;;****************************************************************************
;;;; General Polynomial Equations
;;;;****************************************************************************
(defun-gsl complex-workspace-alloc (n)
"gsl_poly_complex_workspace_alloc" ((n :size))
:check-null-pointers
((:creturn :ENOMEM (format nil "No complex workspace")))
:c-return :pointer
:export nil
:index with-poly-complex-workspace)
(defun-gsl complex-workspace-free (ws)
"gsl_poly_complex_workspace_free" ((ws :pointer))
:c-return :void
:export nil
:index with-poly-complex-workspace)
(export '(with-poly-complex-workspace))
(defmacro with-poly-complex-workspace ((workspace size) &body body)
"Macro to create and cleanup workspace for polynomial root solver."
`(let ((,workspace (complex-workspace-alloc ,size)))
(unwind-protect
(progn ,@body)
(complex-workspace-free ,workspace))))
(defun-gsl polynomial-solve-ws (coefficients workspace answer-pd)
"gsl_poly_complex_solve"
(((gsl-array coefficients) :pointer) ((dim0 coefficients) :size)
(workspace :pointer) ((gsl-array answer-pd) :pointer))
:return
((loop for i from 0 below (dim0 answer-pd) by 2
collect (complex (gsl-aref answer-pd i)
(gsl-aref answer-pd (1+ i)))))
:documentation
"Arguments are:
a GSL array of coefficients, a workspace, a gsl-array of doubles."
:export nil
:index polynomial-solve)
(defun polynomial-solve (coefficients)
"The roots of the general polynomial
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
@math{P(x) = a_0 + a_1 x + a_2 x^2 + ... + a_@{n-1@} x^@{n-1@}} using
balanced-QR reduction of the companion matrix. The parameter @var{n}
specifies the length of the coefficient array. The coefficient of the
highest order term must be non-zero. The function requires a workspace
@var{w} of the appropriate size. The @math{n-1} roots are returned in
the packed complex array @var{z} of length @math{2(n-1)}, alternating
real and imaginary parts."
(let ((len (length coefficients)))
(with-data (coef vector-double len)
(setf (data coef) coefficients)
(with-data (answer vector-double ((* 2 (1- len))))
(with-poly-complex-workspace (ws len)
(values-list (polynomial-solve-ws coef ws answer)))))))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(lisp-unit:define-test polynomial
(lisp-unit:assert-first-fp-equal
"0.200000000000d+01"
(with-data (vec vector-double 3)
(setf (data vec) #(1.0d0 2.0d0 3.0d0))
(polynomial-eval vec -1.0d0)))
(lisp-unit:assert-equal
'(NIL NIL)
(multiple-value-list (solve-quadratic 1.0d0 0.0d0 1.0d0)))
(lisp-unit:assert-equal
'(1.0d0 1.0d0)
(multiple-value-list (solve-quadratic 1.0d0 -2.0d0 1.0d0)))
(lisp-unit:assert-equal
'(#C(1.0d0 0.0d0) #C(1.0d0 0.0d0))
(multiple-value-list (solve-quadratic-complex 1.0d0 -2.0d0 1.0d0)))
(lisp-unit:assert-equal
'("-0.300000000000d+01" "0.200000000000d+01" "0.700000000000d+01")
(lisp-unit:fp-values (solve-cubic -6.0d0 -13.0d0 42.0d0)))
;; This should use double-float-unequal
(lisp-unit:assert-equal
'(("-0.902598308594d-17" "-1.000000000000d+00")
("-0.902598308594d-17" "1.000000000000d+00")
("0.100000000000d+01" "0.000000000000d+01"))
(lisp-unit:fp-values (solve-cubic-complex -1.0d0 1.0d0 -1.0d0)))
(lisp-unit:assert-equal
'(("-0.809016994375d+00" "0.587785252292d+00")
("-0.809016994375d+00" "-0.587785252292d+00")
("0.309016994375d+00" "0.951056516295d+00")
("0.309016994375d+00" "-0.951056516295d+00")
("0.100000000000d+01" "0.000000000000d+01"))
;; Example from GSL manual
(lisp-unit:fp-values (polynomial-solve #(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
))