polynomial.lisp 6.56 KB
 liam committed Feb 04, 2008 1 2 ``````;; Polynomials ;; Liam Healy, Tue Mar 21 2006 - 18:33 `````` Liam Healy committed Dec 30, 2008 3 ``````;; Time-stamp: <2008-12-30 21:33:37EST polynomial.lisp> `````` lhealy committed Jul 24, 2008 4 ``````;; \$Id\$ `````` liam committed Mar 22, 2006 5 6 7 `````` (in-package :gsl) `````` 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 24, 2008 19 `````` (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) (x :double)) `````` Liam Healy committed Dec 30, 2008 20 21 `````` :inputs (coefficients) :c-return :double `````` liam committed Feb 04, 2008 22 `````` :documentation ; FDL `````` Liam Healy committed Dec 30, 2008 23 `````` "Evaluate the polyonomial with coefficients at the point x.") `````` liam committed May 30, 2006 24 `````` `````` liam committed Mar 22, 2006 25 26 27 28 ``````;;;;**************************************************************************** ;;;; Divided Difference Representation of Polynomials ;;;;**************************************************************************** `````` lhealy committed Jul 24, 2008 29 ``````(defmfun divided-difference (dd xa ya) `````` liam committed Jun 08, 2006 30 `````` "gsl_poly_dd_init" `````` lhealy committed Jul 24, 2008 31 32 33 34 35 36 37 38 39 40 `````` (((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 `````` Liam Healy committed Dec 30, 2008 41 `````` divided-differences of (xa,ya) are stored in the array `````` lhealy committed Jul 24, 2008 42 `````` dd, of the same length.") `````` liam committed Jun 08, 2006 43 `````` `````` liam committed Feb 18, 2008 44 ``````(defmfun polynomial-eval-divided-difference (dd xa x) `````` liam committed Jun 08, 2006 45 `````` "gsl_poly_dd_eval" `````` lhealy committed Jul 24, 2008 46 47 48 `````` (((c-pointer dd) :pointer) ((c-pointer xa) :pointer) ((dim0 xa) sizet) `````` liam committed Jun 08, 2006 49 `````` (x :double)) `````` lhealy committed Jul 24, 2008 50 `````` :inputs (dd xa) `````` liam committed Jun 08, 2006 51 `````` :c-return :double `````` liam committed Feb 04, 2008 52 `````` :documentation ; FDL `````` liam committed Mar 22, 2006 53 `````` "Evaluate the polynomial stored in divided-difference form `````` liam committed Feb 04, 2008 54 `````` in the arrays dd and xa at the point x.") `````` liam committed Jun 08, 2006 55 `````` `````` liam committed Feb 18, 2008 56 ``````(defmfun taylor-divided-difference (coefs xp dd xa workspace) `````` liam committed Jun 08, 2006 57 `````` "gsl_poly_dd_taylor" `````` lhealy committed Jul 24, 2008 58 `````` (((c-pointer coefs) :pointer) `````` liam committed Jun 14, 2006 59 `````` (xp :double) `````` lhealy committed Jul 24, 2008 60 61 62 63 64 65 `````` ((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 66 `````` :documentation ; FDL `````` liam committed Jun 14, 2006 67 68 `````` "Convert the divided-difference representation of a polynomial to a Taylor expansion. The divided-difference representation `````` liam committed Feb 04, 2008 69 `````` is supplied in the arrays dd and xa of the same length. `````` liam committed Jun 14, 2006 70 `````` On output the Taylor coefficients of the polynomial expanded about the `````` liam committed Feb 04, 2008 71 `````` point xp are stored in the array coefs which has the same length `````` lhealy committed Jul 24, 2008 72 `````` as xa and dd. A workspace of that length must be provided.") `````` liam committed Mar 22, 2006 73 74 75 76 77 `````` ;;;;**************************************************************************** ;;;; Quadratic Equations ;;;;**************************************************************************** `````` liam committed Feb 18, 2008 78 ``````(defmfun solve-quadratic (a b c) `````` liam committed Mar 22, 2006 79 `````` "gsl_poly_solve_quadratic" `````` liam committed Jun 08, 2006 80 `````` ((a :double) (b :double) (c :double) (root1 :double) (root2 :double)) `````` liam committed Feb 04, 2008 81 82 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 83 `````` "The real roots of the quadratic equation a x^2 + b x + c = 0. `````` liam committed Mar 22, 2006 84 `````` Two values are always returned; if the roots are not real, these `````` liam committed Feb 04, 2008 85 `````` values are NIL.") `````` liam committed Mar 22, 2006 86 `````` `````` liam committed Feb 18, 2008 87 ``````(defmfun solve-quadratic-complex (a b c) `````` liam committed Mar 22, 2006 88 `````` "gsl_poly_complex_solve_quadratic" `````` lhealy committed Jul 24, 2008 89 90 `````` ((a :double) (b :double) (c :double) (root1 complex-double-c) (root2 complex-double-c)) `````` liam committed Feb 04, 2008 91 92 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 93 94 `````` "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 95 `````` value returned will be NIL.") `````` liam committed Mar 22, 2006 96 97 98 99 100 `````` ;;;;**************************************************************************** ;;;; Cubic Equations ;;;;**************************************************************************** `````` liam committed Feb 18, 2008 101 ``````(defmfun solve-cubic (a b c) `````` liam committed Mar 22, 2006 102 `````` "gsl_poly_solve_cubic" `````` liam committed Jun 08, 2006 103 104 `````` ((a :double) (b :double) (c :double) (root1 :double) (root2 :double) (root3 :double)) `````` liam committed Feb 04, 2008 105 106 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 107 `````` "Find the real roots of the cubic equation, x^3 + a x^2 + b x + c = 0 `````` liam committed Mar 22, 2006 108 109 `````` with a leading coefficient of unity. The roots are given in ascending order. Three values are always returned; `````` liam committed Feb 04, 2008 110 `````` if a root is not real, the value returned for it will be NIL.") `````` liam committed Mar 22, 2006 111 `````` `````` liam committed Feb 18, 2008 112 ``````(defmfun solve-cubic-complex (a b c) `````` liam committed Mar 22, 2006 113 `````` "gsl_poly_complex_solve_cubic" `````` liam committed Jun 08, 2006 114 `````` ((a :double) (b :double) (c :double) `````` lhealy committed Jul 24, 2008 115 `````` (root1 complex-double-c) (root2 complex-double-c) (root3 complex-double-c)) `````` liam committed Feb 04, 2008 116 117 `````` :c-return :number-of-answers :documentation ; FDL `````` liam committed Mar 22, 2006 118 `````` "Find the complex roots of the cubic equation, x^3 + a x^2 + b x + c = 0 `````` liam committed Mar 22, 2006 119 `````` with a leading coefficient of unity. Three values are always returned; `````` liam committed Feb 04, 2008 120 `````` if a root does not exist, the value returned for it will be NIL.") `````` liam committed Mar 23, 2006 121 122 123 124 125 `````` ;;;;**************************************************************************** ;;;; General Polynomial Equations ;;;;**************************************************************************** `````` Liam Healy committed Dec 25, 2008 126 127 128 ``````(defmobject polynomial-complex-workspace "gsl_poly_complex_workspace" ((n sizet)) "complex workspace for polynomials") `````` liam committed Mar 26, 2006 129 `````` `````` Liam Healy committed Dec 30, 2008 130 131 132 133 134 135 ``````(defmfun polynomial-solve (coefficients &optional (answer (make-marray '(complex double-float) :dimensions (1- (total-size coefficients)))) (workspace (make-polynomial-complex-workspace (total-size coefficients)))) `````` liam committed Jan 28, 2007 136 `````` "gsl_poly_complex_solve" `````` lhealy committed Jul 24, 2008 137 `````` (((c-pointer coefficients) :pointer) ((dim0 coefficients) sizet) `````` Liam Healy committed Dec 30, 2008 138 `````` ((mpointer workspace) :pointer) ((c-pointer answer) :pointer)) `````` Liam Healy committed Dec 30, 2008 139 `````` :inputs (coefficients) `````` Liam Healy committed Dec 30, 2008 140 141 `````` :outputs (answer) :return (answer) `````` liam committed Feb 04, 2008 142 `````` :documentation ; FDL `````` Liam Healy committed Dec 30, 2008 143 144 145 146 147 148 149 `````` "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.") `````` liam committed Jan 28, 2007 150 `````` `````` liam committed Jun 08, 2006 151 152 153 154 ``````;;;;**************************************************************************** ;;;; Examples and unit test ;;;;**************************************************************************** `````` Liam Healy committed Oct 25, 2008 155 ``````(save-test polynomial `````` 156 157 `````` (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 158 `````` (dd (make-marray 'double-float :dimensions 4))) `````` lhealy committed Jul 24, 2008 159 160 161 162 163 164 `````` (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))) `````` 165 `````` (let ((vec #m(1.0d0 2.0d0 3.0d0))) `````` lhealy committed Jul 24, 2008 166 167 168 169 170 171 172 `````` (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 30, 2008 173 `````` (cl-array (polynomial-solve #m(-1.0d0 0.0d0 0.0d0 0.0d0 0.0d0 1.0d0)))) `````` Liam Healy committed Oct 25, 2008 174