1 | head 1.7;
|
---|
2 | access;
|
---|
3 | symbols;
|
---|
4 | locks; strict;
|
---|
5 | comment @;;; @;
|
---|
6 |
|
---|
7 |
|
---|
8 | 1.7
|
---|
9 | date 2009.01.23.10.49.32; author marek; state Exp;
|
---|
10 | branches;
|
---|
11 | next 1.6;
|
---|
12 |
|
---|
13 | 1.6
|
---|
14 | date 2009.01.23.10.45.47; author marek; state Exp;
|
---|
15 | branches;
|
---|
16 | next 1.5;
|
---|
17 |
|
---|
18 | 1.5
|
---|
19 | date 2009.01.23.10.43.33; author marek; state Exp;
|
---|
20 | branches;
|
---|
21 | next 1.4;
|
---|
22 |
|
---|
23 | 1.4
|
---|
24 | date 2009.01.22.04.01.19; author marek; state Exp;
|
---|
25 | branches;
|
---|
26 | next 1.3;
|
---|
27 |
|
---|
28 | 1.3
|
---|
29 | date 2009.01.19.09.25.14; author marek; state Exp;
|
---|
30 | branches;
|
---|
31 | next 1.2;
|
---|
32 |
|
---|
33 | 1.2
|
---|
34 | date 2009.01.19.07.49.33; author marek; state Exp;
|
---|
35 | branches;
|
---|
36 | next 1.1;
|
---|
37 |
|
---|
38 | 1.1
|
---|
39 | date 2009.01.19.06.48.12; author marek; state Exp;
|
---|
40 | branches;
|
---|
41 | next ;
|
---|
42 |
|
---|
43 |
|
---|
44 | desc
|
---|
45 | @@
|
---|
46 |
|
---|
47 |
|
---|
48 | 1.7
|
---|
49 | log
|
---|
50 | @*** empty log message ***
|
---|
51 | @
|
---|
52 | text
|
---|
53 | @;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
|
---|
54 | #|
|
---|
55 | $Id$
|
---|
56 | *--------------------------------------------------------------------------*
|
---|
57 | | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@@math.arizona.edu) |
|
---|
58 | | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
|
---|
59 | | |
|
---|
60 | | Everyone is permitted to copy, distribute and modify the code in this |
|
---|
61 | | directory, as long as this copyright note is preserved verbatim. |
|
---|
62 | *--------------------------------------------------------------------------*
|
---|
63 | |#
|
---|
64 |
|
---|
65 | (defpackage "DYNAMICS"
|
---|
66 | (:use "ORDER" "MONOM" "COEFFICIENT-RING" "GROBNER" "MAKELIST" "PRINTER" "TERM" "POLY" "COMMON-LISP")
|
---|
67 | (:export poly-composition poly-scalar-composition poly-dynamic-power
|
---|
68 | poly-evaluate poly-scalar-evaluate factorial
|
---|
69 | poly-scalar-diff poly-diff standard-vector
|
---|
70 | scalar-partial partial drop-elt drop-row drop-column
|
---|
71 | determinant minor matrix- poly-list-
|
---|
72 | characteristic-combination
|
---|
73 | characteristic-matrix
|
---|
74 | characteristic-polynomial
|
---|
75 | identity-matrix
|
---|
76 | print-matrix
|
---|
77 | jacobi-matrix
|
---|
78 | jacobian))
|
---|
79 |
|
---|
80 | (in-package "DYNAMICS")
|
---|
81 |
|
---|
82 | #+debug(proclaim '(optimize (speed 0) (debug 3)))
|
---|
83 | #-debug(proclaim '(optimize (speed 3) (debug 0)))
|
---|
84 |
|
---|
85 | (defun poly-scalar-composition (f G &optional (order #'lex>))
|
---|
86 | "Returns a polynomial obtained by substituting a list of polynomials G=(G1,G2,...,GN)
|
---|
87 | into a polynomial F(X1,X2,...,XN). All polynomials are assumed to be in the internal form,
|
---|
88 | so variables do not explicitly apprear in the calculation."
|
---|
89 | (reduce #'(lambda (x y) (poly+ x y order))
|
---|
90 | (mapcar #'(lambda (x)
|
---|
91 | (scalar-times-poly
|
---|
92 | (cdr x)
|
---|
93 | (poly-mexpt G (car x) order)))
|
---|
94 | f)))
|
---|
95 |
|
---|
96 | (defun poly-composition (F G &optional (order #'lex>))
|
---|
97 | "Return the composition of a polynomial map F with a polynomial map
|
---|
98 | G. Both maps are represented as lists of polynomials, and each polynomial
|
---|
99 | is in the internal alist representation. The restriction is that the length
|
---|
100 | of the list G must be the number of variables in the list F."
|
---|
101 | (mapcar #'(lambda (x) (poly-scalar-composition x G order)) F))
|
---|
102 |
|
---|
103 |
|
---|
104 | (defun poly-dynamic-power (F n &optional (order #'lex>))
|
---|
105 | "Calculate the composition FoFo...oF (n times), where
|
---|
106 | F is a polynomial map represented as a list of polynomials."
|
---|
107 | (reduce #'(lambda (x y) (poly-composition x y order))
|
---|
108 | (make-list n :initial-element F)))
|
---|
109 |
|
---|
110 | (defun poly-scalar-evaluate (f x &optional (order #'lex>))
|
---|
111 | "Evaluate a polynomial F at a point X. This operation is implemented
|
---|
112 | through POLY-SCALAR-COMPOSITION."
|
---|
113 | (if (endp f)
|
---|
114 | 0
|
---|
115 | (let ((p (poly-scalar-composition
|
---|
116 | f
|
---|
117 | (mapcar #'(lambda (u)
|
---|
118 | (if (zerop u)
|
---|
119 | nil
|
---|
120 | (list (cons (make-list (length (caar f)) :initial-element 0)
|
---|
121 | u))))
|
---|
122 | x)
|
---|
123 | order)))
|
---|
124 | (if (endp p) 0
|
---|
125 | (cdar p)))))
|
---|
126 |
|
---|
127 | (defun poly-evaluate (F x &optional (order #'lex>))
|
---|
128 | "Evaluate a polynomial map F, represented as list of polynomials, at a point X."
|
---|
129 | (mapcar #'(lambda (h) (poly-scalar-evaluate h x order)) F))
|
---|
130 |
|
---|
131 | (defun factorial (n &optional (k n) &aux (result 1))
|
---|
132 | "Return N!/(N-K)!=N(N-1)(N-K+1)."
|
---|
133 | (dotimes (j k result)
|
---|
134 | (setf result (* result (- n j)))))
|
---|
135 |
|
---|
136 | (defun poly-scalar-diff (f m)
|
---|
137 | "Return the partial derivative of a polynomial F over multiple
|
---|
138 | variables according to multiindex M."
|
---|
139 | (setf f (remove-if #'(lambda (x) (some #'< (car x) m)) f))
|
---|
140 | (mapcar #'(lambda (x) (cons (monom/ (car x) m)
|
---|
141 | (* (cdr x) (apply #'* (mapcar #'factorial (car x) m)))))
|
---|
142 | f))
|
---|
143 |
|
---|
144 | (defun poly-diff (F m)
|
---|
145 | "Return the partial derivative of a polynomial map F, represented as a list of polynomials,
|
---|
146 | with respect to several variables, according to multi-index M."
|
---|
147 | (mapcar #'(lambda (h) (poly-scalar-diff h m)) F))
|
---|
148 |
|
---|
149 |
|
---|
150 | (defun scalar-partial (f k &optional (l 1))
|
---|
151 | "Returns the L-th partial derivative of a polynomial F over the K-th variable."
|
---|
152 | (when f
|
---|
153 | (poly-scalar-diff f (standard-vector (length (caar f)) k l))))
|
---|
154 |
|
---|
155 | (defun partial (F k &optional (l 1))
|
---|
156 | "Returns the L-th partial derivative over the K-th variable, of a polynomial
|
---|
157 | map F, represented as a list of polynomials."
|
---|
158 | (when (and F (car F) (caar F))
|
---|
159 | (poly-diff F (standard-vector (length (caaar F)) k l))))
|
---|
160 |
|
---|
161 |
|
---|
162 | (defun determinant (F &optional (order #'lex>) &aux (result nil))
|
---|
163 | "Returns the determinant of a polynomial matrix F, which is a list of
|
---|
164 | rows of the matrix, each row is a list of polynomials. The algorithm
|
---|
165 | is recursive expansion along columns."
|
---|
166 | (cond
|
---|
167 | ((= (length F) 1) (setf result (caar F)))
|
---|
168 | (t
|
---|
169 | (dotimes (i (length F) result)
|
---|
170 | (setf result (poly+ result
|
---|
171 | (scalar-times-poly
|
---|
172 | (expt -1 i)
|
---|
173 | (poly* (car (elt F i)) (minor F i 0 order) order))
|
---|
174 | order))))))
|
---|
175 |
|
---|
176 |
|
---|
177 | (defun minor (F i j &optional (order #'lex>))
|
---|
178 | "Calculate the minor of a polynomial matrix F with respect to entry (I,J)."
|
---|
179 | (determinant (drop-row (drop-column F j) i) order))
|
---|
180 |
|
---|
181 | (defun drop-row (F i)
|
---|
182 | "Discards the I-th row from a polynomial matrix F."
|
---|
183 | (drop-elt F i))
|
---|
184 |
|
---|
185 | (defun drop-column (F j)
|
---|
186 | "Discards the J-th column from a polynomial matrix F."
|
---|
187 | (mapcar #'(lambda (x) (drop-elt x j)) F))
|
---|
188 |
|
---|
189 | (defun drop-elt (lst j)
|
---|
190 | "Discards the J-th element from a list LST."
|
---|
191 | (append (subseq lst 0 j) (subseq lst (1+ j))))
|
---|
192 |
|
---|
193 | ;; Matrix operations
|
---|
194 |
|
---|
195 | (defun matrix- (F G &optional (order #'lex>))
|
---|
196 | "Returns difference of two polynomial matrices F and G."
|
---|
197 | (mapcar #'(lambda (x y) (poly-list- x y order)) F G))
|
---|
198 |
|
---|
199 | (defun scalar-times-matrix (s F)
|
---|
200 | "Returns a product of a polynomial S by a polynomial matrix F."
|
---|
201 | (mapcar #'(lambda (x) (scalar-times-poly-list s x)) F))
|
---|
202 |
|
---|
203 | (defun monom-times-matrix (m F)
|
---|
204 | "Returns a product of a monomial M by a polynomial matrix F."
|
---|
205 | (mapcar #'(lambda (x) (monom-times-poly-list m x)) F))
|
---|
206 |
|
---|
207 | (defun term-times-matrix (term F)
|
---|
208 | "Returns a product of a term TERM by a polynomial matrix F."
|
---|
209 | (mapcar #'(lambda (x) (term-times-poly-list term x)) F))
|
---|
210 |
|
---|
211 |
|
---|
212 | ;; Polynomial list operations
|
---|
213 | (defun poly-list- (F G &optional (order #'lex>))
|
---|
214 | "Returns the list of differences of two lists of polynomials
|
---|
215 | F and G (polynomial maps)."
|
---|
216 | (mapcar #'(lambda (x y) (poly- x y order)) F G))
|
---|
217 |
|
---|
218 | (defun scalar-times-poly-list (s F)
|
---|
219 | "Returns the list of products of a polynomial S by the
|
---|
220 | list of polynomials F."
|
---|
221 | (mapcar #'(lambda (x) (scalar-times-poly s x)) F))
|
---|
222 |
|
---|
223 | (defun monom-times-poly-list (m f)
|
---|
224 | "Returns the list of products of a monomial M by the
|
---|
225 | list of polynomials F."
|
---|
226 | (mapcar #'(lambda (x) (monom-times-poly m x)) F))
|
---|
227 |
|
---|
228 | (defun term-times-poly-list (term f)
|
---|
229 | "Returns the list of products of a term TERM by the
|
---|
230 | list of polynomials F."
|
---|
231 | (mapcar #'(lambda (x) (term-times-poly term x)) F))
|
---|
232 |
|
---|
233 | ;; Generalized Characteristic polynomial operations
|
---|
234 | ;; det(A - u1*B1-u2*B2-...-um*Bm)
|
---|
235 |
|
---|
236 | (defun characteristic-combination (A B &optional (order #'lex>)
|
---|
237 | &aux (n (length B)))
|
---|
238 | "Returns A - U1 * B1 - U2 * B2 - ... - UM * BM where A is a polynomial
|
---|
239 | and B=(B1,B2,...,BM) is a polynomial list, where U1, U2, ... , UM are
|
---|
240 | new variables. These variables will be added to every polynomial
|
---|
241 | A and BI as the last M variables."
|
---|
242 | (setf A (poly-extend-end A (make-list n :initial-element 0)))
|
---|
243 | (dotimes (i (length B) A)
|
---|
244 | (setf A (poly- A (poly-extend-end (elt B i) (standard-vector n i)) order))))
|
---|
245 |
|
---|
246 | ;; A is a list of polynomials; B is a list of lists of polynomials
|
---|
247 | (defun characteristic-combination-poly-list (A B &optional (order #'lex>))
|
---|
248 | "Returns A - U1 * B1 - U2 * B2 - ... - UM * BM where A is a polynomial list
|
---|
249 | and B=(B1, B2, ... , BM) is a list of polynomial lists, where U1, U2, ... ,UM are
|
---|
250 | new variables. These variables will be added to every polynomial
|
---|
251 | A and BI as the last M variables. Se also CHARACTERISTIC-COMBINATION."
|
---|
252 | (apply #'mapcar
|
---|
253 | (cons #'(lambda (&rest x)
|
---|
254 | (funcall #'characteristic-combination (car x) (rest x) order))
|
---|
255 | (cons A B))))
|
---|
256 |
|
---|
257 | ;; Finally, the case of matrices
|
---|
258 | (defun characteristic-matrix (A &optional
|
---|
259 | (order #'lex>)
|
---|
260 | (B (list (identity-matrix (length A) (length (caaaar A))))))
|
---|
261 | "Returns A - U1*B1 - U2*B2 - ... - UM * BM where A is a polynomial matrix
|
---|
262 | and B=(B1,B2,...,BM) is a list of polynomial matrices, where U1, U2, .., UM are
|
---|
263 | new variables. These variables will be added to every polynomial
|
---|
264 | A and BI as the last M variables. Se also CHARACTERISTIC-COMBINATION."
|
---|
265 | (apply #'mapcar
|
---|
266 | (cons #'(lambda (&rest x)
|
---|
267 | (funcall #'characteristic-combination-poly-list (car x) (rest x) order))
|
---|
268 | (cons A B))))
|
---|
269 |
|
---|
270 | ;; Characteristic polynomial
|
---|
271 | (defun characteristic-polynomial (A &optional
|
---|
272 | (order #'lex>)
|
---|
273 | (B (list (identity-matrix (length A) (length (caaaar A))))))
|
---|
274 | "Returns the generalized characteristic polynomial, i.e. the determinant
|
---|
275 | DET(A - U1 * B1 - U2 * B2 - ... - UM * BM), where A and BI are square polynomial matrices in
|
---|
276 | N variables. The resulting polynomial will have N+M variables, with U1, U2, ..., UM added
|
---|
277 | as the last M variables."
|
---|
278 | (determinant (characteristic-matrix A order B) order))
|
---|
279 |
|
---|
280 | (defun identity-matrix (dim nvars)
|
---|
281 | "Return the polynomial matrix which is the identity matrix. DIM is the requested
|
---|
282 | dimension and NVARS is the number of variables of each entry."
|
---|
283 | (labels ((zero-monom () (make-list nvars :initial-element 0))
|
---|
284 | (entry (i j) (if (= i j)
|
---|
285 | (list (cons (zero-monom) 1))
|
---|
286 | nil)))
|
---|
287 | (makelist (makelist (entry i j) (i 1 dim)) (j 1 dim))))
|
---|
288 |
|
---|
289 | (defun print-matrix (F vars)
|
---|
290 | "Prints a polynomial matrix F, using a list of symbols VARS as variable names."
|
---|
291 | (mapcar #'(lambda (x) (poly-print (cons '[ x) vars) (terpri)) F)
|
---|
292 | F)
|
---|
293 |
|
---|
294 | (defun jacobi-matrix (F &optional (m (length F)) (n (length (caaaar F))))
|
---|
295 | "Returns the Jacobi matrix of a polynomial list F over the first N variables."
|
---|
296 | (makelist (makelist (scalar-partial (elt F i) j)
|
---|
297 | (j 0 (1- n)))
|
---|
298 | (i 0 (1- m))))
|
---|
299 |
|
---|
300 | (defun jacobian (F &optional (order #'lex>) (m (length F)) (n (length (caaaar F))))
|
---|
301 | "Returns the Jacobian (determinant) of a polynomial list F over the first N variables."
|
---|
302 | (determinant (jacobi-matrix F m n) order))
|
---|
303 | @
|
---|
304 |
|
---|
305 |
|
---|
306 | 1.6
|
---|
307 | log
|
---|
308 | @*** empty log message ***
|
---|
309 | @
|
---|
310 | text
|
---|
311 | @d3 1
|
---|
312 | a3 1
|
---|
313 | $Id: dynamics.lisp,v 1.5 2009/01/23 10:43:33 marek Exp marek $
|
---|
314 | a97 7
|
---|
315 | (defun standard-vector (n k &optional (coeff 1)
|
---|
316 | &aux (v (make-list n :initial-element 0)))
|
---|
317 | "Returns vector (0 0 ... 1 ... 0 0) of length N, where 1 appears on K-th place."
|
---|
318 | (setf (elt v k) coeff)
|
---|
319 | v)
|
---|
320 |
|
---|
321 |
|
---|
322 | @
|
---|
323 |
|
---|
324 |
|
---|
325 | 1.5
|
---|
326 | log
|
---|
327 | @*** empty log message ***
|
---|
328 | @
|
---|
329 | text
|
---|
330 | @d3 1
|
---|
331 | a3 1
|
---|
332 | $Id: dynamics.lisp,v 1.4 2009/01/22 04:01:19 marek Exp marek $
|
---|
333 | d17 1
|
---|
334 | a17 1
|
---|
335 | poly-scalar-diff poly-diff std-vector
|
---|
336 | d98 1
|
---|
337 | a98 1
|
---|
338 | (defun std-vector (n k &optional (coeff 1)
|
---|
339 | d108 1
|
---|
340 | a108 1
|
---|
341 | (poly-scalar-diff f (std-vector (length (caar f)) k l))))
|
---|
342 | d114 1
|
---|
343 | a114 1
|
---|
344 | (poly-diff F (std-vector (length (caaar F)) k l))))
|
---|
345 | d199 1
|
---|
346 | a199 1
|
---|
347 | (setf A (poly- A (poly-extend-end (elt B i) (std-vector n i)) order))))
|
---|
348 | @
|
---|
349 |
|
---|
350 |
|
---|
351 | 1.4
|
---|
352 | log
|
---|
353 | @*** empty log message ***
|
---|
354 | @
|
---|
355 | text
|
---|
356 | @d3 1
|
---|
357 | a3 1
|
---|
358 | $Id$
|
---|
359 | d17 1
|
---|
360 | a17 1
|
---|
361 | poly-scalar-diff poly-diff standard-vector
|
---|
362 | d98 1
|
---|
363 | a98 1
|
---|
364 | (defun standard-vector (n k &optional (coeff 1)
|
---|
365 | d108 1
|
---|
366 | a108 1
|
---|
367 | (poly-scalar-diff f (standard-vector (length (caar f)) k l))))
|
---|
368 | d114 1
|
---|
369 | a114 1
|
---|
370 | (poly-diff F (standard-vector (length (caaar F)) k l))))
|
---|
371 | d199 1
|
---|
372 | a199 1
|
---|
373 | (setf A (poly- A (poly-extend-end (elt B i) (standard-vector n i)) order))))
|
---|
374 | @
|
---|
375 |
|
---|
376 |
|
---|
377 | 1.3
|
---|
378 | log
|
---|
379 | @*** empty log message ***
|
---|
380 | @
|
---|
381 | text
|
---|
382 | @d30 2
|
---|
383 | a31 2
|
---|
384 | ;;(proclaim '(optimize (speed 0) (debug 3)))
|
---|
385 | (proclaim '(optimize (speed 3) (debug 0)))
|
---|
386 | @
|
---|
387 |
|
---|
388 |
|
---|
389 | 1.2
|
---|
390 | log
|
---|
391 | @*** empty log message ***
|
---|
392 | @
|
---|
393 | text
|
---|
394 | @d30 2
|
---|
395 | a31 1
|
---|
396 | (proclaim '(optimize (speed 0) (debug 3)))
|
---|
397 | @
|
---|
398 |
|
---|
399 |
|
---|
400 | 1.1
|
---|
401 | log
|
---|
402 | @Initial revision
|
---|
403 | @
|
---|
404 | text
|
---|
405 | @d3 1
|
---|
406 | a3 1
|
---|
407 | $Id: dynamics.lisp,v 1.18 1997/12/13 07:10:00 marek Exp $
|
---|
408 | d30 1
|
---|
409 | @
|
---|