Fix bug #4307: partswitch affects op and operatorp
[maxima.git] / src / rat3e.lisp
blob358ab89731e091f4835063589adceb7a116de1a5
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 rat3e)
15 ;; This is the rational function package part 5.
16 ;; It includes the conversion and top-level routines used
17 ;; by the rest of the functions.
19 (declare-top (special intbs* alflag var alc
20 vlist radlist expsumsplit *ratsimp*))
22 (defvar *ratsimp* nil)
24 (defmvar factorresimp nil "If `t' resimplifies factor(x-y) to x-y")
26 (defun mrateval (x)
27 (let ((varlist (caddar x)))
28 (cond ((and evp $infeval) (meval ($ratdisrep x)))
29 ((or evp
30 (and $float $keepfloat)
31 (not (alike varlist (mapcar #'meval varlist))))
32 (ratf (meval ($ratdisrep x))))
33 (t x))))
35 (defprop mrat mrateval mfexpr*)
37 (defmfun $ratnumer (x)
38 (cond ((mbagp x)
39 (cons (car x) (mapcar '$ratnumer (cdr x))))
41 (setq x (taychk2rat x))
42 (cons (car x) (cons (cadr x) 1)))))
44 (defmfun $ratdenom (x)
45 (cond ((mbagp x)
46 (cons (car x) (mapcar '$ratdenom (cdr x))))
48 (setq x (taychk2rat x))
49 (cons (car x) (cons (cddr x) 1)))))
51 (defun taychk2rat (x)
52 (cond ((and ($ratp x)
53 (member 'trunc (cdar x) :test #'eq))
54 ($taytorat x))
55 (t ($rat x))))
57 (defun tellratdisp (x)
58 (pdisrep+ (trdisp1 (cdr x) (car x))))
60 (defun trdisp1 (p var)
61 (cond ((null p) nil)
62 (t (cons (pdisrep* (if (mtimesp (cadr p))
63 (copy-list (cadr p))
64 (cadr p)) ;prevents clobbering p
65 (pdisrep! (car p) var))
66 (trdisp1 (cddr p) var)))))
68 (defmfun $untellrat (&rest args)
69 (dolist (x args)
70 (if (setq x (assol x tellratlist))
71 (setq tellratlist (remove x tellratlist :test #'equal))))
72 (cons '(mlist) (mapcar #'tellratdisp tellratlist)))
74 (defmfun $tellrat (&rest args)
75 (mapc #'tellrat1 args)
76 (unless (null args) (add2lnc 'tellratlist $myoptions))
77 (cons '(mlist) (mapcar #'tellratdisp tellratlist)))
79 (defun tellrat1 (x &aux varlist genvar $algebraic $ratfac algvar)
80 (setq x ($totaldisrep x))
81 (and (not (atom x)) (eq (caar x) 'mequal) (newvar (cadr x)))
82 (newvar (setq x (meqhk x)))
83 (unless varlist (merror (intl:gettext "tellrat: argument must be a polynomial; found: ~M") x))
84 (setq algvar (car (last varlist)))
85 (setq x (p-terms (primpart (cadr (ratrep* x)))))
86 (unless (equal (pt-lc x) 1) (merror (intl:gettext "tellrat: minimal polynomial must be monic.")))
87 (do ((p (pt-red x) (pt-red p)))
88 ((ptzerop p))
89 (setf (pt-lc p) (pdis (pt-lc p))))
90 (setq algvar (cons algvar x))
91 (if (setq x (assol (car algvar) tellratlist))
92 (setq tellratlist (remove x tellratlist :test #'equal)))
93 (push algvar tellratlist))
96 (defmfun $printvarlist ()
97 (cons '(mlist) (copy-tree varlist)))
99 (defmfun $showratvars (e)
100 (cons '(mlist simp)
101 (cond (($ratp e)
102 (if (member 'trunc (cdar e) :test #'eq) (setq e ($taytorat e)))
103 (caddar (minimize-varlist e)))
104 (t (let (varlist) (lnewvar e) varlist)))))
106 (defmfun $ratvars (&rest args)
107 (add2lnc '$ratvars $myoptions)
108 (setq $ratvars (cons '(mlist simp) (setq varlist (mapfr1 args varlist)))))
110 (defun mapfr1 (l varlist)
111 (mapcar #'(lambda (z) (fr1 z varlist)) l))
113 (defmvar inratsimp nil)
115 (defmfun ($fullratsimp :properties ((evfun t))) (exp &rest argl)
116 (prog (exp1)
117 loop (setq exp1 (simplify (apply #'$ratsimp (cons exp argl))))
118 (when (alike1 exp exp1) (return exp))
119 (setq exp exp1)
120 (go loop)))
122 (defun fullratsimp (l)
123 (let (($expop 0) ($expon 0) (inratsimp t) $ratsimpexpons)
124 (when (not ($ratp l))
125 ;; Not a mrat expression. Remove the special representation.
126 (setq l (specrepcheck l)))
127 (setq l ($totaldisrep l))
128 (fr1 l varlist)))
130 (defmfun $totaldisrep (l)
131 (cond ((atom l) l)
132 ((not (among 'mrat l)) l)
133 ((eq (caar l) 'mrat) (ratdisrep l))
134 (t (cons (delete 'ratsimp (car l) :test #'eq) (mapcar '$totaldisrep (cdr l))))))
136 ;;;VARLIST HAS MAIN VARIABLE AT END
138 (defun joinvarlist (cdrl)
139 (mapc #'(lambda (z) (unless (memalike z varlist) (push z varlist)))
140 (reverse (mapfr1 cdrl nil))))
142 (defmfun $rat (e &rest vars)
143 (cond ((not (null vars))
144 (let (varlist)
145 (joinvarlist vars)
146 (lnewvar e)
147 (rat0 e)))
149 (lnewvar e)
150 (rat0 e))))
152 (defun rat0 (exp) ;SIMP FLAGS?
153 (if (mbagp exp)
154 (cons (car exp) (mapcar #'rat0 (cdr exp)))
155 (ratf exp)))
157 (defmfun ($ratsimp :properties ((evfun t))) (e &rest vars)
158 (cond ((not (null vars))
159 (let (varlist)
160 (joinvarlist vars)
161 (fullratsimp e)))
162 (t (fullratsimp e))))
164 ;; $RATSIMP, $FULLRATSIMP and $RAT allow for optional extra
165 ;; arguments specifying the VARLIST.
167 ;;;PSQFR HAS NOT BEEN CHANGED TO MAKE USE OF THE SQFR FLAGS YET
169 (defmfun $sqfr (x)
170 (let ((varlist (cdr $ratvars)) genvar $keepfloat $ratfac)
171 (sublis '((factored . sqfred) (irreducible . sqfr))
172 (ffactor x #'psqfr))))
174 (declare-top (special fn))
176 (defun whichfn (p)
177 (cond ((and (mexptp p) (integerp (caddr p)))
178 (list '(mexpt) (whichfn (cadr p)) (caddr p)))
179 ((mtimesp p)
180 (cons '(mtimes) (mapcar #'whichfn (cdr p))))
181 (fn (ffactor p #'pfactor))
182 (t (factoralg p))))
184 (declare-top (special var))
186 (defun factoralg (p)
187 (prog (alc ans adn* $gcd)
188 (setq $gcd '$algebraic)
189 (when (or (atom p) (numberp p)) (return p))
190 (setq adn* 1)
191 (when (and (not $nalgfac) (not intbs*))
192 (setq intbs* (findibase minpoly*)))
193 (setq algfac* t)
194 (setq ans (ffactor p #'pfactor))
195 (cond ((eq (caar ans) 'mplus)
196 (return p))
197 (mplc*
198 (setq ans (albk ans))))
199 (if (and (not alc) (equal 1 adn*)) (return ans))
200 (setq ans (partition ans (car (last varlist)) 1))
201 (return (mul (let ((dosimp t))
202 (mul `((rat) 1 ,adn*)
203 (car ans)
204 (if alc (pdis alc) 1)))
205 (cdr ans)))))
207 (defun albk (p) ;to undo monicizing subst
208 (let ((*alpha* (pdis *alpha*)) ($ratfac t))
209 (declare (special *alpha*))
210 ;; don't multiply them back out
211 (maxima-substitute (list '(mtimes simp) mplc* *alpha*) ;assumes mplc* is int
212 *alpha* p)))
215 (defmfun $gfactor (p &aux (gauss t))
216 (when ($ratp p) (setq p ($ratdisrep p)))
217 (setq p ($factor (subst '%i '$%i p) '((mplus) 1 ((mexpt) %i 2))))
218 (setq p (sublis '((factored . gfactored) (irreducible . irreducibleg)) p))
219 (let (($expop 0) ($expon 0) $negdistrib)
220 (maxima-substitute '$%i '%i p)))
222 (defmfun ($factor :properties ((evfun t))) (e &optional (mp nil mp?))
223 (let ($intfaclim (varlist (cdr $ratvars)) genvar ans)
224 (setq ans (if mp? (factor e mp) (factor e)))
225 (if (and factorresimp $negdistrib
226 (mtimesp ans) (null (cdddr ans))
227 (equal (cadr ans) -1) (mplusp (caddr ans)))
228 (let (($expop 0) ($expon 0))
229 ($multthru ans))
230 ans)))
232 (defun factor (e &optional (mp nil mp?))
233 (let ((tellratlist nil)
234 (varlist varlist)
235 (genvar nil)
236 ($gcd $gcd)
237 ($negdistrib $negdistrib))
238 (prog (fn var mm* mplc* intbs* alflag minpoly* *alpha* p algfac*
239 $keepfloat $algebraic cargs)
240 (declare (special cargs fn *alpha*))
241 (unless (member $gcd *gcdl* :test #'eq)
242 (setq $gcd (car *gcdl*)))
243 (let ($ratfac)
244 (setq p e
245 mm* 1
246 cargs (if mp? (list mp) nil))
247 (when (symbolp p) (return p))
248 (when ($numberp p)
249 (return (let (($factorflag (not scanmapp)))
250 (factornumber p))))
251 (when (mbagp p)
252 (return (cons (car p) (mapcar #'(lambda (x) (apply #'factor (cons x cargs))) (cdr p)))))
253 (cond (mp?
254 (setq *alpha* (meqhk mp))
255 (newvar *alpha*)
256 (setq minpoly* (cadr (ratrep* *alpha*)))
257 (when (or (pcoefp minpoly*)
258 (not (univar (cdr minpoly*)))
259 (< (cadr minpoly*) 2))
260 (merror (intl:gettext "factor: second argument must be a nonlinear, univariate polynomial; found: ~M") *alpha*))
261 (setq *alpha* (pdis (list (car minpoly*) 1 1))
262 mm* (cadr minpoly*))
263 (unless (equal (caddr minpoly*) 1)
264 (setq mplc* (caddr minpoly*))
265 (setq minpoly* (pmonz minpoly*))
266 (setq p (maxima-substitute (div *alpha* mplc*) *alpha* p)))
267 (setq $algebraic t)
268 ($tellrat(pdis minpoly*))
269 (setq algfac* t))
271 (setq fn t)))
272 (unless scanmapp (setq p (let (($ratfac t)) (sratsimp p))))
273 (newvar p)
274 (when (symbolp p) (return p))
275 (when (numberp p)
276 (return (let (($factorflag (not scanmapp)))
277 (factornumber p))))
278 (setq $negdistrib nil)
279 (setq p (let ($factorflag ($ratexpand $facexpand))
280 (whichfn p))))
282 (setq p (let (($expop 0) ($expon 0))
283 (simplify p)))
284 (cond ((mnump p) (return (factornumber p)))
285 ((atom p) (return p)))
286 (and $ratfac (not $factorflag) ($ratp e) (return ($rat p)))
287 (and $factorflag (mtimesp p) (mnump (cadr p))
288 (setq *alpha* (factornumber (cadr p)))
289 (cond ((alike1 *alpha* (cadr p)))
290 ((mtimesp *alpha*)
291 (setq p (cons (car p) (append (cdr *alpha*) (cddr p)))))
293 (setq p (cons (car p) (cons *alpha* (cddr p)))))))
294 (when (null (member 'factored (car p) :test #'eq))
295 (setq p (cons (append (car p) '(factored)) (cdr p))))
296 (return p))))
298 (defun factornumber (n)
299 (setq n (nretfactor1 (nratfact (cdr ($rat n)))))
300 (cond ((cdr n)
301 (cons '(mtimes simp factored)
302 (if (equal (car n) -1)
303 (cons (car n) (nreverse (cdr n)))
304 (nreverse n))))
305 ((atom (car n))
306 (car n))
308 (cons (cons (caaar n) '(simp factored)) (cdar n)))))
310 (defun nratfact (x)
311 (cond ((equal (cdr x) 1) (cfactor (car x)))
312 ((equal (car x) 1) (revsign (cfactor (cdr x))))
313 (t (nconc (cfactor (car x)) (revsign (cfactor (cdr x)))))))
315 ;;; FOR LISTS OF JUST NUMBERS
316 (defun nretfactor1 (l)
317 (cond ((null l) nil)
318 ((equal (cadr l) 1) (cons (car l) (nretfactor1 (cddr l))))
319 (t (cons (if (equal (cadr l) -1)
320 (list '(rat simp) 1 (car l))
321 (list '(mexpt simp) (car l) (cadr l)))
322 (nretfactor1 (cddr l))))))
324 (declare-top (unspecial var))
326 (defmfun $polymod (p &optional (m 0 m?))
327 (let ((modulus modulus))
328 (when m?
329 (setq modulus m)
330 (when (or (not (integerp modulus)) (zerop modulus))
331 (merror (intl:gettext "polymod: modulus must be a nonzero integer; found: ~M") modulus)))
332 (when (minusp modulus)
333 (setq modulus (abs modulus)))
334 (mod1 p)))
336 (defun mod1 (e)
337 (if (mbagp e) (cons (car e) (mapcar 'mod1 (cdr e)))
338 (let (formflag)
339 (newvar e)
340 (setq formflag ($ratp e) e (ratrep* e))
341 (setq e (cons (car e) (ratreduce (pmod (cadr e)) (pmod (cddr e)))))
342 (cond (formflag e) (t (ratdisrep e))))))
344 (defmfun $divide (x y &rest vars)
345 (prog (h varlist tt ty formflag $ratfac)
346 (when (and ($ratp x) (setq formflag t) (integerp (cadr x)) (equal (cddr x) 1))
347 (setq x (cadr x)))
348 (when (and ($ratp y) (setq formflag t) (integerp (cadr y)) (equal (cddr y) 1))
349 (setq y (cadr y)))
350 (when (and (integerp x) (integerp y))
351 (when (zerop y)
352 (merror (intl:gettext "divide: Quotient by zero")))
353 (return (cons '(mlist) (multiple-value-list (truncate x y)))))
354 (setq varlist vars)
355 (mapc #'newvar (reverse (cdr $ratvars)))
356 (newvar y)
357 (newvar x)
358 (setq x (ratrep* x))
359 (setq h (car x))
360 (setq x (cdr x))
361 (setq y (cdr (ratrep* y)))
362 (cond ((and (equal (setq tt (cdr x)) 1) (equal (cdr y) 1))
363 (setq x (pdivide (car x) (car y))))
364 (t (setq ty (cdr y))
365 (setq x (ptimes (car x) (cdr y)))
366 (setq x (pdivide x (car y)))
367 (setq x (list
368 (ratqu (car x) tt)
369 (ratqu (cadr x) (ptimes tt ty))))))
370 (setq h (list '(mlist) (cons h (car x)) (cons h (cadr x))))
371 (return (if formflag h ($totaldisrep h)))))
373 (defmfun $quotient (&rest args)
374 (cadr (apply '$divide args)))
376 (defmfun $remainder (&rest args)
377 (caddr (apply '$divide args)))
379 (defmfun $gcd (x y &rest vars)
380 (prog (h varlist genvar $keepfloat formflag)
381 (setq formflag ($ratp x))
382 (and ($ratp y) (setq formflag t))
383 (setq varlist vars)
384 (dolist (v varlist)
385 (when (numberp v) (improper-arg-err v '$gcd)))
386 (newvar x)
387 (newvar y)
388 (when (and ($ratp x) ($ratp y) (equal (car x) (car y)))
389 (setq genvar (car (last (car x))) h (car x) x (cdr x) y (cdr y))
390 (go on))
391 (setq x (ratrep* x))
392 (setq h (car x))
393 (setq x (cdr x))
394 (setq y (cdr (ratrep* y)))
395 on (setq x (cons (pgcd (car x) (car y)) (plcm (cdr x) (cdr y))))
396 (setq h (cons h x))
397 (return (if formflag h (ratdisrep h)))))
399 (defmfun $content (x &rest vars)
400 (prog (y h varlist formflag)
401 (setq formflag ($ratp x))
402 (setq varlist vars)
403 (newvar x)
404 (desetq (h x . y) (ratrep* x))
405 (unless (atom x)
406 ;; (CAR X) => gensym corresponding to apparent main var.
407 ;; MAIN-GENVAR => gensym corresponding to the genuine main var.
408 (let ((main-genvar (nth (1- (length varlist)) genvar)))
409 (unless (eq (car x) main-genvar)
410 (setq x `(,main-genvar 0 ,x)))))
411 (setq x (rcontent x)
412 y (cons 1 y))
413 (setq h (list '(mlist)
414 (cons h (rattimes (car x) y nil))
415 (cons h (cadr x))))
416 (return (if formflag h ($totaldisrep h)))))
418 (defun pget (gen)
419 (cons gen '(1 1)))
421 (defun m$exp? (x)
422 (and (mexptp x) (eq (cadr x) '$%e)))
424 (defun algp ($x)
425 (algpchk $x nil))
427 (defun algpget ($x)
428 (algpchk $x t))
430 (defun algpchk ($x mpflag &aux temp)
431 (cond ((eq $x '$%i) '(2 -1))
432 ((eq $x '$%phi) '(2 1 1 -1 0 -1))
433 ((radfunp $x nil)
434 (if (not mpflag) t
435 (let ((r (prep1 (cadr $x))))
436 (cond ((onep1 (cdr r)) ;INTEGRAL ALG. QUANT.?
437 (list (caddr (caddr $x))
438 (car r)))
439 (*ratsimp* (setq radlist (cons $x radlist)) nil)))))
440 ((not $algebraic) nil)
441 ((and (m$exp? $x) (mtimesp (setq temp (caddr $x)))
442 (equal (cddr temp) '($%i $%pi))
443 (ratnump (setq temp (cadr temp))))
444 (if mpflag (primcyclo (* 2 (caddr temp))) t))
445 ((not mpflag) (assolike $x tellratlist))
446 ((setq temp (copy-list (assolike $x tellratlist)))
447 (do ((p temp (cddr p))) ((null p))
448 (rplaca (cdr p) (car (prep1 (cadr p)))))
449 (setq temp
450 (cond ((ptzerop (pt-red temp)) (list (pt-le temp) (pzero)))
451 ((zerop (pt-le (pt-red temp)))
452 (list (pt-le temp) (pminus (pt-lc (pt-red temp)))))
453 (t temp)))
454 (if (and (= (pt-le temp) 1) (setq $x (assol $x genpairs)))
455 (rplacd $x (cons (cadr temp) 1)))
456 temp)))
458 (defun radfunp (x funcflag) ;FUNCFLAG -> TEST FOR ALG FUNCTION NOT NUMBER
459 (cond ((atom x) nil)
460 ((not (eq (caar x) 'mexpt)) nil)
461 ((not (ratnump (caddr x))) nil)
462 (funcflag (not (numberp (cadr x))))
463 (t t)))
465 (defun ratsetup (vl gl)
466 (ratsetup1 vl gl) (ratsetup2 vl gl))
468 (defun ratsetup1 (vl gl)
469 (and $ratwtlvl
470 (mapc #'(lambda (v g)
471 (setq v (assolike v *ratweights))
472 (if v (putprop g v '$ratweight) (remprop g '$ratweight)))
473 vl gl)))
475 (defun ratsetup2 (vl gl)
476 (when $algebraic
477 (mapc #'(lambda (g) (remprop g 'algord)) gl)
478 (mapl #'(lambda (v lg)
479 (cond ((setq v (algpget (car v)))
480 (algordset v lg) (putprop (car lg) v 'tellrat))
481 (t (remprop (car lg) 'tellrat))))
482 vl gl))
483 (and $ratfac (let ($ratfac)
484 (mapc #'(lambda (v g)
485 (if (mplusp v)
486 (putprop g (car (prep1 v)) 'unhacked)
487 (remprop g 'unhacked)))
488 vl gl))))
490 (defun porder (p)
491 (if (pcoefp p) 0 (valget (car p))))
493 (defun algordset (x gl)
494 (do ((p x (cddr p))
495 (mv 0))
496 ((null p)
497 (do ((l gl (cdr l)))
498 ((or (null l) (> (valget (car l)) mv)))
499 (putprop (car l) t 'algord)))
500 (setq mv (max mv (porder (cadr p))))))
502 (defun gensym-readable (name)
503 (cond ((symbolp name)
504 (gensym (string-trim "$" (string name))))
506 (setq name (aformat nil "~:M" name))
507 (if name (gensym name) (gensym)))))
509 (defun orderpointer (l)
510 (loop for v in l
511 for i below (- (length l) (length genvar))
512 collecting (gensym-readable v) into tem
513 finally (setq genvar (nconc tem genvar))
514 (return (prenumber genvar 1))))
516 (defun creatsym (n)
517 (dotimes (i n)
518 (push (gensym) genvar)))
520 (defun prenumber (v n)
521 (do ((vl v (cdr vl))
522 (i n (1+ i)))
523 ((null vl) nil)
524 (setf (symbol-value (car vl)) i)))
526 (defun rget (genv)
527 (cons (if (and $ratwtlvl
528 (or (fixnump $ratwtlvl)
529 (merror (intl:gettext "rat: 'ratwtlvl' must be an integer; found: ~M") $ratwtlvl))
530 (> (or (get genv '$ratweight) -1) $ratwtlvl))
531 (pzero)
532 (pget genv))
535 (defun ratrep (x varl)
536 (setq varlist varl)
537 (ratrep* x))
539 (defun ratrep* (x)
540 (let (genpairs)
541 (orderpointer varlist)
542 (ratsetup1 varlist genvar)
543 (mapc #'(lambda (y z) (push (cons y (rget z)) genpairs)) varlist genvar)
544 (ratsetup2 varlist genvar)
545 (xcons (prep1 x) ; PREP1 changes VARLIST
546 (list* 'mrat 'simp varlist genvar ; when $RATFAC is T.
547 (if (and (not (atom x)) (member 'irreducible (cdar x) :test #'eq))
548 '(irreducible))))))
550 (defvar *withinratf* nil)
552 (defun ratf (l)
553 (prog (u *withinratf*)
554 (setq *withinratf* t)
555 (when (eq '%% (catch 'ratf (newvar l)))
556 (setq *withinratf* nil)
557 (return (srf l)))
558 (setq u (catch 'ratf (ratrep* l))) ; for truncation routines
559 (return (or u (prog2 (setq *withinratf* nil) (srf l))))))
562 (defun prep1 (x &aux temp)
563 (cond ((floatp x)
564 (cond ($keepfloat (cons x 1.0))
565 ((prepfloat x))))
566 ((integerp x) (cons (cmod x) 1))
567 ((rationalp x)
568 (if (null modulus)
569 (cons (numerator x) (denominator x))
570 (cquotient (numerator x) (denominator x))))
571 ((atom x) (cond ((assolike x genpairs))
572 (t (newsym x))))
573 ((and $ratfac (assolike x genpairs)))
574 ((eq (caar x) 'mplus)
575 (cond ($ratfac
576 (setq x (mapcar #'prep1 (cdr x)))
577 (cond ((every #'frpoly? x)
578 (cons (mfacpplus (mapl #'(lambda (x) (rplaca x (caar x))) x)) 1))
579 (t (do ((a (car x) (facrplus a (car l)))
580 (l (cdr x) (cdr l)))
581 ((null l) a)))))
582 (t (do ((a (prep1 (cadr x)) (ratplus a (prep1 (car l))))
583 (l (cddr x) (cdr l)))
584 ((null l) a)))))
585 ((eq (caar x) 'mtimes)
586 (do ((a (savefactors (prep1 (cadr x)))
587 (rattimes a (savefactors (prep1 (car l))) sw))
588 (l (cddr x) (cdr l))
589 (sw (not (and $norepeat (member 'ratsimp (cdar x) :test #'eq)))))
590 ((null l) a)))
591 ((eq (caar x) 'mexpt)
592 (newvarmexpt x (caddr x) t))
593 ((eq (caar x) 'mquotient)
594 (ratquotient (savefactors (prep1 (cadr x)))
595 (savefactors (prep1 (caddr x)))))
596 ((eq (caar x) 'mminus)
597 (ratminus (prep1 (cadr x))))
598 ((eq (caar x) 'rat)
599 (cond (modulus (cons (cquotient (cmod (cadr x)) (cmod (caddr x))) 1))
600 (t (cons (cadr x) (caddr x)))))
601 ((eq (caar x) 'bigfloat)(bigfloat2rat x))
602 ((eq (caar x) 'mrat)
603 (cond ((and *withinratf* (member 'trunc (car x) :test #'eq))
604 (throw 'ratf nil))
605 ((catch 'compatvl
606 (progn
607 (setq temp (compatvarl (caddar x) varlist (cadddr (car x)) genvar))
609 (cond ((member 'trunc (car x) :test #'eq)
610 (cdr ($taytorat x)))
611 ((and (not $keepfloat)
612 (or (pfloatp (cadr x)) (pfloatp (cddr x))))
613 (cdr (ratrep* ($ratdisrep x))))
614 ((sublis temp (cdr x)))))
615 (t (cdr (ratrep* ($ratdisrep x))))))
616 ((assolike x genpairs))
617 (t (setq x (littlefr1 x))
618 (cond ((assolike x genpairs))
619 (t (newsym x))))))
622 (defun putonvlist (x)
623 (push x vlist)
624 (and $algebraic
625 (setq x (assolike x tellratlist))
626 (mapc 'newvar1 x)))
628 (setq expsumsplit t) ;CONTROLS SPLITTING SUMS IN EXPONS
630 (defun newvarmexpt (x e flag)
631 ;; WHEN FLAG IS T, CALL RETURNS RATFORM
632 (prog (topexp)
633 (when (and (integerp e) (not flag))
634 (return (newvar1 (cadr x))))
635 (setq topexp 1)
636 top (cond
638 ;; X=B^N FOR N A NUMBER
639 ((integerp e)
640 (setq topexp (* topexp e))
641 (setq x (cadr x)))
642 ((atom e) nil)
644 ;; X=B^(P/Q) FOR P AND Q INTEGERS
645 ((eq (caar e) 'rat)
646 (cond ((or (minusp (cadr e)) (> (cadr e) 1))
647 (setq topexp (* topexp (cadr e)))
648 (setq x (list '(mexpt)
649 (cadr x)
650 (list '(rat) 1 (caddr e))))))
651 (cond ((or flag (numberp (cadr x)) ))
652 (*ratsimp*
653 (cond ((memalike x radlist) (return nil))
654 (t (setq radlist (cons x radlist))
655 (return (newvar1 (cadr x))))) )
656 ($algebraic (newvar1 (cadr x)))))
657 ;; X=B^(A*C)
658 ((eq (caar e) 'mtimes)
659 (cond
660 ((or
662 ;; X=B^(N *C)
663 (and (atom (cadr e))
664 (integerp (cadr e))
665 (setq topexp (* topexp (cadr e)))
666 (setq e (cddr e)))
668 ;; X=B^(P/Q *C)
669 (and (not (atom (cadr e)))
670 (eq (caaadr e) 'rat)
671 (not (equal 1 (cadadr e)))
672 (setq topexp (* topexp (cadadr e)))
673 (setq e (cons (list '(rat)
675 (caddr (cadr e)))
676 (cddr e)))))
677 (setq x
678 (list '(mexpt)
679 (cadr x)
680 (setq e (simplify (cons '(mtimes)
681 e)))))
682 (go top))))
684 ;; X=B^(A+C)
685 ((and (eq (caar e) 'mplus) expsumsplit) ;SWITCH CONTROLS
686 (setq ;SPLITTING EXPONENT
687 x ;SUMS
688 (cons
689 '(mtimes)
690 (mapcar
691 #'(lambda (ll)
692 (list '(mexpt)
693 (cadr x)
694 (simplify (list '(mtimes)
695 topexp
696 ll))))
697 (cdr e))))
698 (cond (flag (return (prep1 x)))
699 (t (return (newvar1 x))))))
700 (cond (flag nil)
701 ((equal 1 topexp)
702 (cond ((or (atom x)
703 (not (eq (caar x) 'mexpt)))
704 (newvar1 x))
705 ((or (memalike x varlist) (memalike x vlist))
706 nil)
707 (t (cond ((or (atom x) (null *fnewvarsw))
708 (putonvlist x))
709 (t (setq x (littlefr1 x))
710 (mapc #'newvar1 (cdr x))
711 (or (memalike x vlist)
712 (memalike x varlist)
713 (putonvlist x)))))))
714 (t (newvar1 x)))
715 (return
716 (cond
717 ((null flag) nil)
718 ((equal 1 topexp)
719 (cond
720 ((and (not (atom x)) (eq (caar x) 'mexpt))
721 (cond ((assolike x genpairs))
722 ;; *** SHOULD ONLY GET HERE IF CALLED FROM FR1. *FNEWVARSW=NIL
723 (t (setq x (littlefr1 x))
724 (cond ((assolike x genpairs))
725 (t (newsym x))))))
726 (t (prep1 x))))
727 (t (ratexpt (prep1 x) topexp))))))
729 (defun newvar1 (x)
730 (cond ((numberp x) nil)
731 ((memalike x varlist) nil)
732 ((memalike x vlist) nil)
733 ((atom x) (putonvlist x))
734 ((member (caar x)
735 '(mplus mtimes rat mdifference
736 mquotient mminus bigfloat) :test #'eq)
737 (mapc #'newvar1 (cdr x)))
738 ((eq (caar x) 'mexpt)
739 (newvarmexpt x (caddr x) nil))
740 ((eq (caar x) 'mrat)
741 (and *withinratf* (member 'trunc (cdddar x) :test #'eq) (throw 'ratf '%%))
742 (cond ($ratfac (mapc 'newvar3 (caddar x)))
743 (t (mapc #'newvar1 (reverse (caddar x))))))
744 (t (cond (*fnewvarsw (setq x (littlefr1 x))
745 (mapc #'newvar1 (cdr x))
746 (or (memalike x vlist)
747 (memalike x varlist)
748 (putonvlist x)))
749 (t (putonvlist x))))))
751 (defun newvar3 (x)
752 (or (memalike x vlist)
753 (memalike x varlist)
754 (putonvlist x)))
756 (defun fr1 (x varlist) ;put radicands on initial varlist?
757 (prog (genvar $norepeat *ratsimp* radlist vlist nvarlist ovarlist genpairs)
758 (newvar1 x)
759 (setq nvarlist (mapcar #'fr-args vlist))
760 (cond ((not *ratsimp*) ;*ratsimp* not set for initial varlist
761 (setq varlist (nconc (sortgreat vlist) varlist))
762 (return (rdis (cdr (ratrep* x))))))
763 (setq ovarlist (nconc vlist varlist)
764 vlist nil)
765 (mapc #'newvar1 nvarlist) ;*RATSIMP*=T PUTS RADICANDS ON VLIST
766 (setq nvarlist (nconc nvarlist varlist) ; RADICALS ON RADLIST
767 varlist (nconc (sortgreat vlist) (radsort radlist) varlist))
768 (orderpointer varlist)
769 (setq genpairs (mapcar #'(lambda (x y) (cons x (rget y))) varlist genvar))
770 (let (($algebraic $algebraic) ($ratalgdenom $ratalgdenom) radlist)
771 (and (not $algebraic)
772 (some #'algpget varlist) ;NEEDS *RATSIMP*=T
773 (setq $algebraic t $ratalgdenom nil))
774 (ratsetup varlist genvar)
775 (setq genpairs
776 (mapcar #'(lambda (x y) (cons x (prep1 y))) ovarlist nvarlist))
777 (setq x (rdis (prep1 x)))
778 (cond (radlist ;rational radicands
779 (setq *ratsimp* nil)
780 (setq x (ratsimp (simplify x) nil nil)))))
781 (return x)))
783 (defun ratsimp (x varlist genvar) ($ratdisrep (ratf x)))
785 (defun littlefr1 (x)
786 (cons (remove 'simp (car x) :test #'eq) (mapfr1 (cdr x) nil)))
788 ;;IF T RATSIMP FACTORS RADICANDS AND LOGANDS
789 (defmvar fr-factor nil)
791 (defun fr-args (x) ;SIMP (A/B)^N TO A^N/B^N ?
792 (cond ((atom x)
793 (when (eq x '$%i) (setq *ratsimp* t)) ;indicates algebraic present
795 (t (setq *ratsimp* t) ;FLAG TO CHANGED ELMT.
796 (simplify (zp (cons (remove 'simp (car x) :test #'eq)
797 (if (or (radfunp x nil) (eq (caar x) '%log))
798 (cons (if fr-factor (factor (cadr x))
799 (fr1 (cadr x) varlist))
800 (cddr x))
801 (let (modulus)
802 (mapfr1 (cdr x) varlist)))))))))
804 ;;SIMPLIFY MEXPT'S & RATEXPAND EXPONENT
806 (defun zp (x)
807 (if (and (mexptp x) (not (atom (caddr x))))
808 (list (car x) (cadr x)
809 (let ((varlist varlist) *ratsimp*)
810 ($ratexpand (caddr x))))
813 (defun newsym (e)
814 (prog (g p)
815 (when (setq g (assolike e genpairs))
816 (return g))
817 (setq g (gensym-readable e))
818 (putprop g e 'disrep)
819 (push e varlist)
820 (push (cons e (rget g)) genpairs)
821 (valput g (if genvar (1- (valget (car genvar))) 1))
822 (push g genvar)
823 (when (setq p (and $algebraic (algpget e)))
824 (algordset p genvar)
825 (putprop g p 'tellrat))
826 (return (rget g))))
828 ;; This control of conversion from float to rational appears to be explained
829 ;; nowhere. - RJF
831 (defun maxima-rationalize (x)
832 (cond ((not (floatp x)) x)
833 ;; Handle denormalized floats -- see maxima-discuss 2021-04-27 and bug 3777
834 ((and (not (= x 0.0)) (< (abs x) COMMON-LISP:LEAST-POSITIVE-NORMALIZED-DOUBLE-FLOAT))
835 (let ((r ($rationalize x)))
836 (cons (cadr r) (caddr r))))
837 ((< x 0.0)
838 (setq x (ration1 (* -1.0 x)))
839 (rplaca x (* -1 (car x))))
840 (t (ration1 x))))
842 (defun ration1 (x)
843 (let ((rateps (if (not (floatp $ratepsilon))
844 ($float $ratepsilon)
845 $ratepsilon)))
846 (or (and (zerop x) (cons 0 1))
847 (prog (y a)
848 (return
849 ;; I (rtoy) think this is computing a continued fraction
850 ;; expansion of the given float.
852 ;; FIXME? CMUCL used to use this routine for its
853 ;; RATIONALIZE function, but there were known bugs in
854 ;; that implementation, where the result was not very
855 ;; accurate. Unfortunately, I can't find the example
856 ;; that demonstrates this. In any case, CMUCL replaced
857 ;; it with an algorithm based on the code in Clisp, which
858 ;; was much better.
859 (do ((xx x (setq y (/ 1.0 (- xx (float a x)))))
860 (num (setq a (floor x)) (+ (* (setq a (floor y)) num) onum))
861 (den 1 (+ (* a den) oden))
862 (onum 1 num)
863 (oden 0 den))
864 ((or (zerop (- xx (float a x)))
865 (and (not (zerop den))
866 (not (> (abs (/ (- x (/ (float num x) (float den x))) x)) rateps))))
867 (cons num den))))))))
869 (defun prepfloat (f)
870 (cond (modulus (merror (intl:gettext "rat: can't rationalize ~M when modulus = ~M") f modulus))
871 ($ratprint (mtell (intl:gettext "~&rat: replaced ~A by") f)))
872 (setq f (maxima-rationalize f))
873 (when $ratprint
874 (mtell " ~A/~A = ~A~%" (car f) (cdr f) (fpcofrat1 (car f) (cdr f))))
877 (defun pdisrep (p)
878 (if (atom p)
880 (pdisrep+ (pdisrep2 (cdr p) (get (car p) 'disrep)))))
882 (defun pdisrep! (n var)
883 (cond ((zerop n) 1)
884 ((equal n 1) (cond ((atom var) var)
885 ((or (eq (caar var) 'mtimes)
886 (eq (caar var) 'mplus))
887 (copy-list var))
888 (t var)))
889 ((eql var 1) 1)
890 (t (list '(mexpt ratsimp) var n))))
892 (defun pdisrep+ (p)
893 (cond ((null (cdr p)) (car p))
894 (t (let ((a (last p)))
895 (cond ((mplusp (car a))
896 (rplacd a (cddar a))
897 (rplaca a (cadar a))))
898 (cons '(mplus ratsimp) p)))))
900 (defun pdisrep* (a b)
901 (cond ((equal a 1) b)
902 ((equal b 1) a)
903 (t (cons '(mtimes ratsimp) (nconc (pdisrep*chk a) (pdisrep*chk b))))))
905 (defun pdisrep*chk (a)
906 (if (mtimesp a) (cdr a) (ncons a)))
908 (defun pdisrep2 (p var)
909 (cond ((null p) nil)
910 ($ratexpand (pdisrep2expand p var))
911 (t (do ((l () (cons (pdisrep* (pdisrep (cadr p)) (pdisrep! (car p) var)) l))
912 (p p (cddr p)))
913 ((null p) (nreverse l))))))
915 (defmfun ($ratexpand :properties ((evfun t))) (x)
916 (if (mbagp x)
917 (cons (car x) (mapcar '$ratexpand (cdr x)))
918 (let (($ratexpand t) ($ratfac nil))
919 (ratdisrep (ratf x)))))
921 (defun pdisrep*expand (a b)
922 (cond ((equal a 1) (list b))
923 ((equal b 1) (list a))
924 ((or (atom a) (not (eq (caar a) 'mplus)))
925 (list (cons (quote (mtimes ratsimp))
926 (nconc (pdisrep*chk a) (pdisrep*chk b)))))
927 (t (mapcar #'(lambda (z) (if (equal z 1) b
928 (cons '(mtimes ratsimp)
929 (nconc (pdisrep*chk z)
930 (pdisrep*chk b)))))
931 (cdr a)))))
933 (defun pdisrep2expand (p var)
934 (cond ((null p) nil)
935 (t (nconc (pdisrep*expand (pdisrep (cadr p)) (pdisrep! (car p) var))
936 (pdisrep2expand (cddr p) var)))))
939 (defmvar $ratdenomdivide t)
941 (defmfun $ratdisrep (x)
942 (cond ((mbagp x)
943 ;; Distribute over lists, equations, and matrices.
944 (cons (car x) (mapcar #'$ratdisrep (cdr x))))
945 ((not ($ratp x)) x)
947 (setq x (ratdisrepd x))
948 (if (and (not (atom x))
949 (member 'trunc (cdar x) :test #'eq))
950 (cons (delete 'trunc (copy-list (car x)) :count 1 :test #'eq)
951 (cdr x))
952 x))))
954 ;; RATDISREPD is needed by DISPLA. - JPG
955 (defun ratdisrepd (x)
956 (mapc #'(lambda (y z) (putprop y z (quote disrep)))
957 (cadddr (car x))
958 (caddar x))
959 (let ((varlist (caddar x)))
960 (if (member 'trunc (car x) :test #'eq)
961 (srdisrep x)
962 (cdisrep (cdr x)))))
964 (defun cdisrep (x &aux n d sign)
965 (cond ((pzerop (car x)) 0)
966 ((or (equal 1 (cdr x)) (floatp (cdr x))) (pdisrep (car x)))
967 (t (setq sign (cond ($ratexpand (setq n (pdisrep (car x))) 1)
968 ((pminusp (car x))
969 (setq n (pdisrep (pminus (car x)))) -1)
970 (t (setq n (pdisrep (car x))) 1)))
971 (setq d (pdisrep (cdr x)))
972 (cond ((and (numberp n) (numberp d))
973 (list '(rat) (* sign n) d))
974 ((and $ratdenomdivide $ratexpand
975 (not (atom n))
976 (eq (caar n) 'mplus))
977 (fancydis n d))
978 ((numberp d)
979 (list '(mtimes ratsimp)
980 (list '(rat) sign d) n))
981 ((equal sign -1)
982 (cons '(mtimes ratsimp)
983 (cond ((numberp n)
984 (list (* n -1)
985 (list '(mexpt ratsimp) d -1)))
986 (t (list sign n (list '(mexpt ratsimp) d -1))))))
987 ((equal n 1)
988 (list '(mexpt ratsimp) d -1))
989 (t (list '(mtimes ratsimp) n
990 (list '(mexpt ratsimp) d -1)))))))
993 ;; FANCYDIS GOES THROUGH EACH TERM AND DIVIDES IT BY THE DENOMINATOR.
995 (defun fancydis (n d)
996 (setq d (simplify (list '(mexpt) d -1)))
997 (simplify (cons '(mplus)
998 (mapcar #'(lambda (z) ($ratdisrep (ratf (list '(mtimes) z d))))
999 (cdr n)))))
1002 (defun compatvarl (a b c d)
1003 (cond ((null a) nil)
1004 ((or (null b) (null c) (null d)) (throw 'compatvl nil))
1005 ((alike1 (car a) (car b))
1006 (setq a (compatvarl (cdr a) (cdr b) (cdr c) (cdr d)))
1007 (if (eq (car c) (car d))
1009 (cons (cons (car c) (car d)) a)))
1010 (t (compatvarl a (cdr b) c (cdr d)))))
1012 (defun newvar (l &aux vlist)
1013 (newvar1 l)
1014 (setq varlist (nconc (sortgreat vlist) varlist)))
1016 (defun sortgreat (l)
1017 (and l (nreverse (stable-sort l 'great))));FIXME consider a total order function with #'sort
1019 (defun fnewvar (l &aux (*fnewvarsw t))
1020 (newvar l))
1022 (defun nestlev (exp)
1023 (if (atom exp)
1025 (do ((m (nestlev (cadr exp)) (max m (nestlev (car l))))
1026 (l (cddr exp) (cdr l)))
1027 ((null l) (1+ m)))))
1029 (defun radsort (l)
1030 (sort l #'(lambda (a b)
1031 (let ((na (nestlev a)) (nb (nestlev b)))
1032 (cond ((< na nb) t)
1033 ((> na nb) nil)
1034 (t (great b a)))))))
1036 ;; THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 5
1037 ;; IT INCLUDES THE CONVERSION AND TOP-LEVEL ROUTINES USED
1038 ;; BY THE REST OF THE FUNCTIONS.