[1] | 1 | #|
|
---|
| 2 | $Id: division.lisp,v 1.4 2009/01/22 04:00:56 marek Exp $
|
---|
| 3 | *--------------------------------------------------------------------------*
|
---|
| 4 | | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu) |
|
---|
| 5 | | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
|
---|
| 6 | | |
|
---|
| 7 | | Everyone is permitted to copy, distribute and modify the code in this |
|
---|
| 8 | | directory, as long as this copyright note is preserved verbatim. |
|
---|
| 9 | *--------------------------------------------------------------------------*
|
---|
| 10 | |#
|
---|
| 11 | (defpackage "DIVISION"
|
---|
| 12 | (:use "MONOM" "ORDER" "TERM" "POLY" "COEFFICIENT-RING" "COMMON-LISP")
|
---|
| 13 | (:export divide poly-exact-divide))
|
---|
| 14 |
|
---|
| 15 | (in-package "DIVISION")
|
---|
| 16 |
|
---|
| 17 | #+debug(proclaim '(optimize (speed 0) (debug 3)))
|
---|
| 18 | #-debug(proclaim '(optimize (speed 3) (debug 0)))
|
---|
| 19 |
|
---|
| 20 | (defun divide (f fl &optional
|
---|
| 21 | (pred #'lex>)
|
---|
| 22 | (ring *coefficient-ring*)
|
---|
| 23 | &aux (s (length fl)))
|
---|
| 24 | "Divide polynomial F by a list of polynomials FL; use predicate PRED
|
---|
| 25 | to sort monomials; assumes that the polynomials have initially been
|
---|
| 26 | sorted according to PRED. It returnes multiple values. The first value
|
---|
| 27 | is a list of quotients A. The second value is the remainder R. These
|
---|
| 28 | object satisfy the quation F = SUM A[J]*FL[I] + R."
|
---|
| 29 | (do ((a (make-list s))
|
---|
| 30 | r
|
---|
| 31 | (p f)
|
---|
| 32 | (division-occurred nil nil))
|
---|
| 33 | ((endp p) (values (mapcar #'reverse a) (reverse r)))
|
---|
| 34 | (declare (list a r p))
|
---|
| 35 | (do ((fl fl (rest fl))
|
---|
| 36 | (a a (rest a)))
|
---|
| 37 | ((or (endp fl) division-occurred))
|
---|
| 38 | (when (term-divides-p (car (first fl)) (first p))
|
---|
| 39 | (let ((quot (term/ (first p) (car (first fl)) ring)))
|
---|
| 40 | (push quot (car a))
|
---|
| 41 | (setf p (poly-op (rest p) quot (rest (first fl)) pred ring)
|
---|
| 42 | division-occurred t))))
|
---|
| 43 | (when (not division-occurred)
|
---|
| 44 | (setf r (cons (first p) r)
|
---|
| 45 | p (rest p)))))
|
---|
| 46 |
|
---|
| 47 | (defun poly-exact-divide (f g &optional (order #'lex>) (ring *coefficient-ring*))
|
---|
| 48 | "Divide a polynomial F by another polynomial G. Assume that exact division
|
---|
| 49 | with no remainder is possible. Returns the quotient."
|
---|
| 50 | (multiple-value-bind (q r)
|
---|
| 51 | (divide f (list g) order ring)
|
---|
| 52 | (unless (endp r) (error "Exact division failed."))
|
---|
| 53 | (car q)))
|
---|