src/rat3e.lisp

   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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  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 alc
  20                       vlist radlist expsumsplit *ratsimp*))
  21
  22 (defvar *ratsimp* nil)
  23
  24 (defmvar factorresimp nil "If `t' resimplifies factor(x-y) to x-y")
  25
  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))))
  34
  35 (defprop mrat mrateval mfexpr*)
  36
  37 (defmfun $ratnumer (x)
  38   (cond ((mbagp x)
  39          (cons (car x) (mapcar '$ratnumer (cdr x))))
  40         (t
  41          (setq x (taychk2rat x))
  42          (cons (car x) (cons (cadr x) 1)))))
  43
  44 (defmfun $ratdenom (x)
  45   (cond ((mbagp x)
  46          (cons (car x) (mapcar '$ratdenom (cdr x))))
  47         (t
  48          (setq x (taychk2rat x))
  49          (cons (car x) (cons (cddr x) 1)))))
  50
  51 (defun taychk2rat (x)
  52   (cond ((and ($ratp x)
  53               (member 'trunc (cdar x) :test #'eq))
  54          ($taytorat x))
  55         (t ($rat x))))
  56
  57 (defun tellratdisp (x)
  58   (pdisrep+ (trdisp1 (cdr x) (car x))))
  59
  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)))))
  67
  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)))
  73
  74 (defmfun $tellrat (&rest args)
  75   (mapc #'tellrat1 args)
  76   (unless (null args) (add2lnc 'tellratlist $myoptions))
  77   (cons '(mlist) (mapcar #'tellratdisp tellratlist)))
  78
  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))
  94
  95
  96 (defmfun $printvarlist ()
  97   (cons '(mlist) (copy-tree varlist)))
  98
  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)))))
 105
 106 (defmfun $ratvars (&rest args)
 107   (add2lnc '$ratvars $myoptions)
 108   (setq $ratvars (cons '(mlist simp) (setq varlist (mapfr1 args varlist)))))
 109
 110 (defun mapfr1 (l varlist)
 111   (mapcar #'(lambda (z) (fr1 z varlist)) l))
 112
 113 (defmvar inratsimp nil)
 114
 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)))
 121
 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)))
 129
 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))))))
 135
 136 ;;;VARLIST HAS MAIN VARIABLE AT END
 137
 138 (defun joinvarlist (cdrl)
 139   (mapc #'(lambda (z) (unless (memalike z varlist) (push z varlist)))
 140         (reverse (mapfr1 cdrl nil))))
 141
 142 (defmfun $rat (e &rest vars)
 143   (cond ((not (null vars))
 144          (let (varlist)
 145            (joinvarlist vars)
 146            (lnewvar e)
 147            (rat0 e)))
 148         (t
 149          (lnewvar e)
 150          (rat0 e))))
 151
 152 (defun rat0 (exp)                       ;SIMP FLAGS?
 153   (if (mbagp exp)
 154       (cons (car exp) (mapcar #'rat0 (cdr exp)))
 155       (ratf exp)))
 156
 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))))
 163
 164 ;; $RATSIMP, $FULLRATSIMP and $RAT allow for optional extra
 165 ;; arguments specifying the VARLIST.
 166
 167 ;;;PSQFR HAS NOT BEEN CHANGED TO MAKE USE OF THE SQFR FLAGS YET
 168
 169 (defmfun $sqfr (x)
 170   (let ((varlist (cdr $ratvars)) genvar $keepfloat $ratfac)
 171     (sublis '((factored . sqfred) (irreducible . sqfr))
 172             (ffactor x #'psqfr))))
 173
 174 (declare-top (special fn))
 175
 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))))
 183
 184 (declare-top (special var))
 185
 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)))))
 206
 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)))
 213
 214
 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)))
 221
 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)))
 231
 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))
 270                (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))))
 281
 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)))))
 292                   (t
 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))))
 297
 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))
 307         (t
 308          (cons (cons (caaar n) '(simp factored)) (cdar n)))))
 309
 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)))))))
 314
 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))))))
 323
 324 (declare-top (unspecial var))
 325
 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)))
 335
 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))))))
 343
 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)))))
 372
 373 (defmfun $quotient (&rest args)
 374   (cadr (apply '$divide args)))
 375
 376 (defmfun $remainder (&rest args)
 377   (caddr (apply '$divide args)))
 378
 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)))))
 398
 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)))))
 417
 418 (defun pget (gen)
 419   (cons gen '(1 1)))
 420
 421 (defun m$exp? (x)
 422   (and (mexptp x) (eq (cadr x) '$%e)))
 423
 424 (defun algp ($x)
 425   (algpchk $x nil))
 426
 427 (defun algpget ($x)
 428   (algpchk $x t))
 429
 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)))
 457
 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)))
 464
 465 (defun ratsetup (vl gl)
 466   (ratsetup1 vl gl) (ratsetup2 vl gl))
 467
 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)))
 474
 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))))
 489
 490 (defun porder (p)
 491   (if (pcoefp p) 0 (valget (car p))))
 492
 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))))))
 501
 502 (defun gensym-readable (name)
 503   (cond ((symbolp name)
 504          (gensym (string-trim "$" (string name))))
 505         (t
 506          (setq name (aformat nil "~:M" name))
 507          (if name (gensym name) (gensym)))))
 508
 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))))
 515
 516 (defun creatsym (n)
 517   (dotimes (i n)
 518     (push (gensym) genvar)))
 519
 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)))
 525
 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))
 533         1))
 534
 535 (defun ratrep (x varl)
 536   (setq varlist varl)
 537   (ratrep* x))
 538
 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))))))
 549
 550 (defvar *withinratf* nil)
 551
 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))))))
 560
 561
 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))
 608                     t))
 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))))))
 620
 621
 622 (defun putonvlist (x)
 623   (push x vlist)
 624   (and $algebraic
 625        (setq x (assolike x tellratlist))
 626        (mapc 'newvar1 x)))
 627
 628 (setq expsumsplit t)                ;CONTROLS SPLITTING SUMS IN EXPONS
 629
 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
 637
 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)
 643
 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
 661
 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)))
 667
 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)
 674                                           1
 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))))
 683
 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))))))
 728
 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))))))
 750
 751 (defun newvar3 (x)
 752   (or (memalike x vlist)
 753       (memalike x varlist)
 754       (putonvlist x)))
 755
 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)))
 782
 783 (defun ratsimp (x varlist genvar) ($ratdisrep (ratf x)))
 784
 785 (defun littlefr1 (x)
 786   (cons (remove 'simp (car x) :test #'eq) (mapfr1 (cdr x) nil)))
 787
 788 ;;IF T RATSIMP FACTORS RADICANDS AND LOGANDS
 789 (defmvar fr-factor nil)
 790
 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
 794          x)
 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)))))))))
 803
 804 ;;SIMPLIFY MEXPT'S & RATEXPAND EXPONENT
 805
 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))))
 811       x))
 812
 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))))
 827
 828 ;; This control of conversion from float to rational appears to be explained
 829 ;; nowhere. - RJF
 830
 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))))
 841
 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.
 851              ;;
 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))))))))
 868
 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))))
 875   f)
 876
 877 (defun pdisrep (p)
 878   (if (atom p)
 879       p
 880       (pdisrep+ (pdisrep2 (cdr p) (get (car p) 'disrep)))))
 881
 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))))
 891
 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)))))
 899
 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))))))
 904
 905 (defun pdisrep*chk (a)
 906   (if (mtimesp a) (cdr a) (ncons a)))
 907
 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))))))
 914
 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)))))
 920
 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)))))
 932
 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)))))
 937
 938
 939 (defmvar $ratdenomdivide t)
 940
 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)
 946         (t
 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))))
 953
 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)))))
 963
 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)))))))
 991
 992
 993 ;; FANCYDIS GOES THROUGH EACH TERM AND DIVIDES IT BY THE DENOMINATOR.
 994
 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)))))
1000
1001
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))
1008              a
1009              (cons (cons (car c) (car d)) a)))
1010         (t (compatvarl a (cdr b) c (cdr d)))))
1011
1012 (defun newvar (l &aux vlist)
1013   (newvar1 l)
1014   (setq varlist (nconc (sortgreat vlist) varlist)))
1015
1016 (defun sortgreat (l)
1017   (and l (nreverse (stable-sort l 'great))));FIXME consider a total order function with #'sort
1018
1019 (defun fnewvar (l &aux (*fnewvarsw t))
1020   (newvar l))
1021
1022 (defun nestlev (exp)
1023   (if (atom exp)
1024       0
1025       (do ((m (nestlev (cadr exp)) (max m (nestlev (car l))))
1026            (l (cddr exp) (cdr l)))
1027           ((null l) (1+ m)))))
1028
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)))))))
1035
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.