polynomial.lisp 6.95 KB
Newer Older
1 2
;; Polynomials
;; Liam Healy, Tue Mar 21 2006 - 18:33
3
;; Time-stamp: <2008-12-26 18:53:14EST polynomial.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8 9
;;; /usr/include/gsl/gsl_poly.h

10 11
;;; Provide autotranslation from CL pure arrays?
;;; Divided differences not complete/tested.
12 13 14 15 16

;;;;****************************************************************************
;;;; Polynomial Evaluation
;;;;****************************************************************************

liam's avatar
liam committed
17
(defmfun polynomial-eval (coefficients x)
18
  "gsl_poly_eval"
19
  (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) (x :double))
20
  :documentation			; FDL
21 22 23
  "Evaluate the polyonomial with coefficients at the point x."
  :c-return :double)

24 25 26 27
;;;;****************************************************************************
;;;; Divided Difference Representation of Polynomials
;;;;****************************************************************************

28
(defmfun divided-difference (dd xa ya)
29
  "gsl_poly_dd_init"
30 31 32 33 34 35 36 37 38 39 40 41
  (((c-pointer dd) :pointer)
   ((c-pointer xa) :pointer) ((c-pointer ya) :pointer)
   ((dim0 xa) sizet))
  :inputs (xa ya)
  :outputs (dd)
  :return (dd)
  :documentation			; FDL
  "Compute a divided-difference representation of the
   interpolating polynomial for the points (xa, ya) stored in
   the arrays of equal length.  On output the
   divided-differences of (@var{xa},@var{ya}) are stored in the array
   dd, of the same length.")
42

liam's avatar
liam committed
43
(defmfun polynomial-eval-divided-difference (dd xa x)
44
  "gsl_poly_dd_eval"
45 46 47
  (((c-pointer dd) :pointer)
   ((c-pointer xa) :pointer)
   ((dim0 xa) sizet)
48
   (x :double))
49
  :inputs (dd xa)
50
  :c-return :double
51
  :documentation			; FDL
52
  "Evaluate the polynomial stored in divided-difference form
53
   in the arrays dd and xa at the point x.")
54

liam's avatar
liam committed
55
(defmfun taylor-divided-difference (coefs xp dd xa workspace)
56
  "gsl_poly_dd_taylor"
57
  (((c-pointer coefs) :pointer)
58
   (xp :double)
59 60 61 62 63 64
   ((c-pointer dd) :pointer)
   ((c-pointer xa) :pointer)
   ((dim0 xa) sizet)
   ((c-pointer workspace) :pointer))
  :inputs (coefs xa)
  :outputs (coefs)
65
  :documentation			; FDL
66 67
  "Convert the divided-difference representation of a
  polynomial to a Taylor expansion.  The divided-difference representation
68
  is supplied in the arrays dd and xa of the same length.
69
  On output the Taylor coefficients of the polynomial expanded about the
70
  point xp are stored in the array coefs which has the same length
71
  as xa and dd.  A workspace of that length must be provided.")
72 73 74 75 76

;;;;****************************************************************************
;;;; Quadratic Equations
;;;;****************************************************************************

liam's avatar
liam committed
77
(defmfun solve-quadratic (a b c)
78
  "gsl_poly_solve_quadratic"
79
  ((a :double) (b :double) (c :double) (root1 :double) (root2 :double))
80 81
  :c-return :number-of-answers
  :documentation			; FDL
82
  "The real roots of the quadratic equation a x^2 + b x + c = 0.
83
   Two values are always returned; if the roots are not real, these
84
   values are NIL.")
85

liam's avatar
liam committed
86
(defmfun solve-quadratic-complex (a b c)
87
  "gsl_poly_complex_solve_quadratic"
88 89
  ((a :double) (b :double) (c :double)
   (root1 complex-double-c) (root2 complex-double-c))
90 91
  :c-return :number-of-answers
  :documentation			; FDL
92 93
  "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
94
   value returned will be NIL.") 
95 96 97 98 99

;;;;****************************************************************************
;;;; Cubic Equations
;;;;****************************************************************************

liam's avatar
liam committed
100
(defmfun solve-cubic (a b c)
101
  "gsl_poly_solve_cubic"
102 103
  ((a :double) (b :double) (c :double)
   (root1 :double) (root2 :double) (root3 :double))
104 105
  :c-return :number-of-answers
  :documentation			; FDL
106
  "Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0
107 108
   with a leading coefficient of unity.  The roots are given
   in ascending order.  Three values are always returned;
109
   if a root is not real, the value returned for it will be NIL.")
110

liam's avatar
liam committed
111
(defmfun solve-cubic-complex (a b c)
112
  "gsl_poly_complex_solve_cubic"
113
  ((a :double) (b :double) (c :double)
114
   (root1 complex-double-c) (root2 complex-double-c) (root3 complex-double-c))
115 116
  :c-return :number-of-answers
  :documentation			; FDL
117
  "Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0
118
   with a leading coefficient of unity.  Three values are always returned;
119
   if a root does not exist, the value returned for it will be NIL.")
120 121 122 123 124

;;;;****************************************************************************
;;;; General Polynomial Equations
;;;;****************************************************************************

125 126 127
(defmobject polynomial-complex-workspace "gsl_poly_complex_workspace"
  ((n sizet))
  "complex workspace for polynomials")
128

129
(export 'polynomial-solve)
130
(defun polynomial-solve (coefficients)
131
  ;; FDL
132
  "The roots of the general polynomial 
133 134
  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 n
135 136
  specifies the length of the coefficient array.  The coefficient of the
  highest order term must be non-zero.  The function requires a workspace
137 138
  w of the appropriate size.  The n-1 roots are returned in
  the packed complex array z of length 2(n-1), alternating
139
  real and imaginary parts."
140 141
  (let ((len (total-size coefficients)))
    ;; Should this be making a complex array?
142 143
    (let ((answer (make-marray 'double-float :dimensions (* 2 (1- len))))
	   (ws (make-polynomial-complex-workspace len)))
144
      (values-list (polynomial-solve-ws coefficients ws answer)))))
145

liam's avatar
liam committed
146
(defmfun polynomial-solve-ws (coefficients workspace answer-pd)
147
  "gsl_poly_complex_solve"
148
  (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet)
149
   ((mpointer workspace) :pointer) ((c-pointer answer-pd) :pointer))
150 151
  :return
  ((loop for i from 0 below (dim0 answer-pd) by 2
152 153
	 collect (complex (maref answer-pd i)
			  (maref answer-pd (1+ i)))))
154 155
  :export nil
  :index polynomial-solve
156
  :documentation			; FDL
157
  "Arguments are:
158
   a GSL array of coefficients, a workspace, a gsl-array of doubles.")
159 160


161 162 163 164
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
165
(save-test polynomial
166 167
 (let ((xa #m(0.0d0 1.0d0 2.0d0 3.0d0))
	(ya #m(2.5d0 7.2d0 32.7d0 91.0d0))
168
	(dd (make-marray 'double-float :dimensions 4)))
169 170 171 172 173 174
   (divided-difference dd xa ya)
   (list
    (polynomial-eval-divided-difference dd xa 0.0d0)
    (polynomial-eval-divided-difference dd xa 1.0d0)
    (polynomial-eval-divided-difference dd xa 2.0d0)
    (polynomial-eval-divided-difference dd xa 3.0d0)))
175
 (let ((vec #m(1.0d0 2.0d0 3.0d0)))
176 177 178 179 180 181 182
   (polynomial-eval vec -1.0d0))
 (solve-quadratic 1.0d0 0.0d0 1.0d0)
 (solve-quadratic 1.0d0 -2.0d0 1.0d0)
 (solve-quadratic-complex 1.0d0 -2.0d0 1.0d0)
 (solve-cubic -6.0d0 -13.0d0 42.0d0)
 (solve-cubic-complex -1.0d0 1.0d0 -1.0d0)
 ;; Example from GSL manual
183
 (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
Liam Healy's avatar
Liam Healy committed
184