source: CGBLisp/trunk/src/division.lisp@ 14

#|
        $Id: division.lisp,v 1.4 2009/01/22 04:00:56 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 "DIVISION"
  (:use "MONOM" "ORDER" "TERM" "POLY" "COEFFICIENT-RING" "COMMON-LISP")
  (:export divide poly-exact-divide))

(in-package  "DIVISION")

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

(defun divide (f fl &optional
                    (pred #'lex>)
                    (ring *coefficient-ring*)
                 &aux (s (length fl)))
  "Divide polynomial F by a list of polynomials FL; use predicate PRED
to sort monomials; assumes that the polynomials have initially been
sorted according to PRED. It returnes multiple values. The first value
is a list of quotients A.  The second value is the remainder R. These
object satisfy the quation F = SUM A[J]*FL[I] + R."
  (do ((a (make-list s))
       r 
       (p f)
       (division-occurred nil nil))
      ((endp p) (values (mapcar #'reverse a) (reverse r)))
    (declare (list a r p))
    (do ((fl fl (rest fl))
         (a a (rest a)))
        ((or (endp fl) division-occurred))
      (when (term-divides-p (car (first fl)) (first p))
        (let ((quot (term/ (first p) (car (first fl)) ring)))
          (push quot (car a))
          (setf p (poly-op (rest p) quot (rest (first fl)) pred ring)
                division-occurred t))))
    (when (not division-occurred)
      (setf r (cons (first p) r)
            p (rest p))))) 
        
(defun poly-exact-divide (f g  &optional (order #'lex>) (ring *coefficient-ring*))
  "Divide a polynomial F by another polynomial G. Assume that exact division
with no remainder is possible. Returns the quotient."
  (multiple-value-bind (q r)
      (divide f (list g) order ring)
    (unless (endp r) (error "Exact division failed."))
    (car q)))