1;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;; 2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 3;;; The data in this file contains enhancments. ;;;;; 4;;; ;;;;; 5;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;; 6;;; All rights reserved ;;;;; 7;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 8 9(in-package :maxima) 10 11;; ** (c) Copyright 1982 Massachusetts Institute of Technology ** 12 13(macsyma-module comm) 14 15(declare-top (special $exptsubst $linechar $nolabels $inflag $piece $dispflag 16 $gradefs $props $dependencies derivflag derivlist 17 $linenum $partswitch *linelabel* nn* dn* 18 $powerdisp atvars $errexp $derivsubst $dotdistrib 19 $opsubst $subnumsimp $transrun in-p substp $sqrtdispflag 20 $pfeformat dummy-variable-operators)) 21 22(defvar *islinp* nil) ; When T, sdiff is called from the function islinear. 23(defvar *atp* nil) ; When T, prevents substitution from applying to vars 24 ; bound by %sum, %product, %integrate, %limit 25 26;; op and opr properties 27 28(defvar *opr-table* (make-hash-table :test #'equal)) 29 30(defun getopr0 (x) 31 (or 32 (and (symbolp x) (get x 'opr)) 33 (and (stringp x) (gethash x *opr-table*)))) 34 35(defun getopr (x) 36 (or (getopr0 x) x)) 37 38(defun putopr (x y) 39 (or 40 (and (symbolp x) (setf (get x 'opr) y)) 41 (and (stringp x) (setf (gethash x *opr-table*) y)))) 42 43(defun remopr (x) 44 (or 45 (and (symbolp x) (remprop x 'opr)) 46 (and (stringp x) (remhash x *opr-table*)))) 47 48;; Store built-in operators, which get additional properties. 49;; These operators aren't killed by the function kill-operator. 50(defvar *mopl* nil) 51 52(mapc #'(lambda (x) (putprop (car x) (cadr x) 'op) (putopr (cadr x) (car x)) (push (cadr x) *mopl*)) 53 '((mplus "+") (mminus "-") (mtimes "*") (mexpt "**") (mexpt "^") 54 (mnctimes ".") (rat "/") (mquotient "/") (mncexpt "^^") 55 (mequal "=") (mgreaterp ">") (mlessp "<") (mleqp "<=") (mgeqp ">=") 56 (mnotequal "#") (mand "and") (mor "or") (mnot "not") (msetq ":") 57 (mdefine ":=") (mdefmacro "::=") (mquote "'") (mlist "[") 58 (mset "::") (mfactorial "!") (marrow "-->") (mprogn "(") 59 (mcond "if") (mdo "do") (mdoin "do_in"))) 60 61(mapc #'(lambda (x) (putprop (car x) (cadr x) 'op)) 62 '((mqapply $subvar) (bigfloat $bfloat))) 63 64(defmvar $exptsubst nil) 65(defmvar $partswitch nil) 66(defmvar $inflag nil) 67(defmvar $derivsubst nil) 68(defmvar $opsubst t) 69(defvar $gradefs '((mlist simp))) 70(defvar $dependencies '((mlist simp))) 71(defvar atvars '($@1 $@2 $@3 $@4)) 72(defvar in-p nil) 73(defvar substp nil) 74 75(defmvar $vect_cross nil 76 "If TRUE allows DIFF(X~Y,T) to work where ~ is defined in 77 SHARE;VECT where VECT_CROSS is set to TRUE.") 78 79(defmfun $substitute (old new &optional (expr nil three-arg?)) 80 (cond (three-arg? (maxima-substitute old new expr)) 81 (t 82 (let ((l old) (z new)) 83 (cond ((and ($listp l) ($listp (cadr l)) (null (cddr l))) 84 ($substitute (cadr l) z)) 85 ((notloreq l) (improper-arg-err l '$substitute)) 86 ((eq (caar l) 'mequal) (maxima-substitute (caddr l) (cadr l) z)) 87 (t (do ((l (cdr l) (cdr l))) 88 ((null l) z) 89 (setq z ($substitute (car l) z))))))))) 90 91;; Define an alias $psubst and a reversealias for $psubstitute 92(defprop $psubst $psubstitute alias) 93(defprop $psubstitute $psubst reversealias) 94 95;; $psubstitute is similar to $substitute. In distinction from $substitute 96;; the function $psubstitute does parallel substitution, if the first argument 97;; is a list of equations. 98(defmfun $psubstitute (old new &optional (expr nil three-arg?)) 99 (cond (three-arg? (maxima-substitute old new expr)) 100 (t 101 (let ((l old) (z new)) 102 (cond ((and ($listp l) 103 ($listp (cadr l)) 104 (null (cddr l))) 105 ;; A nested list. 106 ($psubstitute (cadr l) z)) 107 ((and ($listp l) 108 (eq (caar (cadr l)) 'mequal) 109 (null (cddr l))) 110 ;; A list with one equation. 111 ($psubstitute (cadr l) z)) 112 ((notloreq l) (improper-arg-err l '$psubstitute)) 113 ((eq (caar l) 'mequal) 114 ;; Do a substitution for one equation. 115 (maxima-substitute (caddr l) (cadr l) z)) 116 (t 117 ;; We have a list of equations. We do parallel substitution. 118 (let (gensymbol genlist eqn ($simp nil)) 119 ;; At first substitute a gensym for the expressions of 120 ;; the left hand side of the equations. 121 (do ((l (cdr l) (cdr l))) 122 ((null l) z) 123 (setq eqn (car l)) 124 (when (not (eq 'mequal (caar eqn))) 125 (improper-arg-err old '$substitute)) 126 (setq gensymbol (gensym)) 127 ;; Store the gensym and the new expression into a list. 128 (push (cons gensymbol (caddr eqn)) genlist) 129 ;; Substitute a gensym for the old expression. 130 (setq z (maxima-substitute gensymbol (cadr eqn) z))) 131 ;; Substitute the new expressions for the gensyms. 132 (do ((l genlist (cdr l))) 133 ((null l) 134 ;; Resimplify the result. 135 (let (($simp t)) (resimplify z))) 136 (setq z (maxima-substitute (cdar l) (caar l) z)))))))))) 137 138(defun maxima-substitute (x y z) ; The args to SUBSTITUTE are assumed to be simplified. 139 (let ((in-p t) (substp t)) 140 (if (and (mnump y) (= (signum1 y) 1)) 141 (let ($sqrtdispflag ($pfeformat t)) (setq z (nformat-all z)))) 142 (simplifya 143 (if (atom y) 144 (cond ((equal y -1) 145 (setq y '((mminus) 1)) (subst2 x y (nformat-all z) nil nil)) ;; negxpty and timesp don't matter in this call since (caar y) != 'mexpt 146 (t 147 (cond ((and (not (symbolp x)) 148 (functionp x)) 149 (let ((tem (gensym))) 150 (setf (get tem 'operators) 'application-operator) 151 (setf (symbol-function tem) x) 152 (setq x tem)))) 153 (subst1 x y z))) 154 (let ((negxpty (if (and (eq (caar y) 'mexpt) 155 (= (signum1 (caddr y)) 1)) 156 (mul2 -1 (caddr y)))) 157 (timesp (if (eq (caar y) 'mtimes) (setq y (nformat y))))) 158 (subst2 x y z negxpty timesp))) 159 nil))) 160 161;;Remainder of page is update from F302 --gsb 162 163;;Used in COMM2 (AT), limit, and below. 164(defvar dummy-variable-operators '(%product %sum %laplace %integrate %limit %at)) 165 166(defun subst1 (x y z) ; Y is an atom 167 (cond ((atom z) (if (equal y z) x z)) 168 ((specrepp z) (subst1 x y (specdisrep z))) 169 ((eq (caar z) 'bigfloat) z) 170 ((and (eq (caar z) 'rat) (or (equal y (cadr z)) (equal y (caddr z)))) 171 (div (subst1 x y (cadr z)) (subst1 x y (caddr z)))) 172 ((at-substp z) (subst-except-second-arg x y z)) 173 ((and (eq y t) (eq (caar z) 'mcond)) 174 (list (cons (caar z) nil) (subst1 x y (cadr z)) (subst1 x y (caddr z)) 175 (cadddr z) (subst1 x y (car (cddddr z))))) 176 (t (let ((margs (mapcar #'(lambda (z1) (subst1 x y z1)) (cdr z))) 177 (oprx (getopr x)) (opry (getopr y))) 178 (if (and $opsubst 179 (or (eq opry (caar z)) 180 (and (eq (caar z) 'rat) (eq opry 'mquotient)))) 181 (if (or (numberp x) 182 (member x '(t nil $%e $%pi $%i) :test #'eq) 183 (and (not (atom x)) 184 (not (or (eq (car x) 'lambda) 185 (eq (caar x) 'lambda))))) 186 (if (or (and (member 'array (cdar z) :test #'eq) 187 (or (and (mnump x) $subnumsimp) 188 (and (not (mnump x)) (not (atom x))))) 189 ($subvarp x)) 190 (let ((substp 'mqapply)) 191 (subst0 (list* '(mqapply) x margs) z)) 192 (merror (intl:gettext "subst: cannot substitute ~M for operator ~M in expression ~M") x y z)) 193 (subst0 (cons (cons oprx nil) margs) z)) 194 (subst0 (cons (cons (caar z) nil) margs) z)))))) 195 196(defun subst2 (x y z negxpty timesp) 197 (let (newexpt) 198 (cond ((atom z) z) 199 ((specrepp z) (subst2 x y (specdisrep z) negxpty timesp)) 200 ((at-substp z) z) ;; IS SUBST-EXCEPT-SECOND-ARG APPROPRIATE HERE ?? !! 201 ((alike1 y z) x) 202 ((and timesp (eq (caar z) 'mtimes) (alike1 y (setq z (nformat z)))) x) 203 ((and (eq (caar y) 'mexpt) (eq (caar z) 'mexpt) (alike1 (cadr y) (cadr z)) 204 (setq newexpt (cond ((alike1 negxpty (caddr z)) -1) 205 ($exptsubst (expthack (caddr y) (caddr z)))))) 206 (list '(mexpt) x newexpt)) 207 ((and $derivsubst (eq (caar y) '%derivative) (eq (caar z) '%derivative) 208 (alike1 (cadr y) (cadr z))) 209 (let ((tail (subst-diff-match (cddr y) (cdr z)))) 210 (cond ((null tail) z) 211 (t (cons (cons (caar z) nil) (cons x (cdr tail))))))) 212 (t (recur-apply #'(lambda (z1) (subst2 x y z1 negxpty timesp)) z))))) 213 214;; replace y with x in z, but leave z's second arg unchanged. 215;; This is for cases like at(integrate(x, x, a, b), [x=3]) 216;; where second arg of integrate binds a new variable x, 217;; and we do not wish to subst 3 for x inside integrand. 218(defun subst-except-second-arg (x y z) 219 (cond 220 ((member (caar z) '(%integrate %sum %product %limit %laplace)) 221 (append 222 (list (remove 'simp (car z)) ; ensure resimplification after substitution 223 (if (eq y (third z)) ; if (third z) is new var that shadows y 224 (second z) ; leave (second z) unchanged 225 (subst1 x y (second z))) ; otherwise replace y with x in (second z) 226 (third z)) ; never change integration var 227 (mapcar (lambda (z) (subst1 x y z)); do subst in limits of integral 228 (cdddr z)))) 229 ((eq (caar z) '%at) 230 ;; similar considerations here, but different structure of expression. 231 (let* 232 ((at-eqn-or-eqns (third z)) 233 (at-eqns-list (if (eq (caar at-eqn-or-eqns) 'mlist) (rest at-eqn-or-eqns) (list at-eqn-or-eqns)))) 234 (list 235 (remove 'simp (car z)) ;; ensure resimplification after substitution 236 (if (member y (mapcar #'(lambda (e) (second e)) at-eqns-list)) 237 (second z) 238 (subst1 x y (second z))) 239 `((mlist) ,@(mapcar #'(lambda (e) (list (first e) (second e) (subst1 x y (third e)))) at-eqns-list))))) 240 ((eq (caar z) '%derivative) 241 ;; and again, with yet a different structure. 242 (let* 243 ((vars-and-degrees (rest (rest z))) 244 (diff-vars (loop for v in vars-and-degrees by #'cddr collect v)) 245 (diff-degrees (loop for n in (rest vars-and-degrees) by #'cddr collect n))) 246 (append 247 (list 248 (remove 'simp (car z)) ;; ensure resimplification after substitution 249 (if (member y diff-vars) 250 (second z) 251 (subst1 x y (second z)))) 252 (apply #'append (loop for v in diff-vars for n in diff-degrees collect (list v (subst1 x y n))))))) 253 (t z))) 254 255(defun subst0 (new old) 256 (cond ((atom new) new) 257 ((alike (cdr new) (cdr old)) 258 (cond ((eq (caar new) (caar old)) old) 259 (t (simplifya (cons (cons (caar new) (member 'array (cdar old) :test #'eq)) (cdr old)) 260 nil)))) 261 ((member 'array (cdar old) :test #'eq) 262 (simplifya (cons (cons (caar new) '(array)) (cdr new)) nil)) 263 (t (simplifya new nil)))) 264 265(defun expthack (y z) 266 (prog (nn* dn* yn yd zn zd qd) 267 (cond ((and (mnump y) (mnump z)) 268 (return (if (numberp (setq y (div* z y))) y))) 269 ((atom z) (if (not (mnump y)) (return nil))) 270 ((or (ratnump z) (eq (caar z) 'mplus)) (return nil))) 271 (numden y) ; (CSIMP) sets NN* and DN* 272 (setq yn nn* yd dn*) 273 (numden z) 274 (setq zn nn* zd dn*) 275 (setq qd (cond ((and (equal zd 1) (equal yd 1)) 1) 276 ((prog2 (numden (div* zd yd)) 277 (and (equal dn* 1) (equal nn* 1))) 278 1) 279 ((equal nn* 1) (div* 1 dn*)) 280 ((equal dn* 1) nn*) 281 (t (return nil)))) 282 (numden (div* zn yn)) 283 (if (equal dn* 1) (return (div* nn* qd))))) 284 285(defun subst-diff-match (l1 l2) 286 (do ((l l1 (cddr l)) (l2 (copy-list l2)) (failed nil nil)) 287 ((null l) l2) 288 (do ((l2 l2 (cddr l2))) 289 ((null (cdr l2)) (setq failed t)) 290 (if (alike1 (car l) (cadr l2)) 291 (if (and (fixnump (cadr l)) (fixnump (caddr l2))) 292 (cond ((< (cadr l) (caddr l2)) 293 (return (rplacd (cdr l2) 294 (cons (- (caddr l2) (cadr l)) 295 (cdddr l2))))) 296 ((= (cadr l) (caddr l2)) 297 (return (rplacd l2 (cdddr l2)))) 298 (t (return (setq failed t)))) 299 (return (setq failed t))))) 300 (if failed (return nil)))) 301 302;;This probably should be a subst or macro. 303(defun at-substp (z) 304 (and *atp* 305 (or (member (caar z) '(%derivative %del) :test #'eq) 306 (member (caar z) dummy-variable-operators :test #'eq)))) 307 308(defun recur-apply (fun e) 309 (cond ((eq (caar e) 'bigfloat) e) 310 ((specrepp e) (funcall fun (specdisrep e))) 311 (t (let ((newargs (mapcar fun (cdr e)))) 312 (if (alike newargs (cdr e)) 313 e 314 (simplifya (cons (cons (caar e) (member 'array (cdar e) :test #'eq)) newargs) 315 nil)))))) 316 317(defmfun $depends (&rest args) 318 (when (oddp (length args)) 319 (merror (intl:gettext "depends: number of arguments must be even."))) 320 (do ((args args (cddr args)) 321 (l)) 322 ((null args) (i-$dependencies (nreverse l))) 323 (if ($listp (first args)) 324 (mapc #'(lambda (e) (push (depends1 e (second args)) l)) 325 (cdr (first args))) 326 (push (depends1 (first args) (second args)) l)))) 327 328(defun depends1 (x y) 329 (nonsymchk x '$depends) 330 (cons (cons x nil) (if ($listp y) (cdr y) (cons y nil)))) 331 332(defmspec $dependencies (form) 333 (i-$dependencies (cdr form))) 334 335(defun i-$dependencies (l &aux res) 336 (dolist (z l) 337 (cond 338 ((atom z) 339 (merror 340 (intl:gettext 341 "depends: argument must be a non-atomic expression; found ~M") z)) 342 (t 343 (do ((zz z (cdr zz)) 344 (y nil)) 345 ((null zz) 346 (mputprop (caar z) (setq y (reverse y)) 'depends) 347 (setq res (push (cons (ncons (caar z)) y) res)) 348 (unless (cdr $dependencies) 349 (setq $dependencies (copy-list '((mlist simp))))) 350 (add2lnc (cons (cons (caar z) nil) y) $dependencies)) 351 (cond ((and ($subvarp (cadr zz)) 352 (not (member (caar (cadr zz)) y))) 353 (setq y (push (cadr zz) y))) 354 ((not (symbolp (cadr zz))) 355 (merror 356 (intl:gettext "depends: argument must be a symbol; found ~M") 357 (cadr zz))) 358 ((and (cadr zz) 359 (not (member (cadr zz) y))) 360 (setq y (push (cadr zz) y)))))))) 361 (cons '(mlist simp) (reverse res))) 362 363(defmspec $gradef (l) 364 (setq l (cdr l)) 365 (let ((z (car l)) (n 0)) 366 (cond ((atom z) 367 (if (not (= (length l) 3)) (merror (intl:gettext "gradef: expected exactly three arguments."))) 368 (mputprop z 369 (cons (cons (cadr l) (meval (caddr l))) 370 (mget z '$atomgrad)) 371 '$atomgrad) 372 (i-$dependencies (cons (cons (ncons z) 373 ;; Append existing dependencies 374 (cons (cadr l) (mget z 'depends))) 375 nil)) 376 (add2lnc z $props) 377 z) 378 ((or (mopp1 (caar z)) (member 'array (cdar z) :test #'eq)) 379 (merror (intl:gettext "gradef: argument cannot be a built-in operator or subscripted expression; found ~M") z)) 380 ((prog2 (setq n (- (length z) (length l))) (minusp n)) 381 (wna-err '$gradef)) 382 (t (do ((zl (cdr z) (cdr zl))) ((null zl)) 383 (if (not (symbolp (car zl))) 384 (merror (intl:gettext "gradef: argument must be a symbol; found ~M") (car zl)))) 385 (setq l (nconc (mapcar #'(lambda (x) (remsimp (meval x))) 386 (cdr l)) 387 (mapcar #'(lambda (x) (list '(%derivative) z x 1)) 388 (nthcdr (- (length z) n) z)))) 389 (putprop (caar z) 390 (sublis (mapcar #'cons (cdr z) (mapcar #'stripdollar (cdr z))) 391 (cons (cdr z) l)) 392 'grad) 393 (or (cdr $gradefs) (setq $gradefs (copy-list '((mlist simp))))) 394 (add2lnc (cons (cons (caar z) nil) (cdr z)) $gradefs) z)))) 395 396(defmfun $diff (&rest args) 397 #-gcl 398 (declare (dynamic-extent args)) 399 (let (derivlist) 400 (deriv args))) 401 402(defmfun $del (e) 403 (stotaldiff e)) 404 405(defun deriv (e) 406 (prog (exp z count) 407 (cond ((null e) (wna-err '$diff)) 408 ((null (cdr e)) (return (stotaldiff (car e)))) 409 ((null (cddr e)) (nconc e '(1)))) 410 (setq exp (car e) z (setq e (copy-list e))) 411 loop (if (or (null derivlist) (member (cadr z) derivlist :test #'equal)) (go doit)) 412 ; DERIVLIST is set by $EV 413 (setq z (cdr z)) 414 loop2(cond ((cdr z) (go loop)) 415 ((null (cdr e)) (return exp)) 416 (t (go noun))) 417 doit (cond ((nonvarcheck (cadr z) '$diff)) 418 ((null (cddr z)) (wna-err '$diff)) 419 ((not (fixnump (caddr z))) (go noun)) 420 ((minusp (setq count (caddr z))) 421 (merror (intl:gettext "diff: order of derivative must be a nonnegative integer; found ~M") count))) 422 loop1(cond ((zerop count) (rplacd z (cdddr z)) (go loop2)) 423 ((equal (setq exp (sdiff exp (cadr z))) 0) (return 0))) 424 (setq count (1- count)) 425 (go loop1) 426 noun (return (diff%deriv (cons exp (cdr e)))))) 427 428(defun chainrule (e x) 429 (let (w) 430 (cond (*islinp* 431 ;; sdiff is called from the function islinear. 432 (if (and (not (atom e)) 433 (eq (caar e) '%derivative) 434 (not (freel (cdr e) x))) 435 (diff%deriv (list e x 1)) 436 0)) 437 ((atomgrad e x)) 438 ((not (setq w (mget (cond ((atom e) e) 439 ((member 'array (cdar e) :test #'eq) (caar e)) 440 ((atom (cadr e)) (cadr e)) 441 (t (caaadr e))) 442 'depends))) 443 0) 444 (t (let (derivflag) 445 (addn (mapcar 446 #'(lambda (u) 447 (let ((y (sdiff u x))) 448 (if (equal y 0) 449 0 450 (list '(mtimes) 451 (or (atomgrad e u) 452 (list '(%derivative) e u 1)) 453 y)))) 454 w) 455 nil)))))) 456 457(defun atomgrad (e x) 458 (let (y) 459 (and (atom e) (setq y (mget e '$atomgrad)) (assolike x y)))) 460 461(defun depends (e x &aux l) 462 (setq e (specrepcheck e)) 463 (cond ((alike1 e x) t) 464 ((mnump e) nil) 465 ((and (symbolp e) (setq l (mget e 'depends))) 466 ;; Go recursively through the list of dependencies. 467 ;; This code detects indirect dependencies like a(x) and x(t). 468 (dependsl l x)) 469 ((atom e) nil) 470 (t (or (depends (caar e) x) 471 (dependsl (cdr e) x))))) 472 473(defun dependsl (l x) 474 (dolist (u l) 475 (if (depends u x) (return t)))) 476 477(defun sdiff (e x) ; The args to SDIFF are assumed to be simplified. 478 ;; Remove a special representation from the variable of differentiation 479 (setq x (specrepcheck x)) 480 (cond ((alike1 e x) 1) 481 ((mnump e) 0) 482 ((or (atom e) (member 'array (cdar e) :test #'eq)) (chainrule e x)) 483 ((eq (caar e) 'mrat) (ratdx e x)) 484 ((eq (caar e) 'mpois) ($poisdiff e x)) ; Poisson series 485 ((eq (caar e) 'mplus) (addn (sdiffmap (cdr e) x) t)) 486 ((mbagp e) (cons (car e) (sdiffmap (cdr e) x))) 487 ((member (caar e) '(%sum %product) :test #'eq) (diffsumprod e x)) 488 ((eq (caar e) '%at) (diff-%at e x)) 489 ((not (depends e x)) 0) 490 ((eq (caar e) 'mtimes) (addn (sdifftimes (cdr e) x) t)) 491 ((eq (caar e) 'mexpt) (diffexpt e x)) 492 ((eq (caar e) 'mnctimes) 493 (let (($dotdistrib t)) 494 (add2 (ncmuln (cons (sdiff (cadr e) x) (cddr e)) t) 495 (ncmul2 (cadr e) (sdiff (cons '(mnctimes) (cddr e)) x))))) 496 ((and $vect_cross (eq (caar e) '|$~|)) 497 (add2* `((|$~|) ,(cadr e) ,(sdiff (caddr e) x)) 498 `((|$~|) ,(sdiff (cadr e) x) ,(caddr e)))) 499 ((eq (caar e) 'mncexpt) (diffncexpt e x)) 500 ((member (caar e) '(%log %plog) :test #'eq) 501 (sdiffgrad (cond ((and (not (atom (cadr e))) (eq (caaadr e) 'mabs)) 502 (cons (car e) (cdadr e))) 503 (t e)) 504 x)) 505 ((eq (caar e) '%derivative) 506 (cond ((or (atom (cadr e)) (member 'array (cdaadr e) :test #'eq)) (chainrule e x)) 507 ((freel (cddr e) x) (diff%deriv (cons (sdiff (cadr e) x) (cddr e)))) 508 (t (diff%deriv (list e x 1))))) 509 ((member (caar e) '(%binomial $beta) :test #'eq) 510 (let ((efact ($makefact e))) 511 (mul2 (factor (sdiff efact x)) (div e efact)))) 512 ((eq (caar e) '%integrate) (diffint e x)) 513 ((eq (caar e) '%laplace) (difflaplace e x)) 514 ((eq (caar e) '%at) (diff-%at e x)) 515; This rule is not correct. We cut it out. 516; ((member (caar e) '(%realpart %imagpart) :test #'eq) 517; (list (cons (caar e) nil) (sdiff (cadr e) x))) 518 ((and (eq (caar e) 'mqapply) 519 (eq (caaadr e) '$%f)) 520 ;; Handle %f, hypergeometric function 521 ;; 522 ;; The derivative of %f[p,q]([a1,...,ap],[b1,...,bq],z) is 523 ;; 524 ;; a1*a2*...*ap/(b1*b2*...*bq) 525 ;; *%f[p,q]([a1+1,a2+1,...,ap+1],[b1+1,b2+1,...,bq+1],z) 526 (let* ((arg1 (cdr (third e))) 527 (arg2 (cdr (fourth e))) 528 (v (fifth e))) 529 (mul (sdiff v x) 530 (div (mull arg1) (mull arg2)) 531 `((mqapply) (($%f array) ,(length arg1) ,(length arg2)) 532 ((mlist) ,@(incr1 arg1)) 533 ((mlist) ,@(incr1 arg2)) 534 ,v)))) 535 (t (sdiffgrad e x)))) 536 537(defun sdiffgrad (e x) 538 (let ((fun (caar e)) grad args result) 539 (cond ((and (eq fun 'mqapply) (zl-get (caaadr e) 'grad)) 540 ;; Change the array function f[n](x) to f(n,x), call sdiffgrad again. 541 (setq result 542 (sdiffgrad (cons (cons (caaadr e) nil) 543 (append (cdadr e) (cddr e))) 544 x)) 545 ;; If noun form for f(n,x), adjust the noun form for f[n](x) 546 (if (isinop result '%derivative) 547 (if (not (depends e x)) 548 0 549 (diff%deriv (list e x 1))) 550 result)) 551 552 ;; extension for pdiff. 553 ((and (get '$pderivop 'operators) (funcall 'sdiffgrad-pdiff e x))) 554 555 ;; two line extension for hypergeometric. 556 ((and (equal fun '$hypergeometric) (get '$hypergeometric 'operators)) 557 (funcall 'diff-hypergeometric (second e) (third e) (fourth e) x)) 558 559 ((or (eq fun 'mqapply) (null (setq grad (zl-get fun 'grad)))) 560 (if (not (depends e x)) 0 (diff%deriv (list e x 1)))) 561 ((not (= (length (cdr e)) (length (car grad)))) 562 (merror (intl:gettext "~:M: expected exactly ~M arguments.") 563 fun 564 (length (car grad)))) 565 (t 566 (setq args (sdiffmap (cdr e) x)) 567 (setq result 568 (addn 569 (mapcar 570 #'mul2 571 (cdr 572 (substitutel 573 (cdr e) 574 (car grad) 575 (do ((l1 (cdr grad) (cdr l1)) 576 (args args (cdr args)) 577 (l2)) 578 ((null l1) (cons '(mlist) (nreverse l2))) 579 (setq l2 580 (cons (cond ((equal (car args) 0) 0) 581 ((functionp (car l1)) 582 ;; Evaluate a lambda expression 583 ;; given as a derivative. 584 (apply (car l1) (cdr e))) 585 (t (car l1))) 586 l2))))) 587 args) 588 t)) 589 (if (or (null result) (not (freeof nil result))) 590 ;; A derivative has returned NIL. Return a noun form. 591 (if (not (depends e x)) 592 0 593 (diff%deriv (list e x 1))) 594 result))))) 595 596(defun sdiffmap (e x) 597 (mapcar #'(lambda (term) (sdiff term x)) e)) 598 599(defun sdifftimes (l x) 600 (prog (term left out) 601 loop (setq term (car l) l (cdr l)) 602 (setq out (cons (muln (cons (sdiff term x) (append left l)) t) out)) 603 (if (null l) (return out)) 604 (setq left (cons term left)) 605 (go loop))) 606 607(defun diffexpt (e x) 608 (if (mnump (caddr e)) 609 (mul3 (caddr e) (power (cadr e) (addk (caddr e) -1)) (sdiff (cadr e) x)) 610 (mul2 e (add2 (mul3 (power (cadr e) -1) (caddr e) (sdiff (cadr e) x)) 611 (mul2 (simplifya (list '(%log) (cadr e)) t) 612 (sdiff (caddr e) x)))))) 613 614(defun diff%deriv (e) 615 (let (derivflag) 616 (simplifya (cons '(%derivative) e) t))) 617 618 619;; grad properties 620 621(let ((header '(x))) 622 (mapc #'(lambda (z) (putprop (car z) (cons header (cdr z)) 'grad)) 623 ;; All these GRAD templates have been simplified and then the SIMP flags 624 ;; (which are unnecessary) have been removed to save core space. 625 '((%log ((mexpt) x -1)) (%plog ((mexpt) x -1)) 626 (%gamma ((mtimes) ((mqapply) (($psi array) 0) x) ((%gamma) x))) 627 (mfactorial ((mtimes) ((mqapply) (($psi array) 0) ((mplus) 1 x)) ((mfactorial) x))) 628 (%sin ((%cos) x)) 629 (%cos ((mtimes) -1 ((%sin) x))) 630 (%tan ((mexpt) ((%sec) x) 2)) 631 (%cot ((mtimes) -1 ((mexpt) ((%csc) x) 2))) 632 (%sec ((mtimes) ((%sec) x) ((%tan) x))) 633 (%csc ((mtimes) -1 ((%cot) x) ((%csc) x))) 634 (%asin ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) ((rat) -1 2))) 635 (%acos ((mtimes) -1 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) 636 ((rat) -1 2)))) 637 (%atan ((mexpt) ((mplus) 1 ((mexpt) x 2)) -1)) 638 (%acot ((mtimes) -1 ((mexpt) ((mplus) 1 ((mexpt) x 2)) -1))) 639 (%acsc ((mtimes) -1 640 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x -2))) 641 ((rat) -1 2)) 642 ((mexpt) x -2))) 643 (%asec ((mtimes) 644 ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x -2))) 645 ((rat) -1 2)) 646 ((mexpt) x -2))) 647 (%sinh ((%cosh) x)) 648 (%cosh ((%sinh) x)) 649 (%tanh ((mexpt) ((%sech) x) 2)) 650 (%coth ((mtimes) -1 ((mexpt) ((%csch) x) 2))) 651 (%sech ((mtimes) -1 ((%sech) x) ((%tanh) x))) 652 (%csch ((mtimes) -1 ((%coth) x) ((%csch) x))) 653 (%asinh ((mexpt) ((mplus) 1 ((mexpt) x 2)) ((rat) -1 2))) 654 (%acosh ((mexpt) ((mplus) -1 ((mexpt) x 2)) ((rat) -1 2))) 655 (%atanh ((mexpt) ((mplus) 1 ((mtimes) -1 ((mexpt) x 2))) -1)) 656 (%acoth ((mtimes) -1 ((mexpt) ((mplus) -1 ((mexpt) x 2)) -1))) 657 (%asech ((mtimes) -1 658 ((mexpt) ((mplus) -1 ((mexpt) x -2)) ((rat) -1 2)) 659 ((mexpt) x -2))) 660 (%acsch ((mtimes) -1 661 ((mexpt) ((mplus) 1 ((mexpt) x -2)) ((rat) -1 2)) 662 ((mexpt) x -2))) 663 (mabs ((mtimes) x ((mexpt) ((mabs) x) -1))) 664 (%erf ((mtimes) 2 ((mexpt) $%pi ((rat) -1 2)) 665 ((mexpt) $%e ((mtimes) -1 ((mexpt) x 2))))) 666 ))) 667 668(defprop $atan2 ((x y) ((mtimes) y ((mexpt) ((mplus) ((mexpt) x 2) ((mexpt) y 2)) -1)) 669 ((mtimes) -1 x ((mexpt) ((mplus) ((mexpt) x 2) ((mexpt) y 2)) -1))) 670 grad) 671 672(defprop $li 673 ((n x) 674; Do not put a noun form on the property list, but NIL. 675; SDIFFGRAD generates the noun form. 676; ((%derivative) ((mqapply) (($li array) n) x) n 1) 677 nil 678 ((mtimes) ((mqapply) (($li array) ((mplus) -1 n)) x) ((mexpt) x -1))) 679 grad) 680 681(defprop $psi 682 ((n x) 683; Do not put a noun form on the property list, but NIL. 684; SDIFFGRAD generates the noun form. 685 nil 686 ((mqapply) (($psi array) ((mplus) 1 n)) x)) 687 grad) 688 689(defun atvarschk (argl) 690 (do ((largl (length argl) (1- largl)) 691 (latvrs (length atvars)) 692 (l)) 693 ((not (< latvrs largl)) (nconc atvars l)) 694 (setq l (cons (implode (cons '$ (cons '@ (mexploden largl)))) l)))) 695 696(defun notloreq (x) 697 (or (atom x) 698 (not (member (caar x) '(mlist mequal) :test #'eq)) 699 (and (eq (caar x) 'mlist) 700 (dolist (u (cdr x)) 701 (if (not (mequalp u)) (return t)))))) 702 703(defun substitutel (l1 l2 e) 704 "l1 is a list of expressions. l2 is a list of variables. For each 705 element in list l2, substitute corresponding element of l1 into e" 706 (do ((l1 l1 (cdr l1)) 707 (l2 l2 (cdr l2))) 708 ((null l1) e) 709 (setq e (maxima-substitute (car l1) (car l2) e)))) 710 711(defun union* (a b) 712 (do ((a a (cdr a)) 713 (x b)) 714 ((null a) x) 715 (if (not (memalike (car a) b)) (setq x (cons (car a) x))))) 716 717(defun intersect* (a b) 718 (do ((a a (cdr a)) 719 (x)) 720 ((null a) x) 721 (if (memalike (car a) b) (setq x (cons (car a) x))))) 722 723(defun nthelem (n e) 724 (car (nthcdr (1- n) e))) 725 726(defun delsimp (e) 727 (delete 'simp (copy-list e) :count 1 :test #'eq)) 728 729(defun remsimp (e) 730 (if (atom e) e (cons (delsimp (car e)) (mapcar #'remsimp (cdr e))))) 731 732(defmfun $trunc (e) 733 (cond ((atom e) e) 734 ((eq (caar e) 'mplus) (cons (append (car e) '(trunc)) (cdr e))) 735 ((mbagp e) (cons (car e) (mapcar #'$trunc (cdr e)))) 736 ((specrepp e) ($trunc (specdisrep e))) 737 (t e))) 738 739(defun nonvarcheck (e fn) 740 (if (or (mnump e) 741 (maxima-integerp e) 742 ($constantp e) 743 (and (not (atom e)) (not (eq (caar e) 'mqapply)) (mopp1 (caar e)))) 744 (merror (intl:gettext "~:M: second argument must be a variable; found ~M") fn e))) 745 746(defmspec $ldisplay (form) 747 (disp1 (cdr form) t t)) 748 749(defmfun $ldisp (&rest args) 750 #-gcl 751 (declare (dynamic-extent args)) 752 (disp1 args t nil)) 753 754(defmspec $display (form) 755 (disp1 (cdr form) nil t)) 756 757(defmfun $disp (&rest args) 758 #-gcl 759 (declare (dynamic-extent args)) 760 (disp1 args nil nil)) 761 762(defun disp1 (ll lablist eqnsp) 763 (if lablist (setq lablist (cons '(mlist simp) nil))) 764 (do ((ll ll (cdr ll)) 765 (l) 766 (ans) 767 ($dispflag t) 768 (tim 0)) 769 ((null ll) (or lablist '$done)) 770 (setq l (car ll) ans (if eqnsp (meval l) l)) 771 (if (and eqnsp (not (mequalp ans))) 772 (setq ans (list '(mequal simp) (disp2 l) ans))) 773 (if lablist (nconc lablist (cons (elabel ans) nil))) 774 (setq tim (get-internal-run-time)) 775 (let ((*display-labels-p* (not (null lablist)))) 776 (declare (special *display-labels-p*)) 777 (displa (list '(mlabel) (if lablist *linelabel*) ans))) 778 (mterpri) 779 (timeorg tim))) 780 781(defun disp2 (e) 782 (cond ((atom e) e) 783 ((eq (caar e) 'mqapply) 784 (cons '(mqapply) (cons (cons (caadr e) (mapcar #'meval (cdadr e))) 785 (mapcar #'meval (cddr e))))) 786 ((eq (caar e) 'msetq) (disp2 (cadr e))) 787 ((eq (caar e) 'mset) (disp2 (meval (cadr e)))) 788 ((eq (caar e) 'mlist) (cons (car e) (mapcar #'disp2 (cdr e)))) 789 ((mspecfunp (caar e)) e) 790 (t (cons (car e) (mapcar #'meval (cdr e)))))) 791 792; Construct a new intermediate result label, 793; and bind it to the expression e. 794; The global flag $NOLABELS is ignored; the label is always bound. 795; Otherwise (if ELABEL were to observe $NOLABELS) it would be 796; impossible to programmatically refer to intermediate result expression. 797 798(defun elabel (e) 799 (if (not (checklabel $linechar)) (setq $linenum (1+ $linenum))) 800 (let (($nolabels nil)) ; <-- This is pretty ugly. MAKELABEL should take another argument. 801 (makelabel $linechar)) 802 (setf (symbol-value *linelabel*) e) 803 *linelabel*) 804 805(defmfun $dispterms (e) 806 (cond ((or (atom e) (eq (caar e) 'bigfloat)) (displa e)) 807 ((specrepp e) ($dispterms (specdisrep e))) 808 (t (let (($dispflag t)) 809 (mterpri) 810 (displa (getop (mop e))) 811 (do ((e (if (and (eq (caar e) 'mplus) (not $powerdisp)) 812 (reverse (cdr e)) 813 (margs e)) 814 (cdr e))) ((null e)) (mterpri) (displa (car e)) (mterpri))) 815 (mterpri))) 816 '$done) 817 818(defmfun $dispform (e &optional (flag nil flag?)) 819 (when (and flag? (not (eq flag '$all))) 820 (merror (intl:gettext "dispform: second argument, if present, must be 'all'; found ~M") flag)) 821 (if (or (atom e) 822 (atom (setq e (if flag? (nformat-all e) (nformat e)))) 823 (member 'simp (cdar e) :test #'eq)) 824 e 825 (cons (cons (caar e) (cons 'simp (cdar e))) 826 (if (and (eq (caar e) 'mplus) (not $powerdisp)) 827 (reverse (cdr e)) 828 (cdr e))))) 829 830;;; These functions implement the Macsyma functions $op and $operatorp. 831;;; Dan Stanger 832(defmfun $op (expr) 833 ($part expr 0)) 834 835(defmfun $operatorp (expr oplist) 836 (if ($listp oplist) 837 ($member ($op expr) oplist) 838 (equal ($op expr) oplist))) 839 840(defmfun $part (&rest args) 841 #-gcl 842 (declare (dynamic-extent args)) 843 (mpart args nil nil $inflag '$part)) 844 845(defmfun $inpart (&rest args) 846 #-gcl 847 (declare (dynamic-extent args)) 848 (mpart args nil nil t '$inpart)) 849 850(defmspec $substpart (l) 851 (let ((substp t)) 852 (mpart (cdr l) t nil $inflag '$substpart))) 853 854(defmspec $substinpart (l) 855 (let ((substp t)) 856 (mpart (cdr l) t nil t '$substinpart))) 857 858(defun part1 (arglist substflag dispflag inflag) ; called only by TRANSLATE 859 (let ((substp t)) 860 (mpart arglist substflag dispflag inflag '$substpart))) 861 862(defun mpart (arglist substflag dispflag inflag fn) 863 (prog (substitem arg arg1 exp exp1 exp* sevlist count prevcount n specp 864 lastelem lastcount) 865 (setq specp (or substflag dispflag)) 866 (if substflag (setq substitem (car arglist) arglist (cdr arglist))) 867 (if (null arglist) (wna-err '$part)) 868 (setq exp (if substflag (meval (car arglist)) (car arglist))) 869 (when (null (setq arglist (cdr arglist))) 870 (setq $piece exp) 871 (return (cond (substflag (meval substitem)) 872 (dispflag (box exp dispflag)) 873 (t exp)))) 874 (cond ((not inflag) 875 (cond ((or (and ($listp exp) (null (cdr arglist))) 876 (and ($matrixp exp) 877 (or (null (cdr arglist)) (null (cddr arglist))))) 878 (setq inflag t)) 879 ((not specp) (setq exp (nformat exp))) 880 (t (setq exp (nformat-all exp))))) 881 ((specrepp exp) (setq exp (specdisrep exp)))) 882 (when (and (atom exp) (null $partswitch)) 883 (merror (intl:gettext "~:M: argument must be a non-atomic expression; found ~:M") fn exp)) 884 (when (and inflag specp) 885 (setq exp (copy-tree exp))) 886 (when substflag 887 ;; Replace all occurrences of 'rat with 'mquotient when in subst. 888 (setq exp (let (($simp nil)) (maxima-substitute 'mquotient 'rat exp)))) 889 (setq exp* exp) 890 start (cond ((or (atom exp) (eq (caar exp) 'bigfloat)) (go err)) 891 ((equal (setq arg (if substflag (meval (car arglist)) (car arglist))) 892 0) 893 (setq arglist (cdr arglist)) 894 (cond ((mnump substitem) 895 (merror (intl:gettext "~:M: argument cannot be a number; found ~M") fn substitem)) 896 ((and specp arglist) 897 (if (eq (caar exp) 'mqapply) 898 (prog2 (setq exp (cadr exp)) (go start)) 899 ;; NOT CLEAR WHAT IS INVALID HERE. OH WELL. 900 (merror (intl:gettext "~:M: invalid operator.") fn))) 901 (t (setq $piece (getop (mop exp))) 902 (return 903 (cond (substflag 904 (setq substitem (getopr (meval substitem))) 905 (cond ((mnump substitem) 906 (merror (intl:gettext "~:M: argument cannot be a number; found ~M") fn substitem)) 907 ((not (atom substitem)) 908 (if (not (eq (caar exp) 'mqapply)) 909 (rplaca (rplacd exp (cons (car exp) 910 (cdr exp))) 911 '(mqapply))) 912 (rplaca (cdr exp) substitem) 913 (return (resimplify exp*))) 914 ((eq (caar exp) 'mqapply) 915 (rplacd exp (cddr exp)))) 916 (rplaca exp (cons substitem 917 (if (and (member 'array (cdar exp) :test #'eq) 918 (not (mopp substitem))) 919 '(array)))) 920 (resimplify exp*)) 921 (dispflag 922 (rplacd exp (cdr (box (copy-tree exp) dispflag))) 923 (rplaca exp (if (eq dispflag t) 924 '(mbox) 925 '(mlabox))) 926 (resimplify exp*)) 927 (t (when arglist (setq exp $piece) (go a)) 928 $piece)))))) 929 ((not (atom arg)) (go several)) 930 ((not (fixnump arg)) 931 (merror (intl:gettext "~:M: argument must be an integer; found ~M") fn arg)) 932 ((< arg 0) (go bad))) 933 (if (eq (caar exp) 'mqapply) (setq exp (cdr exp))) 934 loop (cond ((not (zerop arg)) (setq arg (1- arg) exp (cdr exp)) 935 (if (null exp) (go err)) 936 (go loop)) 937 ((null (setq arglist (cdr arglist))) 938 (return (cond (substflag (setq $piece (resimplify (car exp))) 939 (rplaca exp (meval substitem)) 940 (resimplify exp*)) 941 (dispflag (setq $piece (resimplify (car exp))) 942 (rplaca exp (box (car exp) dispflag)) 943 (resimplify exp*)) 944 (inflag (setq $piece (car exp))) 945 (t (setq $piece (simplify (car exp)))))))) 946 (setq exp (car exp)) 947 a (cond ((and (not inflag) (not specp)) (setq exp (nformat exp))) 948 ((specrepp exp) (setq exp (specdisrep exp)))) 949 (go start) 950 err (cond ((eq $partswitch 'mapply) 951 (merror (intl:gettext "~M: invalid index ~M of list or matrix.") fn (car arglist))) 952 ($partswitch (return (setq $piece '$end))) 953 (t (merror (intl:gettext "~:M: fell off the end.") fn))) 954 bad (improper-arg-err arg fn) 955 several 956 (if (or (not (member (caar arg) '(mlist $allbut) :test #'eq)) (cdr arglist)) 957 (go bad)) 958 (setq exp1 (cons (caar exp) (if (member 'array (cdar exp) :test #'eq) '(array)))) 959 (if (eq (caar exp) 'mqapply) 960 (setq sevlist (list (cadr exp) exp1) exp (cddr exp)) 961 (setq sevlist (ncons exp1) exp (cdr exp))) 962 (setq arg1 (cdr arg) prevcount 0 exp1 exp) 963 (dolist (arg* arg1) 964 (if (not (fixnump arg*)) 965 (merror (intl:gettext "~:M: argument must be an integer; found ~M") fn arg*))) 966 (when (and specp (eq (caar arg) 'mlist)) 967 (if substflag (setq lastelem (car (last arg1)))) 968 (setq arg1 (sort (copy-list arg1) #'<))) 969 (when (eq (caar arg) '$allbut) 970 (setq n (length exp)) 971 (dolist (i arg1) 972 (if (or (< i 1) (> i n)) 973 (merror (intl:gettext "~:M: index must be in range 1 to ~M, inclusive; found ~M") fn n i))) 974 (do ((i n (1- i)) (arg2)) 975 ((= i 0) (setq arg1 arg2)) 976 (if (not (member i arg1 :test #'equal)) (setq arg2 (cons i arg2)))) 977 (if substflag (setq lastelem (car (last arg1))))) 978 (if (null arg1) (if specp (go bad) (go end))) 979 (if substflag (setq lastcount lastelem)) 980 sevloop 981 (if specp 982 (setq count (- (car arg1) prevcount) prevcount (car arg1)) 983 (setq count (car arg1))) 984 (if (< count 1) (go bad)) 985 (if (and substflag (< (car arg1) lastelem)) 986 (setq lastcount (1- lastcount))) 987 count(cond ((null exp) (go err)) 988 ((not (= count 1)) (setq count (1- count) exp (cdr exp)) (go count))) 989 (setq sevlist (cons (car exp) sevlist)) 990 (setq arg1 (cdr arg1)) 991 end (cond ((null arg1) 992 (setq sevlist (nreverse sevlist)) 993 (setq $piece (if (or inflag (not specp)) 994 (simplify sevlist) 995 (resimplify sevlist))) 996 (return (cond (substflag (rplaca (nthcdr (1- lastcount) exp1) 997 (meval substitem)) 998 (resimplify exp*)) 999 (dispflag (rplaca exp (box (car exp) dispflag)) 1000 (resimplify exp*)) 1001 (t $piece)))) 1002 (substflag (if (null (cdr exp)) (go err)) 1003 (rplaca exp (cadr exp)) (rplacd exp (cddr exp))) 1004 (dispflag (rplaca exp (box (car exp) dispflag)) 1005 (setq exp (cdr exp))) 1006 (t (setq exp exp1))) 1007 (go sevloop))) 1008 1009(defun getop (x) 1010 (or (and (symbolp x) (get x 'op)) x)) 1011 1012(defmfun $listp (x) 1013 (and (not (atom x)) 1014 (not (atom (car x))) 1015 (eq (caar x) 'mlist))) 1016 1017(defmfun $cons (x e) 1018 (atomchk (setq e (format1 e)) '$cons t) 1019 (mcons-exp-args e (cons x (margs e)))) 1020 1021(defmfun $endcons (x e) 1022 (atomchk (setq e (format1 e)) '$endcons t) 1023 (mcons-exp-args e (append (margs e) (ncons x)))) 1024 1025(defmfun $reverse (e) 1026 (atomchk (setq e (format1 e)) '$reverse nil) 1027 (mcons-exp-args e (reverse (margs e)))) 1028 1029(defmfun $append (&rest args) 1030 (if (null args) 1031 '((mlist simp)) 1032 (let ((arg1 (specrepcheck (first args))) op arrp) 1033 (atomchk arg1 '$append nil) 1034 (setq op (mop arg1) 1035 arrp (if (member 'array (cdar arg1) :test #'eq) t)) 1036 (mcons-exp-args 1037 arg1 1038 (apply #'append 1039 (mapcar #'(lambda (u) 1040 (atomchk (setq u (specrepcheck u)) '$append nil) 1041 (unless (and (alike1 op (mop u)) 1042 (eq arrp (if (member 'array (cdar u) :test #'eq) t))) 1043 (merror (intl:gettext "append: operators of arguments must all be the same."))) 1044 (margs u)) 1045 args)))))) 1046 1047(defun mcons-exp-args (e args) 1048 (if (eq (caar e) 'mqapply) 1049 (list* (delsimp (car e)) (cadr e) args) 1050 (cons (delsimp (car e)) args))) 1051 1052(defmfun $member (x e) 1053 (atomchk (setq e ($totaldisrep e)) '$member t) 1054 (if (memalike ($totaldisrep x) (margs e)) t)) 1055 1056(defun atomchk (e fun 2ndp) 1057 (if (or (atom e) (eq (caar e) 'bigfloat)) 1058 (merror (intl:gettext "~:M: ~Margument must be a non-atomic expression; found ~M") fun (if 2ndp "2nd " "") e))) 1059 1060(defun format1 (e) 1061 (cond (($listp e) e) 1062 ($inflag (specrepcheck e)) 1063 (t (nformat e)))) 1064 1065(defmfun $first (e) 1066 (atomchk (setq e (format1 e)) '$first nil) 1067 (if (null (cdr e)) (merror (intl:gettext "first: empty argument."))) 1068 (car (margs e))) 1069 1070;; This macro is used to create functions second thru tenth. 1071;; Rather than try to modify mformat for ~:R, use the quoted symbol 1072 1073(macrolet ((make-nth (si i) 1074 (let ((sim (intern (concatenate 'string "$" (symbol-name si))))) 1075 `(defun ,sim (e) 1076 (atomchk (setq e (format1 e)) ',sim nil) 1077 (if (< (length (margs e)) ,i) 1078 (merror (intl:gettext "~:M: no such element in ~M") ',sim e)) 1079 (,si (margs e)))))) 1080 1081 (make-nth second 2) 1082 (make-nth third 3) 1083 (make-nth fourth 4) 1084 (make-nth fifth 5) 1085 (make-nth sixth 6) 1086 (make-nth seventh 7) 1087 (make-nth eighth 8) 1088 (make-nth ninth 9) 1089 (make-nth tenth 10)) 1090 1091(defmfun $rest (e &optional (n 1 n?)) 1092 (prog (m fun fun1 revp) 1093 (when (and n? (equal n 0)) 1094 (return e)) 1095 (atomchk (setq m (format1 e)) '$rest nil) 1096 (cond ((and n? (not (fixnump n))) 1097 (merror (intl:gettext "rest: second argument, if present, must be an integer; found ~M") n)) 1098 ((minusp n) 1099 (setq n (- n) revp t))) 1100 (if (< (length (margs m)) n) 1101 (if $partswitch 1102 (return '$end) 1103 (merror (intl:gettext "rest: fell off the end.")))) 1104 (setq fun (car m)) 1105 (when (eq (car fun) 'mqapply) 1106 (setq fun1 (cadr m) 1107 m (cdr m))) 1108 (setq m (cdr m)) 1109 (when revp (setq m (reverse m))) 1110 (setq m (nthcdr n m)) 1111 (setq m (cons (if (eq (car fun) 'mlist) fun (delsimp fun)) 1112 (if revp (nreverse m) m))) 1113 (when (eq (car fun) 'mqapply) 1114 (return (cons (car m) (cons fun1 (cdr m))))) 1115 (return m))) 1116 1117(defmfun $last (e) 1118 (atomchk (setq e (format1 e)) '$last nil) 1119 (when (null (cdr e)) 1120 (merror (intl:gettext "last: empty argument."))) 1121 (car (last e))) 1122 1123(defmfun $firstn (e n) 1124 (atomchk (setq e (format1 e)) '$firstn nil) 1125 (if (or (not (fixnump n)) (minusp n)) 1126 (merror (intl:gettext "firstn: second argument must be a nonnegative integer; found: ~M") n)) 1127 (let ((m ($length e))) 1128 (if (> n m) 1129 ($rest e 0) 1130 ($rest e (- n m))))) 1131 1132(defmfun $lastn (e n) 1133 (atomchk (setq e (format1 e)) '$lastn nil) 1134 (if (or (not (fixnump n)) (minusp n)) 1135 (merror (intl:gettext "lastn: second argument must be a nonnegative integer; found: ~M") n)) 1136 (let ((m ($length e))) 1137 (if (> n m) 1138 ($rest e 0) 1139 ($rest e (- m n))))) 1140 1141(defmfun $args (e) 1142 (atomchk (setq e (format1 e)) '$args nil) 1143 (cons '(mlist) (margs e))) 1144 1145(defmfun $delete (x l &optional (n -1 n?)) 1146 (when (and n? (or (not (fixnump n)) (minusp n))) ; if n is set, it must be a nonneg fixnum 1147 (merror (intl:gettext "delete: third argument, if present, must be a nonnegative integer; found ~M") n)) 1148 (atomchk (setq l (specrepcheck l)) '$delete t) 1149 (setq x (specrepcheck x) 1150 l (cons (delsimp (car l)) (copy-list (cdr l)))) 1151 (do ((l1 (if (eq (caar l) 'mqapply) (cdr l) l))) 1152 ((or (null (cdr l1)) (zerop n)) l) 1153 (if (alike1 x (specrepcheck (cadr l1))) 1154 (progn 1155 (decf n) 1156 (rplacd l1 (cddr l1))) 1157 (setq l1 (cdr l1))))) 1158 1159(defmfun $length (e) 1160 (setq e (cond (($listp e) e) 1161 ((or $inflag (not ($ratp e))) (specrepcheck e)) 1162 (t ($ratdisrep e)))) 1163 (cond ((symbolp e) (merror (intl:gettext "length: argument cannot be a symbol; found ~:M") e)) 1164 ((or (numberp e) (eq (caar e) 'bigfloat)) 1165 (if (and (not $inflag) (mnegp e)) 1166 1 1167 (merror (intl:gettext "length: argument cannot be a number; found ~:M") e))) 1168 ((or $inflag (not (member (caar e) '(mtimes mexpt) :test #'eq))) (length (margs e))) 1169 ((eq (caar e) 'mexpt) 1170 (if (and (alike1 (caddr e) '((rat simp) 1 2)) $sqrtdispflag) 1 2)) 1171 (t (length (cdr (nformat e)))))) 1172 1173(defmfun $atom (x) 1174 (setq x (specrepcheck x)) (or (atom x) (eq (caar x) 'bigfloat))) 1175 1176(defmfun $symbolp (x) 1177 (setq x (specrepcheck x)) (symbolp x)) 1178 1179(defmfun $num (e) 1180 (let (x) 1181 (cond ((atom e) e) 1182 ((eq (caar e) 'mrat) ($ratnumer e)) 1183 ((eq (caar e) 'rat) (cadr e)) 1184 ((eq (caar (setq x (nformat e))) 'mquotient) (simplify (cadr x))) 1185 ((and (eq (caar x) 'mminus) (not (atom (setq x (cadr x)))) 1186 (eq (caar x) 'mquotient)) 1187 (simplify (list '(mtimes) -1 (cadr x)))) 1188 (t e)))) 1189 1190(defmfun $denom (e) 1191 (cond ((atom e) 1) 1192 ((eq (caar e) 'mrat) ($ratdenom e)) 1193 ((eq (caar e) 'rat) (caddr e)) 1194 ((or (eq (caar (setq e (nformat e))) 'mquotient) 1195 (and (eq (caar e) 'mminus) (not (atom (setq e (cadr e)))) 1196 (eq (caar e) 'mquotient))) 1197 (simplify (caddr e))) 1198 (t 1))) 1199 1200(defmfun $entier (e) (take '($floor) e)) 1201 1202(defmfun $fix (e) (take '($floor) e)) 1203 1204;; Evaluate THUNK ignoring any error that might occur. If the THUNK 1205;; returns a good number (not infinity or NaN), then return the 1206;; number. Otherwise, return NIL to indicate that we the computation 1207;; failed. This is a pretty brute-force approach. 1208(defun try-float-computation (thunk) 1209 (let ((errcatch (cons bindlist loclist)) 1210 (*mdebug* nil)) 1211 (declare (special errcatch)) 1212 (let ((result (errset (funcall thunk)))) 1213 (labels ((bad-number-p (num) 1214 (if (complexp num) 1215 (or (bad-number-p (realpart num)) 1216 (bad-number-p (imagpart num))) 1217 (or (float-inf-p num) 1218 (float-nan-p num))))) 1219 (if (and result (bad-number-p (car result))) 1220 nil 1221 (car result)))))) 1222 1223(defmfun $float (e) 1224 (cond ((numberp e) 1225 (let ((e1 (float e))) (if (float-inf-p e1) (signal 'floating-point-overflow) e1))) 1226 ((and (symbolp e) (mget e '$numer))) 1227 ((or (atom e) (member 'array (cdar e) :test #'eq)) e) 1228 ((eq (caar e) 'rat) (fpcofrat e)) 1229 ((eq (caar e) 'bigfloat) (fp2flo e)) 1230 ((member (caar e) '(mexpt mncexpt) :test #'eq) 1231 ;; avoid x^2 -> x^2.0, allow %e^%pi -> 23.14 1232 (let ((res (recur-apply #'$float e))) 1233 (if (floatp res) 1234 res 1235 (list (ncons (caar e)) ($float (cadr e)) (caddr e))))) 1236 ((and (eq (caar e) '%log) 1237 (complex-number-p (second e) '$ratnump)) 1238 ;; Basically we try to compute float(log(x)) as directly as 1239 ;; possible, expecting Lisp to return some error if it can't. 1240 ;; Then we do a more complicated approach to compute the 1241 ;; result. However, gcl and ecl don't signal errors in these 1242 ;; cases, so we always use the complicated approach for these lisps. 1243 (let ((n (second e))) 1244 (cond ((integerp n) 1245 ;; float(log(int)). First try to compute (log 1246 ;; (float n)). If that works, we're done. 1247 ;; Otherwise we need to do more. 1248 (to (or (try-float-computation #'(lambda () 1249 (log (float n)))) 1250 (let ((m (integer-length n))) 1251 ;; Write n as (n/2^m)*2^m where m is the number of 1252 ;; bits in n. Then log(n) = log(2^m) + log(n/2^m). 1253 ;; n/2^m is approximately 1, so converting that to a 1254 ;; float is no problem. log(2^m) = m * log(2). 1255 (+ (* m (log 2e0)) 1256 (log (float (/ n (ash 1 m))))))))) 1257 (($ratnump n) 1258 ;; float(log(n/m)) where n and m are integers. Try computing 1259 ;; it first. If it fails, compute as log(n) - log(m). 1260 (let ((try (try-float-computation #'(lambda() 1261 (log (fpcofrat n)))))) 1262 (if try 1263 (to try) 1264 (sub ($float `((%log) ,(second n))) 1265 ($float `((%log) ,(third n))))))) 1266 ((complex-number-p n 'integerp) 1267 ;; float(log(n+m*%i)). 1268 (let ((re ($realpart n)) 1269 (im ($imagpart n))) 1270 (to (or (try-float-computation #'(lambda () 1271 (log (complex (float re) 1272 (float im))))) 1273 (let* ((size (max (integer-length re) 1274 (integer-length im))) 1275 (scale (ash 1 size))) 1276 (+ (* size (log 2e0)) 1277 (log (complex (float (/ re scale)) 1278 (float (/ im scale)))))))))) 1279 (t 1280 ;; log(n1/d1 + n2/d2*%i) = 1281 ;; log(s*(n+m*%i)) = log(s) + log(n+m*%i) 1282 ;; 1283 ;; where s = lcm(d1, d2), n and m are integers 1284 ;; 1285 (let* ((s (lcm ($denom ($realpart n)) 1286 ($denom ($imagpart n)))) 1287 (p ($expand (mul s n)))) 1288 (add ($float `((%log) ,s)) 1289 ($float `((%log) ,p)))))))) 1290 ((and (eq (caar e) '%erf) 1291 (eq (second e) '$%i)) 1292 ;; Handle like erf(%i). float(%i) (via recur-apply below) 1293 ;; just returns %i, so we never numerically evaluate it. 1294 (complexify (complex-erf (complex 0 1d0)))) 1295 ((or (eq (caar e) '%derivative) (eq (caar e) '$diff)) 1296 (append (list (remove 'simp (first e)) ($float (second e))) (rest (rest e)))) 1297 (t (recur-apply #'$float e)))) 1298 1299(defmfun $coeff (e x &optional (n 1)) 1300 (if (equal n 0) 1301 (coeff e x 0) 1302 (coeff e (power x n) 1))) 1303 1304(defun coeff (e var pow) 1305 (simplify 1306 (cond ((alike1 e var) (if (equal pow 1) 1 0)) 1307 ((atom e) (if (equal pow 0) e 0)) 1308 ((eq (caar e) 'mexpt) 1309 (cond ((alike1 (cadr e) var) 1310 (if (or (equal pow 0) (not (alike1 (caddr e) pow))) 0 1)) 1311 ((equal pow 0) e) 1312 (t 0))) 1313 ((or (eq (caar e) 'mplus) (mbagp e)) 1314 (cons (if (eq (caar e) 'mplus) '(mplus) (car e)) 1315 (mapcar #'(lambda (e) (coeff e var pow)) (cdr e)))) 1316 ((eq (caar e) 'mrat) (ratcoeff e var pow)) 1317 ((equal pow 0) (if (coeff-contains-powers e var) 0 e)) 1318 ((eq (caar e) 'mtimes) 1319 (let ((term (if (equal pow 1) var (power var pow)))) 1320 (if (memalike term (cdr e)) ($delete term e 1) 0))) 1321 (t 0)))) 1322 1323(defun coeff-contains-powers (e var) 1324 (cond ((alike1 e var) t) 1325 ((atom e) nil) 1326 ((eq (caar e) 'mexpt) 1327 (alike1 (cadr e) var)) 1328 ((eq (caar e) 'mtimes) 1329 (member t (mapcar #'(lambda (e) 1330 (coeff-contains-powers e var)) (cdr e)))) 1331 (t nil))) 1332 1333(let (my-powers my-num my-flag) 1334 (declare (special my-powers my-num my-flag)) 1335 1336 (defmfun $hipow (e var) 1337 (findpowers e t var)) 1338 1339 ;; These work best on expanded "simple" expressions. 1340 1341 (defmfun $lopow (e var) 1342 (findpowers e nil var)) 1343 1344 (defun findpowers (e hiflg var) 1345 (let (my-powers my-num my-flag) 1346 (declare (special my-powers my-num my-flag)) 1347 (findpowers1 e hiflg var) 1348 (cond ((null my-powers) (if (null my-num) 0 my-num)) 1349 (t (when my-num (setq my-powers (cons my-num my-powers))) 1350 (maximin my-powers (if hiflg '$max '$min)))))) 1351 1352 (defun findpowers1 (e hiflg var) 1353 (cond ((alike1 e var) (checkpow 1 hiflg)) 1354 ((atom e)) 1355 ((eq (caar e) 'mplus) 1356 (cond ((not (freel (cdr e) var)) 1357 (do ((e (cdr e) (cdr e))) ((null e)) 1358 (setq my-flag nil) (findpowers1 (car e) hiflg var) 1359 (if (null my-flag) (checkpow 0 hiflg)))))) 1360 ((and (eq (caar e) 'mexpt) (alike1 (cadr e) var)) (checkpow (caddr e) hiflg)) 1361 ((specrepp e) (findpowers1 (specdisrep e) hiflg var)) 1362 (t (mapc #'(lambda (x) (findpowers1 x hiflg var)) (cdr e))))) 1363 1364 (defun checkpow (pow hiflg) 1365 (setq my-flag t) 1366 (cond ((not (numberp pow)) (setq my-powers (cons pow my-powers))) 1367 ((null my-num) (setq my-num pow)) 1368 (hiflg (if (> pow my-num) (setq my-num pow))) 1369 ((< pow my-num) (setq my-num pow))))) 1370 1371;; push(x,l): The argument l must be a mapatom that is bound to a list. Prepend x to 1372;; the value of l and return this list. The arguments to push are evaluated from left to right. 1373 1374;; We assume that if each member of a list is simplified, the list members are still simplified 1375;; after pushing onto (or popping) the list. Just before returning, the code calls simplifya 1376;; with second argument true. Without this call to simplifya, the general simplifier would simplify 1377;; every list member after returning. (Barton Willis, author of $push and $pop) 1378 1379(defmspec $push (z) 1380 (let* ((o (car (pop z))) 1381 (x (if z (pop z) (wna-err o))) 1382 (l (if z (pop z) (wna-err o))) 1383 (ll) (lo)) 1384 1385 (if z (wna-err o)) 1386 (setq x (meval x)) 1387 (cond (($subvarp l) 1388 (setq lo (simplifya (cons (list (meval (mop l)) 'array) (mevalargs (margs l))) t)) 1389 (setq ll (let ((noevalargs t)) (meval lo)))) 1390 (t 1391 (setq lo l) 1392 (setq ll (meval l)))) 1393 (if (and ($mapatom lo) ($listp ll)) 1394 (progn 1395 (setq ll (cons x (cdr ll))) 1396 (mset lo (simplifya (cons '(mlist) ll) t))) 1397 (merror "Second argument to push must be a mapatom that is bound to a list")))) 1398 1399;; pop(l): The argument l must be a mapatom that is bound to a list. Remove the 1400;; first member of the value of l and return it. 1401 1402;; We assume that if each member of a list is simplified, the list members are still simplified 1403;; after popping (or pushing onto) the list. Just before returning, the code calls simplifya 1404;; with second argument true. Without this call to simplifya, the general simplifier would simplify 1405;; every list member after returning. 1406 1407(defmspec $pop (z) 1408 (let* ((o (car (pop z))) 1409 (l (if z (pop z) (wna-err o))) 1410 (ll) (lo)) 1411 (if z (wna-err o)) 1412 (cond (($subvarp l) 1413 (setq lo (simplifya (cons (list (meval (mop l)) 'array) (mevalargs (margs l))) t)) 1414 (setq ll (let ((noevalargs t)) (meval lo)))) 1415 (t 1416 (setq lo l) 1417 (setq ll (meval lo)))) 1418 (setq ll (meval l)) 1419 (cond ((and ($mapatom lo) ($listp ll)) 1420 (if ($emptyp ll) (merror "Pop called on an empty list") 1421 (prog2 1422 (setq ll (cdr ll)) 1423 (pop ll) 1424 (mset lo (simplifya (cons '(mlist) ll) t))))) 1425 (t (merror "Argument to pop must be a mapatom that is bound to a list"))))) 1426 1427