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


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

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

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

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


desc
@@


1.4
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 "SVPOLY"
  (:export))

(in-package "SVPOLY")

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

(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)))

|#@


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

@


1.2
log
@*** empty log message ***
@
text
@d16 2
a17 2
;;(proclaim '(optimize (speed 0) (debug 3)))
(proclaim '(optimize (speed 3) (debug 0)))
@


1.1
log
@Initial revision
@
text
@d2 1
a2 1
	$Id: svpoly.lisp,v 1.6 1997/12/25 02:00:04 marek Exp $	
d16 2
a17 1
(proclaim '(optimize (speed 0) (debug 3)))
@