source: CGBLisp/src/svpoly.lisp@ 8

#|
        $Id: svpoly.lisp,v 1.4 2009/01/23 10:37:28 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 "SVPOLY"
  (:export))

(in-package "SVPOLY")

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

(defstruct (svpoly (:constructor make-svpoly-raw))
  (nvars "Number of variables." (:type fixnum))
  (coefficient-type "Type of the coefficient.")
  (order "Monomial order.")
  (terms "The array of terms."))

(defun make-svpoly (alist order
                    &aux (nterms (length alist))
                         (nvars (length (caar alist)))
                         (coefficient-type (type-of (cdar alist))))
  "Construct svpolynomial from ALIST."
  (let ((svp (make-svpoly-raw
              :nvars nvars
              :coefficient-type (type-of (cdar alist))
              :order order
              :terms (make-array nterms :element-type 'cons))))
    (dotimes (i nterms)
      (let ((term (nth i alist)))
        (setf (svref (svpoly-terms svp) i)
          (cons (make-array nvars :initial-contents (car term))
                (cdr term)))))
    (svpoly-sort svp)))

(defun svpoly-sort (svp)
  "Destructively sorts an sv-polynomial SVP."
  (setf (svpoly-terms svp) (sort (svpoly-terms svp) (svpoly-order svp) :key #'car)) svp)

(defun make-order (nvars order)
  "Returns an order function with two parameters P and Q such that if called
on a pair of monomials with exactly NVARS variables, this function
will return T if P is greater than Q and NIL otherwise. The keyword
ORDER indicates one of several standard orders (:LEX, etc)."
  (declare (fixnum nvars))
  (ecase order
    (:lex #'(lambda (p q)
              (dotimes (i nvars (values NIL T))
                (declare (fixnum i))
                (cond
                 ((> (svref p i) (svref q i)) 
                  (return (values t nil)))
                 ((< (svref p i) (svref q i))
                  (return (values nil nil)))))))))

(defun svpoly-add (p q)
  "Adds polynomials P and Q, where P and Q are assumed to be ordered
by the same monomial order. Destructive to P and Q."
  (setf (svpoly-terms p) (add-terms (svpoly-terms p) (svpoly-terms q)
                                    (svpoly-order p)))
  p)

(defun add-terms (p q pred &aux (lp (length p)) (lq (length q)))
  (do ((r (make-array (+ lp lq)  :element-type 'cons))
       (i 0) (j 0) (k 0) done)
      (done (coerce (adjust-array r k) 'simple-vector))
    (cond 
     ((= i lp)
      (do nil
          ((>= j lq) (setf done t))
        (setf (aref r k) (svref q j))
        (incf j) (incf k)))
     ((= j lq)
      (do nil
          ((>= i lp) (setf done t))
        (setf (aref r k) (svref p i))
        (incf i) (incf k))) 
     (t
      (multiple-value-bind 
          (mgreater mequal)
          (funcall pred (car (svref p i)) (car (svref q j)))
        (cond
         (mequal
          (let ((s (+ (cdr (svref p i)) (cdr (svref q j)))))
            (unless (zerop s)           ;check for cancellation
              (setf (aref r k) (cons (car (svref p i)) s))
              (incf k))
            (incf i)
            (incf j)))
         (mgreater
          (setf (aref r k) (svref p i))
          (incf i)
          (incf k))
         (t
          (setf (aref r k) (svref q j))
          (incf j) 
          (incf k))))))))

#|

(defun scalar-times-poly (c p)
  "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 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 done)
      (done r)
    (cond 
      ((endp p) (setf r (append (nreverse r) q) done t))
      ((endp q) (setf r (append (nreverse r) p) done t))
      (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) (setf r (append (nreverse r) (minus-poly q ring)) done t))
      ((endp q) (setf r (append (nreverse r) p) done t))
      (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)))

|#