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