source: CGBLisp/src/RCS/poly.lisp,v@ 1

head    1.6;
access;
symbols;
locks; strict;
comment @;;; @;


1.6
date    2009.01.22.04.05.52;    author marek;   state Exp;
branches;
next    1.5;

1.5
date    2009.01.19.09.28.37;    author marek;   state Exp;
branches;
next    1.4;

1.4
date    2009.01.19.07.56.24;    author marek;   state Exp;
branches;
next    1.3;

1.3
date    2009.01.19.07.51.39;    author marek;   state Exp;
branches;
next    1.2;

1.2
date    2009.01.19.07.37.58;    author marek;   state Exp;
branches;
next    1.1;

1.1
date    2009.01.19.06.48.25;    author marek;   state Exp;
branches;
next    ;


desc
@@


1.6
log
@*** empty log message ***
@
text
@#|
        $Id$    
  *--------------------------------------------------------------------------*
  |  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)))

@


1.5
log
@*** empty log message ***
@
text
@d20 2
a21 2
;;(proclaim '(optimize (speed 0) (debug 3)))
(proclaim '(optimize (speed 3) (debug 0)))
@


1.4
log
@*** empty log message ***
@
text
@d20 2
a21 1
(proclaim '(optimize (speed 0) (debug 3)))
@


1.3
log
@*** empty log message ***
@
text
@d80 1
a80 1
      ()
@


1.2
log
@*** empty log message ***
@
text
@a187 3
(eval-when (compile)
  (proclaim '(inline lt lm lc poly-zerop scalar-times-poly)))

@


1.1
log
@Initial revision
@
text
@d2 1
a2 1
        $Id: poly.lisp,v 1.33 1998/01/18 09:05:28 marek Exp $   
d20 1
a20 4
(eval-when (compile)
  (proclaim '(optimize (speed 3) (safety 0)))
  (proclaim '(inline monom-times-poly scalar-times-poly 
              term-times-poly minus-poly sort-poly poly-op)))
@