src/comm2.lisp

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