source: CGBLisp/trunk/src/dynamics.lisp@ 14

;;; -*- Mode: Lisp; Syntax: Common-Lisp; Package: Grobner; Base: 10 -*-
#|
        $Id: dynamics.lisp,v 1.7 2009/01/23 10:49:32 marek Exp $
  *--------------------------------------------------------------------------*
  |  Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu)    |
  |    Department of Mathematics, University of Arizona, Tucson, AZ 85721    |
  |                                                                          |
  | Everyone is permitted to copy, distribute and modify the code in this    |
  | directory, as long as this copyright note is preserved verbatim.         |
  *--------------------------------------------------------------------------*
|#

(defpackage "DYNAMICS"
  (:use "ORDER" "MONOM" "COEFFICIENT-RING" "GROBNER" "MAKELIST" "PRINTER" "TERM" "POLY" "COMMON-LISP")
  (:export poly-composition poly-scalar-composition poly-dynamic-power
           poly-evaluate poly-scalar-evaluate factorial
           poly-scalar-diff poly-diff standard-vector
           scalar-partial partial drop-elt drop-row drop-column
           determinant minor matrix- poly-list-
           characteristic-combination
           characteristic-matrix
           characteristic-polynomial
           identity-matrix
           print-matrix
           jacobi-matrix
           jacobian))

(in-package "DYNAMICS")

(proclaim '(optimize (speed 0) (debug 3)))

(defun poly-scalar-composition (f G &optional (order #'lex>))
  "Returns a polynomial obtained by substituting a list of polynomials G=(G1,G2,...,GN)
 into a polynomial F(X1,X2,...,XN). All polynomials are assumed to be in the internal form,
so variables do not explicitly apprear in the calculation."
  (reduce #'(lambda (x y) (poly+ x y order))
          (mapcar #'(lambda (x)
                      (scalar-times-poly
                       (cdr x)
                       (poly-mexpt G (car x) order)))
                  f)))

(defun poly-composition (F G &optional (order #'lex>))
  "Return the composition of a polynomial map F with a polynomial map
G. Both maps are represented as lists of polynomials, and each polynomial
is in the internal alist representation. The restriction is that the length
of the list G must be the number of variables in the list F."
  (mapcar #'(lambda (x) (poly-scalar-composition x G order)) F))


(defun poly-dynamic-power (F n &optional (order #'lex>))
  "Calculate the composition FoFo...oF (n times), where
F is a polynomial map represented as a list of polynomials."
  (reduce #'(lambda (x y) (poly-composition x y order))
          (make-list n :initial-element F)))

(defun poly-scalar-evaluate (f x &optional (order #'lex>))
  "Evaluate a polynomial F at a point X. This operation is implemented
through POLY-SCALAR-COMPOSITION."
  (if (endp f)
      0
    (let ((p (poly-scalar-composition
              f
              (mapcar #'(lambda (u)
                          (if (zerop u)
                              nil
                            (list (cons (make-list (length (caar f)) :initial-element 0)
                                        u))))
                      x)
              order)))
      (if (endp p) 0
        (cdar p)))))
                        
(defun poly-evaluate (F x &optional (order #'lex>))
  "Evaluate a polynomial map F, represented as list of polynomials, at a point X."
  (mapcar #'(lambda (h) (poly-scalar-evaluate h x order)) F))

(defun factorial (n &optional (k n) &aux (result 1))
  "Return N!/(N-K)!=N(N-1)(N-K+1)."
  (dotimes (j k result)
    (setf result (* result (- n j)))))

(defun poly-scalar-diff (f m)
  "Return the partial derivative of a polynomial F over multiple
  variables according to multiindex M."
  (setf f (remove-if #'(lambda (x) (some #'< (car x) m)) f))
  (mapcar #'(lambda (x) (cons (monom/ (car x) m)
                              (* (cdr x) (apply #'* (mapcar #'factorial (car x) m)))))
          f))
  
(defun poly-diff (F m)
  "Return the partial derivative of a polynomial map F, represented as a list of polynomials,
with respect to several variables, according to multi-index M."
  (mapcar #'(lambda (h) (poly-scalar-diff h m)) F))


(defun scalar-partial (f k &optional (l 1))
  "Returns the L-th partial derivative of a polynomial F over the K-th variable."
  (when f
    (poly-scalar-diff f (standard-vector (length (caar f)) k l))))

(defun partial (F k &optional (l 1))
  "Returns the L-th partial derivative over the K-th variable, of a polynomial
map F, represented as a list of polynomials."
  (when (and F (car F) (caar F))
    (poly-diff F (standard-vector (length (caaar F)) k l))))


(defun determinant (F &optional (order #'lex>) &aux (result nil))
  "Returns the determinant of a polynomial matrix F, which is a list of
rows of the matrix, each row is a list of polynomials.  The algorithm
is recursive expansion along columns."
  (cond 
   ((= (length F) 1) (setf result (caar F)))
   (t
    (dotimes (i (length F) result)
      (setf result (poly+ result
                          (scalar-times-poly
                           (expt -1 i)
                           (poly* (car (elt F i)) (minor F i 0 order) order))
                          order))))))


(defun minor (F i j &optional (order #'lex>))
  "Calculate the minor of a polynomial matrix F with respect to entry (I,J)."
  (determinant (drop-row (drop-column F j) i) order))

(defun drop-row (F i)
  "Discards the I-th row from a polynomial matrix F."
  (drop-elt F i))

(defun drop-column (F j)
  "Discards the J-th column from a polynomial matrix F."
  (mapcar #'(lambda (x) (drop-elt x j)) F))

(defun drop-elt (lst j)
  "Discards the J-th element from a list LST."
  (append (subseq lst 0 j) (subseq lst (1+ j))))

;; Matrix operations

(defun matrix- (F G &optional (order #'lex>))
  "Returns difference of two polynomial matrices F and G."
  (mapcar #'(lambda (x y) (poly-list- x y order)) F G))

(defun scalar-times-matrix (s F)
  "Returns a product of a polynomial S by a polynomial matrix F."
  (mapcar #'(lambda (x) (scalar-times-poly-list s x)) F))

(defun monom-times-matrix (m F)
  "Returns a product of a monomial M by a polynomial matrix F."  
  (mapcar #'(lambda (x) (monom-times-poly-list m x)) F))

(defun term-times-matrix (term F)
  "Returns a product of a term TERM by a polynomial matrix F."  
  (mapcar #'(lambda (x) (term-times-poly-list term x)) F))


;; Polynomial list operations
(defun poly-list- (F G &optional (order #'lex>))
  "Returns the list of differences of two lists of polynomials
F and G (polynomial maps)."
  (mapcar #'(lambda (x y) (poly- x y order)) F G))

(defun scalar-times-poly-list (s F)
  "Returns the list of products of a polynomial S by the
list of polynomials F."
  (mapcar #'(lambda (x) (scalar-times-poly s x)) F))

(defun monom-times-poly-list (m f)
  "Returns the list of products of a monomial M by the
list of polynomials F."
  (mapcar #'(lambda (x) (monom-times-poly m x)) F))
       
(defun term-times-poly-list (term f)
  "Returns the list of products of a term TERM by the
list of polynomials F."
  (mapcar #'(lambda (x) (term-times-poly term x)) F))
       
;; Generalized Characteristic polynomial operations
;; det(A - u1*B1-u2*B2-...-um*Bm)

(defun characteristic-combination (A B &optional (order #'lex>)
                                                 &aux (n (length B)))
  "Returns A - U1 * B1 - U2 * B2 - ... - UM * BM where A is a polynomial
and B=(B1,B2,...,BM) is a polynomial list, where U1, U2, ... , UM are
new variables. These variables will be added to every polynomial
A and BI as the last M variables."
  (setf A (poly-extend-end A (make-list n :initial-element 0)))
  (dotimes (i (length B) A)
    (setf A (poly- A  (poly-extend-end (elt B i) (standard-vector n i)) order))))
                                                 
;; A is a list of polynomials; B is a list of lists of polynomials
(defun characteristic-combination-poly-list (A B &optional (order #'lex>))
  "Returns A - U1 * B1 - U2 * B2 - ... - UM * BM where A is a polynomial list
and B=(B1, B2, ... , BM) is a list of polynomial lists, where U1, U2, ... ,UM are
new variables. These variables will be added to every polynomial
A and BI as the last M variables. Se also CHARACTERISTIC-COMBINATION."
  (apply #'mapcar 
         (cons #'(lambda (&rest x)
                   (funcall #'characteristic-combination (car x) (rest x) order))
               (cons A B))))
  
;; Finally, the case of matrices
(defun characteristic-matrix (A &optional
                                (order #'lex>)
                                (B (list (identity-matrix (length A) (length (caaaar A))))))
  "Returns A - U1*B1 - U2*B2 - ... - UM * BM where A is a polynomial matrix
and B=(B1,B2,...,BM) is a list of polynomial matrices, where U1, U2, .., UM are
new variables. These variables will be added to every polynomial
A and BI as the last M variables. Se also CHARACTERISTIC-COMBINATION."
  (apply #'mapcar 
         (cons #'(lambda (&rest x)
                   (funcall #'characteristic-combination-poly-list (car x) (rest x) order))
               (cons A B))))

;; Characteristic polynomial
(defun characteristic-polynomial (A &optional
                                    (order #'lex>)
                                    (B (list (identity-matrix (length A) (length (caaaar A))))))
  "Returns the generalized characteristic polynomial, i.e. the determinant
DET(A - U1 * B1 - U2 * B2 - ... - UM * BM), where A and BI are square polynomial matrices in
N variables. The resulting polynomial will have N+M variables, with U1, U2, ..., UM added
as the last M variables."
  (determinant (characteristic-matrix A order B) order))

(defun identity-matrix (dim nvars)
  "Return the polynomial matrix which is the identity matrix. DIM is the requested
dimension and NVARS is the number of variables of each entry."
  (labels ((zero-monom () (make-list nvars :initial-element 0))
           (entry (i j) (if (= i j)
                            (list (cons (zero-monom) 1))
                          nil)))
    (makelist (makelist (entry i j) (i 1 dim)) (j 1 dim))))

(defun print-matrix (F vars)
  "Prints a polynomial matrix F, using a list of symbols VARS as variable names."
  (mapcar #'(lambda (x) (poly-print (cons '[ x) vars) (terpri)) F)
  F)

(defun jacobi-matrix (F &optional (m (length F)) (n (length (caaaar F))))
  "Returns the Jacobi matrix of a polynomial list F over the first N variables."
  (makelist (makelist (scalar-partial (elt F i) j)
                      (j 0 (1- n)))
            (i 0 (1- m))))

(defun jacobian (F &optional (order #'lex>) (m (length F)) (n (length (caaaar F))))
  "Returns the Jacobian (determinant) of a polynomial list F over the first N variables."
  (determinant (jacobi-matrix F m n) order))