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