1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 ;; ** (c) Copyright 1982 Massachusetts Institute of Technology **
15 (declare-top (special *linelabel
*))
17 ;; op and opr properties
19 (defvar *opr-table
* (make-hash-table :test
#'equal
))
23 (and (symbolp x
) (get x
'opr
))
24 (and (stringp x
) (gethash x
*opr-table
*))))
31 (and (symbolp x
) (setf (get x
'opr
) y
))
32 (and (stringp x
) (setf (gethash x
*opr-table
*) y
))))
36 (and (symbolp x
) (remprop x
'opr
))
37 (and (stringp x
) (remhash x
*opr-table
*))))
39 ;; This business about operator properties is terrible --
40 ;; this stuff should be in src/nparse.lisp, and it should be split up
41 ;; for each operator. Extra points for making it declarative.
44 (putprop (car x
) (cadr x
) 'op
)
45 (putprop (intern (concatenate 'string
"%" (symbol-name (car x
)))) (cadr x
) 'op
) ;; nounify not yet available
46 (putopr (cadr x
) (car x
))
47 (push (cadr x
) *mopl
*))
49 '((mplus "+") (mminus "-") (mtimes "*") (mexpt "**") (mexpt "^")
50 (mnctimes ".") (rat "/") (mquotient "/") (mncexpt "^^")
51 (mequal "=") (mgreaterp ">") (mlessp "<") (mleqp "<=") (mgeqp ">=")
52 (mnotequal "#") (mand "and") (mor "or") (mnot "not") (msetq ":")
53 (mdefine ":=") (mdefmacro "::=") (mquote "'") (mlist "[")
54 (mset "::") (mfactorial "!") (marrow "-->") (mprogn "(")
55 (mcond "if") (mdo "do") (mdoin "do_in")))
57 (mapc #'(lambda (x) (putprop (car x
) (cadr x
) 'op
))
58 '((mqapply $subvar
) (bigfloat $bfloat
)))
60 (defmvar $vect_cross nil
61 "If TRUE allows DIFF(X~Y,T) to work where ~ is defined in
62 SHARE;VECT where VECT_CROSS is set to TRUE.")
64 (defmfun $substitute
(new old
&optional
(expr nil three-arg?
))
65 (cond (three-arg?
(maxima-substitute new old expr
))
67 (let ((l new
) (z old
))
68 (cond ((and ($listp l
) ($listp
(cadr l
)) (null (cddr l
)))
69 ($substitute
(cadr l
) z
))
70 ((notloreq l
) (improper-arg-err l
'$substitute
))
71 ((eq (caar l
) 'mequal
) (maxima-substitute (caddr l
) (cadr l
) z
))
72 (t (do ((l (cdr l
) (cdr l
)))
74 (setq z
($substitute
(car l
) z
)))))))))
76 ;; Define an alias $psubst and a reversealias for $psubstitute
77 (defprop $psubst $psubstitute alias
)
78 (defprop $psubstitute $psubst reversealias
)
80 ;; $psubstitute is similar to $substitute. In distinction from $substitute
81 ;; the function $psubstitute does parallel substitution, if the first argument
82 ;; is a list of equations.
83 (defmfun $psubstitute
(old new
&optional
(expr nil three-arg?
))
84 (cond (three-arg?
(maxima-substitute old new expr
))
86 (let ((l old
) (z new
))
87 (cond ((and ($listp l
)
91 ($psubstitute
(cadr l
) z
))
93 (eq (caar (cadr l
)) 'mequal
)
95 ;; A list with one equation.
96 ($psubstitute
(cadr l
) z
))
97 ((notloreq l
) (improper-arg-err l
'$psubstitute
))
98 ((eq (caar l
) 'mequal
)
99 ;; Do a substitution for one equation.
100 (maxima-substitute (caddr l
) (cadr l
) z
))
102 ;; We have a list of equations. We do parallel substitution.
103 (let (gensymbol genlist eqn
($simp nil
))
104 ;; At first substitute a gensym for the expressions of
105 ;; the left hand side of the equations.
106 (do ((l (cdr l
) (cdr l
)))
109 (when (not (eq 'mequal
(caar eqn
)))
110 (improper-arg-err old
'$substitute
))
111 (setq gensymbol
(gensym))
112 ;; Store the gensym and the new expression into a list.
113 (push (cons gensymbol
(caddr eqn
)) genlist
)
114 ;; Substitute a gensym for the old expression.
115 (setq z
(maxima-substitute gensymbol
(cadr eqn
) z
)))
116 ;; Substitute the new expressions for the gensyms.
117 (do ((l genlist
(cdr l
)))
119 ;; Resimplify the result.
120 (let (($simp t
)) (resimplify z
)))
121 (setq z
(maxima-substitute (cdar l
) (caar l
) z
))))))))))
123 (defun maxima-substitute (x y z
) ; The args to SUBSTITUTE are assumed to be simplified.
124 ;; Prevent replacing dependent variable with constant in derivative
125 (cond ((and (not (atom z
))
126 (eq (caar z
) '%derivative
)
131 (let ((in-p t
) (substp t
))
132 (if (and (mnump y
) (= (signum1 y
) 1))
133 (let ($sqrtdispflag
($pfeformat t
)) (setq z
(nformat-all z
))))
137 (setq y
'((mminus) 1))
138 (subst2 x y
(nformat-all z
) nil nil
)) ;; negxpty and timesp don't matter in this call since (caar y) != 'mexpt
140 (cond ((and (not (symbolp x
))
142 (let ((tem (gensym)))
143 (setf (get tem
'operators
) 'application-operator
)
144 (setf (symbol-function tem
) x
)
147 (let ((negxpty (if (and (eq (caar y
) 'mexpt
)
148 (= (signum1 (caddr y
)) 1))
149 (mul2 -
1 (caddr y
))))
150 (timesp (if (eq (caar y
) 'mtimes
)
151 (setq y
(nformat y
)))))
152 (subst2 x y z negxpty timesp
)))
155 ;;Remainder of page is update from F302 --gsb
157 (defun subst1 (x y z
) ; Y is an atom
158 (cond ((atom z
) (if (equal y z
) x z
))
159 ((specrepp z
) (subst1 x y
(specdisrep z
)))
160 ((eq (caar z
) 'bigfloat
) z
)
161 ((and (eq (caar z
) 'rat
) (or (equal y
(cadr z
)) (equal y
(caddr z
))))
162 (div (subst1 x y
(cadr z
)) (subst1 x y
(caddr z
))))
163 ((at-substp z
) (subst-except-second-arg x y z
))
164 ((and (eq y t
) (eq (caar z
) 'mcond
))
165 (list (cons (caar z
) nil
) (subst1 x y
(cadr z
)) (subst1 x y
(caddr z
))
166 (cadddr z
) (subst1 x y
(car (cddddr z
)))))
167 (t (let ((margs (mapcar #'(lambda (z1) (subst1 x y z1
)) (cdr z
)))
168 (oprx (getopr x
)) (opry (getopr y
)))
170 (or (eq opry
(caar z
))
171 (and (eq (caar z
) 'rat
) (eq opry
'mquotient
))))
173 (member x
'(t nil $%e $%pi $%i
) :test
#'eq
)
175 (not (or (eq (car x
) 'lambda
)
176 (eq (caar x
) 'lambda
)))))
177 (if (or (and (member 'array
(cdar z
) :test
#'eq
)
178 (or (and (mnump x
) $subnumsimp
)
179 (and (not (mnump x
)) (not (atom x
)))))
181 (let ((substp 'mqapply
))
182 (subst0 (list* '(mqapply) x margs
) z
))
183 (merror (intl:gettext
"subst: cannot substitute ~M for operator ~M in expression ~M") x y z
))
184 (subst0 (cons (cons oprx nil
) margs
) z
))
185 (subst0 (cons (cons (caar z
) nil
) margs
) z
))))))
187 (defun subst2 (x y z negxpty timesp
)
190 ((specrepp z
) (subst2 x y
(specdisrep z
) negxpty timesp
))
191 ((at-substp z
) z
) ;; IS SUBST-EXCEPT-SECOND-ARG APPROPRIATE HERE ?? !!
193 ((and timesp
(eq (caar z
) 'mtimes
) (alike1 y
(setq z
(nformat z
)))) x
)
194 ((and (eq (caar y
) 'mexpt
) (eq (caar z
) 'mexpt
) (alike1 (cadr y
) (cadr z
))
195 (setq newexpt
(cond ((alike1 negxpty
(caddr z
)) -
1)
196 ($exptsubst
(expthack (caddr y
) (caddr z
))))))
197 (list '(mexpt) x newexpt
))
198 ((and $derivsubst
(eq (caar y
) '%derivative
) (eq (caar z
) '%derivative
)
199 (alike1 (cadr y
) (cadr z
)))
200 (let ((tail (subst-diff-match (cddr y
) (cdr z
))))
201 (cond ((null tail
) z
)
202 (t (cons (cons (caar z
) nil
) (cons x
(cdr tail
)))))))
203 (t (recur-apply #'(lambda (z1) (subst2 x y z1 negxpty timesp
)) z
)))))
205 ;; replace y with x in z, but leave z's second arg unchanged.
206 ;; This is for cases like at(integrate(x, x, a, b), [x=3])
207 ;; where second arg of integrate binds a new variable x,
208 ;; and we do not wish to subst 3 for x inside integrand.
209 (defun subst-except-second-arg (x y z
)
211 ((member (caar z
) '(%integrate %sum %product %limit %laplace
))
213 (list (remove 'simp
(car z
)) ; ensure resimplification after substitution
214 (if (eq y
(third z
)) ; if (third z) is new var that shadows y
215 (second z
) ; leave (second z) unchanged
216 (subst1 x y
(second z
))) ; otherwise replace y with x in (second z)
217 (third z
)) ; never change integration var
218 (mapcar (lambda (z) (subst1 x y z
)); do subst in limits of integral
221 ;; similar considerations here, but different structure of expression.
223 ((at-eqn-or-eqns (third z
))
224 (at-eqns-list (if (eq (caar at-eqn-or-eqns
) 'mlist
) (rest at-eqn-or-eqns
) (list at-eqn-or-eqns
))))
226 (remove 'simp
(car z
)) ;; ensure resimplification after substitution
227 (if (member y
(mapcar #'(lambda (e) (second e
)) at-eqns-list
))
229 (subst1 x y
(second z
)))
230 `((mlist) ,@(mapcar #'(lambda (e) (list (first e
) (second e
) (subst1 x y
(third e
)))) at-eqns-list
)))))
231 ((eq (caar z
) '%derivative
)
232 ;; and again, with yet a different structure.
234 ((vars-and-degrees (rest (rest z
)))
235 (diff-vars (loop for v in vars-and-degrees by
#'cddr collect v
))
236 (diff-degrees (loop for n in
(rest vars-and-degrees
) by
#'cddr collect n
)))
239 (remove 'simp
(car z
)) ;; ensure resimplification after substitution
240 (if (member y diff-vars
)
242 (subst1 x y
(second z
))))
243 (apply #'append
(loop for v in diff-vars for n in diff-degrees collect
(list v
(subst1 x y n
)))))))
246 (defun subst0 (new old
)
247 (cond ((atom new
) new
)
248 ((alike (cdr new
) (cdr old
))
249 (cond ((eq (caar new
) (caar old
)) old
)
250 (t (simplifya (cons (cons (caar new
) (member 'array
(cdar old
) :test
#'eq
)) (cdr old
))
252 ((member 'array
(cdar old
) :test
#'eq
)
253 (simplifya (cons (cons (caar new
) '(array)) (cdr new
)) nil
))
254 (t (simplifya new nil
))))
256 (defun expthack (y z
)
257 (prog (nn* dn
* yn yd zn zd qd
)
258 (cond ((and (mnump y
) (mnump z
))
259 (return (if (numberp (setq y
(div* z y
))) y
)))
260 ((atom z
) (if (not (mnump y
)) (return nil
)))
261 ((or (ratnump z
) (eq (caar z
) 'mplus
)) (return nil
)))
262 (numden y
) ; (CSIMP) sets NN* and DN*
266 (setq qd
(cond ((and (equal zd
1) (equal yd
1)) 1)
267 ((prog2 (numden (div* zd yd
))
268 (and (equal dn
* 1) (equal nn
* 1)))
270 ((equal nn
* 1) (div* 1 dn
*))
273 (numden (div* zn yn
))
274 (if (equal dn
* 1) (return (div* nn
* qd
)))))
276 (defun subst-diff-match (l1 l2
)
277 (do ((l l1
(cddr l
)) (l2 (copy-list l2
)) (failed nil nil
))
279 (do ((l2 l2
(cddr l2
)))
280 ((null (cdr l2
)) (setq failed t
))
281 (if (alike1 (car l
) (cadr l2
))
282 (if (and (fixnump (cadr l
)) (fixnump (caddr l2
)))
283 (cond ((< (cadr l
) (caddr l2
))
284 (return (rplacd (cdr l2
)
285 (cons (- (caddr l2
) (cadr l
))
287 ((= (cadr l
) (caddr l2
))
288 (return (rplacd l2
(cdddr l2
))))
289 (t (return (setq failed t
))))
290 (return (setq failed t
)))))
291 (if failed
(return nil
))))
293 ;;This probably should be a subst or macro.
296 (or (member (caar z
) '(%derivative %del
) :test
#'eq
)
297 (member (caar z
) dummy-variable-operators
:test
#'eq
))))
299 (defun recur-apply (fun e
)
300 (cond ((eq (caar e
) 'bigfloat
) e
)
301 ((specrepp e
) (funcall fun
(specdisrep e
)))
302 (t (let ((newargs (mapcar fun
(cdr e
))))
303 (if (alike newargs
(cdr e
))
305 (simplifya (cons (cons (caar e
) (member 'array
(cdar e
) :test
#'eq
)) newargs
)
308 (defmfun $depends
(&rest args
)
309 (when (oddp (length args
))
310 (merror (intl:gettext
"depends: number of arguments must be even.")))
311 (do ((args args
(cddr args
))
313 ((null args
) (i-$dependencies
(nreverse l
)))
314 (if ($listp
(first args
))
315 (mapc #'(lambda (e) (push (depends1 e
(second args
)) l
))
317 (push (depends1 (first args
) (second args
)) l
))))
319 (defun depends1 (x y
)
320 (nonsymchk x
'$depends
)
321 (cons (cons x nil
) (if ($listp y
) (cdr y
) (cons y nil
))))
323 (defmspec $dependencies
(form)
324 (i-$dependencies
(cdr form
)))
326 (defun i-$dependencies
(l &aux res
)
332 "depends: argument must be a non-atomic expression; found ~M") z
))
337 (mputprop (caar z
) (setq y
(reverse y
)) 'depends
)
338 (setq res
(push (cons (ncons (caar z
)) y
) res
))
339 (unless (cdr $dependencies
)
340 (setq $dependencies
(copy-list '((mlist simp
)))))
341 (add2lnc (cons (cons (caar z
) nil
) y
) $dependencies
))
342 (cond ((and ($subvarp
(cadr zz
))
343 (not (member (caar (cadr zz
)) y
)))
344 (setq y
(push (cadr zz
) y
)))
345 ((not (symbolp (cadr zz
)))
347 (intl:gettext
"depends: argument must be a symbol; found ~M")
350 (not (member (cadr zz
) y
)))
351 (setq y
(push (cadr zz
) y
))))))))
352 (cons '(mlist simp
) (reverse res
)))
354 (defmspec $gradef
(l)
356 (let ((z (car l
)) (n 0))
358 (if (not (= (length l
) 3)) (merror (intl:gettext
"gradef: expected exactly three arguments.")))
360 (cons (cons (cadr l
) (meval (caddr l
)))
363 (i-$dependencies
(cons (cons (ncons z
)
364 ;; Append existing dependencies
365 (cons (cadr l
) (mget z
'depends
)))
369 ((or (mopp1 (caar z
)) (member 'array
(cdar z
) :test
#'eq
))
370 (merror (intl:gettext
"gradef: argument cannot be a built-in operator or subscripted expression; found ~M") z
))
371 ((prog2 (setq n
(- (length z
) (length l
))) (minusp n
))
373 (t (do ((zl (cdr z
) (cdr zl
))) ((null zl
))
374 (if (not (symbolp (car zl
)))
375 (merror (intl:gettext
"gradef: argument must be a symbol; found ~M") (car zl
))))
376 (setq l
(nconc (mapcar #'(lambda (x) (remsimp (meval x
)))
378 (mapcar #'(lambda (x) (list '(%derivative
) z x
1))
379 (nthcdr (- (length z
) n
) z
))))
381 (sublis (mapcar #'cons
(cdr z
) (mapcar #'stripdollar
(cdr z
)))
384 (or (cdr $gradefs
) (setq $gradefs
(copy-list '((mlist simp
)))))
385 (add2lnc (cons (cons (caar z
) nil
) (cdr z
)) $gradefs
) z
))))
387 (defmfun $diff
(&rest args
)
388 (declare (dynamic-extent args
))
397 (cond ((null e
) (wna-err '$diff
))
398 ((null (cdr e
)) (return (stotaldiff (car e
))))
399 ((null (cddr e
)) (nconc e
'(1))))
400 (setq exp
(car e
) z
(setq e
(copy-list e
)))
401 loop
(if (or (null derivlist
) (member (cadr z
) derivlist
:test
#'equal
)) (go doit
))
402 ; DERIVLIST is set by $EV
404 loop2
(cond ((cdr z
) (go loop
))
405 ((null (cdr e
)) (return exp
))
407 doit
(cond ((nonvarcheck (cadr z
) '$diff
))
408 ((null (cddr z
)) (wna-err '$diff
))
409 ((not (fixnump (caddr z
))) (go noun
))
410 ((minusp (setq count
(caddr z
)))
411 (merror (intl:gettext
"diff: order of derivative must be a nonnegative integer; found ~M") count
)))
412 loop1
(cond ((zerop count
) (rplacd z
(cdddr z
)) (go loop2
))
413 ((equal (setq exp
(sdiff exp
(cadr z
))) 0) (return 0)))
414 (setq count
(1- count
))
416 noun
(return (diff%deriv
(cons exp
(cdr e
))))))
418 (defun chainrule (e x
)
421 ;; sdiff is called from the function islinear.
422 (if (and (not (atom e
))
423 (eq (caar e
) '%derivative
)
424 (not (freel (cdr e
) x
)))
425 (diff%deriv
(list e x
1))
428 ((not (setq w
(mget (cond ((atom e
) e
)
429 ((member 'array
(cdar e
) :test
#'eq
) (caar e
))
430 ((atom (cadr e
)) (cadr e
))
437 (let ((y (sdiff u x
)))
442 (list '(%derivative
) e u
1))
447 (defun atomgrad (e x
)
449 (and (atom e
) (setq y
(mget e
'$atomgrad
)) (assolike x y
))))
451 (defun depends (e x
&aux l
)
452 (setq e
(specrepcheck e
))
453 (cond ((alike1 e x
) t
)
455 ((and (symbolp e
) (setq l
(mget e
'depends
)))
456 ;; Go recursively through the list of dependencies.
457 ;; This code detects indirect dependencies like a(x) and x(t).
460 (t (or (depends (caar e
) x
)
461 (dependsl (cdr e
) x
)))))
463 (defun dependsl (l x
)
465 (if (depends u x
) (return t
))))
467 (defun sdiff (e x
) ; The args to SDIFF are assumed to be simplified.
468 ;; Remove a special representation from the variable of differentiation
469 (setq x
(specrepcheck x
))
470 (cond ((alike1 e x
) 1)
472 ((or (atom e
) (member 'array
(cdar e
) :test
#'eq
)) (chainrule e x
))
473 ((eq (caar e
) 'mrat
) (ratdx e x
))
474 ((eq (caar e
) 'mpois
) ($poisdiff e x
)) ; Poisson series
475 ((eq (caar e
) 'mplus
) (addn (sdiffmap (cdr e
) x
) t
))
476 ((mbagp e
) (cons (car e
) (sdiffmap (cdr e
) x
)))
477 ((member (caar e
) '(%sum %product
) :test
#'eq
) (diffsumprod e x
))
478 ((eq (caar e
) '%at
) (diff-%at e x
))
479 ((not (depends e x
)) 0)
480 ((eq (caar e
) 'mtimes
) (addn (sdifftimes (cdr e
) x
) t
))
481 ((eq (caar e
) 'mexpt
) (diffexpt e x
))
482 ((eq (caar e
) 'mnctimes
)
483 (let (($dotdistrib t
))
484 (add2 (ncmuln (cons (sdiff (cadr e
) x
) (cddr e
)) t
)
485 (ncmul2 (cadr e
) (sdiff (cons '(mnctimes) (cddr e
)) x
)))))
486 ((and $vect_cross
(eq (caar e
) '|$~|
))
487 (add2* `((|$~|
) ,(cadr e
) ,(sdiff (caddr e
) x
))
488 `((|$~|
) ,(sdiff (cadr e
) x
) ,(caddr e
))))
489 ((eq (caar e
) 'mncexpt
) (diffncexpt e x
))
490 ((member (caar e
) '(%log %plog
) :test
#'eq
)
491 (sdiffgrad (cond ((and (not (atom (cadr e
))) (eq (caaadr e
) 'mabs
))
492 (cons (car e
) (cdadr e
)))
495 ((eq (caar e
) '%derivative
)
496 (cond ((or (atom (cadr e
)) (member 'array
(cdaadr e
) :test
#'eq
)) (chainrule e x
))
497 ((freel (cddr e
) x
) (diff%deriv
(cons (sdiff (cadr e
) x
) (cddr e
))))
498 (t (diff%deriv
(list e x
1)))))
499 ((member (caar e
) '(%binomial %beta
) :test
#'eq
)
500 (let ((efact ($makefact e
)))
501 (mul2 (factor (sdiff efact x
)) (div e efact
))))
502 ((eq (caar e
) '%integrate
) (diffint e x
))
503 ((eq (caar e
) '%laplace
) (difflaplace e x
))
504 ((eq (caar e
) '%at
) (diff-%at e x
))
505 ; This rule is not correct. We cut it out.
506 ; ((member (caar e) '(%realpart %imagpart) :test #'eq)
507 ; (list (cons (caar e) nil) (sdiff (cadr e) x)))
508 ((and (eq (caar e
) 'mqapply
)
509 (eq (caaadr e
) '$%f
))
510 ;; Handle %f, hypergeometric function
512 ;; The derivative of %f[p,q]([a1,...,ap],[b1,...,bq],z) is
514 ;; a1*a2*...*ap/(b1*b2*...*bq)
515 ;; *%f[p,q]([a1+1,a2+1,...,ap+1],[b1+1,b2+1,...,bq+1],z)
516 (let* ((arg1 (cdr (third e
)))
517 (arg2 (cdr (fourth e
)))
520 (div (mull arg1
) (mull arg2
))
521 `((mqapply) (($%f array
) ,(length arg1
) ,(length arg2
))
522 ((mlist) ,@(incr1 arg1
))
523 ((mlist) ,@(incr1 arg2
))
525 (t (sdiffgrad e x
))))
527 (defun sdiffgrad (e x
)
528 (let ((fun (caar e
)) grad args result
)
529 (cond ((and (eq fun
'mqapply
) (zl-get (caaadr e
) 'grad
))
530 ;; Change the array function f[n](x) to f(n,x), call sdiffgrad again.
532 (sdiffgrad (cons (cons (caaadr e
) nil
)
533 (append (cdadr e
) (cddr e
)))
535 ;; If noun form for f(n,x), adjust the noun form for f[n](x)
536 (if (isinop result
'%derivative
)
537 (if (not (depends e x
))
539 (diff%deriv
(list e x
1)))
542 ;; extension for pdiff.
543 ((and (get '$pderivop
'operators
) (funcall 'sdiffgrad-pdiff e x
)))
545 ;; two line extension for hypergeometric.
546 ((and (equal fun
'%hypergeometric
) (get '%hypergeometric
'operators
))
547 (funcall 'diff-hypergeometric
(second e
) (third e
) (fourth e
) x
))
549 ((or (eq fun
'mqapply
) (null (setq grad
(zl-get fun
'grad
))))
550 (if (not (depends e x
)) 0 (diff%deriv
(list e x
1))))
551 ((not (= (length (cdr e
)) (length (car grad
))))
552 (merror (intl:gettext
"~:M: expected exactly ~M arguments.")
554 (length (car grad
))))
556 (setq args
(sdiffmap (cdr e
) x
))
565 (do ((l1 (cdr grad
) (cdr l1
))
566 (args args
(cdr args
))
568 ((null l1
) (cons '(mlist) (nreverse l2
)))
570 (cons (cond ((equal (car args
) 0) 0)
571 ((functionp (car l1
))
572 ;; Evaluate a lambda expression
573 ;; given as a derivative.
574 (apply (car l1
) (cdr e
)))
579 (if (or (null result
) (not (freeof nil result
)))
580 ;; A derivative has returned NIL. Return a noun form.
581 (if (not (depends e x
))
583 (diff%deriv
(list e x
1)))
586 (defun sdiffmap (e x
)
587 (mapcar #'(lambda (term) (sdiff term x
)) e
))
589 (defun sdifftimes (l x
)
590 (prog (term left out
)
591 loop
(setq term
(car l
) l
(cdr l
))
592 (setq out
(cons (muln (cons (sdiff term x
) (append left l
)) t
) out
))
593 (if (null l
) (return out
))
594 (setq left
(cons term left
))
597 (defun diffexpt (e x
)
598 (if (mnump (caddr e
))
599 (mul3 (caddr e
) (power (cadr e
) (addk (caddr e
) -
1)) (sdiff (cadr e
) x
))
600 (mul2 e
(add2 (mul3 (power (cadr e
) -
1) (caddr e
) (sdiff (cadr e
) x
))
601 (mul2 (simplifya (list '(%log
) (cadr e
)) t
)
602 (sdiff (caddr e
) x
))))))
604 (defun diff%deriv
(e)
606 (simplifya (cons '(%derivative
) e
) t
)))
612 (mapc #'(lambda (z) (putprop (car z
) (cons header
(cdr z
)) 'grad
))
613 ;; All these GRAD templates have been simplified and then the SIMP flags
614 ;; (which are unnecessary) have been removed to save core space.
615 '((%log
((mexpt) x -
1)) (%plog
((mexpt) x -
1))
616 (%gamma
((mtimes) ((mqapply) (($psi array
) 0) x
) ((%gamma
) x
)))
617 (mfactorial ((mtimes) ((mqapply) (($psi array
) 0) ((mplus) 1 x
)) ((mfactorial) x
)))
619 (%cos
((mtimes) -
1 ((%sin
) x
)))
620 (%tan
((mexpt) ((%sec
) x
) 2))
621 (%cot
((mtimes) -
1 ((mexpt) ((%csc
) x
) 2)))
622 (%sec
((mtimes) ((%sec
) x
) ((%tan
) x
)))
623 (%csc
((mtimes) -
1 ((%cot
) x
) ((%csc
) x
)))
624 (%asin
((mexpt) ((mplus) 1 ((mtimes) -
1 ((mexpt) x
2))) ((rat) -
1 2)))
625 (%acos
((mtimes) -
1 ((mexpt) ((mplus) 1 ((mtimes) -
1 ((mexpt) x
2)))
627 (%atan
((mexpt) ((mplus) 1 ((mexpt) x
2)) -
1))
628 (%acot
((mtimes) -
1 ((mexpt) ((mplus) 1 ((mexpt) x
2)) -
1)))
630 ((mexpt) ((mplus) 1 ((mtimes) -
1 ((mexpt) x -
2)))
634 ((mexpt) ((mplus) 1 ((mtimes) -
1 ((mexpt) x -
2)))
639 (%tanh
((mexpt) ((%sech
) x
) 2))
640 (%coth
((mtimes) -
1 ((mexpt) ((%csch
) x
) 2)))
641 (%sech
((mtimes) -
1 ((%sech
) x
) ((%tanh
) x
)))
642 (%csch
((mtimes) -
1 ((%coth
) x
) ((%csch
) x
)))
643 (%asinh
((mexpt) ((mplus) 1 ((mexpt) x
2)) ((rat) -
1 2)))
644 (%acosh
((mexpt) ((mplus) -
1 ((mexpt) x
2)) ((rat) -
1 2)))
645 (%atanh
((mexpt) ((mplus) 1 ((mtimes) -
1 ((mexpt) x
2))) -
1))
646 (%acoth
((mtimes) -
1 ((mexpt) ((mplus) -
1 ((mexpt) x
2)) -
1)))
648 ((mexpt) ((mplus) -
1 ((mexpt) x -
2)) ((rat) -
1 2))
651 ((mexpt) ((mplus) 1 ((mexpt) x -
2)) ((rat) -
1 2))
653 (mabs ((mtimes) x
((mexpt) ((mabs) x
) -
1)))
654 (%erf
((mtimes) 2 ((mexpt) $%pi
((rat) -
1 2))
655 ((mexpt) $%e
((mtimes) -
1 ((mexpt) x
2)))))
658 (defprop %atan2
((x y
) ((mtimes) y
((mexpt) ((mplus) ((mexpt) x
2) ((mexpt) y
2)) -
1))
659 ((mtimes) -
1 x
((mexpt) ((mplus) ((mexpt) x
2) ((mexpt) y
2)) -
1)))
664 ; Do not put a noun form on the property list, but NIL.
665 ; SDIFFGRAD generates the noun form.
666 ; ((%derivative) ((mqapply) (($li array) n) x) n 1)
668 ((mtimes) ((mqapply) (($li array
) ((mplus) -
1 n
)) x
) ((mexpt) x -
1)))
673 ; Do not put a noun form on the property list, but NIL.
674 ; SDIFFGRAD generates the noun form.
676 ((mqapply) (($psi array
) ((mplus) 1 n
)) x
))
679 (defun atvarschk (argl)
680 (do ((largl (length argl
) (1- largl
))
681 (latvrs (length atvars
))
683 ((not (< latvrs largl
)) (nconc atvars l
))
684 (setq l
(cons (implode (cons '$
(cons '@ (mexploden largl
)))) l
))))
688 (not (member (caar x
) '(mlist mequal
) :test
#'eq
))
689 (and (eq (caar x
) 'mlist
)
691 (if (not (mequalp u
)) (return t
))))))
693 (defun substitutel (l1 l2 e
)
694 "l1 is a list of expressions. l2 is a list of variables. For each
695 element in list l2, substitute corresponding element of l1 into e"
696 (do ((l1 l1
(cdr l1
))
699 (setq e
(maxima-substitute (car l1
) (car l2
) e
))))
705 (if (not (memalike (car a
) b
)) (setq x
(cons (car a
) x
)))))
707 (defun intersect* (a b
)
711 (if (memalike (car a
) b
) (setq x
(cons (car a
) x
)))))
714 (car (nthcdr (1- n
) e
)))
716 (defun delsimp (e) (remove 'simp e
))
719 (if (atom e
) e
(cons (delsimp (car e
)) (mapcar #'remsimp
(cdr e
)))))
723 ((eq (caar e
) 'mplus
) (cons (append (car e
) '(trunc)) (cdr e
)))
724 ((mbagp e
) (cons (car e
) (mapcar #'$trunc
(cdr e
))))
725 ((specrepp e
) ($trunc
(specdisrep e
)))
728 (defun nonvarcheck (e fn
)
732 (and (not (atom e
)) (not (eq (caar e
) 'mqapply
)) (mopp1 (caar e
))))
733 (merror (intl:gettext
"~:M: second argument must be a variable; found ~M") fn e
)))
735 (defmspec $ldisplay
(form)
736 (disp1 (cdr form
) t t
))
738 (defmfun $ldisp
(&rest args
)
739 (declare (dynamic-extent args
))
742 (defmspec $display
(form)
743 (disp1 (cdr form
) nil t
))
745 (defmfun $disp
(&rest args
)
746 (declare (dynamic-extent args
))
747 (disp1 args nil nil
))
749 (defun disp1 (ll lablist eqnsp
)
750 (if lablist
(setq lablist
(cons '(mlist simp
) nil
)))
751 (do ((ll ll
(cdr ll
))
756 ((null ll
) (or lablist
'$done
))
757 (setq l
(car ll
) ans
(if eqnsp
(meval l
) l
))
758 (if (and eqnsp
(not (mequalp ans
)))
759 (setq ans
(list '(mequal simp
) (disp2 l
) ans
)))
760 (if lablist
(nconc lablist
(cons (elabel ans
) nil
)))
761 (setq tim
(get-internal-run-time))
762 (let ((*display-labels-p
* (not (null lablist
))))
763 (displa (list '(mlabel) (if lablist
*linelabel
*) ans
)))
769 ((eq (caar e
) 'mqapply
)
770 (cons '(mqapply) (cons (cons (caadr e
) (mapcar #'meval
(cdadr e
)))
771 (mapcar #'meval
(cddr e
)))))
772 ((eq (caar e
) 'msetq
) (disp2 (cadr e
)))
773 ((eq (caar e
) 'mset
) (disp2 (meval (cadr e
))))
774 ((eq (caar e
) 'mlist
) (cons (car e
) (mapcar #'disp2
(cdr e
))))
775 ((mspecfunp (caar e
)) e
)
776 (t (cons (car e
) (mapcar #'meval
(cdr e
))))))
778 ; Construct a new intermediate result label,
779 ; and bind it to the expression e.
780 ; The global flag $NOLABELS is ignored; the label is always bound.
781 ; Otherwise (if ELABEL were to observe $NOLABELS) it would be
782 ; impossible to programmatically refer to intermediate result expression.
785 (if (not (checklabel $linechar
)) (setq $linenum
(1+ $linenum
)))
786 (let (($nolabels nil
)) ; <-- This is pretty ugly. MAKELABEL should take another argument.
787 (makelabel $linechar
))
788 (setf (symbol-value *linelabel
*) e
)
791 (defmfun $dispterms
(e)
792 (cond ((or (atom e
) (eq (caar e
) 'bigfloat
)) (displa e
))
793 ((specrepp e
) ($dispterms
(specdisrep e
)))
794 (t (let (($dispflag t
))
796 (displa (getop (mop e
)))
797 (do ((e (if (and (eq (caar e
) 'mplus
) (not $powerdisp
))
800 (cdr e
))) ((null e
)) (mterpri) (displa (car e
)) (mterpri)))
804 (defmfun $dispform
(e &optional
(flag nil flag?
))
805 (when (and flag?
(not (eq flag
'$all
)))
806 (merror (intl:gettext
"dispform: second argument, if present, must be 'all'; found ~M") flag
))
808 (atom (setq e
(if flag?
(nformat-all e
) (nformat e
))))
809 (member 'simp
(cdar e
) :test
#'eq
))
811 (cons (cons (caar e
) (cons 'simp
(cdar e
)))
812 (if (and (eq (caar e
) 'mplus
) (not $powerdisp
))
816 ;;; These functions implement the Macsyma functions $op and $operatorp.
819 (let (($partswitch nil
))
822 (defmfun $operatorp
(expr oplist
)
824 ($member
($op expr
) oplist
)
825 (alike1 ($op expr
) oplist
)))
827 (defmfun $part
(&rest args
)
828 (declare (dynamic-extent args
))
829 (mpart args nil nil $inflag
'$part
))
831 (defmfun $inpart
(&rest args
)
832 (declare (dynamic-extent args
))
833 (mpart args nil nil t
'$inpart
))
835 (defmspec $substpart
(l)
837 (mpart (cdr l
) t nil $inflag
'$substpart
)))
839 (defmspec $substinpart
(l)
841 (mpart (cdr l
) t nil t
'$substinpart
)))
843 (defun part1 (arglist substflag dispflag inflag fn
) ; called only by TRANSLATE
845 (mpart arglist substflag dispflag inflag fn
)))
847 (defun mpart (arglist substflag dispflag inflag fn
)
848 (prog (substitem arg arg1 exp exp1 exp
* sevlist count prevcount n specp
850 (setq specp
(or substflag dispflag
))
851 (if substflag
(setq substitem
(car arglist
) arglist
(cdr arglist
)))
852 (if (null arglist
) (wna-err fn
))
853 (setq exp
(if substflag
(meval (car arglist
)) (car arglist
)))
854 (when (null (setq arglist
(cdr arglist
)))
856 (return (cond (substflag (meval substitem
))
857 (dispflag (box exp dispflag
))
860 (cond ((or (and ($listp exp
) (null (cdr arglist
)))
862 (or (null (cdr arglist
)) (null (cddr arglist
)))))
864 ((not specp
) (setq exp
(nformat exp
)))
865 (t (setq exp
(nformat-all exp
)))))
866 ((specrepp exp
) (setq exp
(specdisrep exp
))))
867 (when (and (atom exp
) (null $partswitch
))
868 (merror (intl:gettext
"~:M: argument must be a non-atomic expression; found ~:M") fn exp
))
869 (when (and inflag specp
)
870 (setq exp
(copy-tree exp
)))
872 ;; Replace all occurrences of 'rat with 'mquotient when in subst.
873 (setq exp
(let (($simp nil
)) (maxima-substitute 'mquotient
'rat exp
))))
875 start
(cond ((or (atom exp
) (eq (caar exp
) 'bigfloat
)) (go err
))
876 ((equal (setq arg
(if substflag
(meval (car arglist
)) (car arglist
)))
878 (setq arglist
(cdr arglist
))
879 (cond ((mnump substitem
)
880 (merror (intl:gettext
"~:M: argument cannot be a number; found ~M") fn substitem
))
882 (if (eq (caar exp
) 'mqapply
)
883 (prog2 (setq exp
(cadr exp
)) (go start
))
884 ;; NOT CLEAR WHAT IS INVALID HERE. OH WELL.
885 (merror (intl:gettext
"~:M: invalid operator.") fn
)))
886 (t (setq $piece
(getop (mop exp
)))
889 (setq substitem
(getopr (meval substitem
)))
890 (cond ((mnump substitem
)
891 (merror (intl:gettext
"~:M: argument cannot be a number; found ~M") fn substitem
))
892 ((not (atom substitem
))
893 (if (not (eq (caar exp
) 'mqapply
))
894 (rplaca (rplacd exp
(cons (car exp
)
897 (rplaca (cdr exp
) substitem
)
898 (return (resimplify exp
*)))
899 ((eq (caar exp
) 'mqapply
)
900 (rplacd exp
(cddr exp
))))
901 (rplaca exp
(cons substitem
902 (if (and (member 'array
(cdar exp
) :test
#'eq
)
903 (not (mopp substitem
)))
907 (rplacd exp
(cdr (box (copy-tree exp
) dispflag
)))
908 (rplaca exp
(if (eq dispflag t
)
912 (t (when arglist
(setq exp $piece
) (go a
))
914 ((not (atom arg
)) (go several
))
916 (merror (intl:gettext
"~:M: argument must be an integer; found ~M") fn arg
))
917 ((< arg
0) (go bad
)))
918 (if (eq (caar exp
) 'mqapply
) (setq exp
(cdr exp
)))
919 loop
(cond ((not (zerop arg
)) (setq arg
(1- arg
) exp
(cdr exp
))
920 (if (null exp
) (go err
))
922 ((null (setq arglist
(cdr arglist
)))
923 (return (cond (substflag (setq $piece
(resimplify (car exp
)))
924 (rplaca exp
(meval substitem
))
926 (dispflag (setq $piece
(resimplify (car exp
)))
927 (rplaca exp
(box (car exp
) dispflag
))
929 (inflag (setq $piece
(car exp
)))
930 (t (setq $piece
(simplify (car exp
))))))))
932 a
(cond ((and (not inflag
) (not specp
)) (setq exp
(nformat exp
)))
933 ((specrepp exp
) (setq exp
(specdisrep exp
))))
935 err
(cond ((eq $partswitch
'mapply
)
936 (merror (intl:gettext
"~M: invalid index ~M of list or matrix.") fn
(car arglist
)))
937 ($partswitch
(return (setq $piece
'$end
)))
938 (t (merror (intl:gettext
"~:M: fell off the end.") fn
)))
939 bad
(improper-arg-err arg fn
)
941 (if (or (not (member (caar arg
) '(mlist $allbut
) :test
#'eq
)) (cdr arglist
))
943 (setq exp1
(cons (caar exp
) (if (member 'array
(cdar exp
) :test
#'eq
) '(array))))
944 (if (eq (caar exp
) 'mqapply
)
945 (setq sevlist
(list (cadr exp
) exp1
) exp
(cddr exp
))
946 (setq sevlist
(ncons exp1
) exp
(cdr exp
)))
947 (setq arg1
(cdr arg
) prevcount
0 exp1 exp
)
949 (if (not (fixnump arg
*))
950 (merror (intl:gettext
"~:M: argument must be an integer; found ~M") fn arg
*)))
951 (when (and specp
(eq (caar arg
) 'mlist
))
952 (if substflag
(setq lastelem
(car (last arg1
))))
953 (setq arg1
(sort (copy-list arg1
) #'<)))
954 (when (eq (caar arg
) '$allbut
)
955 (setq n
(length exp
))
957 (if (or (< i
1) (> i n
))
958 (merror (intl:gettext
"~:M: index must be in range 1 to ~M, inclusive; found ~M") fn n i
)))
959 (do ((i n
(1- i
)) (arg2))
960 ((= i
0) (setq arg1 arg2
))
961 (if (not (member i arg1
:test
#'equal
)) (setq arg2
(cons i arg2
))))
962 (if substflag
(setq lastelem
(car (last arg1
)))))
963 (if (null arg1
) (if specp
(go bad
) (go end
)))
964 (if substflag
(setq lastcount lastelem
))
967 (setq count
(- (car arg1
) prevcount
) prevcount
(car arg1
))
968 (setq count
(car arg1
)))
969 (if (< count
1) (go bad
))
970 (if (and substflag
(< (car arg1
) lastelem
))
971 (setq lastcount
(1- lastcount
)))
972 count
(cond ((null exp
) (go err
))
973 ((not (= count
1)) (setq count
(1- count
) exp
(cdr exp
)) (go count
)))
974 (setq sevlist
(cons (car exp
) sevlist
))
975 (setq arg1
(cdr arg1
))
976 end
(cond ((null arg1
)
977 (setq sevlist
(nreverse sevlist
))
978 (setq $piece
(if (or inflag
(not specp
))
980 (resimplify sevlist
)))
981 (return (cond (substflag (rplaca (nthcdr (1- lastcount
) exp1
)
984 (dispflag (rplaca exp
(box (car exp
) dispflag
))
987 (substflag (if (null (cdr exp
)) (go err
))
988 (rplaca exp
(cadr exp
)) (rplacd exp
(cddr exp
)))
989 (dispflag (rplaca exp
(box (car exp
) dispflag
))
990 (setq exp
(cdr exp
)))
995 (or (and (symbolp x
) (get x
'op
)) x
))
1000 (eq (caar x
) 'mlist
)))
1002 (defmfun $cons
(x e
)
1003 (atomchk (setq e
(format1 e
)) '$cons t
)
1004 (mcons-exp-args e
(cons x
(margs e
))))
1006 (defmfun $endcons
(x e
)
1007 (atomchk (setq e
(format1 e
)) '$endcons t
)
1008 (mcons-exp-args e
(append (margs e
) (ncons x
))))
1010 (defmfun $reverse
(e)
1011 (atomchk (setq e
(format1 e
)) '$reverse nil
)
1012 (mcons-exp-args e
(reverse (margs e
))))
1014 (defmfun $append
(&rest args
)
1017 (let ((arg1 (specrepcheck (first args
))) op arrp
)
1018 (atomchk arg1
'$append nil
)
1020 arrp
(if (member 'array
(cdar arg1
) :test
#'eq
) t
))
1024 (mapcar #'(lambda (u)
1025 (atomchk (setq u
(specrepcheck u
)) '$append nil
)
1026 (unless (and (alike1 op
(mop u
))
1027 (eq arrp
(if (member 'array
(cdar u
) :test
#'eq
) t
)))
1028 (merror (intl:gettext
"append: operators of arguments must all be the same.")))
1032 (defun mcons-exp-args (e args
)
1033 (if (eq (caar e
) 'mqapply
)
1034 (list* (delsimp (car e
)) (cadr e
) args
)
1035 (cons (delsimp (car e
)) args
)))
1037 (defmfun $member
(x e
)
1038 (atomchk (setq e
($totaldisrep e
)) '$member t
)
1039 (if (memalike ($totaldisrep x
) (margs e
)) t
))
1041 (defun atomchk (e fun
2ndp
)
1042 (if (or (atom e
) (eq (caar e
) 'bigfloat
))
1043 (merror (intl:gettext
"~:M: ~Margument must be a non-atomic expression; found ~M") fun
(if 2ndp
"2nd " "") e
)))
1046 (cond (($listp e
) e
)
1047 ($inflag
(specrepcheck e
))
1051 (atomchk (setq e
(format1 e
)) '$first nil
)
1052 (if (null (cdr e
)) (merror (intl:gettext
"first: empty argument.")))
1055 ;; This macro is used to create functions second thru tenth.
1056 ;; Rather than try to modify mformat for ~:R, use the quoted symbol
1058 (macrolet ((make-nth (si i
)
1059 (let ((sim (intern (concatenate 'string
"$" (symbol-name si
)))))
1061 (atomchk (setq e
(format1 e
)) ',sim nil
)
1062 (if (< (length (margs e
)) ,i
)
1063 (merror (intl:gettext
"~:M: no such element in ~M") ',sim e
))
1071 (make-nth seventh
7)
1074 (make-nth tenth
10))
1076 (defmfun $rest
(e &optional
(n 1 n?
))
1077 (prog (m fun fun1 revp
)
1078 (when (and n?
(equal n
0))
1080 (atomchk (setq m
(format1 e
)) '$rest nil
)
1081 (cond ((and n?
(not (fixnump n
)))
1082 (merror (intl:gettext
"rest: second argument, if present, must be an integer; found ~M") n
))
1084 (setq n
(- n
) revp t
)))
1085 (if (< (length (margs m
)) n
)
1088 (merror (intl:gettext
"rest: fell off the end."))))
1090 (when (eq (car fun
) 'mqapply
)
1094 (when revp
(setq m
(reverse m
)))
1095 (setq m
(nthcdr n m
))
1096 (setq m
(cons (if (eq (car fun
) 'mlist
) fun
(delsimp fun
))
1097 (if revp
(nreverse m
) m
)))
1098 (when (eq (car fun
) 'mqapply
)
1099 (return (cons (car m
) (cons fun1
(cdr m
)))))
1103 (atomchk (setq e
(format1 e
)) '$last nil
)
1104 (when (null (cdr e
))
1105 (merror (intl:gettext
"last: empty argument.")))
1108 (defmfun $firstn
(e n
)
1109 (atomchk (setq e
(format1 e
)) '$firstn nil
)
1110 (if (or (not (fixnump n
)) (minusp n
))
1111 (merror (intl:gettext
"firstn: second argument must be a nonnegative integer; found: ~M") n
))
1112 (let ((m ($length e
)))
1115 ($rest e
(- n m
)))))
1117 (defmfun $lastn
(e n
)
1118 (atomchk (setq e
(format1 e
)) '$lastn nil
)
1119 (if (or (not (fixnump n
)) (minusp n
))
1120 (merror (intl:gettext
"lastn: second argument must be a nonnegative integer; found: ~M") n
))
1121 (let ((m ($length e
)))
1124 ($rest e
(- m n
)))))
1127 (atomchk (setq e
(format1 e
)) '$args nil
)
1128 (cons '(mlist) (margs e
)))
1130 (defmfun $delete
(x l
&optional
(n -
1 n?
))
1131 (when (and n?
(or (not (fixnump n
)) (minusp n
))) ; if n is set, it must be a nonneg fixnum
1132 (merror (intl:gettext
"delete: third argument, if present, must be a nonnegative integer; found ~M") n
))
1133 (atomchk (setq l
(specrepcheck l
)) '$delete t
)
1134 (setq x
(specrepcheck x
)
1135 l
(cons (delsimp (car l
)) (copy-list (cdr l
))))
1136 (do ((l1 (if (eq (caar l
) 'mqapply
) (cdr l
) l
)))
1137 ((or (null (cdr l1
)) (zerop n
)) l
)
1138 (if (alike1 x
(specrepcheck (cadr l1
)))
1141 (rplacd l1
(cddr l1
)))
1142 (setq l1
(cdr l1
)))))
1144 (defmfun $length
(e)
1145 (setq e
(cond (($listp e
) e
)
1146 ((or $inflag
(not ($ratp e
))) (specrepcheck e
))
1147 (t ($ratdisrep e
))))
1148 (cond ((symbolp e
) (merror (intl:gettext
"length: argument cannot be a symbol; found ~:M") e
))
1149 ((or (numberp e
) (eq (caar e
) 'bigfloat
))
1150 (if (and (not $inflag
) (mnegp e
))
1152 (merror (intl:gettext
"length: argument cannot be a number; found ~:M") e
)))
1153 ((or $inflag
(not (member (caar e
) '(mtimes mexpt
) :test
#'eq
))) (length (margs e
)))
1154 ((eq (caar e
) 'mexpt
)
1155 (if (and (alike1 (caddr e
) '((rat simp
) 1 2)) $sqrtdispflag
) 1 2))
1156 (t (length (cdr (nformat e
))))))
1159 (setq x
(specrepcheck x
)) (or (atom x
) (eq (caar x
) 'bigfloat
)))
1161 (defmfun $symbolp
(x)
1162 (setq x
(specrepcheck x
)) (symbolp x
))
1167 ((eq (caar e
) 'mrat
) ($ratnumer e
))
1168 ((eq (caar e
) 'rat
) (cadr e
))
1169 ((eq (caar (setq x
(nformat e
))) 'mquotient
) (simplify (cadr x
)))
1170 ((and (eq (caar x
) 'mminus
) (not (atom (setq x
(cadr x
))))
1171 (eq (caar x
) 'mquotient
))
1172 (simplify (list '(mtimes) -
1 (cadr x
))))
1177 ((eq (caar e
) 'mrat
) ($ratdenom e
))
1178 ((eq (caar e
) 'rat
) (caddr e
))
1179 ((or (eq (caar (setq e
(nformat e
))) 'mquotient
)
1180 (and (eq (caar e
) 'mminus
) (not (atom (setq e
(cadr e
))))
1181 (eq (caar e
) 'mquotient
)))
1182 (simplify (caddr e
)))
1185 (defmfun $entier
(e) (take '($floor
) e
))
1187 (defmfun $fix
(e) (take '($floor
) e
))
1189 ;; Evaluate THUNK ignoring any error that might occur. If the THUNK
1190 ;; returns a good number (not infinity or NaN), then return the
1191 ;; number. Otherwise, return NIL to indicate that we the computation
1192 ;; failed. This is a pretty brute-force approach.
1193 (defun try-float-computation (thunk)
1195 (errcatch (cons bindlist loclist
))
1197 (let ((result (errset (funcall thunk
))))
1198 (labels ((bad-number-p (num)
1200 (or (bad-number-p (realpart num
))
1201 (bad-number-p (imagpart num
)))
1202 (or (float-inf-p num
)
1203 (float-nan-p num
)))))
1204 (if (and result
(bad-number-p (car result
)))
1210 (let ((e1 (float e
))) (if (float-inf-p e1
) (signal 'floating-point-overflow
) e1
)))
1212 ;; Handle %i specially.
1214 ((and (symbolp e
) (mget e
'$numer
)))
1215 ((or (atom e
) (member 'array
(cdar e
) :test
#'eq
)) e
)
1216 ((eq (caar e
) 'rat
) (fpcofrat e
))
1217 ((eq (caar e
) 'bigfloat
) (fp2flo e
))
1218 ((member (caar e
) '(mexpt mncexpt
) :test
#'eq
)
1219 ;; avoid x^2 -> x^2.0, allow %e^%pi -> 23.14
1220 (let ((res (recur-apply #'$float e
)))
1223 (list (ncons (caar e
)) ($float
(cadr e
)) (caddr e
)))))
1224 ((and (eq (caar e
) '%log
)
1225 (complex-number-p (second e
) '$ratnump
))
1226 ;; Basically we try to compute float(log(x)) as directly as
1227 ;; possible, expecting Lisp to return some error if it can't.
1228 ;; Then we do a more complicated approach to compute the
1229 ;; result. However, gcl and ecl don't signal errors in these
1230 ;; cases, so we always use the complicated approach for these lisps.
1231 (let ((n (second e
)))
1233 ;; float(log(int)). First try to compute (log
1234 ;; (float n)). If that works, we're done.
1235 ;; Otherwise we need to do more.
1236 (to (or (try-float-computation #'(lambda ()
1238 (let ((m (integer-length n
)))
1239 ;; Write n as (n/2^m)*2^m where m is the number of
1240 ;; bits in n. Then log(n) = log(2^m) + log(n/2^m).
1241 ;; n/2^m is approximately 1, so converting that to a
1242 ;; float is no problem. log(2^m) = m * log(2).
1244 (log (float (/ n
(ash 1 m
)))))))))
1246 ;; float(log(n/m)) where n and m are integers. Try computing
1247 ;; it first. If it fails, compute as log(n) - log(m).
1248 (let ((try (try-float-computation #'(lambda()
1249 (log (fpcofrat n
))))))
1252 (sub ($float
`((%log
) ,(second n
)))
1253 ($float
`((%log
) ,(third n
)))))))
1254 ((complex-number-p n
'integerp
)
1255 ;; float(log(n+m*%i)).
1256 (let ((re ($realpart n
))
1258 (to (or (try-float-computation #'(lambda ()
1259 (log (complex (float re
)
1261 (let* ((size (max (integer-length re
)
1262 (integer-length im
)))
1263 (scale (ash 1 size
)))
1264 (+ (* size
(log 2e0
))
1265 (log (complex (float (/ re scale
))
1266 (float (/ im scale
))))))))))
1268 ;; log(n1/d1 + n2/d2*%i) =
1269 ;; log(s*(n+m*%i)) = log(s) + log(n+m*%i)
1271 ;; where s = lcm(d1, d2), n and m are integers
1273 (let* ((s (lcm ($denom
($realpart n
))
1274 ($denom
($imagpart n
))))
1275 (p ($expand
(mul s n
))))
1276 (add ($float
`((%log
) ,s
))
1277 ($float
`((%log
) ,p
))))))))
1278 ((and (eq (caar e
) '%erf
)
1279 (eq (second e
) '$%i
))
1280 ;; Handle like erf(%i). float(%i) (via recur-apply below)
1281 ;; just returns %i, so we never numerically evaluate it.
1282 (complexify (complex-erf (complex 0 1d0
))))
1283 ((or (eq (caar e
) '%derivative
) (eq (caar e
) '$diff
))
1284 (append (list (remove 'simp
(first e
)) ($float
(second e
))) (rest (rest e
))))
1285 (t (recur-apply #'$float e
))))
1287 (defmfun $coeff
(e x
&optional
(n 1))
1290 (coeff e
(power x n
) 1)))
1292 (defun coeff (e var pow
)
1294 (cond ((alike1 e var
) (if (equal pow
1) 1 0))
1295 ((atom e
) (if (equal pow
0) e
0))
1296 ((eq (caar e
) 'mexpt
)
1297 (cond ((alike1 (cadr e
) var
)
1298 (if (or (equal pow
0) (not (alike1 (caddr e
) pow
))) 0 1))
1301 ((or (eq (caar e
) 'mplus
) (mbagp e
))
1302 (cons (if (eq (caar e
) 'mplus
) '(mplus) (car e
))
1303 (mapcar #'(lambda (e) (coeff e var pow
)) (cdr e
))))
1304 ((eq (caar e
) 'mrat
) (ratcoeff e var pow
))
1305 ((equal pow
0) (if (coeff-contains-powers e var
) 0 e
))
1306 ((eq (caar e
) 'mtimes
)
1307 (let ((term (if (equal pow
1) var
(power var pow
))))
1308 (if (memalike term
(cdr e
)) ($delete term e
1) 0)))
1311 (defun coeff-contains-powers (e var
)
1312 (cond ((alike1 e var
) t
)
1314 ((eq (caar e
) 'mexpt
)
1315 (alike1 (cadr e
) var
))
1316 ((eq (caar e
) 'mtimes
)
1317 (member t
(mapcar #'(lambda (e)
1318 (coeff-contains-powers e var
)) (cdr e
))))
1321 (let (my-powers my-num my-flag
)
1322 (declare (special my-powers my-num my-flag
))
1324 (defmfun $hipow
(e var
)
1325 (findpowers e t var
))
1327 ;; These work best on expanded "simple" expressions.
1329 (defmfun $lopow
(e var
)
1330 (findpowers e nil var
))
1332 (defun findpowers (e hiflg var
)
1333 (let (my-powers my-num my-flag
)
1334 (declare (special my-powers my-num my-flag
))
1335 (findpowers1 e hiflg var
)
1336 (cond ((null my-powers
) (if (null my-num
) 0 my-num
))
1337 (t (when my-num
(setq my-powers
(cons my-num my-powers
)))
1338 (maximin my-powers
(if hiflg
'$max
'$min
))))))
1340 (defun findpowers1 (e hiflg var
)
1341 (cond ((alike1 e var
) (checkpow 1 hiflg
))
1343 ((eq (caar e
) 'mplus
)
1344 (cond ((not (freel (cdr e
) var
))
1345 (do ((e (cdr e
) (cdr e
))) ((null e
))
1346 (setq my-flag nil
) (findpowers1 (car e
) hiflg var
)
1347 (if (null my-flag
) (checkpow 0 hiflg
))))))
1348 ((and (eq (caar e
) 'mexpt
) (alike1 (cadr e
) var
)) (checkpow (caddr e
) hiflg
))
1349 ((specrepp e
) (findpowers1 (specdisrep e
) hiflg var
))
1350 (t (mapc #'(lambda (x) (findpowers1 x hiflg var
)) (cdr e
)))))
1352 (defun checkpow (pow hiflg
)
1354 (cond ((not (numberp pow
)) (setq my-powers
(cons pow my-powers
)))
1355 ((null my-num
) (setq my-num pow
))
1356 (hiflg (if (> pow my-num
) (setq my-num pow
)))
1357 ((< pow my-num
) (setq my-num pow
)))))
1359 ;; push(x,l): The argument l must be a mapatom that is bound to a list. Prepend x to
1360 ;; the value of l and return this list. The arguments to push are evaluated from left to right.
1362 ;; We assume that if each member of a list is simplified, the list members are still simplified
1363 ;; after pushing onto (or popping) the list. Just before returning, the code calls simplifya
1364 ;; with second argument true. Without this call to simplifya, the general simplifier would simplify
1365 ;; every list member after returning. (Barton Willis, author of $push and $pop)
1368 (let* ((o (car (pop z
)))
1369 (x (if z
(pop z
) (wna-err o
)))
1370 (l (if z
(pop z
) (wna-err o
)))
1376 (setq lo
(simplifya (cons (list (meval (mop l
)) 'array
) (mevalargs (margs l
))) t
))
1377 (setq ll
(let ((noevalargs t
)) (meval lo
))))
1380 (setq ll
(meval l
))))
1381 (if (and ($mapatom lo
) ($listp ll
))
1383 (setq ll
(cons x
(cdr ll
)))
1384 (mset lo
(simplifya (cons '(mlist) ll
) t
)))
1385 (merror "Second argument to push must be a mapatom that is bound to a list"))))
1387 ;; pop(l): The argument l must be a mapatom that is bound to a list. Remove the
1388 ;; first member of the value of l and return it.
1390 ;; We assume that if each member of a list is simplified, the list members are still simplified
1391 ;; after popping (or pushing onto) the list. Just before returning, the code calls simplifya
1392 ;; with second argument true. Without this call to simplifya, the general simplifier would simplify
1393 ;; every list member after returning.
1396 (let* ((o (car (pop z
)))
1397 (l (if z
(pop z
) (wna-err o
)))
1401 (setq lo
(simplifya (cons (list (meval (mop l
)) 'array
) (mevalargs (margs l
))) t
))
1402 (setq ll
(let ((noevalargs t
)) (meval lo
))))
1405 (setq ll
(meval lo
))))
1407 (cond ((and ($mapatom lo
) ($listp ll
))
1408 (if ($emptyp ll
) (merror "Pop called on an empty list")
1412 (mset lo
(simplifya (cons '(mlist) ll
) t
)))))
1413 (t (merror "Argument to pop must be a mapatom that is bound to a list")))))
