polynomial.lisp 6.95 KB
 liam committed Feb 04, 2008 1 2 ``````;; Polynomials ;; Liam Healy, Tue Mar 21 2006 - 18:33 `````` Liam Healy committed Dec 27, 2008 3 ``````;; Time-stamp: <2008-12-26 18:53:14EST polynomial.lisp> `````` lhealy committed Jul 25, 2008 4 ``````;; \$Id\$ `````` liam committed Mar 22, 2006 5 6 7 `````` (in-package :gsl) `````` Liam Healy committed Dec 27, 2008 8 9 ``````;;; /usr/include/gsl/gsl_poly.h `````` liam committed Jun 14, 2006 10 11 ``````;;; Provide autotranslation from CL pure arrays? ;;; Divided differences not complete/tested. `````` liam committed Mar 22, 2006 12 13 14 15 16 `````` ;;;;**************************************************************************** ;;;; Polynomial Evaluation ;;;;**************************************************************************** `````` liam committed Feb 18, 2008 17 ``````(defmfun polynomial-eval (coefficients x) `````` liam committed Jun 08, 2006 18 `````` "gsl_poly_eval" `````` lhealy committed Jul 25, 2008 19 `````` (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) (x :double)) `````` liam committed Feb 04, 2008 20 `````` :documentation ; FDL `````` liam committed May 30, 2006 21 22 23 `````` "Evaluate the polyonomial with coefficients at the point x." :c-return :double) `````` liam committed Mar 22, 2006 24 25 26 27 ``````;;;;**************************************************************************** ;;;; Divided Difference Representation of Polynomials ;;;;**************************************************************************** `````` lhealy committed Jul 25, 2008 28 ``````(defmfun divided-difference (dd xa ya) `````` liam committed Jun 08, 2006 29 `````` "gsl_poly_dd_init" `````` lhealy committed Jul 25, 2008 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.") `````` liam committed Jun 08, 2006 42 `````` `````` liam committed Feb 18, 2008 43 ``````(defmfun polynomial-eval-divided-difference (dd xa x) `````` liam committed Jun 08, 2006 44 `````` "gsl_poly_dd_eval" `````` lhealy committed Jul 25, 2008 45 46 47 `````` (((c-pointer dd) :pointer) ((c-pointer xa) :pointer) ((dim0 xa) sizet) `````` liam committed Jun 08, 2006 48 `````` (x :double)) `````` lhealy committed Jul 25, 2008 49 `````` :inputs (dd xa) `````` liam committed Jun 08, 2006 50 `````` :c-return :double `````` liam committed Feb 04, 2008 51 `````` :documentation ; FDL `````` liam committed Mar 22, 2006 52 `````` "Evaluate the polynomial stored in divided-difference form `````` liam committed Feb 04, 2008 53 `````` in the arrays dd and xa at the point x.") `````` liam committed Jun 08, 2006 54 `````` `````` liam committed Feb 18, 2008 55 ``````(defmfun taylor-divided-difference (coefs xp dd xa workspace) `````` liam committed Jun 08, 2006 56 `````` "gsl_poly_dd_taylor" `````` lhealy committed Jul 25, 2008 57 `````` (((c-pointer coefs) :pointer) `````` liam committed Jun 14, 2006 58 `````` (xp :double) `````` lhealy committed Jul 25, 2008 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) `````` liam committed Feb 04, 2008 65 `````` :documentation ; FDL `````` liam committed Jun 14, 2006 66 67 `````` "Convert the divided-difference representation of a polynomial to a Taylor expansion. The divided-difference representation `````` liam committed Feb 04, 2008 68 `````` is supplied in the arrays dd and xa of the same length. `````` liam committed Jun 14, 2006 69 `````` On output the Taylor coefficients of the polynomial expanded about the `````` liam committed Feb 04, 2008 70 `````` point xp are stored in the array coefs which has the same length `````` lhealy committed Jul 25, 2008 71 `````` as xa and dd. A workspace of that length must be provided.") `````` liam committed Mar 22, 2006 72 73 74 75 76 `````` ;;;;**************************************************************************** ;;;; Quadratic Equations ;;;;**************************************************************************** `````` liam committed Feb 18, 2008 77 ``````(defmfun solve-quadratic (a b c) `````` liam committed Mar 22, 2006 78 `````` "gsl_poly_solve_quadratic" `````` liam committed Jun 08, 2006 79 `````` ((a :double) (b :double) (c :double) (root1 :double) (root2 :double)) `````` liam committed Feb 04, 2008 80 81 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 82 `````` "The real roots of the quadratic equation a x^2 + b x + c = 0. `````` liam committed Mar 22, 2006 83 `````` Two values are always returned; if the roots are not real, these `````` liam committed Feb 04, 2008 84 `````` values are NIL.") `````` liam committed Mar 22, 2006 85 `````` `````` liam committed Feb 18, 2008 86 ``````(defmfun solve-quadratic-complex (a b c) `````` liam committed Mar 22, 2006 87 `````` "gsl_poly_complex_solve_quadratic" `````` lhealy committed Jul 25, 2008 88 89 `````` ((a :double) (b :double) (c :double) (root1 complex-double-c) (root2 complex-double-c)) `````` liam committed Feb 04, 2008 90 91 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 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 `````` liam committed Feb 04, 2008 94 `````` value returned will be NIL.") `````` liam committed Mar 22, 2006 95 96 97 98 99 `````` ;;;;**************************************************************************** ;;;; Cubic Equations ;;;;**************************************************************************** `````` liam committed Feb 18, 2008 100 ``````(defmfun solve-cubic (a b c) `````` liam committed Mar 22, 2006 101 `````` "gsl_poly_solve_cubic" `````` liam committed Jun 08, 2006 102 103 `````` ((a :double) (b :double) (c :double) (root1 :double) (root2 :double) (root3 :double)) `````` liam committed Feb 04, 2008 104 105 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 106 `````` "Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0 `````` liam committed Mar 22, 2006 107 108 `````` with a leading coefficient of unity. The roots are given in ascending order. Three values are always returned; `````` liam committed Feb 04, 2008 109 `````` if a root is not real, the value returned for it will be NIL.") `````` liam committed Mar 22, 2006 110 `````` `````` liam committed Feb 18, 2008 111 ``````(defmfun solve-cubic-complex (a b c) `````` liam committed Mar 22, 2006 112 `````` "gsl_poly_complex_solve_cubic" `````` liam committed Jun 08, 2006 113 `````` ((a :double) (b :double) (c :double) `````` lhealy committed Jul 25, 2008 114 `````` (root1 complex-double-c) (root2 complex-double-c) (root3 complex-double-c)) `````` liam committed Feb 04, 2008 115 116 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 117 `````` "Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0 `````` liam committed Mar 22, 2006 118 `````` with a leading coefficient of unity. Three values are always returned; `````` liam committed Feb 04, 2008 119 `````` if a root does not exist, the value returned for it will be NIL.") `````` liam committed Mar 23, 2006 120 121 122 123 124 `````` ;;;;**************************************************************************** ;;;; General Polynomial Equations ;;;;**************************************************************************** `````` Liam Healy committed Dec 25, 2008 125 126 127 ``````(defmobject polynomial-complex-workspace "gsl_poly_complex_workspace" ((n sizet)) "complex workspace for polynomials") `````` liam committed Mar 26, 2006 128 `````` `````` lhealy committed Jul 25, 2008 129 ``````(export 'polynomial-solve) `````` liam committed Jun 08, 2006 130 ``````(defun polynomial-solve (coefficients) `````` liam committed Feb 04, 2008 131 `````` ;; FDL `````` liam committed Mar 23, 2006 132 `````` "The roots of the general polynomial `````` liam committed Feb 04, 2008 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 `````` liam committed Jun 08, 2006 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 `````` liam committed Feb 04, 2008 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 `````` liam committed Jun 08, 2006 139 `````` real and imaginary parts." `````` Liam Healy committed Dec 01, 2008 140 141 `````` (let ((len (total-size coefficients))) ;; Should this be making a complex array? `````` Liam Healy committed Dec 26, 2008 142 143 `````` (let ((answer (make-marray 'double-float :dimensions (* 2 (1- len)))) (ws (make-polynomial-complex-workspace len))) `````` Liam Healy committed Dec 01, 2008 144 `````` (values-list (polynomial-solve-ws coefficients ws answer))))) `````` liam committed Jun 08, 2006 145 `````` `````` liam committed Feb 18, 2008 146 ``````(defmfun polynomial-solve-ws (coefficients workspace answer-pd) `````` liam committed Jan 28, 2007 147 `````` "gsl_poly_complex_solve" `````` lhealy committed Jul 25, 2008 148 `````` (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) `````` Liam Healy committed Dec 26, 2008 149 `````` ((mpointer workspace) :pointer) ((c-pointer answer-pd) :pointer)) `````` liam committed Jan 28, 2007 150 151 `````` :return ((loop for i from 0 below (dim0 answer-pd) by 2 `````` liam committed Feb 24, 2008 152 153 `````` collect (complex (maref answer-pd i) (maref answer-pd (1+ i))))) `````` lhealy committed Jul 25, 2008 154 155 `````` :export nil :index polynomial-solve `````` liam committed Feb 04, 2008 156 `````` :documentation ; FDL `````` liam committed Jan 28, 2007 157 `````` "Arguments are: `````` lhealy committed Jul 25, 2008 158 `````` a GSL array of coefficients, a workspace, a gsl-array of doubles.") `````` liam committed Jan 28, 2007 159 160 `````` `````` liam committed Jun 08, 2006 161 162 163 164 ``````;;;;**************************************************************************** ;;;; Examples and unit test ;;;;**************************************************************************** `````` Liam Healy committed Oct 25, 2008 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)) `````` Liam Healy committed Dec 26, 2008 168 `````` (dd (make-marray 'double-float :dimensions 4))) `````` lhealy committed Jul 25, 2008 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))) `````` lhealy committed Jul 25, 2008 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 `````` Liam Healy committed Dec 01, 2008 183 `````` (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0))) `````` Liam Healy committed Oct 25, 2008 184