polynomial.lisp 7.2 KB
Newer Older
1 2
;; Polynomials
;; Liam Healy, Tue Mar 21 2006 - 18:33
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2008-10-25 11:47:01EDT polynomial.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

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

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

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

22 23 24 25
;;;;****************************************************************************
;;;; Divided Difference Representation of Polynomials
;;;;****************************************************************************

26
(defmfun divided-difference (dd xa ya)
27
  "gsl_poly_dd_init"
28 29 30 31 32 33 34 35 36 37 38 39
  (((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.")
40

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

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

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

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

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

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

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

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

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

123
(defgo-s (complex-workspace n)
124
	 allocate-complex-workspace free-complex-workspace)
125

liam's avatar
liam committed
126
(defmfun allocate-complex-workspace (n)
127
  "gsl_poly_complex_workspace_alloc" ((n sizet))
128 129
  :c-return :pointer
  :export nil
130
  :index (letm complex-workspace))
131

liam's avatar
liam committed
132
(defmfun free-complex-workspace (ws)
133 134 135
  "gsl_poly_complex_workspace_free" ((ws :pointer))
  :c-return :void
  :export nil
136
  :index (letm complex-workspace))
137

138
(export 'polynomial-solve)
139
(defun polynomial-solve (coefficients)
140
  ;; FDL
141
  "The roots of the general polynomial 
142 143
  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
144 145
  specifies the length of the coefficient array.  The coefficient of the
  highest order term must be non-zero.  The function requires a workspace
146 147
  w of the appropriate size.  The n-1 roots are returned in
  the packed complex array z of length 2(n-1), alternating
148 149
  real and imaginary parts."
  (let ((len (length coefficients)))
150 151
    (letm ((coef (vector-double-float coefficients))
	   (answer (vector-double-float (* 2 (1- len))))
152 153
	   (ws (complex-workspace len)))
      (values-list (polynomial-solve-ws coef ws answer)))))
154

liam's avatar
liam committed
155
(defmfun polynomial-solve-ws (coefficients workspace answer-pd)
156
  "gsl_poly_complex_solve"
157 158
  (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet)
   (workspace :pointer) ((c-pointer answer-pd) :pointer))
159 160
  :return
  ((loop for i from 0 below (dim0 answer-pd) by 2
161 162
	 collect (complex (maref answer-pd i)
			  (maref answer-pd (1+ i)))))
163 164
  :export nil
  :index polynomial-solve
165
  :documentation			; FDL
166
  "Arguments are:
167
   a GSL array of coefficients, a workspace, a gsl-array of doubles.")
168 169


170 171 172 173
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
174
(save-test polynomial
Liam Healy's avatar
Liam Healy committed
175 176
 (letm ((xa (vector-double-float (a 0.0d0 1.0d0 2.0d0 3.0d0)))
	(ya (vector-double-float (a 2.5d0 7.2d0 32.7d0 91.0d0)))
177 178 179 180 181 182 183
	(dd (vector-double-float 4)))
   (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)))
Liam Healy's avatar
Liam Healy committed
184
 (letm ((vec (vector-double-float (a 1.0d0 2.0d0 3.0d0))))
185 186 187 188 189 190 191
   (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
Liam Healy's avatar
Liam Healy committed
192
 (polynomial-solve (a double-float -1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
Liam Healy's avatar
Liam Healy committed
193