1 | #|
|
---|
2 | $Id: ratpoly.lisp,v 1.4 2009/01/22 04:07:33 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 | ;; Operations in the ring k(x2,x3,...,xn)[x1]
|
---|
12 | ;; The representation is (n1 a1 n2 a2 ... nk ak)
|
---|
13 | ;; where ak are rational functions in an arbitrary number of variables
|
---|
14 | ;; as in the "rat" package; this representation is referred to as RATPOLY
|
---|
15 | ;;
|
---|
16 | ;; Another representation is POLY1 and this corresponds to reppresenting
|
---|
17 | ;; polynomials as elements of k[x2,x3,...,xn][x1]
|
---|
18 | ;;
|
---|
19 | ;; In order to facilitate input, a function called
|
---|
20 | ;; poly-to-ratpoly is provided; it converts a polynomial in alist form
|
---|
21 | ;; to a form like above (plist form) with ai being rational
|
---|
22 | ;; Output is done by function ratpoly-print
|
---|
23 | ;;
|
---|
24 |
|
---|
25 | (defpackage "RATPOLY"
|
---|
26 | (:export ratpoly+ ratpoly- ratpoly* ratpoly-uminus scalar-times-ratpoly
|
---|
27 | rat-times-ratpoly ratpoly-divide ratpoly-remainder
|
---|
28 | ratpoly-gcd ratpoly-diff ratpoly-square-free
|
---|
29 | ratpoly-normalize ratpoly-resultant
|
---|
30 | deg lead ratpoly-discriminant ratpoly-print
|
---|
31 | poly-to-ratpoly ratpoly-to-poly poly-resultant)
|
---|
32 | (:use "RAT" "PRINTER" "MONOM" "DIVISION" "TERM" "POLY" "COMMON-LISP"))
|
---|
33 |
|
---|
34 | (in-package "RATPOLY")
|
---|
35 |
|
---|
36 | (proclaim '(optimize (speed 3) (space 0) (safety 0) (debug 0)))
|
---|
37 |
|
---|
38 | ;; Arithmetic on polynomials in one variable
|
---|
39 | ;; The polynomial a1 * x^n1 + a2 * x^n2 + ... + ak * x^nk
|
---|
40 | ;; where ai are nonzero numbers and n1>n2>...>nk is represented as list
|
---|
41 | ;; (n1 a1 n2 a2 ... nk ak)
|
---|
42 | ;; This is the sparse recursive representation
|
---|
43 |
|
---|
44 | (defun ratpoly+ (p q)
|
---|
45 | "Add polynomials P and Q."
|
---|
46 | (cond ((endp p) q)
|
---|
47 | ((endp q) p)
|
---|
48 | ((= (car p) (car q))
|
---|
49 | (let ((s (rat+ (cadr p) (cadr q))))
|
---|
50 | (if (rat-zerop s) ;check for cancellation
|
---|
51 | (ratpoly+ (cddr p) (cddr q))
|
---|
52 | (append (list (car p) s)
|
---|
53 | (ratpoly+ (cddr p) (cddr q))))))
|
---|
54 | ((> (car p) (car q))
|
---|
55 | (append (list (car p) (cadr p))
|
---|
56 | (ratpoly+ (cddr p) q)))
|
---|
57 | (t (append (list (car q) (cadr q))
|
---|
58 | (ratpoly+ p (cddr q))))))
|
---|
59 |
|
---|
60 | (defun ratpoly- (p q)
|
---|
61 | ;;"Subtract polynomial Q from P."
|
---|
62 | (cond ((endp p) (ratpoly-uminus q))
|
---|
63 | ((endp q) p)
|
---|
64 | ((= (car p) (car q))
|
---|
65 | (let ((s (rat- (cadr p) (cadr q))))
|
---|
66 | (if (rat-zerop s) ;check for cancellation
|
---|
67 | (ratpoly- (cddr p) (cddr q))
|
---|
68 | (append (list (car p) s)
|
---|
69 | (ratpoly- (cddr p) (cddr q))))))
|
---|
70 | ((> (car p) (car q))
|
---|
71 | (append (list (car p) (cadr p))
|
---|
72 | (ratpoly- (cddr p) q)))
|
---|
73 | (t (append (list (car q) (rat-uminus (cadr q)))
|
---|
74 | (ratpoly- p (cddr q))))))
|
---|
75 |
|
---|
76 | (defun ratpoly-uminus (p)
|
---|
77 | (cond
|
---|
78 | ((endp p) nil)
|
---|
79 | (t (cons (car p) (cons (rat-uminus (cadr p)) (ratpoly-uminus (cddr p)))))))
|
---|
80 |
|
---|
81 | (defun ratpoly* (p q)
|
---|
82 | "Multiply polynomials P and Q."
|
---|
83 | (cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)
|
---|
84 | ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
|
---|
85 | (t (append (list (+ (car p) (car q)) (rat* (cadr p) (cadr q)))
|
---|
86 | (ratpoly+ (ratpoly* (list (car p) (cadr p)) (cddr q))
|
---|
87 | (ratpoly* (cddr p) q))))))
|
---|
88 |
|
---|
89 | (defun scalar-times-ratpoly (scalar p)
|
---|
90 | "Multiply scalar SCALAR by a polynomial P."
|
---|
91 | (cond ((endp p) nil)
|
---|
92 | (t (cons (car p)
|
---|
93 | (cons (scalar-times-rat scalar (cadr p))
|
---|
94 | (scalar-times-ratpoly scalar (cddr p)))))))
|
---|
95 |
|
---|
96 | (defun rat-times-ratpoly (scalar p)
|
---|
97 | "Multiply rational function SCALAR by a polynomial P."
|
---|
98 | (cond ((endp p) nil)
|
---|
99 | (t (cons (car p)
|
---|
100 | (cons (rat* scalar (cadr p))
|
---|
101 | (rat-times-ratpoly scalar (cddr p)))))))
|
---|
102 |
|
---|
103 | (defun ratpoly-divide (f g)
|
---|
104 | "Divide polynomial F by G. Return quotient and remainder as multiple values."
|
---|
105 | (do* ((lp nil (- (car r) (car g)))
|
---|
106 | (lc nil (rat/ (cadr r) (cadr g)))
|
---|
107 | (q nil (nconc q (list lp lc)))
|
---|
108 | (r f (ratpoly- (cddr r) (ratpoly* (list lp lc) (cddr g)))))
|
---|
109 | ((or (endp r) (< (car r) (car g)))
|
---|
110 | (values q r))))
|
---|
111 |
|
---|
112 | (defun ratpoly-remainder (f g)
|
---|
113 | "The remainder of the division of a polynomial F by G."
|
---|
114 | (second (multiple-value-list (ratpoly-divide f g))))
|
---|
115 |
|
---|
116 | (defun ratpoly-gcd (f g)
|
---|
117 | "Return GCD of polynomials F and G."
|
---|
118 | (do ((h f s)
|
---|
119 | (s g (ratpoly-remainder h s)))
|
---|
120 | ((endp s) h)))
|
---|
121 |
|
---|
122 | (defun ratpoly-diff (f)
|
---|
123 | "Differentiate a polynomial."
|
---|
124 | (cond ((endp f) nil)
|
---|
125 | ((zerop (car f)) nil) ;degree 0
|
---|
126 | (t (append (list (1- (car f)) (scalar-times-rat (car f) (cadr f)))
|
---|
127 | (ratpoly-diff (cddr f))))))
|
---|
128 |
|
---|
129 | (defun ratpoly-square-free (f)
|
---|
130 | "Return the square-free part of a polynomial F."
|
---|
131 | (ratpoly-divide f (ratpoly-gcd f (ratpoly-diff f))))
|
---|
132 |
|
---|
133 | (defun ratpoly-normalize (f)
|
---|
134 | "Divide a non-zero polynomial by the coefficient at the highest power."
|
---|
135 | (rat-times-ratpoly (cons (denom (cadr f)) (num (cadr f))) f))
|
---|
136 |
|
---|
137 | ;; A modification of Euclidean algorithm
|
---|
138 | (defun ratpoly-resultant (f g)
|
---|
139 | "Return the resultant of polynomials F and G."
|
---|
140 | (do* ((r nil (ratpoly-remainder h s))
|
---|
141 | (res (list 0 (rat-constant 1 (length (car (caaadr f)))))
|
---|
142 | (rat-times-ratpoly
|
---|
143 | (scalar-times-rat
|
---|
144 | (expt -1 (* (deg h) (deg s)))
|
---|
145 | (rat-expt (lead s) (- (deg h) (deg r))))
|
---|
146 | res))
|
---|
147 | (h f s)
|
---|
148 | (s g r))
|
---|
149 | ((<= (deg s) 0)
|
---|
150 | (cond
|
---|
151 | ((or (endp s) (endp h)) nil)
|
---|
152 | ((plusp (deg h))
|
---|
153 | (rat-times-ratpoly (rat-expt (lead s) (deg h))
|
---|
154 | res))))))
|
---|
155 |
|
---|
156 | (defun deg (s)
|
---|
157 | (cond
|
---|
158 | ((endp s)
|
---|
159 | #+debugging(warn "ratpoly::deg: Calculating degree of 0 polynomial")
|
---|
160 | -1)
|
---|
161 | (t (car s))))
|
---|
162 |
|
---|
163 | (defun lead (s) (cadr s))
|
---|
164 |
|
---|
165 | (defun ratpoly-discriminant (p &aux (l (deg p)))
|
---|
166 | "The discriminant of a polynomial P."
|
---|
167 | (rat-times-ratpoly
|
---|
168 | (scalar-div-rat (expt -1 (mod (* l (1- l)) 2)) (lead p))
|
---|
169 | (ratpoly-resultant p (ratpoly-diff p))))
|
---|
170 |
|
---|
171 | (defun ratpoly-print (p vars &optional (stream t) (beg t) (p-orig p))
|
---|
172 | (when (endp p)
|
---|
173 | (when beg (format stream "0"))
|
---|
174 | (return-from ratpoly-print p-orig))
|
---|
175 | (if beg (format stream "((")
|
---|
176 | (format stream " + (("))
|
---|
177 | (poly-print (caadr p) (rest vars) stream)
|
---|
178 | (format stream ") / (")
|
---|
179 | (poly-print (cdadr p) (rest vars) stream)
|
---|
180 | (case (car p)
|
---|
181 | (1 (format stream ")) * ~a" (car vars)))
|
---|
182 | (0 (format stream "))"))
|
---|
183 | (otherwise
|
---|
184 | (format stream ")) * ~a^~d" (car vars) (car p))))
|
---|
185 | (ratpoly-print (cddr p) vars stream nil p-orig))
|
---|
186 |
|
---|
187 | (defun poly-to-ratpoly (p)
|
---|
188 | (poly1-to-ratpoly (poly-to-poly1 p)))
|
---|
189 |
|
---|
190 | ;; Convert a polynomial to a polynomial in k[x2,...,xn][x1]
|
---|
191 | (defun poly-to-poly1 (p &aux (htab (make-hash-table)) q)
|
---|
192 | (dolist (term p)
|
---|
193 | (push (cons (cdar term) (cdr term))
|
---|
194 | (gethash (caar term) htab nil)))
|
---|
195 | (maphash #'(lambda (key val) (push (cons key val) q)) htab)
|
---|
196 | (mapcan #'(lambda (x &aux (deg (car x)) (coef (cdr x)))
|
---|
197 | (list deg
|
---|
198 | (sort-poly coef)))
|
---|
199 | (sort q #'> :key #'car)))
|
---|
200 |
|
---|
201 |
|
---|
202 | ;; Convert poly1 to ratpoly, i.e. add denominators=1
|
---|
203 | (defun poly1-to-ratpoly (p)
|
---|
204 | (unless (endp p)
|
---|
205 | (cons (car p)
|
---|
206 | (cons
|
---|
207 | (cons (cadr p)
|
---|
208 | (list
|
---|
209 | (cons (make-list (length (caaadr p)) :initial-element 0) 1)))
|
---|
210 | (poly1-to-ratpoly (cddr p))))))
|
---|
211 |
|
---|
212 | (defun ratpoly-to-poly1 (p)
|
---|
213 | "Convert every coefficient of ratpoly to polynomial if possible"
|
---|
214 | (cond
|
---|
215 | ((endp p) nil)
|
---|
216 | (t (cons (car p)
|
---|
217 | (cons (rat-to-poly (cadr p))
|
---|
218 | (ratpoly-to-poly1 (cddr p)))))))
|
---|
219 |
|
---|
220 | (defun poly1-to-poly (p)
|
---|
221 | "Convert a ratpoly, whose coeffs have been converted to poly,
|
---|
222 | into a poly structure, i.e. tack in powers of first variable."
|
---|
223 | (cond
|
---|
224 | ((endp p) nil)
|
---|
225 | (t
|
---|
226 | (append
|
---|
227 | (mapcar #'(lambda (x) (cons (cons (car p) (car x)) (cdr x))) (cadr p))
|
---|
228 | (poly1-to-poly (cddr p))))))
|
---|
229 |
|
---|
230 | (defun ratpoly-to-poly (p)
|
---|
231 | (poly1-to-poly (ratpoly-to-poly1 p)))
|
---|
232 |
|
---|
233 |
|
---|
234 | (defun poly-resultant (f g)
|
---|
235 | "Calculate resultant of F and G given in poly i.e. alist representation."
|
---|
236 | (ratpoly-to-poly
|
---|
237 | (ratpoly-resultant (poly-to-ratpoly f)
|
---|
238 | (poly-to-ratpoly g))))
|
---|
239 |
|
---|
240 | #|
|
---|
241 | ;;----------------------------------------------------------------
|
---|
242 | ;; Multi-variable GCD algorithm
|
---|
243 | ;; Roughly p. 134 of Davenport, Siret, Tournier
|
---|
244 | ;;----------------------------------------------------------------
|
---|
245 | (defun poly-gcd (A B)
|
---|
246 | (multiple-value-bind (Ap Ac)
|
---|
247 | (primitive-part (poly-to-poly1 A))
|
---|
248 | (multiple-value-bind (Bp Bc)
|
---|
249 | (primitive-part (poly-to-poly1 B))
|
---|
250 | (poly* (poly1-to-poly (primitive-part (euclid Ap Bp)))
|
---|
251 | (poly1-to-poly
|
---|
252 | (list 0
|
---|
253 | (cond
|
---|
254 | ((numberp Ac) (list (cons nil (gcd Ac Bc))))
|
---|
255 | (t (poly-gcd Ac Bc)))))))))
|
---|
256 |
|
---|
257 | ;; This operates on poly1
|
---|
258 | (defun content1 (A)
|
---|
259 | (cond
|
---|
260 | ((endp A) (error "content1: Content of 0 is not defined."))
|
---|
261 | ((endp (caaadr A)) ;1-variable poly
|
---|
262 | (content1-aux A)) ;a number
|
---|
263 | ((endp (cddr A)) (cadr A))
|
---|
264 | (t (poly-gcd (cadr A) (content1 (cddr A))))))
|
---|
265 |
|
---|
266 | ;; 1-variable case
|
---|
267 | (defun content1-aux (A)
|
---|
268 | (cond ((endp A) (error "content1-aux: Content of 0 is not defined."))
|
---|
269 | ((endp (cddr A)) (cdaadr A))
|
---|
270 | (t (gcd (cdaadr A) (content1-aux (cddr A))))))
|
---|
271 |
|
---|
272 | ;; operates on poly
|
---|
273 | (defun content (A &aux (A1 (poly-to-poly1 A)))
|
---|
274 | (content1 A1))
|
---|
275 |
|
---|
276 | ;; Operates on poly1
|
---|
277 | (defun primitive-part (A)
|
---|
278 | (let ((Ac (content1 A)))
|
---|
279 | (values (divide-coeffs A Ac) Ac)))
|
---|
280 |
|
---|
281 | ;; Operates on A in poly1 form and A in poly form of 1 variable less
|
---|
282 | (defun divide-coeffs (A Ac)
|
---|
283 | (cond
|
---|
284 | ((endp A) nil)
|
---|
285 | ((numberp Ac) ;ground case
|
---|
286 | (cons (car A)
|
---|
287 | (cons (list (cons nil (floor (cdaadr A) Ac)))
|
---|
288 | (divide-coeffs (cddr A) Ac))))
|
---|
289 | (t (cons (car A)
|
---|
290 | (cons
|
---|
291 | (multiple-value-bind (q r)
|
---|
292 | (divide (cadr A) (list Ac))
|
---|
293 | (unless (endp r) (error "divide-coeffs: not divisible."))
|
---|
294 | (car q))
|
---|
295 | (divide-coeffs (cddr A) Ac))))))
|
---|
296 |
|
---|
297 | ;; Euclid operates on poly1
|
---|
298 | (defun euclid (A B)
|
---|
299 | (ratpoly-to-poly1 (ratpoly-gcd (poly1-to-ratpoly A)
|
---|
300 | (poly1-to-ratpoly B))))
|
---|
301 |
|
---|
302 | |#
|
---|
303 |
|
---|