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


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

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

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

1.1
date	2009.01.19.07.00.23;	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 "POLY-GCD"
  (:export poly-gcd poly-content poly-pseudo-divide poly-primitive-part
		   poly-pseudo-remainder)
  (:use "ORDER" "POLY" "DIVISION" "COMMON-LISP"))

(in-package "POLY-GCD")

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

;; This package calculates GCD of polynomials over integers.
;; They are assumed to be ordered lexicographically. 
;;
;; The algorithm is that on p. 57 of Geddes, Czapor, Labahn
;; Given polynomials a(x),b(x) in D[x] where D is a UFD, we
;; compute g(x)=GCD(a(x),b(x))
;; Assume that the poly's are sorted lexicographically

(defun poly-gcd (a b)
  (cond
   ((endp a) b)
   ((endp b) a)
   ((endp (caar a))			;scalar
    (list (cons nil (gcd (cdar a) (cdar b)))))
   (t
    (do* ((r nil (poly-pseudo-remainder c d))
	  (c (poly-primitive-part a) d)
	  (d (poly-primitive-part b) (poly-primitive-part r)))
	((endp d)
	 (poly* (poly-extend 
		 (poly-gcd (poly-content a) (poly-content b)))
		c))))))
	 
;; Perform a pseudo-division in k[x2,....,xn][x1]
;; Assume that the terms are sorted according to decreasing powers of x1
;; p.297 of Cox
(defun poly-pseudo-divide (f g)	
  (multiple-value-bind (lg grest)
      (lpart g)
    (do* ((m (mdeg g))
	  (dm (poly-extend lg))
	  (lpart nil (multiple-value-list (lpart r)))
	  (lr nil (car lpart))
	  (lrest nil (cadr lpart))
	  (term nil (poly-extend lr (list (- (mdeg r) m))))
	  (r f (poly- (poly* dm lrest)
		      (poly* term grest)))
	  (q nil (append (poly* dm q) term)))
	((or (endp r) (< (mdeg r) m))
	 (values q r)))))


(defun poly-pseudo-remainder (f g)
  (second (multiple-value-list (poly-pseudo-divide f g))))

;; Degree in main variable
(defun mdeg (b) (caaar b))

;; Leading coefficient in the first variable; a poly in k[x2,...,xn]
(defun lcoeff (b)
  (first (multiple-value-list (lpart b))))

(defun lrest (b)
  (second (multiple-value-list (lpart b))))

(defun lpart (b)
  (do ((mdeg (mdeg b))
       (b b (rest b))
       (b1 nil (cons  (cons (cdaar b) (cdar b)) b1)))
      ((or (endp b) (/= (caaar b) mdeg))
       (values (reverse b1) b))))

;; Divide f by its content
(defun poly-primitive-part (f)
  (cond
   ((endp (caar f))			;scalar
    f)
   (t
    (poly-exact-divide f (poly-extend (poly-content f))))))


(defun poly-content (f)
  (cond
   ((endp f) (error "Zero argument"))
   (t (multiple-value-bind (lc r)
	  (lpart f)
	(cond
	 ((endp r) lc)
	 ((and (= (length lc) 1)
	       (every #'zerop (caar lc))
	       (or (= (cdar lc) 1) 
		   (= (cdar lc) -1)))	;lc is 1 or -1
	  lc)
	 (t (poly-gcd lc (poly-content r))))))))
@


1.3
log
@*** empty log message ***
@
text
@d18 2
a19 2
;;(proclaim '(optimize (speed 0) (debug 3)))
(proclaim '(optimize (speed 3) (debug 0)))
@


1.2
log
@*** empty log message ***
@
text
@d18 2
a19 1
(proclaim '(optimize (speed 0) (debug 3)))
@


1.1
log
@Initial revision
@
text
@d2 1
a2 1
	$Id: poly-gcd.lisp,v 1.12 1997/12/13 07:11:59 marek Exp $	
d18 1
@