[1] | 1 | ;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
|
---|
| 2 | #|
|
---|
| 3 | $Id: dynamics.lisp,v 1.7 2009/01/23 10:49:32 marek Exp $
|
---|
| 4 | *--------------------------------------------------------------------------*
|
---|
| 5 | | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu) |
|
---|
| 6 | | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
|
---|
| 7 | | |
|
---|
| 8 | | Everyone is permitted to copy, distribute and modify the code in this |
|
---|
| 9 | | directory, as long as this copyright note is preserved verbatim. |
|
---|
| 10 | *--------------------------------------------------------------------------*
|
---|
| 11 | |#
|
---|
| 12 |
|
---|
| 13 | (defpackage "DYNAMICS"
|
---|
| 14 | (:use "ORDER" "MONOM" "COEFFICIENT-RING" "GROBNER" "MAKELIST" "PRINTER" "TERM" "POLY" "COMMON-LISP")
|
---|
| 15 | (:export poly-composition poly-scalar-composition poly-dynamic-power
|
---|
| 16 | poly-evaluate poly-scalar-evaluate factorial
|
---|
| 17 | poly-scalar-diff poly-diff standard-vector
|
---|
| 18 | scalar-partial partial drop-elt drop-row drop-column
|
---|
| 19 | determinant minor matrix- poly-list-
|
---|
| 20 | characteristic-combination
|
---|
| 21 | characteristic-matrix
|
---|
| 22 | characteristic-polynomial
|
---|
| 23 | identity-matrix
|
---|
| 24 | print-matrix
|
---|
| 25 | jacobi-matrix
|
---|
| 26 | jacobian))
|
---|
| 27 |
|
---|
| 28 | (in-package "DYNAMICS")
|
---|
| 29 |
|
---|
[92] | 30 | (proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 3)))
|
---|
[1] | 31 |
|
---|
| 32 | (defun poly-scalar-composition (f G &optional (order #'lex>))
|
---|
| 33 | "Returns a polynomial obtained by substituting a list of polynomials G=(G1,G2,...,GN)
|
---|
| 34 | into a polynomial F(X1,X2,...,XN). All polynomials are assumed to be in the internal form,
|
---|
| 35 | so variables do not explicitly apprear in the calculation."
|
---|
| 36 | (reduce #'(lambda (x y) (poly+ x y order))
|
---|
| 37 | (mapcar #'(lambda (x)
|
---|
| 38 | (scalar-times-poly
|
---|
| 39 | (cdr x)
|
---|
| 40 | (poly-mexpt G (car x) order)))
|
---|
| 41 | f)))
|
---|
| 42 |
|
---|
| 43 | (defun poly-composition (F G &optional (order #'lex>))
|
---|
| 44 | "Return the composition of a polynomial map F with a polynomial map
|
---|
| 45 | G. Both maps are represented as lists of polynomials, and each polynomial
|
---|
| 46 | is in the internal alist representation. The restriction is that the length
|
---|
| 47 | of the list G must be the number of variables in the list F."
|
---|
| 48 | (mapcar #'(lambda (x) (poly-scalar-composition x G order)) F))
|
---|
| 49 |
|
---|
| 50 |
|
---|
| 51 | (defun poly-dynamic-power (F n &optional (order #'lex>))
|
---|
| 52 | "Calculate the composition FoFo...oF (n times), where
|
---|
| 53 | F is a polynomial map represented as a list of polynomials."
|
---|
| 54 | (reduce #'(lambda (x y) (poly-composition x y order))
|
---|
| 55 | (make-list n :initial-element F)))
|
---|
| 56 |
|
---|
| 57 | (defun poly-scalar-evaluate (f x &optional (order #'lex>))
|
---|
| 58 | "Evaluate a polynomial F at a point X. This operation is implemented
|
---|
| 59 | through POLY-SCALAR-COMPOSITION."
|
---|
| 60 | (if (endp f)
|
---|
| 61 | 0
|
---|
| 62 | (let ((p (poly-scalar-composition
|
---|
| 63 | f
|
---|
| 64 | (mapcar #'(lambda (u)
|
---|
| 65 | (if (zerop u)
|
---|
| 66 | nil
|
---|
| 67 | (list (cons (make-list (length (caar f)) :initial-element 0)
|
---|
| 68 | u))))
|
---|
| 69 | x)
|
---|
| 70 | order)))
|
---|
| 71 | (if (endp p) 0
|
---|
| 72 | (cdar p)))))
|
---|
| 73 |
|
---|
| 74 | (defun poly-evaluate (F x &optional (order #'lex>))
|
---|
| 75 | "Evaluate a polynomial map F, represented as list of polynomials, at a point X."
|
---|
| 76 | (mapcar #'(lambda (h) (poly-scalar-evaluate h x order)) F))
|
---|
| 77 |
|
---|
| 78 | (defun factorial (n &optional (k n) &aux (result 1))
|
---|
| 79 | "Return N!/(N-K)!=N(N-1)(N-K+1)."
|
---|
| 80 | (dotimes (j k result)
|
---|
| 81 | (setf result (* result (- n j)))))
|
---|
| 82 |
|
---|
| 83 | (defun poly-scalar-diff (f m)
|
---|
| 84 | "Return the partial derivative of a polynomial F over multiple
|
---|
| 85 | variables according to multiindex M."
|
---|
| 86 | (setf f (remove-if #'(lambda (x) (some #'< (car x) m)) f))
|
---|
| 87 | (mapcar #'(lambda (x) (cons (monom/ (car x) m)
|
---|
| 88 | (* (cdr x) (apply #'* (mapcar #'factorial (car x) m)))))
|
---|
| 89 | f))
|
---|
| 90 |
|
---|
| 91 | (defun poly-diff (F m)
|
---|
| 92 | "Return the partial derivative of a polynomial map F, represented as a list of polynomials,
|
---|
| 93 | with respect to several variables, according to multi-index M."
|
---|
| 94 | (mapcar #'(lambda (h) (poly-scalar-diff h m)) F))
|
---|
| 95 |
|
---|
| 96 |
|
---|
| 97 | (defun scalar-partial (f k &optional (l 1))
|
---|
| 98 | "Returns the L-th partial derivative of a polynomial F over the K-th variable."
|
---|
| 99 | (when f
|
---|
| 100 | (poly-scalar-diff f (standard-vector (length (caar f)) k l))))
|
---|
| 101 |
|
---|
| 102 | (defun partial (F k &optional (l 1))
|
---|
| 103 | "Returns the L-th partial derivative over the K-th variable, of a polynomial
|
---|
| 104 | map F, represented as a list of polynomials."
|
---|
| 105 | (when (and F (car F) (caar F))
|
---|
| 106 | (poly-diff F (standard-vector (length (caaar F)) k l))))
|
---|
| 107 |
|
---|
| 108 |
|
---|
| 109 | (defun determinant (F &optional (order #'lex>) &aux (result nil))
|
---|
| 110 | "Returns the determinant of a polynomial matrix F, which is a list of
|
---|
| 111 | rows of the matrix, each row is a list of polynomials. The algorithm
|
---|
| 112 | is recursive expansion along columns."
|
---|
| 113 | (cond
|
---|
| 114 | ((= (length F) 1) (setf result (caar F)))
|
---|
| 115 | (t
|
---|
| 116 | (dotimes (i (length F) result)
|
---|
| 117 | (setf result (poly+ result
|
---|
| 118 | (scalar-times-poly
|
---|
| 119 | (expt -1 i)
|
---|
| 120 | (poly* (car (elt F i)) (minor F i 0 order) order))
|
---|
| 121 | order))))))
|
---|
| 122 |
|
---|
| 123 |
|
---|
| 124 | (defun minor (F i j &optional (order #'lex>))
|
---|
| 125 | "Calculate the minor of a polynomial matrix F with respect to entry (I,J)."
|
---|
| 126 | (determinant (drop-row (drop-column F j) i) order))
|
---|
| 127 |
|
---|
| 128 | (defun drop-row (F i)
|
---|
| 129 | "Discards the I-th row from a polynomial matrix F."
|
---|
| 130 | (drop-elt F i))
|
---|
| 131 |
|
---|
| 132 | (defun drop-column (F j)
|
---|
| 133 | "Discards the J-th column from a polynomial matrix F."
|
---|
| 134 | (mapcar #'(lambda (x) (drop-elt x j)) F))
|
---|
| 135 |
|
---|
| 136 | (defun drop-elt (lst j)
|
---|
| 137 | "Discards the J-th element from a list LST."
|
---|
| 138 | (append (subseq lst 0 j) (subseq lst (1+ j))))
|
---|
| 139 |
|
---|
| 140 | ;; Matrix operations
|
---|
| 141 |
|
---|
| 142 | (defun matrix- (F G &optional (order #'lex>))
|
---|
| 143 | "Returns difference of two polynomial matrices F and G."
|
---|
| 144 | (mapcar #'(lambda (x y) (poly-list- x y order)) F G))
|
---|
| 145 |
|
---|
| 146 | (defun scalar-times-matrix (s F)
|
---|
| 147 | "Returns a product of a polynomial S by a polynomial matrix F."
|
---|
| 148 | (mapcar #'(lambda (x) (scalar-times-poly-list s x)) F))
|
---|
| 149 |
|
---|
| 150 | (defun monom-times-matrix (m F)
|
---|
| 151 | "Returns a product of a monomial M by a polynomial matrix F."
|
---|
| 152 | (mapcar #'(lambda (x) (monom-times-poly-list m x)) F))
|
---|
| 153 |
|
---|
| 154 | (defun term-times-matrix (term F)
|
---|
| 155 | "Returns a product of a term TERM by a polynomial matrix F."
|
---|
| 156 | (mapcar #'(lambda (x) (term-times-poly-list term x)) F))
|
---|
| 157 |
|
---|
| 158 |
|
---|
| 159 | ;; Polynomial list operations
|
---|
| 160 | (defun poly-list- (F G &optional (order #'lex>))
|
---|
| 161 | "Returns the list of differences of two lists of polynomials
|
---|
| 162 | F and G (polynomial maps)."
|
---|
| 163 | (mapcar #'(lambda (x y) (poly- x y order)) F G))
|
---|
| 164 |
|
---|
| 165 | (defun scalar-times-poly-list (s F)
|
---|
| 166 | "Returns the list of products of a polynomial S by the
|
---|
| 167 | list of polynomials F."
|
---|
| 168 | (mapcar #'(lambda (x) (scalar-times-poly s x)) F))
|
---|
| 169 |
|
---|
| 170 | (defun monom-times-poly-list (m f)
|
---|
| 171 | "Returns the list of products of a monomial M by the
|
---|
| 172 | list of polynomials F."
|
---|
| 173 | (mapcar #'(lambda (x) (monom-times-poly m x)) F))
|
---|
| 174 |
|
---|
| 175 | (defun term-times-poly-list (term f)
|
---|
| 176 | "Returns the list of products of a term TERM by the
|
---|
| 177 | list of polynomials F."
|
---|
| 178 | (mapcar #'(lambda (x) (term-times-poly term x)) F))
|
---|
| 179 |
|
---|
| 180 | ;; Generalized Characteristic polynomial operations
|
---|
| 181 | ;; det(A - u1*B1-u2*B2-...-um*Bm)
|
---|
| 182 |
|
---|
| 183 | (defun characteristic-combination (A B &optional (order #'lex>)
|
---|
| 184 | &aux (n (length B)))
|
---|
| 185 | "Returns A - U1 * B1 - U2 * B2 - ... - UM * BM where A is a polynomial
|
---|
| 186 | and B=(B1,B2,...,BM) is a polynomial list, where U1, U2, ... , UM are
|
---|
| 187 | new variables. These variables will be added to every polynomial
|
---|
| 188 | A and BI as the last M variables."
|
---|
| 189 | (setf A (poly-extend-end A (make-list n :initial-element 0)))
|
---|
| 190 | (dotimes (i (length B) A)
|
---|
| 191 | (setf A (poly- A (poly-extend-end (elt B i) (standard-vector n i)) order))))
|
---|
| 192 |
|
---|
| 193 | ;; A is a list of polynomials; B is a list of lists of polynomials
|
---|
| 194 | (defun characteristic-combination-poly-list (A B &optional (order #'lex>))
|
---|
| 195 | "Returns A - U1 * B1 - U2 * B2 - ... - UM * BM where A is a polynomial list
|
---|
| 196 | and B=(B1, B2, ... , BM) is a list of polynomial lists, where U1, U2, ... ,UM are
|
---|
| 197 | new variables. These variables will be added to every polynomial
|
---|
| 198 | A and BI as the last M variables. Se also CHARACTERISTIC-COMBINATION."
|
---|
| 199 | (apply #'mapcar
|
---|
| 200 | (cons #'(lambda (&rest x)
|
---|
| 201 | (funcall #'characteristic-combination (car x) (rest x) order))
|
---|
| 202 | (cons A B))))
|
---|
| 203 |
|
---|
| 204 | ;; Finally, the case of matrices
|
---|
| 205 | (defun characteristic-matrix (A &optional
|
---|
| 206 | (order #'lex>)
|
---|
| 207 | (B (list (identity-matrix (length A) (length (caaaar A))))))
|
---|
| 208 | "Returns A - U1*B1 - U2*B2 - ... - UM * BM where A is a polynomial matrix
|
---|
| 209 | and B=(B1,B2,...,BM) is a list of polynomial matrices, where U1, U2, .., UM are
|
---|
| 210 | new variables. These variables will be added to every polynomial
|
---|
| 211 | A and BI as the last M variables. Se also CHARACTERISTIC-COMBINATION."
|
---|
| 212 | (apply #'mapcar
|
---|
| 213 | (cons #'(lambda (&rest x)
|
---|
| 214 | (funcall #'characteristic-combination-poly-list (car x) (rest x) order))
|
---|
| 215 | (cons A B))))
|
---|
| 216 |
|
---|
| 217 | ;; Characteristic polynomial
|
---|
| 218 | (defun characteristic-polynomial (A &optional
|
---|
| 219 | (order #'lex>)
|
---|
| 220 | (B (list (identity-matrix (length A) (length (caaaar A))))))
|
---|
| 221 | "Returns the generalized characteristic polynomial, i.e. the determinant
|
---|
| 222 | DET(A - U1 * B1 - U2 * B2 - ... - UM * BM), where A and BI are square polynomial matrices in
|
---|
| 223 | N variables. The resulting polynomial will have N+M variables, with U1, U2, ..., UM added
|
---|
| 224 | as the last M variables."
|
---|
| 225 | (determinant (characteristic-matrix A order B) order))
|
---|
| 226 |
|
---|
| 227 | (defun identity-matrix (dim nvars)
|
---|
| 228 | "Return the polynomial matrix which is the identity matrix. DIM is the requested
|
---|
| 229 | dimension and NVARS is the number of variables of each entry."
|
---|
| 230 | (labels ((zero-monom () (make-list nvars :initial-element 0))
|
---|
| 231 | (entry (i j) (if (= i j)
|
---|
| 232 | (list (cons (zero-monom) 1))
|
---|
| 233 | nil)))
|
---|
| 234 | (makelist (makelist (entry i j) (i 1 dim)) (j 1 dim))))
|
---|
| 235 |
|
---|
| 236 | (defun print-matrix (F vars)
|
---|
| 237 | "Prints a polynomial matrix F, using a list of symbols VARS as variable names."
|
---|
| 238 | (mapcar #'(lambda (x) (poly-print (cons '[ x) vars) (terpri)) F)
|
---|
| 239 | F)
|
---|
| 240 |
|
---|
| 241 | (defun jacobi-matrix (F &optional (m (length F)) (n (length (caaaar F))))
|
---|
| 242 | "Returns the Jacobi matrix of a polynomial list F over the first N variables."
|
---|
| 243 | (makelist (makelist (scalar-partial (elt F i) j)
|
---|
| 244 | (j 0 (1- n)))
|
---|
| 245 | (i 0 (1- m))))
|
---|
| 246 |
|
---|
| 247 | (defun jacobian (F &optional (order #'lex>) (m (length F)) (n (length (caaaar F))))
|
---|
| 248 | "Returns the Jacobian (determinant) of a polynomial list F over the first N variables."
|
---|
| 249 | (determinant (jacobi-matrix F m n) order))
|
---|