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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
17 ;; Switches dealing with matrices and non-commutative multiplication.
19 (defmvar $doscmxplus nil
20 "Causes SCALAR + MATRIX to return a matrix answer. This switch
21 is not subsumed under DOALLMXOPS.")
24 "Causes SCALAR^MATRIX([1,2],[3,4]) to return
25 MATRIX([SCALAR,SCALAR^2],[SCALAR^3,SCALAR^4]). In general, this
26 transformation affects exponentiations where the *print-base* is a
27 scalar and the power is a matrix or list.")
29 (defmvar $domxplus nil
)
31 (defmvar $domxtimes nil
)
35 ;; Lisp level variables
36 (defmvar derivsimp t
"Hack in `simpderiv' for `rwg'")
38 (defmvar $grindswitch nil
)
41 (defmvar $limitdomain
'$complex
)
42 (defmvar $lognegint nil
44 :properties
((evflag t
)))
45 (defmvar $distribute_over t
46 "If T, functions are distributed over bags.")
48 (defvar *expandp
* nil
)
49 (defvar *timesinp
* nil
)
50 (defvar *exptrlsw
* nil
)
51 (defvar *expandflag
* nil
)
53 (defprop mnctimes t associative
)
54 (defprop lambda t lisp-no-simp
)
56 ;; Local functions should not be simplified. Various lisps
57 ;; use various names for the list structure defining these:
62 (dolist (y (list x z
))
65 (setf (get (car y
) 'lisp-no-simp
) t
))))))
67 (dolist (x '(mplus mtimes mnctimes mexpt mncexpt %sum
))
68 (setf (get x
'msimpind
) (cons x
'(simp))))
70 ;; operators properties
72 (mapc #'(lambda (x) (setf (get (first x
) 'operators
) (second x
)))
73 '((mplus simplus
) (mtimes simptimes
) (mncexpt simpncexpt
)
76 (mnctimes simpnct
) (mquotient simpquot
) (mexpt simpexpt
)
77 (%derivative simpderiv
)
79 (%integrate simpinteg
) (%limit simp-limit
)
80 (bigfloat simpbigfloat
) (lambda simplambda
) (mdefine simpmdef
)
84 ($matrix simpmatrix
) (%matrix simpmatrix
)))
86 (defprop $li lisimp specsimp
)
87 (defprop $psi psisimp specsimp
)
89 (defprop $equal t binary
)
90 (defprop $notequal t binary
)
93 "Returns true if X is a bigfloat"
96 (eq (caar x
) 'bigfloat
)))
99 "Returns non-NIL if X is an integer, float, or bfloat that is equal
101 (or (and (integerp x
) (= 0 x
))
102 (and (floatp x
) (= 0.0 x
))
103 (and ($bfloatp x
) (= 0 (second x
)))))
106 "Returns non-NIL if X is an integer, float, or bfloat that is equal
108 (or (and (integerp x
) (= 1 x
))
109 (and (floatp x
) (= 1.0 x
))
110 (and ($bfloatp x
) (zerop1 (sub x
1)))))
113 "Returns non-NIL if X is a Lisp number or if it is a Maxima rational
114 form or a bigfloat form"
116 (and (not (atom x
)) (not (atom (car x
)))
117 (member (caar x
) '(rat bigfloat
)))))
119 ;; Does X or a subexpression match PREDICATE?
121 ;; If X is a tree then we recurse depth-first down its arguments. The specrep
122 ;; check is because rat forms are built rather differently from normal Maxima
123 ;; expressions so we need to unpack them for the recursion to work properly.
124 (defun subexpression-matches-p (predicate x
)
125 (or (funcall predicate x
)
128 (subexpression-matches-p predicate
(specdisrep x
))
129 (some (lambda (arg) (subexpression-matches-p predicate arg
))
132 ;; Is there a bfloat anywhere in X?
133 (defun some-bfloatp (x) (subexpression-matches-p '$bfloatp x
))
135 ;; Is there a float anywhere in X?
136 (defun some-floatp (x) (subexpression-matches-p 'floatp x
))
138 ;; EVEN works for any arbitrary lisp object since it does an integer
139 ;; check first. In other cases, you may want the Lisp EVENP function
140 ;; which only works for integers.
142 (defun even (a) (and (integerp a
) (not (oddp a
))))
144 ;; Predicates to determine if X satisfies some condition.
147 "Determines if X is a Maxima rational form: ((rat ...) a b)"
148 (and (not (atom x
)) (eq (caar x
) 'rat
)))
151 "Determines if X is a Maxima sum form: ((mplus ...) ...)"
152 (and (not (atom x
)) (eq (caar x
) 'mplus
)))
155 "Determines if X is a Maxima product form: ((mtimes ...) ...)"
156 (and (not (atom x
)) (eq (caar x
) 'mtimes
)))
159 "Determines if X is a Maxima exponential form: ((mexpt ...) ...)"
160 (and (not (atom x
)) (eq (caar x
) 'mexpt
)))
162 (defun mnctimesp (x) (and (not (atom x
)) (eq (caar x
) 'mnctimes
)))
164 (defun mncexptp (x) (and (not (atom x
)) (eq (caar x
) 'mncexpt
)))
167 "Determines if X is a Maxima log form: ((%log ...) ...)"
168 (and (not (atom x
)) (eq (caar x
) '%log
)))
171 "Determines if X is a Maxima negative form: ((mminus ...) ...)
173 This generally only happens on input forms like a - b:
174 ((mplus) $a ((mminus) $b)).
175 After simplification a - b becomes
176 ((mplus) $a ((mtimes) -1 $b))"
178 (and (not (atom x
)) (eq (caar x
) 'mminus
)))
181 "Determines if X is negative if X is a Lisp number or a Maxima rat
182 form or bigfloat form"
183 (cond ((realp x
) (minusp x
))
184 ((or (ratnump x
) ($bfloatp x
)) (minusp (cadr x
)))))
186 (defun mqapplyp (e) (and (not (atom e
)) (eq (caar e
) 'mqapply
)))
188 (defun ratdisrep (e) (simplifya ($ratdisrep e
) nil
))
190 (defun sratsimp (e) (simplifya ($ratsimp e
) nil
))
192 (defun simpcheck (e flag
)
193 (cond ((specrepp e
) (specdisrep e
))
195 (t (let (($%enumer $numer
))
196 ;; Switch $%enumer on, when $numer is TRUE to allow
197 ;; simplification of $%e to its numerical value.
198 (simplifya e nil
)))))
200 (defun mratcheck (e) (if ($ratp e
) (ratdisrep e
) e
))
202 (defmfun $numberp
(e)
203 "Returns true if E is a Maxima rational, a float, or a bigfloat number"
204 (or ($ratnump e
) ($floatnump e
) ($bfloatp e
)))
206 (defmfun $integerp
(x)
209 (not (member 'trunc
(car x
)))
211 (equal (cddr x
) 1))))
213 ;; The call to $INTEGERP in the following two functions checks for a CRE
214 ;; rational number with an integral numerator and a unity denominator.
217 (cond ((integerp x
) (oddp x
))
218 (($integerp x
) (oddp (cadr x
)))))
221 (cond ((integerp x
) (evenp x
))
222 (($integerp x
) (not (oddp (cadr x
))))))
224 (defmfun $floatnump
(x)
226 (and ($ratp x
) (floatp (cadr x
)) (onep1 (cddr x
)))))
228 (defmfun $ratnump
(x)
232 (not (member 'trunc
(car x
)))
234 (integerp (cddr x
)))))
239 (eq (caar x
) 'mrat
)))
241 (defmfun $taylorp
(x)
244 (member 'trunc
(cdar x
)) t
))
246 (defun specrepcheck (e) (if (specrepp e
) (specdisrep e
) e
))
248 ;; Note that the following two functions are carefully coupled.
252 (member (caar e
) '(mrat mpois
))))
254 (defun specdisrep (e)
255 (cond ((eq (caar e
) 'mrat
) (ratdisrep e
))
258 (defmfun $polysign
(x)
259 (setq x
(cadr (ratf x
)))
260 (cond ((equal x
0) 0) ((pminusp x
) -
1) (t 1)))
262 ;; These check for the correct number of operands within Macsyma expressions,
263 ;; not arguments in a procedure call as the name may imply.
265 (defun arg-count-check (required-arg-count expr
)
266 (unless (= required-arg-count
(length (rest expr
)))
267 (wna-err expr required-arg-count
)))
269 (defun oneargcheck (expr)
270 (arg-count-check 1 expr
))
272 (defun twoargcheck (expr)
273 (arg-count-check 2 expr
))
275 ;; WNA-ERR: Wrong Number of Arguments error
277 ;; If REQUIRED-ARG-COUNT is non-NIL, then we check that EXPR has the
278 ;; correct number of arguments. A informative error message is shown
279 ;; if the number of arguments is not given.
281 ;; Otherwise, EXPR must be a symbol and a generic message is printed.
282 ;; (This is for backward compatibility for existing uses of WNA-ERR.)
283 (defun wna-err (exprs &optional required-arg-count
)
284 (if required-arg-count
285 (let ((op (caar exprs
))
286 (actual-count (length (rest exprs
))))
287 (merror (intl:gettext
"~M: expected exactly ~M arguments but got ~M: ~M")
288 op required-arg-count actual-count
(list* '(mlist) (rest exprs
))))
289 (merror (intl:gettext
"~:@M: wrong number of arguments.")
292 (defun improper-arg-err (exp fn
)
293 (merror (intl:gettext
"~:M: improper argument: ~M") fn exp
))
295 (defun subargcheck (form subsharp argsharp fun
)
296 (if (or (not (= (length (subfunsubs form
)) subsharp
))
297 (not (= (length (subfunargs form
)) argsharp
)))
298 (merror (intl:gettext
"~:@M: wrong number of arguments or subscripts.") fun
)))
300 ;; Constructor and extractor primitives for subscripted functions, e.g.
301 ;; F[1,2](X,Y). SUBL is (1 2) and ARGL is (X Y).
303 ;; These will be flushed when NOPERS is finished. They will be macros in
304 ;; NOPERS instead of functions, so we have to be careful that they aren't
305 ;; mapped or applied anyplace. What we really want is open-codable routines.
307 (defun subfunmakes (fun subl argl
)
308 `((mqapply simp
) ((,fun simp array
) .
,subl
) .
,argl
))
310 (defun subfunmake (fun subl argl
)
311 `((mqapply) ((,fun simp array
) .
,subl
) .
,argl
))
313 (defun subfunname (exp) (caaadr exp
))
315 (defun subfunsubs (exp) (cdadr exp
))
317 (defun subfunargs (exp) (cddr exp
))
319 (defmfun $numfactor
(x)
320 (setq x
(specrepcheck x
))
323 ((not (eq (caar x
) 'mtimes
)) 1)
324 ((mnump (cadr x
)) (cadr x
))
327 (defun scalar-or-constant-p (x flag
)
328 (if flag
(not ($nonscalarp x
)) ($scalarp x
)))
330 (defmfun $constantp
(x)
331 (cond ((atom x
) (or ($numberp x
) (kindp x
'$constant
)))
332 ((member (caar x
) '(rat bigfloat
)) t
)
333 ((specrepp x
) ($constantp
(specdisrep x
)))
334 ((or (mopp (caar x
)) (kindp (caar x
) '$constant
))
335 (do ((x (cdr x
) (cdr x
))) ((null x
) t
)
336 (if (not ($constantp
(car x
))) (return nil
))))))
339 (cond ((symbolp x
) (kindp x
'$constant
))
341 (and (kindp (caar x
) '$constant
)
342 (do ((x (cdr x
) (cdr x
))) ((null x
) t
)
343 (if (not ($constantp
(car x
))) (return nil
)))))))
345 (defun maxima-constantp (x)
347 (and (symbolp x
) (kindp x
'$constant
))))
349 (defun consttermp (x) (and ($constantp x
) (not ($nonscalarp x
))))
351 (defmfun $scalarp
(x) (or (consttermp x
) (eq (scalarclass x
) '$scalar
)))
353 (defmfun $nonscalarp
(x) (eq (scalarclass x
) '$nonscalar
))
355 (defun scalarclass (exp) ; Returns $SCALAR, $NONSCALAR, or NIL (unknown).
357 ;; Maxima numbers are scalar.
360 (cond ((or (mget exp
'$nonscalar
)
361 (and (not (mget exp
'$scalar
))
362 ;; Arrays are nonscalar, but not if declared scalar.
364 ($member exp $arrays
))))
366 ((or (mget exp
'$scalar
)
367 ;; Include constant atoms which are not declared nonscalar.
370 ((member 'array
(car exp
))
371 (cond ((mget (caar exp
) '$scalar
) '$scalar
)
372 ((mget (caar exp
) '$nonscalar
) '$nonscalar
)
374 ((specrepp exp
) (scalarclass (specdisrep exp
)))
375 ;; If the function is declared scalar or nonscalar, then return. If it
376 ;; isn't explicitly declared, then try to be intelligent by looking at
377 ;; the arguments to the function.
378 ((scalarclass (caar exp
)))
379 ;; <number> + <scalar> is SCALARP because that seems to be useful.
380 ;; This should probably only be true if <number> is a member of the
381 ;; field of scalars. <number> * <scalar> is SCALARP since
382 ;; <scalar> + <scalar> is SCALARP. Also, this has to be done to make
383 ;; <scalar> - <scalar> SCALARP.
384 ((member (caar exp
) '(mplus mtimes
))
385 (do ((l (cdr exp
) (cdr l
))) ((null l
) '$scalar
)
386 (if (not (consttermp (car l
)))
387 (return (scalarclass-list l
)))))
388 ((and (eq (caar exp
) 'mqapply
) (scalarclass (cadr exp
))))
389 ((mxorlistp exp
) '$nonscalar
)
390 ;; If we can't find out anything about the operator, then look at the
391 ;; arguments to the operator. I think NIL should be returned at this
394 (do ((exp (cdr exp
) (cdr exp
)) (l '(1)))
395 ((null exp
) (scalarclass-list l
))
396 (if (not (consttermp (car exp
)))
397 (setq l
(cons (car exp
) l
)))))))
399 ;; Could also do <scalar> +|-|*|/ |^ <declared constant>, but this is not
400 ;; always correct and could screw somebody.
402 ;; SCALARCLASS-LIST takes a list of expressions as its argument. If their
403 ;; scalarclasses all agree, then that scalarclass is returned.
405 (defun scalarclass-list (llist)
406 (cond ((null llist
) nil
)
407 ((null (cdr llist
)) (scalarclass (car llist
)))
408 (t (let ((sc-car (scalarclass (car llist
)))
409 (sc-cdr (scalarclass-list (cdr llist
))))
410 (cond ((or (eq sc-car
'$nonscalar
)
411 (eq sc-cdr
'$nonscalar
))
413 ((and (eq sc-car
'$scalar
) (eq sc-cdr
'$scalar
))
418 (member (caar x
) '(mequal mlist $matrix
))))
420 (defun mequalp (x) (and (not (atom x
)) (eq (caar x
) 'mequal
)))
424 (member (caar x
) '(mlist $matrix
))))
426 (defun mxorlistp1 (x)
428 (or (eq (caar x
) '$matrix
)
429 (and (eq (caar x
) 'mlist
) $listarith
))))
431 (defun constmx (const x
)
432 (simplifya (fmapl1 #'(lambda (ign)
433 (declare (ignore ign
))
438 ;;; ISINOP returns the complete subexpression with the operator OP, when the
439 ;;; operator OP is found in EXPR.
441 (defun isinop (expr op
) ; OP is assumed to be an atom
442 (cond ((atom expr
) nil
)
443 ((and (eq (caar expr
) op
)
444 (not (member 'array
(cdar expr
))))
447 (do ((expr (cdr expr
) (cdr expr
))
450 (when (setq res
(isinop (car expr
) op
))
453 (defun free (exp free-var
)
454 (cond ((alike1 exp free-var
) nil
)
457 (and (listp (car exp
))
458 (free (caar exp
) free-var
)
459 (freel (cdr exp
) free-var
)))))
461 (defun freel (l free-var
)
462 (do ((l l
(cdr l
))) ((null l
) t
)
464 ((atom l
) (return (free l free-var
))) ;; second element of a pair
465 ((not (free (car l
) free-var
)) (return nil
)))))
468 (defun simplifya (x y
)
469 (cond ((not $simp
) x
)
471 (cond ((and $%enumer $numer
(eq x
'$%e
))
472 ;; Replace $%e with its numerical value,
473 ;; when %enumer and $numer TRUE
477 (cond ((and (cdr x
) (atom (cdr x
)))
478 (merror (intl:gettext
"simplifya: malformed expression (atomic cdr).")))
479 ((get (car x
) 'lisp-no-simp
)
480 ;; this feature is to be used with care. it is meant to be
481 ;; used to implement data objects with minimum of consing.
482 ;; forms must not bash the DISPLA package. Only new forms
483 ;; with carefully chosen names should use this feature.
486 (mapcar #'(lambda (x) (simplifya x y
)) (cdr x
))))))
487 ((eq (caar x
) 'rat
) (*red1 x
))
488 ;; Enforced resimplification: Reset dosimp and strip 'simp tags from x.
489 (dosimp (let ((dosimp nil
)) (simplifya (unsimplify x
) y
)))
490 ((member 'simp
(cdar x
)) x
)
491 ((eq (caar x
) 'mrat
) x
)
493 (simplifya (cons (cons ($verbify
(caar x
)) (rest (car x
))) (rest x
)) y
))
494 ((not (atom (caar x
)))
495 (cond ((or (eq (caaar x
) 'lambda
)
496 (and (not (atom (caaar x
))) (eq (caaaar x
) 'lambda
)))
497 (mapply1 (caar x
) (cdr x
) (caar x
) x
))
498 (t (merror (intl:gettext
"simplifya: operator is neither an atom nor a lambda expression: ~S") x
))))
499 ((and $distribute_over
500 (get (caar x
) 'distribute_over
)
501 ;; A function with the property 'distribute_over.
502 ;; Look, if we have a bag as argument to the function.
503 (distribute-over x
)))
504 ((get (caar x
) 'opers
)
505 (let ((opers-list *opers-list
)) (oper-apply x y
)))
506 ((and (eq (caar x
) 'mqapply
)
508 (and (eq substp
'mqapply
)
509 (or (eq (car (cadr x
)) 'lambda
)
510 (eq (caar (cadr x
)) 'lambda
)))))
511 (cond ((or (symbolp (cadr x
)) (not (atom (cadr x
))))
512 (simplifya (cons (cons (cadr x
) (cdar x
)) (cddr x
)) y
))
513 ((or (not (member 'array
(cdar x
))) (not $subnumsimp
))
514 (merror (intl:gettext
"simplifya: I don't know how to simplify this operator: ~M") x
))
516 (t (let ((w (get (caar x
) 'operators
)))
518 (or (not (member 'array
(cdar x
)))
521 (t (simpargs x y
)))))))
523 ;; EQTEST returns an expression which is the same as X
524 ;; except that it is marked with SIMP and maybe other flags from CHECK.
526 ;; Following description is inferred from the code. Dunno why the function is named "EQTEST".
528 ;; (1) if X is already marked with SIMP flag or doesn't need it: return X.
529 ;; (2) if X is pretty much the same as CHECK (same operator and same arguments),
530 ;; then return CHECK after marking it with SIMP flag.
531 ;; (3) if operator of X has the MSIMPIND property, replace it
532 ;; with value of MSIMPIND (something like '(MPLUS SIMP)) and return X.
533 ;; (4) if X or CHECK is an array expression, return X after marking it with SIMP and ARRAY flags.
534 ;; (5) otherwise, return X after marking it with SIMP flag.
536 (defun eqtest (x check
)
541 (member 'simp
(cdar x
)))
543 ((and (eq (caar x
) (caar check
))
544 (equal (cdr x
) (cdr check
)))
545 (cond ((and (null (cdar check
))
546 (setq y
(get (caar check
) 'msimpind
)))
547 (cons y
(cdr check
)))
548 ((member 'simp
(cdar check
))
551 (cons (cons (caar check
)
553 (cons 'simp
(cdar check
))
556 ((setq y
(get (caar x
) 'msimpind
))
558 ((or (member 'array
(cdar x
))
559 (and (eq (caar x
) (caar check
))
560 (member 'array
(cdar check
))))
561 (rplaca x
(cons (caar x
) '(simp array
))))
563 (rplaca x
(cons (caar x
) '(simp)))))))
565 ;; A function, which distributes of bags like a list, matrix, or equation.
566 ;; Check, if we have to distribute of one of the bags or any other operator.
567 (defun distribute-over (expr)
568 (cond ((= 1 (length (cdr expr
)))
569 ;; Distribute over for a function with one argument.
570 (cond ((and (not (atom (cadr expr
)))
571 (member (caaadr expr
) (get (caar expr
) 'distribute_over
)))
574 (mapcar #'(lambda (u) (simplify (list (car expr
) u
)))
578 ;; A function with more than one argument.
579 (do ((args (cdr expr
) (cdr args
))
582 (when (and (not (atom (car args
)))
583 (member (caar (car args
))
584 (get (caar expr
) 'distribute_over
)))
585 ;; Distribute the function over the arguments and simplify again.
588 (cons (ncons (caar (car args
)))
589 (mapcar #'(lambda (u)
593 (cons (ncons (caar expr
))
594 (reverse first-args
))
597 (cdr (car args
)))))))
598 (setq first-args
(cons (car args
) first-args
))))))
600 (defun rulechk (x) (or (mget x
'oldrules
) (get x
'rules
)))
602 (defun resimplify (x) (let ((dosimp t
)) (simplifya x nil
)))
604 (defun unsimplify (x)
605 (if (or (atom x
) (specrepp x
))
607 (cons (remove 'simp
(car x
) :count
1) (mapcar #'unsimplify
(cdr x
)))))
609 (defun simpargs (x y
)
610 (if (or (and (eq (get (caar x
) 'dimension
) 'dimension-infix
)
611 (not (getl (caar x
) '($lassociative $rassociative
))))
612 (get (caar x
) 'binary
))
614 (if (and (member 'array
(cdar x
)) (null (margs x
)))
615 (merror (intl:gettext
"SIMPARGS: subscripted variable found with no subscripts.")))
616 (eqtest (if y x
(let ((flag (member (caar x
) '(mlist mequal
))))
617 (cons (ncons (caar x
))
618 (mapcar #'(lambda (u)
619 (if flag
(simplifya u nil
)
624 ;;;-----------------------------------------------------------------------------
625 ;;; ADDK (X Y) 27.09.2010/DK
627 ;;; Arguments and values:
628 ;;; X - a Maxima number
629 ;;; Y - a Maxima number
630 ;;; result - a simplified Maxima number
633 ;;; ADDK adds two Maxima numbers and returns a simplified Maxima number.
634 ;;; ADDK can be called in Lisp code, whenever the arguments are valid
635 ;;; Maxima numbers, these are integer, float, Maxima rational, or
636 ;;; Maxima bigfloat numbers. The arguments must not be simplified. The
637 ;;; precision of a bigfloat result depends on the setting of the
638 ;;; global variable $FPPREC. If the option variable $FLOAT is T, a
639 ;;; Maxima rational number as a result is converted to a float number.
643 ;;; (addk 2.0 3) -> 5.0
644 ;;; (addk ($bfloat 2) 3)-> ((BIGFLOAT SIMP 56) 45035996273704960 3)
645 ;;; (addk 2 '((rat) 1 2)) -> ((RAT SIMP) 5 2)
646 ;;; (let (($float t)) (addk 2 '((rat) 1 2))) -> 2.5
649 ;;; The option variables $FLOAT and $FPPREC.
652 ;;; TIMESK to multiply and EXPTRL to exponentiate two Maxima numbers.
655 ;;; The routine works for Lisp rational and Lisp complex numbers too.
656 ;;; This feature is not used in Maxima code. If Lisp complex and
657 ;;; rational numbers are mixed with Maxima rational or bigfloat
658 ;;; numbers the result is wrong or a Lisp error is generated.
659 ;;;-----------------------------------------------------------------------------
664 ((and (numberp x
) (numberp y
)) (+ x y
))
665 ((or ($bfloatp x
) ($bfloatp y
)) ($bfloat
(list '(mplus) x y
)))
668 (cond ((floatp x
) (return (+ x
(fpcofrat y
))))
669 (t (setq x
(list '(rat) x
1)))))
671 (cond ((floatp y
) (return (+ y
(fpcofrat x
))))
672 (t (setq y
(list '(rat) y
1))))))
673 (setq g
(gcd (caddr x
) (caddr y
)))
674 (setq a
(truncate (caddr x
) g
)
675 b
(truncate (caddr y
) g
))
676 (return (timeskl (list '(rat) 1 g
)
682 ;;;-----------------------------------------------------------------------------
683 ;;; *RED1 (X) 27.09.2010/DK
686 ;;; Arguments and values:
687 ;;; X - a Maxima rational number (for *RED1)
688 ;;; N - an integer number representing the numerator of a rational
689 ;;; D - an integer number representing the denominator of a rational
690 ;;; result - a simplified Maxima rational number
693 ;;; *RED1 is called from SIMPLIFYA to reduce and simplify a Maxima rational
694 ;;; number. *RED1 checks if the rational number is already simplified. If
695 ;;; the option variable $FLOAT is T, the rational number is converted to a
696 ;;; float number. If the number is not simplified, *RED is called.
698 ;;; *RED reduces the numerator N and the demoniator D and returns a
699 ;;; simplified Maxima rational number. The result is converted to a float
700 ;;; number, if the option variable $FLOAT is T.
703 ;;; The option variable $FLOAT.
704 ;;;-----------------------------------------------------------------------------
707 (cond ((member 'simp
(cdar x
))
708 (if $float
(fpcofrat x
) x
))
709 (t (*red
(cadr x
) (caddr x
)))))
714 (t (let ((u (gcd n d
)))
715 (setq n
(truncate n u
)
717 (if (minusp d
) (setq n
(- n
) d
(- d
)))
718 (cond ((equal d
1) n
)
719 ($float
(fpcofrat1 n d
))
720 (t (list '(rat simp
) n d
)))))))
722 ;;;-----------------------------------------------------------------------------
723 ;;; TIMESK (X Y) 27.09.2010/DK
725 ;;; Arguments and values:
726 ;;; X - a Maxima number
727 ;;; Y - a Maxima number
728 ;;; result - a simplified Maxima number
731 ;;; TIMESK Multiplies two Maxima numbers and returns a simplified Maxima
732 ;;; number. TIMESK can be called in Lisp code, whenever the arguments are
733 ;;; valid Maxima numbers, these are integer, float, Maxima rational, or
734 ;;; Maxima bigfloat numbers. The arguments must not be simplified. The
735 ;;; precision of a bigfloat result depends on the setting of the
736 ;;; global variable $FPPREC. If the option variable $FLOAT is T, a
737 ;;; Maxima rational number as a result is converted to a float number.
739 ;;; TIMESKL is called from TIMESK to multiply two Maxima rational numbers or
740 ;;; a rational number with an integer number.
743 ;;; (timesk 2 3) -> 6
744 ;;; (timesk 2.0 3) -> 6.0
745 ;;; (timesk ($bfloat 2) 3)-> ((BIGFLOAT SIMP 56) 54043195528445952 3)
746 ;;; (timesk 3 '((rat) 1 2)) -> ((RAT SIMP) 3 2)
747 ;;; (let (($float t)) (timesk 3 '((rat) 1 2))) -> 1.5
750 ;;; The option variables $FLOAT and $FPPREC.
753 ;;; ADDK to add and EXPTRL to exponentiate two Maxima numbers.
756 ;;; The routine works for Lisp rational and Lisp complex numbers too.
757 ;;; This feature is not used in Maxima code. If Lisp complex and
758 ;;; rational numbers are mixed with Maxima rational or bigfloat
759 ;;; numbers the result is wrong or a Lisp error is generated.
760 ;;;-----------------------------------------------------------------------------
762 ;; NUM1 and DENOM1 are helper functions for TIMESKL to get the numerator and the
763 ;; denominator of an integer or Maxima rational number. For an integer the
764 ;; denominator is 1. Both functions are used at other places in Maxima code too.
767 (if (numberp a
) a
(cadr a
)))
770 (if (numberp a
) 1 (caddr a
)))
772 (defun timesk (x y
) ; X and Y are assumed to be already reduced
773 (cond ((equal x
1) y
)
775 ((and (numberp x
) (numberp y
)) (* x y
))
776 ((or ($bfloatp x
) ($bfloatp y
)) ($bfloat
(list '(mtimes) x y
)))
777 ((floatp x
) (* x
(fpcofrat y
)))
778 ((floatp y
) (* y
(fpcofrat x
)))
781 ;; TIMESKL takes one or two Maxima rational numbers, one argument can be an
782 ;; integer number. The result is a Maxima rational or an integer number.
783 ;; If the option variable $FLOAT is T, a Maxima rational number in converted
788 (setq u
(*red
(num1 x
) (denom1 y
)))
789 (setq v
(*red
(num1 y
) (denom1 x
)))
790 (setq g
(cond ((or (equal u
0) (equal v
0)) 0)
792 ((and (numberp u
) (numberp v
)) (* u v
))
794 (* (num1 u
) (num1 v
))
795 (* (denom1 u
) (denom1 v
))))))
796 (return (cond ((numberp g
) g
)
797 ((equal (caddr g
) 1) (cadr g
))
798 ($float
(fpcofrat g
))
801 ;;;-----------------------------------------------------------------------------
802 ;;; FPCOFRAT (RATNO) 27.09.2010/DK
805 ;;; Arguments and values:
806 ;;; RATNO - a Maxima rational number (for FPCOFRAT)
807 ;;; NU - an integer number which represents the numerator of a rational
808 ;;; D - an integer number which represents the denominator of a rational
809 ;;; result - floating point approximation of a rational number
812 ;;; Floating Point Conversion OF RATional number routine.
813 ;;; Finds floating point approximation to rational number.
815 ;;; FPCOFRAT1 computes the quotient of NU/D.
817 ;;; Exceptional situations:
818 ;;; A Lisp error is generated, if the rational number does not fit into a
820 ;;;-----------------------------------------------------------------------------
822 ;; This constant is only needed in the file float.lisp.
824 (:compile-toplevel
:load-toplevel
:execute
)
825 (defconstant machine-mantissa-precision
(float-digits 1.0)))
827 (defun fpcofrat (ratno)
828 (fpcofrat1 (cadr ratno
) (caddr ratno
)))
830 (defun fpcofrat1 (nu d
)
833 ;;;-----------------------------------------------------------------------------
834 ;;; EXPTA (X Y) 27.09.2010/DK
836 ;;; Arguments and values:
837 ;;; X - a Maxima number
838 ;;; Y - an integer number
839 ;;; result - a simplified Maxima number
842 ;;; Computes X^Y, where X is Maxima number and Y an integer. The result is
843 ;;; a simplified Maxima number. Y can be a rational Maxima number. For this
844 ;;; case the numerator is taken as the power.
847 ;;; The option variables $FLOAT and $FPPREC.
850 ;;; This routine is not used within the simplifier. There is only one
851 ;;; call from the file hayat.lisp. This call can be replaced with a
852 ;;; call of the function power.
853 ;;;-----------------------------------------------------------------------------
861 ($bfloat
(list '(mexpt) x y
)))
863 (*red
(exptb (caddr x
) (- (num1 y
)))
864 (exptb (cadr x
) (- (num1 y
)))))
866 (*red
(exptb (cadr x
) (num1 y
))
867 (exptb (caddr x
) (num1 y
))))))
869 ;;;-----------------------------------------------------------------------------
870 ;;; EXPTB (A B) 27.09.2010/DK
872 ;;; Arguments and values:
873 ;;; A - a float or integer number
874 ;;; B - an integer number
875 ;;; result - a simplified Maxima number
878 ;;; Computes A^B, where A is a float or an integer number and B is an
879 ;;; integer number. The result is an integer, float, or Maxima
884 ;;; (exptb 3.0 2) -> 9.0
885 ;;; (exptb 3 -2) -> ((RAT SiMP) 1 9)
886 ;;; (let (($float t)) (exptb 3 -2)) -> 0.1111111111111111
889 ;;; The option variable $FLOAT.
892 ;;; EXPTB calls the Lisp functions EXP or EXPT to compute the result.
893 ;;;-----------------------------------------------------------------------------
897 (cond ((equal a %e-val
)
898 ;; Make B a float so we'll get double-precision result.
900 ((or (floatp a
) (not (minusp b
)))
903 (setq b
(expt a
(- b
)))
905 (when (float-inf-p result
) ;; needed for gcl and sbcl - (sometimes) no trap of overflow
906 (error 'floating-point-overflow
))
910 ;;;-----------------------------------------------------------------------------
911 ;;; SIMPLUS (X W Z) 27.09.2010/DK
913 ;;; Arguments and values:
914 ;;; X - a Maxima expression of the form ((mplus) term1 term2 ...)
915 ;;; W - an arbitrary value, the value is ignored
916 ;;; Z - T or NIL, if T the arguments are assumed to be simplified
917 ;;; result - a simplified mplus-expression or an atom
920 ;;; Implementation of the simplifier for the "+" operator.
921 ;;; A general description of SIMPLUS can be found in the paper:
922 ;;; http://www.cs.berkeley.edu/~fateman/papers/simplifier.txt
925 ;;; The addition of matrices and lists is affected by the following option
927 ;;; $DOALLMXOPS, $DOMXMXOPS, $DOMXPLUS, $DOSCMXOPPS, $DOSCMXPLUS, $LISTARITH
930 ;;; This routine should not be called directly. It is called by SIMPLIFYA.
931 ;;; A save access is to call the function ADD.
932 ;;;-----------------------------------------------------------------------------
934 (defun simplus (x w z
)
935 (prog (res check eqnflag matrixflag sumflag
)
936 (if (null (cdr x
)) (return 0))
940 (if (null x
) (go end
))
941 (setq w
(if z
(car x
) (simplifya (car x
) nil
)))
948 (not (member 'trunc
(cdar w
))))
949 (spsimpcases (cdr x
) w
))
950 (setq w
(ratdisrep w
))
955 (nconc (mapcar #'simplify
(cons w
(cdr x
)))
957 ((eq (caar w
) 'mequal
)
962 (add2 (cadr eqnflag
) (cadr w
))
963 (add2 (caddr eqnflag
) (caddr w
)))))
965 ((member (caar w
) '(mlist $matrix
))
967 (cond ((not matrixflag
) w
)
968 ((and (or $doallmxops $domxmxops $domxplus
969 (and (eq (caar w
) 'mlist
)
970 ($listp matrixflag
)))
971 (or (not (eq (caar w
) 'mlist
)) $listarith
))
972 (addmx matrixflag w
))
973 (t (setq res
(pls w res
)) matrixflag
)))
976 (setq sumflag t res
(sumpls w res
))
977 (setq w
(car res
) res
(cdr res
))))
978 (setq res
(pls w res
))
981 (setq res
(testp res
))
984 (cond ((and (or ($listp matrixflag
)
985 $doallmxops $doscmxplus $doscmxops
)
986 (or (not ($listp matrixflag
)) $listarith
))
987 (mxplusc res matrixflag
))
988 (t (testp (pls matrixflag
(pls res nil
)))))))
989 (setq res
(eqtest res check
))
992 (add2 (cadr eqnflag
) res
)
993 (add2 (caddr eqnflag
) res
))
996 ;;;-----------------------------------------------------------------------------
997 ;;; PLS (X OUT) 27.09.2010/DK
999 ;;; Arguments and values:
1000 ;;; X - a Maxima expression or an atom
1001 ;;; OUT - a form ((mplus) <number> term1 term2 ...) or NIL
1002 ;;; result - a form ((mplus) <number> term1 ...), where x is added in.
1005 ;;; Adds the argument X into the form OUT. If OUT is NIL a form
1006 ;;; ((mplus) 0 X) is initialized, if X is an expression or a symbol,
1007 ;;; or ((mplus) X), if X is a number. Numbers are added to the first
1008 ;;; term <number> of the form. Any other symbol or expression is added
1009 ;;; into the canonical ordered list of arguments. The result is in a
1010 ;;; canonical order, but it is not a valid Maxima expression. To get a
1011 ;;; valid Maxima expression the result has to be checked with the
1012 ;;; function TESTP. This is done by the calling routine SIMPLUS.
1014 ;;; PLS checks the global flag *PLUSFLAG*, which is set in PLUSIN to T,
1015 ;;; if a mplus-expression is part of the result.
1018 ;;; (pls 2 nil) -> ((MPLUS) 2)
1019 ;;; (pls '$A nil) -> ((MPLUS) 0 $A)
1020 ;;; (pls '$B '((mplus) 0 $A)) -> ((MPLUS) 0 $A $B)
1021 ;;; (pls '$A '((mplus) 0 $A)) -> ((MPLUS) 0 ((MTIMES SIMP) 2 $A))
1023 ;;; Examples with the option variables $NUMER and $NEGDISTRIB:
1024 ;;; (let (($numer t)) (pls '$%e nil)) -> ((MPLUS) 2.718281828459045)
1025 ;;; (let (($negdistrib t)) (pls '((mtimes) -1 ((mplus) $A $B)) nil))
1026 ;;; -> ((MPLUS) 0 ((MTIMES SIMP) -1 $A) ((MTIMES SIMP) -1 $B))
1027 ;;; (let (($negdistrib nil)) (pls '((mtimes) -1 ((mplus) $A $B)) nil))
1028 ;;; -> ((MPLUS) 0 ((MTIMES) -1 ((MPLUS) $A $B)))
1031 ;;; The option variables $NUMER and $NEGDISTRIB and the global flag
1032 ;;; *PLUSFLAG*, which is set in the routine PLUSIN.
1035 ;;; PLUSIN and ADDK which are called from PLS and SIMPLUS.
1038 ;;; To add an expression into the list (CDR OUT), the list is passed
1039 ;;; to the routine PLUSIN as an argument. PLUSIN adds the argument to
1040 ;;; the list of terms by modifying the list (CDR OUT) destructively.
1041 ;;; The new value of OUT is returned as a result by PLS.
1042 ;;;-----------------------------------------------------------------------------
1044 ;; Set in PLUSIN to T to indicate a nested mplus expression.
1045 (defvar *plusflag
* nil
)
1047 ;; TESTP checks the result of PLS to get a valid Maxima mplus-expression.
1051 ((null (cddr x
)) (cadr x
))
1053 (cond ((null (cdddr x
)) (caddr x
)) (t (rplacd x
(cddr x
)))))
1057 (prog (fm *plusflag
*)
1058 (if (mtimesp x
) (setq x
(testtneg x
)))
1059 (when (and $numer
(atom x
) (eq x
'$%e
))
1060 ;; Replace $%e with its numerical value, when $numer is TRUE
1063 ;; Initialize a form like ((mplus) <number> expr)
1066 (cond ((mnump x
) (ncons x
))
1068 (list 0 (cond ((atom x
) x
) (t (copy-list x
)))))
1069 ((mnump (cadr x
)) (copy-list (cdr x
) ))
1070 (t (cons 0 (copy-list (cdr x
) )))))))
1072 ;; Add a number into the first term of the list out.
1073 (return (cons '(mplus)
1074 (if (mnump (cadr out
))
1075 (cons (addk (cadr out
) x
) (cddr out
))
1076 (cons x
(cdr out
))))))
1077 ((not (mplusp x
)) (plusin x
(cdr out
)) (go end
)))
1078 ;; At this point we have a mplus expression as argument x. The following
1079 ;; code assumes that the argument x is already simplified and the terms
1080 ;; are in a canonical order.
1081 ;; First we add the number to the first term of the list out.
1083 (addk (if (mnump (cadr out
)) (cadr out
) 0)
1084 (cond ((mnump (cadr x
)) (setq x
(cdr x
)) (car x
)) (t 0))))
1085 ;; Initialize fm with the list of terms and start the loop to add the
1086 ;; terms of an mplus expression into the list out.
1089 (if (null (setq x
(cdr x
))) (go end
))
1090 ;; The return value of PLUSIN is a list, where the first element is the
1091 ;; added argument and the rest are the terms which follow the added
1093 (setq fm
(plusin (car x
) fm
))
1096 (if (not *plusflag
*) (return out
))
1097 (setq *plusflag
* nil
) ; *PLUSFLAG* T handles e.g. a+b+3*(a+b)-2*(a+b)
1099 ;; *PLUSFLAG* is set by PLUSIN to indicate that a mplus expression is
1100 ;; part of the result. For this case go again through the terms of the
1101 ;; result and add any term of the mplus expression into the list out.
1104 (when (mplusp (cadr fm
))
1106 (rplacd fm
(cddr fm
))
1110 (if (null (cdr fm
)) (return out
))
1113 ;;;-----------------------------------------------------------------------------
1114 ;;; PLUSIN (X FM) 27.09.2010/DK
1116 ;;; Arguments and values:
1117 ;;; X - a Maxima expression or atom
1118 ;;; FM - a list with the terms of an addition
1119 ;;; result - part of the list fm, which starts at the inserted expression
1122 ;;; Adds X into running list of additive terms FM. The routine modifies
1123 ;;; the argument FM destructively, but does not return the modified list as
1124 ;;; a result. The return value is a part of the list FM, which starts at the
1125 ;;; inserted term. PLUSIN can not handle Maxima numbers. PLUSIN is called
1126 ;;; only from the routine PLS.
1130 ;;; (plusin '$a fm) -> ($A)
1132 ;;; (plusin '$b fm) -> ($B)
1134 ;;; (plusin '$a fm) -> (((MTIMES SIMP) 2 $A) $B)
1135 ;;; fm -> (0 ((MTIMES SIMP) 2 $A) $B)
1138 ;;; Modifies destructively the argument FM, which contains the result of the
1139 ;;; addition of the argument X into the list FM.
1142 ;;; The option variables $doallmxops and $listarith.
1145 ;;; The return value is used in PLS to go in parallel through the list of
1146 ;;; terms, when adding a complete mplus-expression into the list of terms.
1147 ;;; This is triggered by the flag *PLUSFLAG*, which is set in PLUSIN, if
1148 ;;; a mplus-expression is added to the result list.
1149 ;;;-----------------------------------------------------------------------------
1151 (defun plusin (x fm
)
1152 (prog (x1 x2 flag check v w xnew a n m c
)
1157 (if (mnump (cadr x
)) (setq w
(cadr x
) x
(cddr x
))
1159 (t (setq x
(ncons x
))))
1160 (setq x1
(if (null (cdr x
)) (car x
) (cons '(mtimes) x
))
1161 xnew
(list* '(mtimes) w x
))
1163 (cond ((null (cdr fm
)))
1164 ((and (alike1 x1
(cadr fm
)) (null (cdr x
)))
1166 ;; Implement the simplification of
1167 ;; v*a^(c+n)+w*a^(c+m) -> (v*a^n+w*a^m)*a^c
1168 ;; where a, v, w, and (n-m) are integers.
1169 ((and (or (and (mexptp (setq x2
(cadr fm
)))
1172 (not (alike1 x1 x2
))
1174 (integerp (setq v
(cadr x2
)))
1175 (mexptp (setq x2
(caddr x2
)))))
1176 (integerp (setq a
(cadr x2
)))
1179 (integerp (sub (caddr x2
) (caddr x1
))))
1180 (setq n
(if (and (mplusp (caddr x2
))
1181 (mnump (cadr (caddr x2
))))
1183 (if (mnump (caddr x2
))
1186 (setq m
(if (and (mplusp (caddr x1
))
1187 (mnump (cadr (caddr x1
))))
1189 (if (mnump (caddr x1
))
1192 (setq c
(sub (caddr x2
) n
))
1195 ;; n and m are integers and the result is (v*a^n+w*a^m)*a^c.
1196 (setq x1
(mul (addk (timesk v
(exptb a n
))
1197 (timesk w
(exptb a m
)))
1201 ;; n and m are rational numbers: The difference n-m is an
1202 ;; integer. The rational numbers might be improper fractions.
1203 ;; The mixed numbers are: n = n1 + d1/r and m = n2 + d2/r,
1204 ;; where r is the common denominator. We have two cases:
1205 ;; I) d1 = d2: e.g. 2^(1/3+c)+2^(4/3+c)
1206 ;; The result is (v*a^n1+w*a^n2)*a^(c+d1/r)
1207 ;; II) d1 # d2: e.g. 2^(1/2+c)+2^(-1/2+c)
1208 ;; In this case one of the exponents d1 or d2 must
1209 ;; be negative. The negative exponent is factored out.
1210 ;; This guarantees that the factor (v*a^n1+w*a^n2)
1211 ;; is an integer. But the positive exponent has to be
1212 ;; adjusted accordingly. E.g. when we factor out
1213 ;; a^(d2/r) because d2 is negative, then we have to
1214 ;; adjust the positive exponent to n1 -> n1+(d1-d2)/r.
1216 ;; Part of the simplification is done in simptimes. E.g.
1217 ;; this algorithm simplifies the sum sqrt(2)+3*sqrt(2)
1218 ;; to 4*sqrt(2). In simptimes this is further simplified
1220 (multiple-value-bind (n1 d1
)
1221 (truncate (num1 n
) (denom1 n
))
1222 (multiple-value-bind (n2 d2
)
1223 (truncate (num1 m
) (denom1 m
))
1224 (cond ((equal d1 d2
)
1225 ;; Case I: -> (v*a^n1+w*a^n2)*a^(c+d1/r)
1227 (mul (addk (timesk v
(exptb a n1
))
1228 (timesk w
(exptb a n2
)))
1231 (div d1
(denom1 n
))))))
1234 ;; Case II:: d2 is negative, adjust n1.
1235 (setq n1
(add n1
(div (sub d1 d2
) (denom1 n
))))
1237 (mul (addk (timesk v
(exptb a n1
))
1238 (timesk w
(exptb a n2
)))
1241 (div d2
(denom1 n
))))))
1244 ;; Case II: d1 is negative, adjust n2.
1245 (setq n2
(add n2
(div (sub d2 d1
) (denom1 n
))))
1247 (mul (addk (timesk v
(exptb a n1
))
1248 (timesk w
(exptb a n2
)))
1251 (div d1
(denom1 n
))))))
1253 ;; This clause should never be reached.
1254 (t (merror "Internal error in simplus."))))))))
1255 ((mtimesp (cadr fm
))
1256 (cond ((alike1 x1
(cadr fm
))
1258 ((and (mnump (cadadr fm
)) (alike x
(cddadr fm
)))
1259 (setq flag t
) ; found common factor
1261 ((great xnew
(cadr fm
)) (go gr
))))
1262 ((great x1
(cadr fm
)) (go gr
)))
1263 (setq xnew
(eqtest (testt xnew
) (or check
'((foo)))))
1264 (return (cdr (rplacd fm
(cons xnew
(cdr fm
)))))
1271 (list* '(mtimes simp
) 0 x
)
1272 ;; Call muln to get a simplified product.
1273 (if (mtimesp (setq x1
(muln (cons (addk 1 w
) x
) t
)))
1277 (cond ((not (mtimesp (cadr fm
)))
1280 ;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
1281 (rplacd (cadr fm
) (cddadr fm
))
1283 ((not (zerop1 (cadadr fm
)))
1285 ;; Handle the multiplication with a zero.
1286 ((and (or (not $listarith
) (not $doallmxops
))
1287 (mxorlistp (caddr (cadr fm
))))
1289 (cons (constmx 0 (caddr (cadr fm
))) (cddr fm
))))))
1290 ;; (cadadr fm) is zero. If the first term of fm is a number,
1291 ;; add it to preserve the type.
1292 (when (mnump (car fm
))
1293 (rplaca fm
(addk (car fm
) (cadadr fm
))))
1294 (return (rplacd fm
(cddr fm
)))
1296 ;; Call muln to get a simplified product.
1297 (setq x1
(muln (cons (addk w
(if flag
(cadadr fm
) 1)) x
) t
))
1298 ;; Make a mplus expression to guarantee that x1 is added again into the sum
1299 (setq x1
(list '(mplus) x1
))
1303 (list* '(mtimes) x1 x
)
1304 (if (mtimesp x1
) (testtneg x1
) x1
)))
1305 (if (not (mtimesp (cadr fm
))) (go check
))
1306 (when (and (onep (cadadr fm
)) flag
(null (cdddr (cadr fm
))))
1307 ;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
1308 (rplaca (cdr fm
) (caddr (cadr fm
))) (go check
))
1311 (if (mplusp (cadr fm
)) (setq *plusflag
* t
)) ; A nested mplus expression
1314 ;;;-----------------------------------------------------------------------------
1316 ;; Routines to add matrices
1318 (defun mxplusc (sc mx
)
1320 (setq sc
(partition-ns (cdr sc
)))
1321 (cond ((null (car sc
)) (cons '(mplus) (cons mx
(cadr sc
))))
1322 ((not (null (cadr sc
)))
1325 (outermap1 'mplus
(cons '(mplus) (car sc
)) mx
))
1327 (t (simplify (outermap1 'mplus
(cons '(mplus) (car sc
)) mx
)))))
1328 ((not (scalar-or-constant-p sc $assumescalar
))
1329 (testp (pls mx
(pls sc nil
))))
1330 (t (simplify (outermap1 'mplus sc mx
)))))
1332 (defun partition-ns (x)
1333 (let (sp nsp
) ; SP = scalar part, NSP = nonscalar part
1334 (mapc #'(lambda (z) (if (scalar-or-constant-p z $assumescalar
)
1335 (setq sp
(cons z sp
))
1336 (setq nsp
(cons z nsp
))))
1338 (list (nreverse sp
) (nreverse nsp
))))
1340 (defun addmx (x1 x2
)
1341 (let (($doscmxops t
) ($domxmxops t
) ($listarith t
))
1342 (simplify (fmapl1 'mplus x1 x2
))))
1344 ;;; ----------------------------------------------------------------------------
1346 ;;; Simplification of the Log function
1348 ;; The log function distributes over lists, matrices, and equations
1349 (defprop %log
(mlist $matrix mequal
) distribute_over
)
1351 (def-simplifier log
(y)
1352 (cond ((onep1 y
) (addk -
1 y
))
1354 (cond (radcanp (list '(%log simp
) 0))
1356 (merror (intl:gettext
"log: encountered log(0).")))
1357 (t (throw 'errorsw t
))))
1358 ;; Check evaluation in floating point precision.
1359 ((flonum-eval (mop form
) y
))
1360 ;; Check evaluation in bigfloag precision.
1361 ((and (not (member 'simp
(car form
)))
1362 (big-float-eval (mop form
) y
)))
1365 (cond ((or (and $logexpand
(eq $domain
'$real
))
1366 (member $logexpand
'($all $super
))
1367 (and (eq ($csign
(cadr y
)) '$pos
)
1368 (not (member ($csign
(caddr y
))
1369 '($complex $imaginary
)))))
1370 ;; Simplify log(x^a) -> a*log(x), where x > 0 and a is real
1371 (mul (caddr y
) (take '(%log
) (cadr y
))))
1372 ((or (and (ratnump (caddr y
))
1373 (or (eql 1 (cadr (caddr y
)))
1374 (eql -
1 (cadr (caddr y
)))))
1375 (maxima-integerp (inv (caddr y
))))
1376 ;; Simplify log(z^(1/n)) -> log(z)/n, where n is an integer
1378 (take '(%log
) (cadr y
))))
1379 ((and (eq (cadr y
) '$%e
)
1380 (or (not (member ($csign
(caddr y
))
1381 '($complex $imaginary
)))
1382 (not (member ($csign
(mul '$%i
(caddr y
)))
1383 '($complex $imaginary
)))))
1384 ;; Simplify log(exp(x)) and log(exp(%i*x)), where x is a real
1388 ;; Simplify log(n/d)
1389 (cond ((eql (cadr y
) 1)
1390 (mul -
1 (take '(%log
) (caddr y
))))
1391 ((eq $logexpand
'$super
)
1392 (sub (take '(%log
) (cadr y
)) (take '(%log
) (caddr y
))))
1394 ((and (member $logexpand
'($all $super
))
1396 (do ((y (cdr y
) (cdr y
))
1398 ((null y
) (return (addn b t
)))
1399 (setq b
(cons (take '(%log
) (car y
)) b
))))
1400 ((and (member $logexpand
'($all $super
))
1402 (member (caar y
) '(%product $product
)))
1403 (let ((new-op (if (char= (get-first-char (caar y
)) #\%
) '%sum
'$sum
)))
1404 (simplifya `((,new-op
) ((%log
) ,(cadr y
)) ,@(cddr y
)) t
)))
1407 (eq ($sign y
) '$neg
))
1408 (add (mul '$%i
'$%pi
) (take '(%log
) (neg y
))))
1409 ((taylorize (mop form
) (second form
)))
1413 (simplifya (list '(mtimes) (caddr w
)
1414 (simplifya (list '(%log
) (cadr w
)) t
)) t
))
1416 ;;; ----------------------------------------------------------------------------
1418 ;;; Implementation of the Square root function
1420 (defprop $sqrt %sqrt verb
)
1421 (defprop $sqrt %sqrt alias
)
1423 (defprop %sqrt $sqrt noun
)
1424 (defprop %sqrt $sqrt reversealias
)
1426 (defprop %sqrt simp-sqrt operators
)
1429 (simplify (list '(%sqrt
) z
)))
1431 (defun simp-sqrt (x ignored z
)
1432 (declare (ignore ignored
))
1434 (simplifya (list '(mexpt) (cadr x
) '((rat simp
) 1 2)) z
))
1436 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1438 ;;; Simplification of the "/" operator.
1440 (defun simpquot (x y z
)
1442 (cond ((and (integerp (cadr x
)) (integerp (caddr x
)) (not (zerop (caddr x
))))
1443 (*red
(cadr x
) (caddr x
)))
1444 ((and (numberp (cadr x
)) (numberp (caddr x
)) (not (zerop (caddr x
))))
1445 (/ (cadr x
) (caddr x
)))
1446 ((and (floatp (cadr x
)) (floatp (caddr x
)) #-ieee-floating-point
(not (zerop (caddr x
))))
1447 (/ (cadr x
) (caddr x
)))
1448 ((and ($bfloatp
(cadr x
)) ($bfloatp
(caddr x
)) (not (equal bigfloatzero
(caddr x
))))
1449 ;; Call BIGFLOATP to ensure that arguments have same precision.
1450 ;; Otherwise FPQUOTIENT could return a spurious value.
1451 (bcons (fpquotient (cdr (bigfloatp (cadr x
))) (cdr (bigfloatp (caddr x
))))))
1452 (t (setq y
(simplifya (cadr x
) z
))
1453 (setq x
(simplifya (list '(mexpt) (caddr x
) -
1) z
))
1454 (if (equal y
1) x
(simplifya (list '(mtimes) y x
) t
)))))
1456 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1458 ;;; Implementation of the abs function.
1460 ;; Put the properties alias, reversealiases, noun and verb on the property list.
1461 (defprop $abs mabs alias
)
1462 (defprop $abs mabs verb
)
1463 (defprop mabs $abs reversealias
)
1464 (defprop mabs $abs noun
)
1466 ;; The abs function distributes over bags.
1467 (defprop mabs
(mlist $matrix mequal
) distribute_over
)
1469 ;; Define a verb function $abs
1471 (simplify (list '(mabs) x
)))
1473 ;; The abs function is a simplifying function.
1474 (defprop mabs simpabs operators
)
1476 (defun simpabs (e y z
)
1477 (declare (ignore y
))
1480 (x (simpcheck (second e
) z
)))
1482 (cond ((complex-number-p x
#'(lambda (s) (or (floatp s
) ($bfloatp s
))))
1483 (maxima::to
(bigfloat::abs
(bigfloat:to x
))))
1485 ((complex-number-p x
#'mnump
)
1488 ;; nounform for arrays...
1489 ((or (arrayp x
) ($member x $arrays
)) `((mabs simp
) ,x
))
1491 ;; taylor polynomials
1492 ((taylorize 'mabs x
))
1494 ;; values for extended real arguments:
1495 ((member x
'($inf $infinity $minf
)) '$inf
)
1496 ((member x
'($ind $und
)) x
)
1498 ;; abs(abs(expr)) --> abs(expr). Since x is simplified, it's OK to return x.
1499 ((and (consp x
) (consp (car x
)) (eq (caar x
) 'mabs
))
1502 ;; abs(conjugate(expr)) = abs(expr).
1503 ((and (consp x
) (consp (car x
)) (eq (caar x
) '$conjugate
))
1504 (take '(mabs) (cadr x
)))
1507 (setq sgn
($csign x
))
1508 (cond ((member sgn
'($neg $nz
)) (mul -
1 x
))
1509 ((eq '$zero sgn
) (mul 0 x
))
1510 ((member sgn
'($pos $pz
)) x
)
1512 ;; for complex constant expressions, use $cabs
1513 ((and (eq sgn
'$complex
) ($constantp x
))
1516 ;; abs(pos^complex) --> pos^(realpart(complex)).
1517 ((and (eq sgn
'$complex
) (mexptp x
) (eq '$pos
($csign
(second x
))))
1518 (power (second x
) ($realpart
(third x
))))
1520 ;; for abs(neg^z), use cabs.
1521 ((and (mexptp x
) (eq '$neg
($csign
(second x
))))
1524 ;; When x # 0, we have abs(signum(x)) = 1.
1525 ((and (eq '$pn sgn
) (consp x
) (consp (car x
)) (eq (caar x
) '%signum
)) 1)
1527 ;; multiplicative property: abs(x*y) = abs(x) * abs(y). We would like
1528 ;; assume(a*b > 0), abs(a*b) --> a*b. Thus the multiplicative property
1529 ;; is applied after the sign test.
1531 (muln (mapcar #'(lambda (u) (take '(mabs) u
)) (margs x
)) t
))
1533 ;; abs(x^n) = abs(x)^n for integer n. Is the featurep check worthwhile?
1534 ;; Again the sign check is done first because we'd like abs(x^2) --> x^2.
1535 ((and (mexptp x
) ($featurep
(caddr x
) '$integer
))
1536 (power (take '(mabs) (cadr x
)) (caddr x
)))
1538 ;; Reflection rule: abs(-x) --> abs(x)
1539 ((great (neg x
) x
) (take '(mabs) (neg x
)))
1542 (t (eqtest (list '(mabs) x
) e
)))))))
1544 (defun abs-integral (x)
1545 (mul (div 1 2) x
(take '(mabs) x
)))
1547 (putprop 'mabs
`((x) ,'abs-integral
) 'integral
)
1549 ;; I (rtoy) think this does some simple optimizations of x * y.
1554 ;; We have something like ((mtimes) foo). This is the same as foo.
1557 ;; We have 1*foo. Which is the same as foo. This should not
1558 ;; be applied to 1.0 or 1b0!
1559 (cond ((null (cdddr x
))
1561 (t (rplacd x
(cddr x
)))))
1565 ;; This basically converts -(a+b) to -a-b.
1567 (cond ((and (equal (cadr x
) -
1)
1571 ;; If x is exactly of the form -1*(sum), and $negdistrib is
1572 ;; true, we distribute the -1 across the sum.
1573 (addn (mapcar #'(lambda (z)
1579 ;; Simplification of the "-" operator
1580 (defun simpmin (x vestigial z
)
1581 (declare (ignore vestigial
))
1582 (cond ((null (cdr x
)) 0)
1584 (mul -
1 (simplifya (cadr x
) z
)))
1586 ;; ((mminus) a b ...) -> ((mplus) a ((mtimes) -1 b) ...)
1587 (sub (simplifya (cadr x
) z
) (addn (cddr x
) z
)))))
1589 (defun simptimes (x w z
) ; W must be 1
1590 (prog (res check eqnflag matrixflag sumflag
)
1591 (if (null (cdr x
)) (return 1))
1597 (cond ((and matrixflag
(mxorlistp1 matrixflag
))
1598 (return (constmx res matrixflag
)))
1599 (eqnflag (return (list '(mequal simp
)
1600 (mul2 res
(cadr eqnflag
))
1601 (mul2 res
(caddr eqnflag
)))))
1604 (cond ((mxorlistp u
)
1605 (return (setq res
(constmx res u
))))
1607 (mxorlistp1 (cadr u
))
1608 ($numberp
(caddr u
)))
1609 (return (setq res
(constmx res
(cadr u
)))))
1613 (list '(mequal simp
)
1615 (mul2 res
(caddr u
))))))))))))
1617 ((null x
) (go end
)))
1618 (setq w
(if z
(car x
) (simplifya (car x
) nil
)))
1620 (cond ((atom w
) nil
)
1621 ((eq (caar w
) 'mrat
)
1622 (cond ((or eqnflag matrixflag
1624 (not (member 'trunc
(cdar w
))))
1625 (spsimpcases (cdr x
) w
))
1626 (setq w
(ratdisrep w
))
1630 (ratf (cons '(mtimes)
1631 (nconc (mapcar #'simplify
(cons w
(cdr x
)))
1633 ((eq (caar w
) 'mequal
)
1638 (mul2 (cadr eqnflag
) (cadr w
))
1639 (mul2 (caddr eqnflag
) (caddr w
)))))
1641 ((member (caar w
) '(mlist $matrix
))
1643 (cond ((not matrixflag
) w
)
1644 ((and (or $doallmxops $domxmxops $domxtimes
)
1645 (or (not (eq (caar w
) 'mlist
)) $listarith
)
1646 (not (eq *inv
* '$detout
)))
1647 (stimex matrixflag w
))
1648 (t (setq res
(tms (copy-tree w
) 1 (copy-tree res
))) matrixflag
)))
1650 ((and (eq (caar w
) '%sum
) $sumexpand
)
1651 (setq sumflag
(sumtimes sumflag w
))
1653 (setq res
(tms (copy-tree w
) 1 (copy-tree res
)))
1656 (cond ((mtimesp res
) (setq res
(testt res
))))
1657 (cond (sumflag (setq res
(cond ((or (null res
) (equal res
1)) sumflag
)
1658 ((not (mtimesp res
))
1659 (list '(mtimes) res sumflag
))
1660 (t (nconc res
(list sumflag
)))))))
1661 (cond ((or (atom res
)
1662 (not (member (caar res
) '(mexpt mtimes
)))
1663 (and (zerop $expop
) (zerop $expon
))
1665 ((eq (caar res
) 'mtimes
) (setq res
(expandtimes res
)))
1666 ((and (mplusp (cadr res
))
1667 (fixnump (caddr res
))
1668 (not (or (> (caddr res
) $expop
)
1669 (> (- (caddr res
)) $expon
))))
1670 (setq res
(expandexpt (cadr res
) (caddr res
)))))
1673 (cond ((null res
) matrixflag
)
1674 ((and (or ($listp matrixflag
)
1677 (not (member res
'(-1 -
1.0))))
1678 ;; RES should only be -1 here (not = 1)
1680 (member res
'(-1 -
1.0))))
1681 (or (not ($listp matrixflag
)) $listarith
))
1682 (mxtimesc res matrixflag
))
1683 (t (testt (tms matrixflag
1 (tms res
1 nil
))))))))
1684 (if res
(setq res
(eqtest res check
)))
1685 (return (cond (eqnflag
1686 (if (null res
) (setq res
1))
1688 (mul2 (cadr eqnflag
) res
)
1689 (mul2 (caddr eqnflag
) res
)))
1692 (defun spsimpcases (l e
)
1694 (if (or (mbagp u
) (and (not (atom u
))
1696 (not (member 'trunc
(cdar e
)))))
1699 (defun mxtimesc (sc mx
)
1701 (and (mtimesp sc
) (member (cadr sc
) '(-1 -
1.0))
1702 $doscmxops
(not (or $doallmxops $domxmxops $domxtimes
))
1703 (setq sign
(cadr sc
)) (rplaca (cdr sc
) nil
))
1704 (setq out
(let ((scp* (cond ((mtimesp sc
) (partition-ns (cdr sc
)))
1705 ((not (scalar-or-constant-p sc $assumescalar
))
1708 (cond ((null scp
*) (list '(mtimes simp
) sc mx
))
1709 ((and (not (atom scp
*)) (null (car scp
*)))
1710 (append '((mtimes)) (cadr scp
*) (list mx
)))
1711 ((or (atom scp
*) (and (null (cdr scp
*))
1712 (not (null (cdr sc
)))
1713 (setq scp
* (cons '(mtimes) (car scp
*))))
1715 (simplifya (outermap1 'mtimes scp
* mx
) nil
))
1716 (t (append '((mtimes))
1719 (cons '(mtimes) (car scp
*)) mx
)
1722 (cond (sign (if (mtimesp out
)
1723 (rplacd out
(cons sign
(cdr out
)))
1724 (list '(mtimes) sign out
)))
1725 ((mtimesp out
) (testt out
))
1729 (let (($doscmxops t
) ($domxmxops t
) ($listarith t
))
1730 (simplify (fmapl1 'mtimes x y
))))
1732 ;; TMS takes a simplified expression FACTOR and a cumulative
1733 ;; PRODUCT as arguments and modifies the cumulative product so
1734 ;; that the expression is now one of its factors. The
1735 ;; exception to this occurs when a tellsimp rule is triggered.
1736 ;; The second argument is the POWER to which the expression is
1737 ;; to be raised within the product.
1739 (defun tms (factor power product
&aux tem
)
1740 (let ((*rulesw
* nil
)
1742 (when (mplusp product
) (setq product
(list '(mtimes simp
) product
)))
1743 (cond ((zerop1 factor
)
1744 (cond ((mnegp power
)
1747 (merror (intl:gettext
"Division by 0"))))
1749 ((and (null product
)
1750 (or (and (mtimesp factor
) (equal power
1))
1751 (and (setq product
(list '(mtimes) 1)) nil
)))
1752 (setq tem
(append '((mtimes)) (if (mnump (cadr factor
)) nil
'(1))
1754 (if (= (length tem
) 1)
1755 (setq tem
(copy-list tem
))
1758 (do ((factor-list (cdr factor
) (cdr factor-list
)))
1759 ((or (null factor-list
) (zerop1 product
)) product
)
1760 (setq z
(timesin (car factor-list
) (cdr product
) power
))
1763 (setq product
(tms-format-product z
)))))
1765 (setq z
(timesin factor
(cdr product
) power
))
1767 (tms-format-product z
)
1770 (defun tms-format-product (x)
1771 (cond ((zerop1 x
) x
)
1772 ((mnump x
) (list '(mtimes) x
))
1773 ((not (mtimesp x
)) (list '(mtimes) 1 x
))
1774 ((not (mnump (cadr x
))) (cons '(mtimes) (cons 1 (cdr x
))))
1778 (cond ($ratsimpexpons
(sratsimp (list '(mplus) x y
)))
1779 ((and (mnump x
) (mnump y
)) (addk x y
))
1783 (if (and (mnump x
) (mnump y
))
1787 (defun simp-limit (x vestigial z
)
1788 (declare (ignore vestigial
))
1789 (let ((l1 (length x
))
1791 (unless (or (= l1
2) (= l1
4) (= l1
5))
1792 (merror (intl:gettext
"limit: wrong number of arguments.")))
1793 (setq y
(simpmap (cdr x
) z
))
1794 (cond ((and (= l1
5) (not (member (cadddr y
) '($plus $minus
))))
1795 (merror (intl:gettext
"limit: direction must be either 'plus' or 'minus': ~M") (cadddr y
)))
1797 (merror (intl:gettext
"limit: variable must not be a number; found: ~M") (cadr y
)))
1801 (eqtest (cons '(%limit
) y
) x
)))))
1803 (defun simpinteg (x vestigial z
)
1804 (declare (ignore vestigial
))
1805 (let ((l1 (length x
))
1807 (unless (or (= l1
3) (= l1
5))
1808 (merror (intl:gettext
"integrate: wrong number of arguments.")))
1809 (setq y
(simpmap (cdr x
) z
))
1810 (cond ((mnump (cadr y
))
1811 (merror (intl:gettext
"integrate: variable must not be a number; found: ~M") (cadr y
)))
1812 ((and (= l1
5) (alike1 (caddr y
) (cadddr y
)))
1815 (free (setq z
(sub (cadddr y
) (caddr y
))) '$%i
)
1816 (eq ($sign z
) '$neg
))
1817 (neg (simplifya (list '(%integrate
) (car y
) (cadr y
) (cadddr y
) (caddr y
)) t
)))
1821 (if (or (among '$inf z
) (among '$minf z
))
1825 (eqtest (cons '(%integrate
) y
) x
)))))
1827 (defun simpbigfloat (x vestigial simp-flag
)
1828 (declare (ignore vestigial simp-flag
))
1831 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1833 ;;; Implementation of the Exp function.
1835 (defprop $exp %exp verb
)
1836 (defprop $exp %exp alias
)
1838 (defprop %exp $exp noun
)
1839 (defprop %exp $exp reversealias
)
1841 (defprop %exp simp-exp operators
)
1844 (simplify (list '(%exp
) z
)))
1846 ;; Support a function for code,
1847 ;; which depends on an unsimplified noun form.
1848 (defmfun $exp-form
(z)
1849 (list '(mexpt) '$%e z
))
1851 (defun simp-exp (x ignored z
)
1852 (declare (ignore ignored
))
1854 (simplifya (list '(mexpt) '$%e
(cadr x
)) z
))
1856 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1858 (defun simplambda (x vestigial simp-flag
)
1859 (declare (ignore vestigial simp-flag
))
1860 ; Check for malformed lambda expressions.
1861 ; We verify that we have a valid list of parameters and a non-empty body.
1862 (let ((params (cadr x
)))
1863 (unless ($listp params
)
1864 (merror (intl:gettext
"lambda: first argument must be a list; found: ~M") params
))
1865 (do ((params (cdr params
) (cdr params
))
1868 (when (mdeflistp params
)
1869 (setq params
(cdar params
)))
1870 (let ((p (car params
)))
1871 (unless (or (mdefparam p
)
1872 (and (op-equalp p
'mquote
)
1873 (mdefparam (cadr p
))))
1874 (merror (intl:gettext
"lambda: parameter must be a symbol and must not be a system constant; found: ~M") p
))
1876 (when (member p seen-params
:test
#'eq
)
1877 (merror (intl:gettext
"lambda: ~M occurs more than once in the parameter list") p
))
1878 (push p seen-params
))))
1879 (when (null (cddr x
))
1880 (merror (intl:gettext
"lambda: no body present")))
1881 (cons '(lambda simp
) (cdr x
)))
1883 (defun simpmdef (x vestigial simp-flag
)
1884 (declare (ignore vestigial simp-flag
))
1886 (cons '(mdefine simp
) (cdr x
)))
1888 (defun simpmap (e z
)
1889 (mapcar #'(lambda (u) (simpcheck u z
)) e
))
1892 (let ((x ($expand e
1 1)))
1893 (cond ((or (not (free x
'$ind
)) (not (free x
'$und
))
1894 (not (free x
'$zeroa
)) (not (free x
'$zerob
))
1895 (not (free x
'$infinity
))
1898 ((and (free x
'$inf
) (free x
'$minf
)) x
)
1899 (t (infsimp1 x e
)))))
1901 (defun infsimp1 (x e
)
1902 (let ((minf-coef (coeff x
'$minf
1))
1903 (inf-coef (coeff x
'$inf
1)))
1904 (cond ((or (and (equal minf-coef
0)
1906 (and (not (free minf-coef
'$inf
))
1907 (not (free inf-coef
'$minf
)))
1908 (let ((new-exp (sub (add2 (mul2 minf-coef
'$minf
)
1909 (mul2 inf-coef
'$inf
))
1911 (and (not (free new-exp
'$inf
))
1912 (not (free new-exp
'$minf
)))))
1914 (t (let ((sign-minf-coef ($asksign minf-coef
))
1915 (sign-inf-coef ($asksign inf-coef
)))
1916 (cond ((or (and (eq sign-inf-coef
'$zero
)
1917 (eq sign-minf-coef
'$neg
))
1918 (and (eq sign-inf-coef
'$pos
)
1919 (eq sign-minf-coef
'$zero
))
1920 (and (eq sign-inf-coef
'$pos
)
1921 (eq sign-minf-coef
'$neg
))) '$inf
)
1922 ((or (and (eq sign-inf-coef
'$zero
)
1923 (eq sign-minf-coef
'$pos
))
1924 (and (eq sign-inf-coef
'$neg
)
1925 (eq sign-minf-coef
'$zero
))
1926 (and (eq sign-inf-coef
'$neg
)
1927 (eq sign-minf-coef
'$pos
))) '$minf
)
1928 ((or (and (eq sign-inf-coef
'$pos
)
1929 (eq sign-minf-coef
'$pos
))
1930 (and (eq sign-inf-coef
'$neg
)
1931 (eq sign-minf-coef
'$neg
))) '$und
)))))))
1933 (defun infsimp2 (x e
)
1935 (if (isinop x
'%limit
) e x
))
1937 (defun simpderiv (x y z
)
1939 (cond ((not (even (length x
)))
1940 (cond ((and (cdr x
) (null (cdddr x
))) (nconc x
'(1)))
1941 (t (wna-err '%derivative
)))))
1942 (setq w
(cons '(%derivative
) (simpmap (cdr x
) z
)))
1944 (do ((u (cddr w
) (cddr u
))) ((null u
))
1945 (cond ((mnump (car u
))
1946 (merror (intl:gettext
"diff: variable must not be a number; found: ~M") (car u
)))))
1947 (cond ((or (zerop1 y
)
1948 (and (or (mnump y
) (and (atom y
) (constant y
)))
1950 (and (not (alike1 y
(caddr w
)))
1951 (do ((u (cddr w
) (cddr u
))) ((null u
))
1952 (cond ((and (numberp (cadr u
))
1953 (not (zerop (cadr u
))))
1956 ((and (not (atom y
)) (eq (caar y
) '%derivative
) derivsimp
)
1957 (rplacd w
(append (cdr y
) (cddr w
)))))
1959 (return (if (null derivflag
) (list '(%del simp
) y
) (deriv (cdr w
)))))
1962 (cond ((null u
) (go next
))
1963 ((zerop1 (caddr u
)) (rplacd u
(cdddr u
)))
1964 (t (setq u
(cddr u
))))
1967 (cond ((null (cddr w
)) (return y
))
1968 ((and (null (cddddr w
))
1970 (alike1 (cadr w
) (caddr w
)))
1975 (cond ((null (cddr z
)) (go loop
))
1976 ((alike1 (car z
) (caddr z
))
1977 (rplaca (cdddr z
) (add2 (cadr z
) (cadddr z
)))
1978 (rplacd z
(cdddr z
)))
1979 ((great (car z
) (caddr z
))
1980 (let ((u1 (car z
)) (u2 (cadr z
)) (v1 (caddr z
)) (v2 (cadddr z
)))
1983 (rplacd z
(cons v2
(cons u1
(cons u2
(cddddr z
))))))))
1984 (cond ((setq z
(cddr z
)) (go sort
)))
1986 (cond ((null flag
) (return (cond ((null derivflag
) (eqtest w x
))
1987 (t (deriv (cdr w
)))))))
1993 (setq x
(num1 x
)) (cond ((plusp x
) 1) ((minusp x
) -
1) (t 0)))
1995 ((mplusp x
) (if *expandp
* 1 (signum1 (car (last x
)))))
1996 ((mtimesp x
) (if (mplusp (cadr x
)) 1 (signum1 (cadr x
))))
1999 (defprop %signum
(mlist $matrix mequal
) distribute_over
)
2002 (defun simpsignum (e y z
)
2003 (declare (ignore y
))
2005 (let ((x (simpcheck (second e
) z
)) (sgn))
2007 (cond ((complex-number-p x
#'mnump
)
2008 (if (complex-number-p x
#'$ratnump
) ;; nonfloat complex
2009 (if (zerop1 x
) 0 ($rectform
(div x
($cabs x
))))
2010 (maxima::to
(bigfloat::signum
(bigfloat::to x
)))))
2012 ;; idempotent: signum(signum(z)) = signum(z).
2013 ((and (consp x
) (consp (car x
)) (eq '%signum
(mop x
))) x
)
2016 (setq sgn
($csign x
))
2017 (cond ((eq sgn
'$neg
) -
1)
2021 ;; multiplicative: signum(ab) = signum(a) * signum(b).
2023 (muln (mapcar #'(lambda (s) (take '(%signum
) s
)) (margs x
)) t
))
2025 ;; Reflection rule: signum(-x) --> -signum(x).
2026 ((great (neg x
) x
) (neg (take '(%signum
) (neg x
))))
2029 (t (eqtest (list '(%signum
) x
) e
)))))))
2031 (def-simplifier signum
(x)
2033 (cond ((complex-number-p x
#'mnump
)
2034 (if (complex-number-p x
#'$ratnump
) ;; nonfloat complex
2037 ($rectform
(div x
($cabs x
))))
2038 (maxima::to
(bigfloat::signum
(bigfloat::to x
)))))
2040 ;; idempotent: signum(signum(z)) = signum(z).
2043 (eq '%signum
(mop x
)))
2047 (setq sgn
($csign x
))
2048 (cond ((eq sgn
'$neg
) -
1)
2052 ;; multiplicative: signum(ab) = signum(a) * signum(b).
2054 (muln (mapcar #'(lambda (s) (take '(%signum
) s
)) (margs x
))
2057 ;; Reflection rule: signum(-x) --> -signum(x).
2059 (neg (take '(%signum
) (neg x
))))
2064 (defun exptrl (r1 r2
)
2065 (cond ((equal r2
1) r1
)
2067 (cond ((mnump r1
) (addk 0.0 r1
))
2068 ;; Do not simplify the type of the number away.
2069 (t (list '(mexpt simp
) r1
1.0))))
2070 ((equal r2 bigfloatone
)
2071 (cond ((mnump r1
) ($bfloat r1
))
2072 ;; Do not simplify the type of the number away.
2073 (t (list '(mexpt simp
) r1 bigfloatone
))))
2075 (cond ((or (zerop1 r2
) (mnegp r2
))
2077 (merror (intl:gettext
"expt: undefined: ~M") (list '(mexpt) r1 r2
))
2078 (throw 'errorsw t
)))
2079 (t (zerores r1 r2
))))
2080 ((or (zerop1 r2
) (onep1 r1
))
2081 (cond ((or ($bfloatp r1
) ($bfloatp r2
)) bigfloatone
)
2082 ((or (floatp r1
) (floatp r2
)) 1.0)
2084 ((or ($bfloatp r1
) ($bfloatp r2
)) ($bfloat
(list '(mexpt) r1 r2
)))
2085 ((and (numberp r1
) (integerp r2
)) (exptb r1 r2
))
2086 ((and (numberp r1
) (floatp r2
) (equal r2
(float (floor r2
))))
2087 (exptb (float r1
) (floor r2
)))
2088 ((or $numer
(and (floatp r2
) (or (plusp (num1 r1
)) $numer_pbranch
)))
2090 (cond ((minusp (setq r1
(addk 0.0 r1
)))
2091 (cond ((or $numer_pbranch
(eq $domain
'$complex
))
2093 ;; R1^R2 = (-R1)^R2*cos(pi*R2) + i*(-R1)^R2*sin(pi*R2)
2094 (setq r2
(addk 0.0 r2
))
2095 (setq y
(exptrl (- r1
) r2
) r2
(* (mget '$%pi
'$numer
) r2
))
2096 (add2 (* y
(cos r2
))
2097 (list '(mtimes simp
) (* y
(sin r2
)) '$%i
)))
2098 (t (setq y
(let ($numer $float $keepfloat $ratprint
)
2100 (mul2 y
(exptrl (- r1
) r2
)))))
2101 ((equal (setq r2
(addk 0.0 r2
)) (float (floor r2
)))
2102 (exptb r1
(floor r2
)))
2103 ((and (equal (setq y
(* 2.0 r2
)) (float (floor y
)))
2104 (not (equal r1 %e-val
)))
2105 (exptb (sqrt r1
) (floor y
)))
2106 (t (exp (* r2
(log r1
)))))))
2107 ((floatp r2
) (list '(mexpt simp
) r1 r2
))
2110 (exptrl (cond ((equal (abs (cadr r1
)) 1)
2111 (* (cadr r1
) (caddr r1
)))
2112 ;; We set the simp flag at this place. This
2113 ;; changes nothing for an exponent r2 # -1.
2114 ;; exptrl is called again and does not look at
2115 ;; the simp flag. For the case r2 = -1 exptrl
2116 ;; is called with an exponent 1. For this case
2117 ;; the base is immediately returned. Now the
2118 ;; base has the correct simp flag. (DK 02/2010)
2120 (list '(rat simp
) (- (caddr r1
)) (- (cadr r1
))))
2121 (t (list '(rat simp
) (caddr r1
) (cadr r1
))))
2123 (t (list '(rat simp
) (exptb (cadr r1
) r2
) (exptb (caddr r1
) r2
)))))
2124 ((and (floatp r1
) (alike1 r2
'((rat) 1 2)))
2126 (list '(mtimes simp
) (sqrt (- r1
)) '$%i
)
2128 ((and (floatp r1
) (alike1 r2
'((rat) -
1 2)))
2130 (list '(mtimes simp
) (/ -
1.0 (sqrt (- r1
))) '$%i
)
2134 (exptrl r1
(fpcofrat r2
))
2135 (mul2 (exptrl -
1 r2
) ;; (-4.5)^(1/4) -> (4.5)^(1/4) * (-1)^(1/4)
2136 (exptrl (- r1
) r2
))))
2137 (*exptrlsw
* (list '(mexpt simp
) r1 r2
))
2139 (let ((*exptrlsw
* t
))
2140 (simptimes (list '(mtimes)
2141 (exptrl r1
(truncate (cadr r2
) (caddr r2
)))
2142 (let ((y (let ($keepfloat $ratprint
)
2143 (simpnrt r1
(caddr r2
))))
2144 (z (rem (cadr r2
) (caddr r2
))))
2146 (list (car y
) (cadr y
) (mul2 (caddr y
) z
))
2150 (defun simpexpt (x y z
)
2151 (prog (gr pot check res
*rulesw
* w mlpgr mlppot
)
2153 (cond (z (setq gr
(cadr x
) pot
(caddr x
)) (go cont
)))
2155 (setq gr
(simplifya (cadr x
) nil
))
2157 (let (($%enumer $numer
))
2158 ;; Switch $%enumer on, when $numer is TRUE to allow
2159 ;; simplification of $%e to its numerical value.
2160 (simplifya (if $ratsimpexpons
($ratsimp
(caddr x
)) (caddr x
))
2164 (setq pot
(ratdisrep pot
))
2167 (cond ((member 'trunc
(car gr
))
2168 (return (srf (list '(mexpt) gr pot
))))
2170 (let ((varlist (caddar gr
)) (genvar (cadddr (car gr
))))
2171 (return (ratrep* (list '(mexpt) gr pot
)))))
2173 (setq gr
(ratdisrep gr
))
2175 ((or (setq mlpgr
(mxorlistp gr
))
2176 (setq mlppot
(mxorlistp pot
)))
2178 ((onep1 pot
) (go atgr
))
2179 ((or (zerop1 pot
) (onep1 gr
)) (go retno
))
2181 ;; This code tries to handle 0^a more complete.
2182 ;; If the sign of realpart(a) is not known return an unsimplified
2183 ;; expression. The handling of the flag *zexptsimp? is not changed.
2184 ;; Reverting the return of an unsimplified 0^a, because timesin
2185 ;; can not handle such expressions. (DK 02/2010)
2187 (cond ((or (member (setq z
($csign pot
)) '($neg $nz
))
2188 (and *zexptsimp?
(eq ($asksign pot
) '$neg
)))
2189 ;; A negative exponent. Maxima error.
2190 (cond ((not errorsw
) (merror (intl:gettext
"expt: undefined: 0 to a negative exponent.")))
2191 (t (throw 'errorsw t
))))
2192 ((and (member z
'($complex $imaginary
))
2193 ;; A complex exponent. Look at the sign of the realpart.
2194 (member (setq z
($sign
($realpart pot
)))
2196 (cond ((not errorsw
)
2197 (merror (intl:gettext
"expt: undefined: 0 to a complex exponent.")))
2198 (t (throw 'errorsw t
))))
2199 ((and *zexptsimp?
(eq ($asksign pot
) '$zero
))
2200 (cond ((not errorsw
)
2201 (merror (intl:gettext
"expt: undefined: 0^0")))
2202 (t (throw 'errorsw t
))))
2203 ((not (member z
'($pos $pz
)))
2204 ;; The sign of realpart(pot) is not known. We can not return
2205 ;; an unsimplified 0^a expression, because timesin can not
2206 ;; handle it. We return ZERO. That is the old behavior.
2207 ;; Look for the imaginary symbol to be consistent with
2209 (cond ((not (free pot
'$%i
))
2210 (cond ((not errorsw
)
2211 (merror (intl:gettext
"expt: undefined: 0 to a complex exponent.")))
2212 (t (throw 'errorsw t
))))
2214 ;; Return ZERO and not an unsimplified expression.
2215 (return (zerores gr pot
)))))
2216 (t (return (zerores gr pot
)))))
2220 (or (not (ratnump gr
)) (not (ratnump pot
))))
2221 (return (eqtest (exptrl gr pot
) check
)))
2222 ;; Check for numerical evaluation of the sqrt.
2223 ((and (alike1 pot
'((rat) 1 2))
2224 (or (setq res
(flonum-eval '%sqrt gr
))
2225 (and (not (member 'simp
(car x
)))
2226 (setq res
(big-float-eval '%sqrt gr
)))))
2229 (return (%itopot pot
)))
2230 ((and (realp gr
) (minusp gr
) (mevenp pot
))
2233 ((and (realp gr
) (minusp gr
) (moddp pot
))
2234 (return (mul2 -
1 (power (- gr
) pot
))))
2235 ((and (equal gr -
1) (maxima-integerp pot
) (mminusp pot
))
2236 (setq pot
(neg pot
))
2239 (maxima-integerp pot
)
2242 (integerp (cadr pot
))
2244 (maxima-integerp (caddr pot
)))
2245 (setq pot
(caddr pot
))
2247 ((atom gr
) (go atgr
))
2248 ((and (eq (caar gr
) 'mabs
)
2251 (or (and (eq $domain
'$real
) (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
))))
2252 (and (eq $domain
'$complex
) (apparently-real-to-judge-by-$csign-p
(cadr gr
)))))
2253 (return (power (cadr gr
) pot
)))
2254 ((and (eq (caar gr
) 'mabs
)
2257 (not (equal pot -
1))
2258 (or (and (eq $domain
'$real
) (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
))))
2259 (and (eq $domain
'$complex
) (apparently-real-to-judge-by-$csign-p
(cadr gr
)))))
2260 ;; abs(x)^(2*n+1) -> abs(x)*x^(2*n), n an integer number
2262 (return (mul (power (cadr gr
) (add pot -
1))
2264 (return (mul (power (cadr gr
) (add pot
1))
2266 ((eq (caar gr
) 'mequal
)
2267 (return (eqtest (list (ncons (caar gr
))
2268 (power (cadr gr
) pot
)
2269 (power (caddr gr
) pot
))
2271 ((symbolp pot
) (go opp
))
2272 ((eq (caar gr
) 'mexpt
) (go e1
))
2273 ((and (eq (caar gr
) '%sum
)
2278 (return (do ((i (1- pot
) (1- i
))
2279 (an gr
(simptimes (list '(mtimes) an gr
) 1 t
)))
2282 (return (eqtest (testt (tms gr pot nil
)) check
)))
2284 (return (eqtest (cond ((and (mplusp gr
)
2285 (not (or (> pot $expop
)
2286 (> (- pot
) $expon
))))
2287 (expandexpt gr pot
))
2288 (t (simplifya (tms gr pot nil
) t
)))
2292 (cond ((eq (caar gr
) 'mexpt
) (go e1
))
2293 ((eq (caar gr
) 'rat
)
2294 (return (mul2 (power (cadr gr
) pot
)
2295 (power (caddr gr
) (mul2 -
1 pot
)))))
2296 ((not (eq (caar gr
) 'mtimes
)) (go up
))
2297 ((or (eq $radexpand
'$all
) (and $radexpand
(simplexpon pot
)))
2300 ((and (or (not (numberp (cadr gr
)))
2301 (equal (cadr gr
) -
1))
2302 (equal -
1 ($num gr
)) ; only for -1
2303 ;; Do not simplify for a complex base.
2304 (not (member ($csign gr
) '($complex $imaginary
)))
2305 (and (eq $domain
'$real
) $radexpand
))
2306 ;; (-1/x)^a -> 1/(-x)^a for x negative
2307 ;; For all other cases (-1)^a/x^a
2308 (if (eq ($csign
(setq w
($denom gr
))) '$neg
)
2309 (return (inv (power (neg w
) pot
)))
2310 (return (div (power -
1 pot
)
2312 ((or (eq $domain
'$complex
) (not $radexpand
)) (go up
)))
2313 (return (do ((l (cdr gr
) (cdr l
)) (res (ncons 1)) (rad))
2315 (cond ((equal res
'(1))
2316 (eqtest (list '(mexpt) gr pot
) check
))
2318 (testt (cons '(mtimes simp
) res
)))
2320 (setq rad
(power* ; RADEXPAND=()?
2321 (cons '(mtimes) (nreverse rad
)) pot
))
2322 (cond ((not (onep1 rad
))
2324 (testt (tms rad
1 (cons '(mtimes) res
))))
2326 (setq *rulesw
* nil res
(cdr rad
))))))
2327 (eqtest (testt (cons '(mtimes) res
)) check
))))
2328 ;; Check with $csign to be more complete. This prevents wrong
2329 ;; simplifications like sqrt(-z^2)->%i*sqrt(z^2) for z complex.
2330 (setq z
($csign
(car l
)))
2331 (if (member z
'($complex $imaginary
))
2332 (setq z
'$pnz
)) ; if appears complex, unknown sign
2333 (setq w
(cond ((member z
'($neg $nz
))
2334 (setq rad
(cons -
1 rad
))
2338 ((alike1 w gr
) (return (list '(mexpt simp
) gr pot
)))
2339 ((member z
'($pn $pnz
))
2340 (setq rad
(cons w rad
)))
2342 (setq w
(testt (tms (simplifya (list '(mexpt) w pot
) t
)
2343 1 (cons '(mtimes) res
))))))
2344 (cond (*rulesw
* (setq *rulesw
* nil res
(cdr w
))))))
2347 (cond ((and (cdr res
) (onep1 (car res
)) (ratnump (cadr res
)))
2348 (setq res
(cdr res
))))
2349 (cond ((null (setq gr
(cdr gr
)))
2350 (return (eqtest (testt (cons '(mtimes) res
)) check
)))
2352 (setq y
(power (cadar gr
) (mult (caddar gr
) pot
))))
2354 (setq y
(%itopot pot
)))
2356 (setq y
(list '(mexpt) (car gr
) pot
)))
2357 (t (setq y
(list '(mexpt simp
) (car gr
) pot
))))
2358 (setq w
(testt (tms (simplifya y t
) 1 (cons '(mtimes) res
))))
2359 (cond (*rulesw
* (setq *rulesw
* nil res
(cdr w
))))
2363 (return (exptrl gr pot
))
2366 (cond ((zerop1 pot
) (go retno
))
2368 (let ((y (mget gr
'$numer
)))
2369 (if (and y
(floatp y
) (or $numer
(not (equal pot
1))))
2370 ;; A numeric constant like %e, %pi, ... and
2371 ;; exponent is a float or bigfloat value.
2372 (return (if (and (member gr
*builtin-numeric-constants
*)
2373 (equal pot bigfloatone
))
2374 ;; Return a bigfloat value.
2376 ;; Return a float value.
2378 ;; In all other cases exptrl simplifies accordingly.
2379 (return (exptrl gr pot
)))))
2381 ;; Numerically evaluate if the power is a flonum.
2383 (let ((val (flonum-eval '%exp pot
)))
2384 (when (float-inf-p val
)
2385 ;; needed for gcl and sbcl - (sometimes) no trap of overflow
2386 (error 'floating-point-overflow
))
2389 ;; Numerically evaluate if the power is a (complex)
2390 ;; big-float. (This is basically the guts of
2391 ;; big-float-eval, but we can't use big-float-eval.)
2392 (when (and (not (member 'simp
(car x
)))
2393 (complex-number-p pot
'bigfloat-or-number-p
))
2394 (let ((x ($realpart pot
))
2395 (y ($imagpart pot
)))
2396 (cond ((and ($bfloatp x
) (like 0 y
))
2397 (return ($bfloat
`((mexpt simp
) $%e
,pot
))))
2398 ((or ($bfloatp x
) ($bfloatp y
))
2399 (let ((z (add ($bfloat x
) (mul '$%i
($bfloat y
)))))
2400 (setq z
($rectform
`((mexpt simp
) $%e
,z
)))
2401 (return ($bfloat z
))))))))
2402 (cond ((and $logsimp
(among '%log pot
)) (return (%etolog pot
)))
2403 ((and $demoivre
(setq z
(demoivre pot
))) (return z
))
2407 ;; Exponent contains %i and %pi and %emode is TRUE:
2408 ;; Check simplification of exp(%i*%pi*p/q*x)
2409 (setq z
(%especial pot
)))
2411 (($taylorp
(third x
))
2412 ;; taylorize %e^taylor(...)
2413 (return ($taylor x
)))))
2415 (let ((y (mget gr
'$numer
)))
2416 ;; Check for a numeric constant.
2420 ;; The exponent is a bigfloat. Convert base to bigfloat.
2422 (member gr
*builtin-numeric-constants
*)
2423 (setq y
($bfloat gr
)))
2424 (and $numer
(integerp pot
)))
2425 (return (exptrl y pot
))))))
2428 (return (eqtest (list '(mexpt) gr pot
) check
))
2432 (cond ((mxorlistp1 gr
) (return (constmx (addk 1 pot
) gr
)))
2434 ((onep1 pot
) (return gr
))
2435 ((or $doallmxops $doscmxops $domxexpt
)
2436 (cond ((or (and mlpgr
2437 (or (not ($listp gr
)) $listarith
)
2438 (scalar-or-constant-p pot $assumescalar
))
2441 (or (not ($listp pot
)) $listarith
)
2442 (scalar-or-constant-p gr $assumescalar
)))
2443 (return (simplifya (outermap1 'mexpt gr pot
) t
)))
2445 ((and $domxmxops
(member pot
'(-1 -
1.0)))
2446 (return (simplifya (outermap1 'mexpt gr pot
) t
)))
2449 ;; At this point we have an expression: (z^a)^b with gr = z^a and pot = b
2450 (cond ((or (eq $radexpand
'$all
)
2451 ;; b is an integer or an odd rational
2453 (and (eq $domain
'$complex
)
2454 (not (member ($csign
(caddr gr
)) '($complex $imaginary
)))
2455 ;; z >= 0 and a not a complex
2456 (or (member ($csign
(cadr gr
)) '($pos $pz $zero
))
2458 (and (mnump (caddr gr
))
2459 (eq ($sign
(sub 1 (take '(mabs) (caddr gr
))))
2461 (and (eq $domain
'$real
)
2462 (member ($csign
(cadr gr
)) '($pos $pz $zero
)))
2463 ;; (1/z)^a -> 1/z^a when z a constant complex
2464 (and (eql (caddr gr
) -
1)
2466 (eq $domain
'$real
))
2467 (and (eq ($csign
(cadr gr
)) '$complex
)
2468 ($constantp
(cadr gr
)))))
2469 ;; This does (1/z)^a -> 1/z^a. This is in general wrong.
2470 ;; We switch this type of simplification on, when
2471 ;; $ratsimpexpons is T. E.g. radcan sets this flag to T.
2472 ;; radcan hangs for expressions like sqrt(1/(1+x)) without
2473 ;; this simplification.
2475 (equal (caddr gr
) -
1))
2478 (odnump (caddr gr
))))
2479 ;; Simplify (z^a)^b -> z^(a*b)
2480 (setq pot
(mul pot
(caddr gr
))
2482 ((and (eq $domain
'$real
)
2485 (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
)))
2486 (evnump (caddr gr
)))
2487 ;; Simplify (x^a)^b -> abs(x)^(a*b)
2488 (setq pot
(mul pot
(caddr gr
))
2489 gr
(radmabs (cadr gr
))))
2492 (mminusp (caddr gr
)))
2493 ;; Simplify (1/z^a)^b -> 1/(z^a)^b
2495 gr
(power (cadr gr
) (neg (caddr gr
)))))
2499 (defun apparently-complex-to-judge-by-$csign-p
(e)
2500 (let ((s ($csign e
)))
2501 (member s
'($complex $imaginary
))))
2503 (defun apparently-real-to-judge-by-$csign-p
(e)
2504 (let ((s ($csign e
)))
2505 (member s
'($pos $neg $zero $pn $pnz $pz $nz
))))
2507 ;; Basically computes log of m base b. Except if m is not a power
2508 ;; of b, we return nil. m is a positive integer and base an integer
2509 ;; not equal to +/-1.
2510 (defun exponent-of (m base
)
2511 ;; Just compute base^k until base^k >= m. Then check if they're equal.
2512 ;; If so, we have the exponent. Otherwise, give up.
2515 (multiple-value-bind (q r
)
2521 (if (zerop expo
) nil expo
)))
2523 (defun timesin (x y w
) ; Multiply X^W into Y
2524 (prog (fm temp z check u expo
)
2525 (if (mexptp x
) (setq check x
))
2527 ;; Prepare the factor x^w and initialize the work of timesin
2531 (setq temp
(cons '(mexpt) (if check
2532 (list (cadr x
) (mult (caddr x
) w
))
2534 (if (and (not *timesinp
*) (not (eq x
'$%i
)))
2535 (let ((*timesinp
* t
))
2536 (setq temp
(simplifya temp t
))))))
2537 (setq x
(if (mexptp temp
)
2543 ;; Go through the list of terms in fm and look what is to do.
2544 (cond ((null (cdr fm
))
2545 ;; The list of terms is empty. The loop is finished.
2547 ((or (and (mnump temp
)
2548 (not (or (integerp temp
)
2550 (and (integerp temp
)
2552 ;; Stop the loop for a float or bigfloat number, or number -1.
2555 (cond ((alike1 (car x
) (cadadr fm
))
2556 (cond ((zerop1 (setq w
(plsk (caddr (cadr fm
)) w
)))
2561 (rplacd fm
(cddr fm
))
2562 (cond ((mnump (setq x
(if (mnump (car x
))
2564 (power (car x
) w
))))
2565 (return (rplaca y
(timesk (car y
) x
))))
2570 x
(if (mexptp x
) (cdr x
) (list x
1)))
2574 ((maxima-constantp (car x
))
2577 (cond ((mtimesp (car x
))
2578 ;; A base which is a mtimes expression. Remove
2579 ;; the factor from the lists of products.
2580 (rplacd fm
(cddr fm
))
2581 ;; Multiply the factors of the base with
2582 ;; the list of all remaining products.
2584 (return (muln (nconc y
(cdar x
)) t
)))
2585 (t (return (rplaca (cdr fm
) (car x
))))))
2588 ;; At this place we have to add code for a rational number
2589 ;; as a factor to the list of products.
2591 (or (ratnump (car x
))
2592 (and (integerp (car x
))
2593 (not (onep (car x
))))))
2594 ;; Multiplying bas^k * num/den
2595 (let ((num (num1 (car x
)))
2596 (den (denom1 (car x
)))
2597 (bas (second (cadr fm
))))
2598 (cond ((and (integerp bas
)
2599 (not (eql 1 (abs bas
)))
2600 (setq expo
(exponent-of (abs num
) bas
)))
2601 ;; We have bas^m*bas^k = bas^(k+m).
2602 (setq temp
(power bas
2603 (add (third (cadr fm
)) expo
)))
2604 ;; Set fm to have 1/denom term.
2605 (setq x
(mul (car y
)
2609 ((and (integerp bas
)
2610 (not (eql 1 (abs bas
)))
2611 (setq expo
(exponent-of den bas
)))
2612 (setq expo
(- expo
))
2613 ;; We have bas^(-m)*bas^k = bas^(k-m).
2614 (setq temp
(power bas
2615 (add (third (cadr fm
)) expo
)))
2616 ;; Set fm to have the numerator term.
2617 (setq x
(mul (car y
)
2620 (exptrl bas
(- expo
)))))))
2622 ;; Next term in list of products.
2625 ;; Add in the bas^(k+m) term or bas^(k-m)
2626 (setf y
(rplaca y
1))
2627 (rplacd fm
(cddr fm
))
2628 (rplacd fm
(cons temp
(cdr fm
)))
2634 ((and (not (atom (car x
)))
2635 (eq (caar (car x
)) 'mabs
)
2637 (integerp (caddr (cadr fm
)))
2638 (< (caddr (cadr fm
)) -
1)
2639 (alike1 (cadr (car x
)) (cadr (cadr fm
)))
2640 (not (member ($csign
(cadr (car x
)))
2641 '($complex imaginary
))))
2642 ;; 1/x^n*abs(x) -> 1/(x^(n-2)*abs(x)), where n an integer
2643 ;; Replace 1/x^n -> 1/x^(n-2)
2644 (setq temp
(power (cadr (cadr fm
))
2645 (add (caddr (cadr fm
)) 2)))
2646 (rplacd fm
(cddr fm
))
2647 (if (not (equal temp
1))
2648 (rplacd fm
(cons temp
(cdr fm
))))
2649 ;; Multiply factor 1/abs(x) into list of products.
2650 (setq x
(list (car x
) -
1))
2651 (setq temp
(power (car x
) (cadr x
)))
2655 ((and (not (atom (car x
)))
2656 (eq (caar (car x
)) 'mabs
)
2658 (integerp (caddr (cadr fm
)))
2659 (> (caddr (cadr fm
)) 1)
2660 (alike1 (cadr (car x
)) (cadr (cadr fm
)))
2661 (not (member ($csign
(cadr (car x
)))
2662 '($complex imaginary
))))
2663 ;; x^n/abs(x) -> x^(n-2)*abs(x), where n an integer.
2664 ;; Replace x^n -> x^(n-2)
2665 (setq temp
(power (cadr (cadr fm
))
2666 (add (caddr (cadr fm
)) -
2)))
2667 (rplacd fm
(cddr fm
))
2668 (if (not (equal temp
1))
2669 (rplacd fm
(cons temp
(cdr fm
))))
2670 ;; Multiply factor abs(x) into list of products.
2671 (setq x
(list (car x
) 1))
2672 (setq temp
(power (car x
) (cadr x
)))
2676 ((and (not (atom (cadr fm
)))
2677 (not (atom (cadr (cadr fm
))))
2678 (eq (caaadr (cadr fm
)) 'mabs
)
2679 (equal (caddr (cadr fm
)) -
1)
2682 (alike1 (cadadr (cadr fm
)) (car x
))
2683 (not (member ($csign
(cadadr (cadr fm
)))
2684 '($complex imaginary
))))
2685 ;; 1/abs(x)*x^n -> x^(n-2)*abs(x), where n an integer.
2686 ;; Replace 1/abs(x) -> abs(x)
2687 (setq temp
(cadr (cadr fm
)))
2688 (rplacd fm
(cddr fm
))
2689 (rplacd fm
(cons temp
(cdr fm
)))
2690 ;; Multiply factor x^(n-2) into list of products.
2691 (setq x
(list (car x
) (add (cadr x
) -
2)))
2692 (setq temp
(power (car x
) (cadr x
)))
2696 ((or (maxima-constantp (car x
))
2697 (maxima-constantp (cadadr fm
)))
2698 (if (great temp
(cadr fm
))
2700 ((great (car x
) (cadadr fm
))
2703 ((alike1 (car x
) (cadr fm
))
2706 ;; When a number goto start and look in the next term.
2710 ((and (not (atom (cadr fm
)))
2711 (eq (caar (cadr fm
)) 'mabs
)
2714 (alike1 (cadr (cadr fm
)) (car x
))
2715 (not (member ($csign
(cadr (cadr fm
)))
2716 '($complex imaginary
))))
2717 ;; abs(x)/x^n -> 1/(x^(n-2)*abs(x)), where n an integer.
2718 ;; Replace abs(x) -> 1/abs(x).
2719 (setq temp
(power (cadr fm
) -
1))
2720 (rplacd fm
(cddr fm
))
2721 (rplacd fm
(cons temp
(cdr fm
)))
2722 ;; Multiply factor x^(-n+2) into list of products.
2723 (setq x
(list (car x
) (add (cadr x
) 2)))
2724 (setq temp
(power (car x
) (cadr x
)))
2728 ((maxima-constantp (car x
))
2729 (when (great temp
(cadr fm
))
2731 ((great (car x
) (cadr fm
))
2735 ;; Multiply a number into the list of products.
2736 (return (rplaca y
(timesk (car y
) temp
))))
2737 ((and (eq (car x
) '$%i
)
2740 ((and (eq (car x
) '$%e
)
2743 (return (rplaca y
(timesk (car y
) (exp (float w
))))))
2745 (not (constant (car x
))))
2747 ;; At this point we will insert a mexpt expression,
2748 ;; but first we look at the car of the list of products and
2749 ;; modify the expression if we found a rational number.
2751 (not (onep1 (car y
)))
2752 (or (integerp (car y
))
2754 ;; Multiplying bas^k * num/den.
2755 (let ((num (num1 (car y
)))
2756 (den (denom1 (car y
)))
2758 (cond ((and (integerp bas
)
2759 (not (eql 1 (abs bas
)))
2760 (setq expo
(exponent-of (abs num
) bas
)))
2761 ;; We have bas^m*bas^k.
2762 (setq temp
(power bas
(add (cadr x
) expo
)))
2763 ;; Set fm to have 1/denom term.
2764 (setq x
(div (div num
(exptrl bas expo
)) den
)))
2765 ((and (integerp bas
)
2766 (not (eql 1 (abs bas
)))
2767 (setq expo
(exponent-of den bas
)))
2768 (setq expo
(- expo
))
2769 ;; We have bas^(-m)*bas^k.
2770 (setq temp
(power bas
(add (cadr x
) expo
)))
2771 ;; Set fm to have the numerator term.
2772 (setq x
(div num
(div den
(exptrl bas
(- expo
))))))
2774 ;; The rational doesn't contain any (simple) powers of
2775 ;; the exponential term. We're done.
2776 (return (cdr (rplacd fm
(cons temp
(cdr fm
)))))))
2777 ;; Add in the a^(m+k) or a^(k-m) term.
2778 (setf y
(rplaca y
1))
2779 (rplacd fm
(cons temp
(cdr fm
)))
2785 ((and (maxima-constantp (car x
))
2786 (do ((l (cdr fm
) (cdr l
)))
2788 (when (and (mexptp (cadr l
))
2789 (alike1 (car x
) (cadadr l
)))
2793 ((or (and (mnump (car x
))
2795 (and (eq (car x
) '$%e
)
2799 (setq u
(%especial w
))))
2801 ((alike (cdr check
) x
)
2804 (exptrl (car x
) w
))))
2806 (return (rplaca y
(timesk (car y
) x
))))
2810 (return (cdr (rplacd fm
(cons x
(cdr fm
))))))
2820 (setq temp
(list '(mexpt) (car x
) w
))
2821 (setq temp
(eqtest temp
(or check
'((foo)))))
2822 (return (cdr (rplacd fm
(cons temp
(cdr fm
)))))))
2824 (return (cdr (rplacd fm
(cons (car x
) (cdr fm
)))))
2829 (cond ((and (eq (car x
) '$%i
) (equal w
1))
2830 (rplacd fm
(cddr fm
))
2831 (return (rplaca y
(timesk -
1 (car y
)))))
2832 ((zerop1 (setq w
(plsk 1 w
)))
2834 ((and (mnump (car x
)) (mnump w
))
2835 (return (rplaca (cdr fm
) (exptrl (car x
) w
))))
2836 ((maxima-constantp (car x
))
2839 (setq z
(list '(mexpt) (car x
) w
))
2840 (cond ((alike1 (setq x
(simplifya z t
)) z
)
2841 (return (rplaca (cdr fm
) x
)))
2843 (rplacd fm
(cddr fm
))
2845 (return (muln (cons x y
) t
))))
2847 (rplacd fm
(cddr fm
))
2848 (setq x
(car x
) check nil
)
2851 (setq z
(tms x
1 (setq temp
(cons '(mtimes) y
))))
2852 (return (cond ((eq z temp
)
2855 (setq *rulesw
* t
) z
)))
2857 (return (rplacd fm
(cddr fm
)))
2859 (if (minusp (setq w
(rem w
4)))
2861 (return (cond ((zerop w
)
2864 (rplaca y
(timesk -
1 (car y
))))
2866 (rplaca y
(timesk -
1 (car y
)))
2867 (rplacd fm
(cons '$%i
(cdr fm
))))
2869 (rplacd fm
(cons '$%i
(cdr fm
))))))))
2871 (defun simpmatrix (x vestigial z
)
2872 (declare (ignore vestigial
))
2873 (if (and (null (cddr x
))
2875 (or (eq $scalarmatrixp
'$all
) (member 'mult
(cdar x
)))
2876 ($listp
(cadr x
)) (cdadr x
) (null (cddadr x
)))
2877 (simplifya (cadadr x
) z
)
2878 (let ((badp (dolist (row (cdr x
)) (if (not ($listp row
)) (return t
))))
2879 (args (simpmap (cdr x
) z
)))
2880 (if (and args
(not badp
)) (matcheck args
))
2881 (cons (if badp
'(%matrix simp
) '($matrix simp
)) args
))))
2883 (defun %itopot
(pot)
2885 (let ((i (boole boole-and pot
3)))
2889 (t (list '(mtimes simp
) -
1 '$%i
))))
2890 (power -
1 (mul2 pot
'((rat simp
) 1 2)))))
2893 (cond ((eq (caar pot
) '%log
) (simplifya (cadr pot
) nil
))
2894 ((and (eq (caar pot
) 'mtimes
)
2895 (or (maxima-integerp (cadr pot
))
2896 (and $%e_to_numlog
($numberp
(cadr pot
))))
2897 (not (atom (caddr pot
))) (eq (caar (caddr pot
)) '%log
)
2899 (power (cadr (caddr pot
)) (cadr pot
)))))
2903 loop
(cond ((null pot
)
2904 (cond (a (setq a
(cons '(mtimes) a
))))
2905 (cond (c (setq c
(list '(mexpt simp
) '$%e
(addn c nil
)))))
2906 (return (cond ((null c
) (simptimes a
1 nil
))
2908 (t (simptimes (append a
(list c
)) 1 nil
)))))
2909 ((and (among '%log
(car pot
)) (setq b
(mnlogp (car pot
))))
2910 (setq a
(cons b a
)))
2911 (t (setq c
(cons (car pot
) c
))))
2912 (setq pot
(cdr pot
))
2915 (defun %etolog
(pot) (cond ((mnlogp pot
))
2916 ((eq (caar pot
) 'mplus
) (mnlog (cdr pot
)))
2917 (t (list '(mexpt simp
) '$%e pot
))))
2919 (defun zerores (r1 r2
)
2920 (cond ((or ($bfloatp r1
) ($bfloatp r2
)) bigfloatzero
)
2921 ((or (floatp r1
) (floatp r2
)) 0.0)
2924 (defmfun $orderlessp
(a b
)
2925 (setq a
($totaldisrep
(specrepcheck a
))
2926 b
($totaldisrep
(specrepcheck b
)))
2927 (and (not (alike1 a b
)) (great b a
) t
))
2929 (defmfun $ordergreatp
(a b
)
2930 (setq a
($totaldisrep
(specrepcheck a
))
2931 b
($totaldisrep
(specrepcheck b
)))
2932 (and (not (alike1 a b
)) (great a b
) t
))
2934 ;; Test function to order a and b by magnitude. If it is not possible to
2935 ;; order a and b by magnitude they are ordered by great. This function
2936 ;; can be used by sort, e.g. sort([3,1,7,x,sin(1),minf],ordermagnitudep)
2937 (defmfun $ordermagnitudep
(a b
)
2939 (setq a
($totaldisrep
(specrepcheck a
))
2940 b
($totaldisrep
(specrepcheck b
)))
2941 (cond ((and (or (constp a
) (member a
'($inf $minf
)))
2942 (or (constp b
) (member b
'($inf $minf
)))
2943 (member (setq sgn
($csign
(sub b a
))) '($pos $neg $zero
)))
2944 (cond ((eq sgn
'$pos
) t
)
2945 ((eq sgn
'$zero
) (and (not (alike1 a b
)) (great b a
)))
2947 ((or (constp a
) (member a
'($inf $minf
))) t
)
2948 ((or (constp b
) (member b
'($inf $minf
))) nil
)
2949 (t (and (not (alike1 a b
)) (great b a
))))))
2951 (defun evnump (n) (or (even n
) (and (ratnump n
) (even (cadr n
)))))
2952 (defun odnump (n) (or (and (integerp n
) (oddp n
))
2953 (and (ratnump n
) (oddp (cadr n
)))))
2955 (defun simplexpon (e)
2956 (or (maxima-integerp e
)
2957 (and (eq $domain
'$real
) (ratnump e
) (oddp (caddr e
)))))
2959 ;; This function is not called in Maxima core or share code
2960 ;; and can be cut out.
2962 (and (free p
'$%i
) (member ($sign p
) '($pos $pz $zero
))))
2965 (if (and limitp
(free e
'$%i
)) (asksign-p-or-n e
))
2966 (simplifya (list '(mabs) e
) t
))
2968 (defun simpmqapply (exp y z
)
2969 (let ((simpfun (and (not (atom (cadr exp
))) (safe-get (caaadr exp
) 'specsimp
))) u
)
2971 (funcall simpfun exp y z
)
2972 (progn (setq u
(simpargs exp z
))
2973 (if (symbolp (cadr u
))
2974 (simplifya (cons (cons (cadr u
) (cdar u
)) (cddr u
)) z
)
2977 ;; TRUE, if the symbol e is declared to be $complex or $imaginary.
2978 (defun decl-complexp (e)
2980 (kindp e
'$complex
)))
2982 ;; TRUE, if the symbol e is declared to be $real, $rational, $irrational
2984 (defun decl-realp (e)
2986 (or (kindp e
'$real
)
2987 (kindp e
'$rational
)
2988 (kindp e
'$irrational
)
2989 (kindp e
'$integer
))))
2991 ;; WARNING: Exercise extreme caution when modifying this function!
2993 ;; Richard Fateman and Stavros Macrakis both say that changing the
2994 ;; actual ordering relations (as opposed to making them faster to
2995 ;; determine) could have very subtle and wide-ranging effects. Also,
2996 ;; the simplifier spends the vast majority of its time here, so be
2997 ;; very careful about changes that may drastically slow down the
3005 (cond ((zerop y
) (floatp x
)) (t (plusp y
))))))
3007 (cond ((constant y
) (alphalessp y x
)) (t (numberp y
))))
3009 (cond ((mget y
'$scalar
) (alphalessp y x
))
3010 (t (maxima-constantp y
))))
3012 (cond ((mget y
'$mainvar
) (alphalessp y x
)) (t t
)))
3013 (t (or (maxima-constantp y
) (mget y
'$scalar
)
3014 (and (not (mget y
'$mainvar
)) (not (null (alphalessp y x
))))))))
3015 (t (not (ordfna y x
)))))
3016 ((atom y
) (ordfna x y
))
3018 (cond ((eq (caar y
) 'rat
)
3019 (> (* (caddr y
) (cadr x
)) (* (caddr x
) (cadr y
))))))
3020 ((eq (caar y
) 'rat
))
3021 ((or (member (caar x
) '(mtimes mplus mexpt %del
))
3022 (member (caar y
) '(mtimes mplus mexpt %del
)))
3024 ((and (eq (caar x
) 'bigfloat
) (eq (caar y
) 'bigfloat
)) (mgrp x y
))
3025 ((or (eq (caar x
) 'mrat
) (eq (caar y
) 'mrat
))
3026 (error "GREAT: internal error: unexpected MRAT argument"))
3027 (t (do ((x1 (margs x
) (cdr x1
)) (y1 (margs y
) (cdr y1
))) (())
3029 (return (cond (y1 nil
)
3030 ((not (alike1 (mop x
) (mop y
)))
3031 (great (mop x
) (mop y
)))
3032 ((member 'array
(cdar x
)) t
))))
3033 ((null y1
) (return t
))
3034 ((not (alike1 (car x1
) (car y1
)))
3035 (return (great (car x1
) (car y1
)))))))))
3037 ;; Trivial function used only in ALIKE1.
3038 ;; Should be defined as an open-codable subr.
3040 (defmacro memqarr
(l)
3041 `(if (member 'array
,l
) t
))
3043 ;; Compares two Macsyma expressions ignoring SIMP flags and all other
3044 ;; items in the header except for the ARRAY flag.
3047 ;; Clauses are ordered based on frequency of the case
3048 ;; cons, integer, and symbol are very common
3049 ;; everything else is rare
3053 (not (atom (car x
)))
3054 (not (atom (car y
)))
3055 (eq (caar x
) (caar y
)))
3057 ((eq (caar x
) 'mrat
) (like x y
))
3058 ((eq (caar x
) 'mpois
) (equal (cdr x
) (cdr y
)))
3059 ((eq (memqarr (cdar x
)) (memqarr (cdar y
)))
3060 (alike (cdr x
) (cdr y
)))
3062 ;; (foo) and (foo) test non-alike because the car's aren't standard
3064 ((or (symbolp x
) (symbolp y
)) nil
)
3065 ((integerp x
) (and (integerp y
) (= x y
)))
3066 ;; uncommon cases from here down
3067 ((floatp x
) (and (floatp y
) (= x y
)))
3068 ((stringp x
) (and (stringp y
) (string= x y
)))
3069 ((vectorp x
) (and (vectorp y
) (lisp-vector-alike1 x y
)))
3070 ((arrayp x
) (and (arrayp y
) (lisp-array-alike1 x y
)))
3074 (defun lisp-vector-alike1 (x y
)
3075 (let ((lx (length x
)))
3076 (when (eql lx
(length y
))
3077 (lisp-array-elements-alike1 x y lx
))))
3079 (defun lisp-array-alike1 (x y
)
3080 (when (equal (array-dimensions x
) (array-dimensions y
))
3081 (lisp-array-elements-alike1 x y
(array-total-size x
))))
3083 (defun lisp-array-elements-alike1 (x y n
)
3085 (unless (alike1 (row-major-aref x i
) (row-major-aref y i
))
3086 (return-from lisp-array-elements-alike1 nil
))))
3088 ;; Not sure if we want to enable comparison of maxima arrays.
3089 ;; Aside from that, add2lnc calls alike1 (via memalike) and that causes trouble.
3090 ;; Possible code for ALIKE1 to enable such comparisons should we choose to do so:
3092 ;; ((maxima-declared-arrayp x)
3093 ;; (and (maxima-declared-arrayp y) (maxima-declared-array-alike1 x y)))
3094 ;; ((maxima-undeclared-arrayp x)
3095 ;; (and (maxima-undeclared-arrayp y) (maxima-undeclared-array-alike1 x y)))
3097 (defun maxima-declared-array-alike1 (x y
)
3098 (lisp-array-alike1 (get (mget x
'array
) 'array
) (get (mget y
'array
) 'array
)))
3100 (defun maxima-undeclared-array-alike1 (x y
)
3102 (alike1 (mfuncall '$arrayinfo x
) (mfuncall '$arrayinfo y
))
3103 (alike1 ($listarray x
) ($listarray y
))))
3105 ;; Maps ALIKE1 down two lists.
3108 (do ((x x
(cdr x
)) (y y
(cdr y
))) ((atom x
) (equal x y
))
3109 (cond ((or (atom y
) (not (alike1 (car x
) (car y
))))
3112 (defun ordfna (e a
) ; A is an atom
3114 (or (not (eq (caar e
) 'rat
))
3115 (> (cadr e
) (* (caddr e
) a
))))
3117 (not (member (caar e
) '(mplus mtimes mexpt
))))
3118 (not (member (caar e
) '(rat bigfloat
))))
3119 ((eq (caar e
) 'mrat
)) ;; all MRATs succeed all atoms
3120 ((null (margs e
)) nil
)
3121 ((eq (caar e
) 'mexpt
)
3122 (cond ((and (maxima-constantp (cadr e
))
3123 (or (not (constant a
)) (not (maxima-constantp (caddr e
)))))
3124 (or (not (free (caddr e
) a
)) (great (caddr e
) a
)))
3125 ((eq (cadr e
) a
) (great (caddr e
) 1))
3126 (t (great (cadr e
) a
))))
3127 ((member (caar e
) '(mplus mtimes
))
3128 (let ((u (car (last e
))))
3129 (cond ((eq u a
) (not (ordhack e
))) (t (great u a
)))))
3130 ((eq (caar e
) '%del
))
3131 ((prog2 (setq e
(car (margs e
))) ; use first arg of e
3132 (and (not (atom e
)) (member (caar e
) '(mplus mtimes
))))
3133 (let ((u (car (last e
)))) ; and compare using
3134 (cond ((eq u a
) (not (ordhack e
))) ; same procedure as above
3139 ;; compare lists a and b elementwise from back to front
3140 (defun ordlist (a b cx cy
)
3142 (setq l1
(length a
) l2
(length b
))
3143 loop
(cond ((= l1
0)
3144 (return (cond ((= l2
0) (eq cx
'mplus
))
3145 ((and (eq cx cy
) (= l2
1))
3146 (great (cond ((eq cx
'mplus
) 0) (t 1)) (car b
))))))
3147 ((= l2
0) (return (not (ordlist b a cy cx
)))))
3148 (setq c
(nthelem l1 a
) d
(nthelem l2 b
))
3149 (cond ((not (alike1 c d
)) (return (great c d
))))
3150 (setq l1
(1- l1
) l2
(1- l2
))
3153 (defun term-list (x)
3158 (defun factor-list (x)
3163 ;; one of the exprs x or y should be one of:
3164 ;; %del, mexpt, mplus, mtimes
3166 (let ((cx (caar x
)) (cy (caar y
)))
3167 (cond ((eq cx
'%del
) (if (eq cy
'%del
) (great (cadr x
) (cadr y
)) t
))
3169 ((or (eq cx
'mtimes
) (eq cy
'mtimes
))
3170 (ordlist (factor-list x
) (factor-list y
) 'mtimes
'mtimes
))
3171 ((or (eq cx
'mplus
) (eq cy
'mplus
))
3172 (ordlist (term-list x
) (term-list y
) 'mplus
'mplus
))
3173 ((eq cx
'mexpt
) (ordmexpt x y
))
3174 ((eq cy
'mexpt
) (not (ordmexpt y x
))))))
3177 (if (and (cddr x
) (null (cdddr x
)))
3178 (great (if (eq (caar x
) 'mplus
) 0 1) (cadr x
))))
3180 (defun ordmexpt (x y
)
3181 (cond ((eq (caar y
) 'mexpt
)
3182 (cond ((alike1 (cadr x
) (cadr y
)) (great (caddr x
) (caddr y
)))
3183 ((maxima-constantp (cadr x
))
3184 (if (maxima-constantp (cadr y
))
3185 (if (or (alike1 (caddr x
) (caddr y
))
3186 (and (mnump (caddr x
)) (mnump (caddr y
))))
3187 (great (cadr x
) (cadr y
))
3188 (great (caddr x
) (caddr y
)))
3189 (great x
(cadr y
))))
3190 ((maxima-constantp (cadr y
)) (great (cadr x
) y
))
3192 (great (cadr x
) (if (mnump (caddr y
)) (cadr y
) y
)))
3193 ((mnump (caddr y
)) (great x
(cadr y
)))
3194 (t (let ((x1 (simpln1 x
)) (y1 (simpln1 y
)))
3195 (if (alike1 x1 y1
) (great (cadr x
) (cadr y
))
3197 ((alike1 (cadr x
) y
) (great (caddr x
) 1))
3198 ((mnump (caddr x
)) (great (cadr x
) y
))
3199 (t (great (simpln1 x
)
3202 (defmfun $multthru
(e1 &optional e2
)
3204 (cond (e2 ;called with two args
3205 (setq arg1
(specrepcheck e1
)
3206 arg2
(specrepcheck e2
))
3207 (cond ((or (atom arg2
)
3208 (not (member (caar arg2
) '(mplus mequal
))))
3210 ((eq (caar arg2
) 'mequal
)
3211 (list (car arg2
) ($multthru arg1
(cadr arg2
))
3212 ($multthru arg1
(caddr arg2
))))
3213 (t (expandterms arg1
(cdr arg2
)))))
3214 (t ;called with only one arg
3216 (setq arg1
(setq arg2
(specrepcheck e1
)))
3217 (cond ((atom arg1
) (return arg1
))
3218 ((eq (caar arg1
) 'mnctimes
)
3219 (setq arg1
(cdr arg1
)) (go nct
))
3220 ((not (eq (caar arg1
) 'mtimes
)) (return arg1
)))
3221 (setq arg1
(reverse (cdr arg1
)))
3222 times
(when (mplusp (car arg1
))
3223 (setq l1
(nconc l1
(cdr arg1
)))
3224 (return (expandterms (muln l1 t
) (cdar arg1
))))
3225 (setq l1
(cons (car arg1
) l1
))
3226 (setq arg1
(cdr arg1
))
3227 (if (null arg1
) (return arg2
))
3229 nct
(when (mplusp (car arg1
))
3230 (setq l1
(nreverse l1
))
3231 (return (addn (mapcar
3235 (append l1
(ncons u
) (cdr arg1
)))
3239 (setq l1
(cons (car arg1
) l1
))
3240 (setq arg1
(cdr arg1
))
3241 (if (null arg1
) (return arg2
))
3244 ;; EXPANDEXPT computes the expansion of (x1 + x2 + ... + xm)^n
3245 ;; taking a sum and integer power as arguments.
3246 ;; Its theory is to recurse down the binomial expansion of
3247 ;; (x1 + (x2 + x3 + ... + xm))^n using the Binomial Expansion
3248 ;; Thus it does a sigma:
3252 ;; \ / n \ k (n - k)
3253 ;; > | | x1 (x2 + x3 + ... + xm)
3258 ;; The function EXPONENTIATE-SUM does this and recurses through the second
3259 ;; sum raised to a power. It takes a list of terms and a positive integer
3260 ;; power as arguments.
3263 (defun expandexpt (sum power
)
3264 (declare (fixnum power
))
3265 (let ((expansion (exponentiate-sum (cdr sum
) (abs power
))))
3266 (cond ((plusp power
) expansion
)
3267 (t (inv expansion
)))))
3269 (defun exponentiate-sum (terms rpower
)
3270 (declare (fixnum rpower
))
3271 (cond ((= rpower
0) 1)
3272 ((null (cdr terms
)) (power (car terms
) rpower
))
3273 ((= rpower
1) (cons '(mplus simp
) terms
))
3274 (t (do ((i 0 (1+ i
))
3275 (result 0 (add2 result
3276 (muln (list (combination rpower i
)
3277 (exponentiate-sum (cdr terms
)
3279 (power (car terms
) i
)) t
))))
3280 ((> i rpower
) result
)
3281 (declare (fixnum i
))))))
3283 ;; Computes the combination of n elements taken m at a time by the formula
3285 ;; (n * (n-1) * ... * (n - m + 1)) / m! =
3286 ;; (n / 1) * ((n - 1) / 2) * ... * ((n - m + 1) / m)
3288 ;; Checks for the case when m is greater than n/2 and translates
3289 ;; to an equivalent expression.
3291 (defun combination (n m
)
3292 (declare (fixnum n m
))
3293 (cond ((> m
(truncate n
2))
3294 (combination n
(- n m
)))
3296 (do ((result 1 (truncate (* result n1
) m1
))
3300 (declare (fixnum n1 m1
))))))
3302 (defun expandsums (a b
)
3304 (setq a
(fixexpand a
) b
(cdr b
))
3306 (when (null a
) (return c
))
3307 (setq c
(cons (expandterms (car a
) b
) c
))
3312 (defun expandterms (a b
)
3315 (when (null b
) (return c
))
3316 (setq c
(cons (mul2 a
(car b
)) c
))
3321 (defun expandtimes (a)
3322 (let (simp-prods simp-negprods simp-sums simp-negsums expsums expnegsums
)
3329 (setq simp-prods
(nreverse simp-prods
)
3330 simp-negprods
(nreverse simp-negprods
)
3331 simp-sums
(nreverse simp-sums
)
3332 simp-negsums
(nreverse simp-negsums
))
3335 (push (car l
) simp-prods
))
3336 ((eq (caaar l
) 'rat
)
3337 (unless (equal (cadar l
) 1)
3338 (push (cadar l
) simp-prods
))
3339 (push (caddar l
) simp-negprods
))
3340 ((eq (caaar l
) 'mplus
)
3341 (push (car l
) simp-sums
))
3342 ((and (eq (caaar l
) 'mexpt
)
3343 (equal (caddar l
) -
1) (mplusp (cadar l
)))
3344 (push (cadar l
) simp-negsums
))
3345 ((and (eq (caaar l
) 'mexpt
)
3346 (let ((*expandp
* t
))
3347 (mminusp (caddar l
))))
3348 (push (if (equal (caddar l
) -
1)
3350 (list (caar l
) (cadar l
) (neg (caddar l
))))
3353 (push (car l
) simp-prods
)))
3356 (setq simp-prods
(cond ((null simp-prods
) 1)
3357 ((null (cdr simp-prods
)) (car simp-prods
))
3358 (t (cons '(mtimes simp
) simp-prods
))))
3359 (setq simp-negprods
(cond ((null simp-negprods
) 1)
3360 ((null (cdr simp-negprods
)) (car simp-negprods
))
3361 (t (cons '(mtimes simp
) simp-negprods
))))
3364 (cond ((null simp-sums
) (go down
))
3365 (t (setq expsums
(car simp-sums
))
3367 (setq expsums
(expandsums expsums c
)))
3369 (setq simp-prods
(cond ((equal simp-prods
1) expsums
)
3370 (t (expandterms simp-prods
(fixexpand expsums
)))))
3371 down
(cond ((null simp-negsums
)
3372 (cond ((equal 1 simp-negprods
) (return simp-prods
))
3373 ((mplusp simp-prods
)
3374 (return (expandterms (power simp-negprods -
1) (cdr simp-prods
))))
3375 (t (return (let ((*expandflag
* t
))
3376 (mul2 simp-prods
(power simp-negprods -
1)))))))
3378 (setq expnegsums
(car simp-negsums
))
3380 (setq expnegsums
(expandsums expnegsums c
)))
3381 (cdr simp-negsums
))))
3382 (setq expnegsums
(expandterms simp-negprods
(fixexpand expnegsums
)))
3383 (return (if (mplusp simp-prods
)
3384 (expandterms (inv expnegsums
) (cdr simp-prods
))
3385 (let ((*expandflag
* t
))
3386 (mul2 simp-prods
(inv expnegsums
))))))))))
3388 (defun expand1 (exp $expop $expon
)
3389 (unless (and (integerp $expop
) (> $expop -
1))
3390 (merror (intl:gettext
"expand: expop must be a nonnegative integer; found: ~M") $expop
))
3391 (unless (and (integerp $expon
) (> $expon -
1))
3392 (merror (intl:gettext
"expand: expon must be a nonnegative integer; found: ~M") $expon
))
3393 (resimplify (specrepcheck exp
)))
3395 (defmfun $expand
(exp &optional
(expop $maxposex
) (expon $maxnegex
))
3396 (expand1 exp expop expon
))
3398 (defun fixexpand (a)
3399 (if (not (mplusp a
))
3403 (defun nrthk (in nth-root
)
3405 ((nrthk1 (in nth-root
)
3406 ;; Computes the NTH-ROOT of -IN
3408 (mul2 (nrthk2 in nth-root
) (nrthk -
1 nth-root
))
3409 (nrthk2 (mul2* -
1 in
) nth-root
)))
3411 (nrthk2 (in nth-root
)
3412 ;; Computes the NTH-ROOT of IN
3413 (power* in
(list '(rat) 1 nth-root
))))
3417 (cond ((equal nth-root
2)
3419 ((eq $domain
'$real
)
3421 (nrthk2 -
1 nth-root
)
3425 (power* '$%e
(list '(mtimes) (list '(rat) 1 nth-root
) '$%pi
'$%i
))))
3427 (nrthk2 -
1 nth-root
))))
3428 ((or (and wflag
(eq ($asksign in
) '$neg
))
3429 (and (mnump in
) (equal ($sign in
) '$neg
)))
3430 (nrthk1 (mul2* -
1 in
) nth-root
))
3432 (nrthk2 in nth-root
)))))
3434 (defun simpnrt (x nth-root
)
3435 ;; Computes the NTH-ROOT of X
3436 (prog (simp-in simp-out varlist genvar $factorflag $dontfactor
)
3439 (do ((x x
(cddr x
)) (y))
3441 (cond ((not (equal 1 (setq y
(gcd (cadr x
) nth-root
))))
3442 (push (simpnrt (list '(mexpt) (car x
) (quotient (cadr x
) y
))
3443 (quotient nth-root y
))
3445 ((and (equal (cadr x
) 1) (integerp (car x
)) (plusp (car x
))
3446 (setq y
(pnthrootp (car x
) nth-root
)))
3449 (unless (> nth-root
(abs (cadr x
)))
3450 (push (list '(mexpt) (car x
) (quotient (cadr x
) nth-root
)) simp-out
))
3451 (push (list '(mexpt) (car x
) (rem (cadr x
) nth-root
)) simp-in
))))))
3452 (setq $factorflag t
)
3454 (setq x
(ratrep* x
))
3455 (when (equal (cadr x
) 0) (return 0))
3456 (setq x
(ratfact (cdr x
) 'psqfr
))
3457 (simpnrt1 (mapcar #'pdis x
))
3458 (setq simp-out
(if simp-out
(muln simp-out nil
) 1))
3459 (setq simp-in
(cond (simp-in
3460 (setq simp-in
(muln simp-in nil
))
3461 (nrthk simp-in nth-root
))
3463 (return (let (($%emode t
))
3464 (simplifya (list '(mtimes) simp-in simp-out
)
3465 (not (or (atom simp-in
)
3466 (atom (cadr simp-in
))
3467 (member (caaadr simp-in
) '(mplus mtimes rat
))))))))))
3469 ;; The following was formerly in SININT. This code was placed here because
3470 ;; SININT is now an out-of-core file on MC, and this code is needed in-core
3471 ;; because of the various calls to it. - BMT & JPG
3473 (defmfun $integrate
(expr x
&optional lo hi
)
3474 (declare (special *in-risch-p
* context
))
3477 (with-new-context (context)
3478 (if (member '%risch
*nounl
*)
3480 ;; Give up; we're being called from RISCHINT by some path.
3481 (list '(%integrate
) expr x
)
3484 ($defint expr x lo hi
))))
3486 (defun ratp (a ratp-var
)
3488 ((member (caar a
) '(mplus mtimes
))
3489 (do ((l (cdr a
) (cdr l
))) ((null l
) t
)
3490 (or (ratp (car l
) ratp-var
) (return nil
))))
3491 ((eq (caar a
) 'mexpt
)
3492 (if (free (cadr a
) ratp-var
)
3493 (free (caddr a
) ratp-var
)
3494 (and (integerp (caddr a
)) (ratp (cadr a
) ratp-var
))))
3495 (t (free a ratp-var
))))
3497 (defun ratnumerator (r)
3499 ((atom (cdr r
)) (car r
))
3500 ((numberp (cadr r
)) r
)
3503 (defun ratdenominator (r)
3505 ((atom (cdr r
)) (cdr r
))
3506 ((numberp (cadr r
)) 1)
3509 ;; (BPROG U V) appears to return A and B (as ((A1 . A2) B1 . B2) with A = A1/A2, B = B1/B2)
3510 ;; such that B/U + A/V = 1/(U*V), where U, V are polynomials represented as a list of
3511 ;; exponents and coefficients, (<gensym> E1 C1 E2 C2 ...) = C1*<gensym>^E1 + C2*<gensym>^E2 + ....
3513 ;; (%i73) partfrac ((2*x^2-3)/(x^4-3*x^2+2), x);
3514 ;; 1. Trace: (PARTFRAC '((#:X16910 2 2 0 -3) #:X16910 4 1 2 -3 0 2) '#:X16910)
3515 ;; 2. Trace: (BPROG '(#:X16910 2 1 0 -2) '(#:X16910 2 1 0 -1))
3516 ;; 2. Trace: BPROG ==> ((-1 . 1) 1 . 1)
3517 ;; 2. Trace: (BPROG '(#:X16910 1 1 0 1) '(#:X16910 1 1 0 -1))
3518 ;; 2. Trace: BPROG ==> ((1 . 2) -1 . 2)
3519 ;; 2. Trace: (BPROG '(#:X16910 1 1 0 -1) '1)
3520 ;; 2. Trace: BPROG ==> ((0 . 1) 1 . 1)
3521 ;; 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))
3522 ;; (%o73) 1/(x^2-2)-1/(2*(x+1))+1/(2*(x-1))
3524 (defun bprog (r s bprog-var
)
3525 (prog (p1b p2b coef1r coef2r coef1s coef2s f1 f2 a egcd oldalg
)
3526 (setq oldalg $algebraic
)
3529 (setq coef2r
(setq coef1s
0))
3530 (setq coef2s
(setq coef1r
1))
3533 (unless (zerop (pdegree p1b bprog-var
)) (setq egcd
(pgcdexpon p1b
)))
3535 (unless (or (zerop (pdegree p2b bprog-var
)) (= egcd
1))
3536 (setq egcd
(gcd egcd
(pgcdexpon p2b
)))
3537 (setq p1b
(pexpon*// p1b egcd nil
)
3538 p2b
(pexpon*// p2b egcd nil
)))
3539 b1
(cond ((< (pdegree p1b bprog-var
) (pdegree p2b bprog-var
))
3541 (rotatef coef1r coef2r
)
3542 (rotatef coef1s coef2s
)))
3543 (when (zerop (pdegree p2b bprog-var
))
3544 (unless (zerop (pdegree coef2r bprog-var
))
3545 (setq coef2r
(pexpon*// coef2r egcd t
)))
3546 (unless (zerop (pdegree coef2s bprog-var
))
3547 (setq coef2s
(pexpon*// coef2s egcd t
)))
3548 (return (cons (ratreduce (ptimes (cdr r
) coef2r
) p2b
)
3549 (ratreduce (ptimes (cdr s
) coef2s
) p2b
))))
3550 (setq f1
(psquorem1 (cdr p1b
) (cdr p2b
) t
))
3551 (setq $algebraic $false
)
3552 (setq f2
(psimp bprog-var
(cadr f1
)))
3553 (setq p1b
(pquotientchk (psimp bprog-var
(caddr f1
)) a
))
3554 (setq $algebraic oldalg
)
3556 (setq coef1r
(pquotientchk (pdifference (ptimes f1 coef1r
)
3559 (setq coef1s
(pquotientchk (pdifference (ptimes f1 coef1s
)
3565 (defun ratdifference (a b
) (ratplus a
(ratminus b
)))
3567 (defun ratpl (a b
) (ratplus (ratfix a
) (ratfix b
)))
3569 (defun ratti (a b c
) (rattimes (ratfix a
) (ratfix b
) c
))
3571 (defun ratqu (a b
) (ratquotient (ratfix a
) (ratfix b
)))
3573 (defun ratfix (a) (cond ((equal a
(ratnumerator a
)) (cons a
1)) (t a
)))
3575 (defun ratdivide (f g
)
3576 (destructuring-let* (((fnum . fden
) (ratfix f
))
3577 ((gnum . gden
) (ratfix g
))
3578 ((q r
) (pdivide fnum gnum
)))
3579 (cons (ratqu (ratti q gden t
) fden
)
3582 (defun polcoef (l n polcoef-var
)
3583 (cond ((or (atom l
) (pointergp polcoef-var
(car l
)))
3584 (cond ((equal n
0) l
)
3587 (ptterm (cdr l
) n
))))
3589 (defun disrep (l disrep-ratform
)
3590 (cond ((equal (ratnumerator l
) l
)
3591 ($ratdisrep
(cons disrep-ratform
(cons l
1))))
3592 (t ($ratdisrep
(cons disrep-ratform l
)))))
3594 ;; The following was formerly in MATRUN. This code was placed here because
3595 ;; MATRUN is now an out-of-core file on MC, and this code is needed in-core
3596 ;; so that MACSYMA SAVE files will work. - JPG
3601 (defun kar (x) (if (atom x
) (matcherr) (car x
)))
3603 (defun kaar (x) (kar (kar x
)))
3605 (defun kdr (x) (if (atom x
) (matcherr) (cdr x
)))
3607 (defun simpargs1 (a vestigial c
)
3608 (declare (ignore vestigial
))
3612 (if (not (atom x
)) (car x
)))
3614 (defquote retlist
(&rest l
)
3616 (mapcar #'(lambda (z) (list '(mequal simp
) z
(meval z
))) l
)))
3619 (if (zerop c
) x
(nthkdr (kdr x
) (1- c
))))
