src/outmis.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
   9 (in-package :maxima)
  10
  11 ;;      ** (c) Copyright 1982 Massachusetts Institute of Technology **
  12
  13 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  14 ;;;                                                                ;;;
  15 ;;;                Miscellaneous Out-of-core Files                 ;;;
  16 ;;;                                                                ;;;
  17 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  18
  19 (macsyma-module outmis)
  20
  21
  22 (declare-top (special $exptisolate))
  23
  24 (defmvar $exptisolate nil
  25   nil
  26   :properties ((evflag t)))
  27 (defmvar $isolate_wrt_times nil
  28   nil
  29   :properties ((evflag t)))
  30
  31 (defmfun $isolate (e *xvar)
  32   (iso1 e (getopr *xvar)))
  33
  34 (defun iso1 (e *xvar)
  35   (cond ((specrepp e) (iso1 (specdisrep e) *xvar))
  36         ((and (free e 'mplus) (or (null $isolate_wrt_times) (free e 'mtimes))) e)
  37         ((freeof *xvar e) (mgen2 e))
  38         ((alike1 *xvar e) *xvar)
  39         ((member (caar e) '(mplus mtimes) :test #'eq) (iso2 e *xvar))
  40         ((eq (caar e) 'mexpt)
  41          (cond ((null (atom (cadr e))) (list (car e) (iso1 (cadr e) *xvar) (caddr e)))
  42                ((or (alike1 (cadr e) *xvar) (not $exptisolate)) e)
  43                (t (let ((x ($rat (caddr e) *xvar)) (u 0) (h 0))
  44                     (setq u (ratdisrep ($ratnumer x)) x (ratdisrep ($ratdenom x)))
  45                     (if (not (equal x 1))
  46                         (setq u ($multthru (list '(mexpt) x -1) u)))
  47                     (if (mplusp u)
  48                         (setq u ($partition u *xvar) h (cadr u) u (caddr u)))
  49                     (setq u (power* (cadr e) (iso1 u *xvar)))
  50                     (cond ((not (equal h 0))
  51                            (mul2* (mgen2 (power* (cadr e) h)) u))
  52                           (t u))))))
  53         (t (cons (car e) (mapcar #'(lambda (e1) (iso1 e1 *xvar)) (cdr e))))))
  54
  55 (defun iso2 (e *xvar)
  56   (prog (hasit doesnt op)
  57      (setq op (ncons (caar e)))
  58      (do ((i (cdr e) (cdr i))) ((null i))
  59        (cond ((freeof *xvar (car i)) (setq doesnt (cons (car i) doesnt)))
  60              (t (setq hasit (cons (iso1 (car i) *xvar) hasit)))))
  61      (cond ((null doesnt) (go ret))
  62            ((and (null (cdr doesnt)) (atom (car doesnt))) (go ret))
  63            ((prog2 (setq doesnt (simplify (cons op doesnt)))
  64                 (and (free doesnt 'mplus)
  65                      (or (null $isolate_wrt_times)
  66                          (free doesnt 'mtimes)))))
  67            (t (setq doesnt (mgen2 doesnt))))
  68      (setq doesnt (ncons doesnt))
  69      ret  (return (simplifya (cons op (nconc hasit doesnt)) nil))))
  70
  71 (defun mgen2 (h)
  72   (cond ((memsimilarl h (cdr $labels) (getlabcharn $linechar)))
  73         (t (setq h (displine h)) (and $dispflag (mterpri)) h)))
  74
  75 (defun memsimilarl (item list linechar)
  76   (cond ((null list) nil)
  77         ((and (char= (getlabcharn (car list)) linechar)
  78               (boundp (car list))
  79               (memsimilar item (car list) (symbol-value (car list)))))
  80         (t (memsimilarl item (cdr list) linechar))))
  81
  82 (defun memsimilar (item1 item2 item2ev)
  83   (cond ((equal item2ev 0) nil)
  84         ((alike1 item1 item2ev) item2)
  85         (t (let ((errorsw t) r)
  86              (setq r (catch 'errorsw (div item2ev item1)))
  87              (and (mnump r) (not (zerop1 r)) (div item2 r))))))
  88
  89 (defmfun $pickapart (x lev)
  90   (setq x (format1 x))
  91   (cond ((not (fixnump lev))
  92          (merror (intl:gettext "pickapart: second argument must be an integer; found: ~M") lev))
  93         ((or (atom x) (and (eq (caar x) 'mminus) (atom (cadr x)))) x)
  94         ((= lev 0) (mgen2 x))
  95         ((and (atom (cdr x)) (cdr x)) x)
  96         (t (cons (car x) (mapcar #'(lambda (y) ($pickapart y (1- lev))) (cdr x))))))
  97
  98 (defmfun $reveal (e lev)
  99   (setq e (format1 e))
 100   (if (and (fixnump lev) (plusp lev))
 101       (reveal e 1 lev)
 102       (merror (intl:gettext "reveal: second argument must be a positive integer; found: ~M") lev)))
 103
 104 (defun simple (x)
 105   (or (atom x) (member (caar x) '(rat bigfloat) :test #'eq)))
 106
 107 (defun reveal (e nn lev)
 108   (cond ((simple e) e)
 109         ((= nn lev)
 110          (cond ((eq (caar e) 'mplus) (cons '(|$Sum| simp) (ncons (length (cdr e)))))
 111                ((eq (caar e) 'mtimes) (cons '(|$Product| simp) (ncons (length (cdr e)))))
 112                ((eq (caar e) 'mexpt) '|$Expt|)
 113                ((eq (caar e) 'mquotient) '|$Quotient|)
 114                ((eq (caar e) 'mminus) '|$Negterm|)
 115                ((eq (caar e) 'mlist)
 116                 (cons '(|$List| simp) (ncons (length (cdr e)))))
 117                (t (getop (mop e)))))
 118         (t (let ((u (cond ((member 'simp (cdar e) :test #'eq) (car e))
 119                           (t (cons (caar e) (cons 'simp (cdar e))))))
 120                  (v (mapcar #'(lambda (x) (reveal (format1 x) (1+ nn) lev))
 121                             (margs e))))
 122              (cond ((eq (caar e) 'mqapply) (cons u (cons (cadr e) v)))
 123                    ((eq (caar e) 'mplus) (cons u (nreverse v)))
 124                    (t (cons u v)))))))
 125
 126 (defmspec $properties (x)
 127   (setq x (getopr (fexprcheck x)))
 128   (unless (or (symbolp x) (stringp x))
 129     (merror
 130       (intl:gettext "properties: argument must be a symbol or a string.")))
 131   (let ((u (properties x)) (v (or (safe-get x 'noun) (safe-get x 'verb))))
 132     (if v (nconc u (cdr (properties v))) u)))
 133
 134 (defun properties (x)
 135   (if (stringp x)
 136     ; AT THIS POINT WE MIGHT WANT TO TRY TO TEST ALL CHARS IN STRING ...
 137     (if (and (> (length x) 0) (member (char x 0) *alphabet*))
 138       '((mlist) $alphabetic)
 139       '((mlist)))
 140     (do ((y (symbol-plist x) (cddr y))
 141          (l (cons '(mlist simp) (and (boundp x)
 142                                    (if (optionp x) (ncons "system value")
 143                                        (ncons '$value)))))
 144          (prop))
 145         ((null y)
 146          (if (member x (cdr $features) :test #'eq) (nconc l (ncons '$feature)))
 147          (if (member x (cdr $contexts) :test #'eq) (nconc l (ncons '$context)))
 148          (if (member x (cdr $activecontexts) :test #'eq)
 149            (nconc l (ncons '$activecontext)))
 150          (cond  ((null (symbol-plist x))
 151                (if (fboundp x) (nconc l (list "system function")))))
 152          l)
 153
 154       ;; TOP-LEVEL PROPERTIES
 155       (cond ((setq prop
 156                    (assoc (car y)
 157                           `((bindtest . $bindtest)
 158                             (sp2 . $deftaylor)
 159                             (sp2subs . $deftaylor)
 160                             (assign . "assign property")
 161                             (nonarray . $nonarray)
 162                             (grad . $gradef)
 163                             (integral . $integral)
 164                             (distribute_over . "distributes over bags")
 165                             (simplim%function . "limit function")
 166                             (conjugate-function . "conjugate function")
 167                             (commutes-with-conjugate . "mirror symmetry")
 168                             (risplit-function . "complex characteristic")
 169                             (noun . $noun)
 170                             (evfun . $evfun)
 171                             (evflag . $evflag)
 172                             (op . $operator)) :test #'eq))
 173            (nconc l (ncons (cdr prop))))
 174           ((setq prop (member (car y) opers :test #'eq))
 175            (nconc l (list (car prop))))
 176           ((and (eq (car y) 'operators) (not (or (eq (cadr y) 'simpargs1) (eq (cadr y) nil))))
 177            (nconc l (list '$rule)))
 178           ((and (member (car y) '(fexpr fsubr mfexpr*s mfexpr*) :test #'eq)
 179                 (nconc l (ncons "special evaluation form"))
 180                 nil))
 181           ((and (or (get (car y) 'mfexpr*) (fboundp x))
 182                 ;; Do not add more than one entry to the list.
 183                 (not (member '$transfun l))
 184                 (not (member '$rule l))
 185                 (not (member "system function" l :test #'equal)))
 186            (nconc l
 187                   (list (cond ((get x 'translated) '$transfun)
 188                               ((mgetl x '($rule ruleof)) '$rule)
 189                               (t "system function")))))
 190           ((and (eq (car y) 'autoload)
 191                 (not (member "system function" l :test #'equal)))
 192            (nconc l (ncons (if (member x (cdr $props) :test #'eq)
 193                                "user autoload function"
 194                                "system function"))))
 195           ((and (eq (car y) 'reversealias)
 196                 (member (car y) (cdr $aliases) :test #'eq))
 197            (nconc l (ncons '$alias)))
 198           ((eq (car y) 'data)
 199            (nconc l (cons "database info" (cdr ($facts x)))))
 200           ((eq (car y) 'mprops)
 201            ;; PROPS PROPERTIES
 202            (do ((y
 203                  (cdadr y)
 204                  (cddr y)))
 205                ((null y))
 206              (cond ((setq prop (assoc (car y)
 207                                      `((mexpr . $function)
 208                                        (mmacro . $macro)
 209                                        (hashar . "hashed array")
 210                                        (aexpr . "array function")
 211                                        (atvalues . $atvalue)
 212                                        ($atomgrad . $atomgrad)
 213                                        ($numer . $numer)
 214                                        (depends . $dependency)
 215                                        ($nonscalar . $nonscalar)
 216                                        ($scalar . $scalar)
 217                                        (matchdeclare . $matchdeclare)
 218                                        (mode . $modedeclare)) :test #'eq))
 219                     (nconc l (list (cdr prop))))
 220                    ((eq (car y) 'array)
 221                     (nconc l
 222                            (list (cond ((get x 'array) "complete array")
 223                                        (t "declared array")))))
 224                    ((and (eq (car y) '$props) (cdadr y))
 225                     (nconc l
 226                            (do ((y (cdadr y) (cddr y))
 227                                 (l (list '(mlist) "user properties")))
 228                                ((null y) (list l))
 229                              (nconc l (list (car y)))))))))))))
 230
 231 (defmspec $propvars (x)
 232   (setq x (fexprcheck x))
 233   (do ((iteml (cdr $props) (cdr iteml)) (propvars (ncons '(mlist))))
 234       ((null iteml) propvars)
 235     (and (among x (meval (list '($properties) (car iteml))))
 236          (nconc propvars (ncons (car iteml))))))
 237
 238 (defmspec $printprops (r) (setq r (cdr r))
 239           (if (null (cdr r)) (merror (intl:gettext "printprops: requires two arguments.")))
 240           (let ((s (cadr r)))
 241             (setq r (car r))
 242             (setq r (cond ((atom r)
 243                            (cond ((eq r '$all)
 244                                   (cond ((eq s '$gradef) (mapcar 'caar (cdr $gradefs)))
 245                                         (t (cdr (meval (list '($propvars) s))))))
 246                                  (t (ncons r))))
 247                           (t (cdr r))))
 248             (cond ((eq s '$atvalue) (dispatvalues r))
 249                   ((eq s '$atomgrad) (dispatomgrads r))
 250                   ((eq s '$gradef) (dispgradefs r))
 251                   ((eq s '$matchdeclare) (dispmatchdeclares r))
 252                   (t (merror (intl:gettext "printprops: unknown property ~:M") s)))))
 253
 254 (defun dispatvalues (l)
 255   (do ((l l (cdr l)))
 256       ((null l))
 257     (do ((ll (mget (car l) 'atvalues) (cdr ll)))
 258         ((null ll))
 259       (mtell-open "~M~%"
 260                   (list '(mlabel) nil
 261                         (list '(mequal)
 262                               (atdecode (car l) (caar ll) (cadar ll)) (caddar ll))))))
 263   '$done)
 264
 265 (defun atdecode (fun dl vl)
 266   (setq vl (copy-list vl))
 267   (atvarschk vl)
 268   (let ((eqs nil) (nvarl nil))
 269     (cond ((not (member nil (mapcar #'(lambda (x) (signp e x)) dl) :test #'eq))
 270            (do ((vl vl (cdr vl)) (varl atvars (cdr varl)))
 271                ((null vl))
 272              (and (eq (car vl) munbound) (rplaca vl (car varl))))
 273            (cons (list fun) vl))
 274           (t (setq fun (cons (list fun)
 275                              (do ((n (length vl) (1- n))
 276                                   (varl atvars (cdr varl))
 277                                   (l nil (cons (car varl) l)))
 278                                  ((zerop n) (nreverse l)))))
 279              (do ((vl vl (cdr vl)) (varl atvars (cdr varl)))
 280                  ((null vl))
 281                (and (not (eq (car vl) munbound))
 282                     (setq eqs (cons (list '(mequal) (car varl) (car vl)) eqs))))
 283              (setq eqs (cons '(mlist) (nreverse eqs)))
 284              (do ((varl atvars (cdr varl)) (dl dl (cdr dl)))
 285                  ((null dl) (setq nvarl (nreverse nvarl)))
 286                (and (not (zerop (car dl)))
 287                     (setq nvarl (cons (car dl) (cons (car varl) nvarl)))))
 288              (list '(%at) (cons '(%derivative) (cons fun nvarl)) eqs)))))
 289
 290 (defun dispatomgrads (l)
 291   (do ((i l (cdr i)))
 292       ((null i))
 293     (do ((j (mget (car i) '$atomgrad) (cdr j)))
 294         ((null j))
 295       (mtell-open "~M~%"
 296                   (list '(mlabel) nil
 297                         (list '(mequal)
 298                               (list '(%derivative) (car i) (caar j) 1) (cdar j))))))
 299   '$done)
 300
 301 (defun dispgradefs (l)
 302   (do ((i l (cdr i)))
 303       ((null i))
 304     (setq l (get (car i) 'grad))
 305     (do ((j (car l) (cdr j))
 306          (k (cdr l) (cdr k))
 307          (thing (cons (ncons (car i)) (car l))))
 308         ((or (null k) (null j)))
 309       (mtell-open "~M~%"
 310                   (list '(mlabel)
 311                         nil (list '(mequal) (list '(%derivative) thing (car j) 1.) (car k))))))
 312   '$done)
 313
 314 (defun dispmatchdeclares (l)
 315   (do ((i l (cdr i))
 316        (ret))
 317       ((null i) (cons '(mlist) (reverse ret)))
 318     (setq l (car (mget (car i) 'matchdeclare)))
 319     (setq ret (cons (append (cond ((atom l) (ncons (ncons l))) ((eq (caar l) 'lambda) (list '(mqapply) l))  (t l))
 320                             (ncons (car i)))
 321                     ret))))
 322
 323 (declare-top (special *roots *failures))
 324
 325 (defmfun $changevar (expr trans nvar ovar)
 326   (let ($ratfac)
 327     (cond ((or (atom expr) (eq (caar expr) 'rat) (eq (caar expr) 'mrat))
 328            expr)
 329           ((atom trans)
 330            (merror (intl:gettext "changevar: second argument must not be an atom; found: ~M") trans))
 331           ((null (atom nvar))
 332            (merror (intl:gettext "changevar: third argument must be an atom; found: ~M") nvar))
 333           ((null (atom ovar))
 334            (merror (intl:gettext "changevar: fourth argument must be an atom; found: ~M") ovar)))
 335     (changevar expr trans nvar ovar)))
 336
 337 (defun solvable (l var &optional (errswitch nil))
 338   (let (*roots *failures)
 339     (solve l var 1)
 340     (cond (*roots
 341            ;; We arbitrarily pick the first root.  Should we be more careful?
 342            ($rhs (car *roots)))
 343           (errswitch (merror (intl:gettext "changevar: failed to solve for ~M in ~M") var l))
 344           (t nil))))
 345
 346 (defun changevar (expr trans nvar ovar)
 347   (cond ((atom expr) expr)
 348         ((or (not (member (caar expr) '(%integrate %sum %product) :test #'eq))
 349              (not (alike1 (caddr expr) ovar)))
 350          (recur-apply (lambda (e) (changevar e trans nvar ovar)) expr))
 351         (t
 352          ;; TRANS is the expression that relates old var and new var
 353          ;; and is of the form f(ovar, nvar) = 0. Using TRANS, try to
 354          ;; solve for ovar so that ovar = tfun(nvar), if possible.
 355          (let* ((tfun (solvable (setq trans (meqhk trans)) ovar))
 356                 (deriv
 357                  ;; Compute diff(tfun, nvar) = dovar/dnvar if tfun is
 358                  ;; available.  Otherwise, use implicit
 359                  ;; differentiation.
 360                  (if tfun
 361                      (sdiff tfun nvar)
 362                      (neg (div (sdiff trans nvar) ;IMPLICIT DIFF.
 363                                (sdiff trans ovar)))))
 364                 (sum-product-p (member (caar expr) '(%sum %product) :test #'eq)))
 365
 366            #+nil
 367            (progn
 368              (mformat t "tfun = ~M~%" tfun)
 369              (mformat t "deriv = ~M~%" deriv))
 370
 371            ;; For sums and products, we want deriv to be +/-1 because
 372            ;; I think that means that integers will map into integers
 373            ;; (roughly), so that we don't need to express the
 374            ;; summation index or limits in some special way to account
 375            ;; for it.
 376            (when (and (member (caar expr) '(%sum %product) :test #'eq)
 377                       (not (or (equal deriv 1)
 378                                (equal deriv -1))))
 379              (merror (intl:gettext "changevar: illegal change in summation or product")))
 380
 381            (let ((nfun ($radcan ;NIL IF KERNSUBST FAILS
 382                               (if tfun
 383                                   (mul (maxima-substitute tfun ovar (cadr expr))
 384                                        ;; Don't multiply by deriv
 385                                        ;; for sums/products because
 386                                        ;; reversing the order of
 387                                        ;; limits doesn't change the
 388                                        ;; sign of the result.
 389                                        (if sum-product-p 1 deriv))
 390                                   (kernsubst ($ratsimp (mul (cadr expr)
 391                                                             deriv))
 392                                              trans ovar)))))
 393              (cond
 394                (nfun
 395                 ;; nfun is basically the result of substituting ovar
 396                 ;; with tfun in the integrand (summand).
 397                 (cond
 398                   ((cdddr expr)
 399                    ;; Handle definite integral, summation, or product.
 400                    ;; invfun expresses nvar in terms of ovar so that
 401                    ;; we can compute the new lower and upper limits of
 402                    ;; the integral (sum).
 403                    (let* ((invfun (solvable trans nvar t))
 404                           (lo-limit ($limit invfun ovar (cadddr expr) '$plus))
 405                           (hi-limit ($limit invfun
 406                                             ovar
 407                                             (car (cddddr expr))
 408                                             '$minus)))
 409                      ;; If this is a sum or product and deriv = -1, we
 410                      ;; want to reverse the low and high limits.
 411                      (when (and sum-product-p (equal deriv -1))
 412                        (rotatef lo-limit hi-limit))
 413
 414                      ;; Construct the new result.
 415                      (list (ncons (caar expr))
 416                            nfun
 417                            nvar
 418                            lo-limit
 419                            hi-limit)))
 420                   (t
 421                    ;; Indefinite integral
 422                    (list '(%integrate) nfun nvar))))
 423                (t expr)))))))
 424
 425 (defun kernsubst (expr form ovar)
 426   (let (varlist genvar nvarlist)
 427     (newvar expr)
 428     (setq nvarlist (mapcar #'(lambda (x) (if (freeof ovar x) x
 429                                              (solvable form x)))
 430                            varlist))
 431     (if (member nil nvarlist :test #'eq) nil
 432         (prog2 (setq expr (ratrep* expr)
 433                      varlist nvarlist)
 434             (rdis (cdr expr))))))
 435
 436 (declare-top (special facfun))
 437
 438 (defmfun $factorsum (e)
 439   (factorsum0 e '$factor))
 440
 441 (defmfun $gfactorsum (e)
 442   (factorsum0 e '$gfactor))
 443
 444 (defun factorsum0 (e facfun)
 445   (cond ((mplusp (setq e (funcall facfun e)))
 446          (factorsum1 (cdr e)))
 447         (t (factorsum2 e))))
 448
 449 (defun factorsum1 (e)
 450   (prog (f lv llv lex cl lt c)
 451    loop (setq f (car e))
 452    (setq lv (cdr ($showratvars f)))
 453    (cond ((null lv) (setq cl (cons f cl)) (go skip)))
 454    (do ((q llv (cdr q)) (r lex (cdr r)))
 455        ((null q))
 456      (cond ((intersect (car q) lv)
 457             (rplaca q (union* (car q) lv))
 458             (rplaca r (cons f (car r)))
 459             (return (setq lv nil)))))
 460    (or lv (go skip))
 461    (setq llv (cons lv llv) lex (cons (ncons f) lex))
 462    skip (and (setq e (cdr e)) (go loop))
 463    (or cl (go skip2))
 464    (do ((q llv (cdr q)) (r lex (cdr r)))
 465        ((null q))
 466      (cond ((and (null (cdar q)) (cdar r))
 467             (rplaca r (nconc cl (car r)))
 468             (return (setq cl nil)))))
 469    skip2 (setq llv nil lv nil)
 470    (do ((r lex (cdr r)))
 471        ((null r))
 472      (cond ((cdar r)
 473             (setq llv
 474                   (cons (factorsum2 (funcall facfun (cons '(mplus) (car r))))
 475                         llv)))
 476            ((or (not (mtimesp (setq f (caar r))))
 477                 (not (mnump (setq c (cadr f)))))
 478             (setq llv (cons f llv)))
 479            (t (do ((q lt (cdr q)) (s lv (cdr s)))
 480                   ((null q))
 481                 (cond ((alike1 (car s) c)
 482                        (rplaca q (cons (dcon f) (car q)))
 483                        (return (setq f nil)))))
 484               (and f
 485                    (setq lv (cons c lv)
 486                          lt (cons (ncons (dcon f)) lt))))))
 487    (setq lex
 488          (mapcar #'(lambda (s q)
 489                      (simptimes (list '(mtimes) s
 490                                       (cond ((cdr q)
 491                                              (cons '(mplus) q))
 492                                             (t (car q))))
 493                                 1 nil))
 494                  lv lt))
 495    (return (simplus (cons '(mplus) (nconc cl lex llv)) 1 nil))))
 496
 497 (defun dcon (mt)
 498   (cond ((cdddr mt) (cons (car mt) (cddr mt))) (t (caddr mt))))
 499
 500 (defun factorsum2 (e)
 501   (cond ((not (mtimesp e)) e)
 502         (t (cons '(mtimes)
 503                  (mapcar #'(lambda (f)
 504                              (cond ((mplusp f)
 505                                     (factorsum1 (cdr f)))
 506                                    (t f)))
 507                          (cdr e))))))
 508
 509 (declare-top (special $combineflag))
 510
 511 (defmvar $combineflag t)
 512
 513 (defmfun $combine (e)
 514   (cond ((or (atom e) (eq (caar e) 'rat)) e)
 515         ((eq (caar e) 'mplus) (combine (cdr e)))
 516         (t (recur-apply #'$combine e))))
 517
 518 (defun combine (e)
 519   (prog (term r ld sw nnu d ln xl)
 520    again(setq term (car e) e (cdr e))
 521    (when (or (not (or (ratnump term) (mtimesp term) (mexptp term)))
 522              (equal (setq d ($denom term)) 1))
 523      (setq r (cons term r))
 524      (go end))
 525    (setq nnu ($num term))
 526    (and $combineflag (integerp d) (setq xl (cons term xl)) (go end))
 527    (do ((q ld (cdr q)) (p ln (cdr p)))
 528        ((null q))
 529      (cond ((alike1 (car q) d)
 530             (rplaca p (cons nnu (car p)))
 531             (return (setq sw t)))))
 532    (and sw (go skip))
 533    (setq ld (cons d ld) ln (cons (ncons nnu) ln))
 534    skip (setq sw nil)
 535    end  (and e (go again))
 536    (and xl (setq xl (cond ((cdr xl) ($xthru (addn xl t)))
 537                           (t (car xl)))))
 538    (mapc
 539     #'(lambda (nu de)
 540         (setq r (cons (mul2 (addn nu nil) (power* de -1)) r)))
 541     ln ld)
 542    (return (addn (if xl (cons xl r) r) nil))))
 543
 544 (defmfun $factorout (e &rest vl)
 545   (prog (el fl cl l f x)
 546      (when (null vl)
 547        (merror (intl:gettext "factorout: at least two arguments required.")))
 548      (unless (mplusp e)
 549        (return e))
 550      (or (null vl) (mplusp e) (return e))
 551      (setq e (cdr e))
 552      loop       (setq f (car e) e (cdr e))
 553      (unless (mtimesp f)
 554        (setq f (list '(mtimes) 1 f)))
 555      (setq fl nil cl nil)
 556      (do ((i (cdr f) (cdr i)))
 557          ((null i))
 558        (if (and (not (numberp (car i)))
 559                 (apply '$freeof (append vl (ncons (car i)))))
 560            (setq fl (cons (car i) fl))
 561            (setq cl (cons (car i) cl))))
 562      (when (null fl)
 563        (push f el)
 564        (go end))
 565      (setq fl (if (cdr fl)
 566                   (simptimes (cons '(mtimes) fl) 1 nil)
 567                   (car fl)))
 568      (setq cl (cond ((null cl) 1)
 569                     ((cdr cl) (simptimes (cons '(mtimes) cl) 1 t))
 570                     (t (car cl))))
 571      (setq x t)
 572      (do ((i l (cdr i)))
 573          ((null i))
 574        (when (alike1 (caar i) fl)
 575          (rplacd (car i) (cons cl (cdar i)))
 576          (setq i nil x nil)))
 577      (when x
 578        (push (list fl cl) l))
 579      end  (when e (go loop))
 580      (do ((i l (cdr i)))
 581          ((null i))
 582        (push (simptimes (list '(mtimes) (caar i)
 583                               ($factorsum (simplus (cons '(mplus) (cdar i)) 1 nil))) 1 nil) el))
 584      (return (addn el nil))))