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 | (proclaim '(optimize (speed 0) (debug 3)))
|
---|
18 |
|
---|
19 | (defun divide (f fl &optional
|
---|
20 | (pred #'lex>)
|
---|
21 | (ring *coefficient-ring*)
|
---|
22 | &aux (s (length fl)))
|
---|
23 | "Divide polynomial F by a list of polynomials FL; use predicate PRED
|
---|
24 | to sort monomials; assumes that the polynomials have initially been
|
---|
25 | sorted according to PRED. It returnes multiple values. The first value
|
---|
26 | is a list of quotients A. The second value is the remainder R. These
|
---|
27 | object satisfy the quation F = SUM A[J]*FL[I] + R."
|
---|
28 | (do ((a (make-list s))
|
---|
29 | r
|
---|
30 | (p f)
|
---|
31 | (division-occurred nil nil))
|
---|
32 | ((endp p) (values (mapcar #'reverse a) (reverse r)))
|
---|
33 | (declare (list a r p))
|
---|
34 | (do ((fl fl (rest fl))
|
---|
35 | (a a (rest a)))
|
---|
36 | ((or (endp fl) division-occurred))
|
---|
37 | (when (term-divides-p (car (first fl)) (first p))
|
---|
38 | (let ((quot (term/ (first p) (car (first fl)) ring)))
|
---|
39 | (push quot (car a))
|
---|
40 | (setf p (poly-op (rest p) quot (rest (first fl)) pred ring)
|
---|
41 | division-occurred t))))
|
---|
42 | (when (not division-occurred)
|
---|
43 | (setf r (cons (first p) r)
|
---|
44 | p (rest p)))))
|
---|
45 |
|
---|
46 | (defun poly-exact-divide (f g &optional (order #'lex>) (ring *coefficient-ring*))
|
---|
47 | "Divide a polynomial F by another polynomial G. Assume that exact division
|
---|
48 | with no remainder is possible. Returns the quotient."
|
---|
49 | (multiple-value-bind (q r)
|
---|
50 | (divide f (list g) order ring)
|
---|
51 | (unless (endp r) (error "Exact division failed."))
|
---|
52 | (car q)))
|
---|