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