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