1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1981 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
13 (macsyma-module rat3c
)
15 ;; THIS IS THE NEW RATIONAL FUNCTION PACKAGE PART 3.
16 ;; IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS
20 ((and $keepfloat
(or (floatp a
) (floatp b
))) 1)
23 (defun pquotientchk (a b
)
24 (if (equal b
1) a
(pquotient a b
)))
26 ;; divides polynomial x by polynomial y
27 ;; avoids error "quotient by polynomial of higher degree"
28 ;; (returns nil in this case)
29 (defun pquotientchk-safe (x y
)
30 (ignore-rat-err (pquotientchk x y
)))
32 (defun ptimeschk (a b
)
38 (catch 'float
(if (pcoefp x
) (floatp x
) (pfloatp1 x
))))
41 (mapc #'(lambda (q) (cond ((pcoefp q
) (when (floatp q
) (throw 'float t
)))
47 (setq x
(car (pgcda x y nil
)))
48 (cond ((pminusp x
) (pminus x
))
53 (setq x
(pgcdcofacts x y
))
54 (ptimes (car x
) (ptimes (cadr x
) (caddr x
))))
56 (defun plcmcofacts (x y
)
57 (setq x
(pgcdcofacts x y
))
58 (list (ptimes (car x
) (ptimes (cadr x
) (caddr x
)))
61 ; returns list (gcd xx yy alg)
62 ; where x * y = gcd^2 * xx * yy / alg^2
63 ; and alg is non-nil only when $algebraic is true
64 (defun pgcdcofacts (x y
)
65 (let ((a (pgcda x y t
)))
67 ((equal (setq a
(car a
)) 1) (list 1 x y
))
68 ((and $algebraic
(not (pcoefp a
)))
69 (cons a
(prog2 (setq x
(rquotient x a
)
71 a
(pgcdcofacts (cdr x
) (cdr y
)))
72 (list (ptimes (car x
) (caddr a
))
73 (ptimes (car y
) (cadr a
))
74 (ptimes (cadr a
) (cdr y
))))))
75 ((eq a x
) (list x
1 (pquotient y x
)))
76 ((eq a y
) (list a
(pquotient x y
) 1))
77 (t (list a
(pquotient x a
) (pquotient y a
))))))
79 (defun pgcda (x y cofac?
&aux a c
)
80 (cond ((not $gcd
) (list 1 x y
))
81 ((and $keepfloat
(or (pfloatp x
) (pfloatp y
)))
82 (cond ((or (pcoefp x
) (pcoefp y
)
83 (pcoefp (setq a
(car (ptermcont x
))))
84 (pcoefp (setq a
(pgcd a
(car (ptermcont y
))))))
89 (cons (setq a
(cgcd x y
))
91 (list (cquotient x a
) ;(CQUOTIENT 0 0) = 0
93 ((zerop x
) (list y x
1))
94 (t (list (pcontent1 (cdr y
) x
)))))
95 ((pcoefp y
) (cond ((zerop y
) (list x
1 y
))
96 (t (list (pcontent1 (cdr x
) y
)))))
97 ((equal x y
) (list x
1 1))
98 ($ratfac
(multiple-value-list (fpgcdco x y
)))
99 ((not (eq (p-var x
) (p-var y
)))
100 (list (if (pointergp (p-var x
) (p-var y
))
101 (oldcontent1 (p-terms x
) y
)
102 (oldcontent1 (p-terms y
) x
))))
103 ((progn (desetq (a x
) (ptermcont x
))
104 (desetq (c y
) (ptermcont y
))
105 (not (and (equal a
1) (equal c
1))))
106 (mapcar #'ptimes
(monomgcdco a c cofac?
) (pgcda x y cofac?
)))
107 ((and (not $algebraic
) (not modulus
)
108 (desetq (a . c
) (lin-var-find (nreverse (pdegreevector x
))
109 (nreverse (pdegreevector y
))
111 (cond ((= a
1) (linhack x y
(car c
) (cadr c
) cofac?
))
112 (t (setq a
(linhack y x a
(cadr c
) cofac?
))
113 (if (cdr a
) (rplacd a
(nreverse (cdr a
))))
115 ((eq $gcd
'$spmod
) (list (zgcd x y
)))
116 ((eq $gcd
'$subres
) (list (oldgcd x y
)))
117 ((eq $gcd
'$algebraic
)
118 (if (or (palgp x
) (palgp y
))
119 (let (($gcd
'$subres
)) (list (oldgcd x y
)))
120 (let (($gcd
'$spmod
)) (list (zgcd x y
)))))
121 ((eq $gcd
'$ez
) (ezgcd2 x y
))
122 ((eq $gcd
'$red
) (list (oldgcd x y
)))
123 ((eq $gcd
'$mod
) (newgcd x y modulus
))
124 ((not (member $gcd
*gcdl
* :test
#'eq
))
125 (merror (intl:gettext
"gcd: 'gcd' variable must be one of ~M; found: ~M") *gcdl
* $gcd
))
128 (defun monomgcdco (p q cofac?
)
129 (let ((gcd (monomgcd p q
)))
130 (cons gcd
(if cofac?
(list (pquotient p gcd
) (pquotient q gcd
)) ()))))
132 (defun monomgcd (p q
)
133 (cond ((or (pcoefp p
) (pcoefp q
)) 1)
134 ((eq (p-var p
) (p-var q
))
135 (make-poly (p-var p
) (min (p-le p
) (p-le q
))
136 (monomgcd (p-lc p
) (p-lc q
))))
137 ((pointergp (car p
) (car q
)) (monomgcd (p-lc p
) q
))
138 (t (monomgcd p
(p-lc q
)))))
140 (defun linhack (pol1 pol2 nonlindeg var cofac?
)
141 (prog (coeff11 coeff12 gcdab rpol1 rpol2 gcdcd gcdcoef
)
142 (desetq (coeff11 . coeff12
) (bothprodcoef (make-poly var
) pol1
))
143 (setq gcdab
(if (pzerop coeff12
) coeff11
144 (pgcd coeff11 coeff12
)))
145 (cond ((equal gcdab
1)
146 (cond ((setq coeff11
(testdivide pol2 pol1
))
147 (return (list pol1
1 coeff11
)))
148 (t (return (list 1 pol1 pol2
))))))
149 (setq rpol1
(pquotient pol1 gcdab
))
150 (desetq (gcdcd rpol2
) (linhackcontent var pol2 nonlindeg
))
151 (cond ((equal gcdcd
1)
152 (cond ((setq coeff12
(testdivide rpol2 rpol1
))
153 (return (list rpol1 gcdab coeff12
)))
154 (t (return (list 1 pol1 pol2
))))))
155 (cond (cofac?
(desetq (gcdcoef coeff11 coeff12
)
156 (pgcdcofacts gcdab gcdcd
))
157 (cond ((setq gcdcd
(testdivide* rpol2 rpol1
))
158 (return (list (ptimes gcdcoef rpol1
)
160 (ptimes coeff12 gcdcd
))))
161 (t (return (list gcdcoef
162 (ptimes coeff11 rpol1
)
163 (ptimes coeff12 rpol2
))))))
164 (t (setq gcdcoef
(pgcd gcdcd gcdab
))
165 (cond ((testdivide* rpol2 rpol1
)
166 (return (list (ptimes gcdcoef rpol1
))))
167 (t (return (list gcdcoef
))))))))
169 (defun lin-var-find (a b c
)
170 (do ((varl c
(cdr varl
))
171 (degl1 a
(cdr degl1
))
172 (degl2 b
(cdr degl2
)))
173 ((or (null degl1
) (null degl2
)) nil
)
174 (if (equal (min (car degl1
) (car degl2
)) 1)
175 (return (list (car degl1
) (car degl2
) (car varl
))))))
177 (defun linhackcontent (var pol nonlindeg
&aux
(npol pol
) coef gcd
)
178 (do ((i nonlindeg
(1- i
)))
179 ((= i
0) (list (setq gcd
(pgcd gcd npol
)) (pquotient pol gcd
)))
180 (desetq (coef . npol
) (bothprodcoef (make-poly var i
1) npol
))
181 (unless (pzerop coef
)
182 (setq gcd
(if (null gcd
) coef
(pgcd coef gcd
)))
183 (if (equal gcd
1) (return (list 1 pol
))))))
185 ;;*** THIS IS THE REDUCED POLYNOMIAL REMAINDER SEQUENCE GCD (COLLINS')
187 (defun oldgcd (x y
&aux u v s egcd
) ;only called from pgcda
188 (desetq (x u
) (oldcontent x
))
189 (desetq (y v
) (oldcontent y
))
190 (setq egcd
(gcd (pgcdexpon u
) (pgcdexpon v
)))
192 (setq u
(pexpon*// u egcd nil
)
193 v
(pexpon*// v egcd nil
)))
194 (if (> (p-le v
) (p-le u
)) (rotatef u v
))
197 ($subres
(subresgcd u v
))
198 (t (merror "OLDGCD: found gcd = ~M; how did that happen?" $gcd
))))
199 ;; Check for gcd that simplifies to 0. SourceForge bugs 831445 and 1313987
200 (unless (ignore-rat-err (rainv s
))
203 (setq s
(pexpon*// (primpart
205 (pquotient s
(pquotient (p-lc s
)
206 (pgcd (p-lc u
) (p-lc v
))))))
208 (setq s
(ptimeschk s
(pgcd x y
)))
209 (and $algebraic
(not (pcoefp (setq u
(leadalgcoef s
))))
210 (not (equal u s
)) (setq s
(algnormal s
)))
211 (cond (modulus (monize s
))
212 ((pminusp s
) (pminus s
))
217 (do ((d (cadr p
) (gcd d
(car l
)))
218 (l (cdddr p
) (cddr l
)))
219 ((or (null l
) (= d
1)) d
))))
221 (defun pexpon*// (p n
*?
)
222 (if (or (pcoefp p
) (= n
1)) p
223 (do ((ans (list (car p
))
225 (cons (if *?
(* (car l
) n
)
226 (truncate (car l
) n
))
228 (l (cdr p
) (cddr l
)))
229 ((null l
) (nreverse ans
)))))
231 ;;polynomial gcd using reduced prs
233 (defun redgcd (p q
&aux
(d 0))
234 (loop until
(zerop (pdegree q
(p-var p
)))
236 q
(pquotientchk-safe (prem p q
) (pexpt (p-lc p
) d
))
237 d
(+ (p-le p
) 1 (- (p-le q
))))
238 (if (< d
1) (return 1))
239 finally
(return (if (pzerop q
) p
1))))
241 ;;computes gcd's using subresultant prs
242 ;;ACM Transactions On Mathematical Software Sept. 1978
244 (defun subresgcd (p q
)
245 (loop for g
= 1 then
(p-lc p
)
246 for h
= 1 then
(pquotientchk-safe (pexpt g d
) h^
1-d
)
247 for d
= (- (p-le p
) (p-le q
))
248 for h^
1-d
= 1 then
(if (< d
1)
252 q
(pquotientchk-safe (prem p q
) (ptimes g
(ptimes h h^
1-d
))))
253 if
(zerop (pdegree q
(p-var p
))) return
(if (pzerop q
) p
1)))
255 ;;*** THIS COMPUTES PSEUDO REMAINDERS
257 (defun psquorem1 (u v quop
)
258 (prog (k (m 0) lcu lcv quo lc
)
259 (declare (special lcu lcv
))
261 (setq k
(- (pt-le u
) (pt-le v
)))
262 (cond ((minusp k
) (return (list 1 '(0 0) u
))))
263 (if quop
(setq lc
(pexpt (pt-lc v
) (1+ k
))))
264 a
(setq lcu
(pminus (pt-lc u
)))
265 (if quop
(setq quo
(cons (ptimes (pt-lc u
) (pexpt (pt-lc v
) k
))
267 (cond ((null (setq u
(pgcd2 (pt-red u
) (pt-red v
) k
)))
268 (return (list lc
(nreverse quo
) '(0 0))))
269 ((minusp (setq m
(- (pt-le u
) (pt-le v
))))
270 (setq u
(cond ((zerop k
) u
)
271 (t (pctimes1 (pexpt lcv k
) u
))))
272 (return (list lc
(nreverse quo
) u
)))
274 (setq u
(pctimes1 (pexpt lcv
(- (1- k
) m
)) u
))))
281 (if (or modulus
(floatp p
) (floatp q
))
286 (t (psimp (p-var p
) (pgcd1 (p-terms p
) (p-terms q
))))))
288 (defun pgcd1 (u v
) (caddr (psquorem1 u v nil
)))
290 (defun pgcd2 (u v k
&aux
(i 0))
291 (declare (special lcu lcv
) (fixnum i
))
292 (cond ((null u
) (pcetimes1 v k lcu
))
293 ((null v
) (pctimes1 lcv u
))
294 ((zerop (setq i
(+ (pt-le u
) (- k
) (- (car v
)))))
295 (pcoefadd (pt-le u
) (pplus (ptimes lcv
(pt-lc u
))
296 (ptimes lcu
(pt-lc v
)))
297 (pgcd2 (pt-red u
) (pt-red v
) k
)))
299 (list* (+ (pt-le v
) k
) (ptimes lcu
(pt-lc v
)) (pgcd2 u
(pt-red v
) k
)))
300 (t (list* (pt-le u
) (ptimes lcv
(pt-lc u
)) (pgcd2 (pt-red u
) v k
)))))
302 ;;;*** OLDCONTENT REMOVES ALL BUT MAIN VARIABLE AND PUTS THAT IN CONTENT
303 ;;;*** OLDCONTENT OF 3*A*X IS 3*A (WITH MAINVAR=X)
305 (defun rcontent (p) ;RETURNS RAT-FORMS
306 (let ((q (oldcontenta p
)))
307 (list (cons q
1) (cond ($algebraic
(rquotient p q
))
308 (t (cons (pquotient p q
) 1))))))
310 (defun oldcontenta (x)
312 (t (setq x
(contsort (cdr x
)))
313 (oldcontent2 (cdr x
) (car x
)))))
315 (defun oldcontent (x)
316 (cond ((pcoefp x
) (list x
1))
318 (list (p-lc x
) (make-poly (p-var x
) (p-le x
) 1)))
319 (t (let ((u (contsort (cdr x
))) v
)
320 (setq u
(oldcontent2 (cdr u
) (car u
))
321 v
(cond ($algebraic
(car (rquotient x u
)))
322 (t (pcquotient x u
))))
323 (cond ((pminusp v
) (list (pminus u
) (pminus v
)))
326 (defun oldcontent1 (x gcd
)
327 (cond ((equal gcd
1) 1)
329 (t (oldcontent2 (contsort x
) gcd
))))
331 (defun oldcontent2 (x gcd
)
333 (gcd gcd
(pgcd (car x
) gcd
)))
334 ((or (null x
) (equal gcd
1)) gcd
)))
338 (cond ((member 1 x
) '(1))
340 (t (sort x
#'contodr
))))
344 (ans nil
(cons (cadr x
) ans
)))
350 ((eq (car a
) (car b
)) (not (> (cadr a
) (cadr b
))))
351 (t (pointergp (car b
)(car a
)))))
353 ;;;*** PCONTENT COMPUTES INTEGER CONTENT
354 ;;;*** PCONTENT OF 3*A*X IS 3 IF MODULUS = NIL 1 OTHERWISE
357 (cond ((pcoefp x
) (list x
1))
358 (t (let ((u (pcontentz x
)))
359 (if (equal u
1) (list 1 x
)
360 (list u
(pcquotient x u
)))))))
362 (defun pcontent1 (x gcd
)
364 (gcd gcd
(cgcd gcd
(pcontentz (cadr x
)))))
365 ((or (null x
) (equal gcd
1)) gcd
)))
369 (t (pcontent1 (p-red p
) (pcontentz (p-lc p
))))))
371 (defun ucontent (p) ;CONTENT OF UNIV. POLY
372 (cond ((pcoefp p
) (abs p
))
373 (t (setq p
(mapcar #'abs
(coefl (cdr p
))))
374 (let ((m (reduce #'min p
)))
375 (oldcontent2 (delete m p
:test
#'equal
) m
)))))
377 ;;*** PGCDU CORRESPONDS TO BROWN'S ALGORITHM U
379 ;;;PGCDU IS NOT NOW IN RAT;UFACT >
382 (do () ((pzerop q
) (monize p
))
383 (psetq p q q
(pmodrem p q
))))
386 (cond ((null modulus
)
387 (merror "PMODREM: null modulus; how did that happen?"))
388 ((pacoefp y
) (if (pzerop y
) x
0))
390 ((eq (p-var x
) (p-var y
))
391 (psimp (car x
) (pgcdu1 (p-terms x
) (p-terms y
) nil
)))
392 (t (merror "PMODREM: I can't handle this; x = ~M, y = ~M" x y
))))
394 (defun pmodquo (u v
&aux quo
)
395 (declare (special quo
))
396 (cond ((null modulus
)
397 (merror "PMODQUO: null modulus; how did that happen?"))
398 ((pcoefp v
) (cons (ptimes (crecip v
) u
) 0))
399 ((alg v
) (cons (ptimes (painvmod v
) u
) 0))
400 ((pacoefp u
) (cons 0 u
))
401 ((not (eq (p-var u
) (p-var v
)))
402 (merror "PMODQUO: arguments have different variables; how did that happen?"))
403 (t (xcons (psimp (car u
) (pgcdu1 (cdr u
) (cdr v
) t
))
404 (psimp (car u
) quo
)))))
407 (defun pgcdu1 (u v pquo
*)
408 (let ((invv (painvmod (pt-lc v
))) (k 0) q
*)
409 (declare (special k quo q
*))
410 (loop until
(minusp (setq k
(- (pt-le u
) (pt-le v
))))
411 do
(setq q
* (ptimes invv
(pt-lc u
)))
412 if pquo
* do
(setq quo
(nconc quo
(list k q
*)))
413 when
(ptzerop (setq u
(ptpt-subtract-powered-product
414 (pt-red u
) (pt-red v
) q
* k
)))
416 finally
(return u
))))
419 (do ((pl *bigprimes
* (cdr pl
)))
421 (setq p
(next-prime (1- p
) -
1))
422 (setq *bigprimes
* (nconc *bigprimes
* (list p
)))
427 (defun leadcoefficient (p)
428 (if (pcoefp p
) p
(leadcoefficient (caddr p
))))
430 (defun maxcoefficient (p)
431 (if (pcoefp p
) (abs p
) (maxcoef1 (cdr p
))))
434 (if (null p
) 0 (max (maxcoefficient (cadr p
)) (maxcoef1 (cddr p
)))))
436 (defun maxnorm (poly)
437 (if (null poly
) 0 (max (norm (cadr poly
)) (maxnorm (cddr poly
)))))
440 (cond ((null poly
) 0)
441 ((pcoefp poly
) (abs poly
))
442 (t (+ (norm (caddr poly
)) (norm1 (cdddr poly
)) )) ))
445 (if (null poly
) 0 (+ (norm (cadr poly
)) (norm1 (cddr poly
)) )) )
447 (defun pdegree (p var
)
449 ((eq var
(p-var p
)) (p-le p
))
450 ((pointergp var
(p-var p
)) 0)
451 (t (do ((l (p-red p
) (pt-red l
))
452 (e (pdegree (p-lc p
) var
) (max e
(pdegree (pt-lc l
) var
))))
455 (defun poly-in-var (p v
)
456 (cond ((or (pcoefp p
) (pointergp v
(p-var p
))) (list 0 p
))
457 ((eq (p-var p
) v
) (p-terms p
))
459 for
(exp coef
) on
(p-terms p
) by
#'cddr
460 do
(setq ans
(ptptplus ans
461 (everysubst2 (poly-in-var coef v
)
462 (list (p-var p
) exp
1))))
463 finally
(return ans
)))))
466 (or (null x
) (and (pcoefp (pt-lc x
)) (univar (pt-red x
)))))
468 ;;**THE CHINESE REMAINDER ALGORITHM IS A SPECIAL CASE OF LAGRANGE INTERPOLATION
470 (defun lagrange3 (u uk p qk
)
472 (setq uk
(pdifference uk
(pmod u
)))
473 (cond ((pzerop uk
) (setq modulus nil
) u
)
474 (t (setq uk
(pctimes (crecip (cmod qk
)) uk
))
476 (pplus u
(pctimes qk uk
)))))
479 (defun lagrange33 (u uk qk xk
)
480 (declare (special xv
))
481 (setq uk
(pdifference uk
(pcsubst u xk xv
)))
482 (cond ((pzerop uk
) u
)
484 (pctimes (crecip (pcsubst qk xk xv
)) uk
)
488 ;;;*************************************************************
490 ;; THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 3.
491 ;; IT INCLUDES THE GCD ROUTINES AND THEIR SUPPORTING FUNCTIONS
