Add translation for plog (using TRANSLATE-WITH-FLONUM-OP)
[maxima.git] / src / compar.lisp
blobeec30b4002bb3350a088f2cca92d06a412194d62
1 ;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
2 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3 ;;; The data in this file contains enhancements. ;;;;;
4 ;;; ;;;;;
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 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
11 (in-package :maxima)
13 (macsyma-module compar)
15 (load-macsyma-macros mrgmac)
17 (declare-top (special success))
19 (defvar *debug-compar* nil
20 "Enables debugging code for this file.")
22 (defvar %initiallearnflag)
24 (defmvar complexsign nil
25 "If T, COMPAR attempts to work in a complex mode.
26 This scheme is only very partially developed at this time."
27 no-reset)
29 (defmvar $signbfloat t)
30 (defmvar $askexp)
31 (defmvar $assume_pos nil)
32 (defmvar $assume_pos_pred nil)
34 (defmvar factored nil)
36 ;; The *LOCAL-SIGNS* variable contains a list of facts that are local to the
37 ;; current evaluation. These are stored in the assume database (in the global
38 ;; context) by asksign1 when the user answers questions. A "top-level"
39 ;; evaluation is run by MEVAL* and that function calls CLEARSIGN when it
40 ;; finishes to discard them.
41 (defmvar *local-signs* nil)
43 (defmvar sign nil)
44 (defmvar minus nil)
45 (defmvar odds nil)
46 (defmvar evens nil)
48 (defvar $useminmax t)
50 ;; Remove this (nil'ed out) function after a while. We should be
51 ;; using POWER instead of POW.
52 #+nil
53 (defmacro pow (&rest x)
54 `(power ,@x))
56 (defun lmul (l)
57 (simplify (cons '(mtimes) l)))
59 (defun conssize (x)
60 (if (atom x)
62 (do ((x (cdr x) (cdr x))
63 (sz 1))
64 ((null x) sz)
65 (incf sz (1+ (conssize (car x)))))))
67 ;;; Functions for creating, activating, manipulating, and killing contexts
69 ;;; This "turns on" a context, making its facts visible.
71 (defmfun $activate (&rest args)
72 (dolist (c args)
73 (cond ((not (symbolp c)) (nc-err '$activate c))
74 ((member c (cdr $activecontexts) :test #'eq))
75 ((member c (cdr $contexts) :test #'eq)
76 (setq $activecontexts (mcons c $activecontexts))
77 (activate c))
78 (t (merror (intl:gettext "activate: no such context ~:M") c))))
79 '$done)
81 ;;; This "turns off" a context, keeping the facts, but making them invisible
83 (defmfun $deactivate (&rest args)
84 (dolist (c args)
85 (cond ((not (symbolp c)) (nc-err '$deactivate c))
86 ((member c (cdr $contexts) :test #'eq)
87 (setq $activecontexts ($delete c $activecontexts))
88 (deactivate c))
89 (t (merror (intl:gettext "deactivate: no such context ~:M") c))))
90 '$done)
92 ;;; This function prints out a list of the facts in the specified context.
93 ;;; No argument implies the current context.
95 (defmfun $facts (&optional (ctxt $context))
96 (if (member ctxt (cdr $contexts))
97 (facts1 ctxt)
98 (facts2 ctxt)))
100 (defun facts1 (con)
101 (contextmark)
102 (do ((l (zl-get con 'data) (cdr l))
103 (nl)
104 (u))
105 ((null l) (cons '(mlist) nl))
106 (when (visiblep (car l))
107 (setq u (intext (caaar l) (cdaar l)))
108 (unless (memalike u nl)
109 (push u nl)))))
111 ;; Look up facts from the database which contain expr. expr can be a symbol or
112 ;; a more general expression.
113 (defun facts2 (expr)
114 (labels ((among (x l)
115 (cond ((null l) nil)
116 ((atom l) (eq x l))
117 ((alike1 x l) t)
119 (do ((ll (cdr l) (cdr ll)))
120 ((null ll) nil)
121 (if (among x (car ll)) (return t)))))))
122 (do ((facts (cdr ($facts $context)) (cdr facts))
123 (ans))
124 ((null facts) (return (cons '(mlist) (reverse ans))))
125 (when (or (among expr (cadar facts))
126 (among expr (caddar facts)))
127 (push (car facts) ans)))))
129 (defun intext (rel body)
130 (setq body (mapcar #'doutern body))
131 (cond ((eq 'kind rel) (cons '($kind) body))
132 ((eq 'par rel) (cons '($par) body))
133 ((eq 'mgrp rel) (cons '(mgreaterp) body))
134 ((eq 'mgqp rel) (cons '(mgeqp) body))
135 ((eq 'meqp rel) (cons '($equal) body))
136 ((eq 'mnqp rel) (list '(mnot) (cons '($equal) body)))))
138 (defprop $context asscontext assign)
140 ;;; This function switches contexts, creating one if necessary.
142 (defun asscontext (xx y)
143 (declare (ignore xx))
144 (cond ((not (symbolp y)) (nc-err "context assignment" y))
145 ((member y $contexts :test #'eq) (setq context y $context y))
146 (t ($newcontext y))))
148 ;;; This function actually creates a context whose subcontext is $GLOBAL.
149 ;;; It also switches contexts to the newly created one.
150 ;;; If no argument supplied, then invent a name via gensym and use that.
152 (defmfun $newcontext (&rest args)
153 (if (null args)
154 ($newcontext ($gensym "context")) ;; make up a name and try again
155 (if (> (length args) 1)
156 (merror "newcontext: found more than one argument.")
157 (let ((x (first args)))
158 (cond
159 ((not (symbolp x)) (nc-err '$newcontext x))
160 ((member x $contexts :test #'eq)
161 (mtell (intl:gettext "newcontext: context ~M already exists.") x) nil)
163 (setq $contexts (mcons x $contexts))
164 (putprop x '($global) 'subc)
165 (setq context x $context x)))))))
167 ;;; This function creates a supercontext. If given one argument, it
168 ;;; makes the current context be the subcontext of the argument. If
169 ;;; given more than one argument, the first is assumed the name of the
170 ;;; supercontext and the rest are the subcontexts.
171 ;;; If no arguments supplied, then invent a name via gensym and use that.
173 (defmfun $supcontext (&rest x)
174 (cond ((null x) ($supcontext ($gensym "context"))) ;; make up a name and try again
175 ((caddr x) (merror (intl:gettext "supcontext: found more than two arguments.")))
176 ((not (symbolp (car x))) (nc-err '$supcontext (car x)))
177 ((member (car x) $contexts :test #'eq)
178 (merror (intl:gettext "supcontext: context ~M already exists.") (car x)))
179 ((and (cadr x) (not (member (cadr x) $contexts :test #'eq)))
180 (merror (intl:gettext "supcontext: no such context ~M") (cadr x)))
181 (t (setq $contexts (mcons (car x) $contexts))
182 (putprop (car x) (ncons (or (cadr x) $context)) 'subc)
183 (setq context (car x) $context (car x)))))
185 ;;; This function kills a context or a list of contexts
187 (defmfun $killcontext (&rest args)
188 (dolist (c args)
189 (if (symbolp c)
190 (killcontext c)
191 (nc-err '$killcontext c)))
192 (if (and (= (length args) 1) (eq (car args) '$global))
193 '$not_done
194 '$done))
196 (defun killallcontexts ()
197 (mapcar #'killcontext (cdr $contexts))
198 (setq $context '$initial context '$initial current '$initial
199 $contexts '((mlist) $initial $global) dobjects ())
200 ;;The DB variables
201 ;;conmark, conunmrk, conindex, connumber, and contexts
202 ;;concern garbage-collectible contexts, and so we're
203 ;;better off not resetting them.
204 (defprop $global 1 cmark) (defprop $initial 1 cmark)
205 (defprop $initial ($global) subc))
207 (defun killcontext (x)
208 (cond ((not (member x $contexts :test #'eq))
209 (mtell (intl:gettext "killcontext: no such context ~M.") x))
210 ((eq x '$global) '$global)
211 ((eq x '$initial)
212 (mapc #'remov (zl-get '$initial 'data))
213 (remprop '$initial 'data)
214 '$initial)
215 ((and (not (eq $context x)) (contextmark) (< 0 (zl-get x 'cmark)))
216 (mtell (intl:gettext "killcontext: context ~M is currently active.") x))
217 (t (if (member x $activecontexts)
218 ;; Context is on the list of active contexts. The test above
219 ;; checks for active contexts, but it seems not to work in all
220 ;; cases. So deactivate the context at this place to remove it
221 ;; from the list of active contexts before it is deleted.
222 ($deactivate x))
223 (setq $contexts ($delete x $contexts))
224 (cond ((and (eq x $context)
225 (equal ;;replace eq ?? wfs
226 (zl-get x 'subc) '($global)))
227 (setq $context '$initial)
228 (setq context '$initial))
229 ((eq x $context)
230 (setq $context (car (zl-get x 'subc)))
231 (setq context (car (zl-get x 'subc)))))
232 (killc x)
233 x)))
235 (defun nc-err (fn x)
236 (merror (intl:gettext "~M: context name must be a symbol; found ~M") fn x))
238 ;; Simplification and evaluation of boolean expressions
240 ;; Simplification of boolean expressions:
242 ;; and and or are declared nary. The sole effect of this is to allow Maxima to
243 ;; flatten nested expressions, e.g., a and (b and c) => a and b and c
244 ;; (The nary declaration does not make and and or commutative, and and and or
245 ;; are not otherwise declared commutative.)
247 ;; and: if any argument simplifies to false, return false
248 ;; otherwise omit arguments which simplify to true and simplify others
249 ;; if only one argument remains, return it
250 ;; if none remain, return true
252 ;; or: if any argument simplifies to true, return true
253 ;; otherwise omit arguments which simplify to false and simplify others
254 ;; if only one argument remains, return it
255 ;; if none remain, return false
257 ;; not: if argument simplifies to true / false, return false / true
258 ;; otherwise reverse sense of comparisons (if argument is a comparison)
259 ;; otherwise return not <simplified argument>
261 ;; Evaluation (MEVAL) of boolean expressions:
262 ;; same as simplification except evaluating (MEVALP) arguments instead of simplifying
263 ;; When prederror = true, complain if expression evaluates to something other than T / NIL
264 ;; (otherwise return unevaluated boolean expression)
266 ;; Evaluation (MEVALP) of boolean expressions:
267 ;; same as simplification except evaluating (MEVALP) arguments instead of simplifying
268 ;; When prederror = true, complain if expression evaluates to something other than T / NIL
269 ;; (otherwise return unevaluated boolean expression)
271 ;; Simplification of "is" expressions:
272 ;; if argument simplifies to true/false, return true/false
273 ;; otherwise return is (<simplified argument>)
275 ;; Evaluation of "is" expressions:
276 ;; if argument evaluates to true/false, return true/false
277 ;; otherwise return unknown if prederror = false, else trigger an error
279 ;; Simplification of "maybe" expressions:
280 ;; if argument simplifies to true/false, return true/false
281 ;; otherwise return maybe (<simplified expression>)
283 ;; Evaluation of "maybe" expressions:
284 ;; if argument evaluates to true/false, return true/false
285 ;; otherwise return unknown
287 (defprop $is simp-$is operators)
288 (defprop %is simp-$is operators)
289 (defprop $maybe simp-$is operators)
290 (defprop %maybe simp-$is operators)
292 ; I'VE ASSUMED (NULL Z) => SIMPLIFIY ARGUMENTS
293 ; SAME WITH SIMPCHECK (SRC/SIMP.LISP)
294 ; SAME WITH TELLSIMP-GENERATED SIMPLIFICATION FUNCTIONS
295 ; SAME WITH SIMPLIFICATION OF %SIN
296 ; PRETTY SURE I'VE SEEN OTHER EXAMPLES AS WELL
297 ; Z SEEMS TO SIGNIFY "ARE THE ARGUMENTS SIMPLIFIED YET"
299 (defun maybe-simplifya (x z)
300 (if z x (simplifya x z)))
302 (defun maybe-simplifya-protected (x z)
303 (let ((errcatch t) ($errormsg nil))
304 (ignore-errors (maybe-simplifya x z) x)))
306 (defun simp-$is (x yy z)
307 (declare (ignore yy))
308 (let ((a (maybe-simplifya (cadr x) z)))
309 (if (or (eq a t) (eq a nil))
311 `((,(caar x) simp) ,a))))
313 (defmspec $is (form)
314 (unless (= 1 (length (rest form)))
315 (merror (intl:gettext "is() expects a single argument. Found ~A")
316 (length (rest form))))
317 (destructuring-bind (answer patevalled)
318 (mevalp1 (cadr form))
319 (cond ((member answer '(t nil) :test #'eq) answer)
320 ;; I'D RATHER HAVE ($PREDERROR ($THROW `(($PREDERROR) ,PATEVALLED))) HERE
321 ($prederror (pre-err patevalled))
322 (t '$unknown))))
324 (defmspec $maybe (form)
325 (let* ((pat (cadr form))
326 (x (let (($prederror nil)) (mevalp1 pat)))
327 (ans (car x)))
328 (if (member ans '(t nil) :test #'eq)
330 '$unknown)))
332 (defun is (pred)
333 (let (($prederror t))
334 (mevalp pred)))
336 ; The presence of OPERS tells SIMPLIFYA to call OPER-APPLY,
337 ; which calls NARY1 to flatten nested "and" and "or" expressions
338 ; (due to $NARY property of MAND and MOR, declared elsewhere).
340 (putprop 'mand t 'opers)
341 (putprop 'mor t 'opers)
343 (putprop 'mnot 'simp-mnot 'operators)
344 (putprop 'mand 'simp-mand 'operators)
345 (putprop 'mor 'simp-mor 'operators)
347 (defun simp-mand (x yy z)
348 (declare (ignore yy))
349 (do ((l (cdr x) (cdr l))
351 (simplified))
352 ((null l)
353 (cond ((= (length simplified) 0) t)
354 ((= (length simplified) 1) (car simplified))
355 (t (cons '(mand simp) (reverse simplified)))))
356 (setq a (maybe-simplifya (car l) z))
357 (cond ((null a) (return nil))
358 ((eq a '$unknown) (unless (member '$unknown simplified :test #'eq) (push a simplified)))
359 ((not (member a '(t nil) :test #'eq)) (push a simplified)))))
361 (defun simp-mor (x yy z)
362 (declare (ignore yy))
363 (do ((l (cdr x) (cdr l))
365 (simplified))
366 ((null l)
367 (cond ((= (length simplified) 0) nil)
368 ((= (length simplified) 1) (car simplified))
369 (t (cons '(mor simp) (reverse simplified)))))
370 (setq a (maybe-simplifya (car l) z))
371 (cond ((eq a t) (return t))
372 ((eq a '$unknown) (unless (member '$unknown simplified :test #'eq) (push a simplified)))
373 ((not (member a '(t nil) :test #'eq)) (push a simplified)))))
375 ; ALSO CUT STUFF ABOUT NOT EQUAL => NOTEQUAL AT TOP OF ASSUME
377 (defun simp-mnot (x yy z)
378 (declare (ignore yy))
379 (let ((arg (maybe-simplifya (cadr x) z)))
380 (if (atom arg)
381 (cond ((or (eq arg t) (eq arg '$true))
382 nil)
383 ((or (eq arg nil) (eq arg '$false))
385 ((eq arg '$unknown)
386 '$unknown)
387 (t `((mnot simp) ,arg)))
388 (let ((arg-op (caar arg)) (arg-arg (cdr arg)))
389 ;;(setq arg-arg (mapcar #'(lambda (a) (maybe-simplifya a z)) arg-arg))
390 (cond ((eq arg-op 'mlessp)
391 (simplify `((mgeqp) ,@arg-arg)))
392 ((eq arg-op 'mleqp)
393 (simplify `((mgreaterp) ,@arg-arg)))
394 ((eq arg-op 'mequal)
395 (simplify `((mnotequal) ,@arg-arg)))
396 ((eq arg-op '$equal)
397 (simplify `(($notequal) ,@arg-arg)))
398 ((eq arg-op 'mnotequal)
399 (simplify `((mequal) ,@arg-arg)))
400 ((eq arg-op '$notequal)
401 (simplify `(($equal) ,@arg-arg)))
402 ((eq arg-op 'mgeqp)
403 (simplify `((mlessp) ,@arg-arg)))
404 ((eq arg-op 'mgreaterp)
405 (simplify `((mleqp) ,@arg-arg)))
406 ((eq arg-op 'mnot)
407 (car arg-arg))
408 ;; Distribute negation over conjunction and disjunction;
409 ;; analogous to '(- (a + b)) --> - a - b.
410 ((eq arg-op 'mand)
411 (let ((L (mapcar #'(lambda (e) `((mnot) ,e)) arg-arg)))
412 (simplifya `((mor) ,@L) nil)))
413 ((eq arg-op 'mor)
414 (let ((L (mapcar #'(lambda (e) `((mnot) ,e)) arg-arg)))
415 (simplifya `((mand) ,@L) nil)))
416 (t `((mnot simp) ,arg)))))))
418 ;; =>* N.B. *<=
419 ;; The translator depends on some stuff in here.
420 ;; Check it out in the transl module ACALL before proceeding.
422 (defun mevalp (pat)
423 (let* ((x (mevalp1 pat))
424 (ans (car x))
425 (patevalled (cadr x)))
426 (cond ((member ans '(t ()) :test #'eq) ans)
427 ;; I'D RATHER HAVE ($PREDERROR ($THROW `(($PREDERROR) ,PATEVALLED))) HERE
428 ($prederror (pre-err patevalled))
429 (t (or patevalled ans)))))
431 (defun mevalp1 (pat)
432 (let (patevalled ans)
433 (setq ans
434 (cond ((and (not (atom pat))
435 (member (caar pat) '(mnot mand mor) :test #'eq))
436 (cond ((eq 'mnot (caar pat)) (is-mnot #'mevalp (cadr pat)))
437 ((eq 'mand (caar pat)) (is-mand #'mevalp (cdr pat)))
438 (t (is-mor #'mevalp (cdr pat)))))
439 ((atom (setq patevalled (specrepcheck (meval pat))))
440 patevalled)
441 ((member (caar patevalled) '(mnot mand mor) :test #'eq)
442 (return-from mevalp1 (mevalp1 patevalled)))
444 (mevalp2 patevalled
445 (caar patevalled)
446 (cadr patevalled)
447 (caddr patevalled)))))
448 (list ans patevalled)))
450 (defun mevalp2 (patevalled pred arg1 arg2)
451 (cond ((eq 'mequal pred) (like arg1 arg2))
452 ((eq '$equal pred) (meqp arg1 arg2))
453 ((eq 'mnotequal pred) (not (like arg1 arg2)))
454 ((eq '$notequal pred) (mnqp arg1 arg2))
455 ((eq 'mgreaterp pred) (mgrp arg1 arg2))
456 ((eq 'mlessp pred) (mgrp arg2 arg1))
457 ((eq 'mgeqp pred) (mgqp arg1 arg2))
458 ((eq 'mleqp pred) (mgqp arg2 arg1))
459 (t (isp (munformat patevalled)))))
461 (defun pre-err (pat)
462 (merror (intl:gettext "Unable to evaluate predicate ~M") pat))
464 (defun is-mnot (pred-eval pred)
465 (setq pred (funcall pred-eval pred))
466 (cond ((eq t pred) nil)
467 ((not pred))
468 (t (pred-reverse pred))))
470 (defun pred-reverse (pred)
471 (take '(mnot) pred))
473 (defun is-mand (pred-eval pl)
474 (do ((dummy)
475 (npl))
476 ((null pl) (cond ((null npl))
477 ((null (cdr npl)) (car npl))
478 (t (cons '(mand) (nreverse npl)))))
479 (setq dummy (funcall pred-eval (car pl)) pl (cdr pl))
480 (cond ((eq t dummy))
481 ((null dummy) (return nil))
482 (t (push dummy npl)))))
484 (defun is-mor (pred-eval pl)
485 (do ((dummy)
486 (npl))
487 ((null pl) (cond ((null npl) nil)
488 ((null (cdr npl)) (car npl))
489 (t (cons '(mor) (nreverse npl)))))
490 (setq dummy (funcall pred-eval (car pl)) pl (cdr pl))
491 (cond ((eq t dummy) (return t))
492 ((null dummy))
493 (t (push dummy npl)))))
495 (defmspec $assume (x)
496 (setq x (cdr x))
497 (do ((nl)) ((null x) (cons '(mlist) (nreverse nl)))
498 (cond ((atom (car x)) (push (assume (meval (car x))) nl))
499 ((eq 'mand (caaar x))
500 (mapc #'(lambda (l) (push (assume (meval l)) nl))
501 (cdar x)))
502 ((eq 'mnot (caaar x))
503 (push (assume (meval (pred-reverse (cadar x)))) nl))
504 ((eq 'mor (caaar x))
505 (merror (intl:gettext "assume: argument cannot be an 'or' expression; found ~M") (car x)))
506 ((eq (caaar x) 'mequal)
507 (merror (intl:gettext "assume: argument cannot be an '=' expression; found ~M~%assume: maybe you want 'equal'.") (car x)))
508 ((eq (caaar x) 'mnotequal)
509 (merror (intl:gettext "assume: argument cannot be a '#' expression; found ~M~%assume: maybe you want 'not equal'.") (car x)))
510 (t (push (assume (meval (car x))) nl)))
511 (setq x (cdr x))))
513 (defun assume (pat)
514 (if (and (not (atom pat))
515 (eq (caar pat) 'mnot)
516 (eq (caaadr pat) '$equal))
517 (setq pat `(($notequal) ,@(cdadr pat))))
518 (let ((dummy (let ($assume_pos) (car (mevalp1 pat)))))
519 (cond ((eq dummy t) '$redundant)
520 ((null dummy) '$inconsistent)
521 ((atom dummy) '$meaningless)
522 (t (learn pat t)))))
524 (defun learn (pat flag)
525 (cond ((atom pat))
526 ;; Check for abs function in pattern.
527 ((and (not limitp)
528 (learn-abs pat flag)))
529 ;; Check for constant expression in pattern.
530 ((and (not limitp)
531 (learn-numer pat flag)))
532 ((zl-get (caar pat) (if flag 'learn 'unlearn))
533 (funcall (zl-get (caar pat) (if flag 'learn 'unlearn)) pat))
534 ((eq (caar pat) 'mgreaterp) (daddgr flag (sub (cadr pat) (caddr pat))))
535 ((eq (caar pat) 'mgeqp) (daddgq flag (sub (cadr pat) (caddr pat))))
536 ((member (caar pat) '(mequal $equal) :test #'eq)
537 (daddeq flag (sub (cadr pat) (caddr pat))))
538 ((member (caar pat) '(mnotequal $notequal) :test #'eq)
539 (daddnq flag (sub (cadr pat) (caddr pat))))
540 ((eq (caar pat) 'mleqp) (daddgq flag (sub (caddr pat) (cadr pat))))
541 ((eq (caar pat) 'mlessp) (daddgr flag (sub (caddr pat) (cadr pat))))
542 (flag (true* (munformat pat)))
544 (cond
545 ((eq (caar pat) '$kind)
546 (unkind (second pat) (third pat)))
547 (t (untrue (munformat pat))))
548 pat)))
550 ;;; When abs(x)<a is in the pattern, where a is a positive expression,
551 ;;; then learn x<a and -x<a too. The additional facts are put into the context
552 ;;; '$learndata, if the current context is user context 'initial
554 (defun learn-abs (pat flag)
555 (let (tmp)
556 (when (and (setq tmp (isinop pat 'mabs))
557 (or (and (member (caar pat) '(mlessp mleqp))
558 (isinop (cadr pat) 'mabs)
559 (member ($sign (caddr pat)) '($pos $pz)))
560 (and (member (caar pat) '(mgreaterp mgeqp))
561 (member ($sign (cadr pat)) '($pos $pz))
562 (isinop (caddr pat) 'mabs))))
563 (let ((oldcontext context))
564 (if (eq oldcontext '$initial)
565 (asscontext nil '$learndata)) ; switch to context '$learndata
566 ; learn additional facts
567 (learn ($substitute (cadr tmp) tmp pat) flag)
568 (learn ($substitute (mul -1 (cadr tmp)) tmp pat) flag)
569 (when (eq oldcontext '$initial)
570 (asscontext nil oldcontext) ; switch back to context on entry
571 ($activate '$learndata)))) ; context '$learndata is active
572 nil))
574 ;;; The value of a constant expression which can be numerically evaluated is
575 ;;; put into the context '$learndata.
577 (defun learn-numer (pat flag)
578 (let (dum expr patnew)
579 (do ((x (cdr pat) (cdr x)))
580 ((null x) (setq patnew (reverse patnew)))
581 (setq dum (constp (car x))
582 expr (car x))
583 (cond ((or (numberp (car x))
584 (ratnump (car x))))
585 ((eq dum 'bigfloat)
586 (if (prog2
587 (setq dum ($bfloat (car x)))
588 ($bfloatp dum))
589 (setq expr dum)))
590 ((eq dum 'float)
591 (if (and (setq dum (numer (car x)))
592 (numberp dum))
593 (setq expr dum)))
594 ((and (member dum '(numer symbol) :test #'eq)
595 (prog2
596 (setq dum (numer (car x)))
597 (or (null dum)
598 (and (numberp dum)
599 (prog2
600 (setq expr dum)
601 (< (abs dum) 1.0e-6))))))
602 (cond ($signbfloat
603 (and (setq dum ($bfloat (car x)))
604 ($bfloatp dum)
605 (setq expr dum))))))
606 (setq patnew (cons expr patnew)))
607 (setq patnew (cons (car pat) patnew))
608 (when (and (not (alike (cdr pat) (cdr patnew)))
609 (or (not (mnump (cadr patnew))) ; not both sides of the
610 (not (mnump (caddr patnew))))) ; relation can be number
611 (let ((oldcontext $context))
612 (if (eq oldcontext '$initial)
613 (asscontext nil '$learndata)) ; switch to context '$learndata
614 (learn patnew flag) ; learn additional fact
615 (when (eq oldcontext '$initial)
616 (asscontext nil oldcontext) ; switch back to context on entry
617 ($activate '$learndata)))) ; context '$learndata is active
618 nil))
620 (defmspec $forget (x)
621 (setq x (cdr x))
622 (do ((nl))
623 ((null x) (cons '(mlist) (nreverse nl)))
624 (cond ((atom (car x)) (push (forget (meval (car x))) nl))
625 ((eq 'mand (caaar x))
626 (mapc #'(lambda (l) (push (forget (meval l)) nl)) (cdar x)))
627 ((eq 'mnot (caaar x))
628 (push (forget (meval (pred-reverse (cadar x)))) nl))
629 ((eq 'mor (caaar x))
630 (merror (intl:gettext "forget: argument cannot be an 'or' expression; found ~M") (car x)))
631 (t (push (forget (meval (car x))) nl)))
632 (setq x (cdr x))))
634 (defun forget (pat)
635 (cond (($listp pat)
636 (cons '(mlist simp) (mapcar #'forget1 (cdr pat))))
637 (t (forget1 pat))))
639 (defun forget1 (pat)
640 (cond ((and (not (atom pat))
641 (eq (caar pat) 'mnot)
642 (eq (caaadr pat) '$equal))
643 (setq pat `(($notequal) ,@(cdadr pat)))))
644 (learn pat nil))
646 (defun restore-facts (factl) ; used by SAVE
647 (dolist (fact factl)
648 (cond ((eq (caar fact) '$kind)
649 (declarekind (cadr fact) (caddr fact))
650 (add2lnc (getop (cadr fact)) $props))
651 ((eq (caar fact) '$par))
652 (t (assume fact)))))
654 (defmacro compare (a b)
655 `(sign1 (sub* ,a ,b)))
657 (defun maximum (l)
658 (maximin l '$max))
660 (defun minimum (l)
661 (maximin l '$min))
663 (defmspec mand (form)
664 (setq form (cdr form))
665 (do ((l form (cdr l))
667 (unevaluated))
668 ((null l)
669 (cond ((= (length unevaluated) 0) t)
670 ((= (length unevaluated) 1) (car unevaluated))
671 (t (cons '(mand) (reverse unevaluated)))))
672 (setq x (mevalp (car l)))
673 (cond ((null x) (return nil))
674 ((not (member x '(t nil) :test #'eq)) (push x unevaluated)))))
676 (defmspec mor (form)
677 (setq form (cdr form))
678 (do ((l form (cdr l))
680 (unevaluated))
681 ((null l)
682 (cond ((= (length unevaluated) 0) nil)
683 ((= (length unevaluated) 1) (car unevaluated))
684 (t (cons '(mor) (reverse unevaluated)))))
685 (setq x (mevalp (car l)))
686 (cond ((eq x t) (return t))
687 ((not (member x '(t nil) :test #'eq)) (push x unevaluated)))))
689 (defmspec mnot (form)
690 (setq form (cdr form))
691 (let ((x (mevalp (car form))))
692 (if (member x '(t nil) :test #'eq)
693 (not x)
694 `((mnot) ,x))))
696 ;;;Toplevel functions- $askequal, $asksign, and $sign.
697 ;;;Switches- LIMITP If TRUE $ASKSIGN and $SIGN will look for special
698 ;;; symbols such as EPSILON, $INF, $MINF and attempt
699 ;;; to do the correct thing. In addition calls to
700 ;;; $REALPART and $IMAGPART are made to assure that
701 ;;; the expression is real.
703 ;;; if NIL $ASKSIGN and $SIGN assume the expression
704 ;;; given is real unless it contains an $%I, in which
705 ;;; case they call $RECTFORM.
707 (setq limitp nil)
709 (defmfun $askequal (a b)
710 (let ((answer (meqp (sratsimp a) (sratsimp b)))) ; presumably handles mbags and extended reals.
711 (cond ((eq answer t) '$yes)
712 ((eq answer nil) '$no)
714 (setq answer (retrieve `((mtext) ,(intl:gettext "Is ") ,a ,(intl:gettext " equal to ") ,b ,(intl:gettext "?")) nil))
715 (cond ((member answer '($no |$n| |$N|) :test #'eq)
716 (tdpn (sub b a))
717 '$no)
718 ((member answer '($yes |$y| |$Y|) :test #'eq)
719 (tdzero (sub a b))
720 '$yes)
722 (mtell (intl:gettext "Acceptable answers are yes, y, Y, no, n, N. ~%"))
723 ($askequal a b)))))))
725 (defmfun $asksign (exp)
726 (let (sign minus odds evens factored)
727 (asksign01 (cond (limitp (restorelim exp))
728 ((among '$%i exp) ($rectform exp))
729 (t exp)))))
731 (defun asksign-p-or-n (e)
732 (unwind-protect (prog2
733 (assume `(($notequal) ,e 0))
734 ($asksign e))
735 (forget `(($notequal) ,e 0))))
737 (defun asksign01 (a)
738 (let ((e (sign-prep a)))
739 (cond ((eq e '$pnz) '$pnz)
740 ((member (setq e (asksign1 e)) '($pos $neg) :test #'eq) e)
741 (limitp (eps-sign a))
742 (t '$zero))))
744 ;; csign returns t if x appears to be complex.
745 ;; Else, it returns the sign.
746 (defun csign (x)
747 (or (not (free x '$%i))
748 (let (sign-imag-errp limitp) (catch 'sign-imag-err ($sign x)))))
750 ;;; $csign works like $sign but switches the sign-functions into a complex
751 ;;; mode. In complex mode complex and imaginary expressions give the results
752 ;;; imaginary or complex.
754 (defmfun $csign (z)
755 (let ((*complexsign* t)
756 (limitp nil))
757 ($sign z)))
759 (defmfun $sign (x)
760 (let ((x (specrepcheck x))
761 sign minus odds evens factored)
762 (sign01 (cond (limitp (restorelim x))
763 (*complexsign*
764 ;; No rectform in Complex mode. Rectform ask unnecessary
765 ;; questions about complex expressions and can not handle
766 ;; imaginary expressions completely. Thus $csign can not
767 ;; handle something like (1+%i)*(1-%i) which is real.
768 ;; After improving rectform, we can change this. (12/2008)
769 (when *debug-compar*
770 (format t "~&$SIGN with ~A~%" x))
772 ((not (free x '$%i)) ($rectform x))
773 (t x)))))
775 (defun sign01 (a)
776 (let ((e (sign-prep a)))
777 (cond ((eq e '$pnz) '$pnz)
778 (t (setq e (sign1 e))
779 (if (and limitp (eq e '$zero)) (eps-sign a) e)))))
781 ;;; Preparation for asking questions from DEFINT or LIMIT.
782 (defun sign-prep (x)
783 (if limitp
784 (destructuring-let (((rpart . ipart) (trisplit x)))
785 (cond ((and (equal (sratsimp ipart) 0)
786 (free rpart '$infinity))
787 (setq x (nmr (sratsimp rpart)))
788 (if (free x 'prin-inf)
790 ($limit x 'prin-inf '$inf '$minus)))
791 (t '$pnz))) ; Confess ignorance if COMPLEX.
794 ;; don't ask about internal variables created by gruntz
795 (defun has-int-symbols (e)
796 (cond ((and (symbolp e) (get e 'internal))
798 ((atom e) nil)
799 (t (or (has-int-symbols (car e))
800 (has-int-symbols (cdr e))))))
802 ;;; Do substitutions for special symbols.
803 (defun nmr (a)
804 (unless (free a '$zeroa) (setq a ($limit a '$zeroa 0 '$plus)))
805 (unless (free a '$zerob) (setq a ($limit a '$zerob 0 '$minus)))
806 (unless (free a 'z**) (setq a ($limit a 'z** 0 '$plus)))
807 (unless (free a '*z*) (setq a ($limit a '*z* 0 '$plus)))
808 (unless (free a 'epsilon) (setq a ($limit a 'epsilon 0 '$plus)))
809 a) ;;; Give A back.
811 ;;; Get the sign of EPSILON-like terms. Could be made MUCH hairier.
812 (defun eps-sign (b)
813 (let (temp1 temp2 temp3 free1 free2 free3 limitp)
814 ;; unset limitp to prevent infinite recursion
815 (cond ((not (free b '$zeroa))
816 (setq temp1 (eps-coef-sign b '$zeroa)))
817 (t (setq free1 t)))
818 (cond ((not (free b '$zerob))
819 (setq temp2 (eps-coef-sign b '$zerob)))
820 (t (setq free2 t)))
821 (cond ((not (free b 'epsilon))
822 (setq temp3 (eps-coef-sign b 'epsilon)))
823 (t (setq free3 t)))
824 (cond ((and free1 free2 free3) '$zero)
825 ((or (not (null temp1)) (not (null temp2)) (not (null temp3)))
826 (cond ((and (null temp1) (null temp2)) temp3)
827 ((and (null temp2) (null temp3)) temp1)
828 ((and (null temp1) (null temp3)) temp2)
829 (t (merror (intl:gettext "asksign: internal error."))))))))
831 (defun eps-coef-sign (exp epskind)
832 (let ((eps-power ($lopow exp epskind)) eps-coef)
833 (cond ((and (not (equal eps-power 0))
834 (not (equal (setq eps-coef (ratcoeff exp epskind eps-power))
836 (eq (ask-integer eps-power '$integer) '$yes))
837 (cond ((eq (ask-integer eps-power '$even) '$yes)
838 ($sign eps-coef))
839 ((eq (ask-integer eps-power '$odd) '$yes)
840 (setq eps-coef ($sign eps-coef))
841 (cond ((or (and (eq eps-coef '$pos)
842 (or (eq epskind 'epsilon)
843 (eq epskind '$zeroa)))
844 (and (eq eps-coef '$neg)
845 (or (alike epskind (mul2* -1 'epsilon))
846 (eq epskind '$zerob))))
847 '$pos)
848 (t '$neg)))
849 (t (merror (intl:gettext "sign or asksign: insufficient information.")))))
850 (t (let ((deriv (sdiff exp epskind)) deriv-sign)
851 (cond ((not (eq (setq deriv-sign ($sign deriv)) '$zero))
852 (total-sign epskind deriv-sign))
853 ((not
854 (eq (let ((deriv (sdiff deriv epskind)))
855 (setq deriv-sign ($sign deriv)))
856 '$zero))
857 deriv-sign)
858 (t (merror (intl:gettext "sign or asksign: insufficient data.")))))))))
860 ;;; The above code does a partial Taylor series analysis of something
861 ;;; that isn't a polynomial.
863 (defun total-sign (epskind factor-sign)
864 (cond ((or (eq epskind '$zeroa) (eq epskind 'epsilon))
865 (cond ((eq factor-sign '$pos) '$pos)
866 ((eq factor-sign '$neg) '$neg)
867 ((eq factor-sign '$zero) '$zero)))
868 ((eq epskind '$zerob)
869 (cond ((eq factor-sign '$pos) '$neg)
870 ((eq factor-sign '$neg) '$pos)
871 ((eq factor-sign '$zero) '$zero)))))
873 (defun asksign (x)
874 (setq x ($asksign x))
875 (cond ((eq '$pos x) '$positive)
876 ((eq '$neg x) '$negative)
877 ((eq '$pnz x) '$pnz) ;COMPLEX expression encountered here.
878 (t '$zero)))
880 ;; When asksign is called on an expression that involves a symbol with the
881 ;; internal property and Maxima's sign function can only determine that the
882 ;; sign is pnz, nz, or similar, asksign1 returns '$pnz. This causes asksign01,
883 ;; and ultimately $asksign, to return zero. Lesson: be circumspect about calling
884 ;; asksign on an expression that might involve a symbol with the internal property.
885 (defun asksign1 ($askexp)
886 (let ($radexpand)
887 (sign1 $askexp))
888 (cond
889 ((member sign '($pos $neg $zero $imaginary) :test #'eq) sign)
890 ((has-int-symbols $askexp) '$pnz)
892 (let ((domain sign) (squared nil))
893 (cond
894 ((null odds)
895 (setq $askexp (lmul evens)
896 domain '$znz
897 squared t))
899 (if minus (setq sign (flip sign)))
900 (setq $askexp
901 (lmul (nconc odds (mapcar #'(lambda (l) (power l 2)) evens))))))
902 (setq sign (cdr (assol $askexp *local-signs*)))
903 (ensure-sign $askexp domain squared)))))
905 (defun match-sign (sgn domain expression squared)
906 "If SGN makes sense for DOMAIN store the result (see ENSURE-SIGN) and return
907 it. Otherwise, return NIL. If SQUARED is true, we are actually looking for the
908 sign of the square, so any negative results are converted to positive."
909 ;; The entries in BEHAVIOUR are of the form
910 ;; (MATCH DOMAINS REGISTRAR SIGN SIGN-SQ)
912 ;; The algorithm goes as follows:
914 ;; Look for SGN in MATCH. If found, use REGISTRAR to store SIGN for the
915 ;; expression and then return SIGN if SQUARED is false or SIGN-SQ if it is
916 ;; true.
917 (let* ((behaviour
918 '((($pos |$P| |$p| $positive) (nil $znz $pz $pn $pnz) tdpos $pos $pos)
919 (($neg |$N| |$n| $negative) (nil $znz $nz $pn $pnz) tdneg $neg $pos)
920 (($zero |$Z| |$z| 0 0.0) (nil $znz $pz $nz $pnz) tdzero $zero $zero)
921 (($pn $nonzero $nz $nonz $non0) ($znz) tdpn $pn $pos)))
922 (hit (find-if (lambda (bh)
923 (and (member sgn (first bh) :test #'equal)
924 (member domain (second bh) :test #'eq)))
925 behaviour)))
926 (when hit
927 (let ((registrar (third hit))
928 (found-sign (if squared (fifth hit) (fourth hit))))
929 (funcall registrar expression)
930 (setq sign
931 (if minus (flip found-sign) found-sign))))))
933 (defun ensure-sign (expr &optional domain squared)
934 "Try to determine the sign of EXPR. If DOMAIN is not one of the special values
935 described below, we try to tell whether EXPR is positive, negative or zero. It
936 can be more specialised ($pz => positive or zero; $nz => negative or zero; $pn
937 => positive or negative; $znz => zero or nonzero).
939 If SQUARED is true, then we're actually interested in the sign of EXPR^2. As
940 such, a nonzero sign should be regarded as positive.
942 When calling ENSURE-SIGN, set the special variable SIGN to the best current
943 guess for the sign of EXPR. The function returns the sign, calls one of (TDPOS
944 TDNEG TDZERO TDPN) to store it, and also sets SIGN."
945 (loop
946 (let ((new-sign (match-sign sign domain expr squared)))
947 (when new-sign (return new-sign)))
948 (setf sign (retrieve
949 (list '(mtext)
950 "Is " expr
951 (or (second
952 (assoc domain
953 '(($znz " zero or nonzero?")
954 ($pz " positive or zero?")
955 ($nz " negative or zero?")
956 ($pn " positive or negative?"))))
957 " positive, negative or zero?"))
958 nil))))
960 ;; During one evaluation phase asksign writes answers from the user into the
961 ;; global context '$initial. These facts are removed by clearsign after
962 ;; finishing the evaluation phase. clearsign is called from the top-level
963 ;; evaluation function meval*. The facts which have to be removed are stored
964 ;; in the global variable *local-signs*.
966 (defun clearsign ()
967 (let ((context '$initial))
968 (dolist (cons-pair *local-signs*)
969 (destructuring-bind (x . sgn) cons-pair
970 (cond
971 ((eq '$pos sgn) (daddgr nil x))
972 ((eq '$neg sgn) (daddgr nil (neg x)))
973 ((eq '$zero sgn) (daddeq nil x))
974 ((eq '$pn sgn) (daddnq nil x))
975 ((eq '$pz sgn) (daddgq nil x))
976 ((eq '$nz sgn) (daddgq nil (neg x))))))
977 (setf *local-signs* nil)))
979 (defun like (x y)
980 (alike1 (specrepcheck x) (specrepcheck y)))
982 (setf (get '$und 'sysconst) t)
983 (setf (get '$ind 'sysconst) t)
984 (setf (get '$zeroa 'sysconst) t)
985 (setf (get '$zerob 'sysconst) t)
987 ;; There have been some conversations about NaN on the list, but
988 ;; the issue hasn't been settled.
990 (defvar indefinites `($und $ind))
992 ;; Other than sums, products, and lambda forms, meqp knows nothing
993 ;; about dummy variables. Because of the way niceindices chooses names
994 ;; for the sum indices, it's necessary to locally assign a new value to
995 ;; niceindicespref.
997 (defun meqp-by-csign (z a b)
998 (let (($niceindicespref `((mlist) ,(gensym) ,(gensym) ,(gensym))))
999 (setq z ($niceindices z))
1000 (setq z (if ($constantp z) ($rectform z) (sratsimp z)))
1001 (let ((sgn ($csign z))
1002 (dunno `(($equal) ,a ,b)))
1003 (cond ((eq '$zero sgn) t)
1004 ((memq sgn '($pos $neg $pn)) nil)
1006 ;; previously checked also for (linearp z '$%i))
1007 ((memq sgn '($complex $imaginary))
1008 ;; We call trisplit here, which goes back to general evaluation and
1009 ;; could cause an infinite recursion. To make sure that doesn't
1010 ;; happen, use the with-safe-recursion macro.
1011 (handler-case
1012 (with-safe-recursion meqp-by-csign z
1013 (let* ((ri-parts (trisplit z))
1014 (rsgn ($csign (car ri-parts)))
1015 (isgn ($csign (cdr ri-parts))))
1016 (cond ((and (eq '$zero rsgn)
1017 (eq '$zero isgn)) t)
1019 ((or (memq rsgn '($neg $pos $pn))
1020 (memq isgn '($neg $pos $pn))) nil)
1022 (t dunno))))
1023 (unsafe-recursion () dunno)))
1025 (t dunno)))))
1027 ;; For each fact of the form equal(a,b) in the active context, do e : ratsubst(b,a,e).
1029 (defun equal-facts-simp (e)
1030 (let ((f (margs ($facts))))
1031 (dolist (fi f e)
1032 (if (op-equalp fi '$equal)
1033 (setq e ($ratsubst (nth 2 fi) (nth 1 fi) e))))))
1035 (defun maxima-declared-arrayp (x)
1036 (and
1037 (symbolp x)
1038 (mget x 'array)
1039 (get (mget x 'array) 'array)))
1041 (defun maxima-undeclared-arrayp (x)
1042 (and
1043 (symbolp x)
1044 (mget x 'hashar)
1045 (get (mget x 'hashar) 'array)))
1047 (defun meqp (a b)
1048 ;; Check for some particular types before falling into the general case.
1049 (cond ((stringp a)
1050 (and (stringp b) (equal a b)))
1051 ((stringp b) nil)
1052 ((arrayp a)
1053 (and (arrayp b) (array-meqp a b)))
1054 ((arrayp b) nil)
1055 ((maxima-declared-arrayp a)
1056 (and (maxima-declared-arrayp b) (maxima-declared-array-meqp a b)))
1057 ((maxima-declared-arrayp b) nil)
1058 ((maxima-undeclared-arrayp a)
1059 (and (maxima-undeclared-arrayp b) (maxima-undeclared-array-meqp a b)))
1060 ((maxima-undeclared-arrayp b) nil)
1062 (let ((z) (sign))
1063 (setq a (specrepcheck a))
1064 (setq b (specrepcheck b))
1065 (cond ((or (like a b)) (not (member a indefinites)))
1066 ((or (member a indefinites) (member b indefinites)
1067 (member a *infinities*) (member b *infinities*)) nil)
1068 ((and (symbolp a) (or (eq t a) (eq nil a) (get a 'sysconst))
1069 (symbolp b) (or (eq t b) (eq nil b) (get b 'sysconst))) nil)
1070 ((or (mbagp a) (mrelationp a) (mbagp b) (mrelationp b))
1071 (cond ((and (or (and (mbagp a) (mbagp b)) (and (mrelationp a) (mrelationp b)))
1072 (eq (mop a) (mop b)) (= (length (margs a)) (length (margs b))))
1073 (setq z (list-meqp (margs a) (margs b)))
1074 (if (or (eq z t) (eq z nil)) z `(($equal) ,a ,b)))
1075 (t nil)))
1076 ((and (op-equalp a 'lambda) (op-equalp b 'lambda)) (lambda-meqp a b))
1077 (($setp a) (set-meqp a b))
1078 ;; 0 isn't in the range of an exponential function.
1079 ((or (and (mexptp a) (not (eq '$minf (third a))) (zerop1 b) (eq t (mnqp (second a) 0)))
1080 (and (mexptp b) (not (eq '$minf (third b))) (zerop1 a) (eq t (mnqp (second b) 0))))
1081 nil)
1083 ;; DCOMPARE emits new stuff (via DINTERNP) into the assume database.
1084 ;; Let's avoid littering the database with numbers.
1085 ((and (mnump a) (mnump b)) (zerop1 (sub a b)))
1087 ;; lookup in assumption database
1088 ((and (dcompare a b) (eq '$zero sign))) ; dcompare sets sign
1089 ((memq sign '($pos $neg $pn)) nil)
1091 ;; if database lookup failed, apply all equality facts
1092 (t (meqp-by-csign (equal-facts-simp (sratsimp (sub a b))) a b)))))))
1094 ;; Two arrays are equal (according to MEQP)
1095 ;; if (1) they have the same dimensions,
1096 ;; and (2) their elements are MEQP.
1098 (defun array-meqp (p q)
1099 (and
1100 (equal (array-dimensions p) (array-dimensions q))
1101 (progn
1102 (dotimes (i (array-total-size p))
1103 (let ((z (let ($ratprint)
1104 (declare (special $ratprint))
1105 (meqp (row-major-aref p i) (row-major-aref q i)))))
1106 (cond ((eq z nil) (return-from array-meqp nil))
1107 ((eq z t))
1108 (t (return-from array-meqp `(($equal) ,p ,q))))))
1109 t)))
1111 (defun maxima-declared-array-meqp (p q)
1112 (array-meqp (get (mget p 'array) 'array) (get (mget q 'array) 'array)))
1114 (defun maxima-undeclared-array-meqp (p q)
1115 (and
1116 (alike1 (mfuncall '$arrayinfo p) (mfuncall '$arrayinfo q))
1117 (let ($ratprint)
1118 (declare (special $ratprint))
1119 (meqp ($listarray p) ($listarray q)))))
1121 (defun list-meqp (p q)
1122 (let ((z))
1123 (cond ((or (null p) (null q)) (and (null p) (null q)))
1125 (setq z (meqp (car p) (car q)))
1126 (cond ((eq z nil) nil)
1127 ((or (eq z '$unknown) (op-equalp z '$equal)) z)
1128 (t (list-meqp (cdr p) (cdr q))))))))
1130 (defun lambda-meqp (a b)
1131 (let ((z))
1132 (cond ((= (length (second a)) (length (second b)))
1133 (let ((x) (n ($length (second a))))
1134 (dotimes (i n (push '(mlist) x)) (push (gensym) x))
1135 (setq z (meqp (mfuncall '$apply a x) (mfuncall '$apply b x)))
1136 (if (or (eq t z) (eq nil z)) z `(($equal) ,a ,b))))
1137 (t nil))))
1139 (defun set-meqp (a b)
1140 (let ((aa (equal-facts-simp a))
1141 (bb (equal-facts-simp b)))
1142 (cond ((or (not ($setp bb))
1143 (and ($emptyp aa) (not ($emptyp bb)))
1144 (and ($emptyp bb) (not ($emptyp aa))))
1145 nil)
1146 ((and (= (length aa) (length bb))
1147 (every #'(lambda (p q) (eq t (meqp p q))) (margs aa) (margs bb))) t)
1148 ((set-not-eqp (margs aa) (margs bb)) nil)
1149 (t `(($equal ,a ,b))))))
1151 (defun set-not-eqp (a b)
1152 (catch 'done
1153 (dolist (ak a)
1154 (if (every #'(lambda (s) (eq nil (meqp ak s))) b) (throw 'done t)))
1155 (dolist (bk b)
1156 (if (every #'(lambda (s) (eq nil (meqp bk s))) a) (throw 'done t)))
1157 (throw 'done nil)))
1159 ;; MGRP-GENERAL applies when sign of A - B (possibly complex) makes sense.
1161 (defun mgrp-general (a b)
1162 (setq a (sub a b))
1163 (let ((sgn (csign a)))
1164 (cond ((eq sgn '$pos) t)
1165 ((eq sgn t) nil) ;; csign thinks a - b isn't real
1166 ((member sgn '($neg $zero $nz) :test #'eq) nil)
1167 (t `((mgreaterp) ,a 0)))))
1169 (defun mgrp (a b)
1170 (cond
1171 ((or (stringp a) (stringp b))
1172 (if (and (stringp a) (stringp b)) (not (null (string< b a)))
1173 (when $prederror
1174 (merror (intl:gettext "greater than: arguments are incomparable; found: ~:M, ~:M") a b))))
1175 (t (mgrp-general a b))))
1177 (defun mlsp (x y)
1178 (mgrp y x))
1180 ;; MGQP-GENERAL applies when sign of A - B (possibly complex) makes sense.
1182 (defun mgqp-general (a b)
1183 (setq a (sub a b))
1184 (let ((sgn (csign a)))
1185 (cond ((member sgn '($pos $zero $pz) :test #'eq) t)
1186 ((eq sgn t) nil) ;; csign thinks a - b isn't real
1187 ((eq sgn '$neg) nil)
1188 (t `((mgeqp) ,a 0)))))
1190 (defun mgqp (a b)
1191 (cond
1192 ((or (stringp a) (stringp b))
1193 (if (and (stringp a) (stringp b)) (not (null (string<= b a)))
1194 (when $prederror
1195 (merror (intl:gettext "greater than or equal: arguments are incomparable; found: ~:M, ~:M") a b))))
1196 (t (mgqp-general a b))))
1198 (defun mnqp (x y)
1199 (let ((b (meqp x y)))
1200 (cond ((eq b '$unknown) b)
1201 ((or (eq b t) (eq b nil)) (not b))
1202 (t `(($notequal) ,x ,y)))))
1204 (defun c-$pn (o e)
1205 (list '(mnot) (c-$zero o e)))
1207 (defun c-$zero (o e)
1208 (list '($equal) (lmul (nconc o e)) 0))
1210 (defun c-$pos (o e)
1211 (cond ((null o) (list '(mnot) (list '($equal) (lmul e) 0)))
1212 ((null e) (list '(mgreaterp) (lmul o) 0))
1213 (t (setq e (mapcar #'(lambda (l) (power l 2)) e))
1214 (list '(mgreaterp) (lmul (nconc o e)) 0))))
1216 (defun c-$pz (o e)
1217 (cond ((null o) (list '(mnot) (list '($equal) (lmul e) 0)))
1218 ((null e) (list '(mgeqp) (lmul o) 0))
1219 (t (setq e (mapcar #'(lambda (l) (power l 2)) e))
1220 (list '(mgeqp) (lmul (nconc o e)) 0))))
1222 (defun sign* (x)
1223 (let (sign minus odds evens)
1224 (sign1 x)))
1226 (defun infsimp* (e)
1227 (if (or ($mapatom e) (and (free e '$inf) (free e '$minf)))
1229 (infsimp e)))
1231 ;; Like WITH-COMPSPLT, but runs COMPSPLT-EQ instead
1232 (defmacro with-compsplt-eq ((lhs rhs x) &body forms)
1233 `(multiple-value-bind (,lhs ,rhs) (compsplt-eq ,x)
1234 ,@forms))
1236 ;; Call FORMS with LHS and RHS bound to the splitting of EXPR by COMPSPLT.
1237 (defmacro with-compsplt ((lhs rhs expr) &body forms)
1238 `(multiple-value-bind (,lhs ,rhs) (compsplt ,expr)
1239 ,@forms))
1241 (defun sign1 (x)
1242 (setq x (specrepcheck x))
1243 (setq x (infsimp* x))
1244 (when (and *complexsign* (atom x) (eq x '$infinity))
1245 ;; In Complex Mode the sign of infinity is complex.
1246 (when *debug-compar* (format t "~& in sign1 detect $infinity.~%"))
1247 (return-from sign1 '$complex))
1248 (if (member x '($und $ind $infinity) :test #'eq)
1249 (if limitp '$pnz (merror (intl:gettext "sign: sign of ~:M is undefined.") x)))
1250 (prog (dum exp)
1251 (setq dum (constp x) exp x)
1252 (cond ((or (numberp x) (ratnump x)))
1253 ((eq dum 'bigfloat)
1254 (if (prog2 (setq dum ($bfloat x)) ($bfloatp dum))
1255 (setq exp dum)))
1256 ((eq dum 'float)
1257 (if (and (setq dum (numer x)) (numberp dum)) (setq exp dum)))
1258 ((and (member dum '(numer symbol) :test #'eq)
1259 (prog2 (setq dum (numer x))
1260 (or (null dum)
1261 (and (numberp dum)
1262 (prog2 (setq exp dum)
1263 (< (abs dum) 1.0e-6))))))
1264 (cond ($signbfloat
1265 (and (setq dum ($bfloat x)) ($bfloatp dum) (setq exp dum)))
1266 (t (setq sign '$pnz evens nil odds (ncons x) minus nil)
1267 (return sign)))))
1268 (or (and (not (atom x)) (not (mnump x)) (equal x exp)
1269 (let (s o e m)
1270 (with-compsplt (lhs rhs x)
1271 (dcompare lhs rhs)
1272 (cond ((member sign '($pos $neg) :test #'eq))
1273 ((eq sign '$pnz) nil)
1274 (t (setq s sign o odds e evens m minus)
1275 (sign x)
1276 (if (not (strongp sign s))
1277 (if (and (eq sign '$pnz) (eq s '$pn))
1278 (setq sign s)
1279 (setq sign s odds o evens e minus m)))
1280 t)))))
1281 (sign exp))
1282 (return sign)))
1284 (defun numer (x)
1285 (let (($numer t) ; currently, no effect on $float, but proposed to
1286 ($ratprint nil)
1287 result)
1288 ;; Catch a Lisp error, if a floating point overflow occurs.
1289 (setq result (let ((errset nil)) (errset ($float x))))
1290 (if result (car result) nil)))
1292 (defun constp (x)
1293 (cond ((floatp x) 'float)
1294 ((numberp x) 'numer)
1295 ((symbolp x) (if (member x '($%pi $%e $%phi $%gamma $%catalan) :test #'eq) 'symbol))
1296 ((atom x) nil)
1297 ((eq (caar x) 'rat) 'numer)
1298 ((eq (caar x) 'bigfloat) 'bigfloat)
1299 ((specrepp x) (constp (specdisrep x)))
1300 (t (do ((l (cdr x) (cdr l)) (dum) (ans 'numer))
1301 ((null l) ans)
1302 (setq dum (constp (car l)))
1303 (cond ((eq dum 'float) (return 'float))
1304 ((eq dum 'numer))
1305 ((eq dum 'bigfloat) (setq ans 'bigfloat))
1306 ((eq dum 'symbol)
1307 (if (eq ans 'numer) (setq ans 'symbol)))
1308 (t (return nil)))))))
1310 (mapcar #'(lambda (s) (putprop (first s) (second s) 'sign-function))
1311 (list
1312 (list 'mtimes 'sign-mtimes)
1313 (list 'mplus 'sign-mplus)
1314 (list 'mexpt 'sign-mexpt)
1315 (list '%log 'sign-log)
1316 (list 'mabs 'sign-mabs)
1317 (list '$min #'(lambda (x) (sign-minmax (caar x) (cdr x))))
1318 (list '$max #'(lambda (x) (sign-minmax (caar x) (cdr x))))
1319 (list '%csc #'(lambda (x) (sign (inv* (cons (ncons (zl-get (caar x) 'recip)) (cdr x))))))
1320 (list '%csch #'(lambda (x) (sign (inv* (cons (ncons (zl-get (caar x) 'recip)) (cdr x))))))
1321 (list '%acos 'sign-asin/acos/atanh)
1322 (list '%asin 'sign-asin/acos/atanh)
1323 (list '%acosh 'sign-acosh)
1324 (list '%atanh 'sign-asin/acos/atanh)
1326 (list '%signum #'(lambda (x) (sign (cadr x))))
1327 (list '%erf #'(lambda (x) (sign (cadr x))))
1328 (list '$li #'(lambda (x)
1329 (let ((z (first (margs x))) (n (cadadr x)))
1330 (if (and (mnump n) (eq t (mgrp z 0)) (eq t (mgrp 1 z))) (sign z) (sign-any x)))))))
1331 (defun sign (x)
1332 (cond ((mnump x) (setq sign (rgrp x 0) minus nil odds nil evens nil))
1333 ((and *complexsign* (symbolp x) (eq x '$%i))
1334 ;; In Complex Mode the sign of %i is $imaginary.
1335 (setq sign '$imaginary))
1336 ((symbolp x) (if (eq x '$%i) (imag-err x)) (sign-any x))
1337 ((and (consp x) (symbolp (caar x)) (not (specrepp x)) (get (caar x) 'sign-function))
1338 (funcall (get (caar x) 'sign-function) x))
1339 ((and (consp x) (not (specrepp x)) ($subvarp (mop x)) (get (mop (mop x)) 'sign-function))
1340 (funcall (get (mop (mop x)) 'sign-function) x))
1341 ((specrepp x) (sign (specdisrep x)))
1342 ((kindp (caar x) '$posfun) (sign-posfun x))
1343 ((kindp (caar x) '$oddfun) (sign-oddfun x))
1344 (t (sign-any x))))
1346 (defun sign-any (x)
1347 (cond ((and *complexsign*
1348 (symbolp x)
1349 (decl-complexp x))
1350 ;; In Complex Mode look for symbols declared to be complex.
1351 (if ($featurep x '$imaginary)
1352 (setq sign '$imaginary)
1353 (setq sign '$complex)))
1354 ((and *complexsign*
1355 (not (atom x))
1356 (decl-complexp (caar x)))
1357 ;; A function f(x), where f is declared to be imaginary or complex.
1358 (if ($featurep (caar x) '$imaginary)
1359 (setq sign '$imaginary)
1360 (setq sign '$complex)))
1362 (dcompare x 0)
1363 (if (and $assume_pos
1364 (member sign '($pnz $pz $pn) :test #'eq)
1365 (if $assume_pos_pred
1366 (let ((*x* x))
1367 (declare (special *x*))
1368 (is '(($assume_pos_pred) *x*)))
1369 (mapatom x)))
1370 (setq sign '$pos))
1371 (setq minus nil evens nil
1372 odds (if (not (member sign '($pos $neg $zero) :test #'eq))
1373 (ncons x))))))
1375 (defun sign-mtimes (x)
1376 (setq x (cdr x))
1377 (do ((s '$pos) (m) (o) (e)) ((null x) (setq sign s minus m odds o evens e))
1378 (sign1 (car x))
1379 (cond ((eq sign '$zero) (return t))
1380 ((and *complexsign* (eq sign '$complex))
1381 ;; Found a complex factor. We don't return immediately
1382 ;; because another factor could be zero.
1383 (setq s '$complex))
1384 ((and *complexsign* (eq s '$complex))) ; continue the loop
1385 ((and *complexsign* (eq sign '$imaginary))
1386 ;; Found an imaginary factor. Look if we have already one.
1387 (cond ((eq s '$imaginary)
1388 ;; imaginary*imaginary is real. But remember the sign in m.
1389 (setq s (if m '$pos '$neg) m (not m)))
1390 (t (setq s sign))))
1391 ((and *complexsign* (eq s '$imaginary))) ; continue the loop
1392 ((eq sign '$pos))
1393 ((eq sign '$neg) (setq s (flip s) m (not m)))
1394 ((prog2 (setq m (not (eq m minus)) o (nconc odds o) e (nconc evens e))
1395 nil))
1396 ((eq s sign) (when (eq s '$nz) (setq s '$pz)))
1397 ((eq s '$pos) (setq s sign))
1398 ((eq s '$neg) (setq s (flip sign)))
1399 ((or (and (eq s '$pz) (eq sign '$nz))
1400 (and (eq s '$nz) (eq sign '$pz)))
1401 (setq s '$nz))
1402 (t (setq s '$pnz)))
1403 (setq x (cdr x))))
1405 (defun sign-mplus (x &aux s o e m)
1406 (cond ((signdiff x))
1407 ((prog2 (setq s sign e evens o odds m minus) nil))
1408 ((signsum x))
1409 ((prog2 (cond ((strongp s sign))
1410 (t (setq s sign e evens o odds m minus)))
1411 nil))
1412 ((and (not factored) (signfactor x)))
1413 ((strongp sign s))
1414 (t (setq sign s evens e odds o minus m))))
1416 (defun signdiff (x)
1417 (setq sign '$pnz)
1418 (let ((swapped nil) (retval))
1419 (with-compsplt (lhs rhs x)
1420 (if (and (mplusp lhs) (equal rhs 0) (null (cdddr lhs)))
1421 (cond ((and (negp (cadr lhs)) (not (negp (caddr lhs))))
1422 (setq rhs (neg (cadr lhs)) lhs (caddr lhs)))
1423 ;; The following fixes SourceForge bug #3148
1424 ;; where sign(a-b) returned pnz and sign(b-a) returned pos.
1425 ;; Previously, only the case (-a)+b was handled.
1426 ;; Now we also handle a+(-b) by swapping lhs and rhs,
1427 ;; setting a flag "swapped", running through the same code and
1428 ;; then flipping the answer.
1429 ((and (negp (caddr lhs)) (not (negp (cadr lhs))))
1430 (setq rhs (cadr lhs) lhs (neg (caddr lhs)) swapped t))))
1431 (let (dum)
1432 (setq retval
1433 (cond ((or (equal rhs 0) (mplusp lhs)) nil)
1434 ((and (member (constp rhs) '(numer symbol) :test #'eq)
1435 (numberp (setq dum (numer rhs)))
1436 (prog2 (setq rhs dum) nil)))
1437 ((mplusp rhs) nil)
1438 ((and (dcompare lhs rhs) (member sign '($pos $neg $zero) :test #'eq)))
1439 ((and (not (atom lhs)) (not (atom rhs))
1440 (eq (caar lhs) (caar rhs))
1441 (kindp (caar lhs) '$increasing))
1442 (sign (sub (cadr lhs) (cadr rhs)))
1444 ((and (not (atom lhs)) (not (atom rhs))
1445 (eq (caar lhs) (caar rhs))
1446 (kindp (caar lhs) '$decreasing))
1447 (sign (sub (cadr rhs) (cadr lhs)))
1449 ((and (not (atom lhs)) (eq (caar lhs) 'mabs)
1450 (alike1 (cadr lhs) rhs))
1451 (setq sign '$pz minus nil odds nil evens nil) t)
1452 ((signdiff-special lhs rhs))))))
1453 (if swapped
1454 (setq sign (flip sign)))
1455 retval))
1457 (defun signdiff-special (xlhs xrhs)
1458 ;; xlhs may be a constant
1459 (let ((sgn nil))
1460 (when (or (and (realp xrhs) (minusp xrhs)
1461 (not (atom xlhs)) (eq (sign* xlhs) '$pos))
1462 ; e.g. sign(a^3+%pi-1) where a>0
1463 (and (mexptp xlhs)
1464 ;; e.g. sign(%e^x-1) where x>0
1465 (eq (sign* (caddr xlhs)) '$pos)
1466 (or (and
1467 ;; Q^Rpos - S, S<=1, Q>1
1468 (member (sign* (sub 1 xrhs)) '($pos $zero $pz) :test #'eq)
1469 (eq (sign* (sub (cadr xlhs) 1)) '$pos))
1470 (and (not (eq $domain '$complex))
1471 ;; Qpos ^ Rpos - Spos => Qpos - Spos^(1/Rpos)
1472 (eq (sign* (cadr xlhs)) '$pos)
1473 (eq (sign* xrhs) '$pos)
1474 (eq (sign* (sub (cadr xlhs)
1475 (power xrhs (div 1 (caddr xlhs)))))
1476 '$pos))))
1477 (and (mexptp xlhs) (mexptp xrhs)
1478 ;; Q^R - Q^T, Q>1, (R-T) > 0
1479 ;; e.g. sign(2^x-2^y) where x>y
1480 (alike1 (cadr xlhs) (cadr xrhs))
1481 (eq (sign* (sub (cadr xlhs) 1)) '$pos)
1482 (eq (sign* (sub (caddr xlhs) (caddr xrhs))) '$pos)))
1483 (setq sgn '$pos))
1485 ;; sign(sin(x)+c)
1486 (when (and (not (atom xlhs))
1487 (member (caar xlhs) '(%sin %cos))
1488 (zerop1 ($imagpart (cadr xlhs))))
1489 (cond ((eq (sign* (add xrhs 1)) '$neg) ;; c > 1
1490 (setq sgn '$pos))
1491 ((eq (sign* (add xrhs -1)) '$pos) ;; c < -1
1492 (setq sgn '$neg))
1493 ((zerop1 (add xrhs 1)) ;; c = 1
1494 (setq sgn '$pz))
1495 ((zerop1 (add xrhs -1)) ;; c = -1
1496 (setq sgn '$nz))))
1498 ;; sign(signum(x)+c)
1499 (when (and (not (atom xlhs))
1500 (eq (caar xlhs) '%signum)
1501 (zerop1 ($imagpart (cadr xlhs))))
1502 (cond ((eq (sign* (add xrhs 1)) '$neg) ;; c > 1
1503 (setq sgn '$pos))
1504 ((eq (sign* (add xrhs -1)) '$pos) ;; c < -1
1505 (setq sgn '$neg))
1506 ((zerop1 (add xrhs 1)) ;; c = 1
1507 (setq sgn '$pz))
1508 ((zerop1 (add xrhs -1)) ;; c = -1
1509 (setq sgn '$nz))
1510 ((zerop1 xrhs) ;; c = 0 (necessary?)
1511 (setq sgn '$pnz))
1512 (t ;; -1 < c < 1, but c # 0
1513 (setq sgn '$pn))))
1515 (when (and $useminmax (or (minmaxp xlhs) (minmaxp xrhs)))
1516 (setq sgn (signdiff-minmax xlhs xrhs)))
1517 (when sgn (setq sign sgn minus nil odds nil evens nil)
1518 t)))
1520 ;;; Look for symbols with an assumption a > n or a < -n, where n is a number.
1521 ;;; For this case shift the symbol a -> a+n in a summation and multiplication.
1522 ;;; This handles cases like a>1 and b>1 gives sign(a+b-2) -> pos.
1524 (defun sign-shift (expr)
1525 (do ((l (cdr expr) (cdr l))
1526 (fl nil)
1527 (acc nil))
1528 ((null l) (addn acc nil))
1529 (cond ((symbolp (car l))
1530 ;; Get the facts related to the symbol (car l)
1531 ;; Reverse the order to test the newest facts first.
1532 (setq fl (reverse (cdr (facts1 (car l)))))
1533 (push (car l) acc)
1534 (dolist (f fl)
1535 (cond ((and (eq (caar f) 'mgreaterp)
1536 (mnump (caddr f))
1537 (eq ($sign (caddr f)) '$pos))
1538 ;; The case a > n, where a is a symbol and n a number.
1539 ;; Add the number to the list of terms.
1540 (return (push (caddr f) acc)))
1541 ((and (eq (caar f) 'mgreaterp)
1542 (mnump (cadr f))
1543 (eq ($sign (cadr f)) '$neg))
1544 ;; The case a < -n, where a is a symbol and n a number.
1545 ;; Add the number to the list of terms.
1546 (return (push (cadr f) acc))))))
1547 ((mtimesp (car l))
1548 (let ((acctimes) (flag))
1549 ;; Go through the factors of the multiplication.
1550 (dolist (ll (cdar l))
1551 (cond ((symbolp ll)
1552 ;; Get the facts related to the symbol (car l)
1553 ;; Reverse the order to test the newest facts first.
1554 (setq fl (reverse (cdr (facts1 ll))))
1555 (dolist (f fl)
1556 (cond ((and (eq (caar f) 'mgreaterp)
1557 (mnump (caddr f))
1558 (eq ($sign (caddr f)) '$pos))
1559 ;; The case a > n, where a is a symbol and n a
1560 ;; number. Add the number to the list of terms.
1561 (setq flag t)
1562 (return (push (add ll (caddr f)) acctimes)))
1563 ((and (eq (caar f) 'mgreaterp)
1564 (mnump (cadr f))
1565 (eq ($sign (cadr f)) '$neg))
1566 ;; The case a < -n, where a is a symbol and n a
1567 ;; number. Add the number to the list of terms.
1568 (setq flag t)
1569 (return (push (add ll (cadr f)) acctimes)))))
1570 (when (not flag) (push ll acctimes)))
1572 (push ll acctimes))))
1573 (if flag
1574 ;; If a shift has been done expand the factors.
1575 (push ($multthru (muln acctimes nil)) acc)
1576 (push (muln acctimes nil) acc))))
1578 (push (car l) acc)))))
1580 (defun signsum (x)
1581 (setq x (sign-shift x))
1582 ;; x might be simplified to an atom in sign-shift.
1583 (when (atom x) (setq x (cons '(mplus) (list x))))
1584 (do ((l (cdr x) (cdr l)) (s '$zero))
1585 ((null l) (setq sign s minus nil odds (list x) evens nil)
1586 (cond (*complexsign*
1587 ;; Because we have continued the loop in Complex Mode
1588 ;; we have to look for the sign '$pnz and return nil.
1589 (if (eq s '$pnz) nil t))
1590 (t t))) ; in Real Mode return T
1591 ;; Call sign1 and not sign, because sign1 handles constant expressions.
1592 (sign1 (car l))
1593 (cond ((and *complexsign*
1594 (or (eq sign '$complex) (eq sign '$imaginary)))
1595 ;; Found a complex or imaginary expression. The sign is $complex.
1596 (setq sign '$complex odds nil evens nil minus nil)
1597 (return t))
1598 ((or (and (eq sign '$zero)
1599 (setq x (sub x (car l))))
1600 (and (eq s sign) (not (eq s '$pn))) ; $PN + $PN = $PNZ
1601 (and (eq s '$pos) (eq sign '$pz))
1602 (and (eq s '$neg) (eq sign '$nz))))
1603 ((or (and (member sign '($pz $pos) :test #'eq) (member s '($zero $pz) :test #'eq))
1604 (and (member sign '($nz $neg) :test #'eq) (member s '($zero $nz) :test #'eq))
1605 (and (eq sign '$pn) (eq s '$zero)))
1606 (setq s sign))
1608 (cond (*complexsign*
1609 ;; In Complex Mode we have to continue the loop to look further
1610 ;; for a complex or imaginay expression.
1611 (setq s '$pnz))
1613 ;; In Real mode the loop stops when the sign is 'pnz.
1614 (setq sign '$pnz odds (list x) evens nil minus nil)
1615 (return nil)))))))
1617 (defun signfactor (x)
1618 (let (y (factored t))
1619 (setq y (factor-if-small x))
1620 (cond ((or (mplusp y) (> (conssize y) 50.))
1621 (setq sign '$pnz)
1622 nil)
1623 (t (sign y)))))
1625 (defun factor-if-small (x)
1626 (if (< (conssize x) 51.)
1627 (let ($ratprint)
1628 (declare (special $ratprint))
1629 (factor x)) x))
1631 (defun sign-mexpt (x)
1632 (let* ((expt (caddr x)) (base1 (cadr x))
1633 (sign-expt (sign1 expt)) (sign-base (sign1 base1))
1634 (evod (evod expt)))
1635 ;; The variable sign is now equal to sign-base. This is used below
1636 ;; in some places to avoid an assignment operation for sign.
1637 (cond ((and (eq sign-base '$zero)
1638 (member sign-expt '($zero $neg $nz) :test #'eq))
1639 (dbzs-err x))
1640 ((eq sign-expt '$zero) (setq sign '$pos))
1641 ((eq sign-base '$zero))
1643 ((and *complexsign* (or (eq sign-expt '$complex)
1644 (eq sign-expt '$imaginary)
1645 (eq sign-base '$complex)))
1646 ;; Base is complex or exponent is complex or imaginary.
1647 ;; The sign is $complex.
1648 (when *debug-compar*
1649 (format t "~&in SIGN-MEXPT for ~A, sign is complex.~%" x))
1650 (setq sign '$complex))
1652 ((and *complexsign*
1653 (eq sign-base '$neg)
1654 (eq (evod ($expand (mul 2 expt))) '$odd))
1655 ;; Base is negative and the double of the exponent is odd.
1656 ;; Result is imaginary.
1657 (when *debug-compar*
1658 (format t "~&in SIGN-MEXPT for ~A, sign is $imaginary.~%" x))
1659 (setq sign '$imaginary))
1661 ((and *complexsign*
1662 (eq sign-base '$imaginary))
1663 ;; An imaginary base. Look for even or odd exponent.
1664 (when *debug-compar*
1665 (format t "~&in SIGN-MEXPT for ~A, base is $imaginary.~%" x))
1666 (cond
1667 ((and (integerp expt) (eq evod '$even))
1668 (setq sign (if (eql (mod expt 4) 0) '$pz '$nz)))
1669 ((and (integerp expt) (eq evod '$odd))
1670 (setq sign '$imaginary
1671 minus (if (eql (mod (- expt 1) 4) 0) t nil)))
1672 (t (setq sign '$complex))))
1674 ((and *complexsign*
1675 (not evod)
1676 (not (ratnump expt))
1677 (not (member sign-base '($pos $pz $zero) :test #'eq)))
1678 (when *debug-compar*
1679 (format t "~&in SIGN-MEXPT for ~A, base is not $pos, $pz or $zero.~%" x))
1680 (setq sign (if (maxima-integerp expt) '$pnz '$complex)))
1681 ((eq sign-base '$pos))
1682 ((eq evod '$even)
1683 (cond ((eq sign-expt '$neg)
1684 (setq sign '$pos minus nil evens (ncons base1) odds nil))
1685 ((member sign-base '($pn $neg) :test #'eq)
1686 (setq sign '$pos minus nil
1687 evens (nconc odds evens)
1688 odds nil))
1689 (t (setq sign '$pz minus nil
1690 evens (nconc odds evens)
1691 odds nil))))
1692 ((and (member sign-expt '($neg $nz) :test #'eq)
1693 (member sign-base '($nz $pz $pnz) :test #'eq))
1694 (setq sign (if (eq sign-base '$pz)
1695 '$pos
1696 '$pn)))
1697 ((member sign-expt '($pz $nz $pnz) :test #'eq)
1698 (cond ((member sign-base '($neg $nz) :test #'eq)
1699 (setq odds (ncons x) sign (if (eq sign-base '$neg) '$pn '$pnz)))))
1700 ((eq sign-expt '$pn)
1701 (cond ((eq sign-base '$neg)
1702 (setq sign '$pn))
1703 ((eq sign-base '$nz)
1704 (setq sign '$pnz))))
1705 ((ratnump expt)
1706 (cond ((mevenp (cadr expt))
1707 (cond ((member sign-base '($pn $neg) :test #'eq)
1708 (setq sign-base '$pos))
1709 ((member sign-base '($pnz $nz) :test #'eq)
1710 (setq sign-base '$pz)))
1711 (setq evens (nconc odds evens)
1712 odds nil minus nil))
1713 ((mevenp (caddr expt))
1714 (cond (*complexsign*
1715 (when (not (member sign-base '($pos $pz)))
1716 ;; In Complex Mode the sign is $complex.
1717 (setq sign-base (setq sign-expt '$complex))))
1718 (complexsign
1719 ;; The only place the variable complexsign
1720 ;; is used. Unfortunately, one routine in
1721 ;; to_poly.lisp in /share/to_poly_solve depends on
1722 ;; this piece of code. Perhaps we can remove
1723 ;; the dependency. (12/2008)
1724 (setq sign-base (setq sign-expt '$pnz)))
1725 ((eq sign-base '$neg) (imag-err x))
1726 ((eq sign-base '$pn)
1727 (setq sign-base '$pos))
1728 ((eq sign-base '$nz)
1729 (setq sign-base '$zero))
1730 (t (setq sign-base '$pz)))))
1731 (cond ((eq sign-expt '$neg)
1732 (cond ((eq sign-base '$zero) (dbzs-err x))
1733 ((eq sign-base '$pz)
1734 (setq sign-base '$pos))
1735 ((eq sign-base '$nz)
1736 (setq sign-base '$neg))
1737 ((eq sign-base '$pnz)
1738 (setq sign-base '$pn)))))
1739 (setq sign sign-base))
1740 ((member sign-base '($neg $nz) :test #'eq)
1741 (cond ((eq evod '$odd))
1742 (*complexsign*
1743 (setq sign '$complex))
1744 ((eq sign-base '$neg)
1745 (setq sign '$pn))
1746 (t (setq sign '$pnz)))))))
1748 ;;; Determine the sign of log(expr). This function changes the special variable sign.
1750 (defun sign-log (x)
1751 (setq x (cadr x))
1752 (setq sign
1753 (cond ((eq t (mgrp x 0))
1754 (cond ((eq t (mgrp 1 x)) '$neg)
1755 ((eq t (meqp x 1)) '$zero);; log(1) = 0.
1756 ((eq t (mgqp 1 x)) '$nz)
1757 ((eq t (mgrp x 1)) '$pos)
1758 ((eq t (mgqp x 1)) '$pz)
1759 ((eq t (mnqp x 1)) '$pn)
1760 (t '$pnz)))
1761 ((and *complexsign* (eql 1 (cabs x))) '$imaginary)
1762 (*complexsign* '$complex)
1763 (t '$pnz))))
1765 (defun sign-mabs (x)
1766 (let ((*complexsign* t))
1767 (sign (cadr x))
1768 (cond ((member sign '($pos $zero) :test #'eq))
1769 ((member sign '($neg $pn) :test #'eq) (setq sign '$pos))
1770 (t (setq sign '$pz minus nil evens (nconc odds evens) odds nil)))))
1772 (defun sign-asin/acos/atanh (x)
1773 (cond ((and *complexsign*
1774 (or (not (eq t (mgqp (cadr x) -1)))
1775 (not (eq t (mgqp 1 (cadr x)))))) ; x < -1 or x > 1
1776 (setq sign '$complex))
1777 ((not (eq (caar x) '%acos))
1778 (sign-oddfun x))
1779 ((eq t (mlsp (cadr x) 1)) ; x < 1
1780 (sign-posfun x))
1781 (t ; x <= 1
1782 (setq sign '$pz))))
1784 (defun sign-acosh (x)
1785 (cond ((and *complexsign* (not (eq t (mgqp (cadr x) 1)))) ; x < 1
1786 (setq sign '$complex))
1787 ((eq t (mgrp (cadr x) 1)) ; x > 1
1788 (sign-posfun x))
1789 (t ; x >= 1
1790 (setq sign '$pz))))
1792 ;; This code handles sign(sin(x)), where -%pi <= x <= %pi. Of course, at the
1793 ;; expense of a great deal of additional complexity, this code could catch far
1794 ;; more cases.
1796 ;; Like all sign functions, it works by setting the value of the special variable
1797 ;; sign.
1799 ;; The sign of sin(constant expression) is supposed to be handled elsewhere,
1800 ;; so we don't consider that case here. But running the testsuite,
1801 ;; sometimes sign-sin receives a constant expression. This code could
1802 ;; handle this case, but I think that's fixing a bug in the wrong place.
1804 ;; The general simplifier converts sin(%i*a) to %i*sinh(a), so I don't
1805 ;; think this code needs to handle the case of sin(x), where x is purely imaginary.
1807 ;; This code could handle sin({zerob, zeroa, ind}), but it doesn't.
1809 ;; Finally, sometimes the argument to sin is not simplified; for example,
1810 ;; running the testsuite sometimes this code recieves the expression sin(1/(1/x)).
1811 ;; We could simplify it, but I think that it would be better to fix the code
1812 ;; so that the input to sign-sin is always simplified.
1813 (defun sign-sin (e) ; e = sin(x)
1814 (let ((x (cadr e)) (y 0))
1815 ;; When *complexsign* is true, find the rectangular form of
1816 ;; the argument to sin.
1817 (when *complexsign*
1818 (setq x (risplit x))
1819 (setq y (cdr x)
1820 x (car x)))
1821 (cond
1822 ;; When y = 0 and -%pi <= x <= %pi, sign(sin(x)) = sign(x)
1823 ((and (eql y 0)
1824 (eq t (mgqp x (mul -1 '$%pi)))
1825 (eq t (mgqp '$%pi x)))
1826 (sign x))
1827 ;; When *complexsign* is true & y # 0, set sign to complex.
1828 ;; To test y # 0, we'll use (not (eql y 0)))
1829 ((and *complexsign* (not (eql y 0)))
1830 (setf sign '$complex))
1831 (t (setf sign '$pnz))))
1832 nil)
1833 (putprop '%sin 'sign-sin 'sign-function)
1835 (defun sign-cos (e) ; e = cos(x)
1836 (let ((x (cadr e)) (y 0))
1837 ;; When *complexsign* is true, find the rectangular form of
1838 ;; the argument to cos.
1839 (when *complexsign*
1840 (setq x (risplit x))
1841 (setq y (cdr x)
1842 x (car x)))
1843 (cond
1844 ;; When y = 0 and -%pi/2 <= x <= 3 %pi/2, sign(cos(x)) = sign(%pi/2-x)
1845 ((and (eql y 0)
1846 (eq t (mgqp x (div '$%pi -2)))
1847 (eq t (mgqp (div (mul 3 '$%pi) 2) x)))
1848 (sign (sub (div '$%pi 2) x)))
1849 ;; When *complexsign* is true & y # 0, set sign to complex.
1850 ;; To test y # 0, we'll use (not (eql y 0)))
1851 ((and *complexsign* (not (eql y 0)))
1852 (setf sign '$complex))
1853 (t (setf sign '$pnz))))
1854 nil)
1855 (putprop '%cos 'sign-cos 'sign-function)
1857 ;; When Maxima can deduce that when -1<x<=0, then ceiling(x) simplifies to 0. So this
1858 ;; case is simplified away before it gets to this function.
1859 (defun sign-ceiling (x)
1860 (setq x (cadr x)) ; ceiling(x) --> x
1861 (setf sign (cond ((eq t (mgrp x 0)) '$pos) ; ceiling(> 0) = pos
1862 ((eq t (mgrp x -1)) '$pz) ; ceiling(> -1) = pz
1863 ((eq t (mgqp -1 x)) '$neg) ; ceiling(<= -1) = neg
1864 ((eq t (mgrp 0 x)) '$nz) ; ceiling( < 0) = nz
1865 (t '$pnz))))
1866 (putprop '$ceiling 'sign-ceiling 'sign-function)
1868 (defun sign-floor (x)
1869 (setq x (cadr x)) ; floor(x) --> x
1870 (setf sign (cond ((eq t (mgqp x 1)) '$pos) ;floor(>= 1) = pos
1871 ((eq t (mgqp x 0)) '$pz) ;floor(>= 0) = pz
1872 ((eq t (mgrp 0 x)) '$neg) ;floor(< 0) = neg
1873 ((eq t (mgqp 0 x)) '$nz) ;floor(<= 0) = nz
1874 (t '$pnz))))
1875 (putprop '$floor 'sign-floor 'sign-function)
1877 ;;; Compare min/max
1879 ;;; Macros used in simp min/max
1880 ;;; If op is min, use body; if not, negate sign constants in body
1881 ;;; Used to avoid writing min and max code separately: just write the min code
1882 ;;; in such a way that its dual works for max
1883 (defmacro minmaxforms (op &rest body)
1884 `(if (eq ,op '$min)
1885 ,@body
1886 ,@(sublis '(($neg . $pos)
1887 ($nz . $pz)
1888 ($pz . $nz)
1889 ($pos . $neg)
1890 ;;($zero . $zero)
1891 ;;($pn . $pn)
1892 ;;($pnz . $pnz)
1894 ($max . $min)
1895 ($min . $max)
1897 ($inf . $minf)
1898 ($minf . $inf))
1899 body)))
1901 (defun sign-minmax (op args)
1902 (do ((sgn (minmaxforms op '$pos) ;identity element for min
1903 (sminmax op sgn (sign* (car l))))
1904 (end (minmaxforms op '$neg))
1905 (l args (cdr l)))
1906 ((or (null l) (eq sgn end))
1907 (setq minus nil
1908 odds (if (not (member sgn '($pos $neg $zero) :test #'eq))
1909 (ncons (cons (list op) args)))
1910 evens nil
1911 sign sgn))))
1913 ;; sign(op(a,b)) = sminmax(sign(a),sign(b))
1914 ;; op is $min/$max; s1/s2 in neg, nz, zero, pz, pos, pn, pnz
1915 (defun sminmax (op s1 s2)
1916 (minmaxforms
1918 ;; Many of these cases don't come up in simplified expressions,
1919 ;; since e.g. sign(a)=neg and sign(b)=pos implies min(a,b)=a
1920 ;; the order of these clauses is important
1921 (cond ((eq s1 '$pos) s2)
1922 ((eq s2 '$pos) s1)
1923 ((eq s1 s2) s1)
1924 ((or (eq s1 '$neg) (eq s2 '$neg)) '$neg)
1925 ((or (eq s1 '$nz) (eq s2 '$nz)) '$nz)
1926 ((eq s1 '$zero) (if (eq s2 '$pz) '$zero '$nz))
1927 ((eq s2 '$zero) (if (eq s1 '$pz) '$zero '$nz))
1928 (t '$pnz))))
1930 (defun minmaxp (ex)
1931 (cond ((atom ex) nil)
1932 ((member (caar ex) '($min $max) :test #'eq) (caar ex))
1933 (t nil)))
1935 (defun signdiff-minmax (l r)
1936 ;; sign of l-r; nil if unknown (not PNZ)
1937 (let* ((lm (minmaxp l))
1938 (rm (minmaxp r))
1939 (ll (if lm (cdr l)))
1940 (rr (if rm (cdr r)))) ;distinguish between < and <= argument lists of min/max
1941 (minmaxforms
1942 (or rm lm)
1943 (cond ((eq lm rm) ; min(a,...) - min(b,...)
1944 (multiple-value-bind (both onlyl onlyr) (intersect-info ll rr)
1945 (declare (ignore both))
1946 (cond ((null onlyl) '$pz) ; min(a,b) - min(a,b,c)
1947 ((null onlyr) '$nz) ; min(a,b,c) - min(a,b)
1948 ;; TBD: add processing for full onlyl/onlyr case
1949 (t nil))))
1950 ;; TBD: memalike and set-disjointp are crude approx.
1951 ((null lm) (if (memalike l rr) '$pz)) ; a - min(a,b)
1952 ((null rm) (if (memalike r ll) '$nz)) ; min(a,b) - a
1953 (t ; min/max or max/min
1954 (if (not (set-disjointp ll rr)) '$pz)))))) ; max(a,q,r) - min(a,s,t)
1956 (defun intersect-info (a b)
1957 (let ((both nil)
1958 (onlya nil)
1959 (onlyb nil))
1960 (do-merge-asym
1962 #'like
1963 #'$orderlessp
1964 #'(lambda (x) (push x both))
1965 #'(lambda (x) (push x onlya))
1966 #'(lambda (x) (push x onlyb)))
1967 (values
1968 (reverse both)
1969 (reverse onlya)
1970 (reverse onlyb))))
1972 ;;; end compare min/max
1974 (defun sign-posfun (xx)
1975 (declare (ignore xx))
1976 (setq sign '$pos
1977 minus nil
1978 odds nil
1979 evens nil))
1981 (defun sign-oddfun (x)
1982 (cond ((kindp (caar x) '$increasing)
1983 ; Take the sign of the argument
1984 (sign (cadr x)))
1985 ((kindp (caar x) '$decreasing)
1986 ; Take the sign of negative of the argument
1987 (sign (neg (cadr x))))
1989 ; If the sign of the argument is zero, then we're done (the sign of
1990 ; the function value is the same). Otherwise, punt to SIGN-ANY.
1991 (sign (cadr x))
1992 (unless (eq sign '$zero)
1993 (sign-any x)))))
1995 (defun imag-err (x)
1996 (if sign-imag-errp
1997 (merror (intl:gettext "sign: argument cannot be imaginary; found ~M") x)
1998 (throw 'sign-imag-err t)))
2000 (defun dbzs-err (x)
2001 (merror (intl:gettext "sign: division by zero in ~M") x))
2003 ;; Return true iff e is an expression with operator op1, op2,...,or opn.
2005 (defun op-equalp (e &rest op)
2006 (and (consp e) (consp (car e)) (some #'(lambda (s) (equal (caar e) s)) op)))
2008 ;; Return true iff the operator of a is a Maxima relation operator.
2010 (defun mrelationp (a)
2011 (op-equalp a 'mlessp 'mleqp 'mequal 'mnotequal 'mgeqp 'mgreaterp))
2013 ;; This version of featurep applies ratdisrep to the first argument. This
2014 ;; change allows things like featurep(rat(n),integer) --> true when n has
2015 ;; been declared an integer.
2017 (defmfun $featurep (e ind)
2018 (setq e ($ratdisrep e))
2019 (cond ((not (symbolp ind))
2020 (merror
2021 (intl:gettext "featurep: second argument must be a symbol; found ~M")
2022 ind))
2023 ;; Properties not related to the assume database.
2024 ((and (member ind opers) (safe-get e ind)))
2025 ((and (member ind '($evfun $evflag $bindtest $nonarray))
2026 (safe-get e (stripdollar ind))))
2027 ((and (eq ind '$noun)
2028 (safe-get e (stripdollar ind))
2030 ((and (member ind '($scalar $nonscalar $mainvar))
2031 (mget e ind)))
2032 ((and (eq ind '$feature)
2033 (member e $features)
2035 ((eq ind '$alphabetic)
2036 (dolist (l (coerce e 'list) t)
2037 (when (not (member l *alphabet*)) (return nil))))
2038 ;; Properties related to the assume database.
2039 ((eq ind '$integer) (maxima-integerp e))
2040 ((eq ind '$noninteger) (nonintegerp e))
2041 ((eq ind '$even) (mevenp e))
2042 ((eq ind '$odd) (moddp e))
2043 ((eq ind '$real)
2044 (if (atom e)
2045 (or (numberp e) (kindp e '$real) (numberp (numer e)))
2046 (free ($rectform e) '$%i)))
2047 ((symbolp e) (kindp e ind))))
2049 ;; Give a function the maps-integers-to-integers property when it is integer
2050 ;; valued on the integers; give it the integer-valued property when its
2051 ;; range is a subset of the integers. What have I missed?
2053 (setf (get 'mplus 'maps-integers-to-integers) t)
2054 (setf (get 'mtimes 'maps-integers-to-integers) t)
2055 (setf (get 'mabs 'maps-integers-to-integers) t)
2056 (setf (get '$max 'maps-integers-to-integers) t)
2057 (setf (get '$min 'maps-integers-to-integers) t)
2059 (setf (get '$floor 'integer-valued) t)
2060 (setf (get '$ceiling 'integer-valued) t)
2061 (setf (get '$charfun 'integer-valued) t)
2063 (defun maxima-integerp (x)
2064 (cond ((integerp x))
2065 ((mnump x) nil)
2066 ((and (symbolp x)
2067 (or (kindp x '$integer)
2068 (kindp x '$even)
2069 (kindp x '$odd)
2070 (check-integer-facts x))))
2071 (t (let ((x-op (and (consp x) (consp (car x)) (caar x))) ($prederror nil))
2072 (cond ((null x-op) nil)
2073 ((not (symbolp x-op)) nil) ; fix for mqapply at some point?
2074 ((eq x-op 'mrat) (and (integerp (cadr x)) (equal (cddr x) 1)))
2075 ;; mtimes and mplus are generally handled by this clause
2076 ((and (get x-op 'maps-integers-to-integers) (every #'maxima-integerp (margs x))))
2077 ;; Special case for 1/2*...*even
2078 ((eq x-op 'mtimes)
2079 (and (mnump (cadr x))
2080 (integerp (mul 2 (cadr x)))
2081 (every 'maxima-integerp (cddr x))
2082 (some #'(lambda (s) ($featurep s '$even)) (rest (margs x)))))
2083 ((eq x-op 'mexpt)
2084 (and (every #'maxima-integerp (margs x))
2085 (null (mevalp (mlsp (caddr x) 0)))))
2086 ;; ! in Maxima allows real arguments
2087 ((eq x-op 'mfactorial)
2088 (and (maxima-integerp (cadr x))
2089 (not (mevalp (mlsp (cadr x) 0)))))
2090 ((eq x-op '%gamma)
2091 (and (maxima-integerp (cadr x))
2092 (not (mevalp (mlsp (cadr x) 1)))))
2093 ;; other x-ops
2094 ((or ($featurep ($verbify x-op) '$integervalued)
2095 (get x-op 'integer-valued))))))))
2097 ;; When called with mode 'integer look into the database for symbols which are
2098 ;; declared to be equal to an integer or an expression which is an integer.
2099 ;; In mode 'evod look for odd and even expressions.
2100 (defun check-integer-facts (x &optional (mode 'integer))
2101 (do ((factsl (cdr (facts1 x)) (cdr factsl))
2102 fact)
2103 ((null factsl) nil)
2104 (setq fact (car factsl))
2105 (cond ((and (not (atom fact))
2106 (eq (caar fact) '$equal))
2107 (cond ((and (symbolp (cadr fact))
2108 (eq (cadr fact) x))
2109 ;; Case equal(x,expr): Test expr to be an integer.
2110 (cond ((symbolp (caddr fact))
2111 (cond ((and (eq mode 'integer)
2112 (or (kindp (caddr fact) '$integer)
2113 (kindp (caddr fact) '$odd)
2114 (kindp (caddr fact) '$even)))
2115 (return t))
2116 ((eq mode 'evod)
2117 (cond ((kindp (caddr fact) '$odd)
2118 (return '$odd))
2119 ((kindp (caddr fact) '$even)
2120 (return '$even))
2121 (t (return nil))))
2122 (t (return nil))))
2124 (cond ((eq mode 'integer)
2125 (return (maxima-integerp (caddr fact))))
2126 ((eq mode 'evod)
2127 (return (evod (caddr fact))))
2128 (t (return nil))))))
2129 ((and (symbolp (caddr fact))
2130 (eq (caddr fact) x))
2131 ;; Case equal(expr,x): Test expr to be an integer.
2132 (cond ((symbolp (caddr fact))
2133 (cond ((and (eq mode 'integer)
2134 (or (kindp (cadr fact) '$integer)
2135 (kindp (cadr fact) '$odd)
2136 (kindp (cadr fact) '$even)))
2137 (return t))
2138 ((eq mode 'evod)
2139 (cond ((kindp (cadr fact) '$odd)
2140 (return '$odd))
2141 ((kindp (cadr fact) '$even)
2142 (return '$even))
2143 (t (return nil))))
2144 (t (return nil))))
2146 (cond ((eq mode 'integer)
2147 (return (maxima-integerp (cadr fact))))
2148 ((eq mode 'evod)
2149 (return (evod (cadr fact))))
2150 (t (return nil)))))))))))
2152 (defun nonintegerp (e)
2153 (cond ((and (symbolp e) (or (kindp e '$noninteger) (check-noninteger-facts e) (kindp e '$irrational)))) ;declared noninteger
2154 ((mnump e)
2155 (if (integerp e) nil t)) ;all floats are noninteger and integers are not nonintegers
2156 (($ratp e)
2157 (nonintegerp ($ratdisrep e)))
2158 (t (eq t (mgrp e (take '($floor) e))))))
2161 ;; Look into the database for symbols which are declared to be equal
2162 ;; to a noninteger or an expression which is a noninteger.
2163 (defun check-noninteger-facts (x)
2164 (do ((factsl (cdr (facts1 x)) (cdr factsl)))
2165 ((null factsl) nil)
2166 (cond ((and (not (atom (car factsl)))
2167 (eq (caar (car factsl)) '$equal))
2168 (cond ((and (symbolp (cadr (car factsl)))
2169 (eq (cadr (car factsl)) x))
2170 ;; Case equal(x,expr): Test expr to be a noninteger.
2171 (cond ((symbolp (caddr (car factsl)))
2172 (if (kindp (caddr (car factsl)) '$noninteger)
2173 (return t)))
2175 (return (nonintegerp (caddr (car factsl)))))))
2176 ((and (symbolp (caddr (car factsl)))
2177 (eq (caddr (car factsl)) x))
2178 ;; Case equal(expr,x): Test expr to be a noninteger.
2179 (cond ((symbolp (cadr (car factsl)))
2180 (if (kindp (cadr (car factsl)) '$noninteger)
2181 (return t)))
2183 (return (nonintegerp (cadr (car factsl))))))))))))
2185 (defun intp (l)
2186 (every #'maxima-integerp (cdr l)))
2188 (defun mevenp (e)
2189 (cond ((integerp e) (not (oddp e)))
2190 ((mnump e) nil)
2191 (t (eq '$even (evod e)))))
2193 (defun moddp (e)
2194 (cond ((integerp e) (oddp e))
2195 ((mnump e) nil)
2196 (t (eq '$odd (evod e)))))
2198 ;; An extended evod that recognizes that abs(even) is even and
2199 ;; abs(odd) is odd.
2201 (defun evod (e)
2202 (cond ((integerp e) (if (oddp e) '$odd '$even))
2203 ((mnump e) nil)
2204 ((atom e)
2205 (cond ((kindp e '$odd) '$odd)
2206 ((kindp e '$even) '$even)
2207 ;; Check the database for facts.
2208 ((symbolp e) (check-integer-facts e 'evod))))
2209 ((eq 'mtimes (caar e)) (evod-mtimes e))
2210 ((eq 'mplus (caar e)) (evod-mplus e))
2211 ((eq 'mabs (caar e)) (evod (cadr e))) ;; extra code
2212 ((eq 'mexpt (caar e)) (evod-mexpt e))))
2214 (defun evod-mtimes (x)
2215 (do ((l (cdr x) (cdr l)) (flag '$odd))
2216 ((null l) flag)
2217 (setq x (evod (car l)))
2218 (cond ((eq '$odd x))
2219 ((eq '$even x) (setq flag '$even))
2220 ((maxima-integerp (car l)) (cond ((eq '$odd flag) (setq flag nil))))
2221 (t (return nil)))))
2223 (defun evod-mplus (x)
2224 (do ((l (cdr x) (cdr l)) (flag))
2225 ((null l) (cond (flag '$odd) (t '$even)))
2226 (setq x (evod (car l)))
2227 (cond ((eq '$odd x) (setq flag (not flag)))
2228 ((eq '$even x))
2229 (t (return nil)))))
2231 (defun evod-mexpt (x)
2232 (when (and (integerp (caddr x)) (not (minusp (caddr x))))
2233 (evod (cadr x))))
2235 (declare-top (special mgqp mlqp))
2237 (defmode cl ()
2238 (atom (selector +labs) (selector -labs) (selector data)))
2240 (defmacro c-dobj (&rest x)
2241 `(list ,@x))
2243 (defun dcompare (x y)
2244 (setq odds (list (sub x y)) evens nil minus nil
2245 sign (cond ((eq x y) '$zero)
2246 ((or (eq '$inf x) (eq '$minf y)) '$pos)
2247 ((or (eq '$minf x) (eq '$inf y)) '$neg)
2248 (t (dcomp x y)))))
2250 (defun dcomp (x y)
2251 (let (mgqp mlqp)
2252 (setq x (dinternp x) y (dinternp y))
2253 (cond ((or (null x) (null y)) '$pnz)
2254 ((progn (clear) (deq x y) (sel y +labs)))
2255 (t '$pnz))))
2257 (defun deq (x y)
2258 (cond ((dmark x '$zero) nil)
2259 ((eq x y))
2260 (t (do ((l (sel x data) (cdr l))) ((null l))
2261 (if (and (visiblep (car l)) (deqf x y (car l))) (return t))))))
2263 (defun deqf (x y f)
2264 (cond ((eq 'meqp (caar f))
2265 (if (eq x (cadar f)) (deq (caddar f) y) (deq (cadar f) y)))
2266 ((eq 'mgrp (caar f))
2267 (if (eq x (cadar f)) (dgr (caddar f) y) (dls (cadar f) y)))
2268 ((eq 'mgqp (caar f))
2269 (if (eq x (cadar f)) (dgq (caddar f) y) (dlq (cadar f) y)))
2270 ((eq 'mnqp (caar f))
2271 (if (eq x (cadar f)) (dnq (caddar f) y) (dnq (cadar f) y)))))
2273 (defun dgr (x y)
2274 (cond ((dmark x '$pos) nil)
2275 ((eq x y))
2276 (t (do ((l (sel x data) (cdr l)))
2277 ((null l))
2278 (when (or mlqp (and (visiblep (car l)) (dgrf x y (car l))))
2279 (return t))))))
2281 (defun dgrf (x y f)
2282 (cond ((eq 'mgrp (caar f)) (if (eq x (cadar f)) (dgr (caddar f) y)))
2283 ((eq 'mgqp (caar f)) (if (eq x (cadar f)) (dgr (caddar f) y)))
2284 ((eq 'meqp (caar f))
2285 (if (eq x (cadar f))
2286 (dgr (caddar f) y)
2287 (dgr (cadar f) y)))))
2289 (defun dls (x y)
2290 (cond ((dmark x '$neg) nil)
2291 ((eq x y))
2292 (t (do ((l (sel x data) (cdr l)))
2293 ((null l))
2294 (when (or mgqp (and (visiblep (car l)) (dlsf x y (car l))))
2295 (return t))))))
2297 (defun dlsf (x y f)
2298 (cond ((eq 'mgrp (caar f)) (if (eq x (caddar f)) (dls (cadar f) y)))
2299 ((eq 'mgqp (caar f)) (if (eq x (caddar f)) (dls (cadar f) y)))
2300 ((eq 'meqp (caar f))
2301 (if (eq x (cadar f)) (dls (caddar f) y) (dls (cadar f) y)))))
2303 (defun dgq (x y)
2304 (cond ((member (sel x +labs) '($pos $zero) :test #'eq) nil)
2305 ((eq '$nz (sel x +labs)) (deq x y))
2306 ((eq '$pn (sel x +labs)) (dgr x y))
2307 ((dmark x '$pz) nil)
2308 ((eq x y) (setq mgqp t) nil)
2309 (t (do ((l (sel x data) (cdr l))) ((null l))
2310 (if (and (visiblep (car l)) (dgqf x y (car l))) (return t))))))
2312 (defun dgqf (x y f)
2313 (cond ((eq 'mgrp (caar f)) (if (eq x (cadar f)) (dgr (caddar f) y)))
2314 ((eq 'mgqp (caar f)) (if (eq x (cadar f)) (dgq (caddar f) y)))
2315 ((eq 'meqp (caar f))
2316 (if (eq x (cadar f)) (dgq (caddar f) y) (dgq (cadar f) y)))))
2318 (defun dlq (x y)
2319 (cond ((member (sel x +labs) '($neg $zero) :test #'eq) nil)
2320 ((eq '$pz (sel x +labs)) (deq x y))
2321 ((eq '$pn (sel x +labs)) (dls x y))
2322 ((dmark x '$nz) nil)
2323 ((eq x y) (setq mlqp t) nil)
2324 (t (do ((l (sel x data) (cdr l))) ((null l))
2325 (if (and (visiblep (car l)) (dlqf x y (car l))) (return t))))))
2327 (defun dlqf (x y f)
2328 (cond ((eq 'mgrp (caar f)) (if (eq x (caddar f)) (dls (cadar f) y)))
2329 ((eq 'mgqp (caar f)) (if (eq x (caddar f)) (dlq (cadar f) y)))
2330 ((eq 'meqp (caar f))
2331 (if (eq x (cadar f)) (dlq (caddar f) y) (dlq (cadar f) y)))))
2333 (defun dnq (x y)
2334 (cond ((member (sel x +labs) '($pos $neg) :test #'eq) nil)
2335 ((eq '$pz (sel x +labs)) (dgr x y))
2336 ((eq '$nz (sel x +labs)) (dls x y))
2337 ((dmark x '$pn) nil)
2338 ((eq x y) nil)
2339 (t (do ((l (sel x data) (cdr l))) ((null l))
2340 (if (and (visiblep (car l)) (dnqf x y (car l))) (return t))))))
2342 (defun dnqf (x y f)
2343 (cond ((eq 'meqp (caar f))
2344 (if (eq x (cadar f)) (dnq (caddar f) y) (dnq (cadar f) y)))))
2346 ;; mark sign of x to be m, relative to current comparison point for dcomp.
2347 ;; returns true if this fact is already known, nil otherwise.
2348 (defun dmark (x m)
2349 (cond ((eq m (sel x +labs)))
2350 ((and dbtrace (prog1
2352 (mtell (intl:gettext "DMARK: marking ~M ~M") (if (atom x) x (car x)) m))
2353 nil))
2355 (push x +labs)
2356 (push+sto (sel x +labs) m)
2357 nil)))
2359 (defun daddgr (flag x)
2360 (with-compsplt (lhs rhs x)
2361 (mdata flag 'mgrp (dintern lhs) (dintern rhs))
2362 (if (or (mnump lhs) (constant lhs))
2363 (list '(mlessp) rhs lhs)
2364 (list '(mgreaterp) lhs rhs))))
2366 (defun daddgq (flag x)
2367 (with-compsplt (lhs rhs x)
2368 (mdata flag 'mgqp (dintern lhs) (dintern rhs))
2369 (if (or (mnump lhs) (constant lhs))
2370 (list '(mleqp) rhs lhs)
2371 (list '(mgeqp) lhs rhs))))
2373 (defun daddeq (flag x)
2374 (with-compsplt-eq (lhs rhs x)
2375 (mdata flag 'meqp (dintern lhs) (dintern rhs))
2376 (list '($equal) lhs rhs)))
2378 (defun daddnq (flag x)
2379 (with-compsplt-eq (lhs rhs x)
2380 (cond ((and (mtimesp lhs) (equal rhs 0))
2381 (dolist (term (cdr lhs)) (daddnq flag term)))
2382 ((and (mexptp lhs) (mexptp rhs)
2383 (integerp (caddr lhs)) (integerp (caddr rhs))
2384 (equal (caddr lhs) (caddr rhs)))
2385 (mdata flag 'mnqp (dintern (cadr lhs)) (dintern (cadr rhs)))
2386 (cond ((not (oddp (caddr lhs)))
2387 (mdata flag 'mnqp (dintern (cadr lhs))
2388 (dintern (neg (cadr rhs)))))))
2389 (t (mdata flag 'mnqp (dintern lhs) (dintern rhs))))
2390 (list '(mnot) (list '($equal) lhs rhs))))
2392 ;; The following functions are used by asksign to write answers into the
2393 ;; database. We make sure that these answers are written into the global
2394 ;; context '$initial and not in a local context which might be generated during
2395 ;; the evaluation phase and which will be destroyed before the evaluation has
2396 ;; finshed.
2397 ;; The additional facts are removed from the global context '$initial after
2398 ;; finishing the evaluation phase of meval with a call to clearsign.
2400 (defun tdpos (x)
2401 (let ((context '$initial))
2402 (daddgr t x)
2403 (push (cons x '$pos) *local-signs*)))
2405 (defun tdneg (x)
2406 (let ((context '$initial))
2407 (daddgr t (neg x))
2408 (push (cons x '$neg) *local-signs*)))
2410 (defun tdzero (x)
2411 (let ((context '$initial))
2412 (daddeq t x)
2413 (push (cons x '$zero) *local-signs*)))
2415 (defun tdpn (x)
2416 (let ((context '$initial))
2417 (daddnq t x)
2418 (push (cons x '$pn) *local-signs*)))
2420 (defun tdpz (x)
2421 (let ((context '$initial))
2422 (daddgq t x)
2423 (push (cons x '$pz) *local-signs*)))
2425 (defun compsplt-eq (x)
2426 (with-compsplt (lhs rhs x)
2427 (when (equal lhs 0)
2428 (setq lhs rhs rhs 0))
2429 (if (and (equal rhs 0)
2430 (or (mexptp lhs)
2431 (and (not (atom lhs))
2432 (kindp (caar lhs) '$oddfun)
2433 (kindp (caar lhs) '$increasing))))
2434 (setq lhs (cadr lhs)))
2435 (values lhs rhs)))
2437 (defun mdata (flag r x y)
2438 (if flag
2439 (mfact r x y)
2440 (mkill r x y)))
2442 (defun mfact (r x y)
2443 (let ((f (datum (list r x y))))
2444 (cntxt f context)
2445 (addf f x)
2446 (addf f y)))
2448 (defun mkill (r x y)
2449 (let ((f (car (datum (list r x y)))))
2450 (kcntxt f context)
2451 (maxima-remf f x)
2452 (maxima-remf f y)))
2454 (defun mkind (x y)
2455 (kind (dintern x) (dintern y)))
2457 ;; To guess from the previous incarnation of this code,
2458 ;; each argument is assumed to be a float, bigfloat, integer, or Maxima rational.
2459 ;; Convert Maxima rationals to Lisp rationals to make them comparable to others.
2461 (defun rgrp (x y)
2462 (when (or ($bfloatp x) ($bfloatp y))
2463 (setq
2464 x (let (($float2bf t))
2465 (declare (special $float2bf))
2466 (cadr ($bfloat (sub x y))))
2467 y 0))
2468 (if (not (numberp x))
2469 (setq x (/ (cadr x) (caddr x))))
2470 (if (not (numberp y))
2471 (setq y (/ (cadr y) (caddr y))))
2472 (cond
2473 (#-ecl (> x y) #+ecl (> (- x y) 0) '$pos)
2474 (#-ecl (> y x) #+ecl (> (- y x) 0) '$neg)
2475 (t '$zero)))
2477 (defun mcons (x l)
2478 (cons (car l) (cons x (cdr l))))
2480 (defun flip (s)
2481 (cond ((eq '$pos s) '$neg)
2482 ((eq '$neg s) '$pos)
2483 ((eq '$pz s) '$nz)
2484 ((eq '$nz s) '$pz)
2485 (t s)))
2487 (defun strongp (x y)
2488 (cond ((eq '$pnz y))
2489 ((eq '$pnz x) nil)
2490 ((member y '($pz $nz $pn) :test #'eq))))
2492 (defun munformat (form)
2493 (if (atom form)
2494 form
2495 (cons (caar form) (mapcar #'munformat (cdr form)))))
2497 (defun declarekind (var prop) ; This function is for $DECLARE to use.
2498 (let (prop2)
2499 (cond ((truep (list 'kind var prop)) t)
2500 ((or (falsep (list 'kind var prop))
2501 (and (setq prop2 (assoc prop '(($integer . $noninteger)
2502 ($noninteger . $integer)
2503 ($increasing . $decreasing)
2504 ($decreasing . $increasing)
2505 ($symmetric . $antisymmetric)
2506 ($antisymmetric . $symmetric)
2507 ($rational . $irrational)
2508 ($irrational . $rational)
2509 ($oddfun . $evenfun)
2510 ($evenfun . $oddfun)) :test #'eq))
2511 (truep (list 'kind var (cdr prop2)))))
2512 (merror (intl:gettext "declare: inconsistent declaration ~:M") `(($declare) ,var ,prop)))
2513 (t (mkind var prop) t))))
2515 ;;; These functions reformat expressions to be stored in the data base.
2517 ;; Return a list of all the atoms in X that aren't either numbers or constants
2518 ;; whose numerical value we know.
2519 (defun unknown-atoms (x)
2520 (let (($listconstvars t))
2521 (remove-if (lambda (sym) (mget sym '$numer))
2522 (cdr ($listofvars x)))))
2524 ;; Rewrite a^b to a simpler expression that has the same sign:
2525 ;; If b is odd or 1/b is odd, remove the exponent, e.g. x^3 becomes x.
2526 ;; If b has a negative sign, return a^-b, e.g. 1/x^a becomes x^a.
2527 ;; Otherwise, do nothing.
2528 (defun rewrite-mexpt-retaining-sign (x)
2529 (if (mexptp x)
2530 (let ((base (cadr x)) (exponent (caddr x)))
2531 (cond ((or (eq (evod exponent) '$odd) (eq (evod (inv exponent)) '$odd)) base)
2532 ((negp exponent) (inv x))
2533 (t x)))
2536 ;; COMPSPLT
2538 ;; Split X into (values LHS RHS) so that X>0 <=> LHS > RHS. This is supposed to
2539 ;; be a canonical form for X that can be stored in the database and then looked
2540 ;; up in future.
2542 ;; This uses two worker routines: COMPSPLT-SINGLE and COMPSPLT-GENERAL. The
2543 ;; former assumes that X only contains one symbol whose value is not known (eg not %e,
2544 ;; %pi etc.). The latter tries to deal with arbitrary collections of variables.
2545 (defun compsplt (x)
2546 (multiple-value-bind (lhs rhs)
2547 (cond
2548 ((or (atom x) (atom (car x))) (values x 0))
2549 ((null (cdr (unknown-atoms x))) (compsplt-single x))
2550 (t (compsplt-general x)))
2551 (let ((swapped nil))
2552 ;; If lhs is zero, swap lhs and rhs to make the following code simpler.
2553 ;; Remember that they were swapped to swap them back afterwards.
2554 (when (equal lhs 0)
2555 (setq lhs rhs rhs 0 swapped t))
2556 (when (equal rhs 0)
2557 ;; Rewrite mexpts in factors so that e.g. x^3/y>0 becomes x*y>0. */
2558 (setq lhs
2559 (if (mtimesp lhs)
2560 (cons (car lhs) (mapcar #'rewrite-mexpt-retaining-sign (cdr lhs)))
2561 (rewrite-mexpt-retaining-sign lhs))))
2562 ;; Undo swapping lhs and rhs.
2563 (if swapped
2564 (values rhs lhs)
2565 (values lhs rhs)))))
2567 (defun compsplt-single (x)
2568 (do ((exp (list x 0)) (success nil))
2569 ((or success (symbols (cadr exp))) (values (car exp) (cadr exp)))
2570 (cond ((atom (car exp)) (setq success t))
2571 ((eq (caaar exp) 'mplus) (setq exp (splitsum exp)))
2572 ((eq (caaar exp) 'mtimes) (setq exp (splitprod exp)))
2573 (t (setq success t)))))
2575 (defun compsplt-general (x)
2576 ;; Let compsplt-single work on it first to treat constant factors/terms.
2577 (multiple-value-bind (lhs rhs) (compsplt-single x) (setq x (sub lhs rhs)))
2578 (cond
2579 ;; If x is an atom or a single level list then we won't change it any.
2580 ((or (atom x) (atom (car x)))
2581 (values x 0))
2582 ;; If x is a negative expression but not a sum, then get rid of the
2583 ;; negative sign.
2584 ((negp x) (values 0 (neg x)))
2585 ;; If x is not a sum, or is a sum with more than 2 terms, or has some
2586 ;; symbols common to both summands, then do nothing.
2587 ((or (cdddr x)
2588 (not (eq (caar x) 'mplus))
2589 (intersect* (symbols (cadr x)) (symbols (caddr x))))
2590 (values x 0))
2591 ;; -x + y gives (y, x)
2592 ((and (or (negp (cadr x)) (mnump (cadr x)))
2593 (not (negp (caddr x))))
2594 (values (caddr x) (neg (cadr x))))
2595 ;; x - y gives (x, y)
2596 ((and (not (negp (cadr x)))
2597 (or (negp (caddr x)) (mnump (caddr x))))
2598 (values (cadr x) (neg (caddr x))))
2599 ;; - x - y gives (0, x+y)
2600 ((and (negp (cadr x)) (negp (caddr x)))
2601 (values 0 (neg x)))
2602 ;; Give up! (x, 0)
2604 (values x 0))))
2606 (defun negp (x)
2607 (and (mtimesp x) (mnegp (cadr x))))
2609 (defun splitsum (exp)
2610 (do ((llist (cdar exp) (cdr llist))
2611 (lhs1 (car exp))
2612 (rhs1 (cadr exp)))
2613 ((null llist)
2614 (if (mplusp lhs1) (setq success t))
2615 (list lhs1 rhs1))
2616 (cond ((member '$inf llist :test #'eq)
2617 (setq rhs1 (add2 '$inf (sub* rhs1 (addn llist t)))
2618 lhs1 (add2 '$inf (sub* lhs1 (addn llist t)))
2619 llist nil))
2620 ((member '$minf llist :test #'eq)
2621 (setq rhs1 (add2 '$minf (sub* rhs1 (addn llist t)))
2622 lhs1 (add2 '$minf (sub* lhs1 (addn llist t)))
2623 llist nil))
2624 ((null (symbols (car llist)))
2625 (setq lhs1 (sub lhs1 (car llist))
2626 rhs1 (sub rhs1 (car llist)))))))
2628 (defun splitprod (exp)
2629 (do ((flipsign)
2630 (lhs1 (car exp))
2631 (rhs1 (cadr exp))
2632 (llist (cdar exp) (cdr llist))
2633 (sign)
2634 (minus)
2635 (evens)
2636 (odds))
2637 ((null llist)
2638 (if (mtimesp lhs1) (setq success t))
2639 (cond (flipsign
2640 (setq success t)
2641 (list rhs1 lhs1))
2642 (t (list lhs1 rhs1))))
2643 (when (null (symbols (car llist)))
2644 (sign (car llist))
2645 (if (eq sign '$neg) (setq flipsign (not flipsign)))
2646 (if (member sign '($pos $neg) :test #'eq)
2647 (setq lhs1 (div lhs1 (car llist))
2648 rhs1 (div rhs1 (car llist)))))))
2650 (defun symbols (x)
2651 (let (($listconstvars %initiallearnflag))
2652 (cdr ($listofvars x))))
2654 ;; %initiallearnflag is only necessary so that %PI, %E, etc. can be LEARNed.
2656 (defun initialize-numeric-constant (c)
2657 (setq %initiallearnflag t)
2658 (let ((context '$global))
2659 (learn `((mequal) ,c ,(mget c '$numer)) t))
2660 (setq %initiallearnflag nil))
2662 (eval-when (:load-toplevel :execute)
2663 (mapc #'true*
2664 '(;; even and odd are integer
2665 (par ($even $odd) $integer)
2667 ; Cutting out inferences for integer, rational, real, complex (DK 10/2009).
2668 ; (kind $integer $rational)
2669 ; (par ($rational $irrational) $real)
2670 ; (par ($real $imaginary) $complex)
2672 ;; imaginary is complex
2673 (kind $imaginary $complex)
2675 ;; Declarations for constants
2676 (kind $%i $noninteger)
2677 (kind $%i $imaginary)
2678 (kind $%e $noninteger)
2679 (kind $%e $real)
2680 (kind $%pi $noninteger)
2681 (kind $%pi $real)
2682 (kind $%gamma $noninteger)
2683 (kind $%gamma $real)
2684 (kind $%phi $noninteger)
2685 (kind $%phi $real)
2686 (kind $%pi $irrational)
2687 (kind $%e $irrational)
2688 (kind $%phi $irrational)
2689 (kind $%catalan $noninteger)
2690 (kind $%catalan $real)
2692 ;; Declarations for functions
2693 (kind %log $increasing)
2694 (kind %asin $increasing) (kind %asin $oddfun)
2695 (kind %atan $increasing) (kind %atan $oddfun)
2696 (kind $delta $evenfun)
2697 (kind %sinh $increasing) (kind %sinh $oddfun)
2698 (kind %cosh $posfun)
2699 (kind %tanh $increasing) (kind %tanh $oddfun)
2700 (kind %coth $oddfun)
2701 (kind %csch $oddfun)
2702 (kind %sech $posfun)
2703 (kind %asinh $increasing) (kind %asinh $oddfun)
2704 ;; It would be nice to say %acosh is $posfun, but then
2705 ;; assume(xn<0); abs(acosh(xn)) -> acosh(xn), which is wrong
2706 ;; since acosh(xn) is complex.
2707 (kind %acosh $increasing) (kind %acosh $complex)
2708 (kind %atanh $increasing) (kind %atanh $oddfun)
2709 (kind $li $complex)
2710 (kind $lambert_w $complex)
2711 (kind %cabs $complex)))
2713 ;; Create an initial context for the user which is a subcontext of $global.
2714 ($newcontext '$initial))