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 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
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")
27 (let ((varlist (caddar x
)))
28 (cond ((and evp $infeval
) (meval ($ratdisrep x
)))
30 (and $float $keepfloat
)
31 (not (alike varlist
(mapcar #'meval varlist
))))
32 (ratf (meval ($ratdisrep x
))))
35 (defprop mrat mrateval mfexpr
*)
37 (defmfun $ratnumer
(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)
46 (cons (car x
) (mapcar '$ratdenom
(cdr x
))))
48 (setq x
(taychk2rat x
))
49 (cons (car x
) (cons (cddr x
) 1)))))
53 (member 'trunc
(cdar x
) :test
#'eq
))
57 (defun tellratdisp (x)
58 (pdisrep+ (trdisp1 (cdr x
) (car x
))))
60 (defun trdisp1 (p var
)
62 (t (cons (pdisrep* (if (mtimesp (cadr p
))
64 (cadr p
)) ;prevents clobbering p
65 (pdisrep! (car p
) var
))
66 (trdisp1 (cddr p
) var
)))))
68 (defmfun $untellrat
(&rest 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
)))
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)
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
)
117 loop
(setq exp1
(simplify (apply #'$ratsimp
(cons exp argl
))))
118 (when (alike1 exp exp1
) (return exp
))
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
))
130 (defmfun $totaldisrep
(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
))
152 (defun rat0 (exp) ;SIMP FLAGS?
154 (cons (car exp
) (mapcar #'rat0
(cdr exp
)))
157 (defmfun ($ratsimp
:properties
((evfun t
))) (e &rest vars
)
158 (cond ((not (null vars
))
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
170 (let ((varlist (cdr $ratvars
)) genvar $keepfloat $ratfac
)
171 (sublis '((factored . sqfred
) (irreducible . sqfr
))
172 (ffactor x
#'psqfr
))))
174 (declare-top (special fn
))
177 (cond ((and (mexptp p
) (integerp (caddr p
)))
178 (list '(mexpt) (whichfn (cadr p
)) (caddr p
)))
180 (cons '(mtimes) (mapcar #'whichfn
(cdr p
))))
181 (fn (ffactor p
#'pfactor
))
184 (declare-top (special var
))
187 (prog (alc ans adn
* $gcd
)
188 (setq $gcd
'$algebraic
)
189 (when (or (atom p
) (numberp p
)) (return p
))
191 (when (and (not $nalgfac
) (not intbs
*))
192 (setq intbs
* (findibase minpoly
*)))
194 (setq ans
(ffactor p
#'pfactor
))
195 (cond ((eq (caar ans
) 'mplus
)
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
*)
204 (if alc
(pdis alc
) 1)))
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
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))
232 (defun factor (e &optional
(mp nil mp?
))
233 (let ((tellratlist nil
)
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
*)))
246 cargs
(if mp?
(list mp
) nil
))
247 (when (symbolp p
) (return p
))
249 (return (let (($factorflag
(not scanmapp
)))
252 (return (cons (car p
) (mapcar #'(lambda (x) (apply #'factor
(cons x cargs
))) (cdr p
)))))
254 (setq *alpha
* (meqhk mp
))
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))
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
)))
268 ($tellrat
(pdis minpoly
*))
272 (unless scanmapp
(setq p
(let (($ratfac t
)) (sratsimp p
))))
274 (when (symbolp p
) (return p
))
276 (return (let (($factorflag
(not scanmapp
)))
278 (setq $negdistrib nil
)
279 (setq p
(let ($factorflag
($ratexpand $facexpand
))
282 (setq p
(let (($expop
0) ($expon
0))
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
)))
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
))))
298 (defun factornumber (n)
299 (setq n
(nretfactor1 (nratfact (cdr ($rat n
)))))
301 (cons '(mtimes simp factored
)
302 (if (equal (car n
) -
1)
303 (cons (car n
) (nreverse (cdr n
)))
308 (cons (cons (caaar n
) '(simp factored
)) (cdar n
)))))
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)
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
))
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
)))
337 (if (mbagp e
) (cons (car e
) (mapcar 'mod1
(cdr 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))
348 (when (and ($ratp y
) (setq formflag t
) (integerp (cadr y
)) (equal (cddr y
) 1))
350 (when (and (integerp x
) (integerp y
))
352 (merror (intl:gettext
"divide: Quotient by zero")))
353 (return (cons '(mlist) (multiple-value-list (truncate x y
)))))
355 (mapc #'newvar
(reverse (cdr $ratvars
)))
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
))))
365 (setq x
(ptimes (car x
) (cdr y
)))
366 (setq x
(pdivide x
(car y
)))
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
))
385 (when (numberp v
) (improper-arg-err v
'$gcd
)))
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
))
394 (setq y
(cdr (ratrep* y
)))
395 on
(setq x
(cons (pgcd (car x
) (car y
)) (plcm (cdr x
) (cdr y
))))
397 (return (if formflag h
(ratdisrep h
)))))
399 (defmfun $content
(x &rest vars
)
400 (prog (y h varlist formflag
)
401 (setq formflag
($ratp x
))
404 (desetq (h x . y
) (ratrep* 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
)))))
413 (setq h
(list '(mlist)
414 (cons h
(rattimes (car x
) y nil
))
416 (return (if formflag h
($totaldisrep h
)))))
422 (and (mexptp x
) (eq (cadr x
) '$%e
)))
430 (defun algpchk ($x mpflag
&aux temp
)
431 (cond ((eq $x
'$%i
) '(2 -
1))
432 ((eq $x
'$%phi
) '(2 1 1 -
1 0 -
1))
435 (let ((r (prep1 (cadr $x
))))
436 (cond ((onep1 (cdr r
)) ;INTEGRAL ALG. QUANT.?
437 (list (caddr (caddr $x
))
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
)))))
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
)))))
454 (if (and (= (pt-le temp
) 1) (setq $x
(assol $x genpairs
)))
455 (rplacd $x
(cons (cadr temp
) 1)))
458 (defun radfunp (x funcflag
) ;FUNCFLAG -> TEST FOR ALG FUNCTION NOT NUMBER
460 ((not (eq (caar x
) 'mexpt
)) nil
)
461 ((not (ratnump (caddr x
))) nil
)
462 (funcflag (not (numberp (cadr x
))))
465 (defun ratsetup (vl gl
)
466 (ratsetup1 vl gl
) (ratsetup2 vl gl
))
468 (defun ratsetup1 (vl gl
)
470 (mapc #'(lambda (v g
)
471 (setq v
(assolike v
*ratweights
))
472 (if v
(putprop g v
'$ratweight
) (remprop g
'$ratweight
)))
475 (defun ratsetup2 (vl gl
)
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
))))
483 (and $ratfac
(let ($ratfac
)
484 (mapc #'(lambda (v g
)
486 (putprop g
(car (prep1 v
)) 'unhacked
)
487 (remprop g
'unhacked
)))
491 (if (pcoefp p
) 0 (valget (car p
))))
493 (defun algordset (x gl
)
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)
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))))
518 (push (gensym) genvar
)))
520 (defun prenumber (v n
)
524 (setf (symbol-value (car vl
)) i
)))
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
))
535 (defun ratrep (x varl
)
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
))
550 (defvar *withinratf
* nil
)
553 (prog (u *withinratf
*)
554 (setq *withinratf
* t
)
555 (when (eq '%%
(catch 'ratf
(newvar l
)))
556 (setq *withinratf
* nil
)
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
)
564 (cond ($keepfloat
(cons x
1.0))
566 ((integerp x
) (cons (cmod x
) 1))
569 (cons (numerator x
) (denominator x
))
570 (cquotient (numerator x
) (denominator x
))))
571 ((atom x
) (cond ((assolike x genpairs
))
573 ((and $ratfac
(assolike x genpairs
)))
574 ((eq (caar x
) 'mplus
)
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
)))
582 (t (do ((a (prep1 (cadr x
)) (ratplus a
(prep1 (car l
))))
583 (l (cddr x
) (cdr l
)))
585 ((eq (caar x
) 'mtimes
)
586 (do ((a (savefactors (prep1 (cadr x
)))
587 (rattimes a
(savefactors (prep1 (car l
))) sw
))
589 (sw (not (and $norepeat
(member 'ratsimp
(cdar x
) :test
#'eq
)))))
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
))))
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
))
603 (cond ((and *withinratf
* (member 'trunc
(car x
) :test
#'eq
))
607 (setq temp
(compatvarl (caddar x
) varlist
(cadddr (car x
)) genvar
))
609 (cond ((member 'trunc
(car x
) :test
#'eq
)
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
))
622 (defun putonvlist (x)
625 (setq x
(assolike x tellratlist
))
628 (setq expsumsplit t
) ;CONTROLS SPLITTING SUMS IN EXPONS
630 (defun newvarmexpt (x e flag
)
631 ;; WHEN FLAG IS T, CALL RETURNS RATFORM
633 (when (and (integerp e
) (not flag
))
634 (return (newvar1 (cadr x
))))
638 ;; X=B^N FOR N A NUMBER
640 (setq topexp
(* topexp e
))
644 ;; X=B^(P/Q) FOR P AND Q INTEGERS
646 (cond ((or (minusp (cadr e
)) (> (cadr e
) 1))
647 (setq topexp
(* topexp
(cadr e
)))
648 (setq x
(list '(mexpt)
650 (list '(rat) 1 (caddr e
))))))
651 (cond ((or flag
(numberp (cadr x
)) ))
653 (cond ((memalike x radlist
) (return nil
))
654 (t (setq radlist
(cons x radlist
))
655 (return (newvar1 (cadr x
))))) )
656 ($algebraic
(newvar1 (cadr x
)))))
658 ((eq (caar e
) 'mtimes
)
665 (setq topexp
(* topexp
(cadr e
)))
669 (and (not (atom (cadr e
)))
671 (not (equal 1 (cadadr e
)))
672 (setq topexp
(* topexp
(cadadr e
)))
673 (setq e
(cons (list '(rat)
680 (setq e
(simplify (cons '(mtimes)
685 ((and (eq (caar e
) 'mplus
) expsumsplit
) ;SWITCH CONTROLS
686 (setq ;SPLITTING EXPONENT
694 (simplify (list '(mtimes)
698 (cond (flag (return (prep1 x
)))
699 (t (return (newvar1 x
))))))
703 (not (eq (caar x
) 'mexpt
)))
705 ((or (memalike x varlist
) (memalike x vlist
))
707 (t (cond ((or (atom x
) (null *fnewvarsw
))
709 (t (setq x
(littlefr1 x
))
710 (mapc #'newvar1
(cdr x
))
711 (or (memalike x vlist
)
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
))
727 (t (ratexpt (prep1 x
) topexp
))))))
730 (cond ((numberp x
) nil
)
731 ((memalike x varlist
) nil
)
732 ((memalike x vlist
) nil
)
733 ((atom x
) (putonvlist 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
))
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
)
749 (t (putonvlist x
))))))
752 (or (memalike x vlist
)
756 (defun fr1 (x varlist
) ;put radicands on initial varlist?
757 (prog (genvar $norepeat
*ratsimp
* radlist vlist nvarlist ovarlist genpairs
)
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
)
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
)
776 (mapcar #'(lambda (x y
) (cons x
(prep1 y
))) ovarlist nvarlist
))
777 (setq x
(rdis (prep1 x
)))
778 (cond (radlist ;rational radicands
780 (setq x
(ratsimp (simplify x
) nil nil
)))))
783 (defun ratsimp (x varlist genvar
) ($ratdisrep
(ratf 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 ?
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
))
802 (mapfr1 (cdr x
) varlist
)))))))))
804 ;;SIMPLIFY MEXPT'S & RATEXPAND EXPONENT
807 (if (and (mexptp x
) (not (atom (caddr x
))))
808 (list (car x
) (cadr x
)
809 (let ((varlist varlist
) *ratsimp
*)
810 ($ratexpand
(caddr x
))))
815 (when (setq g
(assolike e genpairs
))
817 (setq g
(gensym-readable e
))
818 (putprop g e
'disrep
)
820 (push (cons e
(rget g
)) genpairs
)
821 (valput g
(if genvar
(1- (valget (car genvar
))) 1))
823 (when (setq p
(and $algebraic
(algpget e
)))
825 (putprop g p
'tellrat
))
828 ;; This control of conversion from float to rational appears to be explained
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
))))
838 (setq x
(ration1 (* -
1.0 x
)))
839 (rplaca x
(* -
1 (car x
))))
843 (let ((rateps (if (not (floatp $ratepsilon
))
846 (or (and (zerop x
) (cons 0 1))
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
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
))
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
))))))))
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
))
874 (mtell " ~A/~A = ~A~%" (car f
) (cdr f
) (fpcofrat1 (car f
) (cdr f
))))
880 (pdisrep+ (pdisrep2 (cdr p
) (get (car p
) 'disrep
)))))
882 (defun pdisrep! (n var
)
884 ((equal n
1) (cond ((atom var
) var
)
885 ((or (eq (caar var
) 'mtimes
)
886 (eq (caar var
) 'mplus
))
890 (t (list '(mexpt ratsimp
) var n
))))
893 (cond ((null (cdr p
)) (car p
))
894 (t (let ((a (last p
)))
895 (cond ((mplusp (car a
))
897 (rplaca a
(cadar a
))))
898 (cons '(mplus ratsimp
) p
)))))
900 (defun pdisrep* (a b
)
901 (cond ((equal a
1) b
)
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
)
910 ($ratexpand
(pdisrep2expand p var
))
911 (t (do ((l () (cons (pdisrep* (pdisrep (cadr p
)) (pdisrep! (car p
) var
)) l
))
913 ((null p
) (nreverse l
))))))
915 (defmfun ($ratexpand
:properties
((evfun t
))) (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
)
933 (defun pdisrep2expand (p var
)
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)
943 ;; Distribute over lists, equations, and matrices.
944 (cons (car x
) (mapcar #'$ratdisrep
(cdr 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
)
954 ;; RATDISREPD is needed by DISPLA. - JPG
955 (defun ratdisrepd (x)
956 (mapc #'(lambda (y z
) (putprop y z
(quote disrep
)))
959 (let ((varlist (caddar x
)))
960 (if (member 'trunc
(car x
) :test
#'eq
)
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)
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
976 (eq (caar n
) 'mplus
))
979 (list '(mtimes ratsimp
)
980 (list '(rat) sign d
) n
))
982 (cons '(mtimes ratsimp
)
985 (list '(mexpt ratsimp
) d -
1)))
986 (t (list sign n
(list '(mexpt ratsimp
) d -
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
))))
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
)
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
))
1022 (defun nestlev (exp)
1025 (do ((m (nestlev (cadr exp
)) (max m
(nestlev (car l
))))
1026 (l (cddr exp
) (cdr l
)))
1027 ((null l
) (1+ m
)))))
1030 (sort l
#'(lambda (a b
)
1031 (let ((na (nestlev a
)) (nb (nestlev b
)))
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.