src/comm2.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 ;;      ** (c) Copyright 1982 Massachusetts Institute of Technology **
  11
  12 (macsyma-module comm2)
  13
  14 ;;;; DIFF2
  15
  16 (declare-top (special $props $dotdistrib))
  17
  18 (defun diffint (e x)
  19   (let (a)
  20     (cond ((null (cdddr e))
  21            (cond ((alike1 x (caddr e)) (cadr e))
  22                  ((and (not (atom (caddr e))) (atom x) (not (free (caddr e) x)))
  23                   (mul2 (cadr e) (sdiff (caddr e) x)))
  24                  ((or ($constantp (setq a (sdiff (cadr e) x)))
  25                       (and (atom (caddr e)) (free a (caddr e))))
  26                   (mul2 a (caddr e)))
  27                  (t (simplifya (list '(%integrate) a (caddr e)) t))))
  28           ((alike1 x (caddr e)) (addn (diffint1 (cdr e) x x) t))
  29           (t (addn (cons (if (equal (setq a (sdiff (cadr e) x)) 0)
  30                              0
  31                              (simplifya (list '(%integrate) a (caddr e)
  32                                               (cadddr e) (car (cddddr e)))
  33                                         t))
  34                          (diffint1 (cdr e) x (caddr e)))
  35                    t)))))
  36
  37 (defun diffint1 (e x y)
  38   (let ((u (sdiff (cadddr e) x)) (v (sdiff (caddr e) x)))
  39     (list (if (pzerop u) 0 (mul2 u (maxima-substitute (cadddr e) y (car e))))
  40           (if (pzerop v) 0 (mul3 v (maxima-substitute (caddr e) y (car e)) -1)))))
  41
  42 (defun diffsumprod (e x)
  43   (cond ((or (not ($mapatom x)) (not (free (cadddr e) x)) (not (free (car (cddddr e)) x)))
  44          (diff%deriv (list e x 1)))
  45         ((eq (caddr e) x) 0)
  46         (t (let ((u (sdiff (cadr e) x)))
  47              (setq u (simplifya (list '(%sum)
  48                                       (if (eq (caar e) '%sum) u (div u (cadr e)))
  49                                       (caddr e) (cadddr e) (car (cddddr e)))
  50                                 t))
  51              (if (eq (caar e) '%sum) u (mul2 e u))))))
  52
  53 (defun difflaplace (e x)
  54   (cond ((or (not (atom x)) (eq (cadddr e) x)) (diff%deriv (list e x 1)))
  55         ((eq (caddr e) x) 0)
  56         (t ($laplace (sdiff (cadr e) x) (caddr e) (cadddr e)))))
  57
  58 (defun diff-%at (e x)
  59   (cond ((freeof x e) 0)
  60         ((not (freeofl x (hand-side (caddr e) 'r))) (diff%deriv (list e x 1)))
  61         (t ($at (sdiff (cadr e) x) (caddr e)))))
  62
  63 (defun diffncexpt (e x)
  64   (let ((base* (cadr e))
  65         (pow (caddr e)))
  66     (cond ((and (mnump pow) (or (not (fixnump pow)) (< pow 0))) ; POW cannot be 0
  67            (diff%deriv (list e x 1)))
  68           ((and (atom base*) (eq base* x) (free pow base*))
  69            (mul2* pow (list '(mncexpt) base* (add2 pow -1))))
  70           ((fixnump pow)
  71            (let ((deriv (sdiff base* x))
  72                  (ans nil))
  73              (do ((i 0 (1+ i))) ((= i pow))
  74                (push (list '(mnctimes) (list '(mncexpt) base* i)
  75                            (list '(mnctimes) deriv
  76                                  (list '(mncexpt) base* (- pow 1 i))))
  77                      ans))
  78              (addn ans nil)))
  79           ((and (not (depends pow x)) (or (atom pow) (and (atom base*) (free pow base*))))
  80            (let ((deriv (sdiff base* x))
  81                  (index (gensumindex)))
  82              (simplifya
  83               (list '(%sum)
  84                     (list '(mnctimes) (list '(mncexpt) base* index)
  85                           (list '(mnctimes) deriv
  86                                 (list '(mncexpt) base*
  87                                       (list '(mplus) pow -1 (list '(mtimes) -1 index)))))
  88                     index 0 (list '(mplus) pow -1)) nil)))
  89           (t (diff%deriv (list e x 1))))))
  90
  91 (defun stotaldiff (e)
  92   (cond ((or (mnump e) (constant e)) 0)
  93         ((or (atom e) (member 'array (cdar e) :test #'eq))
  94          (let ((w (mget (if (atom e) e (caar e)) 'depends)))
  95            (if w (cons '(mplus)
  96                        (mapcar #'(lambda (x)
  97                                    (list '(mtimes) (chainrule e x) (list '(%del) x)))
  98                                w))
  99                (list '(%del) e))))
 100         ((specrepp e) (stotaldiff (specdisrep e)))
 101         ((eq (caar e) 'mnctimes)
 102          (let (($dotdistrib t))
 103            (add2 (ncmuln (cons (stotaldiff (cadr e)) (cddr e)) t)
 104                  (ncmul2 (cadr e) (stotaldiff (ncmuln (cddr e) t))))))
 105         ((eq (caar e) 'mncexpt)
 106          (if (and (fixnump (caddr e)) (> (caddr e) 0))
 107              (stotaldiff (list '(mnctimes) (cadr e)
 108                                (ncpower (cadr e) (1- (caddr e)))))
 109              (list '(%derivative) e)))
 110         (t (addn (cons 0 (mapcar #'(lambda (x)
 111                                      (mul2 (sdiff e x) (list '(%del simp) x)))
 112                                  (extractvars (margs e))))
 113                  t))))
 114
 115 (defun extractvars (e &aux vars)
 116   (cond ((null e) nil)
 117         ((atom (car e))
 118          (cond ((not (maxima-constantp (car e)))
 119                 (cond ((setq vars (mget (car e) 'depends))
 120                        ;; The symbol has dependencies. Put the dependencies on
 121                        ;; the list of extracted vars.
 122                        (union* vars (extractvars (cdr e))))
 123                       (t
 124                        ;; Put the symbol on the list of extracted vars.
 125                        (union* (ncons (car e)) (extractvars (cdr e))))))
 126              (t (extractvars (cdr e)))))
 127         ((member 'array (cdaar e) :test #'eq)
 128          (union* (ncons (car e)) (extractvars (cdr e))))
 129         (t (union* (extractvars (cdar e)) (extractvars (cdr e))))))
 130
 131 ;;;; AT
 132
 133 ;;dummy-variable-operators is defined in COMM, which uses it inside of SUBST1.
 134 (declare-top (special atvars *atp* munbound dummy-variable-operators))
 135
 136 (defmfun $atvalue (exp eqs val)
 137   (let (dl vl fun)
 138     (cond ((notloreq eqs) (improper-arg-err eqs '$atvalue))
 139           ((or (atom exp) (and (eq (caar exp) '%derivative) (atom (cadr exp))))
 140            (improper-arg-err exp '$atvalue)))
 141     (cond ((not (eq (caar exp) '%derivative))
 142            (setq fun (caar exp)
 143                  vl (cdr exp)
 144                  dl (make-list (length vl) :initial-element 0)))
 145           (t (setq fun (caaadr exp) vl (cdadr exp))
 146              (dolist (v vl)
 147                (setq dl (nconc dl (ncons (or (getf (cddr exp) v) 0)))))))
 148     (if (or (mopp fun) (eq fun 'mqapply)) (improper-arg-err exp '$atvalue))
 149     (atvarschk vl)
 150     (do ((vl1 vl (cdr vl1)) (l atvars (cdr l))) ((null vl1))
 151       (if (and (symbolp (car vl1)) (not (kindp (car vl1) '$constant)))
 152           (setq val (maxima-substitute (car l) (car vl1) val))
 153           (improper-arg-err (cons '(mlist) vl) '$atvalue)))
 154     (setq eqs (if (eq (caar eqs) 'mequal) (list eqs) (cdr eqs)))
 155     (setq eqs (do ((eqs eqs (cdr eqs)) (l)) ((null eqs) l)
 156                 (if (not (member (cadar eqs) vl :test #'eq))
 157                     (improper-arg-err (car eqs) '$atvalue))
 158                 (setq l (nconc l (ncons (cons (cadar eqs) (caddar eqs)))))))
 159     (setq vl (do ((vl vl (cdr vl)) (l)) ((null vl) l)
 160                (setq l (nconc l (ncons (cdr (or (assoc (car vl) eqs :test #'eq)
 161                                                 (cons nil munbound))))))))
 162     (do ((atvalues (mget fun 'atvalues) (cdr atvalues)))
 163         ((null atvalues)
 164          (mputprop fun (cons (list dl vl val) (mget fun 'atvalues)) 'atvalues))
 165       (when (and (equal (caar atvalues) dl) (equal (cadar atvalues) vl))
 166         (rplaca (cddar atvalues) val) (return nil)))
 167     (add2lnc fun $props)
 168     val))
 169
 170 (defprop %at simp-%at operators)
 171
 172 (defun simp-%at (expr ignored simp-flag)
 173   (declare (ignore ignored))
 174   (twoargcheck expr)
 175   (let* ((arg (simpcheck (cadr expr) simp-flag))
 176          (e (resimplify (caddr expr)))
 177          (eqn (if ($listp e)
 178                   (if (= ($length e) 1) ($first e) (cons '(mlist simp) (cdr ($sort e))))
 179                   e)))
 180     (cond (($constantp arg) arg)
 181           ((alike1 eqn '((mlist))) arg)
 182           ((at-not-dependent eqn arg))
 183           (t (eqtest (list '(%at) arg eqn) expr)))))
 184
 185 ;; Remove any variable from EQN if ARG is not dependent on it.
 186 (defun at-not-dependent (eqn arg)
 187   (if (eq (caar eqn) 'mequal)
 188     (setq eqn (list '(mlist) eqn)))
 189   (multiple-value-bind (e0 e1) (at-not-dependent-find-vars eqn arg)
 190     (if e0
 191       (if e1
 192         (let*
 193           ((e1 (mapcar #'(lambda (x) (list '(mequal) x ($assoc x eqn))) e1))
 194            (eqn1 (if (= (length e1) 1) (first e1) (cons '(mlist) e1))))
 195           (list '(%at) arg eqn1))
 196         arg))))
 197
 198 ;; Test dependence via derivative to account for declared dependencies.
 199 (defun at-not-dependent-find-vars (eqn arg)
 200   (let ((e (mapcar #'second (rest eqn))))
 201     (partition-by #'(lambda (x) (at-not-dependent-find-vars-1 x arg)) e)))
 202
 203 (defun at-not-dependent-find-vars-1 (x arg)
 204   (if ($mapatom x)
 205     (eql (mfuncall '$diff arg x) 0)
 206     ;; We might be called with something like -1*x as the variable.
 207     ;; (That might or might not be a bug in itself, but let it go for the moment.)
 208     ;; Try to extract a variable and test for dependence on that.
 209     ;; If there are 2 or more variables, return NIL (i.e., not at-not-dependent).
 210     (let ((v ($listofvars x)))
 211       (if (eql ($length v) 1)
 212         (at-not-dependent-find-vars-1 ($first v) arg)))))
 213
 214 (defmfun $at (expr ateqs)
 215   (if (notloreq ateqs) (improper-arg-err ateqs '$at))
 216   (atscan (let ((*atp* t)) ($psubstitute ateqs expr)) ateqs))
 217
 218 (defun atscan (expr ateqs)
 219   (cond ((or (atom expr)
 220              (eq (caar expr) 'mrat)
 221              (like ateqs '((mlist))))
 222          expr)
 223         ((eq (caar expr) '%derivative)
 224          (or (and (not (atom (cadr expr)))
 225                   (let ((vl (cdadr expr)) dl)
 226                     (dolist (v vl)
 227                       (setq dl (nconc dl (ncons (or (getf (cddr expr) v) 0)))))
 228                     (atfind (caaadr expr)
 229                             (cdr ($psubstitute ateqs (cons '(mlist) vl)))
 230                             dl)))
 231              (list '(%at) expr ateqs)))
 232         ((member (caar expr) dummy-variable-operators :test #'eq)
 233          (list '(%at) expr ateqs))
 234         ((at1 expr))
 235         (t (recur-apply #'(lambda (x) (atscan x ateqs)) expr))))
 236
 237 (defun at1 (expr)
 238   (atfind (caar expr) (cdr expr) (make-list (length (cdr expr)) :initial-element 0)))
 239
 240 (defun atfind (fun vl dl)
 241   (do ((atvalues (mget fun 'atvalues) (cdr atvalues)))
 242       ((null atvalues))
 243     (and (equal (caar atvalues) dl)
 244          (do ((l (cadar atvalues) (cdr l)) (vl vl (cdr vl)))
 245              ((null l) t)
 246            (if (and (not (equal (car l) (car vl)))
 247                     (not (eq (car l) munbound)))
 248                (return nil)))
 249          (return (prog2
 250                     (atvarschk vl)
 251                     (substitutel vl atvars (caddar atvalues)))))))
 252
 253 (declare-top (special $ratfac genvar varlist $keepfloat))
 254
 255 (defmvar $logconcoeffp nil)
 256
 257 (defmfun $logcontract (e)
 258   (lgcreciprocal (logcon e))) ; E is assumed to be simplified.
 259
 260 (defun logcon (e)
 261   (cond ((atom e) e)
 262         ((member (caar e) '(mplus mtimes) :test #'eq)
 263          (if (not (lgcsimplep e)) (setq e (lgcsort e)))
 264          (cond ((mplusp e) (lgcplus e))
 265                ((mtimesp e) (lgctimes e))
 266                (t (logcon e))))
 267         (t (recur-apply #'logcon e))))
 268
 269 ;; The logcontract algorithm for a sum.
 270 ;;
 271 ;; The function accumulates the arguments of things like log(a)+log(b) into a
 272 ;; list called LOG. It calls out to lgctimes to deal with things like
 273 ;; a*log(b). When all the arguments have been processed, it simplifies all the
 274 ;; logarithmic arguments using sratsimp.
 275 (defun lgcplus (e)
 276   (let ((log) (notlogs))
 277     (dolist (arg (cdr e))
 278       (cond
 279         ((atom arg) (push arg notlogs))
 280         ;; Only gather up log(x), not log[x]. It's not particularly obvious
 281         ;; whether log(x)+log[y] should become log(x*y) or log[x*y], so we just
 282         ;; ignore the fact that log[x] is a logarithm.
 283         ((and (eq (caar arg) '%log)
 284               (not (member 'array (car arg))))
 285          (push (logcon (second arg)) log))
 286         ((eq (caar arg) 'mtimes)
 287          (let ((y (lgctimes arg)))
 288            (if (or (atom y) (not (eq (caar y) '%log)))
 289                (push y notlogs)
 290                (push (cadr y) log))))
 291         (t
 292          (push (logcon arg) notlogs))))
 293     (cond
 294       ((null log)
 295        (subst0 (cons '(mplus) (nreverse notlogs)) e))
 296       (t
 297        (let ((simplified-log (lgcsimp
 298                               (let (($ratfac t))
 299                                 (sratsimp (muln log t))))))
 300          (addn (cons simplified-log notlogs) t))))))
 301
 302 ;; The logcontract algorithm for a product
 303 ;;
 304 ;; The main transformation this does is of the form 3*log(x) => log(x^3). To
 305 ;; make this work, we find the first %log term and insert any coefficients we
 306 ;; find into that. Coefficients are identified by LOGCONCOEFFP, which checks the
 307 ;; $LOGCONCOEFFP user variable.
 308 (defun lgctimes (e)
 309   ;; Apply logcontract to the arguments. It's possible that the subsequent
 310   ;; simplification means that the result isn't a product any more. In that
 311   ;; case, just return it.
 312   (setq e (subst0 (cons '(mtimes) (mapcar 'logcon (cdr e))) e))
 313   (if (not (mtimesp e))
 314       e
 315       (let ((log) (notlogs) (decints))
 316         (dolist (arg (cdr e))
 317           (cond ((and (null log) (not (atom arg))
 318                       (eq (caar arg) '%log) (not (equal (cadr arg) -1)))
 319                  (setq log (cadr arg)))
 320                 ((logconcoeffp arg) (push arg decints))
 321                 (t (setq notlogs (push arg notlogs)))))
 322         (cond
 323           ((or (null log) (null decints)) e)
 324           (t (muln (cons (lgcsimp (power log (muln decints t)))
 325                          notlogs)
 326                    t))))))
 327
 328 (defun lgcsimp (e)
 329   (cond ((atom e)
 330          ;; e.g. log(1) -> 0, or log(%e) -> 1
 331          (simplify (list '(%log) e)))
 332         ((and (mexptp e) (eq (cadr e) '$%e))
 333          ;; log(%e^expr) -> expr
 334          (simplify (list '(%log) e)))
 335         (t
 336          (list '(%log simp) e))))
 337
 338 ;; Tests that its argument is a sum of terms that are "simple".
 339 ;;
 340 ;; A "simple" term is either completely free of logarithms, is a logarithm
 341 ;; itself, or is a number times a logarithm.
 342 ;;
 343 ;; This function assumes that its argument is not an atom.
 344 (defun lgcsimplep (e)
 345   (flet ((lgc-nonsimple-arg-p (arg)
 346            (not (or (atom arg)
 347                     (eq (caar arg) '%log)
 348                     (not (isinop arg '%log))
 349                     ;; Product of a number with a logarithm e.g. 3*log(x)
 350                     (and (eq (caar arg) 'mtimes)
 351                          (null (cdddr arg))
 352                          (mnump (cadr arg))
 353                          (not (atom (caddr arg)))
 354                          (eq (caar (caddr arg)) '%log))))))
 355     (and (eq (caar e) 'mplus)
 356          (not (find-if #'lgc-nonsimple-arg-p (cdr e))))))
 357
 358 ;; Sort the argument so that coefficients come before logarithms and logarithms
 359 ;; come before everything else.
 360 (defun lgcsort (e)
 361   (let ((logs) (notlogs) (decints) (varlist))
 362     ;; Split the variables in E into logs, notlogs and coefficients. The list of
 363     ;; variables is calculated by NEWVAR (and stored in the special variable
 364     ;; VARLIST, which is why we have to bind it above).
 365     (dolist (var (newvar e))
 366       (cond
 367         ((and (not (atom var)) (eq (caar var) '%log)) (push var logs))
 368         ((logconcoeffp var) (push var decints))
 369         (t (push var notlogs))))
 370     (let* ((vl (nreconc decints (nconc (sort logs #'great)
 371                                        (nreverse notlogs))))
 372            (e1 (ratdisrep (ratrep e vl))))
 373       (if (alike1 e e1) e e1))))
 374
 375 ;; lgcreciprocal performs the transformation log(1/x) => -log(x)
 376 (defun lgcreciprocal (e)
 377   (let (num denom)
 378     (cond
 379       ((atom e) e)
 380       ((and (eq (caar e) '%log)
 381             (setq num (member ($num (cadr e)) '(1 -1) :test #'equal))
 382             (not (equal (setq denom ($denom (cadr e))) 1)))
 383        (list '(mtimes simp) -1
 384              (list '(%log simp) (if (= (car num) 1) denom (neg denom)))))
 385       (t (recur-apply #'lgcreciprocal e)))))
 386
 387 (defun logconcoeffp (e)
 388   (if $logconcoeffp
 389       (is `(($logconcoeffp) ,e))
 390       (maxima-integerp e)))
 391
 392 ;;;; RTCON
 393
 394 (declare-top (special $radexpand $domain))
 395
 396 (defmvar $rootsconmode t)
 397
 398 (defmfun $rootscontract (e)            ; E is assumed to be simplified
 399   (let ((radpe (and $radexpand (not (eq $radexpand '$all)) (eq $domain '$real)))
 400         ($radexpand nil))
 401     (rtcon e radpe)))
 402
 403 (defun rtcon (e radpe)
 404   (cond ((atom e) e)
 405         ((eq (caar e) 'mtimes)
 406          (do ((x (cdr e) (cdr x)) (roots) (notroots) (y))
 407              ((null x)
 408               (cond ((null roots) (subst0 (cons '(mtimes) (nreverse notroots)) e))
 409                     (t (if $rootsconmode
 410                            (multiple-value-bind (min gcd lcm)
 411                                (rtc-getinfo roots)
 412                              (cond ((and (= min gcd) (not (= gcd 1))
 413                                        (not (= min lcm))
 414                                        (not (eq $rootsconmode '$all)))
 415                                     (setq roots
 416                                           (rt-separ
 417                                            (list gcd
 418                                                  (rtcon
 419                                                   (rtc-fixitup
 420                                                    (rtc-divide-by-gcd roots gcd)
 421                                                    nil) radpe)
 422                                                  1)
 423                                            nil)))
 424                                    ((eq $rootsconmode '$all)
 425                                     (setq roots
 426                                           (rt-separ (simp-roots lcm roots)
 427                                                     nil))))))
 428                        (rtc-fixitup roots notroots))))
 429            (cond ((atom (car x))
 430                   (cond ((eq (car x) '$%i) (setq roots (rt-separ (list 2 -1) roots)))
 431                         (t (setq notroots (cons (car x) notroots)))))
 432                  ((and (eq (caaar x) 'mexpt) (ratnump (setq y (caddar x))))
 433                   (setq roots (rt-separ (list (caddr y)
 434                                               (list '(mexpt)
 435                                                     (rtcon (cadar x) radpe) (cadr y)))
 436                                         roots)))
 437
 438                  ((and radpe (eq (caaar x) 'mabs))
 439                   (setq roots (rt-separ (list 2 `((mexpt) ,(rtcon (cadar x) radpe) 2) 1)
 440                                         roots)))
 441                  (t (setq notroots (cons (rtcon (car x) radpe) notroots))))))
 442         ((and radpe (eq (caar e) 'mabs))
 443          (power (power (rtcon (cadr e) radpe) 2) '((rat simp) 1 2)))
 444         (t (recur-apply #'(lambda (x) (rtcon x radpe)) e))))
 445
 446 ;; RT-SEPAR separates like roots into their appropriate "buckets",
 447 ;; where a bucket looks like:
 448 ;; ((<denom of power> (<term to be raised> <numer of power>)
 449 ;;                   (<term> <numer>)) etc)
 450
 451 (defun rt-separ (a roots)
 452   (let ((u (assoc (car a) roots :test #'equal)))
 453     (cond (u (nconc u (cdr a))) (t (setq roots (cons a roots)))))
 454   roots)
 455
 456 (defun simp-roots (lcm root-list)
 457   (let (root1)
 458     (do ((x root-list (cdr x)))
 459         ((null x) (push lcm root1))
 460       (push (list '(mexpt) (muln (cdar x) nil) (quotient lcm (caar x)))
 461             root1))))
 462
 463 (defun rtc-getinfo (list)
 464   (let ((m (caar list))
 465         (g (caar list))
 466         (l (caar list)))
 467     (dolist (x (cdr list) (values m g l))
 468       (setq m (min m (car x))
 469             g (gcd g (car x))
 470             l (lcm l (car x))))))
 471
 472 (defun rtc-fixitup (roots notroots)
 473   (mapcar #'(lambda (x) (rplacd x (list (sratsimp (muln (cdr x) (not $rootsconmode))))))
 474           roots)
 475   (muln (nconc (mapcar #'(lambda (x) (power* (cadr x) `((rat) 1 ,(car x))))
 476                        roots)
 477                notroots)
 478         (not $rootsconmode)))
 479
 480 (defun rtc-divide-by-gcd (llist gcd)
 481   (mapcar #'(lambda (x) (rplaca x (quotient (car x) gcd))) llist)
 482   llist)
 483
 484 (defmfun $nterms (e)
 485   (cond ((zerop1 e) 0)
 486         ((atom e) 1)
 487         ((eq (caar e) 'mtimes)
 488          (if (equal -1 (cadr e)) (setq e (cdr e)))
 489          (do ((l (cdr e) (cdr l)) (c 1 (* c ($nterms (car l)))))
 490              ((null l) c)))
 491         ((eq (caar e) 'mplus)
 492          (do ((l (cdr e) (cdr l)) (c 0 (+ c ($nterms (car l)))))
 493              ((null l) c)))
 494         ((and (eq (caar e) 'mexpt) (integerp (caddr e)) (plusp (caddr e)))
 495          ($binomial (+ (caddr e) ($nterms (cadr e)) -1) (caddr e)))
 496         ((specrepp e) ($nterms (specdisrep e)))
 497         (t 1)))
 498
 499 ;;;; ATAN2
 500
 501 (declare-top (special $numer $logarc $trigsign))
 502
 503 ;; atan2 distributes over lists, matrices, and equations
 504 (defprop $atan2 (mlist $matrix mequal) distribute_over)
 505
 506 (defun simpatan2 (expr vestigial z)     ; atan2(y,x) ~ atan(y/x)
 507   (declare (ignore vestigial))
 508   (twoargcheck expr)
 509   (let (y x signy signx)
 510     (setq y (simpcheck (cadr expr) z)
 511           x (simpcheck (caddr expr) z))
 512     (cond ((and (zerop1 y) (zerop1 x))
 513            (merror (intl:gettext "atan2: atan2(0,0) is undefined.")))
 514           ( ;; float contagion
 515            (and (or (numberp x) (ratnump x))       ; both numbers
 516                 (or (numberp y) (ratnump y))       ; ... but not bigfloats
 517                 (or $numer (floatp x) (floatp y))) ; at least one float
 518            (atan ($float y) ($float x)))
 519           ( ;; bfloat contagion
 520            (and (mnump x)
 521                 (mnump y)
 522                 (or ($bfloatp x) ($bfloatp y)))    ; at least one bfloat
 523            (setq x ($bfloat x)
 524                  y ($bfloat y))
 525            (*fpatan y (list x)))
 526           ;; Simplifify infinities
 527           ((or (eq x '$inf)
 528                (alike1 x '((mtimes) -1 $minf)))
 529            ;; Simplify atan2(y,inf) -> 0
 530            0)
 531           ((or (eq x '$minf)
 532                (alike1 x '((mtimes) -1 $inf)))
 533            ;; Simplify atan2(y,minf) -> %pi for realpart(y)>=0 or
 534            ;; -%pi for realpart(y)<0. When sign of y unknwon, return noun form.
 535            (cond ((member (setq signy ($sign ($realpart x))) '($pos $pz $zero))
 536                   '$%pi)
 537                  ((eq signy '$neg) (mul -1 '$%pi))
 538                  (t (eqtest (list '($atan2) y x) expr))))
 539           ((or (eq y '$inf)
 540                (alike1 y '((mtimes) -1 $minf)))
 541            ;; Simplify atan2(inf,x) -> %pi/2
 542            (div '$%pi 2))
 543           ((or (eq y '$minf)
 544                (alike1 y '((mtimes -1 $inf))))
 545            ;; Simplify atan2(minf,x) -> -%pi/2
 546            (div '$%pi -2))
 547           ((and (free x '$%i) (setq signx ($sign x))
 548                 (free y '$%i) (setq signy ($sign y))
 549                 (cond ((zerop1 y)
 550                        (cond ((eq signx '$neg) '$%pi)
 551                              ((member signx '($pos $pz)) 0)))
 552                       ((zerop1 x)
 553                        (cond ((eq signy '$neg) (div '$%pi -2))
 554                              ((member signy '($pos $pz)) (div '$%pi 2))))
 555                       ((alike1 y x)
 556                        (cond ((eq signx '$neg) (mul -3 (div '$%pi 4)))
 557                              ((member signx '($pos $pz)) (div '$%pi 4))))
 558                       ((alike1 y (mul -1 x))
 559                        (cond ((eq signx '$neg) (mul 3 (div '$%pi 4)))
 560                              ((member signx '($pos $pz)) (div '$%pi -4)))))))
 561           ($logarc
 562            (logarc '%atan2 (list ($logarc y) ($logarc x))))
 563           ((and $trigsign (eq t (mminusp y)))
 564            (neg (take '($atan2) (neg y) x)))
 565           ;; atan2(y,x) = atan(y/x) + pi sign(y) (1-sign(x))/2
 566           ((eq signx '$pos)
 567            (take '(%atan) (div y x)))
 568           ((and (eq signx '$neg)
 569                 (member (setq signy ($csign y)) '($pos $neg) :test #'eq))
 570            (add (take '(%atan) (div y x))
 571                 (porm (eq signy '$pos) '$%pi)))
 572           ((and (eq signx '$zero) (eq signy '$zero))
 573            ;; Unfortunately, we'll rarely get here.  For example,
 574            ;; assume(equal(x,0)) atan2(x,x) simplifies via the alike1 case above
 575            (merror (intl:gettext "atan2: atan2(0,0) is undefined.")))
 576           (t (eqtest (list '($atan2) y x) expr)))))
 577
 578 ;;;; ARITHF
 579
 580 (defmfun $fibtophi (e &optional (lnorecurse nil))
 581   (cond ((atom e) e)
 582         ((eq (caar e) '$fib)
 583          (setq e (cond (lnorecurse (cadr e)) (t ($fibtophi (cadr e) lnorecurse))))
 584          (let ((phi (meval '$%phi)))
 585            (div (add2 (power phi e) (neg (power (add2 1 (neg phi)) e)))
 586                 (add2 -1 (mul2 2 phi)))))
 587         (t (recur-apply #'(lambda (x) ($fibtophi x lnorecurse)) e))))
 588
 589 (defmspec $numerval (l) (setq l (cdr l))
 590           (do ((l l (cddr l)) (x (ncons '(mlist simp)))) ((null l) x)
 591             (cond ((null (cdr l)) (merror (intl:gettext "numerval: expected an even number of arguments.")))
 592                   ((not (symbolp (car l)))
 593                    (merror (intl:gettext "numerval: expected a symbol; found ~M") (car l)))
 594                   ((boundp (car l))
 595                    (merror (intl:gettext "numerval: cannot declare a value because ~M is bound.") (car l))))
 596             (mputprop (car l) (cadr l) '$numer)
 597             (add2lnc (car l) $props)
 598             (nconc x (ncons (car l)))))
 599
 600 (let (my-powers)
 601   (declare (special my-powers))
 602
 603   (defmfun $derivdegree (e depvar var)
 604     (let (my-powers) (declare (special my-powers)) (derivdeg1 e depvar var) (if (null my-powers) 0 (maximin my-powers '$max))))
 605
 606   (defun derivdeg1 (e depvar var)
 607     (cond ((or (atom e) (specrepp e)))
 608           ((eq (caar e) '%derivative)
 609            (cond ((alike1 (cadr e) depvar)
 610                   (do ((l (cddr e) (cddr l))) ((null l))
 611                     (cond ((alike1 (car l) var)
 612                            (return (setq my-powers (cons (cadr l) my-powers)))))))))
 613           (t (mapc #'(lambda (x) (derivdeg1 x depvar var)) (cdr e))))))
 614
 615 ;;;; BOX
 616
 617 ;; Set the the property reversealias
 618 (defprop mbox $box reversealias)
 619 (defprop mlabox $box reversealias)
 620
 621 (defmfun $dpart (&rest args)
 622   (mpart args nil t nil '$dpart))
 623
 624 (defmfun $lpart (e &rest args)
 625   (mpart args nil (list e) nil '$lpart))
 626
 627 (defmfun $box (e &optional (l nil l?))
 628   (if l?
 629       (list '(mlabox) e (box-label l))
 630       (list '(mbox) e)))
 631
 632 (defun box (e label)
 633   (if (eq label t)
 634       (list '(mbox) e)
 635       ($box e (car label))))
 636
 637 (defun box-label (x)
 638   (if (atom x)
 639       x
 640       (coerce (mstring x) 'string)))
 641
 642 (defmfun $rembox (e &optional (l nil l?))
 643   (let ((label (if l? (box-label l) '(nil))))
 644     (rembox1 e label)))
 645
 646 (defun rembox1 (e label)
 647   (cond ((atom e) e)
 648         ((or (and (eq (caar e) 'mbox)
 649                   (or (equal label '(nil)) (member label '($unlabelled $unlabeled) :test #'eq)))
 650              (and (eq (caar e) 'mlabox)
 651                   (or (equal label '(nil)) (equal label (caddr e)))))
 652          (rembox1 (cadr e) label))
 653         (t (recur-apply #'(lambda (x) (rembox1 x label)) e))))
 654
 655 ;;;; MAPF
 656
 657 (declare-top (special scanmapp))
 658
 659 (defmspec $scanmap (l)
 660   (let ((scanmapp t))
 661     (resimplify (apply #'scanmap1 (mmapev l)))))
 662
 663 (defun scanmap1 (func e &optional (flag nil flag?))
 664   (let ((arg2 (specrepcheck e)) newarg2)
 665     (cond ((eq func '$rat)
 666            (merror (intl:gettext "scanmap: cannot apply 'rat'.")))
 667           (flag?
 668            (unless (eq flag '$bottomup)
 669              (merror (intl:gettext "scanmap: third argument must be 'bottomup', if present; found ~M") flag))
 670            (if (mapatom arg2)
 671                (funcer func (ncons arg2))
 672                (subst0 (funcer func
 673                                (ncons (mcons-op-args (mop arg2)
 674                                                      (mapcar #'(lambda (u)
 675                                                                  (scanmap1 func u '$bottomup))
 676                                                              (margs arg2)))))
 677                        arg2)))
 678           ((mapatom arg2)
 679            (funcer func (ncons arg2)))
 680           (t
 681            (setq newarg2 (specrepcheck (funcer func (ncons arg2))))
 682              (cond ((mapatom newarg2)
 683                     newarg2)
 684                    ((and (alike1 (cadr newarg2) arg2) (null (cddr newarg2)))
 685                     (subst0 (cons (ncons (caar newarg2))
 686                                   (ncons (subst0
 687                                           (mcons-op-args (mop arg2)
 688                                                          (mapcar #'(lambda (u) (scanmap1 func u))
 689                                                                  (margs arg2)))
 690                                           arg2)))
 691                             newarg2))
 692                    (t
 693                     (subst0 (mcons-op-args (mop newarg2)
 694                                            (mapcar #'(lambda (u) (scanmap1 func u))
 695                                                    (margs newarg2)))
 696                             newarg2)))))))
 697
 698 (defun subgen (form)       ; This function does mapping of subscripts.
 699   (do ((ds (if (eq (caar form) 'mqapply) (list (car form) (cadr form))
 700                (ncons (car form)))
 701            (outermap1 #'dsfunc1 (simplify (car sub)) ds))
 702        (sub (reverse (or (and (eq 'mqapply (caar form)) (cddr form))
 703                          (cdr form)))
 704             (cdr sub)))
 705       ((null sub) ds)))
 706
 707 (defun dsfunc1 (dsn dso)
 708   (cond ((or (atom dso) (atom (car dso))) dso)
 709         ((member 'array (car dso) :test #'eq)
 710          (cond ((eq 'mqapply (caar dso))
 711                 (nconc (list (car dso) (cadr dso) dsn) (cddr dso)))
 712                (t (nconc (list (car dso) dsn) (cdr dso)))))
 713         (t (mapcar #'(lambda (d) (dsfunc1 dsn d)) dso))))
 714
 715 ;;;; GENMAT
 716
 717 (defmfun $genmatrix (a i2 &optional (j2 i2) (i1 1) (j1 i1))
 718   (let ((f) (l (ncons '($matrix))))
 719     (setq f (if (or (symbolp a) (hash-table-p a) (arrayp a))
 720                 #'(lambda (i j) (meval (list (list a 'array) i j)))
 721               #'(lambda (i j) (mfuncall a i j))))
 722
 723     (if (notevery #'fixnump (list i2 j2 i1 j1))
 724       (merror (intl:gettext "genmatrix: bounds must be integers; found ~M, ~M, ~M, ~M") i2 j2 i1 j1))
 725
 726     (if (or (> i1 i2) (> j1 j2))
 727       (merror (intl:gettext "genmatrix: upper bounds must be greater than or equal to lower bounds; found ~M, ~M, ~M, ~M") i2 j2 i1 j1))
 728
 729     (dotimes (i (1+ (- i2 i1)))
 730       (nconc l (ncons (ncons '(mlist)))))
 731     (do ((i i1 (1+ i))
 732          (l (cdr l) (cdr l)))
 733         ((> i i2))
 734       (do ((j j1 (1+ j)))
 735           ((> j j2))
 736         (nconc (car l) (ncons (funcall f i j)))))
 737     l))
 738
 739 ; Execute deep copy for copymatrix and copylist.
 740 ; Resolves SF bug report [ 1224960 ] sideeffect with copylist.
 741 ; An optimization would be to call COPY-TREE only on mutable expressions.
 742
 743 (defmfun $copymatrix (x)
 744   (unless ($matrixp x)
 745     (merror (intl:gettext "copymatrix: argument must be a matrix; found ~M") x))
 746   (copy-tree x))
 747
 748 (defmfun $copylist (x)
 749   (unless ($listp x)
 750     (merror (intl:gettext "copylist: argument must be a list; found ~M") x))
 751   (copy-tree x))
 752
 753 (defmfun $copy (x)
 754   (copy-tree x))
 755
 756 ;;;; ADDROW
 757
 758 (defmfun $addrow (m &rest rows)
 759   (declare (dynamic-extent rows))
 760   (cond ((not ($matrixp m))
 761          (merror
 762            (intl:gettext "addrow: first argument must be a matrix; found ~M")
 763           m))
 764         ((null rows) m)
 765         (t
 766          (let ((m (copy-tree m)))
 767            (dolist (r rows m)
 768              (setq m (addrow m r)))))))
 769
 770 (defmfun $addcol (m &rest cols)
 771   (declare (dynamic-extent cols))
 772   (cond ((not ($matrixp m)) (merror (intl:gettext "addcol: first argument must be a matrix; found ~M") m))
 773         ((null cols) m)
 774         ((null (cdr m))
 775          (apply '$addcol (cons (ensure-matrix-column (first cols)) (rest cols))))
 776         (t (let ((m ($transpose m)))
 777              (dolist (c cols ($transpose m))
 778                (setq m (addrow m ($transpose c))))))))
 779
 780 (defun ensure-matrix-column (a)
 781   (if ($matrixp a) a
 782     ;; otherwise must be a MLIST.
 783     `(($matrix) ,@(mapcar #'(lambda (e) `((mlist) ,e)) (cdr a)))))
 784
 785 (defun addrow (m r)
 786   (cond ((not (mxorlistp r)) (merror (intl:gettext "addrow or addcol: argument must be a matrix or list; found ~M") r))
 787         ((and (cdr m)
 788               (or (and (eq (caar r) 'mlist) (not (= (length (cadr m)) (length r))))
 789                   (and (eq (caar r) '$matrix)
 790                        (not (= (length (cadr m)) (length (cadr r))))
 791                        (prog2 (setq r ($transpose r))
 792                            (not (= (length (cadr m)) (length (cadr r))))))))
 793          (merror (intl:gettext "addrow or addcol: incompatible structure."))))
 794   (append m (if (eq (caar r) '$matrix) (cdr r) (ncons r))))
 795
 796 ;;;; ARRAYF
 797
 798 (defun my-nonatomic-expr-p (e)
 799   (and (consp e) (consp (car e)) (symbolp (caar e))))
 800
 801 (defun my-lambda-expr-p (e)
 802   (and (consp e) (consp (car e)) (eq 'lambda (caar e))))
 803
 804 (defmfun $arraymake (ary subs)
 805   (cond
 806     ;; We go through some gyrations here to allow as wide a range of inputs as possible.
 807     ;; Previously $ARRAYMAKE didn't check the first argument at all;
 808     ;; this is an attempt at a minimally-restrictive change.
 809         ((not (or (symbolp ary) ($subvarp ary) (and (my-nonatomic-expr-p ary) (not (my-lambda-expr-p ary)))))
 810          (merror (intl:gettext "arraymake: first argument must be a symbol, subscripted symbol, or nonatomic expression (but not a lambda expression); found: ~M") ary))
 811         ((or (not ($listp subs)) (null (cdr subs)))
 812          (merror (intl:gettext "arraymake: second argument must be a list of one or more elements; found ~M") subs))
 813         ((symbolp ary)
 814          (cons (cons (getopr ary) '(array)) (cdr subs)))
 815         (t (cons '(mqapply array) (cons ary (cdr subs))))))
 816
 817 (defmspec $arrayinfo (ary)
 818   (setq ary (cdr ary))
 819   (arrayinfo-aux (car ary) (getvalue (car ary))))
 820
 821 (defun arrayinfo-aux (sym val)
 822   (prog (arra ary)
 823      (setq arra val)
 824      (setq ary sym)
 825      (if (and arra
 826               (or (hash-table-p arra)
 827                   (arrayp arra)
 828                   (eq (marray-type arra) '$functional)))
 829          (cond ((hash-table-p arra)
 830                 (let ((dim1 (gethash 'dim1 arra)))
 831                   (return (list* '(mlist) '$hash_table (if dim1 1 t)
 832                                  (loop for u being the hash-keys in arra
 833                                     unless (eq u 'dim1)
 834                                     collect
 835                                     (if dim1
 836                                         u
 837                                         (cons '(mlist simp) u)))))))
 838                ((arrayp arra)
 839                 (return (let ((dims (array-dimensions arra)))
 840                           (list '(mlist) '$declared
 841                                 ;; they don't want more info (array-type arra)
 842                                 (length dims)
 843                                 (cons '(mlist) (mapcar #'1- dims))))))
 844                ((eq (marray-type arra) '$functional)
 845                 (return (arrayinfo-aux sym (mgenarray-content arra)))))
 846          (let ((gen (safe-mgetl sym '(hashar array))) ary1)
 847            (when (null gen)
 848              (merror (intl:gettext "arrayinfo: ~M is not an array.") ary))
 849            (setq ary1 (cadr gen))
 850            (cond ((eq (car gen) 'hashar)
 851                   (setq ary1 (symbol-array ary1))
 852                   (return (append '((mlist simp) $hashed)
 853                                   (cons (aref ary1 2)
 854                                         (do ((i 3 (1+ i)) (l)
 855                                              (n (cadr (arraydims ary1))))
 856                                             ((= i n) (sort l #'(lambda (x y) (great y x))))
 857                                           (do ((l1 (aref ary1 i) (cdr l1)))
 858                                               ((null l1))
 859                                             (push (cons '(mlist simp) (caar l1)) l)))))))
 860                  (t (setq ary1 (arraydims ary1))
 861                     (return (list '(mlist simp)
 862                                   (cond ((safe-get ary 'array)
 863                                          (cdr (assoc (car ary1)
 864                                                      '((t . $complete) (fixnum . $integer)
 865                                                        (flonum . $float)) :test #'eq)))
 866                                         (t '$declared))
 867                                   (length (cdr ary1))
 868                                   (cons '(mlist simp) (mapcar #'1- (cdr ary1)))))))))))
 869
 870 ;;;; ALIAS
 871
 872 (declare-top (special greatorder lessorder))
 873
 874 (defmspec $ordergreat (l)
 875   (if greatorder (merror (intl:gettext "ordergreat: reordering is not allowed.")))
 876   (makorder (setq greatorder (reverse (cdr l))) '_))
 877
 878 (defmspec $orderless (l)
 879   (if lessorder (merror (intl:gettext "orderless: reordering is not allowed.")))
 880   (makorder (setq lessorder (cdr l)) '|#|))
 881
 882 (defun makorder (l char)
 883   (do ((l l (cdr l))
 884        (n 101 (1+ n)))
 885       ((null l) '$done)
 886     (alias (car l)
 887            (implode (nconc (ncons char) (mexploden n)
 888                            (exploden (stripdollar (car l))))))))
 889
 890 (defmfun $unorder ()
 891   (let ((l (delete nil
 892                  (cons '(mlist simp)
 893                        (nconc (mapcar #'(lambda (x) (remalias (getalias x))) lessorder)
 894                               (mapcar #'(lambda (x) (remalias (getalias x))) greatorder)))
 895                  :test #'eq)))
 896     (setq lessorder nil greatorder nil)
 897     l))
 898
 899 ;;;; CONCAT
 900
 901 (defmfun $concat (&rest l)
 902   "Concatenates its arguments.
 903 The arguments must evaluate to atoms. The return value is a symbol if
 904 the first argument is a symbol and a string otherwise."
 905   (when (null l)
 906     (merror (intl:gettext "concat: there must be at least one argument.")))
 907   (let ((result-is-a-string (or (numberp (car l)) (stringp (car l)))))
 908     (setq l (mapcan #'(lambda (x) (unless (atom x) (merror (intl:gettext "concat: argument must be an atom; found ~M") x)) (string* x)) l))
 909     (if result-is-a-string
 910       (coerce l 'string)
 911       (getalias (implode (cons '#\$ l))))))