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