1;; Author: Barton Willis with help from Richard Fateman 2 3#| 4To simplify a sum with n terms, the standard simplus function calls 5great O(n^2) times. By using sorting more effectively, this code 6reduces the calls to great to O(n log_2(n)). 7 8Also, this code tries to be "infinity correct" for addition. By this I 9mean inf + inf --> inf, inf + number --> inf, and inf + minf --> und. 10Since -1 * inf doesn't simplify to minf, this code doesn't simplify 11inf - inf to und; consequently, this portion of the code is largely 12untested. There are other problems too. For one, this code does 13f(inf) - f(inf) --> 0 (comment from Stavros Macrakis). I don't know 14how far we can go with such things without making a mess. You could 15argue that Maxima should do 16 17 f(x) - f(x) --> if finitep(f(x)) then 0 else und 18 19instead of f(x) - f(x) --> 0. 20 21There is a great deal more that could be done. We could tag each 22infinity with a unique identifier (idea from RJF). That way we could 23have x : inf, x - x --> 0 and x : inf, y : inf, x - y --> und, for 24example. 25 26In short, this code is highly experimental; you should not use it for 27anything that is important. I place it in /share/contrib because 28others encouraged me too; also I thought it might be useful to others 29who would like to experiment with similar code. Since a great deal of 30work has gone into the current simplus code, I'm not sure that a total 31re-write is the best route. 32 33Maybe the special case dispatch part of this code makes the task of 34extending simplus (to intervals, for example) easier. The basic design 35of this code is due to Barton Willis. 36|# 37 38#| Fixed bugs: 39 40(1) ceiling(asin(-107) -42) <--- bug! Gets stuck. I think -1.0 * (complex) should 41 expand, but it doesn't. I fixed this by changing the condition for a "do-over" from 42 (and (or (equalp cf 1) (equalp cf -1)) (mplusp x) ...) to (and (or (eq cf 1) (eq cf -1)) (mplusp x) ...) 43 44(2) rat(x) + taylor(x^42,x,0,1) --> error. Fixed by adding taylor terms separately from mrat terms. 45 46Maxima 5.17.0 bugs: I think the rtest16 bugs 74 and 121 are related to the fact that 47-1 * inf doesn't simplify to minf. Fixing this requires a new simptimes, I think. 48 49 Errors found in rtest15.mac, problems: (189 222) <-- correct, but differ from expected 50 Errors found in rtest16.mac, problems: (74 121) <-- wrong and do asksign 51 Error found in rtestsum.mac, problem: (226) <-- not wrong but differs from expected 52 Errors found in rtest_expintegral.mac, problems: (133 134) <-- small differences in big float 53 54Unfixed: 55 56 (1) sqrt(3) + sqrt(3) + sqrt(3) --> 3^(3/2), but altsimp does 3 * sqrt(3). I'm not so sure 57 we want sqrt(3) + sqrt(3) + sqrt(3) --> 3^(3/2). 58 59|# 60 61(in-package :maxima) 62(declaim (optimize (speed 3)(safety 0))) 63 64(define-modify-macro mincf (&optional (i 1)) addk) 65 66(defmacro opcons (op &rest args) 67 `(simplify (list (list ,op) ,@args))) 68 69(defmacro opapply (op args) 70 `(simplify (cons (list ,op) ,args))) 71 72(defun mzerop (z) 73 (and (mnump z) 74 (or (and (numberp z)(= z 0)) 75 (and (bigfloatp z)(= (cadr z) 0))))) ;bigfloat zeros may be diff precisions 76 77(defun convert-to-coeff-form (x) 78 (let ((c)) 79 (cond ((mnump x) (cons 1 x)) 80 ((mtimesp x) 81 (pop x) ;remove (car x) which is (mtimes ..) 82 (cond ((mnump (setf c (car x))) ;set c to numeric coeff. 83 (pop x) ; remove numeric coeff. 84 (if (null (cdr x));; if only one more item, that's it. 85 (cons (car x) c) 86 (cons `((mtimes simp) ,@x) c))) 87 (t (cons `((mtimes simp) ,@x) 1)))) 88 (t (cons x 1))))) 89 90;; The expression e must be simplified (ok) 91;; (a) 1 * x --> x, 92;; (b) 0 * x --> 0, 0.0 * x --> 0.0, 0.0b0 * x --> 0.0b0 93;; (c) cf * e --> timesk(ck,e) when e is a maxima number, 94;; (d) -1 * (a + b) --> -a - b, 95;; (e) cf * (* a b c) --> (* (* cf a) b c ...) when a is a number; otherwise (* cf a b ...) 96;; (f) (* cf e) (default) 97 98(defun number-times-expr (cf e) 99 (cond ((eql cf 1) e) 100 ((mzerop cf) cf) 101 ((mnump e) (timesk cf e)) ; didn't think this should happen 102 ((and (onep1 (neg cf)) (mplusp e)) 103 (opapply 'mplus (mapcar 'neg (cdr e)))) 104 ((mtimesp e) 105 (if (mnump (cadr e)) 106 `((mtimes simp) ,@(cons (timesk cf (cadr e)) (cddr e))) 107 `((mtimes simp) ,@(cons cf (cdr e))))) 108 (t `((mtimes simp) ,cf ,e)))) 109 110;; Add an expression x to a list of equalities l. 111 112(defun add-expr-mequal (x l) 113 (setq l (mapcar 'cdr l)) 114 (push (list x x) l) 115 (setq l (list (reduce 'add (mapcar 'first l)) (reduce 'add (mapcar 'second l)))) 116 (simplifya (cons '(mequal) l) t)) 117 118(defun add-expr-mrat (x l) 119 (ratf (cons '(mplus) (cons (ratf x) l)))) 120 121(defun add-expr-taylor (x l) 122 ($taylor (cons '(mplus) (cons x l)))) 123 124(defun add-expr-mlist (x l) 125 (setq l (if (cdr l) (reduce 'addmx l) (car l))) 126 (opapply 'mlist (mapcar #'(lambda (s) (add x s)) (cdr l)))) 127 128;; Simple demo showing how to define addition for a new object. 129;; We could append simplification rules for intervals: 130 131;; (a) interval(a,a) --> a, 132;; (b) if p > q then interval(p,q) --> standardized empty interval? 133 134(defun add-expr-interval (x l) 135 (setq l (mapcar #'(lambda (s) `((mlist) ,@(cdr s))) l)) 136 (setq l (if (cdr l) (reduce 'addmx l) (car l))) 137 (opapply '$interval (mapcar #'(lambda (s) (add x s)) (cdr l)))) 138 139;; Add an expression x to a list of matrices l. 140 141(defun add-expr-matrix (x l) 142 (mxplusc x (if (cdr l) (reduce 'addmx l) (car l)))) 143 144;; Return a + b, where a, b in {minf, inf, ind, und, infinity}. I should 145;; extend this to allow zeroa and zerob (but I'm not sure zeroa and zerob 146;; are supposed to be allowed outside the limit code). 147 148(defun add-extended-real (a b) 149 (cond ((eq a '$minf) 150 (cond ((memq b '($minf $ind)) '$minf) 151 ((memq b '($und $inf)) '$und) 152 ((eq b '$infinity) '$infinity))) 153 ((eq a '$ind) 154 (cond ((eq b '$minf) '$minf) 155 ((eq b '$ind) '$ind) 156 ((eq b '$und) '$und) 157 ((eq b '$inf) '$inf) 158 ((eq b '$infinity) '$infinity))) 159 ((eq a '$und) '$und) 160 ((eq a '$inf) 161 (cond ((memq b '($minf $und)) '$und) 162 ((memq b '($inf $ind)) '$inf) 163 ((eq b '$infinity) '$infinity))) 164 ((eq a '$infinity) (if (eq b '$und) '$und '$infinity)))) 165 166;; Add an expression x to a list of infinities. 167 168(defun add-expr-infinities (x l) 169 (setq l (if l (reduce 'add-extended-real l) (car l))) 170 (if (mnump x) l `((mplus simp) ,x ,l))) 171 172;; I assumed that if a list of distinct members is sorted using great, 173;; then it's still sorted after multiplying each list member by a nonzero 174;; maxima number. I'm not sure this is true. 175 176;; If l has n summands, simplus calls great O(n log_2(n)) times. All 177;; other spendy functions are called O(n) times. The standard simplus 178;; function calls great O(n^2) times, I think. 179 180;(defvar *calls-to-simplus* 0) 181;(defvar *simplus-length* 0) 182;(defvar *its-an-atom* 0) 183;(defvar *not-an-atom* 0) 184 185 186(defun simplus (l w z) 187 (declare (ignore w)) 188 ;;(incf *calls-to-simplus*) 189 ;;(if (> 8 (length l)) (incf *simplus-length*)) 190 (let ((acc nil) (cf) (x) (num-sum 0) (do-over nil) (mequal-terms nil) (mrat-terms nil) 191 (inf-terms nil) (matrix-terms nil) (mlist-terms nil) (taylor-terms nil) (interval-terms nil) (op) 192 (atom-hash (make-hash-table :test #'eq :size 8))) 193 194 (setq l (margs l)) 195 196 ;; simplfy and flatten 197 (let (($%enumer $numer)) ;; convert %e --> 2.718...Why not %pi too? See simpcheck in simp.lisp. 198 (dolist (li l) 199 (setq li (simplifya li z)) 200 (if (mplusp li) (setq acc (append acc (cdr li))) (push li acc)))) 201 (setq l acc) 202 (setq acc nil) 203 (dolist (li l) 204 ;;(if (atom li) (incf *its-an-atom*) (incf *not-an-atom*)) 205 (cond ((mnump li) (mincf num-sum li)) 206 ;; factor out infrequent cases. 207 ((and (consp li) (consp (car li)) (memq (caar li) '(mequal mrat $matrix mlist $interval))) 208 (setq op (caar li)) 209 (cond ((eq op 'mequal) 210 (push li mequal-terms)) 211 (($taylorp li) 212 (push li taylor-terms)) 213 ((eq op 'mrat) 214 (push li mrat-terms)) 215 ((eq op '$matrix) 216 (push li matrix-terms)) 217 ((eq op '$interval) 218 (push li interval-terms)) 219 ((eq op 'mlist) 220 (if $listarith (push li mlist-terms) (push (convert-to-coeff-form li) acc))))) 221 222 ;; Put non-infinite atoms into a hashtable; push infinite atoms into inf-terms. 223 ((atom li) 224 (if (memq li '($minf $inf $infinity $und $ind)) 225 (push li inf-terms) 226 (progn 227 (setq cf (gethash li atom-hash)) 228 (setf (gethash li atom-hash) (if cf (1+ cf) 1))))) 229 230 (t (push (convert-to-coeff-form li) acc)))) 231 232 ;; push atoms in the hashtable into the accumulator acc; sort acc. 233 (maphash #'(lambda (cf a) (push (cons cf a) acc)) atom-hash) 234 (setq l (sort acc 'great :key 'car)) 235 236 ;; common term crunch: when the new coefficient is -1 or 1 (for example, 5*a - 4*a), 237 ;; set the "do-over" flag to true. In this case, the sum needs to be re-simplified. 238 ;; Without the do over flag, a + 5*a - 4*a --> a + a. Last I checked, the testsuite 239 ;; does not test the do-over scheme. 240 241 (setq acc nil) 242 (while l 243 (setq x (pop l)) 244 (setq cf (cdr x)) 245 (setq x (car x)) 246 (while (and l (like x (caar l))) 247 (mincf cf (cdr (pop l)))) 248 (if (and (or (eql cf 1) (eql cf -1)) (mplusp x)) (setq do-over t)) 249 (setq x (number-times-expr cf x)) 250 (cond ((mnump x) (mincf num-sum x)) 251 ((not (mzerop x)) (push x acc)))) 252 253 ;;(setq acc (sort acc '$orderlessp)) ;;<-- not sure this is needed. 254 255 ;; I think we want x + 0.0 --> x + 0.0, not x + 0.0 --> x. 256 ;; If float and bfloat were simplifying functions we could do 257 ;; x + 0.0 --> float(x) and 0.0b0 + x --> bfloat(x). Changing this 258 ;; test from mzerop to (eq 0 num-sum) causes problems with the test suite. 259 ;; For example, if x + 0.0 --> x + 0.0, we get an asksign for 260 ;; tlimit((x*atan(x))/(1+x),x,inf). That's due to the (bogus) floating point 261 ;; calculations done by the limit code. 262 263 ;;(if (not (eq 0 num-sum)) (push num-sum acc)) 264 (if (not (mzerop num-sum)) (push num-sum acc)) 265 266 ;;(if do-over (incf *do-over*)) ;; never happens for testsuite! 267 (setq acc 268 (cond (do-over (simplifya `((mplus) ,@acc) nil)) 269 ((null acc) num-sum) 270 ((null (cdr acc)) (car acc)) 271 (t (cons '(mplus simp) acc)))) 272 273 ;; special case dispatch 274 (if mequal-terms 275 (setq acc (add-expr-mequal acc mequal-terms))) 276 (if taylor-terms 277 (setq acc (add-expr-taylor acc taylor-terms))) 278 (if mrat-terms 279 (setq acc (add-expr-mrat acc mrat-terms))) 280 (if mlist-terms 281 (setq acc (add-expr-mlist acc mlist-terms))) 282 (if interval-terms 283 (setq acc (add-expr-interval acc interval-terms))) 284 (if matrix-terms 285 (setq acc (add-expr-matrix acc matrix-terms))) 286 (if inf-terms 287 (setq acc (add-expr-infinities acc inf-terms))) 288 289 acc)) 290 291 292