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