source: branches/f4grobner/polynomial.lisp@ 52

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; Polynomials
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defstruct (poly
             ;;
             ;; BOA constructor, by default constructs zero polynomial
             (:constructor make-poly-from-termlist (termlist &optional (sugar (termlist-sugar termlist))))
             (:constructor make-poly-zero (&aux (termlist nil) (sugar -1)))
             ;; Constructor of polynomials representing a variable
             (:constructor make-variable (ring nvars pos &optional (power 1)
                                              &aux
                                              (termlist (list
                                                         (make-term-variable ring nvars pos power)))
                                              (sugar power)))
            (:constructor poly-unit (ring dimension
                                     &aux
                                     (termlist (termlist-unit ring dimension))
                                     (sugar 0))))
  (termlist nil :type list)
  (sugar -1 :type fixnum))

;; Leading term
(defmacro poly-lt (p) `(car (poly-termlist ,p)))

;; Second term
(defmacro poly-second-lt (p) `(cadar (poly-termlist ,p)))

;; Leading monomial
(defun poly-lm (p) (term-monom (poly-lt p)))

;; Second monomial
(defun poly-second-lm (p) (term-monom (poly-second-lt p)))

;; Leading coefficient
(defun poly-lc (p) (term-coeff (poly-lt p)))

;; Second coefficient
(defun poly-second-lc (p) (term-coeff (poly-second-lt p)))

;; Testing for a zero polynomial
(defun poly-zerop (p) (null (poly-termlist p)))

;; The number of terms
(defun poly-length (p) (length (poly-termlist p)))

(defun scalar-times-poly (ring c p)
  (make-poly-from-termlist (scalar-times-termlist ring c (poly-termlist p)) (poly-sugar p)))

;; The scalar product omitting the head term
(defun scalar-times-poly-1 (ring c p)
  (make-poly-from-termlist (scalar-times-termlist ring c (cdr (poly-termlist p))) (poly-sugar p)))
    
(defun monom-times-poly (m p)
  (make-poly-from-termlist (monom-times-termlist m (poly-termlist p)) (+ (poly-sugar p) (monom-sugar m))))

(defun term-times-poly (ring term p)
  (make-poly-from-termlist (term-times-termlist ring term (poly-termlist p)) (+ (poly-sugar p) (term-sugar term))))

(defun poly-add (ring p q)
  (make-poly-from-termlist (termlist-add ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))

(defun poly-sub (ring p q)
  (make-poly-from-termlist (termlist-sub ring (poly-termlist p) (poly-termlist q)) (max (poly-sugar p) (poly-sugar q))))

(defun poly-uminus (ring p)
  (make-poly-from-termlist (termlist-uminus ring (poly-termlist p)) (poly-sugar p)))

(defun poly-mul (ring p q)
  (make-poly-from-termlist (termlist-mul ring (poly-termlist p) (poly-termlist q)) (+ (poly-sugar p) (poly-sugar q))))

(defun poly-expt (ring p n)
  (make-poly-from-termlist (termlist-expt ring (poly-termlist p) n) (* n (poly-sugar p))))

(defun poly-append (&rest plist)
  (make-poly-from-termlist (apply #'append (mapcar #'poly-termlist plist))
             (apply #'max (mapcar #'poly-sugar plist))))

(defun poly-nreverse (p)
  (setf (poly-termlist p) (nreverse (poly-termlist p)))
  p)

(defun poly-contract (p &optional (k 1))
  (make-poly-from-termlist (termlist-contract (poly-termlist p) k)
             (poly-sugar p)))

(defun poly-extend (p &optional (m (make-monom 1 :initial-element 0)))
  (make-poly-from-termlist
   (termlist-extend (poly-termlist p) m)
   (+ (poly-sugar p) (monom-sugar m))))

(defun poly-add-variables (p k)
  (setf (poly-termlist p) (termlist-add-variables (poly-termlist p) k))
  p)

(defun poly-list-add-variables (plist k)
  (mapcar #'(lambda (p) (poly-add-variables p k)) plist))

(defun poly-standard-extension (plist &aux (k (length plist)))
  "Calculate [U1*P1,U2*P2,...,UK*PK], where PLIST=[P1,P2,...,PK]."
  (declare (list plist) (fixnum k))
  (labels ((incf-power (g i)
             (dolist (x (poly-termlist g))
               (incf (monom-elt (term-monom x) i)))
             (incf (poly-sugar g))))
    (setf plist (poly-list-add-variables plist k))
    (dotimes (i k plist)
      (incf-power (nth i plist) i))))

(defun saturation-extension (ring f plist &aux (k (length plist)) (d (monom-dimension (poly-lm (car plist)))))
  "Calculate [F, U1*P1-1,U2*P2-1,...,UK*PK-1], where PLIST=[P1,P2,...,PK]."
  (setf f (poly-list-add-variables f k)
        plist (mapcar #'(lambda (x)
                          (setf (poly-termlist x) (nconc (poly-termlist x)
                                                         (list (make-term (make-monom d :initial-element 0)
                                                                          (funcall (ring-uminus ring) (funcall (ring-unit ring)))))))
                          x)
                      (poly-standard-extension plist)))
  (append f plist))


(defun polysaturation-extension (ring f plist &aux (k (length plist))
                                                   (d (+ k (length (poly-lm (car plist))))))
  "Calculate [F, U1*P1+U2*P2+...+UK*PK-1], where PLIST=[P1,P2,...,PK]."
  (setf f (poly-list-add-variables f k)
        plist (apply #'poly-append (poly-standard-extension plist))
        (cdr (last (poly-termlist plist))) (list (make-term (make-monom d :initial-element 0)
                                                            (funcall (ring-uminus ring) (funcall (ring-unit ring))))))
  (append f (list plist)))

(defun saturation-extension-1 (ring f p) (polysaturation-extension ring f (list p)))