1 ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancments. ;;;;;
5 ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
6 ;;; All rights reserved ;;;;;
7 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
8 ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
9 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
15 (declare-top (special rulesw
*inv
* substp limitp
16 prods negprods sums negsums
17 $scalarmatrixp
*nounl
*
20 bigfloatzero bigfloatone $assumescalar
21 opers-list
*opers-list $dontfactor
*n
22 *out
*in varlist genvar $factorflag radcanp
23 *builtin-numeric-constants
*))
25 ;; General purpose simplification and conversion switches.
27 (defmvar $negdistrib t
28 "Causes negations to be distributed over sums, e.g. -(A+B) is
32 "Causes SOME mathematical functions (including exponentiation)
33 with numerical arguments to be evaluated in floating point.
34 It causes variables in an expression which have been given
35 NUMERVALs to be replaced by their values. It also turns
37 see-also
($numerval $float
))
39 (defmvar $simp t
"Enables simplification.")
41 (defmvar $sumexpand nil
42 "If TRUE, products of sums and exponentiated sums go into nested
45 (defmvar $numer_pbranch nil
)
47 ;; Switches dealing with matrices and non-commutative multiplication.
49 (defmvar $doscmxplus nil
50 "Causes SCALAR + MATRIX to return a matrix answer. This switch
51 is not subsumed under DOALLMXOPS.")
54 "Causes SCALAR^MATRIX([1,2],[3,4]) to return
55 MATRIX([SCALAR,SCALAR^2],[SCALAR^3,SCALAR^4]). In general, this
56 transformation affects exponentiations where the *print-base* is a
57 scalar and the power is a matrix or list.")
59 (defmvar $domxplus nil
)
61 (defmvar $domxtimes nil
)
65 ;; Switches dealing with expansion.
68 "The largest positive exponent which will be automatically
69 expanded. (X+1)^3 will be automatically expanded if
70 EXPOP is greater than or equal to 3."
72 see-also
($expon $maxposex $expand
))
75 "The largest negative exponent which will be automatically
76 expanded. (X+1)^(-3) will be automatically expanded if
77 EXPON is greater than or equal to 3."
79 see-also
($expop $maxnegex $expand
))
81 (defmvar $maxposex
1000.
82 "The largest positive exponent which will be expanded by
85 see-also
($maxnegex $expop $expand
))
87 ;; Check assignment to be a positive integer
88 (putprop '$maxposex
'posintegerset
'assign
)
90 (defmvar $maxnegex
1000.
91 "The largest negative exponent which will be expanded by
94 see-also
($maxposex $expon $expand
))
96 ;; Check assignment to be a positive integer
97 (putprop '$maxnegex
'posintegerset
'assign
)
99 ;; Lisp level variables
102 "Causes SIMP flags to be ignored. $EXPAND works by binding
103 $EXPOP to $MAXPOSEX, $EXPON to $MAXNEGEX, and DOSIMP to T.")
106 "Causes a throw to the tag ERRORSW when certain errors occur
107 rather than the printing of a message. Kludgy MAXIMA-SUBSTITUTE for
108 MAXIMA-ERROR signalling.")
110 (defmvar derivsimp t
"Hack in `simpderiv' for `rwg'")
112 (defmvar $rootsepsilon
#+gcl
(float 1/10000000) #-gcl
1e-7)
113 (defmvar $grindswitch nil
)
114 (defmvar $algepsilon
100000000)
119 (defmvar $logabs nil
)
120 (defmvar $limitdomain
'$complex
)
121 (defmvar $listarith t
)
122 (defmvar $domain
'$real
)
123 (defmvar $m1pbranch nil
)
124 (defmvar $%e_to_numlog nil
)
126 (defmvar $lognegint nil
)
127 (defmvar $ratsimpexpons nil
)
128 (defmvar $logexpand t
) ; Possible values are T, $ALL and $SUPER
129 (defmvar $radexpand t
)
130 (defmvar $subnumsimp nil
)
132 (defmvar $distribute_over t
) ; If T, functions are distributed over bags.
135 (defvar rp-polylogp nil
)
138 (defvar timesinp nil
)
139 (defvar %e-val
(mget '$%e
'$numer
))
140 (defvar %pi-val
(mget '$%pi
'$numer
))
141 (defvar derivflag nil
)
142 (defvar exptrlsw nil
)
143 (defvar expandflag nil
)
144 (defvar *zexptsimp? nil
)
147 (defprop mnctimes t associative
)
148 (defprop lambda t lisp-no-simp
)
150 ;; Local functions should not be simplified. Various lisps
151 ;; use various names for the list structure defining these:
154 #-gcl
(:load-toplevel
)
157 (dolist (y (list x z
))
160 (setf (get (car y
) 'lisp-no-simp
) t
))))))
162 (dolist (x '(mplus mtimes mnctimes mexpt mncexpt %sum
))
163 (setf (get x
'msimpind
) (cons x
'(simp))))
165 ;; operators properties
167 (mapc #'(lambda (x) (setf (get (first x
) 'operators
) (second x
)))
168 '((mplus simplus
) (mtimes simptimes
) (mncexpt simpncexpt
)
169 (mminus simpmin
) (%gamma simpgamma
) (mfactorial simpfact
)
170 (mnctimes simpnct
) (mquotient simpquot
) (mexpt simpexpt
)
172 (%derivative simpderiv
)
174 (%integrate simpinteg
) (%limit simp-limit
)
175 (bigfloat simpbigfloat
) (lambda simplambda
) (mdefine simpmdef
)
176 (mqapply simpmqapply
) (%gamma simpgamma
)
177 ($beta simpbeta
) (%sum simpsum
) (%binomial simpbinocoef
)
178 (%plog simpplog
) (%product simpprod
) (%genfact simpgfact
)
179 ($atan2 simpatan2
) ($matrix simpmatrix
) (%matrix simpmatrix
)
180 ($bern simpbern
) ($euler simpeuler
)))
182 (defprop $li lisimp specsimp
)
183 (defprop $psi psisimp specsimp
)
185 (defprop $equal t binary
)
186 (defprop $notequal t binary
)
188 (defmfun $bfloatp
(x)
189 "Returns true if X is a bigfloat"
192 (eq (caar x
) 'bigfloat
)))
195 "Returns non-NIL if X is an integer, float, or bfloat that is equal
197 (or (and (integerp x
) (= 0 x
))
198 (and (floatp x
) (= 0.0 x
))
199 (and ($bfloatp x
) (= 0 (second x
)))))
202 "Returns non-NIL if X is an integer, float, or bfloat that is equal
204 (or (and (integerp x
) (= 1 x
))
205 (and (floatp x
) (= 1.0 x
))
206 (and ($bfloatp x
) (zerop1 (sub x
1)))))
209 "Returns non-NIL if X is a Lisp number or if it is a Maxima rational
210 form or a bigfloat form"
212 (and (not (atom x
)) (not (atom (car x
)))
213 (member (caar x
) '(rat bigfloat
)))))
215 ;; Does X or a subexpression match PREDICATE?
217 ;; If X is a tree then we recurse depth-first down its arguments. The specrep
218 ;; check is because rat forms are built rather differently from normal Maxima
219 ;; expressions so we need to unpack them for the recursion to work properly.
220 (defun subexpression-matches-p (predicate x
)
221 (or (funcall predicate x
)
224 (subexpression-matches-p predicate
(specdisrep x
))
225 (some (lambda (arg) (subexpression-matches-p predicate arg
))
228 ;; Is there a bfloat anywhere in X?
229 (defun some-bfloatp (x) (subexpression-matches-p '$bfloatp x
))
231 ;; Is there a float anywhere in X?
232 (defun some-floatp (x) (subexpression-matches-p 'floatp x
))
234 ;; EVEN works for any arbitrary lisp object since it does an integer
235 ;; check first. In other cases, you may want the Lisp EVENP function
236 ;; which only works for integers.
238 (defun even (a) (and (integerp a
) (not (oddp a
))))
240 ;; Predicates to determine if X satisfies some condition.
243 "Determines if X is a Maxima rational form: ((rat ...) a b)"
244 (and (not (atom x
)) (eq (caar x
) 'rat
)))
247 "Determines if X is a Maxima sum form: ((mplus ...) ...)"
248 (and (not (atom x
)) (eq (caar x
) 'mplus
)))
251 "Determines if X is a Maxima product form: ((mtimes ...) ...)"
252 (and (not (atom x
)) (eq (caar x
) 'mtimes
)))
255 "Determines if X is a Maxima exponential form: ((mexpt ...) ...)"
256 (and (not (atom x
)) (eq (caar x
) 'mexpt
)))
258 (defun mnctimesp (x) (and (not (atom x
)) (eq (caar x
) 'mnctimes
)))
260 (defun mncexptp (x) (and (not (atom x
)) (eq (caar x
) 'mncexpt
)))
263 "Determines if X is a Maxima log form: ((%log ...) ...)"
264 (and (not (atom x
)) (eq (caar x
) '%log
)))
267 "Determines if X is a Maxima negative form: ((mminus ...) ...)
269 This generally only happens on input forms like a - b:
270 ((mplus) $a ((mminus) $b)).
271 After simplification a - b becomes
272 ((mplus) $a ((mtimes) -1 $b))"
274 (and (not (atom x
)) (eq (caar x
) 'mminus
)))
277 "Determines if X is negative if X is a Lisp number or a Maxima rat
278 form or bigfloat form"
279 (cond ((realp x
) (minusp x
))
280 ((or (ratnump x
) ($bfloatp x
)) (minusp (cadr x
)))))
282 (defun mqapplyp (e) (and (not (atom e
)) (eq (caar e
) 'mqapply
)))
284 (defun ratdisrep (e) (simplifya ($ratdisrep e
) nil
))
286 (defun sratsimp (e) (simplifya ($ratsimp e
) nil
))
288 (defun simpcheck (e flag
)
289 (cond ((specrepp e
) (specdisrep e
))
291 (t (let (($%enumer $numer
))
292 ;; Switch $%enumer on, when $numer is TRUE to allow
293 ;; simplification of $%e to its numerical value.
294 (simplifya e nil
)))))
296 (defun mratcheck (e) (if ($ratp e
) (ratdisrep e
) e
))
298 (defmfun $numberp
(e)
299 "Returns true if E is a Maxima rational, a float, or a bigfloat number"
300 (or ($ratnump e
) ($floatnump e
) ($bfloatp e
)))
302 (defmfun $integerp
(x)
305 (not (member 'trunc
(car x
)))
307 (equal (cddr x
) 1))))
309 ;; The call to $INTEGERP in the following two functions checks for a CRE
310 ;; rational number with an integral numerator and a unity denominator.
313 (cond ((integerp x
) (oddp x
))
314 (($integerp x
) (oddp (cadr x
)))))
317 (cond ((integerp x
) (evenp x
))
318 (($integerp x
) (not (oddp (cadr x
))))))
320 (defmfun $floatnump
(x)
322 (and ($ratp x
) (floatp (cadr x
)) (onep1 (cddr x
)))))
324 (defmfun $ratnump
(x)
328 (not (member 'trunc
(car x
)))
330 (integerp (cddr x
)))))
335 (eq (caar x
) 'mrat
)))
337 (defmfun $taylorp
(x)
340 (member 'trunc
(cdar x
)) t
))
342 (defun specrepcheck (e) (if (specrepp e
) (specdisrep e
) e
))
344 ;; Note that the following two functions are carefully coupled.
348 (member (caar e
) '(mrat mpois
))))
350 (defun specdisrep (e)
351 (cond ((eq (caar e
) 'mrat
) (ratdisrep e
))
354 (defmfun $polysign
(x)
355 (setq x
(cadr (ratf x
)))
356 (cond ((equal x
0) 0) ((pminusp x
) -
1) (t 1)))
358 ;; These check for the correct number of operands within Macsyma expressions,
359 ;; not arguments in a procedure call as the name may imply.
361 (defun arg-count-check (required-arg-count expr
)
362 (unless (= required-arg-count
(length (rest expr
)))
363 (wna-err expr required-arg-count
)))
365 (defun oneargcheck (expr)
366 (arg-count-check 1 expr
))
368 (defun twoargcheck (expr)
369 (arg-count-check 2 expr
))
371 ;; WNA-ERR: Wrong Number of Arguments error
373 ;; If REQUIRED-ARG-COUNT is non-NIL, then we check that EXPR has the
374 ;; correct number of arguments. A informative error message is shown
375 ;; if the number of arguments is not given.
377 ;; Otherwise, EXPR must be a symbol and a generic message is printed.
378 ;; (This is for backward compatibility for existing uses of WNA-ERR.)
379 (defun wna-err (exprs &optional required-arg-count
)
380 (if required-arg-count
381 (let ((op (caar exprs
))
382 (actual-count (length (rest exprs
))))
383 (merror (intl:gettext
"~M: expected exactly ~M arguments but got ~M: ~M")
384 op required-arg-count actual-count
(list* '(mlist) (rest exprs
))))
385 (merror (intl:gettext
"~:@M: wrong number of arguments.")
388 (defun improper-arg-err (exp fn
)
389 (merror (intl:gettext
"~:M: improper argument: ~M") fn exp
))
391 (defun subargcheck (form subsharp argsharp fun
)
392 (if (or (not (= (length (subfunsubs form
)) subsharp
))
393 (not (= (length (subfunargs form
)) argsharp
)))
394 (merror (intl:gettext
"~:@M: wrong number of arguments or subscripts.") fun
)))
396 ;; Constructor and extractor primitives for subscripted functions, e.g.
397 ;; F[1,2](X,Y). SUBL is (1 2) and ARGL is (X Y).
399 ;; These will be flushed when NOPERS is finished. They will be macros in
400 ;; NOPERS instead of functions, so we have to be careful that they aren't
401 ;; mapped or applied anyplace. What we really want is open-codable routines.
403 (defun subfunmakes (fun subl argl
)
404 `((mqapply simp
) ((,fun simp array
) .
,subl
) .
,argl
))
406 (defun subfunmake (fun subl argl
)
407 `((mqapply) ((,fun simp array
) .
,subl
) .
,argl
))
409 (defun subfunname (exp) (caaadr exp
))
411 (defun subfunsubs (exp) (cdadr exp
))
413 (defun subfunargs (exp) (cddr exp
))
415 (defmfun $numfactor
(x)
416 (setq x
(specrepcheck x
))
419 ((not (eq (caar x
) 'mtimes
)) 1)
420 ((mnump (cadr x
)) (cadr x
))
423 (defun scalar-or-constant-p (x flag
)
424 (if flag
(not ($nonscalarp x
)) ($scalarp x
)))
426 (defmfun $constantp
(x)
427 (cond ((atom x
) (or ($numberp x
) (kindp x
'$constant
)))
428 ((member (caar x
) '(rat bigfloat
)) t
)
429 ((specrepp x
) ($constantp
(specdisrep x
)))
430 ((or (mopp (caar x
)) (kindp (caar x
) '$constant
))
431 (do ((x (cdr x
) (cdr x
))) ((null x
) t
)
432 (if (not ($constantp
(car x
))) (return nil
))))))
435 (cond ((symbolp x
) (kindp x
'$constant
))
437 (and (kindp (caar x
) '$constant
)
438 (do ((x (cdr x
) (cdr x
))) ((null x
) t
)
439 (if (not ($constantp
(car x
))) (return nil
)))))))
441 (defun maxima-constantp (x)
443 (and (symbolp x
) (kindp x
'$constant
))))
445 (defun consttermp (x) (and ($constantp x
) (not ($nonscalarp x
))))
447 (defmfun $scalarp
(x) (or (consttermp x
) (eq (scalarclass x
) '$scalar
)))
449 (defmfun $nonscalarp
(x) (eq (scalarclass x
) '$nonscalar
))
451 (defun scalarclass (exp) ; Returns $SCALAR, $NONSCALAR, or NIL (unknown).
453 ;; Maxima numbers are scalar.
456 (cond ((or (mget exp
'$nonscalar
)
457 (and (not (mget exp
'$scalar
))
458 ;; Arrays are nonscalar, but not if declared scalar.
460 ($member exp $arrays
))))
462 ((or (mget exp
'$scalar
)
463 ;; Include constant atoms which are not declared nonscalar.
466 ((member 'array
(car exp
))
467 (cond ((mget (caar exp
) '$scalar
) '$scalar
)
468 ((mget (caar exp
) '$nonscalar
) '$nonscalar
)
470 ((specrepp exp
) (scalarclass (specdisrep exp
)))
471 ;; If the function is declared scalar or nonscalar, then return. If it
472 ;; isn't explicitly declared, then try to be intelligent by looking at
473 ;; the arguments to the function.
474 ((scalarclass (caar exp
)))
475 ;; <number> + <scalar> is SCALARP because that seems to be useful.
476 ;; This should probably only be true if <number> is a member of the
477 ;; field of scalars. <number> * <scalar> is SCALARP since
478 ;; <scalar> + <scalar> is SCALARP. Also, this has to be done to make
479 ;; <scalar> - <scalar> SCALARP.
480 ((member (caar exp
) '(mplus mtimes
))
481 (do ((l (cdr exp
) (cdr l
))) ((null l
) '$scalar
)
482 (if (not (consttermp (car l
)))
483 (return (scalarclass-list l
)))))
484 ((and (eq (caar exp
) 'mqapply
) (scalarclass (cadr exp
))))
485 ((mxorlistp exp
) '$nonscalar
)
486 ;; If we can't find out anything about the operator, then look at the
487 ;; arguments to the operator. I think NIL should be returned at this
490 (do ((exp (cdr exp
) (cdr exp
)) (l '(1)))
491 ((null exp
) (scalarclass-list l
))
492 (if (not (consttermp (car exp
)))
493 (setq l
(cons (car exp
) l
)))))))
495 ;; Could also do <scalar> +|-|*|/ |^ <declared constant>, but this is not
496 ;; always correct and could screw somebody.
498 ;; SCALARCLASS-LIST takes a list of expressions as its argument. If their
499 ;; scalarclasses all agree, then that scalarclass is returned.
501 (defun scalarclass-list (llist)
502 (cond ((null llist
) nil
)
503 ((null (cdr llist
)) (scalarclass (car llist
)))
504 (t (let ((sc-car (scalarclass (car llist
)))
505 (sc-cdr (scalarclass-list (cdr llist
))))
506 (cond ((or (eq sc-car
'$nonscalar
)
507 (eq sc-cdr
'$nonscalar
))
509 ((and (eq sc-car
'$scalar
) (eq sc-cdr
'$scalar
))
514 (member (caar x
) '(mequal mlist $matrix
))))
516 (defun mequalp (x) (and (not (atom x
)) (eq (caar x
) 'mequal
)))
520 (member (caar x
) '(mlist $matrix
))))
522 (defun mxorlistp1 (x)
524 (or (eq (caar x
) '$matrix
)
525 (and (eq (caar x
) 'mlist
) $listarith
))))
527 (defun constfun (ign)
528 (declare (ignore ign
)) ; Arg ignored. Function used for mapping down lists.
531 (defun constmx (*const
* x
)
532 (simplifya (fmapl1 'constfun x
) t
))
534 ;;; ISINOP returns the complete subexpression with the operator OP, when the
535 ;;; operator OP is found in EXPR.
537 (defun isinop (expr op
) ; OP is assumed to be an atom
538 (cond ((atom expr
) nil
)
539 ((and (eq (caar expr
) op
)
540 (not (member 'array
(cdar expr
))))
543 (do ((expr (cdr expr
) (cdr expr
))
546 (when (setq res
(isinop (car expr
) op
))
549 (defun free (exp var
)
550 (cond ((alike1 exp var
) nil
)
553 (and (listp (car exp
))
554 (free (caar exp
) var
)
555 (freel (cdr exp
) var
)))))
558 (do ((l l
(cdr l
))) ((null l
) t
)
560 ((atom l
) (return (free l var
))) ;; second element of a pair
561 ((not (free (car l
) var
)) (return nil
)))))
564 (defun freeargs (exp var
)
565 (cond ((alike1 exp var
) nil
)
567 (t (do ((l (margs exp
) (cdr l
))) ((null l
) t
)
568 (cond ((not (freeargs (car l
) var
)) (return nil
)))))))
570 (defun simplifya (x y
)
571 (cond ((not $simp
) x
)
573 (cond ((and $%enumer $numer
(eq x
'$%e
))
574 ;; Replace $%e with its numerical value,
575 ;; when %enumer and $numer TRUE
579 (cond ((and (cdr x
) (atom (cdr x
)))
580 (merror (intl:gettext
"simplifya: malformed expression (atomic cdr).")))
581 ((get (car x
) 'lisp-no-simp
)
582 ;; this feature is to be used with care. it is meant to be
583 ;; used to implement data objects with minimum of consing.
584 ;; forms must not bash the DISPLA package. Only new forms
585 ;; with carefully chosen names should use this feature.
588 (mapcar #'(lambda (x) (simplifya x y
)) (cdr x
))))))
589 ((eq (caar x
) 'rat
) (*red1 x
))
590 ;; Enforced resimplification: Reset dosimp and strip 'simp tags from x.
591 (dosimp (let ((dosimp nil
)) (simplifya (unsimplify x
) y
)))
592 ((member 'simp
(cdar x
)) x
)
593 ((eq (caar x
) 'mrat
) x
)
594 ((not (atom (caar x
)))
595 (cond ((or (eq (caaar x
) 'lambda
)
596 (and (not (atom (caaar x
))) (eq (caaaar x
) 'lambda
)))
597 (mapply1 (caar x
) (cdr x
) (caar x
) x
))
598 (t (merror (intl:gettext
"simplifya: operator is neither an atom nor a lambda expression: ~S") x
))))
599 ((and $distribute_over
600 (get (caar x
) 'distribute_over
)
601 ;; A function with the property 'distribute_over.
602 ;; Look, if we have a bag as argument to the function.
603 (distribute-over x
)))
604 ((get (caar x
) 'opers
)
605 (let ((opers-list *opers-list
)) (oper-apply x y
)))
606 ((and (eq (caar x
) 'mqapply
)
608 (and (eq substp
'mqapply
)
609 (or (eq (car (cadr x
)) 'lambda
)
610 (eq (caar (cadr x
)) 'lambda
)))))
611 (cond ((or (symbolp (cadr x
)) (not (atom (cadr x
))))
612 (simplifya (cons (cons (cadr x
) (cdar x
)) (cddr x
)) y
))
613 ((or (not (member 'array
(cdar x
))) (not $subnumsimp
))
614 (merror (intl:gettext
"simplifya: I don't know how to simplify this operator: ~M") x
))
616 (t (let ((w (get (caar x
) 'operators
)))
618 (or (not (member 'array
(cdar x
)))
621 (t (simpargs x y
)))))))
623 ;; EQTEST returns an expression which is the same as X
624 ;; except that it is marked with SIMP and maybe other flags from CHECK.
626 ;; Following description is inferred from the code. Dunno why the function is named "EQTEST".
628 ;; (1) if X is already marked with SIMP flag or doesn't need it: return X.
629 ;; (2) if X is pretty much the same as CHECK (same operator and same arguments),
630 ;; then return CHECK after marking it with SIMP flag.
631 ;; (3) if operator of X has the MSIMPIND property, replace it
632 ;; with value of MSIMPIND (something like '(MPLUS SIMP)) and return X.
633 ;; (4) if X or CHECK is an array expression, return X after marking it with SIMP and ARRAY flags.
634 ;; (5) otherwise, return X after marking it with SIMP flag.
636 (defun eqtest (x check
)
641 (member 'simp
(cdar x
)))
643 ((and (eq (caar x
) (caar check
))
644 (equal (cdr x
) (cdr check
)))
645 (cond ((and (null (cdar check
))
646 (setq y
(get (caar check
) 'msimpind
)))
647 (cons y
(cdr check
)))
648 ((member 'simp
(cdar check
))
651 (cons (cons (caar check
)
653 (cons 'simp
(cdar check
))
656 ((setq y
(get (caar x
) 'msimpind
))
658 ((or (member 'array
(cdar x
))
659 (and (eq (caar x
) (caar check
))
660 (member 'array
(cdar check
))))
661 (rplaca x
(cons (caar x
) '(simp array
))))
663 (rplaca x
(cons (caar x
) '(simp)))))))
665 ;; A function, which distributes of bags like a list, matrix, or equation.
666 ;; Check, if we have to distribute of one of the bags or any other operator.
667 (defun distribute-over (expr)
668 (cond ((= 1 (length (cdr expr
)))
669 ;; Distribute over for a function with one argument.
670 (cond ((and (not (atom (cadr expr
)))
671 (member (caaadr expr
) (get (caar expr
) 'distribute_over
)))
674 (mapcar #'(lambda (u) (simplify (list (car expr
) u
)))
678 ;; A function with more than one argument.
679 (do ((args (cdr expr
) (cdr args
))
682 (when (and (not (atom (car args
)))
683 (member (caar (car args
))
684 (get (caar expr
) 'distribute_over
)))
685 ;; Distribute the function over the arguments and simplify again.
688 (cons (ncons (caar (car args
)))
689 (mapcar #'(lambda (u)
693 (cons (ncons (caar expr
))
694 (reverse first-args
))
697 (cdr (car args
)))))))
698 (setq first-args
(cons (car args
) first-args
))))))
700 (defun rulechk (x) (or (mget x
'oldrules
) (get x
'rules
)))
702 (defun resimplify (x) (let ((dosimp t
)) (simplifya x nil
)))
704 (defun unsimplify (x)
705 (if (or (atom x
) (specrepp x
))
707 (cons (remove 'simp
(car x
) :count
1) (mapcar #'unsimplify
(cdr x
)))))
709 (defun simpargs (x y
)
710 (if (or (and (eq (get (caar x
) 'dimension
) 'dimension-infix
)
711 (not (getl (caar x
) '($lassociative $rassociative
))))
712 (get (caar x
) 'binary
))
714 (if (and (member 'array
(cdar x
)) (null (margs x
)))
715 (merror (intl:gettext
"SIMPARGS: subscripted variable found with no subscripts.")))
716 (eqtest (if y x
(let ((flag (member (caar x
) '(mlist mequal
))))
717 (cons (ncons (caar x
))
718 (mapcar #'(lambda (u)
719 (if flag
(simplifya u nil
)
724 ;;;-----------------------------------------------------------------------------
725 ;;; ADDK (X Y) 27.09.2010/DK
727 ;;; Arguments and values:
728 ;;; X - a Maxima number
729 ;;; Y - a Maxima number
730 ;;; result - a simplified Maxima number
733 ;;; ADDK adds two Maxima numbers and returns a simplified Maxima number.
734 ;;; ADDK can be called in Lisp code, whenever the arguments are valid
735 ;;; Maxima numbers, these are integer, float, Maxima rational, or
736 ;;; Maxima bigfloat numbers. The arguments must not be simplified. The
737 ;;; precision of a bigfloat result depends on the setting of the
738 ;;; global variable $FPPREC. If the option variable $FLOAT is T, a
739 ;;; Maxima rational number as a result is converted to a float number.
743 ;;; (addk 2.0 3) -> 5.0
744 ;;; (addk ($bfloat 2) 3)-> ((BIGFLOAT SIMP 56) 45035996273704960 3)
745 ;;; (addk 2 '((rat) 1 2)) -> ((RAT SIMP) 5 2)
746 ;;; (let (($float t)) (addk 2 '((rat) 1 2))) -> 2.5
749 ;;; The option variables $FLOAT and $FPPREC.
752 ;;; TIMESK to multiply and EXPTRL to exponentiate two Maxima numbers.
755 ;;; The routine works for Lisp rational and Lisp complex numbers too.
756 ;;; This feature is not used in Maxima code. If Lisp complex and
757 ;;; rational numbers are mixed with Maxima rational or bigfloat
758 ;;; numbers the result is wrong or a Lisp error is generated.
759 ;;;-----------------------------------------------------------------------------
764 ((and (numberp x
) (numberp y
)) (+ x y
))
765 ((or ($bfloatp x
) ($bfloatp y
)) ($bfloat
(list '(mplus) x y
)))
768 (cond ((floatp x
) (return (+ x
(fpcofrat y
))))
769 (t (setq x
(list '(rat) x
1)))))
771 (cond ((floatp y
) (return (+ y
(fpcofrat x
))))
772 (t (setq y
(list '(rat) y
1))))))
773 (setq g
(gcd (caddr x
) (caddr y
)))
774 (setq a
(truncate (caddr x
) g
)
775 b
(truncate (caddr y
) g
))
776 (return (timeskl (list '(rat) 1 g
)
782 ;;;-----------------------------------------------------------------------------
783 ;;; *RED1 (X) 27.09.2010/DK
786 ;;; Arguments and values:
787 ;;; X - a Maxima rational number (for *RED1)
788 ;;; N - an integer number representing the numerator of a rational
789 ;;; D - an integer number representing the denominator of a rational
790 ;;; result - a simplified Maxima rational number
793 ;;; *RED1 is called from SIMPLIFYA to reduce and simplify a Maxima rational
794 ;;; number. *RED1 checks if the rational number is already simplified. If
795 ;;; the option variable $FLOAT is T, the rational number is converted to a
796 ;;; float number. If the number is not simplified, *RED is called.
798 ;;; *RED reduces the numerator N and the demoniator D and returns a
799 ;;; simplified Maxima rational number. The result is converted to a float
800 ;;; number, if the option variable $FLOAT is T.
803 ;;; The option variable $FLOAT.
804 ;;;-----------------------------------------------------------------------------
807 (cond ((member 'simp
(cdar x
))
808 (if $float
(fpcofrat x
) x
))
809 (t (*red
(cadr x
) (caddr x
)))))
814 (t (let ((u (gcd n d
)))
815 (setq n
(truncate n u
)
817 (if (minusp d
) (setq n
(- n
) d
(- d
)))
818 (cond ((equal d
1) n
)
819 ($float
(fpcofrat1 n d
))
820 (t (list '(rat simp
) n d
)))))))
822 ;;;-----------------------------------------------------------------------------
823 ;;; TIMESK (X Y) 27.09.2010/DK
825 ;;; Arguments and values:
826 ;;; X - a Maxima number
827 ;;; Y - a Maxima number
828 ;;; result - a simplified Maxima number
831 ;;; TIMESK Multiplies two Maxima numbers and returns a simplified Maxima
832 ;;; number. TIMESK can be called in Lisp code, whenever the arguments are
833 ;;; valid Maxima numbers, these are integer, float, Maxima rational, or
834 ;;; Maxima bigfloat numbers. The arguments must not be simplified. The
835 ;;; precision of a bigfloat result depends on the setting of the
836 ;;; global variable $FPPREC. If the option variable $FLOAT is T, a
837 ;;; Maxima rational number as a result is converted to a float number.
839 ;;; TIMESKL is called from TIMESK to multiply two Maxima rational numbers or
840 ;;; a rational number with an integer number.
843 ;;; (timesk 2 3) -> 6
844 ;;; (timesk 2.0 3) -> 6.0
845 ;;; (timesk ($bfloat 2) 3)-> ((BIGFLOAT SIMP 56) 54043195528445952 3)
846 ;;; (timesk 3 '((rat) 1 2)) -> ((RAT SIMP) 3 2)
847 ;;; (let (($float t)) (timesk 3 '((rat) 1 2))) -> 1.5
850 ;;; The option variables $FLOAT and $FPPREC.
853 ;;; ADDK to add and EXPTRL to exponentiate two Maxima numbers.
856 ;;; The routine works for Lisp rational and Lisp complex numbers too.
857 ;;; This feature is not used in Maxima code. If Lisp complex and
858 ;;; rational numbers are mixed with Maxima rational or bigfloat
859 ;;; numbers the result is wrong or a Lisp error is generated.
860 ;;;-----------------------------------------------------------------------------
862 ;; NUM1 and DENOM1 are helper functions for TIMESKL to get the numerator and the
863 ;; denominator of an integer or Maxima rational number. For an integer the
864 ;; denominator is 1. Both functions are used at other places in Maxima code too.
867 (if (numberp a
) a
(cadr a
)))
870 (if (numberp a
) 1 (caddr a
)))
872 (defun timesk (x y
) ; X and Y are assumed to be already reduced
873 (cond ((equal x
1) y
)
875 ((and (numberp x
) (numberp y
)) (* x y
))
876 ((or ($bfloatp x
) ($bfloatp y
)) ($bfloat
(list '(mtimes) x y
)))
877 ((floatp x
) (* x
(fpcofrat y
)))
878 ((floatp y
) (* y
(fpcofrat x
)))
881 ;; TIMESKL takes one or two Maxima rational numbers, one argument can be an
882 ;; integer number. The result is a Maxima rational or an integer number.
883 ;; If the option variable $FLOAT is T, a Maxima rational number in converted
888 (setq u
(*red
(num1 x
) (denom1 y
)))
889 (setq v
(*red
(num1 y
) (denom1 x
)))
890 (setq g
(cond ((or (equal u
0) (equal v
0)) 0)
892 ((and (numberp u
) (numberp v
)) (* u v
))
894 (* (num1 u
) (num1 v
))
895 (* (denom1 u
) (denom1 v
))))))
896 (return (cond ((numberp g
) g
)
897 ((equal (caddr g
) 1) (cadr g
))
898 ($float
(fpcofrat g
))
901 ;;;-----------------------------------------------------------------------------
902 ;;; FPCOFRAT (RATNO) 27.09.2010/DK
905 ;;; Arguments and values:
906 ;;; RATNO - a Maxima rational number (for FPCOFRAT)
907 ;;; NU - an integer number which represents the numerator of a rational
908 ;;; D - an integer number which represents the denominator of a rational
909 ;;; result - floating point approximation of a rational number
912 ;;; Floating Point Conversion OF RATional number routine.
913 ;;; Finds floating point approximation to rational number.
915 ;;; FPCOFRAT1 computes the quotient of NU/D.
917 ;;; Exceptional situations:
918 ;;; A Lisp error is generated, if the rational number does not fit into a
920 ;;;-----------------------------------------------------------------------------
922 ;; This constant is only needed in the file float.lisp.
924 #+gcl
(compile load eval
)
925 #-gcl
(:compile-toplevel
:load-toplevel
:execute
)
926 (defconstant machine-mantissa-precision
(float-digits 1.0)))
928 (defun fpcofrat (ratno)
929 (fpcofrat1 (cadr ratno
) (caddr ratno
)))
931 (defun fpcofrat1 (nu d
)
934 ;;;-----------------------------------------------------------------------------
935 ;;; EXPTA (X Y) 27.09.2010/DK
937 ;;; Arguments and values:
938 ;;; X - a Maxima number
939 ;;; Y - an integer number
940 ;;; result - a simplified Maxima number
943 ;;; Computes X^Y, where X is Maxima number and Y an integer. The result is
944 ;;; a simplified Maxima number. Y can be a rational Maxima number. For this
945 ;;; case the numerator is taken as the power.
948 ;;; The option variables $FLOAT and $FPPREC.
951 ;;; This routine is not used within the simplifier. There is only one
952 ;;; call from the file hayat.lisp. This call can be replaced with a
953 ;;; call of the function power.
954 ;;;-----------------------------------------------------------------------------
962 ($bfloat
(list '(mexpt) x y
)))
964 (*red
(exptb (caddr x
) (- (num1 y
)))
965 (exptb (cadr x
) (- (num1 y
)))))
967 (*red
(exptb (cadr x
) (num1 y
))
968 (exptb (caddr x
) (num1 y
))))))
970 ;;;-----------------------------------------------------------------------------
971 ;;; EXPTB (A B) 27.09.2010/DK
973 ;;; Arguments and values:
974 ;;; A - a float or integer number
975 ;;; B - an integer number
976 ;;; result - a simplified Maxima number
979 ;;; Computes A^B, where A is a float or an integer number and B is an
980 ;;; integer number. The result is an integer, float, or Maxima
985 ;;; (exptb 3.0 2) -> 9.0
986 ;;; (exptb 3 -2) -> ((RAT SiMP) 1 9)
987 ;;; (let (($float t)) (exptb 3 -2)) -> 0.1111111111111111
990 ;;; The option variable $FLOAT.
993 ;;; EXPTB calls the Lisp functions EXP or EXPT to compute the result.
994 ;;;-----------------------------------------------------------------------------
998 (cond ((equal a %e-val
)
999 ;; Make B a float so we'll get double-precision result.
1001 ((or (floatp a
) (not (minusp b
)))
1004 (setq b
(expt a
(- b
)))
1006 (if (float-inf-p result
) ;; needed for gcl - no trap of overflow
1007 (signal 'floating-point-overflow
)
1011 ;;;-----------------------------------------------------------------------------
1012 ;;; SIMPLUS (X W Z) 27.09.2010/DK
1014 ;;; Arguments and values:
1015 ;;; X - a Maxima expression of the form ((mplus) term1 term2 ...)
1016 ;;; W - an arbitrary value, the value is ignored
1017 ;;; Z - T or NIL, if T the arguments are assumed to be simplified
1018 ;;; result - a simplified mplus-expression or an atom
1021 ;;; Implementation of the simplifier for the "+" operator.
1022 ;;; A general description of SIMPLUS can be found in the paper:
1023 ;;; http://www.cs.berkeley.edu/~fateman/papers/simplifier.txt
1026 ;;; The addition of matrices and lists is affected by the following option
1028 ;;; $DOALLMXOPS, $DOMXMXOPS, $DOMXPLUS, $DOSCMXOPPS, $DOSCMXPLUS, $LISTARITH
1031 ;;; This routine should not be called directely. It is called by SIMPLIFYA.
1032 ;;; A save access is to call the function ADD.
1033 ;;;-----------------------------------------------------------------------------
1035 (defun simplus (x w z
)
1036 (prog (res check eqnflag matrixflag sumflag
)
1037 (if (null (cdr x
)) (return 0))
1041 (if (null x
) (go end
))
1042 (setq w
(if z
(car x
) (simplifya (car x
) nil
)))
1044 (cond ((atom w
) nil
)
1045 ((eq (caar w
) 'mrat
)
1049 (not (member 'trunc
(cdar w
))))
1050 (spsimpcases (cdr x
) w
))
1051 (setq w
(ratdisrep w
))
1055 (ratf (cons '(mplus)
1056 (nconc (mapcar #'simplify
(cons w
(cdr x
)))
1058 ((eq (caar w
) 'mequal
)
1063 (add2 (cadr eqnflag
) (cadr w
))
1064 (add2 (caddr eqnflag
) (caddr w
)))))
1066 ((member (caar w
) '(mlist $matrix
))
1068 (cond ((not matrixflag
) w
)
1069 ((and (or $doallmxops $domxmxops $domxplus
1070 (and (eq (caar w
) 'mlist
)
1071 ($listp matrixflag
)))
1072 (or (not (eq (caar w
) 'mlist
)) $listarith
))
1073 (addmx matrixflag w
))
1074 (t (setq res
(pls w res
)) matrixflag
)))
1076 ((eq (caar w
) '%sum
)
1077 (setq sumflag t res
(sumpls w res
))
1078 (setq w
(car res
) res
(cdr res
))))
1079 (setq res
(pls w res
))
1082 (setq res
(testp res
))
1085 (cond ((and (or ($listp matrixflag
)
1086 $doallmxops $doscmxplus $doscmxops
)
1087 (or (not ($listp matrixflag
)) $listarith
))
1088 (mxplusc res matrixflag
))
1089 (t (testp (pls matrixflag
(pls res nil
)))))))
1090 (setq res
(eqtest res check
))
1093 (add2 (cadr eqnflag
) res
)
1094 (add2 (caddr eqnflag
) res
))
1097 ;;;-----------------------------------------------------------------------------
1098 ;;; PLS (X OUT) 27.09.2010/DK
1100 ;;; Arguments and values:
1101 ;;; X - a Maxima expression or an atom
1102 ;;; OUT - a form ((mplus) <number> term1 term2 ...) or NIL
1103 ;;; result - a form ((mplus) <number> term1 ...), where x is added in.
1106 ;;; Adds the argument X into the form OUT. If OUT is NIL a form
1107 ;;; ((mplus) 0 X) is initialized, if X is an expression or a symbol,
1108 ;;; or ((mplus) X), if X is a number. Numbers are added to the first
1109 ;;; term <number> of the form. Any other symbol or expression is added
1110 ;;; into the canonical ordered list of arguments. The result is in a
1111 ;;; canonical order, but it is not a valid Maxima expression. To get a
1112 ;;; valid Maxima expression the result has to be checked with the
1113 ;;; function TESTP. This is done by the calling routine SIMPLUS.
1115 ;;; PLS checks the global flag *PLUSFLAG*, which is set in PLUSIN to T,
1116 ;;; if a mplus-expression is part of the result.
1119 ;;; (pls 2 nil) -> ((MPLUS) 2)
1120 ;;; (pls '$A nil) -> ((MPLUS) 0 $A)
1121 ;;; (pls '$B '((mplus) 0 $A)) -> ((MPLUS) 0 $A $B)
1122 ;;; (pls '$A '((mplus) 0 $A)) -> ((MPLUS) 0 ((MTIMES SIMP) 2 $A))
1124 ;;; Examples with the option variables $NUMER and $NEGDISTRIB:
1125 ;;; (let (($numer t)) (pls '$%e nil)) -> ((MPLUS) 2.718281828459045)
1126 ;;; (let (($negdistrib t)) (pls '((mtimes) -1 ((mplus) $A $B)) nil))
1127 ;;; -> ((MPLUS) 0 ((MTIMES SIMP) -1 $A) ((MTIMES SIMP) -1 $B))
1128 ;;; (let (($negdistrib nil)) (pls '((mtimes) -1 ((mplus) $A $B)) nil))
1129 ;;; -> ((MPLUS) 0 ((MTIMES) -1 ((MPLUS) $A $B)))
1132 ;;; The option variables $NUMER and $NEGDISTRIB and the global flag
1133 ;;; *PLUSFLAG*, which is set in the routine PLUSIN.
1136 ;;; PLUSIN and ADDK which are called from PLS and SIMPLUS.
1139 ;;; To add an expression into the list (CDR OUT), the list is passed
1140 ;;; to the routine PLUSIN as an argument. PLUSIN adds the argument to
1141 ;;; the list of terms by modifying the list (CDR OUT) destructively.
1142 ;;; The new value of OUT is returned as a result by PLS.
1143 ;;;-----------------------------------------------------------------------------
1145 ;; Set in PLUSIN to T to indicate a nested mplus expression.
1146 (defvar *plusflag
* nil
)
1148 ;; TESTP checks the result of PLS to get a valid Maxima mplus-expression.
1152 ((null (cddr x
)) (cadr x
))
1154 (cond ((null (cdddr x
)) (caddr x
)) (t (rplacd x
(cddr x
)))))
1158 (prog (fm *plusflag
*)
1159 (if (mtimesp x
) (setq x
(testtneg x
)))
1160 (when (and $numer
(atom x
) (eq x
'$%e
))
1161 ;; Replace $%e with its numerical value, when $numer ist TRUE
1164 ;; Initialize a form like ((mplus) <number> expr)
1167 (cond ((mnump x
) (ncons x
))
1169 (list 0 (cond ((atom x
) x
) (t (copy-list x
)))))
1170 ((mnump (cadr x
)) (copy-list (cdr x
) ))
1171 (t (cons 0 (copy-list (cdr x
) )))))))
1173 ;; Add a number into the first term of the list out.
1174 (return (cons '(mplus)
1175 (if (mnump (cadr out
))
1176 (cons (addk (cadr out
) x
) (cddr out
))
1177 (cons x
(cdr out
))))))
1178 ((not (mplusp x
)) (plusin x
(cdr out
)) (go end
)))
1179 ;; At this point we have a mplus expression as argument x. The following
1180 ;; code assumes that the argument x is already simplified and the terms
1181 ;; are in a canonical order.
1182 ;; First we add the number to the first term of the list out.
1184 (addk (if (mnump (cadr out
)) (cadr out
) 0)
1185 (cond ((mnump (cadr x
)) (setq x
(cdr x
)) (car x
)) (t 0))))
1186 ;; Initialize fm with the list of terms and start the loop to add the
1187 ;; terms of an mplus expression into the list out.
1190 (if (null (setq x
(cdr x
))) (go end
))
1191 ;; The return value of PLUSIN is a list, where the first element is the
1192 ;; added argument and the rest are the terms which follow the added
1194 (setq fm
(plusin (car x
) fm
))
1197 (if (not *plusflag
*) (return out
))
1198 (setq *plusflag
* nil
) ; *PLUSFLAG* T handles e.g. a+b+3*(a+b)-2*(a+b)
1200 ;; *PLUSFLAG* is set by PLUSIN to indicate that a mplus expression is
1201 ;; part of the result. For this case go again through the terms of the
1202 ;; result and add any term of the mplus expression into the list out.
1205 (when (mplusp (cadr fm
))
1207 (rplacd fm
(cddr fm
))
1211 (if (null (cdr fm
)) (return out
))
1214 ;;;-----------------------------------------------------------------------------
1215 ;;; PLUSIN (X FM) 27.09.2010/DK
1217 ;;; Arguments and values:
1218 ;;; X - a Maxima expression or atom
1219 ;;; FM - a list with the terms of an addition
1220 ;;; result - part of the list fm, which starts at the inserted expression
1223 ;;; Adds X into running list of additive terms FM. The routine modifies
1224 ;;; the argument FM destructively, but does not return the modified list as
1225 ;;; a result. The return value is a part of the list FM, which starts at the
1226 ;;; inserted term. PLUSIN can not handle Maxima numbers. PLUSIN is called
1227 ;;; only from the routine PLS.
1231 ;;; (plusin '$a fm) -> ($A)
1233 ;;; (plusin '$b fm) -> ($B)
1235 ;;; (plusin '$a fm) -> (((MTIMES SIMP) 2 $A) $B)
1236 ;;; fm -> (0 ((MTIMES SIMP) 2 $A) $B)
1239 ;;; Modifies destructively the argument FM, which contains the result of the
1240 ;;; addition of the argument X into the list FM.
1243 ;;; The option variables $doallmxops and $listarith.
1246 ;;; The return value is used in PLS to go in parallel through the list of
1247 ;;; terms, when adding a complete mplus-expression into the list of terms.
1248 ;;; This is triggered by the flag *PLUSFLAG*, which is set in PLUSIN, if
1249 ;;; a mplus-expression is added to the result list.
1250 ;;;-----------------------------------------------------------------------------
1252 (defun plusin (x fm
)
1253 (prog (x1 x2 flag check v w xnew a n m c
)
1258 (if (mnump (cadr x
)) (setq w
(cadr x
) x
(cddr x
))
1260 (t (setq x
(ncons x
))))
1261 (setq x1
(if (null (cdr x
)) (car x
) (cons '(mtimes) x
))
1262 xnew
(list* '(mtimes) w x
))
1264 (cond ((null (cdr fm
)))
1265 ((and (alike1 x1
(cadr fm
)) (null (cdr x
)))
1267 ;; Implement the simplification of
1268 ;; v*a^(c+n)+w*a^(c+m) -> (v*a^n+w*a^m)*a^c
1269 ;; where a, v, w, and (n-m) are integers.
1270 ((and (or (and (mexptp (setq x2
(cadr fm
)))
1273 (not (alike1 x1 x2
))
1275 (integerp (setq v
(cadr x2
)))
1276 (mexptp (setq x2
(caddr x2
)))))
1277 (integerp (setq a
(cadr x2
)))
1280 (integerp (sub (caddr x2
) (caddr x1
))))
1281 (setq n
(if (and (mplusp (caddr x2
))
1282 (mnump (cadr (caddr x2
))))
1284 (if (mnump (caddr x2
))
1287 (setq m
(if (and (mplusp (caddr x1
))
1288 (mnump (cadr (caddr x1
))))
1290 (if (mnump (caddr x1
))
1293 (setq c
(sub (caddr x2
) n
))
1296 ;; n and m are integers and the result is (v*a^n+w*a^m)*a^c.
1297 (setq x1
(mul (addk (timesk v
(exptb a n
))
1298 (timesk w
(exptb a m
)))
1302 ;; n and m are rational numbers: The difference n-m is an
1303 ;; integer. The rational numbers might be improper fractions.
1304 ;; The mixed numbers are: n = n1 + d1/r and m = n2 + d2/r,
1305 ;; where r is the common denominator. We have two cases:
1306 ;; I) d1 = d2: e.g. 2^(1/3+c)+2^(4/3+c)
1307 ;; The result is (v*a^n1+w*a^n2)*a^(c+d1/r)
1308 ;; II) d1 # d2: e.g. 2^(1/2+c)+2^(-1/2+c)
1309 ;; In this case one of the exponents d1 or d2 must
1310 ;; be negative. The negative exponent is factored out.
1311 ;; This guarantees that the factor (v*a^n1+w*a^n2)
1312 ;; is an integer. But the positive exponent has to be
1313 ;; adjusted accordingly. E.g. when we factor out
1314 ;; a^(d2/r) because d2 is negative, then we have to
1315 ;; adjust the positive exponent to n1 -> n1+(d1-d2)/r.
1317 ;; Part of the simplification is done in simptimes. E.g.
1318 ;; this algorithm simplifies the sum sqrt(2)+3*sqrt(2)
1319 ;; to 4*sqrt(2). In simptimes this is further simplified
1321 (multiple-value-bind (n1 d1
)
1322 (truncate (num1 n
) (denom1 n
))
1323 (multiple-value-bind (n2 d2
)
1324 (truncate (num1 m
) (denom1 m
))
1325 (cond ((equal d1 d2
)
1326 ;; Case I: -> (v*a^n1+w*a^n2)*a^(c+d1/r)
1328 (mul (addk (timesk v
(exptb a n1
))
1329 (timesk w
(exptb a n2
)))
1332 (div d1
(denom1 n
))))))
1335 ;; Case II:: d2 is negative, adjust n1.
1336 (setq n1
(add n1
(div (sub d1 d2
) (denom1 n
))))
1338 (mul (addk (timesk v
(exptb a n1
))
1339 (timesk w
(exptb a n2
)))
1342 (div d2
(denom1 n
))))))
1345 ;; Case II: d1 is negative, adjust n2.
1346 (setq n2
(add n2
(div (sub d2 d1
) (denom1 n
))))
1348 (mul (addk (timesk v
(exptb a n1
))
1349 (timesk w
(exptb a n2
)))
1352 (div d1
(denom1 n
))))))
1354 ;; This clause should never be reached.
1355 (t (merror "Internal error in simplus."))))))))
1356 ((mtimesp (cadr fm
))
1357 (cond ((alike1 x1
(cadr fm
))
1359 ((and (mnump (cadadr fm
)) (alike x
(cddadr fm
)))
1360 (setq flag t
) ; found common factor
1362 ((great xnew
(cadr fm
)) (go gr
))))
1363 ((great x1
(cadr fm
)) (go gr
)))
1364 (setq xnew
(eqtest (testt xnew
) (or check
'((foo)))))
1365 (return (cdr (rplacd fm
(cons xnew
(cdr fm
)))))
1372 (list* '(mtimes simp
) 0 x
)
1373 ;; Call muln to get a simplified product.
1374 (if (mtimesp (setq x1
(muln (cons (addk 1 w
) x
) t
)))
1378 (cond ((not (mtimesp (cadr fm
)))
1381 ;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
1382 (rplacd (cadr fm
) (cddadr fm
))
1384 ((not (zerop1 (cadadr fm
)))
1386 ;; Handle the multiplication with a zero.
1387 ((and (or (not $listarith
) (not $doallmxops
))
1388 (mxorlistp (caddr (cadr fm
))))
1390 (cons (constmx 0 (caddr (cadr fm
))) (cddr fm
))))))
1391 ;; (cadadr fm) is zero. If the first term of fm is a number,
1392 ;; add it to preserve the type.
1393 (when (mnump (car fm
))
1394 (rplaca fm
(addk (car fm
) (cadadr fm
))))
1395 (return (rplacd fm
(cddr fm
)))
1397 ;; Call muln to get a simplified product.
1398 (setq x1
(muln (cons (addk w
(if flag
(cadadr fm
) 1)) x
) t
))
1399 ;; Make a mplus expression to guarantee that x1 is added again into the sum
1400 (setq x1
(list '(mplus) x1
))
1404 (list* '(mtimes) x1 x
)
1405 (if (mtimesp x1
) (testtneg x1
) x1
)))
1406 (if (not (mtimesp (cadr fm
))) (go check
))
1407 (when (and (onep (cadadr fm
)) flag
(null (cdddr (cadr fm
))))
1408 ;; Do this simplification for an integer 1, not for 1.0 and 1.0b0
1409 (rplaca (cdr fm
) (caddr (cadr fm
))) (go check
))
1412 (if (mplusp (cadr fm
)) (setq *plusflag
* t
)) ; A nested mplus expression
1415 ;;;-----------------------------------------------------------------------------
1417 ;; Routines to add matrices
1419 (defun mxplusc (sc mx
)
1421 (setq sc
(partition-ns (cdr sc
)))
1422 (cond ((null (car sc
)) (cons '(mplus) (cons mx
(cadr sc
))))
1423 ((not (null (cadr sc
)))
1426 (outermap1 'mplus
(cons '(mplus) (car sc
)) mx
))
1428 (t (simplify (outermap1 'mplus
(cons '(mplus) (car sc
)) mx
)))))
1429 ((not (scalar-or-constant-p sc $assumescalar
))
1430 (testp (pls mx
(pls sc nil
))))
1431 (t (simplify (outermap1 'mplus sc mx
)))))
1433 (defun partition-ns (x)
1434 (let (sp nsp
) ; SP = scalar part, NSP = nonscalar part
1435 (mapc #'(lambda (z) (if (scalar-or-constant-p z $assumescalar
)
1436 (setq sp
(cons z sp
))
1437 (setq nsp
(cons z nsp
))))
1439 (list (nreverse sp
) (nreverse nsp
))))
1441 (defun addmx (x1 x2
)
1442 (let (($doscmxops t
) ($domxmxops t
) ($listarith t
))
1443 (simplify (fmapl1 'mplus x1 x2
))))
1445 ;;; ----------------------------------------------------------------------------
1447 ;;; Simplification of the Log function
1449 ;; The log function distributes over lists, matrices, and equations
1450 (defprop %log
(mlist $matrix mequal
) distribute_over
)
1452 (defun simpln (x y z
)
1454 (setq y
(simpcheck (cadr x
) z
))
1455 (cond ((onep1 y
) (addk -
1 y
))
1457 (cond (radcanp (list '(%log simp
) 0))
1459 (merror (intl:gettext
"log: encountered log(0).")))
1460 (t (throw 'errorsw t
))))
1461 ;; Check evaluation in floating point precision.
1462 ((flonum-eval (mop x
) y
))
1463 ;; Check evaluation in bigfloag precision.
1464 ((and (not (member 'simp
(car x
)))
1465 (big-float-eval (mop x
) y
)))
1468 (cond ((or (and $logexpand
(eq $domain
'$real
))
1469 (member $logexpand
'($all $super
))
1470 (and (eq ($csign
(cadr y
)) '$pos
)
1471 (not (member ($csign
(caddr y
))
1472 '($complex $imaginary
)))))
1473 ;; Simplify log(x^a) -> a*log(x), where x > 0 and a is real
1474 (mul (caddr y
) (take '(%log
) (cadr y
))))
1475 ((or (and (ratnump (caddr y
))
1476 (or (eql 1 (cadr (caddr y
)))
1477 (eql -
1 (cadr (caddr y
)))))
1478 (maxima-integerp (inv (caddr y
))))
1479 ;; Simplify log(z^(1/n)) -> log(z)/n, where n is an integer
1481 (take '(%log
) (cadr y
))))
1482 ((and (eq (cadr y
) '$%e
)
1483 (or (not (member ($csign
(caddr y
))
1484 '($complex $imaginary
)))
1485 (not (member ($csign
(mul '$%i
(caddr y
)))
1486 '($complex $imaginary
)))))
1487 ;; Simplify log(exp(x)) and log(exp(%i*x)), where x is a real
1489 (t (eqtest (list '(%log
) y
) x
))))
1491 ;; Simplify log(n/d)
1492 (cond ((eql (cadr y
) 1)
1493 (mul -
1 (take '(%log
) (caddr y
))))
1494 ((eq $logexpand
'$super
)
1495 (sub (take '(%log
) (cadr y
)) (take '(%log
) (caddr y
))))
1496 (t (eqtest (list '(%log
) y
) x
))))
1497 ((and (member $logexpand
'($all $super
))
1499 (do ((y (cdr y
) (cdr y
))
1501 ((null y
) (return (addn b t
)))
1502 (setq b
(cons (take '(%log
) (car y
)) b
))))
1503 ((and (member $logexpand
'($all $super
))
1505 (member (caar y
) '(%product $product
)))
1506 (let ((new-op (if (char= (get-first-char (caar y
)) #\%
) '%sum
'$sum
)))
1507 (simplifya `((,new-op
) ((%log
) ,(cadr y
)) ,@(cddr y
)) t
)))
1510 (eq ($sign y
) '$neg
))
1511 (add (mul '$%i
'$%pi
) (take '(%log
) (neg y
))))
1512 ((taylorize (mop x
) (second x
)))
1513 (t (eqtest (list '(%log
) y
) x
))))
1516 (simplifya (list '(mtimes) (caddr w
)
1517 (simplifya (list '(%log
) (cadr w
)) t
)) t
))
1519 ;;; ----------------------------------------------------------------------------
1521 ;;; Implementation of the Square root function
1523 (defprop $sqrt %sqrt verb
)
1524 (defprop $sqrt %sqrt alias
)
1526 (defprop %sqrt $sqrt noun
)
1527 (defprop %sqrt $sqrt reversealias
)
1529 (defprop %sqrt simp-sqrt operators
)
1532 (simplify (list '(%sqrt
) z
)))
1534 (defun simp-sqrt (x ignored z
)
1535 (declare (ignore ignored
))
1537 (simplifya (list '(mexpt) (cadr x
) '((rat simp
) 1 2)) z
))
1539 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1541 ;;; Simplification of the "/" operator.
1543 (defun simpquot (x y z
)
1545 (cond ((and (integerp (cadr x
)) (integerp (caddr x
)) (not (zerop (caddr x
))))
1546 (*red
(cadr x
) (caddr x
)))
1547 ((and (numberp (cadr x
)) (numberp (caddr x
)) (not (zerop (caddr x
))))
1548 (/ (cadr x
) (caddr x
)))
1549 ((and (floatp (cadr x
)) (floatp (caddr x
)) #-ieee-floating-point
(not (zerop (caddr x
))))
1550 (/ (cadr x
) (caddr x
)))
1551 ((and ($bfloatp
(cadr x
)) ($bfloatp
(caddr x
)) (not (equal bigfloatzero
(caddr x
))))
1552 ;; Call BIGFLOATP to ensure that arguments have same precision.
1553 ;; Otherwise FPQUOTIENT could return a spurious value.
1554 (bcons (fpquotient (cdr (bigfloatp (cadr x
))) (cdr (bigfloatp (caddr x
))))))
1555 (t (setq y
(simplifya (cadr x
) z
))
1556 (setq x
(simplifya (list '(mexpt) (caddr x
) -
1) z
))
1557 (if (equal y
1) x
(simplifya (list '(mtimes) y x
) t
)))))
1559 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1561 ;;; Implementation of the abs function.
1563 ;; Put the properties alias, reversealiases, noun and verb on the property list.
1564 (defprop $abs mabs alias
)
1565 (defprop $abs mabs verb
)
1566 (defprop mabs $abs reversealias
)
1567 (defprop mabs $abs noun
)
1569 ;; The abs function distributes over bags.
1570 (defprop mabs
(mlist $matrix mequal
) distribute_over
)
1572 ;; Define a verb function $abs
1574 (simplify (list '(mabs) x
)))
1576 ;; The abs function is a simplifying function.
1577 (defprop mabs simpabs operators
)
1579 (defun simpabs (e y z
)
1580 (declare (ignore y
))
1583 (x (simpcheck (second e
) z
)))
1585 (cond ((complex-number-p x
#'(lambda (s) (or (floatp s
) ($bfloatp s
))))
1586 (maxima::to
(bigfloat::abs
(bigfloat:to x
))))
1588 ((complex-number-p x
#'mnump
)
1591 ;; nounform for arrays...
1592 ((or (arrayp x
) ($member x $arrays
)) `((mabs simp
) ,x
))
1594 ;; taylor polynomials
1595 ((taylorize 'mabs x
))
1597 ;; values for extended real arguments:
1598 ((member x
'($inf $infinity $minf
)) '$inf
)
1599 ((member x
'($ind $und
)) x
)
1601 ;; abs(abs(expr)) --> abs(expr). Since x is simplified, it's OK to return x.
1602 ((and (consp x
) (consp (car x
)) (eq (caar x
) 'mabs
))
1605 ;; abs(conjugate(expr)) = abs(expr).
1606 ((and (consp x
) (consp (car x
)) (eq (caar x
) '$conjugate
))
1607 (take '(mabs) (cadr x
)))
1610 (setq sgn
($csign x
))
1611 (cond ((member sgn
'($neg $nz
)) (mul -
1 x
))
1612 ((eq '$zero sgn
) (mul 0 x
))
1613 ((member sgn
'($pos $pz
)) x
)
1615 ;; for complex constant expressions, use $cabs
1616 ((and (eq sgn
'$complex
) ($constantp x
))
1619 ;; abs(pos^complex) --> pos^(realpart(complex)).
1620 ((and (eq sgn
'$complex
) (mexptp x
) (eq '$pos
($csign
(second x
))))
1621 (power (second x
) ($realpart
(third x
))))
1623 ;; for abs(neg^z), use cabs.
1624 ((and (mexptp x
) (eq '$neg
($csign
(second x
))))
1627 ;; When x # 0, we have abs(signum(x)) = 1.
1628 ((and (eq '$pn sgn
) (consp x
) (consp (car x
)) (eq (caar x
) '%signum
)) 1)
1630 ;; multiplicative property: abs(x*y) = abs(x) * abs(y). We would like
1631 ;; assume(a*b > 0), abs(a*b) --> a*b. Thus the multiplicative property
1632 ;; is applied after the sign test.
1634 (muln (mapcar #'(lambda (u) (take '(mabs) u
)) (margs x
)) t
))
1636 ;; abs(x^n) = abs(x)^n for integer n. Is the featurep check worthwhile?
1637 ;; Again the sign check is done first because we'd like abs(x^2) --> x^2.
1638 ((and (mexptp x
) ($featurep
(caddr x
) '$integer
))
1639 (power (take '(mabs) (cadr x
)) (caddr x
)))
1641 ;; Reflection rule: abs(-x) --> abs(x)
1642 ((great (neg x
) x
) (take '(mabs) (neg x
)))
1645 (t (eqtest (list '(mabs) x
) e
)))))))
1647 (defun abs-integral (x)
1648 (mul (div 1 2) x
(take '(mabs) x
)))
1650 (putprop 'mabs
`((x) ,#'abs-integral
) 'integral
)
1652 ;; I (rtoy) think this does some simple optimizations of x * y.
1657 ;; We have something like ((mtimes) foo). This is the same as foo.
1660 ;; We have 1*foo. Which is the same as foo. This should not
1661 ;; be applied to 1.0 or 1b0!
1662 (cond ((null (cdddr x
))
1664 (t (rplacd x
(cddr x
)))))
1668 ;; This basically converts -(a+b) to -a-b.
1670 (cond ((and (equal (cadr x
) -
1)
1674 ;; If x is exactly of the form -1*(sum), and $negdistrib is
1675 ;; true, we distribute the -1 across the sum.
1676 (addn (mapcar #'(lambda (z)
1682 ;; Simplification of the "-" operator
1683 (defun simpmin (x vestigial z
)
1684 (declare (ignore vestigial
))
1685 (cond ((null (cdr x
)) 0)
1687 (mul -
1 (simplifya (cadr x
) z
)))
1689 ;; ((mminus) a b ...) -> ((mplus) a ((mtimes) -1 b) ...)
1690 (sub (simplifya (cadr x
) z
) (addn (cddr x
) z
)))))
1692 (defun simptimes (x w z
) ; W must be 1
1693 (prog (res check eqnflag matrixflag sumflag
)
1694 (if (null (cdr x
)) (return 1))
1700 (cond ((and matrixflag
(mxorlistp1 matrixflag
))
1701 (return (constmx res matrixflag
)))
1702 (eqnflag (return (list '(mequal simp
)
1703 (mul2 res
(cadr eqnflag
))
1704 (mul2 res
(caddr eqnflag
)))))
1707 (cond ((mxorlistp u
)
1708 (return (setq res
(constmx res u
))))
1710 (mxorlistp1 (cadr u
))
1711 ($numberp
(caddr u
)))
1712 (return (setq res
(constmx res
(cadr u
)))))
1716 (list '(mequal simp
)
1718 (mul2 res
(caddr u
))))))))))))
1720 ((null x
) (go end
)))
1721 (setq w
(if z
(car x
) (simplifya (car x
) nil
)))
1723 (cond ((atom w
) nil
)
1724 ((eq (caar w
) 'mrat
)
1725 (cond ((or eqnflag matrixflag
1727 (not (member 'trunc
(cdar w
))))
1728 (spsimpcases (cdr x
) w
))
1729 (setq w
(ratdisrep w
))
1733 (ratf (cons '(mtimes)
1734 (nconc (mapcar #'simplify
(cons w
(cdr x
)))
1736 ((eq (caar w
) 'mequal
)
1741 (mul2 (cadr eqnflag
) (cadr w
))
1742 (mul2 (caddr eqnflag
) (caddr w
)))))
1744 ((member (caar w
) '(mlist $matrix
))
1746 (cond ((not matrixflag
) w
)
1747 ((and (or $doallmxops $domxmxops $domxtimes
)
1748 (or (not (eq (caar w
) 'mlist
)) $listarith
)
1749 (not (eq *inv
* '$detout
)))
1750 (stimex matrixflag w
))
1751 (t (setq res
(tms (copy-tree w
) 1 (copy-tree res
))) matrixflag
)))
1753 ((and (eq (caar w
) '%sum
) $sumexpand
)
1754 (setq sumflag
(sumtimes sumflag w
))
1756 (setq res
(tms (copy-tree w
) 1 (copy-tree res
)))
1759 (cond ((mtimesp res
) (setq res
(testt res
))))
1760 (cond (sumflag (setq res
(cond ((or (null res
) (equal res
1)) sumflag
)
1761 ((not (mtimesp res
))
1762 (list '(mtimes) res sumflag
))
1763 (t (nconc res
(list sumflag
)))))))
1764 (cond ((or (atom res
)
1765 (not (member (caar res
) '(mexpt mtimes
)))
1766 (and (zerop $expop
) (zerop $expon
))
1768 ((eq (caar res
) 'mtimes
) (setq res
(expandtimes res
)))
1769 ((and (mplusp (cadr res
))
1770 (fixnump (caddr res
))
1771 (not (or (> (caddr res
) $expop
)
1772 (> (- (caddr res
)) $expon
))))
1773 (setq res
(expandexpt (cadr res
) (caddr res
)))))
1776 (cond ((null res
) matrixflag
)
1777 ((and (or ($listp matrixflag
)
1780 (not (member res
'(-1 -
1.0))))
1781 ;; RES should only be -1 here (not = 1)
1783 (member res
'(-1 -
1.0))))
1784 (or (not ($listp matrixflag
)) $listarith
))
1785 (mxtimesc res matrixflag
))
1786 (t (testt (tms matrixflag
1 (tms res
1 nil
))))))))
1787 (if res
(setq res
(eqtest res check
)))
1788 (return (cond (eqnflag
1789 (if (null res
) (setq res
1))
1791 (mul2 (cadr eqnflag
) res
)
1792 (mul2 (caddr eqnflag
) res
)))
1795 (defun spsimpcases (l e
)
1797 (if (or (mbagp u
) (and (not (atom u
))
1799 (not (member 'trunc
(cdar e
)))))
1802 (defun mxtimesc (sc mx
)
1804 (and (mtimesp sc
) (member (cadr sc
) '(-1 -
1.0))
1805 $doscmxops
(not (or $doallmxops $domxmxops $domxtimes
))
1806 (setq sign
(cadr sc
)) (rplaca (cdr sc
) nil
))
1807 (setq out
(let ((scp* (cond ((mtimesp sc
) (partition-ns (cdr sc
)))
1808 ((not (scalar-or-constant-p sc $assumescalar
))
1811 (cond ((null scp
*) (list '(mtimes simp
) sc mx
))
1812 ((and (not (atom scp
*)) (null (car scp
*)))
1813 (append '((mtimes)) (cadr scp
*) (list mx
)))
1814 ((or (atom scp
*) (and (null (cdr scp
*))
1815 (not (null (cdr sc
)))
1816 (setq scp
* (cons '(mtimes) (car scp
*))))
1818 (simplifya (outermap1 'mtimes scp
* mx
) nil
))
1819 (t (append '((mtimes))
1822 (cons '(mtimes) (car scp
*)) mx
)
1825 (cond (sign (if (mtimesp out
)
1826 (rplacd out
(cons sign
(cdr out
)))
1827 (list '(mtimes) sign out
)))
1828 ((mtimesp out
) (testt out
))
1832 (let (($doscmxops t
) ($domxmxops t
) ($listarith t
))
1833 (simplify (fmapl1 'mtimes x y
))))
1835 ;; TMS takes a simplified expression FACTOR and a cumulative
1836 ;; PRODUCT as arguments and modifies the cumulative product so
1837 ;; that the expression is now one of its factors. The
1838 ;; exception to this occurs when a tellsimp rule is triggered.
1839 ;; The second argument is the POWER to which the expression is
1840 ;; to be raised within the product.
1842 (defun tms (factor power product
&aux tem
)
1845 (when (mplusp product
) (setq product
(list '(mtimes simp
) product
)))
1846 (cond ((zerop1 factor
)
1847 (cond ((mnegp power
)
1850 (merror (intl:gettext
"Division by 0"))))
1852 ((and (null product
)
1853 (or (and (mtimesp factor
) (equal power
1))
1854 (and (setq product
(list '(mtimes) 1)) nil
)))
1855 (setq tem
(append '((mtimes)) (if (mnump (cadr factor
)) nil
'(1))
1857 (if (= (length tem
) 1)
1858 (setq tem
(copy-list tem
))
1861 (do ((factor-list (cdr factor
) (cdr factor-list
)))
1862 ((or (null factor-list
) (zerop1 product
)) product
)
1863 (setq z
(timesin (car factor-list
) (cdr product
) power
))
1866 (setq product
(tms-format-product z
)))))
1868 (setq z
(timesin factor
(cdr product
) power
))
1870 (tms-format-product z
)
1873 (defun tms-format-product (x)
1874 (cond ((zerop1 x
) x
)
1875 ((mnump x
) (list '(mtimes) x
))
1876 ((not (mtimesp x
)) (list '(mtimes) 1 x
))
1877 ((not (mnump (cadr x
))) (cons '(mtimes) (cons 1 (cdr x
))))
1881 (cond ($ratsimpexpons
(sratsimp (list '(mplus) x y
)))
1882 ((and (mnump x
) (mnump y
)) (addk x y
))
1886 (if (and (mnump x
) (mnump y
))
1890 (defun simp-limit (x vestigial z
)
1891 (declare (ignore vestigial
))
1892 (let ((l1 (length x
))
1894 (unless (or (= l1
2) (= l1
4) (= l1
5))
1895 (merror (intl:gettext
"limit: wrong number of arguments.")))
1896 (setq y
(simpmap (cdr x
) z
))
1897 (cond ((and (= l1
5) (not (member (cadddr y
) '($plus $minus
))))
1898 (merror (intl:gettext
"limit: direction must be either 'plus' or 'minus': ~M") (cadddr y
)))
1900 (merror (intl:gettext
"limit: variable must not be a number; found: ~M") (cadr y
)))
1904 (eqtest (cons '(%limit
) y
) x
)))))
1906 (defun simpinteg (x vestigial z
)
1907 (declare (ignore vestigial
))
1908 (let ((l1 (length x
))
1910 (unless (or (= l1
3) (= l1
5))
1911 (merror (intl:gettext
"integrate: wrong number of arguments.")))
1912 (setq y
(simpmap (cdr x
) z
))
1913 (cond ((mnump (cadr y
))
1914 (merror (intl:gettext
"integrate: variable must not be a number; found: ~M") (cadr y
)))
1915 ((and (= l1
5) (alike1 (caddr y
) (cadddr y
)))
1918 (free (setq z
(sub (cadddr y
) (caddr y
))) '$%i
)
1919 (eq ($sign z
) '$neg
))
1920 (neg (simplifya (list '(%integrate
) (car y
) (cadr y
) (cadddr y
) (caddr y
)) t
)))
1924 (if (or (among '$inf z
) (among '$minf z
))
1928 (eqtest (cons '(%integrate
) y
) x
)))))
1930 (defun simpbigfloat (x vestigial simp-flag
)
1931 (declare (ignore vestigial simp-flag
))
1934 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1936 ;;; Implementation of the Exp function.
1938 (defprop $exp %exp verb
)
1939 (defprop $exp %exp alias
)
1941 (defprop %exp $exp noun
)
1942 (defprop %exp $exp reversealias
)
1944 (defprop %exp simp-exp operators
)
1947 (simplify (list '(%exp
) z
)))
1949 ;; Support a function for code,
1950 ;; which depends on an unsimplified noun form.
1951 (defmfun $exp-form
(z)
1952 (list '(mexpt) '$%e z
))
1954 (defun simp-exp (x ignored z
)
1955 (declare (ignore ignored
))
1957 (simplifya (list '(mexpt) '$%e
(cadr x
)) z
))
1959 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
1961 (defun simplambda (x vestigial simp-flag
)
1962 (declare (ignore vestigial simp-flag
))
1963 ; Check for malformed lambda expressions.
1964 ; We verify that we have a valid list of parameters and a non-empty body.
1965 (let ((params (cadr x
)))
1966 (unless ($listp params
)
1967 (merror (intl:gettext
"lambda: first argument must be a list; found: ~M") params
))
1968 (do ((params (cdr params
) (cdr params
))
1971 (when (mdeflistp params
)
1972 (setq params
(cdar params
)))
1973 (let ((p (car params
)))
1974 (unless (or (mdefparam p
)
1975 (and (op-equalp p
'mquote
)
1976 (mdefparam (cadr p
))))
1977 (merror (intl:gettext
"lambda: parameter must be a symbol and must not be a system constant; found: ~M") p
))
1979 (when (member p seen-params
:test
#'eq
)
1980 (merror (intl:gettext
"lambda: ~M occurs more than once in the parameter list") p
))
1981 (push p seen-params
))))
1982 (when (null (cddr x
))
1983 (merror (intl:gettext
"lambda: no body present")))
1984 (cons '(lambda simp
) (cdr x
)))
1986 (defun simpmdef (x vestigial simp-flag
)
1987 (declare (ignore vestigial simp-flag
))
1989 (cons '(mdefine simp
) (cdr x
)))
1991 (defun simpmap (e z
)
1992 (mapcar #'(lambda (u) (simpcheck u z
)) e
))
1995 (let ((x ($expand e
1 1)))
1996 (cond ((or (not (free x
'$ind
)) (not (free x
'$und
))
1997 (not (free x
'$zeroa
)) (not (free x
'$zerob
))
1998 (not (free x
'$infinity
))
2001 ((and (free x
'$inf
) (free x
'$minf
)) x
)
2002 (t (infsimp1 x e
)))))
2004 (defun infsimp1 (x e
)
2005 (let ((minf-coef (coeff x
'$minf
1))
2006 (inf-coef (coeff x
'$inf
1)))
2007 (cond ((or (and (equal minf-coef
0)
2009 (and (not (free minf-coef
'$inf
))
2010 (not (free inf-coef
'$minf
)))
2011 (let ((new-exp (sub (add2 (mul2 minf-coef
'$minf
)
2012 (mul2 inf-coef
'$inf
))
2014 (and (not (free new-exp
'$inf
))
2015 (not (free new-exp
'$minf
)))))
2017 (t (let ((sign-minf-coef ($asksign minf-coef
))
2018 (sign-inf-coef ($asksign inf-coef
)))
2019 (cond ((or (and (eq sign-inf-coef
'$zero
)
2020 (eq sign-minf-coef
'$neg
))
2021 (and (eq sign-inf-coef
'$pos
)
2022 (eq sign-minf-coef
'$zero
))
2023 (and (eq sign-inf-coef
'$pos
)
2024 (eq sign-minf-coef
'$neg
))) '$inf
)
2025 ((or (and (eq sign-inf-coef
'$zero
)
2026 (eq sign-minf-coef
'$pos
))
2027 (and (eq sign-inf-coef
'$neg
)
2028 (eq sign-minf-coef
'$zero
))
2029 (and (eq sign-inf-coef
'$neg
)
2030 (eq sign-minf-coef
'$pos
))) '$minf
)
2031 ((or (and (eq sign-inf-coef
'$pos
)
2032 (eq sign-minf-coef
'$pos
))
2033 (and (eq sign-inf-coef
'$neg
)
2034 (eq sign-minf-coef
'$neg
))) '$und
)))))))
2036 (defun infsimp2 (x e
)
2038 (if (isinop x
'%limit
) e x
))
2040 (defun simpderiv (x y z
)
2042 (cond ((not (even (length x
)))
2043 (cond ((and (cdr x
) (null (cdddr x
))) (nconc x
'(1)))
2044 (t (wna-err '%derivative
)))))
2045 (setq w
(cons '(%derivative
) (simpmap (cdr x
) z
)))
2047 (do ((u (cddr w
) (cddr u
))) ((null u
))
2048 (cond ((mnump (car u
))
2049 (merror (intl:gettext
"diff: variable must not be a number; found: ~M") (car u
)))))
2050 (cond ((or (zerop1 y
)
2051 (and (or (mnump y
) (and (atom y
) (constant y
)))
2053 (and (not (alike1 y
(caddr w
)))
2054 (do ((u (cddr w
) (cddr u
))) ((null u
))
2055 (cond ((and (numberp (cadr u
))
2056 (not (zerop (cadr u
))))
2059 ((and (not (atom y
)) (eq (caar y
) '%derivative
) derivsimp
)
2060 (rplacd w
(append (cdr y
) (cddr w
)))))
2062 (return (if (null derivflag
) (list '(%del simp
) y
) (deriv (cdr w
)))))
2065 (cond ((null u
) (go next
))
2066 ((zerop1 (caddr u
)) (rplacd u
(cdddr u
)))
2067 (t (setq u
(cddr u
))))
2070 (cond ((null (cddr w
)) (return y
))
2071 ((and (null (cddddr w
))
2073 (alike1 (cadr w
) (caddr w
)))
2078 (cond ((null (cddr z
)) (go loop
))
2079 ((alike1 (car z
) (caddr z
))
2080 (rplaca (cdddr z
) (add2 (cadr z
) (cadddr z
)))
2081 (rplacd z
(cdddr z
)))
2082 ((great (car z
) (caddr z
))
2083 (let ((u1 (car z
)) (u2 (cadr z
)) (v1 (caddr z
)) (v2 (cadddr z
)))
2086 (rplacd z
(cons v2
(cons u1
(cons u2
(cddddr z
))))))))
2087 (cond ((setq z
(cddr z
)) (go sort
)))
2089 (cond ((null flag
) (return (cond ((null derivflag
) (eqtest w x
))
2090 (t (deriv (cdr w
)))))))
2096 (setq x
(num1 x
)) (cond ((plusp x
) 1) ((minusp x
) -
1) (t 0)))
2098 ((mplusp x
) (if expandp
1 (signum1 (car (last x
)))))
2099 ((mtimesp x
) (if (mplusp (cadr x
)) 1 (signum1 (cadr x
))))
2102 (defprop %signum
(mlist $matrix mequal
) distribute_over
)
2104 (defun simpsignum (e y z
)
2105 (declare (ignore y
))
2107 (let ((x (simpcheck (second e
) z
)) (sgn))
2109 (cond ((complex-number-p x
#'mnump
)
2110 (if (complex-number-p x
#'$ratnump
) ;; nonfloat complex
2111 (if (zerop1 x
) 0 ($rectform
(div x
($cabs x
))))
2112 (maxima::to
(bigfloat::signum
(bigfloat::to x
)))))
2114 ;; idempotent: signum(signum(z)) = signum(z).
2115 ((and (consp x
) (consp (car x
)) (eq '%signum
(mop x
))) x
)
2118 (setq sgn
($csign x
))
2119 (cond ((eq sgn
'$neg
) -
1)
2123 ;; multiplicative: signum(ab) = signum(a) * signum(b).
2125 (muln (mapcar #'(lambda (s) (take '(%signum
) s
)) (margs x
)) t
))
2127 ;; Reflection rule: signum(-x) --> -signum(x).
2128 ((great (neg x
) x
) (neg (take '(%signum
) (neg x
))))
2131 (t (eqtest (list '(%signum
) x
) e
)))))))
2133 (defun exptrl (r1 r2
)
2134 (cond ((equal r2
1) r1
)
2136 (cond ((mnump r1
) (addk 0.0 r1
))
2137 ;; Do not simplify the type of the number away.
2138 (t (list '(mexpt simp
) r1
1.0))))
2139 ((equal r2 bigfloatone
)
2140 (cond ((mnump r1
) ($bfloat r1
))
2141 ;; Do not simplify the type of the number away.
2142 (t (list '(mexpt simp
) r1 bigfloatone
))))
2144 (cond ((or (zerop1 r2
) (mnegp r2
))
2146 (merror (intl:gettext
"expt: undefined: ~M") (list '(mexpt) r1 r2
))
2147 (throw 'errorsw t
)))
2148 (t (zerores r1 r2
))))
2149 ((or (zerop1 r2
) (onep1 r1
))
2150 (cond ((or ($bfloatp r1
) ($bfloatp r2
)) bigfloatone
)
2151 ((or (floatp r1
) (floatp r2
)) 1.0)
2153 ((or ($bfloatp r1
) ($bfloatp r2
)) ($bfloat
(list '(mexpt) r1 r2
)))
2154 ((and (numberp r1
) (integerp r2
)) (exptb r1 r2
))
2155 ((and (numberp r1
) (floatp r2
) (equal r2
(float (floor r2
))))
2156 (exptb (float r1
) (floor r2
)))
2157 ((or $numer
(and (floatp r2
) (or (plusp (num1 r1
)) $numer_pbranch
)))
2158 (let (y #+kcl
(r1 r1
) #+kcl
(r2 r2
))
2159 (cond ((minusp (setq r1
(addk 0.0 r1
)))
2160 (cond ((or $numer_pbranch
(eq $domain
'$complex
))
2162 ;; R1^R2 = (-R1)^R2*cos(pi*R2) + i*(-R1)^R2*sin(pi*R2)
2163 (setq r2
(addk 0.0 r2
))
2164 (setq y
(exptrl (- r1
) r2
) r2
(* %pi-val r2
))
2165 (add2 (* y
(cos r2
))
2166 (list '(mtimes simp
) (* y
(sin r2
)) '$%i
)))
2167 (t (setq y
(let ($numer $float $keepfloat $ratprint
)
2169 (mul2 y
(exptrl (- r1
) r2
)))))
2170 ((equal (setq r2
(addk 0.0 r2
)) (float (floor r2
)))
2171 (exptb r1
(floor r2
)))
2172 ((and (equal (setq y
(* 2.0 r2
)) (float (floor y
)))
2173 (not (equal r1 %e-val
)))
2174 (exptb (sqrt r1
) (floor y
)))
2175 (t (exp (* r2
(log r1
)))))))
2176 ((floatp r2
) (list '(mexpt simp
) r1 r2
))
2179 (exptrl (cond ((equal (abs (cadr r1
)) 1)
2180 (* (cadr r1
) (caddr r1
)))
2181 ;; We set the simp flag at this place. This
2182 ;; changes nothing for an exponent r2 # -1.
2183 ;; exptrl is called again and does not look at
2184 ;; the simp flag. For the case r2 = -1 exptrl
2185 ;; is called with an exponent 1. For this case
2186 ;; the base is immediately returned. Now the
2187 ;; base has the correct simp flag. (DK 02/2010)
2189 (list '(rat simp
) (- (caddr r1
)) (- (cadr r1
))))
2190 (t (list '(rat simp
) (caddr r1
) (cadr r1
))))
2192 (t (list '(rat simp
) (exptb (cadr r1
) r2
) (exptb (caddr r1
) r2
)))))
2193 ((and (floatp r1
) (alike1 r2
'((rat) 1 2)))
2195 (list '(mtimes simp
) (sqrt (- r1
)) '$%i
)
2197 ((and (floatp r1
) (alike1 r2
'((rat) -
1 2)))
2199 (list '(mtimes simp
) (/ -
1.0 (sqrt (- r1
))) '$%i
)
2203 (exptrl r1
(fpcofrat r2
))
2204 (mul2 (exptrl -
1 r2
) ;; (-4.5)^(1/4) -> (4.5)^(1/4) * (-1)^(1/4)
2205 (exptrl (- r1
) r2
))))
2206 (exptrlsw (list '(mexpt simp
) r1 r2
))
2209 (simptimes (list '(mtimes)
2210 (exptrl r1
(truncate (cadr r2
) (caddr r2
)))
2211 (let ((y (let ($keepfloat $ratprint
)
2212 (simpnrt r1
(caddr r2
))))
2213 (z (rem (cadr r2
) (caddr r2
))))
2215 (list (car y
) (cadr y
) (mul2 (caddr y
) z
))
2219 (defun simpexpt (x y z
)
2220 (prog (gr pot check res rulesw w mlpgr mlppot
)
2222 (cond (z (setq gr
(cadr x
) pot
(caddr x
)) (go cont
)))
2224 (setq gr
(simplifya (cadr x
) nil
))
2226 (let (($%enumer $numer
))
2227 ;; Switch $%enumer on, when $numer is TRUE to allow
2228 ;; simplification of $%e to its numerical value.
2229 (simplifya (if $ratsimpexpons
($ratsimp
(caddr x
)) (caddr x
))
2233 (setq pot
(ratdisrep pot
))
2236 (cond ((member 'trunc
(car gr
))
2237 (return (srf (list '(mexpt) gr pot
))))
2239 (let ((varlist (caddar gr
)) (genvar (cadddr (car gr
))))
2240 (return (ratrep* (list '(mexpt) gr pot
)))))
2242 (setq gr
(ratdisrep gr
))
2244 ((or (setq mlpgr
(mxorlistp gr
))
2245 (setq mlppot
(mxorlistp pot
)))
2247 ((onep1 pot
) (go atgr
))
2248 ((or (zerop1 pot
) (onep1 gr
)) (go retno
))
2250 ;; This code tries to handle 0^a more complete.
2251 ;; If the sign of realpart(a) is not known return an unsimplified
2252 ;; expression. The handling of the flag *zexptsimp? is not changed.
2253 ;; Reverting the return of an unsimplified 0^a, because timesin
2254 ;; can not handle such expressions. (DK 02/2010)
2256 (cond ((or (member (setq z
($csign pot
)) '($neg $nz
))
2257 (and *zexptsimp?
(eq ($asksign pot
) '$neg
)))
2258 ;; A negative exponent. Maxima error.
2259 (cond ((not errorsw
) (merror (intl:gettext
"expt: undefined: 0 to a negative exponent.")))
2260 (t (throw 'errorsw t
))))
2261 ((and (member z
'($complex $imaginary
))
2262 ;; A complex exponent. Look at the sign of the realpart.
2263 (member (setq z
($sign
($realpart pot
)))
2265 (cond ((not errorsw
)
2266 (merror (intl:gettext
"expt: undefined: 0 to a complex exponent.")))
2267 (t (throw 'errorsw t
))))
2268 ((and *zexptsimp?
(eq ($asksign pot
) '$zero
))
2269 (cond ((not errorsw
)
2270 (merror (intl:gettext
"expt: undefined: 0^0")))
2271 (t (throw 'errorsw t
))))
2272 ((not (member z
'($pos $pz
)))
2273 ;; The sign of realpart(pot) is not known. We can not return
2274 ;; an unsimplified 0^a expression, because timesin can not
2275 ;; handle it. We return ZERO. That is the old behavior.
2276 ;; Look for the imaginary symbol to be consistent with
2278 (cond ((not (free pot
'$%i
))
2279 (cond ((not errorsw
)
2280 (merror (intl:gettext
"expt: undefined: 0 to a complex exponent.")))
2281 (t (throw 'errorsw t
))))
2283 ;; Return ZERO and not an unsimplified expression.
2284 (return (zerores gr pot
)))))
2285 (t (return (zerores gr pot
)))))
2289 (or (not (ratnump gr
)) (not (ratnump pot
))))
2290 (return (eqtest (exptrl gr pot
) check
)))
2291 ;; Check for numerical evaluation of the sqrt.
2292 ((and (alike1 pot
'((rat) 1 2))
2293 (or (setq res
(flonum-eval '%sqrt gr
))
2294 (and (not (member 'simp
(car x
)))
2295 (setq res
(big-float-eval '%sqrt gr
)))))
2298 (return (%itopot pot
)))
2299 ((and (realp gr
) (minusp gr
) (mevenp pot
))
2302 ((and (realp gr
) (minusp gr
) (moddp pot
))
2303 (return (mul2 -
1 (power (- gr
) pot
))))
2304 ((and (equal gr -
1) (maxima-integerp pot
) (mminusp pot
))
2305 (setq pot
(neg pot
))
2308 (maxima-integerp pot
)
2311 (integerp (cadr pot
))
2313 (maxima-integerp (caddr pot
)))
2314 (setq pot
(caddr pot
))
2316 ((atom gr
) (go atgr
))
2317 ((and (eq (caar gr
) 'mabs
)
2320 (or (and (eq $domain
'$real
) (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
))))
2321 (and (eq $domain
'$complex
) (apparently-real-to-judge-by-$csign-p
(cadr gr
)))))
2322 (return (power (cadr gr
) pot
)))
2323 ((and (eq (caar gr
) 'mabs
)
2326 (not (equal pot -
1))
2327 (or (and (eq $domain
'$real
) (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
))))
2328 (and (eq $domain
'$complex
) (apparently-real-to-judge-by-$csign-p
(cadr gr
)))))
2329 ;; abs(x)^(2*n+1) -> abs(x)*x^(2*n), n an integer number
2331 (return (mul (power (cadr gr
) (add pot -
1))
2333 (return (mul (power (cadr gr
) (add pot
1))
2335 ((eq (caar gr
) 'mequal
)
2336 (return (eqtest (list (ncons (caar gr
))
2337 (power (cadr gr
) pot
)
2338 (power (caddr gr
) pot
))
2340 ((symbolp pot
) (go opp
))
2341 ((eq (caar gr
) 'mexpt
) (go e1
))
2342 ((and (eq (caar gr
) '%sum
)
2347 (return (do ((i (1- pot
) (1- i
))
2348 (an gr
(simptimes (list '(mtimes) an gr
) 1 t
)))
2351 (return (eqtest (testt (tms gr pot nil
)) check
)))
2353 (return (eqtest (cond ((and (mplusp gr
)
2354 (not (or (> pot $expop
)
2355 (> (- pot
) $expon
))))
2356 (expandexpt gr pot
))
2357 (t (simplifya (tms gr pot nil
) t
)))
2361 (cond ((eq (caar gr
) 'mexpt
) (go e1
))
2362 ((eq (caar gr
) 'rat
)
2363 (return (mul2 (power (cadr gr
) pot
)
2364 (power (caddr gr
) (mul2 -
1 pot
)))))
2365 ((not (eq (caar gr
) 'mtimes
)) (go up
))
2366 ((or (eq $radexpand
'$all
) (and $radexpand
(simplexpon pot
)))
2369 ((and (or (not (numberp (cadr gr
)))
2370 (equal (cadr gr
) -
1))
2371 (equal -
1 ($num gr
)) ; only for -1
2372 ;; Do not simplify for a complex base.
2373 (not (member ($csign gr
) '($complex $imaginary
)))
2374 (and (eq $domain
'$real
) $radexpand
))
2375 ;; (-1/x)^a -> 1/(-x)^a for x negative
2376 ;; For all other cases (-1)^a/x^a
2377 (if (eq ($csign
(setq w
($denom gr
))) '$neg
)
2378 (return (inv (power (neg w
) pot
)))
2379 (return (div (power -
1 pot
)
2381 ((or (eq $domain
'$complex
) (not $radexpand
)) (go up
)))
2382 (return (do ((l (cdr gr
) (cdr l
)) (res (ncons 1)) (rad))
2384 (cond ((equal res
'(1))
2385 (eqtest (list '(mexpt) gr pot
) check
))
2387 (testt (cons '(mtimes simp
) res
)))
2389 (setq rad
(power* ; RADEXPAND=()?
2390 (cons '(mtimes) (nreverse rad
)) pot
))
2391 (cond ((not (onep1 rad
))
2393 (testt (tms rad
1 (cons '(mtimes) res
))))
2395 (setq rulesw nil res
(cdr rad
))))))
2396 (eqtest (testt (cons '(mtimes) res
)) check
))))
2397 ;; Check with $csign to be more complete. This prevents wrong
2398 ;; simplifications like sqrt(-z^2)->%i*sqrt(z^2) for z complex.
2399 (setq z
($csign
(car l
)))
2400 (if (member z
'($complex $imaginary
))
2401 (setq z
'$pnz
)) ; if appears complex, unknown sign
2402 (setq w
(cond ((member z
'($neg $nz
))
2403 (setq rad
(cons -
1 rad
))
2407 ((alike1 w gr
) (return (list '(mexpt simp
) gr pot
)))
2408 ((member z
'($pn $pnz
))
2409 (setq rad
(cons w rad
)))
2411 (setq w
(testt (tms (simplifya (list '(mexpt) w pot
) t
)
2412 1 (cons '(mtimes) res
))))))
2413 (cond (rulesw (setq rulesw nil res
(cdr w
))))))
2416 (cond ((and (cdr res
) (onep1 (car res
)) (ratnump (cadr res
)))
2417 (setq res
(cdr res
))))
2418 (cond ((null (setq gr
(cdr gr
)))
2419 (return (eqtest (testt (cons '(mtimes) res
)) check
)))
2421 (setq y
(power (cadar gr
) (mult (caddar gr
) pot
))))
2423 (setq y
(%itopot pot
)))
2425 (setq y
(list '(mexpt) (car gr
) pot
)))
2426 (t (setq y
(list '(mexpt simp
) (car gr
) pot
))))
2427 (setq w
(testt (tms (simplifya y t
) 1 (cons '(mtimes) res
))))
2428 (cond (rulesw (setq rulesw nil res
(cdr w
))))
2432 (return (exptrl gr pot
))
2435 (cond ((zerop1 pot
) (go retno
))
2437 (let ((y (mget gr
'$numer
)))
2438 (if (and y
(floatp y
) (or $numer
(not (equal pot
1))))
2439 ;; A numeric constant like %e, %pi, ... and
2440 ;; exponent is a float or bigfloat value.
2441 (return (if (and (member gr
*builtin-numeric-constants
*)
2442 (equal pot bigfloatone
))
2443 ;; Return a bigfloat value.
2445 ;; Return a float value.
2447 ;; In all other cases exptrl simplifies accordingly.
2448 (return (exptrl gr pot
)))))
2450 ;; Numerically evaluate if the power is a flonum.
2452 (let ((val (flonum-eval '%exp pot
)))
2453 (if (float-inf-p val
)
2454 ;; needed for gcl - no trap of overflow
2455 (signal 'floating-point-overflow
))
2458 ;; Numerically evaluate if the power is a (complex)
2459 ;; big-float. (This is basically the guts of
2460 ;; big-float-eval, but we can't use big-float-eval.)
2461 (when (and (not (member 'simp
(car x
)))
2462 (complex-number-p pot
'bigfloat-or-number-p
))
2463 (let ((x ($realpart pot
))
2464 (y ($imagpart pot
)))
2465 (cond ((and ($bfloatp x
) (like 0 y
))
2466 (return ($bfloat
`((mexpt simp
) $%e
,pot
))))
2467 ((or ($bfloatp x
) ($bfloatp y
))
2468 (let ((z (add ($bfloat x
) (mul '$%i
($bfloat y
)))))
2469 (setq z
($rectform
`((mexpt simp
) $%e
,z
)))
2470 (return ($bfloat z
))))))))
2471 (cond ((and $logsimp
(among '%log pot
)) (return (%etolog pot
)))
2472 ((and $demoivre
(setq z
(demoivre pot
))) (return z
))
2476 ;; Exponent contains %i and %pi and %emode is TRUE:
2477 ;; Check simplification of exp(%i*%pi*p/q*x)
2478 (setq z
(%especial pot
)))
2480 (($taylorp
(third x
))
2481 ;; taylorize %e^taylor(...)
2482 (return ($taylor x
)))))
2484 (let ((y (mget gr
'$numer
)))
2485 ;; Check for a numeric constant.
2489 ;; The exponent is a bigfloat. Convert base to bigfloat.
2491 (member gr
*builtin-numeric-constants
*)
2492 (setq y
($bfloat gr
)))
2493 (and $numer
(integerp pot
)))
2494 (return (exptrl y pot
))))))
2497 (return (eqtest (list '(mexpt) gr pot
) check
))
2501 (cond ((mxorlistp1 gr
) (return (constmx (addk 1 pot
) gr
)))
2503 ((onep1 pot
) (return gr
))
2504 ((or $doallmxops $doscmxops $domxexpt
)
2505 (cond ((or (and mlpgr
2506 (or (not ($listp gr
)) $listarith
)
2507 (scalar-or-constant-p pot $assumescalar
))
2510 (or (not ($listp pot
)) $listarith
)
2511 (scalar-or-constant-p gr $assumescalar
)))
2512 (return (simplifya (outermap1 'mexpt gr pot
) t
)))
2514 ((and $domxmxops
(member pot
'(-1 -
1.0)))
2515 (return (simplifya (outermap1 'mexpt gr pot
) t
)))
2518 ;; At this point we have an expression: (z^a)^b with gr = z^a and pot = b
2519 (cond ((or (eq $radexpand
'$all
)
2520 ;; b is an integer or an odd rational
2522 (and (eq $domain
'$complex
)
2523 (not (member ($csign
(caddr gr
)) '($complex $imaginary
)))
2524 ;; z >= 0 and a not a complex
2525 (or (member ($csign
(cadr gr
)) '($pos $pz $zero
))
2527 (and (mnump (caddr gr
))
2528 (eq ($sign
(sub 1 (take '(mabs) (caddr gr
))))
2530 (and (eq $domain
'$real
)
2531 (member ($csign
(cadr gr
)) '($pos $pz $zero
)))
2532 ;; (1/z)^a -> 1/z^a when z a constant complex
2533 (and (eql (caddr gr
) -
1)
2535 (eq $domain
'$real
))
2536 (and (eq ($csign
(cadr gr
)) '$complex
)
2537 ($constantp
(cadr gr
)))))
2538 ;; This does (1/z)^a -> 1/z^a. This is in general wrong.
2539 ;; We switch this type of simplification on, when
2540 ;; $ratsimpexpons is T. E.g. radcan sets this flag to T.
2541 ;; radcan hangs for expressions like sqrt(1/(1+x)) without
2542 ;; this simplification.
2544 (equal (caddr gr
) -
1))
2547 (odnump (caddr gr
))))
2548 ;; Simplify (z^a)^b -> z^(a*b)
2549 (setq pot
(mul pot
(caddr gr
))
2551 ((and (eq $domain
'$real
)
2554 (not (apparently-complex-to-judge-by-$csign-p
(cadr gr
)))
2555 (evnump (caddr gr
)))
2556 ;; Simplify (x^a)^b -> abs(x)^(a*b)
2557 (setq pot
(mul pot
(caddr gr
))
2558 gr
(radmabs (cadr gr
))))
2561 (mminusp (caddr gr
)))
2562 ;; Simplify (1/z^a)^b -> 1/(z^a)^b
2564 gr
(power (cadr gr
) (neg (caddr gr
)))))
2568 (defun apparently-complex-to-judge-by-$csign-p
(e)
2569 (let ((s ($csign e
)))
2570 (member s
'($complex $imaginary
))))
2572 (defun apparently-real-to-judge-by-$csign-p
(e)
2573 (let ((s ($csign e
)))
2574 (member s
'($pos $neg $zero $pn $pnz $pz $nz
))))
2576 ;; Basically computes log of m base b. Except if m is not a power
2577 ;; of b, we return nil. m is a positive integer and base an integer
2578 ;; not equal to +/-1.
2579 (defun exponent-of (m base
)
2580 ;; Just compute base^k until base^k >= m. Then check if they're equal.
2581 ;; If so, we have the exponent. Otherwise, give up.
2584 (multiple-value-bind (q r
)
2590 (if (zerop expo
) nil expo
)))
2592 (defun timesin (x y w
) ; Multiply X^W into Y
2593 (prog (fm temp z check u expo
)
2594 (if (mexptp x
) (setq check x
))
2596 ;; Prepare the factor x^w and initialize the work of timesin
2600 (setq temp
(cons '(mexpt) (if check
2601 (list (cadr x
) (mult (caddr x
) w
))
2603 (if (and (not timesinp
) (not (eq x
'$%i
)))
2605 (setq temp
(simplifya temp t
))))))
2606 (setq x
(if (mexptp temp
)
2612 ;; Go through the list of terms in fm and look what is to do.
2613 (cond ((null (cdr fm
))
2614 ;; The list of terms is empty. The loop is finshed.
2616 ((or (and (mnump temp
)
2617 (not (or (integerp temp
)
2619 (and (integerp temp
)
2621 ;; Stop the loop for a float or bigfloat number, or number -1.
2624 (cond ((alike1 (car x
) (cadadr fm
))
2625 (cond ((zerop1 (setq w
(plsk (caddr (cadr fm
)) w
)))
2630 (rplacd fm
(cddr fm
))
2631 (cond ((mnump (setq x
(if (mnump (car x
))
2633 (power (car x
) w
))))
2634 (return (rplaca y
(timesk (car y
) x
))))
2639 x
(if (mexptp x
) (cdr x
) (list x
1)))
2643 ((maxima-constantp (car x
))
2646 (cond ((mtimesp (car x
))
2647 ;; A base which is a mtimes expression. Remove
2648 ;; the factor from the lists of products.
2649 (rplacd fm
(cddr fm
))
2650 ;; Multiply the factors of the base with
2651 ;; the list of all remaining products.
2653 (return (muln (nconc y
(cdar x
)) t
)))
2654 (t (return (rplaca (cdr fm
) (car x
))))))
2657 ;; At this place we have to add code for a rational number
2658 ;; as a factor to the list of products.
2660 (or (ratnump (car x
))
2661 (and (integerp (car x
))
2662 (not (onep (car x
))))))
2663 ;; Multiplying bas^k * num/den
2664 (let ((num (num1 (car x
)))
2665 (den (denom1 (car x
)))
2666 (bas (second (cadr fm
))))
2667 (cond ((and (integerp bas
)
2668 (not (eql 1 (abs bas
)))
2669 (setq expo
(exponent-of (abs num
) bas
)))
2670 ;; We have bas^m*bas^k = bas^(k+m).
2671 (setq temp
(power bas
2672 (add (third (cadr fm
)) expo
)))
2673 ;; Set fm to have 1/denom term.
2674 (setq x
(mul (car y
)
2678 ((and (integerp bas
)
2679 (not (eql 1 (abs bas
)))
2680 (setq expo
(exponent-of den bas
)))
2681 (setq expo
(- expo
))
2682 ;; We have bas^(-m)*bas^k = bas^(k-m).
2683 (setq temp
(power bas
2684 (add (third (cadr fm
)) expo
)))
2685 ;; Set fm to have the numerator term.
2686 (setq x
(mul (car y
)
2689 (exptrl bas
(- expo
)))))))
2691 ;; Next term in list of products.
2694 ;; Add in the bas^(k+m) term or bas^(k-m)
2695 (setf y
(rplaca y
1))
2696 (rplacd fm
(cddr fm
))
2697 (rplacd fm
(cons temp
(cdr fm
)))
2703 ((and (not (atom (car x
)))
2704 (eq (caar (car x
)) 'mabs
)
2706 (integerp (caddr (cadr fm
)))
2707 (< (caddr (cadr fm
)) -
1)
2708 (alike1 (cadr (car x
)) (cadr (cadr fm
)))
2709 (not (member ($csign
(cadr (car x
)))
2710 '($complex imaginary
))))
2711 ;; 1/x^n*abs(x) -> 1/(x^(n-2)*abs(x)), where n an integer
2712 ;; Replace 1/x^n -> 1/x^(n-2)
2713 (setq temp
(power (cadr (cadr fm
))
2714 (add (caddr (cadr fm
)) 2)))
2715 (rplacd fm
(cddr fm
))
2716 (if (not (equal temp
1))
2717 (rplacd fm
(cons temp
(cdr fm
))))
2718 ;; Multiply factor 1/abs(x) into list of products.
2719 (setq x
(list (car x
) -
1))
2720 (setq temp
(power (car x
) (cadr x
)))
2724 ((and (not (atom (car x
)))
2725 (eq (caar (car x
)) 'mabs
)
2727 (integerp (caddr (cadr fm
)))
2728 (> (caddr (cadr fm
)) 1)
2729 (alike1 (cadr (car x
)) (cadr (cadr fm
)))
2730 (not (member ($csign
(cadr (car x
)))
2731 '($complex imaginary
))))
2732 ;; x^n/abs(x) -> x^(n-2)*abs(x), where n an integer.
2733 ;; Replace x^n -> x^(n-2)
2734 (setq temp
(power (cadr (cadr fm
))
2735 (add (caddr (cadr fm
)) -
2)))
2736 (rplacd fm
(cddr fm
))
2737 (if (not (equal temp
1))
2738 (rplacd fm
(cons temp
(cdr fm
))))
2739 ;; Multiply factor abs(x) into list of products.
2740 (setq x
(list (car x
) 1))
2741 (setq temp
(power (car x
) (cadr x
)))
2745 ((and (not (atom (cadr fm
)))
2746 (not (atom (cadr (cadr fm
))))
2747 (eq (caaadr (cadr fm
)) 'mabs
)
2748 (equal (caddr (cadr fm
)) -
1)
2751 (alike1 (cadadr (cadr fm
)) (car x
))
2752 (not (member ($csign
(cadadr (cadr fm
)))
2753 '($complex imaginary
))))
2754 ;; 1/abs(x)*x^n -> x^(n-2)*abs(x), where n an integer.
2755 ;; Replace 1/abs(x) -> abs(x)
2756 (setq temp
(cadr (cadr fm
)))
2757 (rplacd fm
(cddr fm
))
2758 (rplacd fm
(cons temp
(cdr fm
)))
2759 ;; Multiply factor x^(n-2) into list of products.
2760 (setq x
(list (car x
) (add (cadr x
) -
2)))
2761 (setq temp
(power (car x
) (cadr x
)))
2765 ((or (maxima-constantp (car x
))
2766 (maxima-constantp (cadadr fm
)))
2767 (if (great temp
(cadr fm
))
2769 ((great (car x
) (cadadr fm
))
2772 ((alike1 (car x
) (cadr fm
))
2775 ;; When a number goto start and look in the next term.
2779 ((and (not (atom (cadr fm
)))
2780 (eq (caar (cadr fm
)) 'mabs
)
2783 (alike1 (cadr (cadr fm
)) (car x
))
2784 (not (member ($csign
(cadr (cadr fm
)))
2785 '($complex imaginary
))))
2786 ;; abs(x)/x^n -> 1/(x^(n-2)*abs(x)), where n an integer.
2787 ;; Replace abs(x) -> 1/abs(x).
2788 (setq temp
(power (cadr fm
) -
1))
2789 (rplacd fm
(cddr fm
))
2790 (rplacd fm
(cons temp
(cdr fm
)))
2791 ;; Multiply factor x^(-n+2) into list of products.
2792 (setq x
(list (car x
) (add (cadr x
) 2)))
2793 (setq temp
(power (car x
) (cadr x
)))
2797 ((maxima-constantp (car x
))
2798 (when (great temp
(cadr fm
))
2800 ((great (car x
) (cadr fm
))
2804 ;; Multiply a number into the list of products.
2805 (return (rplaca y
(timesk (car y
) temp
))))
2806 ((and (eq (car x
) '$%i
)
2809 ((and (eq (car x
) '$%e
)
2812 (return (rplaca y
(timesk (car y
) (exp (float w
))))))
2814 (not (constant (car x
))))
2816 ;; At this point we will insert a mexpt expression,
2817 ;; but first we look at the car of the list of products and
2818 ;; modify the expression if we found a rational number.
2820 (not (onep1 (car y
)))
2821 (or (integerp (car y
))
2823 ;; Multiplying bas^k * num/den.
2824 (let ((num (num1 (car y
)))
2825 (den (denom1 (car y
)))
2827 (cond ((and (integerp bas
)
2828 (not (eql 1 (abs bas
)))
2829 (setq expo
(exponent-of (abs num
) bas
)))
2830 ;; We have bas^m*bas^k.
2831 (setq temp
(power bas
(add (cadr x
) expo
)))
2832 ;; Set fm to have 1/denom term.
2833 (setq x
(div (div num
(exptrl bas expo
)) den
)))
2834 ((and (integerp bas
)
2835 (not (eql 1 (abs bas
)))
2836 (setq expo
(exponent-of den bas
)))
2837 (setq expo
(- expo
))
2838 ;; We have bas^(-m)*bas^k.
2839 (setq temp
(power bas
(add (cadr x
) expo
)))
2840 ;; Set fm to have the numerator term.
2841 (setq x
(div num
(div den
(exptrl bas
(- expo
))))))
2843 ;; The rational doesn't contain any (simple) powers of
2844 ;; the exponential term. We're done.
2845 (return (cdr (rplacd fm
(cons temp
(cdr fm
)))))))
2846 ;; Add in the a^(m+k) or a^(k-m) term.
2847 (setf y
(rplaca y
1))
2848 (rplacd fm
(cons temp
(cdr fm
)))
2854 ((and (maxima-constantp (car x
))
2855 (do ((l (cdr fm
) (cdr l
)))
2857 (when (and (mexptp (cadr l
))
2858 (alike1 (car x
) (cadadr l
)))
2862 ((or (and (mnump (car x
))
2864 (and (eq (car x
) '$%e
)
2868 (setq u
(%especial w
))))
2870 ((alike (cdr check
) x
)
2873 (exptrl (car x
) w
))))
2875 (return (rplaca y
(timesk (car y
) x
))))
2879 (return (cdr (rplacd fm
(cons x
(cdr fm
))))))
2889 (setq temp
(list '(mexpt) (car x
) w
))
2890 (setq temp
(eqtest temp
(or check
'((foo)))))
2891 (return (cdr (rplacd fm
(cons temp
(cdr fm
)))))))
2893 (return (cdr (rplacd fm
(cons (car x
) (cdr fm
)))))
2898 (cond ((and (eq (car x
) '$%i
) (equal w
1))
2899 (rplacd fm
(cddr fm
))
2900 (return (rplaca y
(timesk -
1 (car y
)))))
2901 ((zerop1 (setq w
(plsk 1 w
)))
2903 ((and (mnump (car x
)) (mnump w
))
2904 (return (rplaca (cdr fm
) (exptrl (car x
) w
))))
2905 ((maxima-constantp (car x
))
2908 (setq z
(list '(mexpt) (car x
) w
))
2909 (cond ((alike1 (setq x
(simplifya z t
)) z
)
2910 (return (rplaca (cdr fm
) x
)))
2912 (rplacd fm
(cddr fm
))
2914 (return (muln (cons x y
) t
))))
2916 (rplacd fm
(cddr fm
))
2917 (setq x
(car x
) check nil
)
2920 (setq z
(tms x
1 (setq temp
(cons '(mtimes) y
))))
2921 (return (cond ((eq z temp
)
2924 (setq rulesw t
) z
)))
2926 (return (rplacd fm
(cddr fm
)))
2928 (if (minusp (setq w
(rem w
4)))
2930 (return (cond ((zerop w
)
2933 (rplaca y
(timesk -
1 (car y
))))
2935 (rplaca y
(timesk -
1 (car y
)))
2936 (rplacd fm
(cons '$%i
(cdr fm
))))
2938 (rplacd fm
(cons '$%i
(cdr fm
))))))))
2940 (defun simpmatrix (x vestigial z
)
2941 (declare (ignore vestigial
))
2942 (if (and (null (cddr x
))
2944 (or (eq $scalarmatrixp
'$all
) (member 'mult
(cdar x
)))
2945 ($listp
(cadr x
)) (cdadr x
) (null (cddadr x
)))
2946 (simplifya (cadadr x
) z
)
2947 (let ((badp (dolist (row (cdr x
)) (if (not ($listp row
)) (return t
))))
2948 (args (simpmap (cdr x
) z
)))
2949 (if (and args
(not badp
)) (matcheck args
))
2950 (cons (if badp
'(%matrix simp
) '($matrix simp
)) args
))))
2952 (defun %itopot
(pot)
2954 (let ((i (boole boole-and pot
3)))
2958 (t (list '(mtimes simp
) -
1 '$%i
))))
2959 (power -
1 (mul2 pot
'((rat simp
) 1 2)))))
2962 (cond ((eq (caar pot
) '%log
) (simplifya (cadr pot
) nil
))
2963 ((and (eq (caar pot
) 'mtimes
)
2964 (or (maxima-integerp (cadr pot
))
2965 (and $%e_to_numlog
($numberp
(cadr pot
))))
2966 (not (atom (caddr pot
))) (eq (caar (caddr pot
)) '%log
)
2968 (power (cadr (caddr pot
)) (cadr pot
)))))
2972 loop
(cond ((null pot
)
2973 (cond (a (setq a
(cons '(mtimes) a
))))
2974 (cond (c (setq c
(list '(mexpt simp
) '$%e
(addn c nil
)))))
2975 (return (cond ((null c
) (simptimes a
1 nil
))
2977 (t (simptimes (append a
(list c
)) 1 nil
)))))
2978 ((and (among '%log
(car pot
)) (setq b
(mnlogp (car pot
))))
2979 (setq a
(cons b a
)))
2980 (t (setq c
(cons (car pot
) c
))))
2981 (setq pot
(cdr pot
))
2984 (defun %etolog
(pot) (cond ((mnlogp pot
))
2985 ((eq (caar pot
) 'mplus
) (mnlog (cdr pot
)))
2986 (t (list '(mexpt simp
) '$%e pot
))))
2988 (defun zerores (r1 r2
)
2989 (cond ((or ($bfloatp r1
) ($bfloatp r2
)) bigfloatzero
)
2990 ((or (floatp r1
) (floatp r2
)) 0.0)
2993 (defmfun $orderlessp
(a b
)
2994 (setq a
($totaldisrep
(specrepcheck a
))
2995 b
($totaldisrep
(specrepcheck b
)))
2996 (and (not (alike1 a b
)) (great b a
) t
))
2998 (defmfun $ordergreatp
(a b
)
2999 (setq a
($totaldisrep
(specrepcheck a
))
3000 b
($totaldisrep
(specrepcheck b
)))
3001 (and (not (alike1 a b
)) (great a b
) t
))
3003 ;; Test function to order a and b by magnitude. If it is not possible to
3004 ;; order a and b by magnitude they are ordered by great. This function
3005 ;; can be used by sort, e.g. sort([3,1,7,x,sin(1),minf],ordermagnitudep)
3006 (defmfun $ordermagnitudep
(a b
)
3008 (setq a
($totaldisrep
(specrepcheck a
))
3009 b
($totaldisrep
(specrepcheck b
)))
3010 (cond ((and (or (constp a
) (member a
'($inf $minf
)))
3011 (or (constp b
) (member b
'($inf $minf
)))
3012 (member (setq sgn
($csign
(sub b a
))) '($pos $neg $zero
)))
3013 (cond ((eq sgn
'$pos
) t
)
3014 ((eq sgn
'$zero
) (and (not (alike1 a b
)) (great b a
)))
3016 ((or (constp a
) (member a
'($inf $minf
))) t
)
3017 ((or (constp b
) (member b
'($inf $minf
))) nil
)
3018 (t (and (not (alike1 a b
)) (great b a
))))))
3020 (defun evnump (n) (or (even n
) (and (ratnump n
) (even (cadr n
)))))
3021 (defun odnump (n) (or (and (integerp n
) (oddp n
))
3022 (and (ratnump n
) (oddp (cadr n
)))))
3024 (defun simplexpon (e)
3025 (or (maxima-integerp e
)
3026 (and (eq $domain
'$real
) (ratnump e
) (oddp (caddr e
)))))
3028 ;; This function is not called in Maxima core or share code
3029 ;; and can be cut out.
3031 (and (free p
'$%i
) (member ($sign p
) '($pos $pz $zero
))))
3034 (if (and limitp
(free e
'$%i
)) (asksign-p-or-n e
))
3035 (simplifya (list '(mabs) e
) t
))
3037 (defun simpmqapply (exp y z
)
3038 (let ((simpfun (and (not (atom (cadr exp
))) (safe-get (caaadr exp
) 'specsimp
))) u
)
3040 (funcall simpfun exp y z
)
3041 (progn (setq u
(simpargs exp z
))
3042 (if (symbolp (cadr u
))
3043 (simplifya (cons (cons (cadr u
) (cdar u
)) (cddr u
)) z
)
3046 ;; TRUE, if the symbol e is declared to be $complex or $imaginary.
3047 (defun decl-complexp (e)
3049 (kindp e
'$complex
)))
3051 ;; TRUE, if the symbol e is declared to be $real, $rational, $irrational
3053 (defun decl-realp (e)
3055 (or (kindp e
'$real
)
3056 (kindp e
'$rational
)
3057 (kindp e
'$irrational
)
3058 (kindp e
'$integer
))))
3060 ;; WARNING: Exercise extreme caution when modifying this function!
3062 ;; Richard Fateman and Stavros Macrakis both say that changing the
3063 ;; actual ordering relations (as opposed to making them faster to
3064 ;; determine) could have very subtle and wide-ranging effects. Also,
3065 ;; the simplifier spends the vast majority of its time here, so be
3066 ;; very careful about changes that may drastically slow down the
3074 (cond ((zerop y
) (floatp x
)) (t (plusp y
))))))
3076 (cond ((constant y
) (alphalessp y x
)) (t (numberp y
))))
3078 (cond ((mget y
'$scalar
) (alphalessp y x
))
3079 (t (maxima-constantp y
))))
3081 (cond ((mget y
'$mainvar
) (alphalessp y x
)) (t t
)))
3082 (t (or (maxima-constantp y
) (mget y
'$scalar
)
3083 (and (not (mget y
'$mainvar
)) (not (null (alphalessp y x
))))))))
3084 (t (not (ordfna y x
)))))
3085 ((atom y
) (ordfna x y
))
3087 (cond ((eq (caar y
) 'rat
)
3088 (> (* (caddr y
) (cadr x
)) (* (caddr x
) (cadr y
))))))
3089 ((eq (caar y
) 'rat
))
3090 ((or (member (caar x
) '(mtimes mplus mexpt %del
))
3091 (member (caar y
) '(mtimes mplus mexpt %del
)))
3093 ((and (eq (caar x
) 'bigfloat
) (eq (caar y
) 'bigfloat
)) (mgrp x y
))
3094 ((or (eq (caar x
) 'mrat
) (eq (caar y
) 'mrat
))
3095 (error "GREAT: internal error: unexpected MRAT argument"))
3096 (t (do ((x1 (margs x
) (cdr x1
)) (y1 (margs y
) (cdr y1
))) (())
3098 (return (cond (y1 nil
)
3099 ((not (alike1 (mop x
) (mop y
)))
3100 (great (mop x
) (mop y
)))
3101 ((member 'array
(cdar x
)) t
))))
3102 ((null y1
) (return t
))
3103 ((not (alike1 (car x1
) (car y1
)))
3104 (return (great (car x1
) (car y1
)))))))))
3106 ;; Trivial function used only in ALIKE1.
3107 ;; Should be defined as an open-codable subr.
3109 (defmacro memqarr
(l)
3110 `(if (member 'array
,l
) t
))
3112 ;; Compares two Macsyma expressions ignoring SIMP flags and all other
3113 ;; items in the header except for the ARRAY flag.
3120 (and (arrayp y
) (lisp-array-alike1 x y
)))
3122 ;; NOT SURE IF WE WANT TO ENABLE COMPARISON OF MAXIMA ARRAYS
3123 ;; ASIDE FROM THAT, ADD2LNC CALLS ALIKE1 (VIA MEMALIKE) AND THAT CAUSES TROUBLE
3124 ;; ((maxima-declared-arrayp x)
3125 ;; (and (maxima-declared-arrayp y) (maxima-declared-array-alike1 x y)))
3126 ;; ((maxima-undeclared-arrayp x)
3127 ;; (and (maxima-undeclared-arrayp y) (maxima-undeclared-array-alike1 x y)))
3132 (not (atom (car x
)))
3133 (not (atom (car y
)))
3134 (eq (caar x
) (caar y
)))
3137 ;; Punt back to LIKE, which handles specreps
3140 (eq (memqarr (cdar x
)) (memqarr (cdar y
)))
3141 (alike (cdr x
) (cdr y
))))))))
3143 (defun lisp-array-alike1 (x y
)
3145 (equal (array-dimensions x
) (array-dimensions y
))
3147 (dotimes (i (array-total-size x
))
3148 (if (not (alike1 (row-major-aref x i
) (row-major-aref y i
)))
3149 (return-from lisp-array-alike1 nil
)))
3152 (defun maxima-declared-array-alike1 (x y
)
3153 (lisp-array-alike1 (get (mget x
'array
) 'array
) (get (mget y
'array
) 'array
)))
3155 (defun maxima-undeclared-array-alike1 (x y
)
3157 (alike1 (mfuncall '$arrayinfo x
) (mfuncall '$arrayinfo y
))
3158 (alike1 ($listarray x
) ($listarray y
))))
3160 ;; Maps ALIKE1 down two lists.
3163 (do ((x x
(cdr x
)) (y y
(cdr y
))) ((atom x
) (equal x y
))
3164 (cond ((or (atom y
) (not (alike1 (car x
) (car y
))))
3167 (defun ordfna (e a
) ; A is an atom
3169 (or (not (eq (caar e
) 'rat
))
3170 (> (cadr e
) (* (caddr e
) a
))))
3172 (not (member (caar e
) '(mplus mtimes mexpt
))))
3173 (not (member (caar e
) '(rat bigfloat
))))
3174 ((eq (caar e
) 'mrat
)) ;; all MRATs succeed all atoms
3175 ((null (margs e
)) nil
)
3176 ((eq (caar e
) 'mexpt
)
3177 (cond ((and (maxima-constantp (cadr e
))
3178 (or (not (constant a
)) (not (maxima-constantp (caddr e
)))))
3179 (or (not (free (caddr e
) a
)) (great (caddr e
) a
)))
3180 ((eq (cadr e
) a
) (great (caddr e
) 1))
3181 (t (great (cadr e
) a
))))
3182 ((member (caar e
) '(mplus mtimes
))
3183 (let ((u (car (last e
))))
3184 (cond ((eq u a
) (not (ordhack e
))) (t (great u a
)))))
3185 ((eq (caar e
) '%del
))
3186 ((prog2 (setq e
(car (margs e
))) ; use first arg of e
3187 (and (not (atom e
)) (member (caar e
) '(mplus mtimes
))))
3188 (let ((u (car (last e
)))) ; and compare using
3189 (cond ((eq u a
) (not (ordhack e
))) ; same procedure as above
3194 ;; compare lists a and b elementwise from back to front
3195 (defun ordlist (a b cx cy
)
3197 (setq l1
(length a
) l2
(length b
))
3198 loop
(cond ((= l1
0)
3199 (return (cond ((= l2
0) (eq cx
'mplus
))
3200 ((and (eq cx cy
) (= l2
1))
3201 (great (cond ((eq cx
'mplus
) 0) (t 1)) (car b
))))))
3202 ((= l2
0) (return (not (ordlist b a cy cx
)))))
3203 (setq c
(nthelem l1 a
) d
(nthelem l2 b
))
3204 (cond ((not (alike1 c d
)) (return (great c d
))))
3205 (setq l1
(1- l1
) l2
(1- l2
))
3208 (defun term-list (x)
3213 (defun factor-list (x)
3218 ;; one of the exprs x or y should be one of:
3219 ;; %del, mexpt, mplus, mtimes
3221 (let ((cx (caar x
)) (cy (caar y
)))
3222 (cond ((eq cx
'%del
) (if (eq cy
'%del
) (great (cadr x
) (cadr y
)) t
))
3224 ((or (eq cx
'mtimes
) (eq cy
'mtimes
))
3225 (ordlist (factor-list x
) (factor-list y
) 'mtimes
'mtimes
))
3226 ((or (eq cx
'mplus
) (eq cy
'mplus
))
3227 (ordlist (term-list x
) (term-list y
) 'mplus
'mplus
))
3228 ((eq cx
'mexpt
) (ordmexpt x y
))
3229 ((eq cy
'mexpt
) (not (ordmexpt y x
))))))
3232 (if (and (cddr x
) (null (cdddr x
)))
3233 (great (if (eq (caar x
) 'mplus
) 0 1) (cadr x
))))
3235 (defun ordmexpt (x y
)
3236 (cond ((eq (caar y
) 'mexpt
)
3237 (cond ((alike1 (cadr x
) (cadr y
)) (great (caddr x
) (caddr y
)))
3238 ((maxima-constantp (cadr x
))
3239 (if (maxima-constantp (cadr y
))
3240 (if (or (alike1 (caddr x
) (caddr y
))
3241 (and (mnump (caddr x
)) (mnump (caddr y
))))
3242 (great (cadr x
) (cadr y
))
3243 (great (caddr x
) (caddr y
)))
3244 (great x
(cadr y
))))
3245 ((maxima-constantp (cadr y
)) (great (cadr x
) y
))
3247 (great (cadr x
) (if (mnump (caddr y
)) (cadr y
) y
)))
3248 ((mnump (caddr y
)) (great x
(cadr y
)))
3249 (t (let ((x1 (simpln1 x
)) (y1 (simpln1 y
)))
3250 (if (alike1 x1 y1
) (great (cadr x
) (cadr y
))
3252 ((alike1 (cadr x
) y
) (great (caddr x
) 1))
3253 ((mnump (caddr x
)) (great (cadr x
) y
))
3254 (t (great (simpln1 x
) (simpln (list '(%log
) y
) 1 t
)))))
3256 (defmfun $multthru
(e1 &optional e2
)
3258 (cond (e2 ;called with two args
3259 (setq arg1
(specrepcheck e1
)
3260 arg2
(specrepcheck e2
))
3261 (cond ((or (atom arg2
)
3262 (not (member (caar arg2
) '(mplus mequal
))))
3264 ((eq (caar arg2
) 'mequal
)
3265 (list (car arg2
) ($multthru arg1
(cadr arg2
))
3266 ($multthru arg1
(caddr arg2
))))
3267 (t (expandterms arg1
(cdr arg2
)))))
3268 (t ;called with only one arg
3270 (setq arg1
(setq arg2
(specrepcheck e1
)))
3271 (cond ((atom arg1
) (return arg1
))
3272 ((eq (caar arg1
) 'mnctimes
)
3273 (setq arg1
(cdr arg1
)) (go nct
))
3274 ((not (eq (caar arg1
) 'mtimes
)) (return arg1
)))
3275 (setq arg1
(reverse (cdr arg1
)))
3276 times
(when (mplusp (car arg1
))
3277 (setq l1
(nconc l1
(cdr arg1
)))
3278 (return (expandterms (muln l1 t
) (cdar arg1
))))
3279 (setq l1
(cons (car arg1
) l1
))
3280 (setq arg1
(cdr arg1
))
3281 (if (null arg1
) (return arg2
))
3283 nct
(when (mplusp (car arg1
))
3284 (setq l1
(nreverse l1
))
3285 (return (addn (mapcar
3289 (append l1
(ncons u
) (cdr arg1
)))
3293 (setq l1
(cons (car arg1
) l1
))
3294 (setq arg1
(cdr arg1
))
3295 (if (null arg1
) (return arg2
))
3298 ;; EXPANDEXPT computes the expansion of (x1 + x2 + ... + xm)^n
3299 ;; taking a sum and integer power as arguments.
3300 ;; Its theory is to recurse down the binomial expansion of
3301 ;; (x1 + (x2 + x3 + ... + xm))^n using the Binomial Expansion
3302 ;; Thus it does a sigma:
3306 ;; \ / n \ k (n - k)
3307 ;; > | | x1 (x2 + x3 + ... + xm)
3312 ;; The function EXPONENTIATE-SUM does this and recurses through the second
3313 ;; sum raised to a power. It takes a list of terms and a positive integer
3314 ;; power as arguments.
3317 (defun expandexpt (sum power
)
3318 (declare (fixnum power
))
3319 (let ((expansion (exponentiate-sum (cdr sum
) (abs power
))))
3320 (cond ((plusp power
) expansion
)
3321 (t (inv expansion
)))))
3323 (defun exponentiate-sum (terms rpower
)
3324 (declare (fixnum rpower
))
3325 (cond ((= rpower
0) 1)
3326 ((null (cdr terms
)) (power (car terms
) rpower
))
3327 ((= rpower
1) (cons '(mplus simp
) terms
))
3328 (t (do ((i 0 (1+ i
))
3329 (result 0 (add2 result
3330 (muln (list (combination rpower i
)
3331 (exponentiate-sum (cdr terms
)
3333 (power (car terms
) i
)) t
))))
3334 ((> i rpower
) result
)
3335 (declare (fixnum i
))))))
3337 ;; Computes the combination of n elements taken m at a time by the formula
3339 ;; (n * (n-1) * ... * (n - m + 1)) / m! =
3340 ;; (n / 1) * ((n - 1) / 2) * ... * ((n - m + 1) / m)
3342 ;; Checks for the case when m is greater than n/2 and translates
3343 ;; to an equivalent expression.
3345 (defun combination (n m
)
3346 (declare (fixnum n m
))
3347 (cond ((> m
(truncate n
2))
3348 (combination n
(- n m
)))
3350 (do ((result 1 (truncate (* result n1
) m1
))
3354 (declare (fixnum n1 m1
))))))
3356 (defun expandsums (a b
)
3358 (setq a
(fixexpand a
) b
(cdr b
))
3360 (when (null a
) (return c
))
3361 (setq c
(cons (expandterms (car a
) b
) c
))
3366 (defun expandterms (a b
)
3369 (when (null b
) (return c
))
3370 (setq c
(cons (mul2 a
(car b
)) c
))
3375 (defun genexpands (l)
3380 (setq prods
(nreverse prods
)
3381 negprods
(nreverse negprods
)
3382 sums
(nreverse sums
)
3383 negsums
(nreverse negsums
))
3386 (push (car l
) prods
))
3387 ((eq (caaar l
) 'rat
)
3388 (unless (equal (cadar l
) 1)
3389 (push (cadar l
) prods
))
3390 (push (caddar l
) negprods
))
3391 ((eq (caaar l
) 'mplus
)
3392 (push (car l
) sums
))
3393 ((and (eq (caaar l
) 'mexpt
)
3394 (equal (caddar l
) -
1) (mplusp (cadar l
)))
3395 (push (cadar l
) negsums
))
3396 ((and (eq (caaar l
) 'mexpt
)
3398 (mminusp (caddar l
))))
3399 (push (if (equal (caddar l
) -
1)
3401 (list (caar l
) (cadar l
) (neg (caddar l
))))
3404 (push (car l
) prods
)))
3407 (defun expandtimes (a)
3408 (prog (prods negprods sums negsums expsums expnegsums
)
3410 (setq prods
(cond ((null prods
) 1)
3411 ((null (cdr prods
)) (car prods
))
3412 (t (cons '(mtimes simp
) prods
))))
3413 (setq negprods
(cond ((null negprods
) 1)
3414 ((null (cdr negprods
)) (car negprods
))
3415 (t (cons '(mtimes simp
) negprods
))))
3416 (cond ((null sums
) (go down
))
3417 (t (setq expsums
(car sums
))
3419 (setq expsums
(expandsums expsums c
)))
3421 (setq prods
(cond ((equal prods
1) expsums
)
3422 (t (expandterms prods
(fixexpand expsums
)))))
3423 down
(cond ((null negsums
)
3424 (cond ((equal 1 negprods
) (return prods
))
3426 (return (expandterms (power negprods -
1) (cdr prods
))))
3427 (t (return (let ((expandflag t
))
3428 (mul2 prods
(power negprods -
1)))))))
3430 (setq expnegsums
(car negsums
))
3432 (setq expnegsums
(expandsums expnegsums c
)))
3434 (setq expnegsums
(expandterms negprods
(fixexpand expnegsums
)))
3435 (return (if (mplusp prods
)
3436 (expandterms (inv expnegsums
) (cdr prods
))
3437 (let ((expandflag t
))
3438 (mul2 prods
(inv expnegsums
)))))))
3440 (defun expand1 (exp $expop $expon
)
3441 (unless (and (integerp $expop
) (> $expop -
1))
3442 (merror (intl:gettext
"expand: expop must be a nonnegative integer; found: ~M") $expop
))
3443 (unless (and (integerp $expon
) (> $expon -
1))
3444 (merror (intl:gettext
"expand: expon must be a nonnegative integer; found: ~M") $expon
))
3445 (resimplify (specrepcheck exp
)))
3447 (defmfun $expand
(exp &optional
(expop $maxposex
) (expon $maxnegex
))
3448 (expand1 exp expop expon
))
3450 (defun fixexpand (a)
3451 (if (not (mplusp a
))
3455 (defun simpnrt (x *n
) ; computes X^(1/*N)
3456 (prog (*in
*out varlist genvar $factorflag $dontfactor
)
3457 (setq $factorflag t
)
3459 (setq x
(ratrep* x
))
3460 (when (equal (cadr x
) 0) (return 0))
3461 (setq x
(ratfact (cdr x
) 'psqfr
))
3462 (simpnrt1 (mapcar #'pdis x
))
3463 (setq *out
(if *out
(muln *out nil
) 1))
3464 (setq *in
(cond (*in
3465 (setq *in
(muln *in nil
))
3468 (return (let (($%emode t
))
3469 (simplifya (list '(mtimes) *in
*out
)
3472 (member (caaadr *in
) '(mplus mtimes rat
)))))))))
3475 (do ((x x
(cddr x
)) (y))
3477 (cond ((not (equal 1 (setq y
(gcd (cadr x
) *n
))))
3478 (push (simpnrt (list '(mexpt) (car x
) (quotient (cadr x
) y
))
3481 ((and (equal (cadr x
) 1) (integerp (car x
)) (plusp (car x
))
3482 (setq y
(pnthrootp (car x
) *n
)))
3485 (unless (> *n
(abs (cadr x
)))
3486 (push (list '(mexpt) (car x
) (quotient (cadr x
) *n
)) *out
))
3487 (push (list '(mexpt) (car x
) (rem (cadr x
) *n
)) *in
)))))
3489 (defun nrthk (in *n
)
3495 ((eq $domain
'$real
)
3501 (power* '$%e
(list '(mtimes) (list '(rat) 1 *n
) '$%pi
'$%i
))))
3504 ((or (and wflag
(eq ($asksign in
) '$neg
))
3505 (and (mnump in
) (equal ($sign in
) '$neg
)))
3506 (nrthk1 (mul2* -
1 in
) *n
))
3510 (defun nrthk1 (in *n
) ; computes (-IN)^(1/*N)
3512 (mul2 (nrthk2 in
*n
) (nrthk -
1 *n
))
3513 (nrthk2 (mul2* -
1 in
) *n
)))
3515 (defun nrthk2 (in *n
)
3516 (power* in
(list '(rat) 1 *n
))) ; computes IN^(1/*N)
3518 ;; The following was formerly in SININT. This code was placed here because
3519 ;; SININT is now an out-of-core file on MC, and this code is needed in-core
3520 ;; because of the various calls to it. - BMT & JPG
3522 (declare-top (special var $ratfac ratform context
))
3524 (defmfun $integrate
(expr x
&optional lo hi
)
3525 (declare (special *in-risch-p
*))
3528 (with-new-context (context)
3529 (if (member '%risch
*nounl
*)
3531 ;; Give up; we're being called from RISCHINT by some path.
3532 (list '(%integrate
) expr x
)
3535 ($defint expr x lo hi
))))
3539 ((member (caar a
) '(mplus mtimes
))
3540 (do ((l (cdr a
) (cdr l
))) ((null l
) t
)
3541 (or (ratp (car l
) var
) (return nil
))))
3542 ((eq (caar a
) 'mexpt
)
3543 (if (free (cadr a
) var
)
3544 (free (caddr a
) var
)
3545 (and (integerp (caddr a
)) (ratp (cadr a
) var
))))
3548 (defun ratnumerator (r)
3550 ((atom (cdr r
)) (car r
))
3551 ((numberp (cadr r
)) r
)
3554 (defun ratdenominator (r)
3556 ((atom (cdr r
)) (cdr r
))
3557 ((numberp (cadr r
)) 1)
3560 (declare-top (special var
))
3562 ;; (BPROG U V) appears to return A and B (as ((A1 . A2) B1 . B2) with A = A1/A2, B = B1/B2)
3563 ;; such that B/U + A/V = 1/(U*V), where U, V are polynomials represented as a list of
3564 ;; exponents and coefficients, (<gensym> E1 C1 E2 C2 ...) = C1*<gensym>^E1 + C2*<gensym>^E2 + ....
3566 ;; (%i73) partfrac ((2*x^2-3)/(x^4-3*x^2+2), x);
3567 ;; 1. Trace: (PARTFRAC '((#:X16910 2 2 0 -3) #:X16910 4 1 2 -3 0 2) '#:X16910)
3568 ;; 2. Trace: (BPROG '(#:X16910 2 1 0 -2) '(#:X16910 2 1 0 -1))
3569 ;; 2. Trace: BPROG ==> ((-1 . 1) 1 . 1)
3570 ;; 2. Trace: (BPROG '(#:X16910 1 1 0 1) '(#:X16910 1 1 0 -1))
3571 ;; 2. Trace: BPROG ==> ((1 . 2) -1 . 2)
3572 ;; 2. Trace: (BPROG '(#:X16910 1 1 0 -1) '1)
3573 ;; 2. Trace: BPROG ==> ((0 . 1) 1 . 1)
3574 ;; 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))
3575 ;; (%o73) 1/(x^2-2)-1/(2*(x+1))+1/(2*(x-1))
3578 (prog (p1b p2b coef1r coef2r coef1s coef2s f1 f2 a egcd oldalg
)
3579 (setq oldalg $algebraic
)
3582 (setq coef2r
(setq coef1s
0))
3583 (setq coef2s
(setq coef1r
1))
3586 (unless (zerop (pdegree p1b var
)) (setq egcd
(pgcdexpon p1b
)))
3588 (unless (or (zerop (pdegree p2b var
)) (= egcd
1))
3589 (setq egcd
(gcd egcd
(pgcdexpon p2b
)))
3590 (setq p1b
(pexpon*// p1b egcd nil
)
3591 p2b
(pexpon*// p2b egcd nil
)))
3592 b1
(cond ((< (pdegree p1b var
) (pdegree p2b var
))
3594 (rotatef coef1r coef2r
)
3595 (rotatef coef1s coef2s
)))
3596 (when (zerop (pdegree p2b var
))
3597 (unless (zerop (pdegree coef2r var
))
3598 (setq coef2r
(pexpon*// coef2r egcd t
)))
3599 (unless (zerop (pdegree coef2s var
))
3600 (setq coef2s
(pexpon*// coef2s egcd t
)))
3601 (return (cons (ratreduce (ptimes (cdr r
) coef2r
) p2b
)
3602 (ratreduce (ptimes (cdr s
) coef2s
) p2b
))))
3603 (setq f1
(psquorem1 (cdr p1b
) (cdr p2b
) t
))
3604 (setq $algebraic $false
)
3605 (setq f2
(psimp var
(cadr f1
)))
3606 (setq p1b
(pquotientchk (psimp var
(caddr f1
)) a
))
3607 (setq $algebraic oldalg
)
3609 (setq coef1r
(pquotientchk (pdifference (ptimes f1 coef1r
)
3612 (setq coef1s
(pquotientchk (pdifference (ptimes f1 coef1s
)
3618 (defun ratdifference (a b
) (ratplus a
(ratminus b
)))
3620 (defun ratpl (a b
) (ratplus (ratfix a
) (ratfix b
)))
3622 (defun ratti (a b c
) (rattimes (ratfix a
) (ratfix b
) c
))
3624 (defun ratqu (a b
) (ratquotient (ratfix a
) (ratfix b
)))
3626 (defun ratfix (a) (cond ((equal a
(ratnumerator a
)) (cons a
1)) (t a
)))
3628 (defun ratdivide (f g
)
3629 (destructuring-let* (((fnum . fden
) (ratfix f
))
3630 ((gnum . gden
) (ratfix g
))
3631 ((q r
) (pdivide fnum gnum
)))
3632 (cons (ratqu (ratti q gden t
) fden
)
3635 (defun polcoef (l n
) (cond ((or (atom l
) (pointergp var
(car l
)))
3636 (cond ((equal n
0) l
) (t 0)))
3637 (t (ptterm (cdr l
) n
))))
3639 (defun disrep (l) (cond ((equal (ratnumerator l
) l
)
3640 ($ratdisrep
(cons ratform
(cons l
1))))
3641 (t ($ratdisrep
(cons ratform l
)))))
3643 (declare-top (unspecial var
))
3646 ;; The following was formerly in MATRUN. This code was placed here because
3647 ;; MATRUN is now an out-of-core file on MC, and this code is needed in-core
3648 ;; so that MACSYMA SAVE files will work. - JPG
3653 (defun kar (x) (if (atom x
) (matcherr) (car x
)))
3655 (defun kaar (x) (kar (kar x
)))
3657 (defun kdr (x) (if (atom x
) (matcherr) (cdr x
)))
3659 (defun simpargs1 (a vestigial c
)
3660 (declare (ignore vestigial
))
3664 (if (not (atom x
)) (car x
)))
3666 (defquote retlist
(&rest l
)
3668 (mapcar #'(lambda (z) (list '(mequal simp
) z
(meval z
))) l
)))
3671 (if (zerop c
) x
(nthkdr (kdr x
) (1- c
))))
