source: CGBLisp/src/poly.lisp@ 1

#|
        $Id: poly.lisp,v 1.6 2009/01/22 04:05:52 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 "POLY"
  (:export monom-times-poly scalar-times-poly term-times-poly
           minus-poly poly+ poly- poly* poly-expt
           sort-poly poly-op poly-constant-p poly-extend
           poly-mexpt poly-extend-end
           lt lm lc poly-zerop)
  (:use "ORDER" "MONOM" "TERM" "COEFFICIENT-RING" "COMMON-LISP"))
(in-package "POLY")

#+debug(proclaim '(optimize (speed 0) (debug 3)))
#-debug(proclaim '(optimize (speed 3) (debug 0)))

;;----------------------------------------------------------------
;; BASIC OPERATIONS ON POLYNOMIALS of many variables
;; Hybrid operations involving polynomials and monomials or terms
;;----------------------------------------------------------------
;; The representations are as follows:
;;
;;      monom:          (n1 n2 ... nk) where ni are non-negative integers
;;      term:           (monom . coefficient)
;;      polynomial:     (term1 term2 ... termk)
;;
;; The terms in a polynomial are assumed to be sorted in descending order by
;; some admissible well-ordering, typically defaulting to the lexicographic
;; order from  the "order" package.  Thus, the polynomials are represented
;; non-recursively as alists.
;;----------------------------------------------------------------
;; EXAMPLES: Suppose that variables are x and y. Then
;; Monom x*y^2 ---> (1 2)
;; Term  3*x*y^2 ---> ((1 2) . 3) 
;; Polynomial x^2+x*y ---> (((2 0) . 1) ((1 1) . 1))
;;----------------------------------------------------------------


(defun scalar-times-poly (c p &optional (ring *coefficient-ring*))
  "Return product of a scalar C by a polynomial P with coefficient ring RING."
  (unless (funcall (ring-zerop ring) c)
    (mapcar
     #'(lambda (term) (cons (car term) (funcall (ring-* ring) c (cdr term))))
     p)))

(defun term-times-poly (term f &optional (ring *coefficient-ring*))
  "Return product of a term TERM by a polynomial F with coefficient ring RING."
  (mapcar #'(lambda (x) (term* term x ring)) f))

(defun monom-times-poly (m f)
  "Return product of a monomial M by a polynomial F with coefficient ring RING."
  (mapcar #'(lambda (x) (monom-times-term m x)) f))

(defun minus-poly (f &optional (ring *coefficient-ring*))
  "Changes the sign of a polynomial F with coefficients in coefficient ring
RING, and returns the result."
  (mapcar #'(lambda (x) (cons (car x) (funcall (ring-- ring) (cdr x)))) f))

(defun sort-poly (poly
                  &optional 
                  (pred #'lex>)
                  (start 0)
                  (end (unless (null poly) ;zero
                         (length (caar poly)))))
  "Destructively Sorts a polynomial POLY by predicate PRED; the predicate is assumed
to take arguments START and END in addition to the pair of
monomials, as the functions in the ORDER package do."
  (sort poly #'(lambda (p q) (funcall pred p q :start start :end end))
        :key #'first))

(defun poly+ (p q &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns the sum of two polynomials P and Q with coefficients in
ring RING, with terms ordered according to monomial order PRED."
  (do (r)
      (nil)
    (cond 
      ((endp p) (return (nreconc r q)))
      ((endp q) (return (nreconc r p)))
      (t 
       (multiple-value-bind 
             (mgreater mequal)
           (funcall pred (caar p) (caar q))
         (cond
           (mequal
            (let ((s (funcall (ring-+ ring) (cdar p) (cdar q))))
              (unless (funcall (ring-zerop ring) s)             ;check for cancellation
                      (setf r (cons (cons (caar p) s) r)))
              (setf p (cdr p) q (cdr q))))
           (mgreater
            (setf r (cons (car p) r)
                  p (cdr p)))
           (t (setf r (cons (car q) r)
                    q (cdr q)))))))))

(defun poly- (p q &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns the difference of two polynomials P and Q with coefficients
in ring RING, with terms ordered according to monomial order PRED."
  (do (r done)
      (done r)
    (cond 
      ((endp p) (return (revappend r (minus-poly q ring))))
      ((endp q) (return (revappend r p)))
      (t 
       (multiple-value-bind 
             (mgreater mequal)
           (funcall pred (caar p) (caar q))
         (cond
           (mequal
            (let ((s (funcall (ring-- ring) (cdar p) (cdar q))))
              (unless (zerop s)         ;check for cancellation
                (setf r (cons (cons (caar p) s) r)))
              (setf p (cdr p) q (cdr q))))
           (mgreater
            (setf r (cons (car p) r)
                  p (cdr p)))
           (t (setf r (cons (cons (caar q) (funcall (ring-- ring) (cdar q))) r)
                    q (cdr q)))))))))

;; Multiplication of polynomials
;; Non-destructive version
(defun poly* (p q &optional (pred #'lex>) (ring *coefficient-ring*))
  "Returns the product of two polynomials P and Q with coefficients in
ring RING, with terms ordered according to monomial order PRED."
  (cond ((or (endp p) (endp q)) nil)    ;p or q is 0 (represented by NIL)
        ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
        (t (cons (cons  (monom* (caar p) (caar q))
                        (funcall (ring-* ring) (cdar p) (cdar q)))
                   (poly+ (term-times-poly (car p) (cdr q) ring)
                          (poly* (cdr p) q pred ring)
                          pred ring)))))

(defun poly-op (f m g pred ring)
  "Returns F-M*G, where F and G are polynomials with coefficients in
ring RING, ordered according to monomial order PRED and M is a
monomial."
  (poly- f (term-times-poly m g ring) pred ring))

(defun poly-expt (poly n &optional (pred #'lex>) (ring *coefficient-ring*))
  "Exponentiate a polynomial POLY to power N. The terms of the
polynomial are assumed to be ordered by monomial order PRED and with
coefficients in ring RING.  Use the Chinese algorithm; assume N>=0 and
POLY is non-zero (not NIL)."
  (labels ((poly-one ()
             (list
              (cons (make-list (length (caar poly)) :initial-element 0) (ring-unit ring)))))
    (cond
     ((minusp n) (error "Negative exponent."))
     ((endp poly) (if (zerop n) (poly-one) nil))
     (t
      (do ((k 1 (ash k 1))
           (q poly (poly* q q pred ring))       ;keep squaring
           (p (poly-one) (if (not (zerop (logand k n))) (poly* p q pred ring) p)))
          ((> k n) p))))))

(defun poly-mexpt (plist monom &optional (pred #'lex>) (ring *coefficient-ring*))
  "Raise a polynomial vector represented ad a list of polynomials
PLIST to power MULTIINDEX.  Every polynomial has its terms ordered by
predicate PRED and coefficients in the ring RING."
  (reduce 
   #'(lambda (u v) (poly* u v pred ring))
   (mapcan #'(lambda (y i)
               (cond ((endp y) (if (zerop i) nil (list nil)))
                     (t (list (poly-expt y i pred ring)))))
           plist monom)))

(defun poly-constant-p (p)
  "Returns T if P is a constant polynomial."
  (and (= (length p) 1)
       (every #'zerop (caar p))))

       
(defun  poly-extend (p &optional (m (list 0)))
  "Given a polynomial P in k[x[r+1],...,xn], it returns the same
polynomial as an element of k[x1,...,xn], optionally multiplying it by
a monomial x1^m1*x2^m2*...*xr^mr, where m=(m1,m2,...,mr) is a
multiindex."
  (mapcar #'(lambda (term) (cons (append m (car term)) (cdr term))) p))

(defun  poly-extend-end (p &optional (m (list 0)))
  "Similar to POLY-EXTEND, but it adds new variables at the end."
  (mapcar #'(lambda (term) (cons (append (car term) m) (cdr term))) p))

(defun poly-zerop (p)
  "Returns T if P is a zero polynomial."
  (null p))

(defun lt (p)
  "Returns the leading term of a polynomial P."
  #+debugging(assert (consp p))
  (first p))

(defun lm (p)
  "Returns the leading monomial of a polynomial P."
  #+debugging(assert (consp (lt p)))
  (car (lt p)))

(defun lc (p)
  "Returns the leading coefficient of a polynomial P."
  #+debugging(assert (consp (lt p)))
  (cdr (lt p)))