src/rat3e.lisp

   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         ;; Handle denormalized floats -- see maxima-discuss 2021-04-27 and bug 3777
 872         ((and (not (= x 0.0)) (< (abs x) COMMON-LISP:LEAST-POSITIVE-NORMALIZED-DOUBLE-FLOAT))
 873          (let ((r ($rationalize x)))
 874            (cons (cadr r) (caddr r))))
 875         ((< x 0.0)
 876          (setq x (ration1 (* -1.0 x)))
 877          (rplaca x (* -1 (car x))))
 878         (t (ration1 x))))
 879
 880 (defun ration1 (x)
 881   (let ((rateps (if (not (floatp $ratepsilon))
 882                     ($float $ratepsilon)
 883                     $ratepsilon)))
 884     (or (and (zerop x) (cons 0 1))
 885         (prog (y a)
 886            (return
 887              ;; I (rtoy) think this is computing a continued fraction
 888              ;; expansion of the given float.
 889              ;;
 890              ;; FIXME?  CMUCL used to use this routine for its
 891              ;; RATIONALIZE function, but there were known bugs in
 892              ;; that implementation, where the result was not very
 893              ;; accurate.  Unfortunately, I can't find the example
 894              ;; that demonstrates this.  In any case, CMUCL replaced
 895              ;; it with an algorithm based on the code in Clisp, which
 896              ;; was much better.
 897              (do ((xx x (setq y (/ 1.0 (- xx (float a x)))))
 898                   (num (setq a (floor x)) (+ (* (setq a (floor y)) num) onum))
 899                   (den 1 (+ (* a den) oden))
 900                   (onum 1 num)
 901                   (oden 0 den))
 902                  ((or (zerop (- xx (float a x)))
 903                       (and (not (zerop den))
 904                            (not (> (abs (/ (- x (/ (float num x) (float den x))) x)) rateps))))
 905                   (cons num den))))))))
 906
 907 (defun prepfloat (f)
 908   (cond (modulus (merror (intl:gettext "rat: can't rationalize ~M when modulus = ~M") f modulus))
 909         ($ratprint (mtell (intl:gettext "~&rat: replaced ~A by") f)))
 910   (setq f (maxima-rationalize f))
 911   (when $ratprint
 912     (mtell " ~A/~A = ~A~%"  (car f) (cdr f) (fpcofrat1 (car f) (cdr f))))
 913   f)
 914
 915 (defun pdisrep (p)
 916   (if (atom p)
 917       p
 918       (pdisrep+ (pdisrep2 (cdr p) (get (car p) 'disrep)))))
 919
 920 (defun pdisrep! (n var)
 921   (cond ((zerop n) 1)
 922         ((equal n 1) (cond ((atom var) var)
 923                            ((or (eq (caar var) 'mtimes)
 924                                 (eq (caar var) 'mplus))
 925                             (copy-list var))
 926                            (t var)))
 927         ((eql var 1) 1)
 928         (t (list '(mexpt ratsimp) var n))))
 929
 930 (defun pdisrep+ (p)
 931   (cond ((null (cdr p)) (car p))
 932         (t (let ((a (last p)))
 933              (cond ((mplusp (car a))
 934                     (rplacd a (cddar a))
 935                     (rplaca a (cadar a))))
 936              (cons '(mplus ratsimp) p)))))
 937
 938 (defun pdisrep* (a b)
 939   (cond ((equal a 1) b)
 940         ((equal b 1) a)
 941         (t (cons '(mtimes ratsimp) (nconc (pdisrep*chk a) (pdisrep*chk b))))))
 942
 943 (defun pdisrep*chk (a)
 944   (if (mtimesp a) (cdr a) (ncons a)))
 945
 946 (defun pdisrep2 (p var)
 947   (cond ((null p) nil)
 948         ($ratexpand (pdisrep2expand p var))
 949         (t (do ((l () (cons (pdisrep* (pdisrep (cadr p)) (pdisrep! (car p) var)) l))
 950                 (p p (cddr p)))
 951                ((null p) (nreverse l))))))
 952
 953 ;; IF $RATEXPAND IS TRUE, (X+1)*(Y+1) WILL DISPLAY AS
 954 ;; XY + Y + X + 1  OTHERWISE, AS (X+1)Y + X + 1
 955 (defmvar $ratexpand nil)
 956
 957 (defmfun $ratexpand (x)
 958   (if (mbagp x)
 959       (cons (car x) (mapcar '$ratexpand (cdr x)))
 960       (let (($ratexpand t) ($ratfac nil))
 961         (ratdisrep (ratf x)))))
 962
 963 (defun pdisrep*expand (a b)
 964   (cond ((equal a 1) (list b))
 965         ((equal b 1) (list a))
 966         ((or (atom a) (not (eq (caar a) 'mplus)))
 967          (list (cons (quote (mtimes ratsimp))
 968                      (nconc (pdisrep*chk a) (pdisrep*chk b)))))
 969         (t (mapcar #'(lambda (z) (if (equal z 1) b
 970                                      (cons '(mtimes ratsimp)
 971                                            (nconc (pdisrep*chk z)
 972                                                   (pdisrep*chk b)))))
 973                    (cdr a)))))
 974
 975 (defun pdisrep2expand (p var)
 976   (cond ((null p) nil)
 977         (t (nconc (pdisrep*expand (pdisrep (cadr p)) (pdisrep! (car p) var))
 978                   (pdisrep2expand (cddr p) var)))))
 979
 980
 981 (defmvar $ratdenomdivide t)
 982
 983 (defmfun $ratdisrep (x)
 984   (cond ((mbagp x)
 985          ;; Distribute over lists, equations, and matrices.
 986          (cons (car x) (mapcar #'$ratdisrep (cdr x))))
 987         ((not ($ratp x)) x)
 988         (t
 989          (setq x (ratdisrepd x))
 990          (if (and (not (atom x))
 991                   (member 'trunc (cdar x) :test #'eq))
 992            (cons (delete 'trunc (copy-list (car x)) :count 1 :test #'eq)
 993                  (cdr x))
 994            x))))
 995
 996 ;; RATDISREPD is needed by DISPLA. - JPG
 997 (defun ratdisrepd (x)
 998   (mapc #'(lambda (y z) (putprop y z (quote disrep)))
 999         (cadddr (car x))
1000         (caddar x))
1001   (let ((varlist (caddar x)))
1002     (if (member 'trunc (car x) :test #'eq)
1003         (srdisrep x)
1004         (cdisrep (cdr x)))))
1005
1006 (defun cdisrep (x &aux n d sign)
1007   (cond ((pzerop (car x)) 0)
1008         ((or (equal 1 (cdr x)) (floatp (cdr x))) (pdisrep (car x)))
1009         (t (setq sign (cond ($ratexpand (setq n (pdisrep (car x))) 1)
1010                             ((pminusp (car x))
1011                              (setq n (pdisrep (pminus (car x)))) -1)
1012                             (t (setq n (pdisrep (car x))) 1)))
1013            (setq d (pdisrep (cdr x)))
1014            (cond ((and (numberp n) (numberp d))
1015                   (list '(rat) (* sign n) d))
1016                  ((and $ratdenomdivide $ratexpand
1017                        (not (atom n))
1018                        (eq (caar n) 'mplus))
1019                   (fancydis n d))
1020                  ((numberp d)
1021                   (list '(mtimes ratsimp)
1022                         (list '(rat) sign d) n))
1023                  ((equal sign -1)
1024                   (cons '(mtimes ratsimp)
1025                         (cond ((numberp n)
1026                                (list (* n -1)
1027                                      (list '(mexpt ratsimp) d -1)))
1028                               (t (list sign n (list '(mexpt ratsimp) d -1))))))
1029                  ((equal n 1)
1030                   (list '(mexpt ratsimp) d -1))
1031                  (t (list '(mtimes ratsimp) n
1032                           (list '(mexpt ratsimp) d -1)))))))
1033
1034
1035 ;; FANCYDIS GOES THROUGH EACH TERM AND DIVIDES IT BY THE DENOMINATOR.
1036
1037 (defun fancydis (n d)
1038   (setq d (simplify (list '(mexpt) d -1)))
1039   (simplify (cons '(mplus)
1040                   (mapcar #'(lambda (z) ($ratdisrep (ratf (list '(mtimes) z d))))
1041                           (cdr n)))))
1042
1043
1044 (defun compatvarl (a b c d)
1045   (cond ((null a) nil)
1046         ((or (null b) (null c) (null d)) (throw 'compatvl nil))
1047         ((alike1 (car a) (car b))
1048          (setq a (compatvarl (cdr a) (cdr b) (cdr c) (cdr d)))
1049          (if (eq (car c) (car d))
1050              a
1051              (cons (cons (car c) (car d)) a)))
1052         (t (compatvarl a (cdr b) c (cdr d)))))
1053
1054 (defun newvar (l &aux vlist)
1055   (newvar1 l)
1056   (setq varlist (nconc (sortgreat vlist) varlist)))
1057
1058 (defun sortgreat (l)
1059   (and l (nreverse (stable-sort l 'great))));FIXME consider a total order function with #'sort
1060
1061 (defun fnewvar (l &aux (*fnewvarsw t))
1062   (newvar l))
1063
1064 (defun nestlev (exp)
1065   (if (atom exp)
1066       0
1067       (do ((m (nestlev (cadr exp)) (max m (nestlev (car l))))
1068            (l (cddr exp) (cdr l)))
1069           ((null l) (1+ m)))))
1070
1071 (defun radsort (l)
1072   (sort l #'(lambda (a b)
1073               (let ((na (nestlev a)) (nb (nestlev b)))
1074                 (cond ((< na nb) t)
1075                       ((> na nb) nil)
1076                       (t (great b a)))))))
1077
1078 ;;      THIS IS THE END OF THE NEW RATIONAL FUNCTION PACKAGE PART 5
1079 ;;      IT INCLUDES THE CONVERSION AND TOP-LEVEL ROUTINES USED
1080 ;;      BY THE REST OF THE FUNCTIONS.