source: CGBLisp/trunk/src/poly-gcd.lisp@ 14

#|
        $Id: poly-gcd.lisp,v 1.4 2009/01/22 04:05:32 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 "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")

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

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