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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 (declare-top (special *rulesw
*
16 *prods
* *negprods
* *sums
* *negsums
*
20 ;; Switches dealing with matrices and non-commutative multiplication.
22 (defmvar $doscmxplus nil
23 "Causes SCALAR + MATRIX to return a matrix answer. This switch
24 is not subsumed under DOALLMXOPS.")
27 "Causes SCALAR^MATRIX([1,2],[3,4]) to return
28 MATRIX([SCALAR,SCALAR^2],[SCALAR^3,SCALAR^4]). In general, this
29 transformation affects exponentiations where the *print-base* is a
30 scalar and the power is a matrix or list.")
32 (defmvar $domxplus nil
)
34 (defmvar $domxtimes nil
)
38 ;; Lisp level variables
39 (defmvar derivsimp t
"Hack in `simpderiv' for `rwg'")
41 (defmvar $grindswitch nil
)
44 (defmvar $limitdomain
'$complex
)
45 (defmvar $lognegint nil
47 :properties
((evflag t
)))
48 (defmvar $distribute_over t
) ; If T, functions are distributed over bags.
52 (defvar %e-val
(mget '$%e
'$numer
))
53 (defvar %pi-val
(mget '$%pi
'$numer
))
55 (defvar expandflag nil
)
58 (defprop mnctimes t associative
)
59 (defprop lambda t lisp-no-simp
)
61 ;; Local functions should not be simplified. Various lisps
62 ;; use various names for the list structure defining these:
67 (dolist (y (list x z
))
70 (setf (get (car y
) 'lisp-no-simp
) t
))))))
72 (dolist (x '(mplus mtimes mnctimes mexpt mncexpt %sum
))
73 (setf (get x
'msimpind
) (cons x
'(simp))))
75 ;; operators properties
77 (mapc #'(lambda (x) (setf (get (first x
) 'operators
) (second x
)))
78 '((mplus simplus
) (mtimes simptimes
) (mncexpt simpncexpt
)
79 (mminus simpmin
) (%gamma simpgamma
) (mfactorial simpfact
)
80 (mnctimes simpnct
) (mquotient simpquot
) (mexpt simpexpt
)
82 (%derivative simpderiv
)
84 (%integrate simpinteg
) (%limit simp-limit
)
85 (bigfloat simpbigfloat
) (lambda simplambda
) (mdefine simpmdef
)
86 (mqapply simpmqapply
) (%gamma simpgamma
)
87 ($beta simpbeta
) (%sum simpsum
) (%binomial simpbinocoef
)
88 (%plog simpplog
) (%product simpprod
) (%genfact simpgfact
)
89 ($atan2 simpatan2
) ($matrix simpmatrix
) (%matrix simpmatrix
)
90 ($bern simpbern
) ($euler simpeuler
)))
92 (defprop $li lisimp specsimp
)
93 (defprop $psi psisimp specsimp
)
95 (defprop $equal t binary
)
96 (defprop $notequal t binary
)
99 "Returns true if X is a bigfloat"
102 (eq (caar x
) 'bigfloat
)))
105 "Returns non-NIL if X is an integer, float, or bfloat that is equal
107 (or (and (integerp x
) (= 0 x
))
108 (and (floatp x
) (= 0.0 x
))
109 (and ($bfloatp x
) (= 0 (second x
)))))
112 "Returns non-NIL if X is an integer, float, or bfloat that is equal
114 (or (and (integerp x
) (= 1 x
))
115 (and (floatp x
) (= 1.0 x
))
116 (and ($bfloatp x
) (zerop1 (sub x
1)))))
119 "Returns non-NIL if X is a Lisp number or if it is a Maxima rational
120 form or a bigfloat form"
122 (and (not (atom x
)) (not (atom (car x
)))
123 (member (caar x
) '(rat bigfloat
)))))
125 ;; Does X or a subexpression match PREDICATE?
127 ;; If X is a tree then we recurse depth-first down its arguments. The specrep
128 ;; check is because rat forms are built rather differently from normal Maxima
129 ;; expressions so we need to unpack them for the recursion to work properly.
130 (defun subexpression-matches-p (predicate x
)
131 (or (funcall predicate x
)
134 (subexpression-matches-p predicate
(specdisrep x
))
135 (some (lambda (arg) (subexpression-matches-p predicate arg
))
138 ;; Is there a bfloat anywhere in X?
139 (defun some-bfloatp (x) (subexpression-matches-p '$bfloatp x
))
141 ;; Is there a float anywhere in X?
142 (defun some-floatp (x) (subexpression-matches-p 'floatp x
))
144 ;; EVEN works for any arbitrary lisp object since it does an integer
145 ;; check first. In other cases, you may want the Lisp EVENP function
146 ;; which only works for integers.
148 (defun even (a) (and (integerp a
) (not (oddp a
))))
150 ;; Predicates to determine if X satisfies some condition.
153 "Determines if X is a Maxima rational form: ((rat ...) a b)"
154 (and (not (atom x
)) (eq (caar x
) 'rat
)))
157 "Determines if X is a Maxima sum form: ((mplus ...) ...)"
158 (and (not (atom x
)) (eq (caar x
) 'mplus
)))
161 "Determines if X is a Maxima product form: ((mtimes ...) ...)"
162 (and (not (atom x
)) (eq (caar x
) 'mtimes
)))
165 "Determines if X is a Maxima exponential form: ((mexpt ...) ...)"
166 (and (not (atom x
)) (eq (caar x
) 'mexpt
)))
168 (defun mnctimesp (x) (and (not (atom x
)) (eq (caar x
) 'mnctimes
)))
170 (defun mncexptp (x) (and (not (atom x
)) (eq (caar x
) 'mncexpt
)))
173 "Determines if X is a Maxima log form: ((%log ...) ...)"
174 (and (not (atom x
)) (eq (caar x
) '%log
)))
177 "Determines if X is a Maxima negative form: ((mminus ...) ...)
179 This generally only happens on input forms like a - b:
180 ((mplus) $a ((mminus) $b)).
181 After simplification a - b becomes
182 ((mplus) $a ((mtimes) -1 $b))"
184 (and (not (atom x
)) (eq (caar x
) 'mminus
)))
187 "Determines if X is negative if X is a Lisp number or a Maxima rat
188 form or bigfloat form"
189 (cond ((realp x
) (minusp x
))
190 ((or (ratnump x
) ($bfloatp x
)) (minusp (cadr x
)))))
192 (defun mqapplyp (e) (and (not (atom e
)) (eq (caar e
) 'mqapply
)))
194 (defun ratdisrep (e) (simplifya ($ratdisrep e
) nil
))
196 (defun sratsimp (e) (simplifya ($ratsimp e
) nil
))
198 (defun simpcheck (e flag
)
199 (cond ((specrepp e
) (specdisrep e
))
201 (t (let (($%enumer $numer
))
202 ;; Switch $%enumer on, when $numer is TRUE to allow
203 ;; simplification of $%e to its numerical value.
204 (simplifya e nil
)))))
206 (defun mratcheck (e) (if ($ratp e
) (ratdisrep e
) e
))
208 (defmfun $numberp
(e)
209 "Returns true if E is a Maxima rational, a float, or a bigfloat number"
210 (or ($ratnump e
) ($floatnump e
) ($bfloatp e
)))
212 (defmfun $integerp
(x)
215 (not (member 'trunc
(car x
)))
217 (equal (cddr x
) 1))))
219 ;; The call to $INTEGERP in the following two functions checks for a CRE
220 ;; rational number with an integral numerator and a unity denominator.
223 (cond ((integerp x
) (oddp x
))
224 (($integerp x
) (oddp (cadr x
)))))
227 (cond ((integerp x
) (evenp x
))
228 (($integerp x
) (not (oddp (cadr x
))))))
230 (defmfun $floatnump
(x)
232 (and ($ratp x
) (floatp (cadr x
)) (onep1 (cddr x
)))))
234 (defmfun $ratnump
(x)
238 (not (member 'trunc
(car x
)))
240 (integerp (cddr x
)))))
245 (eq (caar x
) 'mrat
)))
247 (defmfun $taylorp
(x)
250 (member 'trunc
(cdar x
)) t
))
252 (defun specrepcheck (e) (if (specrepp e
) (specdisrep e
) e
))
254 ;; Note that the following two functions are carefully coupled.
258 (member (caar e
) '(mrat mpois
))))
260 (defun specdisrep (e)
261 (cond ((eq (caar e
) 'mrat
) (ratdisrep e
))
264 (defmfun $polysign
(x)
265 (setq x
(cadr (ratf x
)))
266 (cond ((equal x
0) 0) ((pminusp x
) -
1) (t 1)))
268 ;; These check for the correct number of operands within Macsyma expressions,
269 ;; not arguments in a procedure call as the name may imply.
271 (defun arg-count-check (required-arg-count expr
)
272 (unless (= required-arg-count
(length (rest expr
)))
273 (wna-err expr required-arg-count
)))
275 (defun oneargcheck (expr)
276 (arg-count-check 1 expr
))
278 (defun twoargcheck (expr)
279 (arg-count-check 2 expr
))
281 ;; WNA-ERR: Wrong Number of Arguments error
283 ;; If REQUIRED-ARG-COUNT is non-NIL, then we check that EXPR has the
284 ;; correct number of arguments. A informative error message is shown
285 ;; if the number of arguments is not given.
287 ;; Otherwise, EXPR must be a symbol and a generic message is printed.
288 ;; (This is for backward compatibility for existing uses of WNA-ERR.)
289 (defun wna-err (exprs &optional required-arg-count
)
290 (if required-arg-count
291 (let ((op (caar exprs
))
292 (actual-count (length (rest exprs
))))
293 (merror (intl:gettext
"~M: expected exactly ~M arguments but got ~M: ~M")
294 op required-arg-count actual-count
(list* '(mlist) (rest exprs
))))
295 (merror (intl:gettext
"~:@M: wrong number of arguments.")
298 (defun improper-arg-err (exp fn
)
299 (merror (intl:gettext
"~:M: improper argument: ~M") fn exp
))
301 (defun subargcheck (form subsharp argsharp fun
)
302 (if (or (not (= (length (subfunsubs form
)) subsharp
))
303 (not (= (length (subfunargs form
)) argsharp
)))
304 (merror (intl:gettext
"~:@M: wrong number of arguments or subscripts.") fun
)))
306 ;; Constructor and extractor primitives for subscripted functions, e.g.
307 ;; F[1,2](X,Y). SUBL is (1 2) and ARGL is (X Y).
309 ;; These will be flushed when NOPERS is finished. They will be macros in
310 ;; NOPERS instead of functions, so we have to be careful that they aren't
311 ;; mapped or applied anyplace. What we really want is open-codable routines.
313 (defun subfunmakes (fun subl argl
)
314 `((mqapply simp
) ((,fun simp array
) .
,subl
) .
,argl
))
316 (defun subfunmake (fun subl argl
)
317 `((mqapply) ((,fun simp array
) .
,subl
) .
,argl
))
319 (defun subfunname (exp) (caaadr exp
))
321 (defun subfunsubs (exp) (cdadr exp
))
323 (defun subfunargs (exp) (cddr exp
))
325 (defmfun $numfactor
(x)
326 (setq x
(specrepcheck x
))
329 ((not (eq (caar x
) 'mtimes
)) 1)
330 ((mnump (cadr x
)) (cadr x
))
333 (defun scalar-or-constant-p (x flag
)
334 (if flag
(not ($nonscalarp x
)) ($scalarp x
)))
336 (defmfun $constantp
(x)
337 (cond ((atom x
) (or ($numberp x
) (kindp x
'$constant
)))
338 ((member (caar x
) '(rat bigfloat
)) t
)
339 ((specrepp x
) ($constantp
(specdisrep x
)))
340 ((or (mopp (caar x
)) (kindp (caar x
) '$constant
))
341 (do ((x (cdr x
) (cdr x
))) ((null x
) t
)
342 (if (not ($constantp
(car x
))) (return nil
))))))
345 (cond ((symbolp x
) (kindp x
'$constant
))
347 (and (kindp (caar x
) '$constant
)
348 (do ((x (cdr x
) (cdr x
))) ((null x
) t
)
349 (if (not ($constantp
(car x
))) (return nil
)))))))
351 (defun maxima-constantp (x)
353 (and (symbolp x
) (kindp x
'$constant
))))
355 (defun consttermp (x) (and ($constantp x
) (not ($nonscalarp x
))))
357 (defmfun $scalarp
(x) (or (consttermp x
) (eq (scalarclass x
) '$scalar
)))
359 (defmfun $nonscalarp
(x) (eq (scalarclass x
) '$nonscalar
))
361 (defun scalarclass (exp) ; Returns $SCALAR, $NONSCALAR, or NIL (unknown).
363 ;; Maxima numbers are scalar.
366 (cond ((or (mget exp
'$nonscalar
)
367 (and (not (mget exp
'$scalar
))
368 ;; Arrays are nonscalar, but not if declared scalar.
370 ($member exp $arrays
))))
372 ((or (mget exp
'$scalar
)
373 ;; Include constant atoms which are not declared nonscalar.
376 ((member 'array
(car exp
))
377 (cond ((mget (caar exp
) '$scalar
) '$scalar
)
378 ((mget (caar exp
) '$nonscalar
) '$nonscalar
)
380 ((specrepp exp
) (scalarclass (specdisrep exp
)))
381 ;; If the function is declared scalar or nonscalar, then return. If it
382 ;; isn't explicitly declared, then try to be intelligent by looking at
383 ;; the arguments to the function.
384 ((scalarclass (caar exp
)))
385 ;; <number> + <scalar> is SCALARP because that seems to be useful.
386 ;; This should probably only be true if <number> is a member of the
387 ;; field of scalars. <number> * <scalar> is SCALARP since
388 ;; <scalar> + <scalar> is SCALARP. Also, this has to be done to make
389 ;; <scalar> - <scalar> SCALARP.
390 ((member (caar exp
) '(mplus mtimes
))
391 (do ((l (cdr exp
) (cdr l
))) ((null l
) '$scalar
)
392 (if (not (consttermp (car l
)))
393 (return (scalarclass-list l
)))))
394 ((and (eq (caar exp
) 'mqapply
) (scalarclass (cadr exp
))))
395 ((mxorlistp exp
) '$nonscalar
)
396 ;; If we can't find out anything about the operator, then look at the
397 ;; arguments to the operator. I think NIL should be returned at this
400 (do ((exp (cdr exp
) (cdr exp
)) (l '(1)))
401 ((null exp
) (scalarclass-list l
))
402 (if (not (consttermp (car exp
)))
403 (setq l
(cons (car exp
) l
)))))))
405 ;; Could also do <scalar> +|-|*|/ |^ <declared constant>, but this is not
406 ;; always correct and could screw somebody.
408 ;; SCALARCLASS-LIST takes a list of expressions as its argument. If their
409 ;; scalarclasses all agree, then that scalarclass is returned.
411 (defun scalarclass-list (llist)
412 (cond ((null llist
) nil
)
413 ((null (cdr llist
)) (scalarclass (car llist
)))
414 (t (let ((sc-car (scalarclass (car llist
)))
415 (sc-cdr (scalarclass-list (cdr llist
))))
416 (cond ((or (eq sc-car
'$nonscalar
)
417 (eq sc-cdr
'$nonscalar
))
419 ((and (eq sc-car
'$scalar
) (eq sc-cdr
'$scalar
))
424 (member (caar x
) '(mequal mlist $matrix
))))
426 (defun mequalp (x) (and (not (atom x
)) (eq (caar x
) 'mequal
)))
430 (member (caar x
) '(mlist $matrix
))))
432 (defun mxorlistp1 (x)
434 (or (eq (caar x
) '$matrix
)
435 (and (eq (caar x
) 'mlist
) $listarith
))))
437 (defun constfun (ign)
438 (declare (ignore ign
)) ; Arg ignored. Function used for mapping down lists.
441 (defun constmx (*const
* x
)
442 (simplifya (fmapl1 'constfun x
) t
))
444 ;;; ISINOP returns the complete subexpression with the operator OP, when the
445 ;;; operator OP is found in EXPR.
447 (defun isinop (expr op
) ; OP is assumed to be an atom
448 (cond ((atom expr
) nil
)
449 ((and (eq (caar expr
) op
)
450 (not (member 'array
(cdar expr
))))
453 (do ((expr (cdr expr
) (cdr expr
))
456 (when (setq res
(isinop (car expr
) op
))
459 (defun free (exp var
)
460 (cond ((alike1 exp var
) nil
)
463 (and (listp (car exp
))
464 (free (caar exp
) var
)
465 (freel (cdr exp
) var
)))))
468 (do ((l l
(cdr l
))) ((null l
) t
)
470 ((atom l
) (return (free l var
))) ;; second element of a pair
471 ((not (free (car l
) var
)) (return nil
)))))
474 (defun freeargs (exp var
)
475 (cond ((alike1 exp var
) nil
)
477 (t (do ((l (margs exp
) (cdr l
))) ((null l
) t
)
478 (cond ((not (freeargs (car l
) var
)) (return nil
)))))))
480 (defun simplifya (x y
)
481 (cond ((not $simp
) x
)
483 (cond ((and $%enumer $numer
(eq x
'$%e
))
484 ;; Replace $%e with its numerical value,
485 ;; when %enumer and $numer TRUE
489 (cond ((and (cdr x
) (atom (cdr x
)))
490 (merror (intl:gettext
"simplifya: malformed expression (atomic cdr).")))
491 ((get (car x
) 'lisp-no-simp
)
492 ;; this feature is to be used with care. it is meant to be
493 ;; used to implement data objects with minimum of consing.
494 ;; forms must not bash the DISPLA package. Only new forms
495 ;; with carefully chosen names should use this feature.
498 (mapcar #'(lambda (x) (simplifya x y
)) (cdr x
))))))
499 ((eq (caar x
) 'rat
) (*red1 x
))
500 ;; Enforced resimplification: Reset dosimp and strip 'simp tags from x.
501 (dosimp (let ((dosimp nil
)) (simplifya (unsimplify x
) y
)))
502 ((member 'simp
(cdar x
)) x
)
503 ((eq (caar x
) 'mrat
) x
)
504 ((not (atom (caar x
)))
505 (cond ((or (eq (caaar x
) 'lambda
)
506 (and (not (atom (caaar x
))) (eq (caaaar x
) 'lambda
)))
507 (mapply1 (caar x
) (cdr x
) (caar x
) x
))
508 (t (merror (intl:gettext
"simplifya: operator is neither an atom nor a lambda expression: ~S") x
))))
509 ((and $distribute_over
510 (get (caar x
) 'distribute_over
)
511 ;; A function with the property 'distribute_over.
512 ;; Look, if we have a bag as argument to the function.
513 (distribute-over x
)))
514 ((get (caar x
) 'opers
)
515 (let ((opers-list *opers-list
)) (oper-apply x y
)))
516 ((and (eq (caar x
) 'mqapply
)
518 (and (eq substp
'mqapply
)
519 (or (eq (car (cadr x
)) 'lambda
)
520 (eq (caar (cadr x
)) 'lambda
)))))
521 (cond ((or (symbolp (cadr x
)) (not (atom (cadr x
))))
522 (simplifya (cons (cons (cadr x
) (cdar x
)) (cddr x
)) y
))
523 ((or (not (member 'array
(cdar x
))) (not $subnumsimp
))
524 (merror (intl:gettext
"simplifya: I don't know how to simplify this operator: ~M") x
))
526 (t (let ((w (get (caar x
) 'operators
)))
528 (or (not (member 'array
(cdar x
)))
531 (t (simpargs x y
)))))))
533 ;; EQTEST returns an expression which is the same as X
534 ;; except that it is marked with SIMP and maybe other flags from CHECK.
536 ;; Following description is inferred from the code. Dunno why the function is named "EQTEST".
538 ;; (1) if X is already marked with SIMP flag or doesn't need it: return X.
539 ;; (2) if X is pretty much the same as CHECK (same operator and same arguments),
540 ;; then return CHECK after marking it with SIMP flag.
541 ;; (3) if operator of X has the MSIMPIND property, replace it
542 ;; with value of MSIMPIND (something like '(MPLUS SIMP)) and return X.
543 ;; (4) if X or CHECK is an array expression, return X after marking it with SIMP and ARRAY flags.
544 ;; (5) otherwise, return X after marking it with SIMP flag.
546 (defun eqtest (x check
)
551 (member 'simp
(cdar x
)))
553 ((and (eq (caar x
) (caar check
))
554 (equal (cdr x
) (cdr check
)))
555 (cond ((and (null (cdar check
))
556 (setq y
(get (caar check
) 'msimpind
)))
557 (cons y
(cdr check
)))
558 ((member 'simp
(cdar check
))
561 (cons (cons (caar check
)
563 (cons 'simp
(cdar check
))
566 ((setq y
(get (caar x
) 'msimpind
))
568 ((or (member 'array
(cdar x
))
569 (and (eq (caar x
) (caar check
))
570 (member 'array
(cdar check
))))
571 (rplaca x
(cons (caar x
) '(simp array
))))
573 (rplaca x
(cons (caar x
) '(simp)))))))
575 ;; A function, which distributes of bags like a list, matrix, or equation.
576 ;; Check, if we have to distribute of one of the bags or any other operator.
577 (defun distribute-over (expr)
578 (cond ((= 1 (length (cdr expr
)))
579 ;; Distribute over for a function with one argument.
580 (cond ((and (not (atom (cadr expr
)))
581 (member (caaadr expr
) (get (caar expr
) 'distribute_over
)))
584 (mapcar #'(lambda (u) (simplify (list (car expr
) u
)))
588 ;; A function with more than one argument.
589 (do ((args (cdr expr
) (cdr args
))
592 (when (and (not (atom (car args
)))
593 (member (caar (car args
))
594 (get (caar expr
) 'distribute_over
)))
595 ;; Distribute the function over the arguments and simplify again.
598 (cons (ncons (caar (car args
)))
599 (mapcar #'(lambda (u)
603 (cons (ncons (caar expr
))
604 (reverse first-args
))
607 (cdr (car args
)))))))
608 (setq first-args
(cons (car args
) first-args
))))))
610 (defun rulechk (x) (or (mget x
'oldrules
) (get x
'rules
)))
612 (defun resimplify (x) (let ((dosimp t
)) (simplifya x nil
)))
614 (defun unsimplify (x)
615 (if (or (atom x
) (specrepp x
))
617 (cons (remove 'simp
(car x
) :count
1) (mapcar #'unsimplify
(cdr x
)))))
619 (defun simpargs (x y
)
620 (if (or (and (eq (get (caar x
) 'dimension
) 'dimension-infix
)
621 (not (getl (caar x
) '($lassociative $rassociative
))))
622 (get (caar x
) 'binary
))
624 (if (and (member 'array
(cdar x
)) (null (margs x
)))
625 (merror (intl:gettext
"SIMPARGS: subscripted variable found with no subscripts.")))
626 (eqtest (if y x
(let ((flag (member (caar x
) '(mlist mequal
))))
627 (cons (ncons (caar x
))
628 (mapcar #'(lambda (u)
629 (if flag
(simplifya u nil
)
634 ;;;-----------------------------------------------------------------------------
635 ;;; ADDK (X Y) 27.09.2010/DK
637 ;;; Arguments and values:
638 ;;; X - a Maxima number
639 ;;; Y - a Maxima number
640 ;;; result - a simplified Maxima number
643 ;;; ADDK adds two Maxima numbers and returns a simplified Maxima number.
644 ;;; ADDK can be called in Lisp code, whenever the arguments are valid
645 ;;; Maxima numbers, these are integer, float, Maxima rational, or
646 ;;; Maxima bigfloat numbers. The arguments must not be simplified. The
647 ;;; precision of a bigfloat result depends on the setting of the
648 ;;; global variable $FPPREC. If the option variable $FLOAT is T, a
649 ;;; Maxima rational number as a result is converted to a float number.
653 ;;; (addk 2.0 3) -> 5.0
654 ;;; (addk ($bfloat 2) 3)-> ((BIGFLOAT SIMP 56) 45035996273704960 3)
655 ;;; (addk 2 '((rat) 1 2)) -> ((RAT SIMP) 5 2)
656 ;;; (let (($float t)) (addk 2 '((rat) 1 2))) -> 2.5
659 ;;; The option variables $FLOAT and $FPPREC.
662 ;;; TIMESK to multiply and EXPTRL to exponentiate two Maxima numbers.
665 ;;; The routine works for Lisp rational and Lisp complex numbers too.
666 ;;; This feature is not used in Maxima code. If Lisp complex and
667 ;;; rational numbers are mixed with Maxima rational or bigfloat
668 ;;; numbers the result is wrong or a Lisp error is generated.
669 ;;;-----------------------------------------------------------------------------
674 ((and (numberp x
) (numberp y
)) (+ x y
))
675 ((or ($bfloatp x
) ($bfloatp y
)) ($bfloat
(list '(mplus) x y
)))
678 (cond ((floatp x
) (return (+ x
(fpcofrat y
))))
679 (t (setq x
(list '(rat) x
1)))))
681 (cond ((floatp y
) (return (+ y
(fpcofrat x
))))
682 (t (setq y
(list '(rat) y
1))))))
683 (setq g
(gcd (caddr x
) (caddr y
)))
684 (setq a
(truncate (caddr x
) g
)
685 b
(truncate (caddr y
) g
))
686 (return (timeskl (list '(rat) 1 g
)
692 ;;;-----------------------------------------------------------------------------
693 ;;; *RED1 (X) 27.09.2010/DK
696 ;;; Arguments and values:
697 ;;; X - a Maxima rational number (for *RED1)
698 ;;; N - an integer number representing the numerator of a rational
699 ;;; D - an integer number representing the denominator of a rational
700 ;;; result - a simplified Maxima rational number
703 ;;; *RED1 is called from SIMPLIFYA to reduce and simplify a Maxima rational
704 ;;; number. *RED1 checks if the rational number is already simplified. If
705 ;;; the option variable $FLOAT is T, the rational number is converted to a
706 ;;; float number. If the number is not simplified, *RED is called.
708 ;;; *RED reduces the numerator N and the demoniator D and returns a
709 ;;; simplified Maxima rational number. The result is converted to a float
710 ;;; number, if the option variable $FLOAT is T.
713 ;;; The option variable $FLOAT.
714 ;;;-----------------------------------------------------------------------------
717 (cond ((member 'simp
(cdar x
))
718 (if $float
(fpcofrat x
) x
))
719 (t (*red
(cadr x
) (caddr x
)))))
724 (t (let ((u (gcd n d
)))
725 (setq n
(truncate n u
)
727 (if (minusp d
) (setq n
(- n
) d
(- d
)))
728 (cond ((equal d
1) n
)
729 ($float
(fpcofrat1 n d
))
730 (t (list '(rat simp
) n d
)))))))
732 ;;;-----------------------------------------------------------------------------
733 ;;; TIMESK (X Y) 27.09.2010/DK
735 ;;; Arguments and values:
736 ;;; X - a Maxima number
737 ;;; Y - a Maxima number
738 ;;; result - a simplified Maxima number
741 ;;; TIMESK Multiplies two Maxima numbers and returns a simplified Maxima
742 ;;; number. TIMESK can be called in Lisp code, whenever the arguments are
743 ;;; valid Maxima numbers, these are integer, float, Maxima rational, or
744 ;;; Maxima bigfloat numbers. The arguments must not be simplified. The
745 ;;; precision of a bigfloat result depends on the setting of the
746 ;;; global variable $FPPREC. If the option variable $FLOAT is T, a
747 ;;; Maxima rational number as a result is converted to a float number.
749 ;;; TIMESKL is called from TIMESK to multiply two Maxima rational numbers or
750 ;;; a rational number with an integer number.
753 ;;; (timesk 2 3) -> 6
754 ;;; (timesk 2.0 3) -> 6.0
755 ;;; (timesk ($bfloat 2) 3)-> ((BIGFLOAT SIMP 56) 54043195528445952 3)
756 ;;; (timesk 3 '((rat) 1 2)) -> ((RAT SIMP) 3 2)
757 ;;; (let (($float t)) (timesk 3 '((rat) 1 2))) -> 1.5
760 ;;; The option variables $FLOAT and $FPPREC.
763 ;;; ADDK to add and EXPTRL to exponentiate two Maxima numbers.
766 ;;; The routine works for Lisp rational and Lisp complex numbers too.
767 ;;; This feature is not used in Maxima code. If Lisp complex and
768 ;;; rational numbers are mixed with Maxima rational or bigfloat
769 ;;; numbers the result is wrong or a Lisp error is generated.
770 ;;;-----------------------------------------------------------------------------
772 ;; NUM1 and DENOM1 are helper functions for TIMESKL to get the numerator and the
773 ;; denominator of an integer or Maxima rational number. For an integer the
774 ;; denominator is 1. Both functions are used at other places in Maxima code too.
777 (if (numberp a
) a
(cadr a
)))
780 (if (numberp a
) 1 (caddr a
)))
782 (defun timesk (x y
) ; X and Y are assumed to be already reduced
783 (cond ((equal x
1) y
)
785 ((and (numberp x
) (numberp y
)) (* x y
))
786 ((or ($bfloatp x
) ($bfloatp y
)) ($bfloat
(list '(mtimes) x y
)))
787 ((floatp x
) (* x
(fpcofrat y
)))
788 ((floatp y
) (* y
(fpcofrat x
)))
791 ;; TIMESKL takes one or two Maxima rational numbers, one argument can be an
792 ;; integer number. The result is a Maxima rational or an integer number.
793 ;; If the option variable $FLOAT is T, a Maxima rational number in converted
798 (setq u
(*red
(num1 x
) (denom1 y
)))
799 (setq v
(*red
(num1 y
) (denom1 x
)))
800 (setq g
(cond ((or (equal u
0) (equal v
0)) 0)
802 ((and (numberp u
) (numberp v
)) (* u v
))
804 (* (num1 u
) (num1 v
))
805 (* (denom1 u
) (denom1 v
))))))
806 (return (cond ((numberp g
) g
)
807 ((equal (caddr g
) 1) (cadr g
))
808 ($float
(fpcofrat g
))
811 ;;;-----------------------------------------------------------------------------
812 ;;; FPCOFRAT (RATNO) 27.09.2010/DK
815 ;;; Arguments and values:
816 ;;; RATNO - a Maxima rational number (for FPCOFRAT)
817 ;;; NU - an integer number which represents the numerator of a rational
818 ;;; D - an integer number which represents the denominator of a rational
819 ;;; result - floating point approximation of a rational number
822 ;;; Floating Point Conversion OF RATional number routine.
823 ;;; Finds floating point approximation to rational number.
825 ;;; FPCOFRAT1 computes the quotient of NU/D.
827 ;;; Exceptional situations:
828 ;;; A Lisp error is generated, if the rational number does not fit into a
830 ;;;-----------------------------------------------------------------------------
832 ;; This constant is only needed in the file float.lisp.
834 (:compile-toplevel
:load-toplevel
:execute
)
835 (defconstant machine-mantissa-precision
(float-digits 1.0)))
837 (defun fpcofrat (ratno)
838 (fpcofrat1 (cadr ratno
) (caddr ratno
)))
840 (defun fpcofrat1 (nu d
)
843 ;;;-----------------------------------------------------------------------------
844 ;;; EXPTA (X Y) 27.09.2010/DK
846 ;;; Arguments and values:
847 ;;; X - a Maxima number
848 ;;; Y - an integer number
849 ;;; result - a simplified Maxima number
852 ;;; Computes X^Y, where X is Maxima number and Y an integer. The result is
853 ;;; a simplified Maxima number. Y can be a rational Maxima number. For this
854 ;;; case the numerator is taken as the power.
857 ;;; The option variables $FLOAT and $FPPREC.
860 ;;; This routine is not used within the simplifier. There is only one
861 ;;; call from the file hayat.lisp. This call can be replaced with a
862 ;;; call of the function power.
863 ;;;-----------------------------------------------------------------------------
871 ($bfloat
(list '(mexpt) x y
)))
873 (*red
(exptb (caddr x
) (- (num1 y
)))
874 (exptb (cadr x
) (- (num1 y
)))))
876 (*red
(exptb (cadr x
) (num1 y
))
877 (exptb (caddr x
) (num1 y
))))))
879 ;;;-----------------------------------------------------------------------------
880 ;;; EXPTB (A B) 27.09.2010/DK
882 ;;; Arguments and values:
883 ;;; A - a float or integer number
884 ;;; B - an integer number
885 ;;; result - a simplified Maxima number
888 ;;; Computes A^B, where A is a float or an integer number and B is an
889 ;;; integer number. The result is an integer, float, or Maxima
894 ;;; (exptb 3.0 2) -> 9.0
895 ;;; (exptb 3 -2) -> ((RAT SiMP) 1 9)
896 ;;; (let (($float t)) (exptb 3 -2)) -> 0.1111111111111111
899 ;;; The option variable $FLOAT.
902 ;;; EXPTB calls the Lisp functions EXP or EXPT to compute the result.
903 ;;;-----------------------------------------------------------------------------
907 (cond ((equal a %e-val
)
908 ;; Make B a float so we'll get double-precision result.
910 ((or (floatp a
) (not (minusp b
)))
913 (setq b
(expt a
(- b
)))
915 (when (float-inf-p result
) ;; needed for gcl and sbcl - (sometimes) no trap of overflow
916 (error 'floating-point-overflow
))
920 ;;;-----------------------------------------------------------------------------
921 ;;; SIMPLUS (X W Z) 27.09.2010/DK
923 ;;; Arguments and values:
924 ;;; X - a Maxima expression of the form ((mplus) term1 term2 ...)
925 ;;; W - an arbitrary value, the value is ignored
926 ;;; Z - T or NIL, if T the arguments are assumed to be simplified
927 ;;; result - a simplified mplus-expression or an atom
930 ;;; Implementation of the simplifier for the "+" operator.
931 ;;; A general description of SIMPLUS can be found in the paper:
932 ;;; http://www.cs.berkeley.edu/~fateman/papers/simplifier.txt
935 ;;; The addition of matrices and lists is affected by the following option
937 ;;; $DOALLMXOPS, $DOMXMXOPS, $DOMXPLUS, $DOSCMXOPPS, $DOSCMXPLUS, $LISTARITH
940 ;;; This routine should not be called directely. It is called by SIMPLIFYA.
941 ;;; A save access is to call the function ADD.
942 ;;;-----------------------------------------------------------------------------
944 (defun simplus (x w z
)
945 (prog (res check eqnflag matrixflag sumflag
)
946 (if (null (cdr x
)) (return 0))
950 (if (null x
) (go end
))
951 (setq w
(if z
(car x
) (simplifya (car x
) nil
)))
958 (not (member 'trunc
(cdar w
))))
959 (spsimpcases (cdr x
) w
))
960 (setq w
(ratdisrep w
))
965 (nconc (mapcar #'simplify
(cons w
(cdr x
)))
967 ((eq (caar w
) 'mequal
)
972 (add2 (cadr eqnflag
) (cadr w
))
973 (add2 (caddr eqnflag
) (caddr w
)))))
975 ((member (caar w
) '(mlist $matrix
))
977 (cond ((not matrixflag
) w
)
978 ((and (or $doallmxops $domxmxops $domxplus
979 (and (eq (caar w
) 'mlist
)
980 ($listp matrixflag
)))
981 (or (not (eq (caar w
) 'mlist
)) $listarith
))
982 (addmx matrixflag w
))
983 (t (setq res
(pls w res
)) matrixflag
)))
986 (setq sumflag t res
(sumpls w res
))
987 (setq w
(car res
) res
(cdr res
))))
988 (setq res
(pls w res
))
991 (setq res
(testp res
))
994 (cond ((and (or ($listp matrixflag
)
995 $doallmxops $doscmxplus $doscmxops
)
996 (or (not ($listp matrixflag
)) $listarith
))
997 (mxplusc res matrixflag
))
998 (t (testp (pls matrixflag
(pls res nil
)))))))
999 (setq res
(eqtest res check
))
1002 (add2 (cadr eqnflag
) res
)
1003 (add2 (caddr eqnflag
) res
))
1006 ;;;-----------------------------------------------------------------------------
1007 ;;; PLS (X OUT) 27.09.2010/DK
1009 ;;; Arguments and values:
1010 ;;; X - a Maxima expression or an atom
1011 ;;; OUT - a form ((mplus) <number> term1 term2 ...) or NIL
1012 ;;; result - a form ((mplus) <number> term1 ...), where x is added in.
1015 ;;; Adds the argument X into the form OUT. If OUT is NIL a form
1016 ;;; ((mplus) 0 X) is initialized, if X is an expression or a symbol,
1017 ;;; or ((mplus) X), if X is a number. Numbers are added to the first
1018 ;;; term <number> of the form. Any other symbol or expression is added
1019 ;;; into the canonical ordered list of arguments. The result is in a
1020 ;;; canonical order, but it is not a valid Maxima expression. To get a
1021 ;;; valid Maxima expression the result has to be checked with the
1022 ;;; function TESTP. This is done by the calling routine SIMPLUS.
1024 ;;; PLS checks the global flag *PLUSFLAG*, which is set in PLUSIN to T,
1025 ;;; if a mplus-expression is part of the result.
1028 ;;; (pls 2 nil) -> ((MPLUS) 2)
1029 ;;; (pls '$A nil) -> ((MPLUS) 0 $A)
1030 ;;; (pls '$B '((mplus) 0 $A)) -> ((MPLUS) 0 $A $B)
1031 ;;; (pls '$A '((mplus) 0 $A)) -> ((MPLUS) 0 ((MTIMES SIMP) 2 $A))
1033 ;;; Examples with the option variables $NUMER and $NEGDISTRIB:
1034 ;;; (let (($numer t)) (pls '$%e nil)) -> ((MPLUS) 2.718281828459045)
1035 ;;; (let (($negdistrib t)) (pls '((mtimes) -1 ((mplus) $A $B)) nil))
1036 ;;; -> ((MPLUS) 0 ((MTIMES SIMP) -1 $A) ((MTIMES SIMP) -1 $B))
1037 ;;; (let (($negdistrib nil)) (pls '((mtimes) -1 ((mplus) $A $B)) nil))
1038 ;;; -> ((MPLUS) 0 ((MTIMES) -1 ((MPLUS) $A $B)))
1041 ;;; The option variables $NUMER and $NEGDISTRIB and the global flag
1042 ;;; *PLUSFLAG*, which is set in the routine PLUSIN.
1045 ;;; PLUSIN and ADDK which are called from PLS and SIMPLUS.
1048 ;;; To add an expression into the list (CDR OUT), the list is passed
1049 ;;; to the routine PLUSIN as an argument. PLUSIN adds the argument to
1050 ;;; the list of terms by modifying the list (CDR OUT) destructively.
1051 ;;; The new value of OUT is returned as a result by PLS.
1052 ;;;-----------------------------------------------------------------------------
1054 ;; Set in PLUSIN to T to indicate a nested mplus expression.
1055 (defvar *plusflag
* nil
)
1057 ;; TESTP checks the result of PLS to get a valid Maxima mplus-expression.
1061 ((null (cddr x
)) (cadr x
))
1063 (cond ((null (cdddr x
)) (caddr x
)) (t (rplacd x
(cddr x
)))))
1067 (prog (fm *plusflag
*)
1068 (if (mtimesp x
) (setq x
(testtneg x
)))
1069 (when (and $numer
(atom x
) (eq x
'$%e
))
1070 ;; Replace $%e with its numerical value, when $numer ist TRUE
1073 ;; Initialize a form like ((mplus) <number> expr)
1076 (cond ((mnump x
) (ncons x
))
1078 (list 0 (cond ((atom x
) x
) (t (copy-list x
)))))
1079 ((mnump (cadr x
)) (copy-list (cdr x
) ))
1080 (t (cons 0 (copy-list (cdr x
) )))))))
1082 ;; Add a number into the first term of the list out.
1083 (return (cons '(mplus)
1084 (if (mnump (cadr out
))
1085 (cons (addk (cadr out
) x
) (cddr out
))
1086 (cons x
(cdr out
))))))
1087 ((not (mplusp x
)) (plusin x
(cdr out
)) (go end
)))
1088 ;; At this point we have a mplus expression as argument x. The following
1089 ;; code assumes that the argument x is already simplified and the terms
1090 ;; are in a canonical order.
1091 ;; First we add the number to the first term of the list out.
1093 (addk (if (mnump (cadr out
)) (cadr out
) 0)
1094 (cond ((mnump (cadr x
)) (setq x
(cdr x
)) (car x
)) (t 0))))
1095 ;; Initialize fm with the list of terms and start the loop to add the
1096 ;; terms of an mplus expression into the list out.
1099 (if (null (setq x
(cdr x
))) (go end
))
1100 ;; The return value of PLUSIN is a list, where the first element is the
1101 ;; added argument and the rest are the terms which follow the added
1103 (setq fm
(plusin (car x
) fm
))
1106 (if (not *plusflag
*) (return out
))
1107 (setq *plusflag
* nil
) ; *PLUSFLAG* T handles e.g. a+b+3*(a+b)-2*(a+b)
1109 ;; *PLUSFLAG* is set by PLUSIN to indicate that a mplus expression is
1110 ;; part of the result. For this case go again through the terms of the
1111 ;; result and add any term of the mplus expression into the list out.
1114 (when (mplusp (cadr fm
))
1116 (rplacd fm
(cddr fm
))
1120 (if (null (cdr fm
)) (return out
))
1123 ;;;-----------------------------------------------------------------------------
1124 ;;; PLUSIN (X FM) 27.09.2010/DK
1126 ;;; Arguments and values:
1127 ;;; X - a Maxima expression or atom
1128 ;;; FM - a list with the terms of an addition
1129 ;;; result - part of the list fm, which starts at the inserted expression
1132 ;;; Adds X into running list of additive terms FM. The routine modifies
1133 ;;; the argument FM destructively, but does not return the modified list as
1134 ;;; a result. The return value is a part of the list FM, which starts at the
1135 ;;; inserted term. PLUSIN can not handle Maxima numbers. PLUSIN is called
1136 ;;; only from the routine PLS.
1140 ;;; (plusin '$a fm) -> ($A)
1142 ;;; (plusin '$b fm) -> ($B)
1144 ;;; (plusin '$a fm) -> (((MTIMES SIMP) 2 $A) $B)
1145 ;;; fm -> (0 ((MTIMES SIMP) 2 $A) $B)
1148 ;;; Modifies destructively the argument FM, which contains the result of the
1149 ;;; addition of the argument X into the list FM.
1152 ;;; The option variables $doallmxops and $listarith.
1155 ;;; The return value is used in PLS to go in parallel through the list of
1156 ;;; terms, when adding a complete mplus-expression into the list of terms.
1157 ;;; This is triggered by the flag *PLUSFLAG*, which is set in PLUSIN, if
1158 ;;; a mplus-expression is added to the result list.
1159 ;;;-----------------------------------------------------------------------------
1161 (defun plusin (x fm
)
1162 (prog (x1 x2 flag check v w xnew a n m c
)
1167 (if (mnump (cadr x
)) (setq w
(cadr x
) x
(cddr x
))
1169 (t (setq x
(ncons x
))))
1170 (setq x1
(if (null (cdr x
)) (car x
) (cons '(mtimes) x
))
1171 xnew
(list* '(mtimes) w x
))
1173 (cond ((null (cdr fm
)))
1174 ((and (alike1 x1
(cadr fm
)) (null (cdr x
)))
1176 ;; Implement the simplification of
1177 ;; v*a^(c+n)+w*a^(c+m) -> (v*a^n+w*a^m)*a^c
1178 ;; where a, v, w, and (n-m) are integers.
1179 ((and (or (and (mexptp (setq x2
(cadr fm
)))
1182 (not (alike1 x1 x2
))
1184 (integerp (setq v
(cadr x2
)))
1185 (mexptp (setq x2
(caddr x2
)))))
1186 (integerp (setq a
(cadr x2
)))
1189 (integerp (sub (caddr x2
) (caddr x1
))))
1190 (setq n
(if (and (mplusp (caddr x2
))
1191 (mnump (cadr (caddr x2
))))
1193 (if (mnump (caddr x2
))
1196 (setq m
(if (and (mplusp (caddr x1
))
1197 (mnump (cadr (caddr x1
))))
1199 (if (mnump (caddr x1
))
1202 (setq c
(sub (caddr x2
) n
))
1205 ;; n and m are integers and the result is (v*a^n+w*a^m)*a^c.
1206 (setq x1
(mul (addk (timesk v
(exptb a n
))
1207 (timesk w
(exptb a m
)))
1211 ;; n and m are rational numbers: The difference n-m is an
1212 ;; integer. The rational numbers might be improper fractions.
1213 ;; The mixed numbers are: n = n1 + d1/r and m = n2 + d2/r,
1214 ;; where r is the common denominator. We have two cases:
1215 ;; I) d1 = d2: e.g. 2^(1/3+c)+2^(4/3+c)
1216 ;; The result is (v*a^n1+w*a^n2)*a^(c+d1/r)
1217 ;; II) d1 # d2: e.g. 2^(1/2+c)+2^(-1/2+c)
1218 ;; In this case one of the exponents d1 or d2 must
1219 ;; be negative. The negative exponent is factored out.
1220 ;; This guarantees that the factor (v*a^n1+w*a^n2)
1221 ;; is an integer. But the positive exponent has to be
1222 ;; adjusted accordingly. E.g. when we factor out
1223 ;; a^(d2/r) because d2 is negative, then we have to
1224 ;; adjust the positive exponent to n1 -> n1+(d1-d2)/r.
1226 ;; Part of the simplification is done in simptimes. E.g.
1227 ;; this algorithm simplifies the sum sqrt(2)+3*sqrt(2)
1228 ;; to 4*sqrt(2). In simptimes this is further simplified
1230 (multiple-value-bind (n1 d1
)
1231 (truncate (num1 n
) (denom1 n
))
1232 (multiple-value-bind (n2 d2
)
1233 (truncate (num1 m
) (denom1 m
))
1234 (cond ((equal d1 d2
)
1235 ;; Case I: -> (v*a^n1+w*a^n2)*a^(c+d1/r)
1237 (mul (addk (timesk v
(exptb a n1
))
1238 (timesk w
(exptb a n2
)))
1241 (div d1
(denom1 n
))))))
1244 ;; Case II:: d2 is negative, adjust n1.
1245 (setq n1
(add n1
(div (sub d1 d2
) (denom1 n
))))
1247 (mul (addk (timesk v
(exptb a n1
))
1248 (timesk w
(exptb a n2
)))
1251 (div d2
(denom1 n
))))))
1254 ;; Case II: d1 is negative, adjust n2.
1255 (setq n2
(add n2
(div (sub d2 d1
) (denom1 n
))))
1257 (mul (addk (timesk v
(exptb a n1
))
1258 (timesk w
(exptb a n2
)))
1261 (div d1
(denom1 n
))))))
1263 ;; This clause should never be reached.
1264 (t (merror "Internal error in simplus."))))))))
1265 ((mtimesp (cadr fm
))
1266 (cond ((alike1 x1
(cadr fm
))
1268 ((and (mnump (cadadr fm
)) (alike x
(cddadr fm
)))
1269 (setq flag t
) ; found common factor
1271 ((great xnew
(cadr fm
)) (go gr
))))
1272 ((great x1
(cadr fm
)) (go gr
)))
1273 (setq xnew
(eqtest (testt xnew
) (or check
'((foo)))))
1274 (return (cdr (rplacd fm
(cons xnew
(cdr fm
)))))
1281 (list* '(mtimes simp
) 0 x
)
1282 ;; Call muln to get a simplified product.
1283 (if (mtimesp (setq x1
(muln (cons (addk 1 w
) x
) t
)))
1287 (cond ((not (mtimesp (cadr fm
)))
1290 ;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
1291 (rplacd (cadr fm
) (cddadr fm
))
1293 ((not (zerop1 (cadadr fm
)))
1295 ;; Handle the multiplication with a zero.
1296 ((and (or (not $listarith
) (not $doallmxops
))
1297 (mxorlistp (caddr (cadr fm
))))
1299 (cons (constmx 0 (caddr (cadr fm
))) (cddr fm
))))))
1300 ;; (cadadr fm) is zero. If the first term of fm is a number,
1301 ;; add it to preserve the type.
1302 (when (mnump (car fm
))
1303 (rplaca fm
(addk (car fm
) (cadadr fm
))))
1304 (return (rplacd fm
(cddr fm
)))
1306 ;; Call muln to get a simplified product.
1307 (setq x1
(muln (cons (addk w
(if flag
(cadadr fm
) 1)) x
) t
))
1308 ;; Make a mplus expression to guarantee that x1 is added again into the sum
1309 (setq x1
(list '(mplus) x1
))
1313 (list* '(mtimes) x1 x
)
1314 (if (mtimesp x1
) (testtneg x1
) x1
)))
1315 (if (not (mtimesp (cadr fm
))) (go check
))
1316 (when (and (onep (cadadr fm
)) flag
(null (cdddr (cadr fm
))))
1317 ;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
1318 (rplaca (cdr fm
) (caddr (cadr fm
))) (go check
))
1321 (if (mplusp (cadr fm
)) (setq *plusflag
* t
)) ; A nested mplus expression
1324 ;;;-----------------------------------------------------------------------------
1326 ;; Routines to add matrices
1328 (defun mxplusc (sc mx
)
1330 (setq sc
(partition-ns (cdr sc
)))
1331 (cond ((null (car sc
)) (cons '(mplus) (cons mx
(cadr sc
))))
1332 ((not (null (cadr sc
)))
1335 (outermap1 'mplus
(cons '(mplus) (car sc
)) mx
))
1337 (t (simplify (outermap1 'mplus
(cons '(mplus) (car sc
)) mx
)))))
1338 ((not (scalar-or-constant-p sc $assumescalar
))
1339 (testp (pls mx
(pls sc nil
))))
1340 (t (simplify (outermap1 'mplus sc mx
)))))
1342 (defun partition-ns (x)
1343 (let (sp nsp
) ; SP = scalar part, NSP = nonscalar part
1344 (mapc #'(lambda (z) (if (scalar-or-constant-p z $assumescalar
)
1345 (setq sp
(cons z sp
))
1346 (setq nsp
(cons z nsp
))))
1348 (list (nreverse sp
) (nreverse nsp
))))
1350 (defun addmx (x1 x2
)
1351 (let (($doscmxops t
) ($domxmxops t
) ($listarith t
))
1352 (simplify (fmapl1 'mplus x1 x2
))))
1354 ;;; ----------------------------------------------------------------------------
1356 ;;; Simplification of the Log function
1358 ;; The log function distributes over lists, matrices, and equations
1359 (defprop %log
(mlist $matrix mequal
) distribute_over
)
1361 (defun simpln (x y z
)
1363 (setq y
(simpcheck (cadr x
) z
))
1364 (cond ((onep1 y
) (addk -
1 y
))
1366 (cond (radcanp (list '(%log simp
) 0))
1368 (merror (intl:gettext
"log: encountered log(0).")))
1369 (t (throw 'errorsw t
))))
1370 ;; Check evaluation in floating point precision.
1371 ((flonum-eval (mop x
) y
))
1372 ;; Check evaluation in bigfloag precision.
1373 ((and (not (member 'simp
(car x
)))
1374 (big-float-eval (mop x
) y
)))
1377 (cond ((or (and $logexpand
(eq $domain
'$real
))
1378 (member $logexpand
'($all $super
))
1379 (and (eq ($csign
(cadr y
)) '$pos
)
1380 (not (member ($csign
(caddr y
))
1381 '($complex $imaginary
)))))
1382 ;; Simplify log(x^a) -> a*log(x), where x > 0 and a is real
1383 (mul (caddr y
) (take '(%log
) (cadr y
))))
1384 ((or (and (ratnump (caddr y
))
1385 (or (eql 1 (cadr (caddr y
)))
1386 (eql -
1 (cadr (caddr y
)))))
1387 (maxima-integerp (inv (caddr y
))))
1388 ;; Simplify log(z^(1/n)) -> log(z)/n, where n is an integer
1390 (take '(%log
) (cadr y
))))
1391 ((and (eq (cadr y
) '$%e
)
1392 (or (not (member ($csign
(caddr y
))
1393 '($complex $imaginary
)))
1394 (not (member ($csign
(mul '$%i
(caddr y
)))
1395 '($complex $imaginary
)))))
1396 ;; Simplify log(exp(x)) and log(exp(%i*x)), where x is a real
1398 (t (eqtest (list '(%log
) y
) x
))))
1400 ;; Simplify log(n/d)
1401 (cond ((eql (cadr y
) 1)
1402 (mul -
1 (take '(%log
) (caddr y
))))
1403 ((eq $logexpand
'$super
)
1404 (sub (take '(%log
) (cadr y
)) (take '(%log
) (caddr y
))))
1405 (t (eqtest (list '(%log
) y
) x
))))
1406 ((and (member $logexpand
'($all $super
))
1408 (do ((y (cdr y
) (cdr y
))
1410 ((null y
) (return (addn b t
)))
1411 (setq b
(cons (take '(%log
) (car y
)) b
))))
1412 ((and (member $logexpand
'($all $super
))
1414 (member (caar y
) '(%product $product
)))
1415 (let ((new-op (if (char= (get-first-char (caar y
)) #\%
) '%sum
'$sum
)))
1416 (simplifya `((,new-op
) ((%log
) ,(cadr y
)) ,@(cddr y
)) t
)))
1419 (eq ($sign y
) '$neg
))
1420 (add (mul '$%i
'$%pi
) (take '(%log
) (neg y
))))
1421 ((taylorize (mop x
) (second x
)))
1422 (t (eqtest (list '(%log
) y
) x
))))
1425 (simplifya (list '(mtimes) (caddr w
)
1426 (simplifya (list '(%log
) (cadr w
)) t
)) t
))
1428 ;;; ----------------------------------------------------------------------------
1430 ;;; Implementation of the Square root function
1432 (defprop $sqrt %sqrt verb
)
1433 (defprop $sqrt %sqrt alias
)
1435 (defprop %sqrt $sqrt noun
)
1436 (defprop %sqrt $sqrt reversealias
)
1438 (defprop %sqrt simp-sqrt operators
)
1441 (simplify (list '(%sqrt
) z
)))
1443 (defun simp-sqrt (x ignored z
)
1444 (declare (ignore ignored
))
1446 (simplifya (list '(mexpt) (cadr x
) '((rat simp
) 1 2)) z
))
1448 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1450 ;;; Simplification of the "/" operator.
1452 (defun simpquot (x y z
)
1454 (cond ((and (integerp (cadr x
)) (integerp (caddr x
)) (not (zerop (caddr x
))))
1455 (*red
(cadr x
) (caddr x
)))
1456 ((and (numberp (cadr x
)) (numberp (caddr x
)) (not (zerop (caddr x
))))
1457 (/ (cadr x
) (caddr x
)))
1458 ((and (floatp (cadr x
)) (floatp (caddr x
)) #-ieee-floating-point
(not (zerop (caddr x
))))
1459 (/ (cadr x
) (caddr x
)))
1460 ((and ($bfloatp
(cadr x
)) ($bfloatp
(caddr x
)) (not (equal bigfloatzero
(caddr x
))))
1461 ;; Call BIGFLOATP to ensure that arguments have same precision.
1462 ;; Otherwise FPQUOTIENT could return a spurious value.
1463 (bcons (fpquotient (cdr (bigfloatp (cadr x
))) (cdr (bigfloatp (caddr x
))))))
1464 (t (setq y
(simplifya (cadr x
) z
))
1465 (setq x
(simplifya (list '(mexpt) (caddr x
) -
1) z
))
1466 (if (equal y
1) x
(simplifya (list '(mtimes) y x
) t
)))))
1468 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1470 ;;; Implementation of the abs function.
1472 ;; Put the properties alias, reversealiases, noun and verb on the property list.
1473 (defprop $abs mabs alias
)
1474 (defprop $abs mabs verb
)
1475 (defprop mabs $abs reversealias
)
1476 (defprop mabs $abs noun
)
1478 ;; The abs function distributes over bags.
1479 (defprop mabs
(mlist $matrix mequal
) distribute_over
)
1481 ;; Define a verb function $abs
1483 (simplify (list '(mabs) x
)))
1485 ;; The abs function is a simplifying function.
1486 (defprop mabs simpabs operators
)
1488 (defun simpabs (e y z
)
1489 (declare (ignore y
))
1492 (x (simpcheck (second e
) z
)))
1494 (cond ((complex-number-p x
#'(lambda (s) (or (floatp s
) ($bfloatp s
))))
1495 (maxima::to
(bigfloat::abs
(bigfloat:to x
))))
1497 ((complex-number-p x
#'mnump
)
1500 ;; nounform for arrays...
1501 ((or (arrayp x
) ($member x $arrays
)) `((mabs simp
) ,x
))
1503 ;; taylor polynomials
1504 ((taylorize 'mabs x
))
1506 ;; values for extended real arguments:
1507 ((member x
'($inf $infinity $minf
)) '$inf
)
1508 ((member x
'($ind $und
)) x
)
1510 ;; abs(abs(expr)) --> abs(expr). Since x is simplified, it's OK to return x.
1511 ((and (consp x
) (consp (car x
)) (eq (caar x
) 'mabs
))
1514 ;; abs(conjugate(expr)) = abs(expr).
1515 ((and (consp x
) (consp (car x
)) (eq (caar x
) '$conjugate
))
1516 (take '(mabs) (cadr x
)))
1519 (setq sgn
($csign x
))
1520 (cond ((member sgn
'($neg $nz
)) (mul -
1 x
))
1521 ((eq '$zero sgn
) (mul 0 x
))
1522 ((member sgn
'($pos $pz
)) x
)
1524 ;; for complex constant expressions, use $cabs
1525 ((and (eq sgn
'$complex
) ($constantp x
))
1528 ;; abs(pos^complex) --> pos^(realpart(complex)).
1529 ((and (eq sgn
'$complex
) (mexptp x
) (eq '$pos
($csign
(second x
))))
1530 (power (second x
) ($realpart
(third x
))))
1532 ;; for abs(neg^z), use cabs.
1533 ((and (mexptp x
) (eq '$neg
($csign
(second x
))))
1536 ;; When x # 0, we have abs(signum(x)) = 1.
1537 ((and (eq '$pn sgn
) (consp x
) (consp (car x
)) (eq (caar x
) '%signum
)) 1)
1539 ;; multiplicative property: abs(x*y) = abs(x) * abs(y). We would like
1540 ;; assume(a*b > 0), abs(a*b) --> a*b. Thus the multiplicative property
1541 ;; is applied after the sign test.
1543 (muln (mapcar #'(lambda (u) (take '(mabs) u
)) (margs x
)) t
))
1545 ;; abs(x^n) = abs(x)^n for integer n. Is the featurep check worthwhile?
1546 ;; Again the sign check is done first because we'd like abs(x^2) --> x^2.
1547 ((and (mexptp x
) ($featurep
(caddr x
) '$integer
))
1548 (power (take '(mabs) (cadr x
)) (caddr x
)))
1550 ;; Reflection rule: abs(-x) --> abs(x)
1551 ((great (neg x
) x
) (take '(mabs) (neg x
)))
1554 (t (eqtest (list '(mabs) x
) e
)))))))
1556 (defun abs-integral (x)
1557 (mul (div 1 2) x
(take '(mabs) x
)))
1559 (putprop 'mabs
`((x) ,'abs-integral
) 'integral
)
1561 ;; I (rtoy) think this does some simple optimizations of x * y.
1566 ;; We have something like ((mtimes) foo). This is the same as foo.
1569 ;; We have 1*foo. Which is the same as foo. This should not
1570 ;; be applied to 1.0 or 1b0!
1571 (cond ((null (cdddr x
))
1573 (t (rplacd x
(cddr x
)))))
1577 ;; This basically converts -(a+b) to -a-b.
1579 (cond ((and (equal (cadr x
) -
1)
1583 ;; If x is exactly of the form -1*(sum), and $negdistrib is
1584 ;; true, we distribute the -1 across the sum.
1585 (addn (mapcar #'(lambda (z)
1591 ;; Simplification of the "-" operator
1592 (defun simpmin (x vestigial z
)
1593 (declare (ignore vestigial
))
1594 (cond ((null (cdr x
)) 0)
1596 (mul -
1 (simplifya (cadr x
) z
)))
1598 ;; ((mminus) a b ...) -> ((mplus) a ((mtimes) -1 b) ...)
1599 (sub (simplifya (cadr x
) z
) (addn (cddr x
) z
)))))
1601 (defun simptimes (x w z
) ; W must be 1
1602 (prog (res check eqnflag matrixflag sumflag
)
1603 (if (null (cdr x
)) (return 1))
1609 (cond ((and matrixflag
(mxorlistp1 matrixflag
))
1610 (return (constmx res matrixflag
)))
1611 (eqnflag (return (list '(mequal simp
)
1612 (mul2 res
(cadr eqnflag
))
1613 (mul2 res
(caddr eqnflag
)))))
1616 (cond ((mxorlistp u
)
1617 (return (setq res
(constmx res u
))))
1619 (mxorlistp1 (cadr u
))
1620 ($numberp
(caddr u
)))
1621 (return (setq res
(constmx res
(cadr u
)))))
1625 (list '(mequal simp
)
1627 (mul2 res
(caddr u
))))))))))))
1629 ((null x
) (go end
)))
1630 (setq w
(if z
(car x
) (simplifya (car x
) nil
)))
1632 (cond ((atom w
) nil
)
1633 ((eq (caar w
) 'mrat
)
1634 (cond ((or eqnflag matrixflag
1636 (not (member 'trunc
(cdar w
))))
1637 (spsimpcases (cdr x
) w
))
1638 (setq w
(ratdisrep w
))
1642 (ratf (cons '(mtimes)
1643 (nconc (mapcar #'simplify
(cons w
(cdr x
)))
1645 ((eq (caar w
) 'mequal
)
1650 (mul2 (cadr eqnflag
) (cadr w
))
1651 (mul2 (caddr eqnflag
) (caddr w
)))))
1653 ((member (caar w
) '(mlist $matrix
))
1655 (cond ((not matrixflag
) w
)
1656 ((and (or $doallmxops $domxmxops $domxtimes
)
1657 (or (not (eq (caar w
) 'mlist
)) $listarith
)
1658 (not (eq *inv
* '$detout
)))
1659 (stimex matrixflag w
))
1660 (t (setq res
(tms (copy-tree w
) 1 (copy-tree res
))) matrixflag
)))
1662 ((and (eq (caar w
) '%sum
) $sumexpand
)
1663 (setq sumflag
(sumtimes sumflag w
))
1665 (setq res
(tms (copy-tree w
) 1 (copy-tree res
)))
1668 (cond ((mtimesp res
) (setq res
(testt res
))))
1669 (cond (sumflag (setq res
(cond ((or (null res
) (equal res
1)) sumflag
)
1670 ((not (mtimesp res
))
1671 (list '(mtimes) res sumflag
))
1672 (t (nconc res
(list sumflag
)))))))
1673 (cond ((or (atom res
)
1674 (not (member (caar res
) '(mexpt mtimes
)))
1675 (and (zerop $expop
) (zerop $expon
))
1677 ((eq (caar res
) 'mtimes
) (setq res
(expandtimes res
)))
1678 ((and (mplusp (cadr res
))
1679 (fixnump (caddr res
))
1680 (not (or (> (caddr res
) $expop
)
1681 (> (- (caddr res
)) $expon
))))
1682 (setq res
(expandexpt (cadr res
) (caddr res
)))))
1685 (cond ((null res
) matrixflag
)
1686 ((and (or ($listp matrixflag
)
1689 (not (member res
'(-1 -
1.0))))
1690 ;; RES should only be -1 here (not = 1)
1692 (member res
'(-1 -
1.0))))
1693 (or (not ($listp matrixflag
)) $listarith
))
1694 (mxtimesc res matrixflag
))
1695 (t (testt (tms matrixflag
1 (tms res
1 nil
))))))))
1696 (if res
(setq res
(eqtest res check
)))
1697 (return (cond (eqnflag
1698 (if (null res
) (setq res
1))
1700 (mul2 (cadr eqnflag
) res
)
1701 (mul2 (caddr eqnflag
) res
)))
1704 (defun spsimpcases (l e
)
1706 (if (or (mbagp u
) (and (not (atom u
))
1708 (not (member 'trunc
(cdar e
)))))
1711 (defun mxtimesc (sc mx
)
1713 (and (mtimesp sc
) (member (cadr sc
) '(-1 -
1.0))
1714 $doscmxops
(not (or $doallmxops $domxmxops $domxtimes
))
1715 (setq sign
(cadr sc
)) (rplaca (cdr sc
) nil
))
1716 (setq out
(let ((scp* (cond ((mtimesp sc
) (partition-ns (cdr sc
)))
1717 ((not (scalar-or-constant-p sc $assumescalar
))
1720 (cond ((null scp
*) (list '(mtimes simp
) sc mx
))
1721 ((and (not (atom scp
*)) (null (car scp
*)))
1722 (append '((mtimes)) (cadr scp
*) (list mx
)))
1723 ((or (atom scp
*) (and (null (cdr scp
*))
1724 (not (null (cdr sc
)))
1725 (setq scp
* (cons '(mtimes) (car scp
*))))
1727 (simplifya (outermap1 'mtimes scp
* mx
) nil
))
1728 (t (append '((mtimes))
1731 (cons '(mtimes) (car scp
*)) mx
)
1734 (cond (sign (if (mtimesp out
)
1735 (rplacd out
(cons sign
(cdr out
)))
1736 (list '(mtimes) sign out
)))
1737 ((mtimesp out
) (testt out
))
1741 (let (($doscmxops t
) ($domxmxops t
) ($listarith t
))
1742 (simplify (fmapl1 'mtimes x y
))))
1744 ;; TMS takes a simplified expression FACTOR and a cumulative
1745 ;; PRODUCT as arguments and modifies the cumulative product so
1746 ;; that the expression is now one of its factors. The
1747 ;; exception to this occurs when a tellsimp rule is triggered.
1748 ;; The second argument is the POWER to which the expression is
1749 ;; to be raised within the product.
1751 (defun tms (factor power product
&aux tem
)
1752 (let ((*rulesw
* nil
)
1754 (when (mplusp product
) (setq product
(list '(mtimes simp
) product
)))
1755 (cond ((zerop1 factor
)
1756 (cond ((mnegp power
)
1759 (merror (intl:gettext
"Division by 0"))))
1761 ((and (null product
)
1762 (or (and (mtimesp factor
) (equal power
1))
1763 (and (setq product
(list '(mtimes) 1)) nil
)))
1764 (setq tem
(append '((mtimes)) (if (mnump (cadr factor
)) nil
'(1))
1766 (if (= (length tem
) 1)
1767 (setq tem
(copy-list tem
))
1770 (do ((factor-list (cdr factor
) (cdr factor-list
)))
1771 ((or (null factor-list
) (zerop1 product
)) product
)
1772 (setq z
(timesin (car factor-list
) (cdr product
) power
))
1775 (setq product
(tms-format-product z
)))))
1777 (setq z
(timesin factor
(cdr product
) power
))
1779 (tms-format-product z
)
1782 (defun tms-format-product (x)
1783 (cond ((zerop1 x
) x
)
1784 ((mnump x
) (list '(mtimes) x
))
1785 ((not (mtimesp x
)) (list '(mtimes) 1 x
))
1786 ((not (mnump (cadr x
))) (cons '(mtimes) (cons 1 (cdr x
))))
1790 (cond ($ratsimpexpons
(sratsimp (list '(mplus) x y
)))
1791 ((and (mnump x
) (mnump y
)) (addk x y
))
1795 (if (and (mnump x
) (mnump y
))
1799 (defun simp-limit (x vestigial z
)
1800 (declare (ignore vestigial
))
1801 (let ((l1 (length x
))
1803 (unless (or (= l1
2) (= l1
4) (= l1
5))
1804 (merror (intl:gettext
"limit: wrong number of arguments.")))
1805 (setq y
(simpmap (cdr x
) z
))
1806 (cond ((and (= l1
5) (not (member (cadddr y
) '($plus $minus
))))
1807 (merror (intl:gettext
"limit: direction must be either 'plus' or 'minus': ~M") (cadddr y
)))
1809 (merror (intl:gettext
"limit: variable must not be a number; found: ~M") (cadr y
)))
1813 (eqtest (cons '(%limit
) y
) x
)))))
1815 (defun simpinteg (x vestigial z
)
1816 (declare (ignore vestigial
))
1817 (let ((l1 (length x
))
1819 (unless (or (= l1
3) (= l1
5))
1820 (merror (intl:gettext
"integrate: wrong number of arguments.")))
1821 (setq y
(simpmap (cdr x
) z
))
1822 (cond ((mnump (cadr y
))
1823 (merror (intl:gettext
"integrate: variable must not be a number; found: ~M") (cadr y
)))
1824 ((and (= l1
5) (alike1 (caddr y
) (cadddr y
)))
1827 (free (setq z
(sub (cadddr y
) (caddr y
))) '$%i
)
1828 (eq ($sign z
) '$neg
))
1829 (neg (simplifya (list '(%integrate
) (car y
) (cadr y
) (cadddr y
) (caddr y
)) t
)))
1833 (if (or (among '$inf z
) (among '$minf z
))
1837 (eqtest (cons '(%integrate
) y
) x
)))))
1839 (defun simpbigfloat (x vestigial simp-flag
)
1840 (declare (ignore vestigial simp-flag
))
1843 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1845 ;;; Implementation of the Exp function.
1847 (defprop $exp %exp verb
)
1848 (defprop $exp %exp alias
)
1850 (defprop %exp $exp noun
)
1851 (defprop %exp $exp reversealias
)
1853 (defprop %exp simp-exp operators
)
1856 (simplify (list '(%exp
) z
)))
1858 ;; Support a function for code,
1859 ;; which depends on an unsimplified noun form.
1860 (defmfun $exp-form
(z)
1861 (list '(mexpt) '$%e z
))
1863 (defun simp-exp (x ignored z
)
1864 (declare (ignore ignored
))
1866 (simplifya (list '(mexpt) '$%e
(cadr x
)) z
))
1868 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1870 (defun simplambda (x vestigial simp-flag
)
1871 (declare (ignore vestigial simp-flag
))
1872 ; Check for malformed lambda expressions.
1873 ; We verify that we have a valid list of parameters and a non-empty body.
1874 (let ((params (cadr x
)))
1875 (unless ($listp params
)
1876 (merror (intl:gettext
"lambda: first argument must be a list; found: ~M") params
))
1877 (do ((params (cdr params
) (cdr params
))
1880 (when (mdeflistp params
)
1881 (setq params
(cdar params
)))
1882 (let ((p (car params
)))
1883 (unless (or (mdefparam p
)
1884 (and (op-equalp p
'mquote
)
1885 (mdefparam (cadr p
))))
1886 (merror (intl:gettext
"lambda: parameter must be a symbol and must not be a system constant; found: ~M") p
))
1888 (when (member p seen-params
:test
#'eq
)
1889 (merror (intl:gettext
"lambda: ~M occurs more than once in the parameter list") p
))
1890 (push p seen-params
))))
1891 (when (null (cddr x
))
1892 (merror (intl:gettext
"lambda: no body present")))
1893 (cons '(lambda simp
) (cdr x
)))
1895 (defun simpmdef (x vestigial simp-flag
)
1896 (declare (ignore vestigial simp-flag
))
1898 (cons '(mdefine simp
) (cdr x
)))
1900 (defun simpmap (e z
)
1901 (mapcar #'(lambda (u) (simpcheck u z
)) e
))
1904 (let ((x ($expand e
1 1)))
1905 (cond ((or (not (free x
'$ind
)) (not (free x
'$und
))
1906 (not (free x
'$zeroa
)) (not (free x
'$zerob
))
1907 (not (free x
'$infinity
))
1910 ((and (free x
'$inf
) (free x
'$minf
)) x
)
1911 (t (infsimp1 x e
)))))
1913 (defun infsimp1 (x e
)
1914 (let ((minf-coef (coeff x
'$minf
1))
1915 (inf-coef (coeff x
'$inf
1)))
1916 (cond ((or (and (equal minf-coef
0)
1918 (and (not (free minf-coef
'$inf
))
1919 (not (free inf-coef
'$minf
)))
1920 (let ((new-exp (sub (add2 (mul2 minf-coef
'$minf
)
1921 (mul2 inf-coef
'$inf
))
1923 (and (not (free new-exp
'$inf
))
1924 (not (free new-exp
'$minf
)))))
1926 (t (let ((sign-minf-coef ($asksign minf-coef
))
1927 (sign-inf-coef ($asksign inf-coef
)))
1928 (cond ((or (and (eq sign-inf-coef
'$zero
)
1929 (eq sign-minf-coef
'$neg
))
1930 (and (eq sign-inf-coef
'$pos
)
1931 (eq sign-minf-coef
'$zero
))
1932 (and (eq sign-inf-coef
'$pos
)
1933 (eq sign-minf-coef
'$neg
))) '$inf
)
1934 ((or (and (eq sign-inf-coef
'$zero
)
1935 (eq sign-minf-coef
'$pos
))
1936 (and (eq sign-inf-coef
'$neg
)
1937 (eq sign-minf-coef
'$zero
))
1938 (and (eq sign-inf-coef
'$neg
)
1939 (eq sign-minf-coef
'$pos
))) '$minf
)
1940 ((or (and (eq sign-inf-coef
'$pos
)
1941 (eq sign-minf-coef
'$pos
))
1942 (and (eq sign-inf-coef
'$neg
)
1943 (eq sign-minf-coef
'$neg
))) '$und
)))))))
1945 (defun infsimp2 (x e
)
1947 (if (isinop x
'%limit
) e x
))
1949 (defun simpderiv (x y z
)
1951 (cond ((not (even (length x
)))
1952 (cond ((and (cdr x
) (null (cdddr x
))) (nconc x
'(1)))
1953 (t (wna-err '%derivative
)))))
1954 (setq w
(cons '(%derivative
) (simpmap (cdr x
) z
)))
1956 (do ((u (cddr w
) (cddr u
))) ((null u
))
1957 (cond ((mnump (car u
))
1958 (merror (intl:gettext
"diff: variable must not be a number; found: ~M") (car u
)))))
1959 (cond ((or (zerop1 y
)
1960 (and (or (mnump y
) (and (atom y
) (constant y
)))
1962 (and (not (alike1 y
(caddr w
)))
1963 (do ((u (cddr w
) (cddr u
))) ((null u
))
1964 (cond ((and (numberp (cadr u
))
1965 (not (zerop (cadr u
))))
1968 ((and (not (atom y
)) (eq (caar y
) '%derivative
) derivsimp
)
1969 (rplacd w
(append (cdr y
) (cddr w
)))))
1971 (return (if (null derivflag
) (list '(%del simp
) y
) (deriv (cdr w
)))))
1974 (cond ((null u
) (go next
))
1975 ((zerop1 (caddr u
)) (rplacd u
(cdddr u
)))
1976 (t (setq u
(cddr u
))))
1979 (cond ((null (cddr w
)) (return y
))
1980 ((and (null (cddddr w
))
1982 (alike1 (cadr w
) (caddr w
)))
1987 (cond ((null (cddr z
)) (go loop
))
1988 ((alike1 (car z
) (caddr z
))
1989 (rplaca (cdddr z
) (add2 (cadr z
) (cadddr z
)))
1990 (rplacd z
(cdddr z
)))
1991 ((great (car z
) (caddr z
))
1992 (let ((u1 (car z
)) (u2 (cadr z
)) (v1 (caddr z
)) (v2 (cadddr z
)))
1995 (rplacd z
(cons v2
(cons u1
(cons u2
(cddddr z
))))))))
1996 (cond ((setq z
(cddr z
)) (go sort
)))
1998 (cond ((null flag
) (return (cond ((null derivflag
) (eqtest w x
))
1999 (t (deriv (cdr w
)))))))
2005 (setq x
(num1 x
)) (cond ((plusp x
) 1) ((minusp x
) -
1) (t 0)))
2007 ((mplusp x
) (if expandp
1 (signum1 (car (last x
)))))
2008 ((mtimesp x
) (if (mplusp (cadr x
)) 1 (signum1 (cadr x
))))
2011 (defprop %signum
(mlist $matrix mequal
) distribute_over
)
2013 (defun simpsignum (e y z
)
2014 (declare (ignore y
))
2016 (let ((x (simpcheck (second e
) z
)) (sgn))
2018 (cond ((complex-number-p x
#'mnump
)
2019 (if (complex-number-p x
#'$ratnump
) ;; nonfloat complex
2020 (if (zerop1 x
) 0 ($rectform
(div x
($cabs x
))))
2021 (maxima::to
(bigfloat::signum
(bigfloat::to x
)))))
2023 ;; idempotent: signum(signum(z)) = signum(z).
2024 ((and (consp x
) (consp (car x
)) (eq '%signum
(mop x
))) x
)
2027 (setq sgn
($csign x
))
2028 (cond ((eq sgn
'$neg
) -
1)
2032 ;; multiplicative: signum(ab) = signum(a) * signum(b).
2034 (muln (mapcar #'(lambda (s) (take '(%signum
) s
)) (margs x
)) t
))
2036 ;; Reflection rule: signum(-x) --> -signum(x).
2037 ((great (neg x
) x
) (neg (take '(%signum
) (neg x
))))
2040 (t (eqtest (list '(%signum
) x
) e
)))))))
2042 (defun exptrl (r1 r2
)
2043 (cond ((equal r2
1) r1
)
2045 (cond ((mnump r1
) (addk 0.0 r1
))
2046 ;; Do not simplify the type of the number away.
2047 (t (list '(mexpt simp
) r1
1.0))))
2048 ((equal r2 bigfloatone
)
2049 (cond ((mnump r1
) ($bfloat r1
))
2050 ;; Do not simplify the type of the number away.
2051 (t (list '(mexpt simp
) r1 bigfloatone
))))
2053 (cond ((or (zerop1 r2
) (mnegp r2
))
2055 (merror (intl:gettext
"expt: undefined: ~M") (list '(mexpt) r1 r2
))
2056 (throw 'errorsw t
)))
2057 (t (zerores r1 r2
))))
2058 ((or (zerop1 r2
) (onep1 r1
))
2059 (cond ((or ($bfloatp r1
) ($bfloatp r2
)) bigfloatone
)
2060 ((or (floatp r1
) (floatp r2
)) 1.0)
2062 ((or ($bfloatp r1
) ($bfloatp r2
)) ($bfloat
(list '(mexpt) r1 r2
)))
2063 ((and (numberp r1
) (integerp r2
)) (exptb r1 r2
))
2064 ((and (numberp r1
) (floatp r2
) (equal r2
(float (floor r2
))))
2065 (exptb (float r1
) (floor r2
)))
2066 ((or $numer
(and (floatp r2
) (or (plusp (num1 r1
)) $numer_pbranch
)))
2067 (let (y #+kcl
(r1 r1
) #+kcl
(r2 r2
))
2068 (cond ((minusp (setq r1
(addk 0.0 r1
)))
2069 (cond ((or $numer_pbranch
(eq $domain
'$complex
))
2071 ;; R1^R2 = (-R1)^R2*cos(pi*R2) + i*(-R1)^R2*sin(pi*R2)
2072 (setq r2
(addk 0.0 r2
))
2073 (setq y
(exptrl (- r1
) r2
) r2
(* %pi-val r2
))
2074 (add2 (* y
(cos r2
))
2075 (list '(mtimes simp
) (* y
(sin r2
)) '$%i
)))
2076 (t (setq y
(let ($numer $float $keepfloat $ratprint
)
2078 (mul2 y
(exptrl (- r1
) r2
)))))
2079 ((equal (setq r2
(addk 0.0 r2
)) (float (floor r2
)))
2080 (exptb r1
(floor r2
)))
2081 ((and (equal (setq y
(* 2.0 r2
)) (float (floor y
)))
2082 (not (equal r1 %e-val
)))
2083 (exptb (sqrt r1
) (floor y
)))
2084 (t (exp (* r2
(log r1
)))))))
2085 ((floatp r2
) (list '(mexpt simp
) r1 r2
))
2088 (exptrl (cond ((equal (abs (cadr r1
)) 1)
2089 (* (cadr r1
) (caddr r1
)))
2090 ;; We set the simp flag at this place. This
2091 ;; changes nothing for an exponent r2 # -1.
2092 ;; exptrl is called again and does not look at
2093 ;; the simp flag. For the case r2 = -1 exptrl
2094 ;; is called with an exponent 1. For this case
2095 ;; the base is immediately returned. Now the
2096 ;; base has the correct simp flag. (DK 02/2010)
2098 (list '(rat simp
) (- (caddr r1
)) (- (cadr r1
))))
2099 (t (list '(rat simp
) (caddr r1
) (cadr r1
))))
2101 (t (list '(rat simp
) (exptb (cadr r1
) r2
) (exptb (caddr r1
) r2
)))))
2102 ((and (floatp r1
) (alike1 r2
'((rat) 1 2)))
2104 (list '(mtimes simp
) (sqrt (- r1
)) '$%i
)
2106 ((and (floatp r1
) (alike1 r2
'((rat) -
1 2)))
2108 (list '(mtimes simp
) (/ -
1.0 (sqrt (- r1
))) '$%i
)
2112 (exptrl r1
(fpcofrat r2
))
2113 (mul2 (exptrl -
1 r2
) ;; (-4.5)^(1/4) -> (4.5)^(1/4) * (-1)^(1/4)
2114 (exptrl (- r1
) r2
))))
2115 (exptrlsw (list '(mexpt simp
) r1 r2
))
2118 (simptimes (list '(mtimes)
2119 (exptrl r1
(truncate (cadr r2
) (caddr r2
)))
2120 (let ((y (let ($keepfloat $ratprint
)
2121 (simpnrt r1
(caddr r2
))))
2122 (z (rem (cadr r2
) (caddr r2
))))
2124 (list (car y
) (cadr y
) (mul2 (caddr y
) z
))
2128 (defun simpexpt (x y z
)
2129 (prog (gr pot check res
*rulesw
* w mlpgr mlppot
)
2131 (cond (z (setq gr
(cadr x
) pot
(caddr x
)) (go cont
)))
2133 (setq gr
(simplifya (cadr x
) nil
))
2135 (let (($%enumer $numer
))
2136 ;; Switch $%enumer on, when $numer is TRUE to allow
2137 ;; simplification of $%e to its numerical value.
2138 (simplifya (if $ratsimpexpons
($ratsimp
(caddr x
)) (caddr x
))
2142 (setq pot
(ratdisrep pot
))
2145 (cond ((member 'trunc
(car gr
))
2146 (return (srf (list '(mexpt) gr pot
))))
2148 (let ((varlist (caddar gr
)) (genvar (cadddr (car gr
))))
2149 (return (ratrep* (list '(mexpt) gr pot
)))))
2151 (setq gr
(ratdisrep gr
))
2153 ((or (setq mlpgr
(mxorlistp gr
))
2154 (setq mlppot
(mxorlistp pot
)))
2156 ((onep1 pot
) (go atgr
))
2157 ((or (zerop1 pot
) (onep1 gr
)) (go retno
))
2159 ;; This code tries to handle 0^a more complete.
2160 ;; If the sign of realpart(a) is not known return an unsimplified
2161 ;; expression. The handling of the flag *zexptsimp? is not changed.
2162 ;; Reverting the return of an unsimplified 0^a, because timesin
2163 ;; can not handle such expressions. (DK 02/2010)
2165 (cond ((or (member (setq z
($csign pot
)) '($neg $nz
))
2166 (and *zexptsimp?
(eq ($asksign pot
) '$neg
)))
2167 ;; A negative exponent. Maxima error.
2168 (cond ((not errorsw
) (merror (intl:gettext
"expt: undefined: 0 to a negative exponent.")))
2169 (t (throw 'errorsw t
))))
2170 ((and (member z
'($complex $imaginary
))
2171 ;; A complex exponent. Look at the sign of the realpart.
2172 (member (setq z
($sign
($realpart pot
)))
2174 (cond ((not errorsw
)
2175 (merror (intl:gettext
"expt: undefined: 0 to a complex exponent.")))
2176 (t (throw 'errorsw t
))))
2177 ((and *zexptsimp?
(eq ($asksign pot
) '$zero
))
2178 (cond ((not errorsw
)
2179 (merror (intl:gettext
"expt: undefined: 0^0")))
2180 (t (throw 'errorsw t
))))
2181 ((not (member z
'($pos $pz
)))
2182 ;; The sign of realpart(pot) is not known. We can not return
2183 ;; an unsimplified 0^a expression, because timesin can not
2184 ;; handle it. We return ZERO. That is the old behavior.
2185 ;; Look for the imaginary symbol to be consistent with
2187 (cond ((not (free pot
'$%i
))
2188 (cond ((not errorsw
)
2189 (merror (intl:gettext
"expt: undefined: 0 to a complex exponent.")))
2190 (t (throw 'errorsw t
))))
2192 ;; Return ZERO and not an unsimplified expression.
2193 (return (zerores gr pot
)))))
2194 (t (return (zerores gr pot
)))))
2198 (or (not (ratnump gr
)) (not (ratnump pot
))))
2199 (return (eqtest (exptrl gr pot
) check
)))
2200 ;; Check for numerical evaluation of the sqrt.
2201 ((and (alike1 pot
'((rat) 1 2))
2202 (or (setq res
(flonum-eval '%sqrt gr
))
2203 (and (not (member 'simp
(car x
)))
2204 (setq res
(big-float-eval '%sqrt gr
)))))
2207 (return (%itopot pot
)))
2208 ((and (realp gr
) (minusp gr
) (mevenp pot
))
2211 ((and (realp gr
) (minusp gr
) (moddp pot
))
2212 (return (mul2 -
1 (power (- gr
) pot
))))
2213 ((and (equal gr -
1) (maxima-integerp pot
) (mminusp pot
))
2214 (setq pot
(neg pot
))
2217 (maxima-integerp pot
)
2220 (integerp (cadr pot
))
2222 (maxima-integerp (caddr pot
)))
2223 (setq pot
(caddr pot
))
2225 ((atom gr
) (go atgr
))
2226 ((and (eq (caar gr
) 'mabs
)
2229 (or (and (eq $domain
'$real
) (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
))))
2230 (and (eq $domain
'$complex
) (apparently-real-to-judge-by-$csign-p
(cadr gr
)))))
2231 (return (power (cadr gr
) pot
)))
2232 ((and (eq (caar gr
) 'mabs
)
2235 (not (equal pot -
1))
2236 (or (and (eq $domain
'$real
) (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
))))
2237 (and (eq $domain
'$complex
) (apparently-real-to-judge-by-$csign-p
(cadr gr
)))))
2238 ;; abs(x)^(2*n+1) -> abs(x)*x^(2*n), n an integer number
2240 (return (mul (power (cadr gr
) (add pot -
1))
2242 (return (mul (power (cadr gr
) (add pot
1))
2244 ((eq (caar gr
) 'mequal
)
2245 (return (eqtest (list (ncons (caar gr
))
2246 (power (cadr gr
) pot
)
2247 (power (caddr gr
) pot
))
2249 ((symbolp pot
) (go opp
))
2250 ((eq (caar gr
) 'mexpt
) (go e1
))
2251 ((and (eq (caar gr
) '%sum
)
2256 (return (do ((i (1- pot
) (1- i
))
2257 (an gr
(simptimes (list '(mtimes) an gr
) 1 t
)))
2260 (return (eqtest (testt (tms gr pot nil
)) check
)))
2262 (return (eqtest (cond ((and (mplusp gr
)
2263 (not (or (> pot $expop
)
2264 (> (- pot
) $expon
))))
2265 (expandexpt gr pot
))
2266 (t (simplifya (tms gr pot nil
) t
)))
2270 (cond ((eq (caar gr
) 'mexpt
) (go e1
))
2271 ((eq (caar gr
) 'rat
)
2272 (return (mul2 (power (cadr gr
) pot
)
2273 (power (caddr gr
) (mul2 -
1 pot
)))))
2274 ((not (eq (caar gr
) 'mtimes
)) (go up
))
2275 ((or (eq $radexpand
'$all
) (and $radexpand
(simplexpon pot
)))
2278 ((and (or (not (numberp (cadr gr
)))
2279 (equal (cadr gr
) -
1))
2280 (equal -
1 ($num gr
)) ; only for -1
2281 ;; Do not simplify for a complex base.
2282 (not (member ($csign gr
) '($complex $imaginary
)))
2283 (and (eq $domain
'$real
) $radexpand
))
2284 ;; (-1/x)^a -> 1/(-x)^a for x negative
2285 ;; For all other cases (-1)^a/x^a
2286 (if (eq ($csign
(setq w
($denom gr
))) '$neg
)
2287 (return (inv (power (neg w
) pot
)))
2288 (return (div (power -
1 pot
)
2290 ((or (eq $domain
'$complex
) (not $radexpand
)) (go up
)))
2291 (return (do ((l (cdr gr
) (cdr l
)) (res (ncons 1)) (rad))
2293 (cond ((equal res
'(1))
2294 (eqtest (list '(mexpt) gr pot
) check
))
2296 (testt (cons '(mtimes simp
) res
)))
2298 (setq rad
(power* ; RADEXPAND=()?
2299 (cons '(mtimes) (nreverse rad
)) pot
))
2300 (cond ((not (onep1 rad
))
2302 (testt (tms rad
1 (cons '(mtimes) res
))))
2304 (setq *rulesw
* nil res
(cdr rad
))))))
2305 (eqtest (testt (cons '(mtimes) res
)) check
))))
2306 ;; Check with $csign to be more complete. This prevents wrong
2307 ;; simplifications like sqrt(-z^2)->%i*sqrt(z^2) for z complex.
2308 (setq z
($csign
(car l
)))
2309 (if (member z
'($complex $imaginary
))
2310 (setq z
'$pnz
)) ; if appears complex, unknown sign
2311 (setq w
(cond ((member z
'($neg $nz
))
2312 (setq rad
(cons -
1 rad
))
2316 ((alike1 w gr
) (return (list '(mexpt simp
) gr pot
)))
2317 ((member z
'($pn $pnz
))
2318 (setq rad
(cons w rad
)))
2320 (setq w
(testt (tms (simplifya (list '(mexpt) w pot
) t
)
2321 1 (cons '(mtimes) res
))))))
2322 (cond (*rulesw
* (setq *rulesw
* nil res
(cdr w
))))))
2325 (cond ((and (cdr res
) (onep1 (car res
)) (ratnump (cadr res
)))
2326 (setq res
(cdr res
))))
2327 (cond ((null (setq gr
(cdr gr
)))
2328 (return (eqtest (testt (cons '(mtimes) res
)) check
)))
2330 (setq y
(power (cadar gr
) (mult (caddar gr
) pot
))))
2332 (setq y
(%itopot pot
)))
2334 (setq y
(list '(mexpt) (car gr
) pot
)))
2335 (t (setq y
(list '(mexpt simp
) (car gr
) pot
))))
2336 (setq w
(testt (tms (simplifya y t
) 1 (cons '(mtimes) res
))))
2337 (cond (*rulesw
* (setq *rulesw
* nil res
(cdr w
))))
2341 (return (exptrl gr pot
))
2344 (cond ((zerop1 pot
) (go retno
))
2346 (let ((y (mget gr
'$numer
)))
2347 (if (and y
(floatp y
) (or $numer
(not (equal pot
1))))
2348 ;; A numeric constant like %e, %pi, ... and
2349 ;; exponent is a float or bigfloat value.
2350 (return (if (and (member gr
*builtin-numeric-constants
*)
2351 (equal pot bigfloatone
))
2352 ;; Return a bigfloat value.
2354 ;; Return a float value.
2356 ;; In all other cases exptrl simplifies accordingly.
2357 (return (exptrl gr pot
)))))
2359 ;; Numerically evaluate if the power is a flonum.
2361 (let ((val (flonum-eval '%exp pot
)))
2362 (if (float-inf-p val
)
2363 ;; needed for gcl and sbcl - (sometimes) no trap of overflow
2364 (error 'floating-point-overflow
))
2367 ;; Numerically evaluate if the power is a (complex)
2368 ;; big-float. (This is basically the guts of
2369 ;; big-float-eval, but we can't use big-float-eval.)
2370 (when (and (not (member 'simp
(car x
)))
2371 (complex-number-p pot
'bigfloat-or-number-p
))
2372 (let ((x ($realpart pot
))
2373 (y ($imagpart pot
)))
2374 (cond ((and ($bfloatp x
) (like 0 y
))
2375 (return ($bfloat
`((mexpt simp
) $%e
,pot
))))
2376 ((or ($bfloatp x
) ($bfloatp y
))
2377 (let ((z (add ($bfloat x
) (mul '$%i
($bfloat y
)))))
2378 (setq z
($rectform
`((mexpt simp
) $%e
,z
)))
2379 (return ($bfloat z
))))))))
2380 (cond ((and $logsimp
(among '%log pot
)) (return (%etolog pot
)))
2381 ((and $demoivre
(setq z
(demoivre pot
))) (return z
))
2385 ;; Exponent contains %i and %pi and %emode is TRUE:
2386 ;; Check simplification of exp(%i*%pi*p/q*x)
2387 (setq z
(%especial pot
)))
2389 (($taylorp
(third x
))
2390 ;; taylorize %e^taylor(...)
2391 (return ($taylor x
)))))
2393 (let ((y (mget gr
'$numer
)))
2394 ;; Check for a numeric constant.
2398 ;; The exponent is a bigfloat. Convert base to bigfloat.
2400 (member gr
*builtin-numeric-constants
*)
2401 (setq y
($bfloat gr
)))
2402 (and $numer
(integerp pot
)))
2403 (return (exptrl y pot
))))))
2406 (return (eqtest (list '(mexpt) gr pot
) check
))
2410 (cond ((mxorlistp1 gr
) (return (constmx (addk 1 pot
) gr
)))
2412 ((onep1 pot
) (return gr
))
2413 ((or $doallmxops $doscmxops $domxexpt
)
2414 (cond ((or (and mlpgr
2415 (or (not ($listp gr
)) $listarith
)
2416 (scalar-or-constant-p pot $assumescalar
))
2419 (or (not ($listp pot
)) $listarith
)
2420 (scalar-or-constant-p gr $assumescalar
)))
2421 (return (simplifya (outermap1 'mexpt gr pot
) t
)))
2423 ((and $domxmxops
(member pot
'(-1 -
1.0)))
2424 (return (simplifya (outermap1 'mexpt gr pot
) t
)))
2427 ;; At this point we have an expression: (z^a)^b with gr = z^a and pot = b
2428 (cond ((or (eq $radexpand
'$all
)
2429 ;; b is an integer or an odd rational
2431 (and (eq $domain
'$complex
)
2432 (not (member ($csign
(caddr gr
)) '($complex $imaginary
)))
2433 ;; z >= 0 and a not a complex
2434 (or (member ($csign
(cadr gr
)) '($pos $pz $zero
))
2436 (and (mnump (caddr gr
))
2437 (eq ($sign
(sub 1 (take '(mabs) (caddr gr
))))
2439 (and (eq $domain
'$real
)
2440 (member ($csign
(cadr gr
)) '($pos $pz $zero
)))
2441 ;; (1/z)^a -> 1/z^a when z a constant complex
2442 (and (eql (caddr gr
) -
1)
2444 (eq $domain
'$real
))
2445 (and (eq ($csign
(cadr gr
)) '$complex
)
2446 ($constantp
(cadr gr
)))))
2447 ;; This does (1/z)^a -> 1/z^a. This is in general wrong.
2448 ;; We switch this type of simplification on, when
2449 ;; $ratsimpexpons is T. E.g. radcan sets this flag to T.
2450 ;; radcan hangs for expressions like sqrt(1/(1+x)) without
2451 ;; this simplification.
2453 (equal (caddr gr
) -
1))
2456 (odnump (caddr gr
))))
2457 ;; Simplify (z^a)^b -> z^(a*b)
2458 (setq pot
(mul pot
(caddr gr
))
2460 ((and (eq $domain
'$real
)
2463 (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
)))
2464 (evnump (caddr gr
)))
2465 ;; Simplify (x^a)^b -> abs(x)^(a*b)
2466 (setq pot
(mul pot
(caddr gr
))
2467 gr
(radmabs (cadr gr
))))
2470 (mminusp (caddr gr
)))
2471 ;; Simplify (1/z^a)^b -> 1/(z^a)^b
2473 gr
(power (cadr gr
) (neg (caddr gr
)))))
2477 (defun apparently-complex-to-judge-by-$csign-p
(e)
2478 (let ((s ($csign e
)))
2479 (member s
'($complex $imaginary
))))
2481 (defun apparently-real-to-judge-by-$csign-p
(e)
2482 (let ((s ($csign e
)))
2483 (member s
'($pos $neg $zero $pn $pnz $pz $nz
))))
2485 ;; Basically computes log of m base b. Except if m is not a power
2486 ;; of b, we return nil. m is a positive integer and base an integer
2487 ;; not equal to +/-1.
2488 (defun exponent-of (m base
)
2489 ;; Just compute base^k until base^k >= m. Then check if they're equal.
2490 ;; If so, we have the exponent. Otherwise, give up.
2493 (multiple-value-bind (q r
)
2499 (if (zerop expo
) nil expo
)))
2501 (defun timesin (x y w
) ; Multiply X^W into Y
2502 (prog (fm temp z check u expo
)
2503 (if (mexptp x
) (setq check x
))
2505 ;; Prepare the factor x^w and initialize the work of timesin
2509 (setq temp
(cons '(mexpt) (if check
2510 (list (cadr x
) (mult (caddr x
) w
))
2512 (if (and (not timesinp
) (not (eq x
'$%i
)))
2514 (setq temp
(simplifya temp t
))))))
2515 (setq x
(if (mexptp temp
)
2521 ;; Go through the list of terms in fm and look what is to do.
2522 (cond ((null (cdr fm
))
2523 ;; The list of terms is empty. The loop is finshed.
2525 ((or (and (mnump temp
)
2526 (not (or (integerp temp
)
2528 (and (integerp temp
)
2530 ;; Stop the loop for a float or bigfloat number, or number -1.
2533 (cond ((alike1 (car x
) (cadadr fm
))
2534 (cond ((zerop1 (setq w
(plsk (caddr (cadr fm
)) w
)))
2539 (rplacd fm
(cddr fm
))
2540 (cond ((mnump (setq x
(if (mnump (car x
))
2542 (power (car x
) w
))))
2543 (return (rplaca y
(timesk (car y
) x
))))
2548 x
(if (mexptp x
) (cdr x
) (list x
1)))
2552 ((maxima-constantp (car x
))
2555 (cond ((mtimesp (car x
))
2556 ;; A base which is a mtimes expression. Remove
2557 ;; the factor from the lists of products.
2558 (rplacd fm
(cddr fm
))
2559 ;; Multiply the factors of the base with
2560 ;; the list of all remaining products.
2562 (return (muln (nconc y
(cdar x
)) t
)))
2563 (t (return (rplaca (cdr fm
) (car x
))))))
2566 ;; At this place we have to add code for a rational number
2567 ;; as a factor to the list of products.
2569 (or (ratnump (car x
))
2570 (and (integerp (car x
))
2571 (not (onep (car x
))))))
2572 ;; Multiplying bas^k * num/den
2573 (let ((num (num1 (car x
)))
2574 (den (denom1 (car x
)))
2575 (bas (second (cadr fm
))))
2576 (cond ((and (integerp bas
)
2577 (not (eql 1 (abs bas
)))
2578 (setq expo
(exponent-of (abs num
) bas
)))
2579 ;; We have bas^m*bas^k = bas^(k+m).
2580 (setq temp
(power bas
2581 (add (third (cadr fm
)) expo
)))
2582 ;; Set fm to have 1/denom term.
2583 (setq x
(mul (car y
)
2587 ((and (integerp bas
)
2588 (not (eql 1 (abs bas
)))
2589 (setq expo
(exponent-of den bas
)))
2590 (setq expo
(- expo
))
2591 ;; We have bas^(-m)*bas^k = bas^(k-m).
2592 (setq temp
(power bas
2593 (add (third (cadr fm
)) expo
)))
2594 ;; Set fm to have the numerator term.
2595 (setq x
(mul (car y
)
2598 (exptrl bas
(- expo
)))))))
2600 ;; Next term in list of products.
2603 ;; Add in the bas^(k+m) term or bas^(k-m)
2604 (setf y
(rplaca y
1))
2605 (rplacd fm
(cddr fm
))
2606 (rplacd fm
(cons temp
(cdr fm
)))
2612 ((and (not (atom (car x
)))
2613 (eq (caar (car x
)) 'mabs
)
2615 (integerp (caddr (cadr fm
)))
2616 (< (caddr (cadr fm
)) -
1)
2617 (alike1 (cadr (car x
)) (cadr (cadr fm
)))
2618 (not (member ($csign
(cadr (car x
)))
2619 '($complex imaginary
))))
2620 ;; 1/x^n*abs(x) -> 1/(x^(n-2)*abs(x)), where n an integer
2621 ;; Replace 1/x^n -> 1/x^(n-2)
2622 (setq temp
(power (cadr (cadr fm
))
2623 (add (caddr (cadr fm
)) 2)))
2624 (rplacd fm
(cddr fm
))
2625 (if (not (equal temp
1))
2626 (rplacd fm
(cons temp
(cdr fm
))))
2627 ;; Multiply factor 1/abs(x) into list of products.
2628 (setq x
(list (car x
) -
1))
2629 (setq temp
(power (car x
) (cadr x
)))
2633 ((and (not (atom (car x
)))
2634 (eq (caar (car x
)) 'mabs
)
2636 (integerp (caddr (cadr fm
)))
2637 (> (caddr (cadr fm
)) 1)
2638 (alike1 (cadr (car x
)) (cadr (cadr fm
)))
2639 (not (member ($csign
(cadr (car x
)))
2640 '($complex imaginary
))))
2641 ;; x^n/abs(x) -> x^(n-2)*abs(x), where n an integer.
2642 ;; Replace x^n -> x^(n-2)
2643 (setq temp
(power (cadr (cadr fm
))
2644 (add (caddr (cadr fm
)) -
2)))
2645 (rplacd fm
(cddr fm
))
2646 (if (not (equal temp
1))
2647 (rplacd fm
(cons temp
(cdr fm
))))
2648 ;; Multiply factor abs(x) into list of products.
2649 (setq x
(list (car x
) 1))
2650 (setq temp
(power (car x
) (cadr x
)))
2654 ((and (not (atom (cadr fm
)))
2655 (not (atom (cadr (cadr fm
))))
2656 (eq (caaadr (cadr fm
)) 'mabs
)
2657 (equal (caddr (cadr fm
)) -
1)
2660 (alike1 (cadadr (cadr fm
)) (car x
))
2661 (not (member ($csign
(cadadr (cadr fm
)))
2662 '($complex imaginary
))))
2663 ;; 1/abs(x)*x^n -> x^(n-2)*abs(x), where n an integer.
2664 ;; Replace 1/abs(x) -> abs(x)
2665 (setq temp
(cadr (cadr fm
)))
2666 (rplacd fm
(cddr fm
))
2667 (rplacd fm
(cons temp
(cdr fm
)))
2668 ;; Multiply factor x^(n-2) into list of products.
2669 (setq x
(list (car x
) (add (cadr x
) -
2)))
2670 (setq temp
(power (car x
) (cadr x
)))
2674 ((or (maxima-constantp (car x
))
2675 (maxima-constantp (cadadr fm
)))
2676 (if (great temp
(cadr fm
))
2678 ((great (car x
) (cadadr fm
))
2681 ((alike1 (car x
) (cadr fm
))
2684 ;; When a number goto start and look in the next term.
2688 ((and (not (atom (cadr fm
)))
2689 (eq (caar (cadr fm
)) 'mabs
)
2692 (alike1 (cadr (cadr fm
)) (car x
))
2693 (not (member ($csign
(cadr (cadr fm
)))
2694 '($complex imaginary
))))
2695 ;; abs(x)/x^n -> 1/(x^(n-2)*abs(x)), where n an integer.
2696 ;; Replace abs(x) -> 1/abs(x).
2697 (setq temp
(power (cadr fm
) -
1))
2698 (rplacd fm
(cddr fm
))
2699 (rplacd fm
(cons temp
(cdr fm
)))
2700 ;; Multiply factor x^(-n+2) into list of products.
2701 (setq x
(list (car x
) (add (cadr x
) 2)))
2702 (setq temp
(power (car x
) (cadr x
)))
2706 ((maxima-constantp (car x
))
2707 (when (great temp
(cadr fm
))
2709 ((great (car x
) (cadr fm
))
2713 ;; Multiply a number into the list of products.
2714 (return (rplaca y
(timesk (car y
) temp
))))
2715 ((and (eq (car x
) '$%i
)
2718 ((and (eq (car x
) '$%e
)
2721 (return (rplaca y
(timesk (car y
) (exp (float w
))))))
2723 (not (constant (car x
))))
2725 ;; At this point we will insert a mexpt expression,
2726 ;; but first we look at the car of the list of products and
2727 ;; modify the expression if we found a rational number.
2729 (not (onep1 (car y
)))
2730 (or (integerp (car y
))
2732 ;; Multiplying bas^k * num/den.
2733 (let ((num (num1 (car y
)))
2734 (den (denom1 (car y
)))
2736 (cond ((and (integerp bas
)
2737 (not (eql 1 (abs bas
)))
2738 (setq expo
(exponent-of (abs num
) bas
)))
2739 ;; We have bas^m*bas^k.
2740 (setq temp
(power bas
(add (cadr x
) expo
)))
2741 ;; Set fm to have 1/denom term.
2742 (setq x
(div (div num
(exptrl bas expo
)) den
)))
2743 ((and (integerp bas
)
2744 (not (eql 1 (abs bas
)))
2745 (setq expo
(exponent-of den bas
)))
2746 (setq expo
(- expo
))
2747 ;; We have bas^(-m)*bas^k.
2748 (setq temp
(power bas
(add (cadr x
) expo
)))
2749 ;; Set fm to have the numerator term.
2750 (setq x
(div num
(div den
(exptrl bas
(- expo
))))))
2752 ;; The rational doesn't contain any (simple) powers of
2753 ;; the exponential term. We're done.
2754 (return (cdr (rplacd fm
(cons temp
(cdr fm
)))))))
2755 ;; Add in the a^(m+k) or a^(k-m) term.
2756 (setf y
(rplaca y
1))
2757 (rplacd fm
(cons temp
(cdr fm
)))
2763 ((and (maxima-constantp (car x
))
2764 (do ((l (cdr fm
) (cdr l
)))
2766 (when (and (mexptp (cadr l
))
2767 (alike1 (car x
) (cadadr l
)))
2771 ((or (and (mnump (car x
))
2773 (and (eq (car x
) '$%e
)
2777 (setq u
(%especial w
))))
2779 ((alike (cdr check
) x
)
2782 (exptrl (car x
) w
))))
2784 (return (rplaca y
(timesk (car y
) x
))))
2788 (return (cdr (rplacd fm
(cons x
(cdr fm
))))))
2798 (setq temp
(list '(mexpt) (car x
) w
))
2799 (setq temp
(eqtest temp
(or check
'((foo)))))
2800 (return (cdr (rplacd fm
(cons temp
(cdr fm
)))))))
2802 (return (cdr (rplacd fm
(cons (car x
) (cdr fm
)))))
2807 (cond ((and (eq (car x
) '$%i
) (equal w
1))
2808 (rplacd fm
(cddr fm
))
2809 (return (rplaca y
(timesk -
1 (car y
)))))
2810 ((zerop1 (setq w
(plsk 1 w
)))
2812 ((and (mnump (car x
)) (mnump w
))
2813 (return (rplaca (cdr fm
) (exptrl (car x
) w
))))
2814 ((maxima-constantp (car x
))
2817 (setq z
(list '(mexpt) (car x
) w
))
2818 (cond ((alike1 (setq x
(simplifya z t
)) z
)
2819 (return (rplaca (cdr fm
) x
)))
2821 (rplacd fm
(cddr fm
))
2823 (return (muln (cons x y
) t
))))
2825 (rplacd fm
(cddr fm
))
2826 (setq x
(car x
) check nil
)
2829 (setq z
(tms x
1 (setq temp
(cons '(mtimes) y
))))
2830 (return (cond ((eq z temp
)
2833 (setq *rulesw
* t
) z
)))
2835 (return (rplacd fm
(cddr fm
)))
2837 (if (minusp (setq w
(rem w
4)))
2839 (return (cond ((zerop w
)
2842 (rplaca y
(timesk -
1 (car y
))))
2844 (rplaca y
(timesk -
1 (car y
)))
2845 (rplacd fm
(cons '$%i
(cdr fm
))))
2847 (rplacd fm
(cons '$%i
(cdr fm
))))))))
2849 (defun simpmatrix (x vestigial z
)
2850 (declare (ignore vestigial
))
2851 (if (and (null (cddr x
))
2853 (or (eq $scalarmatrixp
'$all
) (member 'mult
(cdar x
)))
2854 ($listp
(cadr x
)) (cdadr x
) (null (cddadr x
)))
2855 (simplifya (cadadr x
) z
)
2856 (let ((badp (dolist (row (cdr x
)) (if (not ($listp row
)) (return t
))))
2857 (args (simpmap (cdr x
) z
)))
2858 (if (and args
(not badp
)) (matcheck args
))
2859 (cons (if badp
'(%matrix simp
) '($matrix simp
)) args
))))
2861 (defun %itopot
(pot)
2863 (let ((i (boole boole-and pot
3)))
2867 (t (list '(mtimes simp
) -
1 '$%i
))))
2868 (power -
1 (mul2 pot
'((rat simp
) 1 2)))))
2871 (cond ((eq (caar pot
) '%log
) (simplifya (cadr pot
) nil
))
2872 ((and (eq (caar pot
) 'mtimes
)
2873 (or (maxima-integerp (cadr pot
))
2874 (and $%e_to_numlog
($numberp
(cadr pot
))))
2875 (not (atom (caddr pot
))) (eq (caar (caddr pot
)) '%log
)
2877 (power (cadr (caddr pot
)) (cadr pot
)))))
2881 loop
(cond ((null pot
)
2882 (cond (a (setq a
(cons '(mtimes) a
))))
2883 (cond (c (setq c
(list '(mexpt simp
) '$%e
(addn c nil
)))))
2884 (return (cond ((null c
) (simptimes a
1 nil
))
2886 (t (simptimes (append a
(list c
)) 1 nil
)))))
2887 ((and (among '%log
(car pot
)) (setq b
(mnlogp (car pot
))))
2888 (setq a
(cons b a
)))
2889 (t (setq c
(cons (car pot
) c
))))
2890 (setq pot
(cdr pot
))
2893 (defun %etolog
(pot) (cond ((mnlogp pot
))
2894 ((eq (caar pot
) 'mplus
) (mnlog (cdr pot
)))
2895 (t (list '(mexpt simp
) '$%e pot
))))
2897 (defun zerores (r1 r2
)
2898 (cond ((or ($bfloatp r1
) ($bfloatp r2
)) bigfloatzero
)
2899 ((or (floatp r1
) (floatp r2
)) 0.0)
2902 (defmfun $orderlessp
(a b
)
2903 (setq a
($totaldisrep
(specrepcheck a
))
2904 b
($totaldisrep
(specrepcheck b
)))
2905 (and (not (alike1 a b
)) (great b a
) t
))
2907 (defmfun $ordergreatp
(a b
)
2908 (setq a
($totaldisrep
(specrepcheck a
))
2909 b
($totaldisrep
(specrepcheck b
)))
2910 (and (not (alike1 a b
)) (great a b
) t
))
2912 ;; Test function to order a and b by magnitude. If it is not possible to
2913 ;; order a and b by magnitude they are ordered by great. This function
2914 ;; can be used by sort, e.g. sort([3,1,7,x,sin(1),minf],ordermagnitudep)
2915 (defmfun $ordermagnitudep
(a b
)
2917 (setq a
($totaldisrep
(specrepcheck a
))
2918 b
($totaldisrep
(specrepcheck b
)))
2919 (cond ((and (or (constp a
) (member a
'($inf $minf
)))
2920 (or (constp b
) (member b
'($inf $minf
)))
2921 (member (setq sgn
($csign
(sub b a
))) '($pos $neg $zero
)))
2922 (cond ((eq sgn
'$pos
) t
)
2923 ((eq sgn
'$zero
) (and (not (alike1 a b
)) (great b a
)))
2925 ((or (constp a
) (member a
'($inf $minf
))) t
)
2926 ((or (constp b
) (member b
'($inf $minf
))) nil
)
2927 (t (and (not (alike1 a b
)) (great b a
))))))
2929 (defun evnump (n) (or (even n
) (and (ratnump n
) (even (cadr n
)))))
2930 (defun odnump (n) (or (and (integerp n
) (oddp n
))
2931 (and (ratnump n
) (oddp (cadr n
)))))
2933 (defun simplexpon (e)
2934 (or (maxima-integerp e
)
2935 (and (eq $domain
'$real
) (ratnump e
) (oddp (caddr e
)))))
2937 ;; This function is not called in Maxima core or share code
2938 ;; and can be cut out.
2940 (and (free p
'$%i
) (member ($sign p
) '($pos $pz $zero
))))
2943 (if (and limitp
(free e
'$%i
)) (asksign-p-or-n e
))
2944 (simplifya (list '(mabs) e
) t
))
2946 (defun simpmqapply (exp y z
)
2947 (let ((simpfun (and (not (atom (cadr exp
))) (safe-get (caaadr exp
) 'specsimp
))) u
)
2949 (funcall simpfun exp y z
)
2950 (progn (setq u
(simpargs exp z
))
2951 (if (symbolp (cadr u
))
2952 (simplifya (cons (cons (cadr u
) (cdar u
)) (cddr u
)) z
)
2955 ;; TRUE, if the symbol e is declared to be $complex or $imaginary.
2956 (defun decl-complexp (e)
2958 (kindp e
'$complex
)))
2960 ;; TRUE, if the symbol e is declared to be $real, $rational, $irrational
2962 (defun decl-realp (e)
2964 (or (kindp e
'$real
)
2965 (kindp e
'$rational
)
2966 (kindp e
'$irrational
)
2967 (kindp e
'$integer
))))
2969 ;; WARNING: Exercise extreme caution when modifying this function!
2971 ;; Richard Fateman and Stavros Macrakis both say that changing the
2972 ;; actual ordering relations (as opposed to making them faster to
2973 ;; determine) could have very subtle and wide-ranging effects. Also,
2974 ;; the simplifier spends the vast majority of its time here, so be
2975 ;; very careful about changes that may drastically slow down the
2983 (cond ((zerop y
) (floatp x
)) (t (plusp y
))))))
2985 (cond ((constant y
) (alphalessp y x
)) (t (numberp y
))))
2987 (cond ((mget y
'$scalar
) (alphalessp y x
))
2988 (t (maxima-constantp y
))))
2990 (cond ((mget y
'$mainvar
) (alphalessp y x
)) (t t
)))
2991 (t (or (maxima-constantp y
) (mget y
'$scalar
)
2992 (and (not (mget y
'$mainvar
)) (not (null (alphalessp y x
))))))))
2993 (t (not (ordfna y x
)))))
2994 ((atom y
) (ordfna x y
))
2996 (cond ((eq (caar y
) 'rat
)
2997 (> (* (caddr y
) (cadr x
)) (* (caddr x
) (cadr y
))))))
2998 ((eq (caar y
) 'rat
))
2999 ((or (member (caar x
) '(mtimes mplus mexpt %del
))
3000 (member (caar y
) '(mtimes mplus mexpt %del
)))
3002 ((and (eq (caar x
) 'bigfloat
) (eq (caar y
) 'bigfloat
)) (mgrp x y
))
3003 ((or (eq (caar x
) 'mrat
) (eq (caar y
) 'mrat
))
3004 (error "GREAT: internal error: unexpected MRAT argument"))
3005 (t (do ((x1 (margs x
) (cdr x1
)) (y1 (margs y
) (cdr y1
))) (())
3007 (return (cond (y1 nil
)
3008 ((not (alike1 (mop x
) (mop y
)))
3009 (great (mop x
) (mop y
)))
3010 ((member 'array
(cdar x
)) t
))))
3011 ((null y1
) (return t
))
3012 ((not (alike1 (car x1
) (car y1
)))
3013 (return (great (car x1
) (car y1
)))))))))
3015 ;; Trivial function used only in ALIKE1.
3016 ;; Should be defined as an open-codable subr.
3018 (defmacro memqarr
(l)
3019 `(if (member 'array
,l
) t
))
3021 ;; Compares two Macsyma expressions ignoring SIMP flags and all other
3022 ;; items in the header except for the ARRAY flag.
3029 (and (arrayp y
) (lisp-array-alike1 x y
)))
3031 ;; NOT SURE IF WE WANT TO ENABLE COMPARISON OF MAXIMA ARRAYS
3032 ;; ASIDE FROM THAT, ADD2LNC CALLS ALIKE1 (VIA MEMALIKE) AND THAT CAUSES TROUBLE
3033 ;; ((maxima-declared-arrayp x)
3034 ;; (and (maxima-declared-arrayp y) (maxima-declared-array-alike1 x y)))
3035 ;; ((maxima-undeclared-arrayp x)
3036 ;; (and (maxima-undeclared-arrayp y) (maxima-undeclared-array-alike1 x y)))
3041 (not (atom (car x
)))
3042 (not (atom (car y
)))
3043 (eq (caar x
) (caar y
)))
3046 ;; Punt back to LIKE, which handles specreps
3049 (eq (memqarr (cdar x
)) (memqarr (cdar y
)))
3050 (alike (cdr x
) (cdr y
))))))))
3052 (defun lisp-array-alike1 (x y
)
3054 (equal (array-dimensions x
) (array-dimensions y
))
3056 (dotimes (i (array-total-size x
))
3057 (if (not (alike1 (row-major-aref x i
) (row-major-aref y i
)))
3058 (return-from lisp-array-alike1 nil
)))
3061 (defun maxima-declared-array-alike1 (x y
)
3062 (lisp-array-alike1 (get (mget x
'array
) 'array
) (get (mget y
'array
) 'array
)))
3064 (defun maxima-undeclared-array-alike1 (x y
)
3066 (alike1 (mfuncall '$arrayinfo x
) (mfuncall '$arrayinfo y
))
3067 (alike1 ($listarray x
) ($listarray y
))))
3069 ;; Maps ALIKE1 down two lists.
3072 (do ((x x
(cdr x
)) (y y
(cdr y
))) ((atom x
) (equal x y
))
3073 (cond ((or (atom y
) (not (alike1 (car x
) (car y
))))
3076 (defun ordfna (e a
) ; A is an atom
3078 (or (not (eq (caar e
) 'rat
))
3079 (> (cadr e
) (* (caddr e
) a
))))
3081 (not (member (caar e
) '(mplus mtimes mexpt
))))
3082 (not (member (caar e
) '(rat bigfloat
))))
3083 ((eq (caar e
) 'mrat
)) ;; all MRATs succeed all atoms
3084 ((null (margs e
)) nil
)
3085 ((eq (caar e
) 'mexpt
)
3086 (cond ((and (maxima-constantp (cadr e
))
3087 (or (not (constant a
)) (not (maxima-constantp (caddr e
)))))
3088 (or (not (free (caddr e
) a
)) (great (caddr e
) a
)))
3089 ((eq (cadr e
) a
) (great (caddr e
) 1))
3090 (t (great (cadr e
) a
))))
3091 ((member (caar e
) '(mplus mtimes
))
3092 (let ((u (car (last e
))))
3093 (cond ((eq u a
) (not (ordhack e
))) (t (great u a
)))))
3094 ((eq (caar e
) '%del
))
3095 ((prog2 (setq e
(car (margs e
))) ; use first arg of e
3096 (and (not (atom e
)) (member (caar e
) '(mplus mtimes
))))
3097 (let ((u (car (last e
)))) ; and compare using
3098 (cond ((eq u a
) (not (ordhack e
))) ; same procedure as above
3103 ;; compare lists a and b elementwise from back to front
3104 (defun ordlist (a b cx cy
)
3106 (setq l1
(length a
) l2
(length b
))
3107 loop
(cond ((= l1
0)
3108 (return (cond ((= l2
0) (eq cx
'mplus
))
3109 ((and (eq cx cy
) (= l2
1))
3110 (great (cond ((eq cx
'mplus
) 0) (t 1)) (car b
))))))
3111 ((= l2
0) (return (not (ordlist b a cy cx
)))))
3112 (setq c
(nthelem l1 a
) d
(nthelem l2 b
))
3113 (cond ((not (alike1 c d
)) (return (great c d
))))
3114 (setq l1
(1- l1
) l2
(1- l2
))
3117 (defun term-list (x)
3122 (defun factor-list (x)
3127 ;; one of the exprs x or y should be one of:
3128 ;; %del, mexpt, mplus, mtimes
3130 (let ((cx (caar x
)) (cy (caar y
)))
3131 (cond ((eq cx
'%del
) (if (eq cy
'%del
) (great (cadr x
) (cadr y
)) t
))
3133 ((or (eq cx
'mtimes
) (eq cy
'mtimes
))
3134 (ordlist (factor-list x
) (factor-list y
) 'mtimes
'mtimes
))
3135 ((or (eq cx
'mplus
) (eq cy
'mplus
))
3136 (ordlist (term-list x
) (term-list y
) 'mplus
'mplus
))
3137 ((eq cx
'mexpt
) (ordmexpt x y
))
3138 ((eq cy
'mexpt
) (not (ordmexpt y x
))))))
3141 (if (and (cddr x
) (null (cdddr x
)))
3142 (great (if (eq (caar x
) 'mplus
) 0 1) (cadr x
))))
3144 (defun ordmexpt (x y
)
3145 (cond ((eq (caar y
) 'mexpt
)
3146 (cond ((alike1 (cadr x
) (cadr y
)) (great (caddr x
) (caddr y
)))
3147 ((maxima-constantp (cadr x
))
3148 (if (maxima-constantp (cadr y
))
3149 (if (or (alike1 (caddr x
) (caddr y
))
3150 (and (mnump (caddr x
)) (mnump (caddr y
))))
3151 (great (cadr x
) (cadr y
))
3152 (great (caddr x
) (caddr y
)))
3153 (great x
(cadr y
))))
3154 ((maxima-constantp (cadr y
)) (great (cadr x
) y
))
3156 (great (cadr x
) (if (mnump (caddr y
)) (cadr y
) y
)))
3157 ((mnump (caddr y
)) (great x
(cadr y
)))
3158 (t (let ((x1 (simpln1 x
)) (y1 (simpln1 y
)))
3159 (if (alike1 x1 y1
) (great (cadr x
) (cadr y
))
3161 ((alike1 (cadr x
) y
) (great (caddr x
) 1))
3162 ((mnump (caddr x
)) (great (cadr x
) y
))
3163 (t (great (simpln1 x
) (simpln (list '(%log
) y
) 1 t
)))))
3165 (defmfun $multthru
(e1 &optional e2
)
3167 (cond (e2 ;called with two args
3168 (setq arg1
(specrepcheck e1
)
3169 arg2
(specrepcheck e2
))
3170 (cond ((or (atom arg2
)
3171 (not (member (caar arg2
) '(mplus mequal
))))
3173 ((eq (caar arg2
) 'mequal
)
3174 (list (car arg2
) ($multthru arg1
(cadr arg2
))
3175 ($multthru arg1
(caddr arg2
))))
3176 (t (expandterms arg1
(cdr arg2
)))))
3177 (t ;called with only one arg
3179 (setq arg1
(setq arg2
(specrepcheck e1
)))
3180 (cond ((atom arg1
) (return arg1
))
3181 ((eq (caar arg1
) 'mnctimes
)
3182 (setq arg1
(cdr arg1
)) (go nct
))
3183 ((not (eq (caar arg1
) 'mtimes
)) (return arg1
)))
3184 (setq arg1
(reverse (cdr arg1
)))
3185 times
(when (mplusp (car arg1
))
3186 (setq l1
(nconc l1
(cdr arg1
)))
3187 (return (expandterms (muln l1 t
) (cdar arg1
))))
3188 (setq l1
(cons (car arg1
) l1
))
3189 (setq arg1
(cdr arg1
))
3190 (if (null arg1
) (return arg2
))
3192 nct
(when (mplusp (car arg1
))
3193 (setq l1
(nreverse l1
))
3194 (return (addn (mapcar
3198 (append l1
(ncons u
) (cdr arg1
)))
3202 (setq l1
(cons (car arg1
) l1
))
3203 (setq arg1
(cdr arg1
))
3204 (if (null arg1
) (return arg2
))
3207 ;; EXPANDEXPT computes the expansion of (x1 + x2 + ... + xm)^n
3208 ;; taking a sum and integer power as arguments.
3209 ;; Its theory is to recurse down the binomial expansion of
3210 ;; (x1 + (x2 + x3 + ... + xm))^n using the Binomial Expansion
3211 ;; Thus it does a sigma:
3215 ;; \ / n \ k (n - k)
3216 ;; > | | x1 (x2 + x3 + ... + xm)
3221 ;; The function EXPONENTIATE-SUM does this and recurses through the second
3222 ;; sum raised to a power. It takes a list of terms and a positive integer
3223 ;; power as arguments.
3226 (defun expandexpt (sum power
)
3227 (declare (fixnum power
))
3228 (let ((expansion (exponentiate-sum (cdr sum
) (abs power
))))
3229 (cond ((plusp power
) expansion
)
3230 (t (inv expansion
)))))
3232 (defun exponentiate-sum (terms rpower
)
3233 (declare (fixnum rpower
))
3234 (cond ((= rpower
0) 1)
3235 ((null (cdr terms
)) (power (car terms
) rpower
))
3236 ((= rpower
1) (cons '(mplus simp
) terms
))
3237 (t (do ((i 0 (1+ i
))
3238 (result 0 (add2 result
3239 (muln (list (combination rpower i
)
3240 (exponentiate-sum (cdr terms
)
3242 (power (car terms
) i
)) t
))))
3243 ((> i rpower
) result
)
3244 (declare (fixnum i
))))))
3246 ;; Computes the combination of n elements taken m at a time by the formula
3248 ;; (n * (n-1) * ... * (n - m + 1)) / m! =
3249 ;; (n / 1) * ((n - 1) / 2) * ... * ((n - m + 1) / m)
3251 ;; Checks for the case when m is greater than n/2 and translates
3252 ;; to an equivalent expression.
3254 (defun combination (n m
)
3255 (declare (fixnum n m
))
3256 (cond ((> m
(truncate n
2))
3257 (combination n
(- n m
)))
3259 (do ((result 1 (truncate (* result n1
) m1
))
3263 (declare (fixnum n1 m1
))))))
3265 (defun expandsums (a b
)
3267 (setq a
(fixexpand a
) b
(cdr b
))
3269 (when (null a
) (return c
))
3270 (setq c
(cons (expandterms (car a
) b
) c
))
3275 (defun expandterms (a b
)
3278 (when (null b
) (return c
))
3279 (setq c
(cons (mul2 a
(car b
)) c
))
3284 (defun genexpands (l)
3289 (setq *prods
* (nreverse *prods
*)
3290 *negprods
* (nreverse *negprods
*)
3291 *sums
* (nreverse *sums
*)
3292 *negsums
* (nreverse *negsums
*))
3295 (push (car l
) *prods
*))
3296 ((eq (caaar l
) 'rat
)
3297 (unless (equal (cadar l
) 1)
3298 (push (cadar l
) *prods
*))
3299 (push (caddar l
) *negprods
*))
3300 ((eq (caaar l
) 'mplus
)
3301 (push (car l
) *sums
*))
3302 ((and (eq (caaar l
) 'mexpt
)
3303 (equal (caddar l
) -
1) (mplusp (cadar l
)))
3304 (push (cadar l
) *negsums
*))
3305 ((and (eq (caaar l
) 'mexpt
)
3307 (mminusp (caddar l
))))
3308 (push (if (equal (caddar l
) -
1)
3310 (list (caar l
) (cadar l
) (neg (caddar l
))))
3313 (push (car l
) *prods
*)))
3316 (defun expandtimes (a)
3317 (prog (*prods
* *negprods
* *sums
* *negsums
* expsums expnegsums
)
3319 (setq *prods
* (cond ((null *prods
*) 1)
3320 ((null (cdr *prods
*)) (car *prods
*))
3321 (t (cons '(mtimes simp
) *prods
*))))
3322 (setq *negprods
* (cond ((null *negprods
*) 1)
3323 ((null (cdr *negprods
*)) (car *negprods
*))
3324 (t (cons '(mtimes simp
) *negprods
*))))
3325 (cond ((null *sums
*) (go down
))
3326 (t (setq expsums
(car *sums
*))
3328 (setq expsums
(expandsums expsums c
)))
3330 (setq *prods
* (cond ((equal *prods
* 1) expsums
)
3331 (t (expandterms *prods
* (fixexpand expsums
)))))
3332 down
(cond ((null *negsums
*)
3333 (cond ((equal 1 *negprods
*) (return *prods
*))
3335 (return (expandterms (power *negprods
* -
1) (cdr *prods
*))))
3336 (t (return (let ((expandflag t
))
3337 (mul2 *prods
* (power *negprods
* -
1)))))))
3339 (setq expnegsums
(car *negsums
*))
3341 (setq expnegsums
(expandsums expnegsums c
)))
3343 (setq expnegsums
(expandterms *negprods
* (fixexpand expnegsums
)))
3344 (return (if (mplusp *prods
*)
3345 (expandterms (inv expnegsums
) (cdr *prods
*))
3346 (let ((expandflag t
))
3347 (mul2 *prods
* (inv expnegsums
)))))))
3349 (defun expand1 (exp $expop $expon
)
3350 (unless (and (integerp $expop
) (> $expop -
1))
3351 (merror (intl:gettext
"expand: expop must be a nonnegative integer; found: ~M") $expop
))
3352 (unless (and (integerp $expon
) (> $expon -
1))
3353 (merror (intl:gettext
"expand: expon must be a nonnegative integer; found: ~M") $expon
))
3354 (resimplify (specrepcheck exp
)))
3356 (defmfun $expand
(exp &optional
(expop $maxposex
) (expon $maxnegex
))
3357 (expand1 exp expop expon
))
3359 (defun fixexpand (a)
3360 (if (not (mplusp a
))
3364 (defun simpnrt (x *n
) ; computes X^(1/*N)
3365 (prog (*in
*out varlist genvar $factorflag $dontfactor
)
3366 (setq $factorflag t
)
3368 (setq x
(ratrep* x
))
3369 (when (equal (cadr x
) 0) (return 0))
3370 (setq x
(ratfact (cdr x
) 'psqfr
))
3371 (simpnrt1 (mapcar #'pdis x
))
3372 (setq *out
(if *out
(muln *out nil
) 1))
3373 (setq *in
(cond (*in
3374 (setq *in
(muln *in nil
))
3377 (return (let (($%emode t
))
3378 (simplifya (list '(mtimes) *in
*out
)
3381 (member (caaadr *in
) '(mplus mtimes rat
)))))))))
3384 (do ((x x
(cddr x
)) (y))
3386 (cond ((not (equal 1 (setq y
(gcd (cadr x
) *n
))))
3387 (push (simpnrt (list '(mexpt) (car x
) (quotient (cadr x
) y
))
3390 ((and (equal (cadr x
) 1) (integerp (car x
)) (plusp (car x
))
3391 (setq y
(pnthrootp (car x
) *n
)))
3394 (unless (> *n
(abs (cadr x
)))
3395 (push (list '(mexpt) (car x
) (quotient (cadr x
) *n
)) *out
))
3396 (push (list '(mexpt) (car x
) (rem (cadr x
) *n
)) *in
)))))
3398 (defun nrthk (in *n
)
3404 ((eq $domain
'$real
)
3410 (power* '$%e
(list '(mtimes) (list '(rat) 1 *n
) '$%pi
'$%i
))))
3413 ((or (and wflag
(eq ($asksign in
) '$neg
))
3414 (and (mnump in
) (equal ($sign in
) '$neg
)))
3415 (nrthk1 (mul2* -
1 in
) *n
))
3419 (defun nrthk1 (in *n
) ; computes (-IN)^(1/*N)
3421 (mul2 (nrthk2 in
*n
) (nrthk -
1 *n
))
3422 (nrthk2 (mul2* -
1 in
) *n
)))
3424 (defun nrthk2 (in *n
)
3425 (power* in
(list '(rat) 1 *n
))) ; computes IN^(1/*N)
3427 ;; The following was formerly in SININT. This code was placed here because
3428 ;; SININT is now an out-of-core file on MC, and this code is needed in-core
3429 ;; because of the various calls to it. - BMT & JPG
3431 (declare-top (special var ratform context
))
3433 (defmfun $integrate
(expr x
&optional lo hi
)
3434 (declare (special *in-risch-p
*))
3437 (with-new-context (context)
3438 (if (member '%risch
*nounl
*)
3440 ;; Give up; we're being called from RISCHINT by some path.
3441 (list '(%integrate
) expr x
)
3444 ($defint expr x lo hi
))))
3448 ((member (caar a
) '(mplus mtimes
))
3449 (do ((l (cdr a
) (cdr l
))) ((null l
) t
)
3450 (or (ratp (car l
) var
) (return nil
))))
3451 ((eq (caar a
) 'mexpt
)
3452 (if (free (cadr a
) var
)
3453 (free (caddr a
) var
)
3454 (and (integerp (caddr a
)) (ratp (cadr a
) var
))))
3457 (defun ratnumerator (r)
3459 ((atom (cdr r
)) (car r
))
3460 ((numberp (cadr r
)) r
)
3463 (defun ratdenominator (r)
3465 ((atom (cdr r
)) (cdr r
))
3466 ((numberp (cadr r
)) 1)
3469 (declare-top (special var
))
3471 ;; (BPROG U V) appears to return A and B (as ((A1 . A2) B1 . B2) with A = A1/A2, B = B1/B2)
3472 ;; such that B/U + A/V = 1/(U*V), where U, V are polynomials represented as a list of
3473 ;; exponents and coefficients, (<gensym> E1 C1 E2 C2 ...) = C1*<gensym>^E1 + C2*<gensym>^E2 + ....
3475 ;; (%i73) partfrac ((2*x^2-3)/(x^4-3*x^2+2), x);
3476 ;; 1. Trace: (PARTFRAC '((#:X16910 2 2 0 -3) #:X16910 4 1 2 -3 0 2) '#:X16910)
3477 ;; 2. Trace: (BPROG '(#:X16910 2 1 0 -2) '(#:X16910 2 1 0 -1))
3478 ;; 2. Trace: BPROG ==> ((-1 . 1) 1 . 1)
3479 ;; 2. Trace: (BPROG '(#:X16910 1 1 0 1) '(#:X16910 1 1 0 -1))
3480 ;; 2. Trace: BPROG ==> ((1 . 2) -1 . 2)
3481 ;; 2. Trace: (BPROG '(#:X16910 1 1 0 -1) '1)
3482 ;; 2. Trace: BPROG ==> ((0 . 1) 1 . 1)
3483 ;; 1. Trace: PARTFRAC ==> ((0 . 1) ((1 . 2) (#:X16910 1 1 0 -1) 1) ((-1 . 2) (#:X16910 1 1 0 1) 1) ((1 . 1) (#:X16910 2 1 0 -2) 1))
3484 ;; (%o73) 1/(x^2-2)-1/(2*(x+1))+1/(2*(x-1))
3487 (prog (p1b p2b coef1r coef2r coef1s coef2s f1 f2 a egcd oldalg
)
3488 (setq oldalg $algebraic
)
3491 (setq coef2r
(setq coef1s
0))
3492 (setq coef2s
(setq coef1r
1))
3495 (unless (zerop (pdegree p1b var
)) (setq egcd
(pgcdexpon p1b
)))
3497 (unless (or (zerop (pdegree p2b var
)) (= egcd
1))
3498 (setq egcd
(gcd egcd
(pgcdexpon p2b
)))
3499 (setq p1b
(pexpon*// p1b egcd nil
)
3500 p2b
(pexpon*// p2b egcd nil
)))
3501 b1
(cond ((< (pdegree p1b var
) (pdegree p2b var
))
3503 (rotatef coef1r coef2r
)
3504 (rotatef coef1s coef2s
)))
3505 (when (zerop (pdegree p2b var
))
3506 (unless (zerop (pdegree coef2r var
))
3507 (setq coef2r
(pexpon*// coef2r egcd t
)))
3508 (unless (zerop (pdegree coef2s var
))
3509 (setq coef2s
(pexpon*// coef2s egcd t
)))
3510 (return (cons (ratreduce (ptimes (cdr r
) coef2r
) p2b
)
3511 (ratreduce (ptimes (cdr s
) coef2s
) p2b
))))
3512 (setq f1
(psquorem1 (cdr p1b
) (cdr p2b
) t
))
3513 (setq $algebraic $false
)
3514 (setq f2
(psimp var
(cadr f1
)))
3515 (setq p1b
(pquotientchk (psimp var
(caddr f1
)) a
))
3516 (setq $algebraic oldalg
)
3518 (setq coef1r
(pquotientchk (pdifference (ptimes f1 coef1r
)
3521 (setq coef1s
(pquotientchk (pdifference (ptimes f1 coef1s
)
3527 (defun ratdifference (a b
) (ratplus a
(ratminus b
)))
3529 (defun ratpl (a b
) (ratplus (ratfix a
) (ratfix b
)))
3531 (defun ratti (a b c
) (rattimes (ratfix a
) (ratfix b
) c
))
3533 (defun ratqu (a b
) (ratquotient (ratfix a
) (ratfix b
)))
3535 (defun ratfix (a) (cond ((equal a
(ratnumerator a
)) (cons a
1)) (t a
)))
3537 (defun ratdivide (f g
)
3538 (destructuring-let* (((fnum . fden
) (ratfix f
))
3539 ((gnum . gden
) (ratfix g
))
3540 ((q r
) (pdivide fnum gnum
)))
3541 (cons (ratqu (ratti q gden t
) fden
)
3544 (defun polcoef (l n
) (cond ((or (atom l
) (pointergp var
(car l
)))
3545 (cond ((equal n
0) l
) (t 0)))
3546 (t (ptterm (cdr l
) n
))))
3548 (defun disrep (l) (cond ((equal (ratnumerator l
) l
)
3549 ($ratdisrep
(cons ratform
(cons l
1))))
3550 (t ($ratdisrep
(cons ratform l
)))))
3552 (declare-top (unspecial var
))
3555 ;; The following was formerly in MATRUN. This code was placed here because
3556 ;; MATRUN is now an out-of-core file on MC, and this code is needed in-core
3557 ;; so that MACSYMA SAVE files will work. - JPG
3562 (defun kar (x) (if (atom x
) (matcherr) (car x
)))
3564 (defun kaar (x) (kar (kar x
)))
3566 (defun kdr (x) (if (atom x
) (matcherr) (cdr x
)))
3568 (defun simpargs1 (a vestigial c
)
3569 (declare (ignore vestigial
))
3573 (if (not (atom x
)) (car x
)))
3575 (defquote retlist
(&rest l
)
3577 (mapcar #'(lambda (z) (list '(mequal simp
) z
(meval z
))) l
)))
3580 (if (zerop c
) x
(nthkdr (kdr x
) (1- c
))))
