I have finished some work on the Complex components. The aim was to get correct
results for Complex expressions using the functions abs, cabs, carg, rectform,
realpart, imagpart, signum and conjugate.
I needed several attemps to find a way so that the complete testsuite has no
problems and we get the expected results for all Complex components.
I have done the following steps:
1. Indroducing a Complex mode for the $sign
I used the global variable complexsign and added code to the subroutines of
$sign to detect Complex symbols and expressions.
With this Complex mode of the function $sign it is possible to define a function
which works for Complex values very nice:
(defun $csign (x)
(let ((complexsign t))
($sign x)))
In addition I have implemented a routine which gives correct answers for
something like sign(a+b-2) is positive when a>1 and b>1 (I called the algorithm
sign-shift).
When complexsign is set to false Maxima don't regocnize the Complex mode and all
run as known including the testsuite. But switching the Complex mode on, csign
gives the following results:
(%i2) declare(z,complex);
(%o2) done
(%i3) csign(z);
(%o3) complex
(%i4) csign(10*z*x+y);
(%o4) complex
(%i5) csign((10*z*x+y)^2);
(%o5) complex
(%i6) csign(%i);
(%o6) imaginary
(%i7) csign(2*%i*x);
(%o7) imaginary
(%i8) csign(2*%i*x+y);
(%o8) complex
(%i9) csign((2*%i*x+y)^3);
(%o9) complex
(%i10) csign(sqrt(-1)+10);
(%o10) complex
(%i11) csign(sqrt(-1));
(%o11) imaginary
(%i12) csign(sqrt(-1)^3);
(%o12) imaginary
2. Adding step by step code at 17 places in the following routines:
rpart.lisp: absarg, risplit-exp, risplit-noun, $cabs, argnum
comm2.lisp: simpatan2
simp.lisp: simpabs, simpsignum
Furthermore I have written a testfile with 245 tests containing a lot of Complex
expressions. A special point is that the functions abs, carg, realpart,
imagpart, signum and conjugate now works together and gives the expected results
for the Complex characteristics of the Complex components.
Some remarks:
signum is extended and returns z/abs(z) for Complex expressions.
carg simplifies to expressions with atan2. The noun carg is not introduced. The
code does not give an error for atan2(0,0) but returns a noun form.
cabs(z) does not simplify to sqrt(realpart(z)^2+imagpart(z)^2) for z a Complex
symbol but returns the noun form abs(z). We often get much more simple results
when we return the noun form abs(z).
The following test were wrong and gives now correct answers:
in rtest_abs Problem 42:
abs(exp(z));
%e^(realpart(z));
in rtest_abs Problem 43:
subst(z?%i,abs(z^2));
1;
in rtest_sign: Problem 84
sign(-2+b+a);
pos;
The testsuite has no problems. (Tested with GCL 2.6.8 and CLISP 2.46).
There are also some bugs related to the above changes.
At the weekend Dan Gildea has commited code to absarg:
o absarg: handle symbols declared complex
fixes [ 620246 ] carg(complex)
[ 1269020 ] rectform(log(z)) with z declared complex
The more general code I have implemented handle these problems too. But returns
the following results:
rectform(log(u));
log(abs(u)) + %i*atan2(imagpart(u),realpart(u));
and
cabs(u);
abs(u);
Both is mathematically equivalent.
The following diffs show all changes to rpart.lisp, simp.lisp and comm2.lisp.
Some code is only for debugging and can be removed.
Index: rpart.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/rpart.lisp,v
retrieving revision 1.12
diff -u -r1.12 rpart.lisp
--- rpart.lisp 30 Aug 2008 19:10:32 -0000 1.12
+++ rpart.lisp 31 Aug 2008 11:18:45 -0000
@@ -10,6 +10,8 @@
(in-package :maxima)
+(defvar *debug-rpart* nil)
+
(macsyma-module rpart)
;;; Complex variable utilities
@@ -52,7 +54,12 @@
;;; be syntactically real without being real (e.g. sqrt(x), x<0). Thus
;;; Cabs must lead an independent existence from Abs.
-(defmfun $cabs (xx) (cabs ($rectform xx)))
+;;; CHANGE 11:
+;;; The user function $cabs calls in addition to the internal function cabs
+;;; the function rectform. That causes problems. The functions don't behave
+;;; identically. e.g. ($cabs '$infinity) doesn' work, but (cabs '$infinity)
+
+(defmfun $cabs (xx) (cabs xx))
;;; Carg gives the complex argument.
(defmfun $carg (xx)
@@ -151,6 +158,7 @@
(risplit (car l))))))
(defun risplit-expt (l)
+ (when *debug-rpart* (format t "~&RISPLIT-EXPT with ~A~%" l))
((lambda (pow $radexpand ris) ; Don't want 'simplifications' like
(cond ((fixnump pow) ; Sqrt(-x) -> %i*sqrt(x)
((lambda (sp)
@@ -160,6 +168,8 @@
(mul -1 (div (cdr sp) a2+b2))))
(spabs sp)))
((> (abs pow) $maxposex)
+ (when *debug-rpart*
+ (format t "~&RISPLIT-EXPT with ~A and ~A~%" sp pow))
(cond ((=0 (cdr sp)) (cons (powers (car sp) pow) 0))
(t ((lambda (abs^n natan)
(cons (mul abs^n
@@ -182,11 +192,17 @@
(not (> (cadr pow) $maxposex))
(prog2 (setq ris (risplit (cadr l)))
(or (= (caddr pow) 2) (=0 (cdr ris)))))
+ (when *debug-rpart*
+ (format t "~&RISPLIT-EXPT in COND ratnump with ~A and ~A~%" ris
pow))
(cond ((=0 (cdr ris))
- (case (cond ((mnegp (car ris)) '$negative)
- (implicit-real '$positive)
- (t (asksign (car ris))))
- ($negative (risplit (mul2 (power -1 pow) (power (neg (car ris))
pow))))
+ ;; The base is real.
+ (case (cond ((mnegp (car ris)) '$neg)
+ (implicit-real '$pos)
+ ;; CHANGE 10:
+ ;; We don't ask for the sign but use $sign.
+ ;; Is this code general enough?
+ (t ($sign (car ris))))
+ ($neg (risplit (mul2 (power -1 pow) (power (neg (car ris)) pow))))
($zero (cons (power 0 pow) 0))
(t (cons (power (car ris) pow) 0))))
(t ((lambda (abs2 n pos?)
@@ -206,6 +222,10 @@
(risplit ((lambda ($numer) (exptrl ris pow)) t)))
(t ((lambda (sp aa)
;;If all else fails, we use the trigonometric form.
+ (when *debug-rpart*
+ (format t "~&RISPLIT-EXPT the general case~%")
+ (format t "~& : base ~A~%" aa)
+ (format t "~& : power ~A~%" sp))
((lambda (pre post)
(cons (mul pre (take '(%cos) post))
(mul pre (take '(%sin) post))))
@@ -213,11 +233,27 @@
(powers (car aa) (car sp)))
(add (mul (cdr sp) (take '(%log) (car aa)))
(mul (car sp) (cdr aa)))))
- (risplit (caddr l)) (absarg1 (cadr l))))))
+ ;; CHANGE 6:
+ ;; absarg1 tries to give an answer on the assumption that
+ ;; the base is real valued. The problem is that we have to
+ ;; take into account a phase factor for negative real values.
+ ;; We change absarg1 to the more general function absarg.
+ (risplit (caddr l)) (absarg (cadr l))))))
(caddr l) nil nil))
(defun risplit-noun (l)
- (cons (simplify (list '(%realpart) l)) (simplify (list '(%imagpart) l))))
+ (when *debug-rpart* (format t "~&RISPLIT-NOUN with ~A~%" l))
+ ;; CHANGE 8:
+ ;; We are more complete and distinguish real, imaginary and complex nouns.
+ ;; We call risplit-noun in this file only for symbols which are
+ ;; declared complex. So we do here more as needed.
+ (let ((sgn ($csign l)))
+ (cond ((eq sgn '$imaginary)
+ (cons 0 (simplify (list '(%imagpart) l))))
+ ((eq sgn '$complex)
+ (cons (simplify (list '(%realpart) l))
+ (simplify (list '(%imagpart) l))))
+ (t (cons l 0)))))
(defun absarg1 (arg)
(let ((arg1 arg) ($keepfloat t))
@@ -238,6 +274,7 @@
;;; (<Real part> . <imaginary part>).
(defun risplit (l)
+ (when *debug-rpart* (format t "~&RISPLIT with ~A~%" l))
(let (($domain '$complex) ($m1pbranch t) $logarc op)
(cond ((atom l)
(cond ((eq l '$%i) (cons 0 1))
@@ -428,9 +465,15 @@
;;; Rectify into polar form; Arguments similar to risplit
(defun argnum (n)
- (cond ((minusp n) (simplify '$%pi)) (t 0)))
+ (cond ((minusp n) (simplify '$%pi))
+ ;; CHANGE 16:
+ ;; Simplify 0 to atan2(0,0). Mathematically atan2(0,0) is not
+ ;; definied. If we simplify to 0 we can get wrong results.
+ ((zerop1 n) (list '(atan2 simp) 0 0))
+ (t 0)))
(defun absarg (l)
+ (when *debug-rpart* (format t "~&ABSARG with ~A~%" l))
(setq l ($expand l))
(cond ((atom l)
(cond ((eq l '$%i)
@@ -439,16 +482,31 @@
(cons (abs l) (argnum l)))
((member l '($%e $%pi) :test #'eq) (cons l 0))
((eq l '$infinity) (cons '$inf '$ind))
- ((kindp l '$complex) (cons (list '($cabs) l)
- (list '($carg) l)))
- (absflag (cons (take '(mabs) l) 0))
- (t ((lambda (gs)
- (cond ((eq gs '$positive) (cons l 0))
- ((eq gs '$zero) (cons 0 0))
- ((eq gs '$negative)
- (cons (neg l) (simplify '$%pi)))
- (t (cons (take '(mabs) l) 0))))
- (cond ((eq rischp l) '$positive) (t (asksign l)))))))
+
+ ;; Change 14:
+ ;; We simplify the real infinities.
+ ((member l '($inf $minf)) (cons '$inf (genatan 0 l)))
+
+; ((kindp l '$complex) (cons (list '($cabs) l)
+; (list '($carg) l)))
+;;; CHANGE 13: We don't do this simplification because we get wrong results.
+;;;
+; (absflag (cons (take '(mabs) l) 0))
+ (t
+ ;; CHANGE 1: Testing with $csign and not asksign
+ (let ((gs (cond ((eq rischp l) '$positive) (t ($csign l)))))
+ (cond ;; CHANGE 3: Adding the test $pz
+ ((member gs '($pos $pz)) (cons l 0))
+ ((eq gs '$zero) (cons 0 0))
+ ((eq gs '$neg) (cons (neg l) '$%pi))
+ ;; CHANGE 2: Answer for a complex symbol.
+ ((member gs '($complex $imaginary))
+ (cons (list '(mabs simp) l)
+ (genatan ($imagpart l) ($realpart l))))
+ ;; CHANGE 7:
+ ;; We have a real symbol but don't know the sign.
+ ;; We return a more general result.
+ (t (cons (take '(mabs) l) (genatan 0 l))))))))
((member (caar l) '(rat bigfloat) :test #'eq)
(cons (list (car l) (abs (cadr l)) (caddr l))
(argnum (cadr l))))
@@ -462,10 +520,15 @@
(return (cons (muln absl t)
(2pistrip (addn argl t))))))
(setq abars (absarg (car n)))))
- ((eq (caar l) 'mexpt)
- (let ((aa (absarg (cadr l)))
+ ((eq (caar l) 'mexpt)
+ ;; CHANGE 4: Setting absflag nil
+ (let ((aa (let ((absflag nil)) (absarg (cadr l))))
(sp (risplit (caddr l)))
($radexpand nil))
+ (when *debug-rpart*
+ (format t "~&ABSARG in MEXPT with ~A~%" l)
+ (format t "~& : base = ~A~%" aa)
+ (format t "~& : pow = ~A~%" sp))
(cons (mul (powers (car aa) (car sp))
(powers '$%e (neg (mul (cdr aa) (cdr sp)))))
(add (mul (cdr aa) (car sp))
@@ -496,6 +559,8 @@
(and foot (not (=0 (cdr (risplit (cadr l))))) (absarg foot)))
(coversinemyfoot l)))
(t (let ((ris (trisplit l)))
+ (when *debug-rpart*
+ (format t "~&ABSARG the general case with ~A~%" ris))
(xcons
;;; Arguments must be in this order so that the side-effect of the Atan2,
;;; that is, determining the Asksign of the argument, can happen before
@@ -515,6 +580,7 @@
(take '(%atan) (m// num den)))))
(defun absarg-mabs (l)
+ (when *debug-rpart* (format t "~&ABSARG-MABS with ~A~%" l))
(if (eq (csign l) t)
(if (member (caar l) '(mabs %cabs) :test #'eq) l (list '(%cabs simp) l))
(take '(mabs) l)))
Success, CVS operation completed
Index: simp.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/simp.lisp,v
retrieving revision 1.58
diff -u -r1.58 simp.lisp
--- simp.lisp 24 Aug 2008 16:11:33 -0000 1.58
+++ simp.lisp 31 Aug 2008 15:10:19 -0000
@@ -880,7 +880,34 @@
((taylorize 'mabs (second x)))
((member y '($inf $infinity $minf) :test #'eq) '$inf)
((member y '($ind $und) :test #'eq) y)
- ((eq (setq z (csign y)) t) (cabs y))
+
+; Problem: if do the multiplication at this point we get no more
+; the same results for abs(z*(1+%i)) and abs(z+z*%i)
+; ((mtimesp y)
+; (muln (mapcar #'(lambda (u) (simplify (list '(mabs) u))) (cdr y)) t))
+
+ ;; CHANGE 12:
+ ;; We test for complex expressions with $csign
+ ((member (setq z ($csign y)) '($complex $imaginary))
+ ;; But we don't call cabs for complex symbols
+ ;; because we would get an endless loop.
+ (cond ((symbolp y)
+ (cond ((eq y '$%i) 1)
+ (t (eqtest (list '(mabs) y) x))))
+ (t
+ ;; call cabs if y is not a symbol
+ (cabs y))))
+
+ ;; We have a problem with Problem 57 in rtest16.mac:
+ ;; There is one test in the testsuite which calls simpabs with an
+ ;; abs function containing an imaginary argument. $csign give the
+ ;; correct answer $PZ because abs is real valued. But the code
+ ;; of the test depends on the wrong answer that the expression is
+ ;; complex. We add the old part of the code. So we enter cabs for
+ ;; some more cases and the test no longer fails.
+
+ ((eq (setq z (csign y)) t) (cabs y))
+
((member z '($pos $pz) :test #'eq) y)
((member z '($neg $nz) :test #'eq) (neg y))
((eq z '$zero) 0)
@@ -1370,7 +1397,17 @@
(setq y (simpcheck (cadr x) z))
(cond ((mnump y)
(setq y (num1 y)) (cond ((plusp y) 1) ((minusp y) -1) (t 0)))
- ((eq (setq z (csign y)) t) (eqtest (list '(%signum) y) x))
+
+ ;; Change 17:
+ ;; We handle infinity, ind, und
+ ((and (atom y) (eq y '$infinity)) '$und)
+ ((and (atom y) (member y '($und $ind))) y)
+
+ ;; Change 15:
+ ;; For a complex expression z the definition is z/abs(z)
+ ((member (setq z ($csign y)) '($complex $imaginary))
+ (div y (simplify (list '(mabs) y))))
+
((eq z '$pos) 1)
((eq z '$neg) -1)
((eq z '$zero) 0)
Success, CVS operation completed
Index: comm2.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/comm2.lisp,v
retrieving revision 1.21
diff -u -r1.21 comm2.lisp
--- comm2.lisp 12 Jul 2008 11:59:03 -0000 1.21
+++ comm2.lisp 31 Aug 2008 11:24:15 -0000
@@ -399,13 +399,19 @@
(declare-top (special $numer $%piargs $logarc $trigsign half%pi fourth%pi))
+(defvar *debug-simpatan2* nil)
+
(defun simpatan2 (e vestigial z) ; atan2(y,x) ~ atan(y/x)
(declare (ignore vestigial))
(twoargcheck e)
(let (y x signy signx)
(setq y (simpcheck (cadr e) z) x (simpcheck (caddr e) z))
(cond ((and (zerop1 y) (zerop1 x))
- (merror "atan2(0,0) has been generated."))
+ ;; CHANGE 17:
+ ;; Don't generate an error but return a noun form.
+ ;; So the calculation can go on with a noun form.
+; (merror "atan2(0,0) has been generated."))
+ (list '(atan2 simp) 0 0))
( ;; float contagion
(and (or (numberp x) (ratnump x)) ; both numbers
(or (numberp y) (ratnump y)) ; ...but not bigfloats
@@ -420,9 +426,27 @@
(*fpatan y (list x)))
((and $%piargs (free x '$%i) (free y '$%i)
;; Only use asksign if %piargs is on.
- (cond ((zerop1 y) (if (atan2negp x) (simplify '$%pi) 0))
+ (cond ((zerop1 y)
+ ;; CHANGE 6:
+ ;; The imaginary part is zero. We try to simplify
+ ;; further for the realpart but don't use asksign.
+ ;; The testsuite execute this code 4 times.
+; (if (atan2negp x) (simplify '$%pi) 0)
+ (when *debug-simpatan2*
+ (format t "~&SIMPATAN in COND ZEROP1 y with x=~A" x))
+ (cond ((eq ($sign x) '$neg) '$%pi)
+ ((eq ($sign x) '$pos ) 0)
+ ;; The general case: We return a noun form.
+ (t (list '(atan2 simp) y x))))
((zerop1 x)
- (if (atan2negp y) (mul2* -1 half%pi) (simplify half%pi)))
+ ;; CHANGE 9:
+ ;; The realpart is zero. We try to simplify
+ ;; further for the imagpart but don't use asksign.
+; (if (atan2negp y) (mul2* -1 half%pi) (simplify half%pi)))
+ (cond ((eq ($sign y) '$neg) (mul -1 (div '$%pi 2)))
+ ((eq ($sign x) '$pos ) (div '$%pi 2))
+ ;; The general case: We return a noun form.
+ (t (list '(atan2 simp) y x))))
((alike1 y x)
;; Should we check if ($sign x) is $zero here?
(if (atan2negp x) (mul2* -3 fourth%pi) (simplify fourth%pi)))
Success, CVS operation completed
At last the code for the function $csign:
Index: compar.lisp
===================================================================
RCS file: /cvsroot/maxima/maxima/src/compar.lisp,v
retrieving revision 1.38
diff -u -r1.38 compar.lisp
--- compar.lisp 27 Jul 2008 07:04:06 -0000 1.38
+++ compar.lisp 26 Aug 2008 21:08:35 -0000
@@ -15,10 +15,11 @@
(load-macsyma-macros mrgmac)
(declare-top (special $float2bf $radexpand $ratprint $ratsimpexpons
$listconstvars
- success $props *x* $%enumer)
+ success $props *x* $%enumer complexsign)
;; Variables defined in DB
(special context current dobjects dbtrace +labs))
+(defvar *debug-compar* nil)
(defvar %initiallearnflag)
(defmvar $context '$initial
@@ -44,10 +45,9 @@
If NIL THROWs in that case."
no-reset)
-(defmvar complexsign nil
+(defvar complexsign nil
"If T, COMPAR attempts to work in a complex mode.
- This scheme is only very partially developed at this time."
- no-reset)
+ This scheme is only very partially developed at this time.")
(defmvar $prederror nil)
(defmvar $signbfloat t)
@@ -656,14 +656,27 @@
(or (not (free x '$%i))
(let (sign-imag-errp limitp) (catch 'sign-imag-err ($sign x)))))
+(defvar *debug-$csign* nil)
+
+(defmfun $csign (x)
+ (let ((complexsign t))
+ (when *debug-$csign*
+ (format t "~&$CSIGN called with ~A~%" x))
+ ($sign x)))
+
(defmfun $sign (x)
+ (when (or *debug-compar* *debug-$csign*)
+ (format t "~&$SIGN called with ~A and complexsign ~A~%" x complexsign))
(let ((x (specrepcheck x))
sign minus odds evens factored)
(sign01 (cond (limitp (restorelim x))
- ((not (free x '$%i)) ($rectform x))
+ ;; We don't do rectform in Complex Mode
+ ;; In Real mode the code stops finding a %i
+ ((and (not complexsign) (not (free x '$%i))) ($rectform x))
(t x)))))
(defun sign01 (a)
+ (when *debug-compar* (format t "~&SIGN01 with ~A~%" a))
(let ((e (sign-prep a)))
(cond ((eq e '$pnz) '$pnz)
(t (setq e (sign1 e))
@@ -671,6 +684,7 @@
;;; Preparation for asking questions from DEFINT or LIMIT.
(defun sign-prep (x)
+ (when *debug-compar* (format t "~&SIGN-PREP with ~A~%" x))
(if limitp
(destructuring-let (((rpart . ipart) (trisplit x)))
(cond ((and (equal (sratsimp ipart) 0)
@@ -1025,6 +1039,7 @@
(defun sign1 (x)
(setq x (specrepcheck x))
(setq x (infsimp* x))
+ (when *debug-compar* (format t "~&SIGN1 with ~A~%" x))
(if (member x '($und $ind $infinity) :test #'eq)
(if limitp '$pnz (merror "The sign of ~:M is undefined" x)))
(prog (dum exp)
@@ -1082,8 +1097,12 @@
(if (eq ans 'numer) (setq ans 'symbol)))
(t (return nil)))))))
+
(defmfun sign (x)
+ (when *debug-compar* (format t "~&SIGN with ~A~%" x))
(cond ((mnump x) (setq sign (rgrp x 0) minus nil odds nil evens nil))
+ ;; In Complex mode return $imaginary
+ ((and complexsign (atom x) (eq x '$%i)) (setq sign '$imaginary))
((atom x) (if (eq x '$%i) (imag-err x)) (sign-any x))
((eq (caar x) 'mtimes) (sign-mtimes x))
((eq (caar x) 'mplus) (sign-mplus x))
@@ -1100,21 +1119,51 @@
(sign-oddinc x))
(t (sign-any x))))
+;;; Is the symbol declared to be Complex
+;;; Maxima has more definitions of similar functions.
+
+(defun declared-complex-p (sym)
+ (and (symbolp sym)
+ (kindp sym '$complex)
+ (not (or (kindp sym '$real)
+ (kindp sym '$rational)
+ (kindp sym '$integer)))))
+
(defun sign-any (x)
- (dcompare x 0)
- (if (and $assume_pos
- (member sign '($pnz $pz $pn) :test #'eq)
- (if $assume_pos_pred (let ((*x* x)) (is '(($assume_pos_pred) *x*)))
- (mapatom x)))
- (setq sign '$pos))
- (setq minus nil evens nil
- odds (if (not (member sign '($pos $neg $zero) :test #'eq)) (ncons x))))
+ (when *debug-compar* (format t "~&SIGN-ANY with ~A~%" x))
+ (cond
+ ((and complexsign (or (and (atom x) (declared-complex-p x))
+ (and (listp x) (declared-complex-p (caar x)))))
+ (when *debug-compar* (format t "~& : Symbol declared complex found.~%"))
+ (setq minus nil evens nil odds nil)
+ (if ($featurep x '$imaginary)
+ (setq sign '$imaginary)
+ (setq sign '$complex)))
+
+ (t
+ (dcompare x 0)
+ (if (and $assume_pos
+ (member sign '($pnz $pz $pn) :test #'eq)
+ (if $assume_pos_pred
+ (let ((*x* x)) (is '(($assume_pos_pred) *x*)))
+ (mapatom x)))
+ (setq sign '$pos))
+ (setq minus nil evens nil
+ odds (if (not (member sign '($pos $neg $zero) :test #'eq))
+ (ncons x))))))
(defun sign-mtimes (x)
+ (when *debug-compar* (format t "~&SIGN-MTIMES with ~A~%" x))
(setq x (cdr x))
(do ((s '$pos) (m) (o) (e)) ((null x) (setq sign s minus m odds o evens e))
(sign1 (car x))
+ (when *debug-compar* (format t "~&SIGN-MTIMES test for ~A~%" x))
(cond ((eq sign '$zero) (return t))
+
+ ((and complexsign (eq sign '$complex)) (return t))
+ ((and complexsign (eq sign '$imaginary)) (setq s sign))
+ ((and complexsign (eq s '$imaginary)))
+
((eq sign '$pos))
((eq sign '$neg) (setq s (flip s) m (not m)))
((prog2 (setq m (not (eq m minus)) o (nconc odds o) e (nconc evens e))
@@ -1129,6 +1178,7 @@
(setq x (cdr x))))
(defun sign-mplus (x &aux s o e m)
+ (when *debug-compar* (format t "~&SIGN-MPLUS with ~A~%" x))
(cond ((signdiff x))
((prog2 (setq s sign e evens o odds m minus) nil))
((signsum x))
@@ -1146,6 +1196,10 @@
(null (cdddr lhs))
(negp (cadr lhs)) (not (negp (caddr lhs))))
(setq rhs (neg (cadr lhs)) lhs (caddr lhs)))
+ (when *debug-compar*
+ (format t "~&SIGNDIFF with ~A~%" x)
+ (format t "~& : lhs = ~A~%" lhs)
+ (format t "~& : rhs = ~A~%" rhs))
(let (dum)
(cond ((or (equal rhs 0) (mplusp lhs)) nil)
((and (member (constp rhs) '(numer symbol) :test #'eq)
@@ -1195,21 +1249,89 @@
(when sgn (setq sign sgn minus nil odds nil evens nil)
t)))
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
+;;; If a symbol has a bound e.g. a > 1 or b < -1 the interval for the tests
+;;; in signsum is accordingly shifted by substituting the shifted expression
+;;; e. g. (a+1) or (b-1) into the expression. This shift of the values allows
+;;; correct answers for expressions like is(a+b>2) when a>1 and b>1
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
+
+(defun sign-shift (expr)
+ (when *debug-compar* (format t "~&SIGN-SHIFT with ~A~%" expr))
+ (do ((l (cdr ($facts)) (cdr l))
+ (e expr)
+ (flag) (fact) (num))
+ ((null l) (if flag ($expand e) expr))
+ (setq fact (car l))
+ (when *debug-compar* (format t "~& : Testing fact ~A" fact))
+ (cond
+ ((member (caar fact) '(mgreaterp mgeqp))
+ (cond
+ ((and (symbolp (cadr fact))
+ (not (member (cadr fact) '($%pi $%e $%gamma $%phi)))
+ (mnump (setq num (caddr fact)))
+ (not (zerop1 num))
+ (not ($freeof (cadr fact) e)))
+ (when *debug-compar*
+ (format t "~&Substitute ~A with ~A~%"
+ (cadr fact)
+ (add (cadr fact) (caddr fact))))
+ (setq flag t)
+ (if (mminusp num) (setq num (mul -1 num)))
+ (setq e ($substitute (add (cadr fact) num) (cadr fact) e)))
+ ((and (symbolp (caddr fact))
+ (not (member (caddr fact) '($%pi $%e $%gamma $%phi)))
+ (mnump (setq num (cadr fact)))
+ (not (zerop1 num))
+ (not ($freeof (caddr fact) e)))
+ (when *debug-compar*
+ (format t "~&Substitute ~A with ~A~%"
+ (caddr fact)
+ (sub (caddr fact) (cadr fact))))
+ (setq flag t)
+ (if (mminusp num) (setq num (mul -1 num)))
+ (setq e ($substitute (sub (caddr fact) num) (caddr fact) e))))))))
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;
+
(defun signsum (x)
+ (setq x (sign-shift x)) ; Check if we should shift for symbolys with a rule.
+ ;; x might be simplified to an atom in sign-shift
+ (when (atom x) (setq x (cons '(mplus) (list x))))
+ (when *debug-compar* (format t "~&SIGNSUM with ~A~%" x))
(do ((l (cdr x) (cdr l)) (s '$zero))
- ((null l) (setq sign s minus nil odds (list x) evens nil) t)
+ ((null l)
+ ;; We have to return nil if don't find a sign
+ (setq sign s minus nil odds (list x) evens nil)
+ (if (eq sign '$pnz) nil t))
(sign (car l))
- (cond ((or (and (eq sign '$zero)
+ (when *debug-compar* (format t "~&SIGNSUM with sign ~A~%" sign))
+ (cond ((and complexsign (or (eq sign '$complex) (eq sign '$imaginary)))
+ ;; We have found a complex or imaginary value.
+ ;; the result is complex. We are finish.
+ (when *debug-compar* (format t "~&Complex in SIGNSUM found.~%"))
+ (setq sign '$complex odds nil evens nil minus nil)
+ (return t))
+
+ ((or (and (eq sign '$zero)
(setq x (sub x (car l))))
(and (eq s sign) (not (eq s '$pn))) ; $PN + $PN = $PNZ
(and (eq s '$pos) (eq sign '$pz))
(and (eq s '$neg) (eq sign '$nz))))
- ((or (and (member sign '($pz $pos) :test #'eq) (member s '($zero $pz)
:test #'eq))
- (and (member sign '($nz $neg) :test #'eq) (member s '($zero $nz)
:test #'eq))
+ ((or (and (member sign '($pz $pos) :test #'eq)
+ (member s '($zero $pz) :test #'eq))
+ (and (member sign '($nz $neg) :test #'eq)
+ (member s '($zero $nz) :test #'eq))
(and (eq sign '$pn) (eq s '$zero)))
(setq s sign))
- (t (setq sign '$pnz odds (list x) evens nil minus nil)
- (return nil)))))
+ (t
+ (cond
+ (complexsign
+ ;; In Complex Mode we have to look further for a Complex.
+ (setq s '$pnz))
+ (t
+ (setq sign '$pnz odds (list x) evens nil minus nil)
+ (return nil)))))))
(defun signfactor (x)
(let (y (factored t))
@@ -1225,13 +1347,34 @@
(let ($ratprint)
(factor x)) x))
-(defmvar complexsign nil)
-
(defun sign-mexpt (x)
(let* ((expt (caddr x)) (base1 (cadr x))
(sign-expt (sign1 expt)) (sign-base (sign1 base1))
(evod (evod expt)))
- (cond ((and (eq sign-base '$zero)
+
+ (when *debug-compar*
+ (format t "~&SIGN-MEXPT with ~A~%" x)
+ (format t "~& : expt = ~A~%" expt)
+ (format t "~& : sign-expt = ~A~%" sign-expt)
+ (format t "~& : base = ~A~%" base1)
+ (format t "~& : sign-base = ~A~%" sign-base)
+ (format t "~& : evod = ~A~%" evod))
+
+ (cond ((and complexsign (or (eq sign-expt '$complex)
+ (eq sign-expt '$imaginary)
+ (eq sign-base '$complex)
+ (eq sign-base '$imaginary)))
+ (setq sign '$complex))
+
+ ((and complexsign
+ (eq sign-base '$neg)
+ (eq (evod ($expand (mul 2 expt))) '$odd))
+ ;; When the double of the exponent is odd we get %i
+ ;; if the base is negative.
+ (when *debug-compar* (format t "half integral Exponent detected."))
+ (setq sign '$imaginary))
+
+ ((and (eq sign-base '$zero)
(member sign-expt '($zero $neg) :test #'eq))
(dbzs-err x))
((eq sign-expt '$zero) (setq sign '$pos) (tdzero (sub x 1)))
@@ -1268,8 +1411,10 @@
(setq evens (nconc odds evens)
odds nil minus nil))
((mevenp (caddr expt))
- (cond (complexsign
- (setq sign-base (setq sign-expt '$pnz)))
+ (cond ;; at this point we have a problem
+ ;; the results differ for ...
+ ((and complexsign (eq sign-base '$neg))
+ (setq sign-base (setq sign-expt '$complex)))
((eq sign-base '$neg) (imag-err x))
((eq sign-base '$pn)
(setq sign-base '$pos)
@@ -1564,10 +1709,18 @@
`(list , at x))
(defun dcompare (x y)
+ (when *debug-compar*
+ (format t "~&DCOMPARE:~%")
+ (format t "~& : x = ~A~%" x)
+ (format t "~& : y = ~A~%" y))
(setq odds (list (sub x y)) evens nil minus nil
sign (cond ((eq x y) '$zero)
((or (eq '$inf x) (eq '$minf y)) '$pos)
((or (eq '$minf x) (eq '$inf y)) '$neg)
+ ((and complexsign
+ (or (declared-complex-p x) (declared-complex-p y)))
+ (when *debug-compar* (format t "& : Complex found.~%"))
+ '$complex)
(t (dcomp x y)))))
(defun dcomp (x y)
Success, CVS operation completed
What do you think? Should I commit the changes?
Dieter Kaiser