polynomial.lisp 11.5 KB
Newer Older
1 2
;; Polynomials
;; Liam Healy, Tue Mar 21 2006 - 18:33
3
;; Time-stamp: <2016-06-15 22:31:17EDT polynomial.lisp>
4
;;
5
;; Copyright 2009, 2011, 2012, 2014, 2016 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21

(in-package :gsl)
22
(named-readtables:in-readtable :antik)
23

24 25
;;; /usr/include/gsl/gsl_poly.h

26 27 28 29
;;;;****************************************************************************
;;;; Polynomial Evaluation
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
30
(defmfun evaluate
31
    ((coefficients grid:vector-double-float) (x float) &key divided-difference)
Liam Healy's avatar
Liam Healy committed
32
  ("gsl_poly_eval" "gsl_poly_dd_eval")
33
  ((((grid:foreign-pointer coefficients) :pointer) ((dim0 coefficients) :sizet)
Liam Healy's avatar
Liam Healy committed
34
    (x :double))
35 36
   (((grid:foreign-pointer divided-difference) :pointer)
    ((grid:foreign-pointer coefficients) :pointer)
37
    ((dim0 coefficients) :sizet)
Liam Healy's avatar
Liam Healy committed
38 39 40
    (x :double)))
  :definition :method
  :inputs (coefficients divided-difference)
41
  :c-return :double
42
  :documentation			; FDL
43
  "Evaluate the polyonomial with coefficients at the point x.")
44

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
(defmfun evaluate-with-derivatives
    (coefficients x &optional (derivatives 2)
		  &aux (result (grid:ensure-foreign-array derivatives derivatives 'double-float)))
  "gsl_poly_eval_derivs"
  (((grid:foreign-pointer coefficients) :pointer)
   ((grid:dim0 coefficients) :sizet)
   (x :double)
   ((grid:foreign-pointer result) :pointer)
   ((grid:dim0 result) :sizet))
  :gsl-version (1 13)
  :inputs (coefficients x)
  :return (result)
  :documentation
  "Evaluates a polynomial and its derivatives and stores the results in the array @var{res} of size @var{lenres}.  The output array contains the values of @math{d^k P/d x^k} for the specified value of @var{x} starting with @math{k = 0}. The optional argument 'derivatives may be either a vector-double-float, or a non-negative integer. If the former, the function value and derivatives are put in the vector supplied; if the latter, a new vector-double-float is created with the specified length.")

Liam Healy's avatar
Liam Healy committed
60
#+fsbv
Liam Healy's avatar
Liam Healy committed
61
(defmfun evaluate
62
    ((coefficients grid:vector-double-float) (x complex)
63
     &key)
Liam Healy's avatar
Liam Healy committed
64
  "gsl_poly_complex_eval"
65
  (((grid:foreign-pointer coefficients) :pointer) ((dim0 coefficients) :sizet)
66
   (x (:struct grid:complex-double-c)))
Liam Healy's avatar
Liam Healy committed
67 68 69
  :definition :method
  :gsl-version (1 11)
  :inputs (coefficients)
70
  :c-return (:struct grid:complex-double-c)
Liam Healy's avatar
Liam Healy committed
71 72 73
  :documentation			; FDL
  "Evaluate the polyonomial with coefficients at the complex value x.")

Liam Healy's avatar
Liam Healy committed
74
#+fsbv
Liam Healy's avatar
Liam Healy committed
75
(defmfun evaluate
76
    ((coefficients grid:vector-complex-double-float) (x complex)
77
     &key)
Liam Healy's avatar
Liam Healy committed
78
  "gsl_complex_poly_complex_eval"
79
  (((grid:foreign-pointer coefficients) :pointer) ((dim0 coefficients) :sizet)
80
   (x (:struct grid:complex-double-c)))
Liam Healy's avatar
Liam Healy committed
81 82 83
  :definition :method
  :gsl-version (1 11)
  :inputs (coefficients)
84
  :c-return (:struct grid:complex-double-c)
Liam Healy's avatar
Liam Healy committed
85 86 87
  :documentation			; FDL
  "Evaluate the polyonomial with coefficients at the complex value x.")

88 89 90 91
;;;;****************************************************************************
;;;; Divided Difference Representation of Polynomials
;;;;****************************************************************************

92
(defmfun divided-difference
93
    (xa ya &optional (dd (grid:make-foreign-array 'double-float :dimensions (dim0 xa))))
94
  "gsl_poly_dd_init"
95 96
  (((grid:foreign-pointer dd) :pointer)
   ((grid:foreign-pointer xa) :pointer) ((grid:foreign-pointer ya) :pointer)
97
   ((dim0 xa) :sizet))
98 99 100 101 102 103 104
  :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
105
   divided-differences of (xa,ya) are stored in the array
106
   dd, of the same length.")
107

108 109 110
(defmfun taylor-divided-difference
    (xp dd xa
	&optional
111 112
	(coefficients (grid:make-foreign-array 'double-float :dimensions (dim0 xa)))
	(workspace (grid:make-foreign-array 'double-float :dimensions (dim0 xa))))
113
  "gsl_poly_dd_taylor"
114
  (((grid:foreign-pointer coefficients) :pointer)
115
   (xp :double)
116 117
   ((grid:foreign-pointer dd) :pointer)
   ((grid:foreign-pointer xa) :pointer)
118
   ((dim0 xa) :sizet)
119
   ((grid:foreign-pointer workspace) :pointer))
120 121
  :inputs (dd xa)
  :outputs (coefficients)
122
  :documentation			; FDL
123 124
  "Convert the divided-difference representation of a
  polynomial to a Taylor expansion.  The divided-difference representation
125
  is supplied in the arrays dd and xa of the same length.
126
  On output the Taylor coefficients of the polynomial expanded about the
127 128
  point xp are stored in the array coefficients which has the same length
  as xa and dd.")
129 130 131 132 133

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

liam's avatar
liam committed
134
(defmfun solve-quadratic (a b c)
135
  "gsl_poly_solve_quadratic"
136 137
  ((a :double) (b :double) (c :double)
   (root1 (:pointer :double)) (root2 (:pointer :double)))
138 139
  :c-return :number-of-answers
  :documentation			; FDL
140
  "The real roots of the quadratic equation a x^2 + b x + c = 0.
141
   Two values are always returned; if the roots are not real, these
142
   values are NIL.")
143

liam's avatar
liam committed
144
(defmfun solve-quadratic-complex (a b c)
145
  "gsl_poly_complex_solve_quadratic"
146
  ((a :double) (b :double) (c :double)
147 148
   (root1 (:pointer (:struct grid:complex-double-c)))
   (root2 (:pointer (:struct grid:complex-double-c))))
149 150
  :c-return :number-of-answers
  :documentation			; FDL
151 152
  "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
153
   value returned will be NIL.") 
154 155 156 157 158

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

liam's avatar
liam committed
159
(defmfun solve-cubic (a b c)
160
  "gsl_poly_solve_cubic"
161
  ((a :double) (b :double) (c :double)
162 163
   (root1 (:pointer :double))
   (root2 (:pointer :double)) (root3 (:pointer :double)))
164 165
  :c-return :number-of-answers
  :documentation			; FDL
166
  "Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0
167 168
   with a leading coefficient of unity.  The roots are given
   in ascending order.  Three values are always returned;
169
   if a root is not real, the value returned for it will be NIL.")
170

liam's avatar
liam committed
171
(defmfun solve-cubic-complex (a b c)
172
  "gsl_poly_complex_solve_cubic"
173
  ((a :double) (b :double) (c :double)
174 175 176
   (root1 (:pointer (:struct grid:complex-double-c)))
   (root2 (:pointer (:struct grid:complex-double-c)))
   (root3 (:pointer (:struct grid:complex-double-c))))
177 178
  :c-return :number-of-answers
  :documentation			; FDL
179
  "Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0
180
   with a leading coefficient of unity.  Three values are always returned;
181
   if a root does not exist, the value returned for it will be NIL.")
182 183 184 185 186

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

187
(defmobject polynomial-complex-workspace "gsl_poly_complex_workspace"
188
  ((n :sizet))
189
  "complex workspace for polynomials")
190

191 192 193
(defmfun polynomial-solve
    (coefficients
     &optional
194
     (answer (grid:make-foreign-array '(complex double-float)
195 196
			  :dimensions (1- (size coefficients))))
     (workspace (make-polynomial-complex-workspace (size coefficients))))
197
  "gsl_poly_complex_solve"
198
  (((grid:foreign-pointer coefficients) :pointer) ((dim0 coefficients) :sizet)
199
   ((mpointer workspace) :pointer) ((grid:foreign-pointer answer) :pointer))
200
  :inputs (coefficients)
201 202
  :outputs (answer)
  :return (answer)
203
  :documentation			; FDL
204 205 206 207 208 209 210
  "Arguments are: a vector-double-float of coefficients, a complex
   vector of length one less than coefficients that will hold the
   answer, and a workspace made by make-polynomial-complex-workspace.
   The roots of the general polynomial 
   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 coefficient of the
   highest order term must be non-zero.")
211

212 213 214 215
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
216
(save-test polynomial
217 218
 (let* ((xa #m(0.0d0 1.0d0 2.0d0 3.0d0))
	(dd (divided-difference xa #m(2.5d0 7.2d0 32.7d0 91.0d0))))
219
   (list
Liam Healy's avatar
Liam Healy committed
220 221 222 223
    (evaluate xa 0.0d0 :divided-difference dd)
    (evaluate xa 1.0d0 :divided-difference dd)
    (evaluate xa 2.0d0 :divided-difference dd)
    (evaluate xa 3.0d0 :divided-difference dd)))
224
 (let ((vec #m(1.0d0 2.0d0 3.0d0)))
Liam Healy's avatar
Liam Healy committed
225
   (evaluate vec -1.0d0))
226 227 228 229 230 231
 (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
232
 (grid:copy-to (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
233
 ;; tests from gsl-1.11/poly/test.c
234 235 236 237
 (evaluate #m(1.0d0 0.5d0 0.3d0) 0.5d0)
 (evaluate
  #m(1.0d0 -1.0d0 1.0d0 -1.0d0 1.0d0 -1.0d0 1.0d0 -1.0d0 1.0d0 -1.0d0 1.0d0)
  1.0d0)
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
 (solve-quadratic 4.0d0 -20.0d0 26.0d0)		  ; no roots
 (solve-quadratic 4.0d0 -20.0d0 25.0d0)		  ; one root
 (solve-quadratic 4.0d0 -20.0d0 21.0d0)		  ; two roots
 (solve-quadratic 4.0d0 7.0d0 0.0d0)		  ; two roots
 (solve-quadratic 5.0d0 0.0d0 -20.0d0)		  ; two roots
 (solve-quadratic 0.0d0 3.0d0 -21.0d0)		  ; one root (linear)
 (solve-quadratic 0.0d0 0.0d0 1.0d0)
 (solve-cubic 0.0d0 0.0d0 -27.0d0)		      ; one root
 (solve-cubic -51.0d0 867.0d0 -4913.0d0)	      ; three roots
 (solve-cubic -57.0d0 1071.0d0 -6647.0d0)	      ; three roots
 (solve-cubic -11.0d0 -493.0d0 +6647.0d0)	      ; three roots
 (solve-cubic -143.0d0 5087.0d0 -50065.0d0)	      ; three roots
 (solve-cubic -109.0d0 803.0d0 50065.0d0)	      ; three roots
 (solve-quadratic-complex 4.0d0 -20.0d0 26.0d0)
 (solve-quadratic-complex 4.0d0 -20.0d0 25.0d0)
 (solve-quadratic-complex 4.0d0 -20.0d0 21.0d0)
 (solve-quadratic-complex 4.0d0 7.0d0 0.0d0)
 (solve-quadratic-complex 5.0d0 0.0d0 -20.0d0)
 (solve-quadratic-complex 5.0d0 0.0d0 20.0d0)
 (solve-quadratic-complex 0.0d0 3.0d0 -21.0d0)
 (solve-quadratic-complex 0.0d0 0.0d0 1.0d0)
 (solve-cubic-complex 0.0d0 0.0d0 -27.0d0)
 (solve-cubic-complex -1.0d0 1.0d0 39.0d0)
 (solve-cubic-complex -51.0d0 867.0d0 -4913.0d0)
 (solve-cubic-complex -57.0d0 1071.0d0 -6647.0d0)
 (solve-cubic-complex -11.0d0 -493.0d0 +6647.0d0)
 (solve-cubic-complex -143.0d0 5087.0d0 -50065.0d0)
265 266 267 268
 (grid:copy-to (polynomial-solve #m(-120.0d0 274.0d0 -225.0d0 85.0d0 -15.0d0 1.0d0)))
 (grid:copy-to (polynomial-solve #m(1.0d0 0.0d0 0.0d0 0.0d0 1.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))
 (let* ((xa #m(0.16d0 0.97d0 1.94d0 2.74d0 3.58d0 3.73d0 4.70d0))
	(ya #m(0.73d0 1.11d0 1.49d0 1.84d0 2.30d0 2.41d0 3.07d0))
269 270
	(dd (divided-difference xa ya)))
   (list
271
    (grid:copy-to dd)
272
    (map 'vector (lambda (x) (evaluate xa x :divided-difference dd))
273
	 (grid:copy-to xa))
274 275
    (map 'vector
	 (lambda (x) (evaluate (taylor-divided-difference 1.5d0 dd xa) (- x 1.5d0)))
276
	 (grid:copy-to xa)))))