Rename *ll* and *ul* to ll and ul in poles-in-interval
[maxima.git] / src / nrat4.lisp
blobe5383cd89352db649dc42b2bcbbbf3c6bc56d585
1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
4 ;;; ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
13 (macsyma-module nrat4)
15 (declare-top (special *exp *exp2 *radsubst *loglist $radsubstflag
16 *v *var))
18 (defmvar $radsubstflag nil
19 "`radsubstflag' `t' makes `ratsubs' call `radcan' when it appears useful")
22 (defun pdis (x) ($ratdisrep (pdis* x)))
24 (defun pdis* (x) `((mrat simp ,varlist ,genvar) ,x . 1))
26 (defun rdis (x) ($ratdisrep (rdis* x)))
28 (defun rdis* (x) `((mrat simp ,varlist ,genvar) . ,x))
30 (defun rform (x) (cdr (ratf x)))
32 (defmfun $ratcoef (e x &optional (n 1))
33 (ratcoeff e x n)) ; The spelling "ratcoeff" is nicer.
35 (defun ratcoeff (a b c)
36 (let* ((formflag ($ratp a))
37 (taylorform (and formflag (member 'trunc (cdar a) :test #'eq))))
38 (cond ((zerop1 b) (improper-arg-err b '$ratcoeff))
39 ((mbagp a) (cons (car a)
40 (mapcar #'(lambda (a) (ratcoeff a b c))
41 (cdr a))))
42 ((and taylorform (mnump c) (assolike b (cadddr (cdar a))))
43 (pscoeff1 a b c))
44 ((and taylorform (mexptp b) (mnump c) (mnump (caddr b))
45 (assolike (cadr b) (cadddr (cdar a))))
46 (pscoeff1 a (cadr b) (mul2 c (caddr b))))
47 ((and taylorform (equal c 0)) a)
48 (t (if taylorform (setq a (ratdisrep a)))
49 (setq a (let ($ratwtlvl)
50 (if (equal c 0)
51 (ratcoef (mul2* a b) b)
52 (ratcoef a (if (equal c 1) b (list '(mexpt) b c))))))
53 (if (and formflag (not taylorform))
54 (minimize-varlist a)
55 (ratdisrep a))))))
57 (defun minimize-varlist (ratfun)
58 (if (not ($ratp ratfun)) (setq ratfun (ratf ratfun)))
59 (minvarlist-mrat (caddr (car ratfun)) (cadddr (car ratfun))
60 (cdr ratfun)))
62 (defun minvarlist-mrat (vars gens ratform)
63 (let ((newgens (union* (listovars (car ratform))
64 (listovars (cdr ratform)))))
65 (do ((lv vars (cdr lv))
66 (lg gens (cdr lg))
67 (nlv ())
68 (nlg ()))
69 ((null lg)
70 (cons (list 'mrat 'simp (nreverse nlv) (nreverse nlg))
71 ratform))
72 (cond ((member (car lg) newgens :test #'eq)
73 (push (car lg) nlg)
74 (push (car lv) nlv))))))
76 (defun ratcoef (exp var)
77 (prog (varlist genvar $ratfac $algebraic $ratwtlvl bas minvar)
78 (setq var (ratdisrep var))
79 (setq bas (if (and (mexptp var) (mnump (caddr var))) (cadr var) var))
80 (newvar var)
81 (newvar bas)
82 (setq minvar (car varlist))
83 (newvar exp)
84 (setq exp (cdr (ratrep* exp)))
85 (setq var (cdr (ratrep* var)))
86 (setq bas (cadr (ratrep* bas)))
87 (if (and (onep1 (cdr exp)) (onep1 (cdr var)) (pureprod (car var)))
88 (return (pdis* (prodcoef (car var) (car exp)))))
89 (setq exp (ratquotient exp var))
90 (if (null minvar) (return (pdis* (prodcoef (cdr exp) (car exp)))))
91 (setq minvar (caadr (ratrep* minvar)))
92 loop (if (or (pcoefp (cdr exp)) (pointergp minvar (cadr exp)))
93 (return (rdis* (cdr (ratdivide exp bas)))))
94 (setq exp (ratcoef1 (car exp) (cdr exp)))
95 (go loop)))
97 (defun ratcoef1 (num den)
98 (cond ((pcoefp num) (rzero))
99 ((eq (car num) (car den)) (car (pdivide num den)))
100 ((pointergp (car den) (car num)) (rzero))
101 (t (ratcoef1 (constcoef (cdr num)) den))))
103 (defun constcoef (p)
104 (cond ((null p) 0)
105 ((zerop (car p)) (cadr p))
106 (t (constcoef (cddr p)))))
108 (setq *radsubst nil)
110 (defmfun $ratsubst (a b c) ; NEEDS CODE FOR FAC. FORM
111 (prog (varlist newvarlist dontdisrepit $ratfac genvar $keepfloat $float $numer)
112 ;; hard to maintain user ordering info.
113 (if ($ratp c) (setq dontdisrepit t))
114 (if (and $radsubstflag
115 (prog2 (newvar b) (some #'mexptp varlist)))
116 (let (($factorflag t) *exp *exp2 *radsubst)
117 (setq a (fullratsimp a))
118 (setq b (fullratsimp b))
119 (setq c (fullratsimp c))
120 (setq varlist nil)
121 (fnewvar b)
122 (fnewvar c)
123 (setq *exp (cdr (ratrep* b)))
124 (setq *exp2 (cdr (ratrep* c)))
125 ;; since *radsubst is t, both *exp and *exp2 will be radcan simplified
126 (setq *radsubst t)
127 (spc0 *var)
128 (setq b (rdis *exp) c (rdis *exp2))
129 (setq varlist nil))
130 (progn
131 (setq a ($rat a))
132 (setq b ($rat b))
133 (setq c ($rat c))
134 (setq varlist nil)))
135 (setq a ($ratdisrep a) b ($ratdisrep b) c ($ratdisrep c))
136 (cond ((integerp b) (setq c (ratf (maxima-substitute a b c)))
137 (return (cond (dontdisrepit c) (t ($ratdisrep c))))))
138 (newvar c)
139 (setq
140 newvarlist
141 (mapcar
142 #'(lambda (z)
143 (cond ((atom z) z)
144 (t (resimplify
145 (cons (car z)
146 (mapcar #'(lambda (zz)
147 (cond ((alike1 zz b) a)
148 ((atom zz) zz)
149 (t ($ratdisrep
150 ($ratsubst a b zz)))))
151 (cdr z)))))))
152 varlist))
153 (newvar a) (newvar b)
154 (setq newvarlist (reverse (pairoff (reverse varlist)
155 (reverse newvarlist))))
156 (setq a (cdr (ratrep* a)))
157 (setq b (cdr (ratrep* b)))
158 (setq c (cdr (ratrep* c)))
159 (when (pminusp (car b))
160 (setq b (ratminus b))
161 (setq a (ratminus a)))
162 (when (and (equal 1 (car b))
163 (not (equal 1 (cdr b)))
164 (not (equal 0 (car a))))
165 (setq a (ratinvert a))
166 (setq b (ratinvert b)))
167 (cond ((not (equal 1 (cdr b)))
168 (setq a (rattimes a (cons (cdr b) 1) t))
169 (setq b (cons (car b) 1))))
170 (setq c
171 (cond ((member (car b) '(0 1) :test #'equal)
172 (ratf (maxima-substitute (rdis a) b (rdis c))))
173 (t (cons (list 'mrat 'simp varlist genvar)
174 (if (equal (cdr a) 1)
175 (ratreduce (everysubst0 (car a) (car b) (car c))
176 (everysubst0 (car a) (car b) (cdr c)))
177 (allsubst00 a b c))))))
178 (unless (alike newvarlist varlist)
179 (setq varlist newvarlist
180 c (rdis (cdr c))
181 varlist nil
182 c (ratf c)))
183 (return (cond (dontdisrepit c) (t ($ratdisrep c))))))
185 (defun xptimes (x y) (if $ratwtlvl (wtptimes x y 0) (ptimes x y)))
187 (defun allsubst00 (a b c)
188 (cond ((equal a b) c)
189 ((not (equal (cdr b) 1)) c)
190 (t (ratquotient (everysubst00 a (car b) (car c))
191 (everysubst00 a (car b) (cdr c))))))
193 (defun everysubst00 (x i z)
194 (loop with ans = (rzero)
195 for (exp coef) on (everysubst i z *alpha) by #'cddr
196 do (setq ans (ratplus ans (rattimes (cons coef 1) (ratexpt x exp) t)))
197 finally (return ans)))
199 (defun everysubst0 (x i z)
200 (loop with ans = (pzero)
201 for (exp coef) on (everysubst i z *alpha) by #'cddr
202 do (setq ans (pplus ans (xptimes coef (pexpt x exp))))
203 finally (return ans)))
205 (defun everysubst1 (a b maxpow)
206 (loop for (exp coef) on (p-terms b) by #'cddr
207 for part = (everysubst a coef maxpow)
208 nconc (if (= 0 exp) part
209 (everysubst2 part (make-poly (p-var b) exp 1)))))
211 (defun everysubst2 (l h)
212 (do ((ptr l (cddr ptr)))
213 ((null ptr) l)
214 (setf (cadr ptr) (ptimes h (cadr ptr)))))
217 (defun pairoff (l m)
218 (cond ((null m) l) (t (cons (car m) (pairoff (cdr l) (cdr m))))))
220 ;;(DEFUN PAIROFF (L M)
221 ;; ;(COND ((NULL M) L) (T (CONS (CAR M) (PAIROFF (CDR L) (CDR M)))))
222 ;; (let ((ans nil))
223 ;; (dolist (x m (nreconc ans l))
224 ;; (push x ans) (setq l (cdr l)))))
226 (defun everysubst (a b maxpow)
227 (cond ((pcoefp a)
228 (cond ((equal a 1) (list maxpow b))
229 ((pcoefp b)
230 (list (setq maxpow
231 (do ((b b (quotient b a))
232 (ans 0 (1+ ans)))
233 ((or (> (abs a) (abs b))
234 (equal maxpow ans))
235 ans)))
236 (quotient b (setq maxpow (expt a maxpow)))
238 (rem b maxpow)))
239 (t (everysubst1 a b maxpow))))
240 ((or (pcoefp b) (pointergp (car a) (car b))) (list 0 b))
241 ((eq (car a) (car b))
242 (cond ((null (cdddr a)) (everypterms b (caddr a) (cadr a) maxpow))
243 (t (substforsum a b maxpow))))
244 (t (everysubst1 a b maxpow))))
246 (defun everypterms (x p n maxpow)
247 (if (< (cadr x) n)
248 (list 0 x)
249 (prog (k ans q part)
250 (setq k (car x))
251 (setq x (cdr x))
252 l (setq q (min maxpow (quotient (car x) n)))
253 m (when (equal q 0)
254 (return (if (null x)
256 (cons 0 (cons (psimp k x) ans)))))
257 (setq part (everysubst p (cadr x) q))
258 (setq ans (nconc (everypterms1 part k n (car x)) ans))
259 (setq x (cddr x))
260 (when (null x)
261 (setq q 0)
262 (go m))
263 (go l))))
265 (defun everypterms1 (l k n j)
266 (do ((ptr l (cddr ptr)))
267 ((null ptr) l)
268 (setf (cadr ptr)
269 (ptimes (psimp k (list (- j (* n (car ptr))) 1))
270 (cadr ptr)))))
272 (defun substforsum (a b maxpow)
273 (do ((pow 0 (1+ pow))
274 (quot) (zl-rem) (ans))
275 ((not (< pow maxpow)) (list* maxpow b ans))
276 (desetq (quot zl-rem) (pdivide b a))
277 (unless (and (equal (cdr quot) 1)
278 (not (pzerop (car quot)))
279 (equal (cdr zl-rem) 1))
280 (return (cons pow (cons b ans))))
281 (unless (pzerop (car zl-rem))
282 (setq ans (cons pow (cons (car zl-rem) ans))))
283 (setq b (car quot))))
285 (defun prodcoef (a b)
286 (cond ((pcoefp a)
287 (cond ((pcoefp b) (quotient b a)) (t (prodcoef1 a b))))
288 ((pcoefp b) (pzero))
289 ((pointergp (car a) (car b)) (pzero))
290 ((eq (car a) (car b))
291 (cond ((null (cdddr a))
292 (prodcoef (caddr a) (ptterm (cdr b) (cadr a))))
293 (t (sumcoef a b))))
294 (t (prodcoef1 a b))))
296 (defun sumcoef (a b)
297 (desetq (a b) (pdivide b a))
298 (if (and (equal (cdr a) 1) (equal (cdr b) 1))
299 (car a)
300 (pzero)))
302 (defun prodcoef1 (a b)
303 (loop with ans = (pzero)
304 for (bexp bcoef) on (p-terms b) by #'cddr
305 for part = (prodcoef a bcoef)
306 unless (pzerop part)
307 do (setq ans (pplus ans (psimp (p-var b) (list bexp part))))
308 finally (return ans)))
310 (defun pureprod (x)
311 (or (atom x)
312 (and (not (atom (cdr x)))
313 (null (cdddr x))
314 (pureprod (caddr x)))))
316 (defmfun $bothcoef (r var)
317 (prog (*var h varlist genvar $ratfac)
318 (unless ($ratp r)
319 (return `((mlist)
320 ,(setq h (coeff r var 1.))
321 ((mplus) ,r ((mtimes) -1 ,h ,var)))))
322 (newvar var)
323 (setq h (and varlist (car varlist)))
324 (newvar r)
325 (setq var (cdr (ratrep* var)))
326 (setq r (cdr (ratrep* r)))
327 (and h (setq h (caadr (ratrep* h))))
328 (cond ((and h (or (pcoefp (cdr r)) (pointergp h (cadr r)))
329 (equal 1 (cdr var)))
330 (setq var (bothprodcoef (car var) (car r)))
331 (return (list '(mlist)
332 (rdis* (ratreduce (car var) (cdr r)))
333 (rdis* (ratreduce (cdr var) (cdr r))))))
335 ;; CAN'T TELL WHAT BROUGHT US TO THIS POINT, SORRY
336 (merror (intl:gettext "bothcoef: invalid arguments."))))))
338 ;;COEFF OF A IN B
340 (defun bothprodcoef (a b)
341 (let ((c (prodcoef a b)))
342 (if (pzerop c) (cons (pzero) b) (cons c (pdifference b (ptimes c a))))))
344 (defvar argsfreeofp nil)
346 (defun argsfreeof (var e)
347 (let ((argsfreeofp t)) (freeof var e)))
349 ;;; This is a version of freeof for a list first argument
350 (defmfun $lfreeof (l e) "`freeof' for a list first argument"
351 (unless ($listp l)
352 (merror (intl:gettext "lfreeof: first argument must be a list; found: ~M") l))
353 (let ((exp ($totaldisrep e)))
354 (dolist (var (margs l) t)
355 (unless (freeof ($totaldisrep var) exp) (return nil)))))
357 (defmfun $freeof (&rest args)
358 (prog (l e)
359 (setq l (mapcar #'$totaldisrep (nreverse args))
360 e (car l))
361 loop (or (setq l (cdr l)) (return t))
362 (if (freeof (getopr (car l)) e) (go loop))
363 (return nil)))
365 (defun freeof (var e)
366 (cond ((alike1 var e) nil)
367 ((atom e) t)
368 ((and (not argsfreeofp)
369 (or (alike1 var ($verbify (caar e)))
370 (alike1 var ($nounify (caar e)))))
371 nil)
372 ((and (or (member (caar e) '(%product %sum %laplace) :test #'eq)
373 (and (eq (caar e) '%integrate) (cdddr e))
374 (and (eq (caar e) '%limit) (cddr e)))
375 (alike1 var (caddr e)))
376 (freeofl var (cdddr e)))
377 ((eq (caar e) '%at)
378 (cond ((not (freeofl var (hand-side (caddr e) 'r))) nil)
379 ((not (freeofl var (hand-side (caddr e) 'l))) t)
380 (t (freeof var (cadr e)))))
381 ((and (eq (caar e) 'lambda)
382 (not (member 'array (cdar e) :test #'eq))
383 ($listp (cadr e))
384 ; Check if var appears in the lambda list in any of the
385 ; following ways: var, 'var, [var] or ['var].
386 (some (lambda (v)
387 (or (eq v var)
388 (alike1 v `((mquote) ,var))
389 (alike1 v `((mlist) ,var))
390 (alike1 v `((mlist) ((mquote) ,var)))))
391 (cdadr e)))
393 ;; Check for a local variable in a block.
394 ((and (eq (caar e) 'mprog)
395 ($listp (cadr e))
396 ; Check if var appears in the variable list alone or
397 ; in an assignment
398 (some (lambda (v)
399 (or (eq v var)
400 (and (msetqp v)
401 (eq (cadr v) var))))
402 (cdadr e)))
404 ;; Check for a loop variable.
405 ((and (member (caar e) '(mdo mdoin) :test #'eq)
406 (alike1 var (cadr e)))
408 (argsfreeofp (freeofl var (margs e)))
409 (t (freeofl var (cdr e)))))
411 (defun freeofl (var l) (loop for x in l always (freeof var x)))
413 (defun hand-side (e flag)
414 (setq e (if (eq (caar e) 'mequal) (ncons e) (cdr e)))
415 (mapcar #'(lambda (u) (if (eq flag 'l) (cadr u) (caddr u))) e))
417 ;; subtitle radcan
419 (defmfun ($radcan :properties ((evfun t))) (exp)
420 (cond ((mbagp exp) (cons (car exp) (mapcar '$radcan (cdr exp))))
421 (t (let (($ratsimpexpons t))
422 (simplify (let (($expop 0) ($expon 0))
423 (radcan1 (fr1 exp nil) *var)))))))
425 (defun radcan1 (*exp *var)
426 (cond ((atom *exp) *exp)
427 (t (let (($factorflag t) varlist genvar $ratfac $norepeat
428 ($gcd (or $gcd (car *gcdl*)))
429 (radcanp t))
430 (newvar *exp)
431 (setq *exp (cdr (ratrep* *exp)))
432 (setq varlist
433 (mapcar
434 #'(lambda (x) (cond
435 ((atom x) x)
436 (t (cons (car x)
437 (mapcar #'(lambda (e)
438 (radcan1 e *var))
439 (cdr x))))))
440 varlist))
441 (spc0 *var)
442 (fr1 (rdis *exp) nil)))))
444 (defun spc0 (*var)
445 (prog (*v *loglist)
446 (if (allatoms varlist) (return nil))
447 (setq varlist (mapcar #'spc1 varlist)) ;make list of logs
448 (setq *loglist (factorlogs *loglist))
449 (mapc #'spc2 *loglist) ;subst log factorizations
450 (mapc #'spc3 varlist genvar) ;expand exponents
451 (mapc #'spc4 varlist) ;make exponent list
452 (desetq (varlist . genvar) (spc5 *v varlist genvar))
453 ;find expon dependencies
454 (setq varlist (mapcar #'(lambda (x)
455 (rjfsimp x *var))
456 varlist)) ;restore radicals
457 (mapc #'spc7 varlist))) ;simplify radicals
459 (defun allatoms (l)
460 (loop for x in l always (atom x)))
462 (defun rjfsimp (x *var &aux expon)
463 (cond ((and *radsubst $radsubstflag) x)
464 ((not (m$exp? (setq x (let ($logsimp) (resimplify x))))) x)
465 ((mlogp (setq expon (caddr x))) (cadr expon))
466 ((not (and (mtimesp expon) (or $logsimp *var))) x)
467 (t (do ((rischflag (and *var (not $logsimp) (not (freeof *var x))))
468 (power (cdr expon) (cdr power))) ;POWER IS A PRODUCT
469 ((null power) x)
470 (cond ((numberp (car power)))
471 ((mlogp (car power))
472 (and rischflag (cdr power) (return x))
473 (return
474 `((mexpt) ,(cadar power)
475 ,(muln (remove (car power) (cdr expon) :count 1 :test #'equal)
476 nil))))
477 (rischflag (return x)))))))
479 (defun dsubsta (x y zl)
480 (cond ((null zl) zl)
481 (t (cond ((alike1 y (car zl)) (rplaca zl x))
482 ((not (atom (car zl))) (dsubsta x y (cdar zl))))
483 (dsubsta x y (cdr zl))
484 zl)))
486 (defun radsubst (a b)
487 (setq *exp (allsubst00 a b *exp))
488 (if *radsubst (setq *exp2 (allsubst00 a b *exp2))))
490 (setq *var nil)
492 (defun spc1 (x)
493 (cond ((mlogp x) (putonloglist x))
494 ((and (mexptp x) (not (eq (cadr x) '$%e)))
495 ($exp-form (list '(mtimes)
496 (caddr x)
497 (putonloglist (list '(%log simp ratsimp)
498 (cadr x))))))
499 (t x)))
501 (defun putonloglist (l)
502 (unless (memalike l *loglist) (push l *loglist))
505 (defun spc2 (p)
506 (radsubst (rform (cdr p)) (rform (car p)))
507 (dsubsta (cdr p) (car p) varlist))
509 (defun spc2a (x) ;CONVERTS FACTORED
510 (let ((sum (mapcar #'spc2b x))) ;RFORM LOGAND TO SUM
511 (if (cdr sum) ;OF LOGS
512 (cons '(mplus) sum)
513 (car sum))))
515 (defun spc2b (x)
516 (let ((log `((%log simp ratsimp irreducible) ,(pdis (car x)))))
517 (if (equal 1 (cdr x)) log
518 (list '(mtimes) (cdr x) log))))
520 (defun spc3 (x v &aux y)
521 (when (and (m$exp? x)
522 (not (atom (setq y (caddr x))))
523 (mplusp (setq y (expand1 (if *var ($partfrac y *var) y) 10 10))))
524 (setq y (cons '(mtimes)
525 (mapcar #'(lambda (z) ($ratsimp ($exp-form z))) (cdr y))))
526 (radsubst (rform y) (rget v))
527 (dsubsta y x varlist)))
529 (defun spc4 (x)
530 (if (and (m$exp? x)
531 (not (memalike (caddr x) *v)))
532 (push (caddr x) *v)))
534 (defun rzcontent (r)
535 (destructuring-let (((c1 p) (pcontent (car r)))
536 ((c2 q) (pcontent (cdr r))))
537 (if (pminusp p) (setq p (pminus p) c1 (cminus c1)))
538 (cons (cons c1 c2) (cons p q))))
540 ;;The GCDLIST looks like (( GCM1pair occurrencepair11 occurrencepair12 ...) ...
541 ;;(GCMnpair occurrencepairn1 occurrencepairn2 ...))
542 ;;where GCMpairs are lists of ratforms and prefix forms for the greatest common
543 ;;multiple of the occurrencepairs. Each of these pairs is a list of a ratform
544 ;;and a prefix form. The prefix form is a pointer into the varlist.
545 ;;The occurrences are exponents of the base %E.
547 (defun spc5 (vl oldvarlist oldgenvar &aux gcdlist varlist genvar)
548 (dolist (v vl)
549 (destructuring-let* ((((c1 . c) . r) (rzcontent (rform v)))
550 (g (assoc r gcdlist :test #'equal)))
551 (cond (g (setf (cadr g) (plcm c (cadr g)))
552 (push (list ($exp-form (div* v c1)) c) (cddr g)))
553 (t (push (list r c (list ($exp-form (div* v c1)) c)) gcdlist)))))
554 (dolist (g gcdlist)
555 (let ((rd (rdis (car g))))
556 (when (and (mlogp rd) (memalike (cadr rd) oldvarlist))
557 (push (list (cadr rd) 1) (cddr g)))
558 (rplaca g ($exp-form (div rd (cadr g))))))
559 (spc5b gcdlist oldvarlist oldgenvar))
561 ;;(DEFUN SPC5B (V VARLIST GENVAR)
562 ;; (DOLIST (L V)
563 ;; (DOLIST (X (CDDR L))
564 ;; (UNLESS (EQUAL (CADR L) (CADR X))
565 ;; (RADSUBST (RATEXPT (RFORM (CAR L))
566 ;; (CAR (QUOTIENT (CADR X) (CADR L))))
567 ;; (RFORM (CAR X))))))
568 ;; (CONS VARLIST GENVAR))
571 (defun spc5b (v varlist genvar)
572 (dolist (l v)
573 (dolist (x (cddr l))
574 (unless (equal (cadr l) (cadr x))
575 (radsubst (ratexpt (rform (car l))
576 (quotient (cadr l) (cadr x)))
577 (rform (car x))))))
578 (cons varlist genvar))
580 (defun spc7 (x)
581 (if (eq x '$%i) (setq x '((mexpt) -1 ((rat) 1 2))))
582 (when (and (mexptp x)
583 (ratnump (caddr x)))
584 (let ((rad (rform x))
585 (rbase (rform (cadr x)))
586 (expon (caddr x)))
587 (radsubst (ratexpt rbase (cadr expon))
588 (ratexpt rad (caddr expon))))))
591 (defun goodform (l) ;;bad -> good
592 (loop for (exp coef) on l by #'cddr
593 collect (cons exp coef)))
595 (defun factorlogs (l)
596 (prog (negl posl maxpl maxnl maxn)
597 (dolist (log l)
598 (setq log
599 (cons log (goodform
600 (ratfact (rform (radcan1 (cadr log) *var))
601 #'pfactor))))
602 (cond ((equal (caadr log) -1) (push log negl))
603 (t (push log posl))))
604 (setq negl (flsort negl) posl (flsort posl) l (append negl posl))
605 (setq negl (mapcar #'cdr negl)
606 posl (mapcar #'cdr posl))
607 a (setq negl (delete '((-1 . 1)) negl :test #'equal))
608 (or negl
609 (return (mapc #'(lambda (x) (rplacd x (spc2a (cdr x)))) l)))
610 (setq maxnl (flmaxl negl)
611 maxn (caaar maxnl))
612 b (setq maxpl (flmaxl posl))
613 (cond ((and maxpl (flgreat (caaar maxpl) maxn))
614 (setq posl (flred posl (caaar maxpl)))
615 (go b))
616 ((and maxpl
617 (not (equal (caaar maxpl) maxn)))
618 (setq maxpl nil)))
619 (cond ((and (flevenp maxpl) (not (flevenp maxnl)))
620 (mapc #'(lambda (fp) (rplaca (car fp) (pminus (caar fp)))
621 (cond ((oddp (cdar fp))
622 (setq fp (delete '(-1 . 1) fp :test #'equal))
623 (setq negl (delete fp negl :test #'equal))
624 (and (cdr fp) (push (cdr fp) posl)))))
625 maxnl)
626 (go a))
627 (t (setq posl (flred posl maxn)
628 negl (flred negl maxn))
629 (go a)))))
631 (defun flevenp (pl)
632 (loop for l in pl never (oddp (cdar l))))
634 (defun flred (pl p)
635 (mapl #'(lambda (x) (if (equal p (caaar x))
636 (rplaca x (cdar x))))
638 (delete nil pl :test #'equal))
640 (defun flmaxl (fpl) ;lists of fac. polys
641 (cond ((null fpl) nil)
642 (t (do ((maxl (list (car fpl))
643 (cond ((equal (caaar maxl) (caaar ll))
644 (cons (car ll) maxl))
645 ((flgreat (caaar maxl) (caaar ll)) maxl)
646 (t (list (car ll)))))
647 (ll (cdr fpl) (cdr ll)))
648 ((null ll) maxl)))))
650 (defun flsort (fpl)
651 (mapc #'(lambda (x) (rplacd x (sort (cdr x) #'flgreat :key #'car)))
652 fpl))
654 (defun nmt (p any)
655 (cond ((pcoefp p)
656 (if (or any (cminusp p)) 1 0))
657 (t (loop for lp on (p-terms p) by #'cddr
658 sum (nmt (cadr lp) any)))))
660 (defun nmterms (p)
661 (cond ((equal p -1) (cons 0 0))
662 (t (cons (nmt p nil) (nmt p t)))))
664 (defun flgreat (p q)
665 (let ((pn (nmterms p)) (qn (nmterms q)))
666 (cond ((> (car pn) (car qn)) t)
667 ((< (car pn) (car qn)) nil)
668 ((> (cdr pn) (cdr qn)) t)
669 ((< (cdr pn) (cdr qn)) nil)
670 (t (flgreat1 p q)))))
672 (defun flgreat1 (p q)
673 (cond ((numberp p)
674 (cond ((numberp q) (> p q))
675 (t nil)))
676 ((numberp q) t)
677 ((pointergp (car p) (car q)) t)
678 ((pointergp (car q) (car p)) nil)
679 ((> (cadr p) (cadr q)) t)
680 ((< (cadr p) (cadr q)) nil)
681 (t (flgreat1 (caddr p) (caddr q)))))