;;; changes to simp.lisp  to make it faster by RJF and Mark Hayden, Nov, 1991.
;;; revised by RJF on 10/14/03 to make it sort of work with sourceforge
;;; and again 7/17/06
;;; and again 7/22/06; not quite working.

;;; copied to www.cs.berkeley.edu/~fateman/lisp/newsimp.lisp  on 3/26/2008
;;; RJF

;;(in-package :climax)  ;; to read it into commercial macsyma
(in-package :maxima)
(eval-when (compile)
   (declaim (optimize (speed 3) (safety 1) (space 0) (compilation-speed 0))))
  ;;(use-package climax)

(define-modify-macro mincf (&optional (i 1)) addk)

(eval-when (compile load) (proclaim '(special simplifya $simp dosimp *opers-list 
		  substp $subnumsimp opers-list
		  $doallmxops $domxmxops $domxplus $listarith $float $%i
		  $doscmxplus $doscmxops check bigfloatzero bigfloatone
		  $sumhack $prodhack uniqtab $sqrtdispflag $%edispflag 
		  $vect_cross sign lind $exptdispflag
		  %e-val $ratfac genvar varlist *withinratf* modulus $norepeat
		  genpairs $keepfloat vlist $numer $numer_pbranch exptrlsw 
	    $factorflag le derivlist *fnewvarsw *errsw minus odds 
	    evens *great-hash* *great-etree* errorsw $dotdistrib
		  *mts* *mps* *mes*)))
;;bugs 

;;  expand is a no-op 11/23/91

;;  diff is slow.. need better uniquifying; maybe hashlists are too slow
;; for small lists or great is slow for common atom case?
;;    fixed ... not doing  simplification of 2*(sum) properly.
;;  still a bug, repeated diff(f[5],x,3) puts new things on uniqtab
;;  12/4/91
;;  there seems to be some program that is changing stuff that is
;; hashed.  saving a pointer to the result of hcons shows that after
;; it is set, someone is changing it.  rplaca/rplacd/nconc, when traced,
;; don't show anything.

;;  rat stuff may not work

;;  ev(foo,x=y)  doesn't work if foo is hashlist.

;;; ret2list, if nil, leaves the result of mplus simplification as 
;;; a hash table.  if t, it changes the hash table to a list each
;; time.   It is set to nil, which should be the 'fast' setting
;; but may break things.

;;(defvar ret2list nil)
(defvar ret2list t) ;; for now, easier to interface with rest of world

;;; (setf ret2list t) makes most everything work, I think.

;;; many programs from simp.lisp don't need to be changed.
;;; some of these probably should be changed to use ucons etc

;;; revised by RJF  Nov 13--18, 1991

;; changes:  
;; 
;; 1. hash-coding of mplus data form in plshash
;;
;; 2. introduce use of CL rational numbers.
;;  e.g. 1/2 instead of ((rat simp) 1 2). This affects many parts of
;; the code from the parser, through the simplifier (etc), to the display.
;;  We try to maintain the old way as a transition.
;; 

;; Regarding mplus as a hashtable:

"proposal a.
 instead of ((mplus simp) x y ((mtimes) 3 z)) we {sometimes}
 use ht=<hashtable> with
 (get-hash 'x ht) = 1
 (get-hash 'y ht) = 1
 (get-hash 'z ht) = 3

proposal b.
 instead of just ht  we use ((mplus simp) ht) or perhaps
 ((mplus simp ht) ht) {two pointers to ht if it is a convenience for old code}

proposal c.
 ((mplus simp ht) x y ((mtimes) 3 z))  .. both forms stored. Handy for
 computing $part, $inpart, output to disk, etc. They must both be up to
 date, or else one of them is missing:
 ((mplus simp) x y ((mtimes) 3 z))    ;; ht is missing
 ((mplus simp ht) ht)  ;; list is missing

This way if you repeatedly add stuff to the expression, you just keep
updating the hash-table.  When some program requires the list version,
you re-convert and save it. 

proposal d.

Use (only) a hashtable but insert into it one
element with key = *tabletype* and val = mplus

proposal e.

Use (only) a hashtable but insert into it one
element with key = *tabletype* and val = mplus,
and one element with key *listversion* and val = ((mplus simp ) ...)
Question: how do you know if the listversion is up to date?
 possibilities: every time you change the hashtable you clear the *listversion*
 value; or you store a time on the hashtable (last changed time) and
 on the value of listversion.  Check time-stamps. Bad news if you
 use rplacx. 

proposal f.

Use a hash-table only internal to simplus. On conversion out of a 
hash table, leave the hash table as the cadddr of the result.

***** THIS IS WHAT IS USED HERE if ret2list is set to true (default)
 but if ret2list is nil, proposal  a   ***
"

;;  Here are the programs...

;; Simplifya simplifies the expression X
;; y is t or nil depending on whether args are simplified or not resp.
(defun simplifya (x y)
  (cond
   ;;screwy inefficiency. if numer=true,let meval bind $%pi to a number. -rjf 11/91
   ;;((atom x) (cond ((and (eq x '$%pi) $numer) %pi-val) (t x)))
   ((atom x)
    (cond((hash-table-p x)
	  (if ret2list (hash2list+ x) x))
	 (t x)))
   ;;most atoms are already simplified.  check for hashtable
   ((not $simp) x)			;if simplify is turned off, do nothing
   ((atom (car x))			;an atomic car is not usual --check for error, mostly
    (cond ((and (cdr x) (atom (cdr x)))
	   (merror "~%~s is a cons with an atomic cdr - simplifya" x))
	  ((get (car x) 'lisp-no-simp)
	   ;; this feature is to be used with care. it is meant to be
	   ;; used to implement data objects with minimum of consing.
	   ;; forms must not bash the displa package. only new forms
	   ;; with carefully chosen names should use this feature.
	   x)
	  (t (cons (car x)		;simplify the arguments
		   (mapcar (function (lambda (x) (simplifya x y)))
			   (cdr x))))))
   ;; this is the most usual case, e.g. (($foo) bar)
   (t 
    (uniq (let* ((head (caar x))
		 (w (if (atom head)(get head 'operators) nil))
		 (headcdr (cdar x)));; if x=(($foo simp) bar), headcdr =(simp)
	    (cond 
	     ((and (not dosimp) (memq 'simp headcdr)) x);already simp'd
	     ((eq head 'mrat) x)	; canonical rational form -> leave alone
	     ((and w (not (memq 'array headcdr)))
	      ;; old code has rulechk to  see if simplifier has been
	      ;; modified by tellsimp etc.
	      ;; we will change tellsimp etc to remove operators property.
	      ;; (when?   -- rjf 11/13/91)
	      ;;	  (format t "funcall ~s~%" w)
	      (funcall w x 1 y))
	     ((not (atom head));;various peculiar and rare cases
	      (cond ((or (eq (car head) 'lambda)
			 (and (not (atom (car head)))
			      (eq(caar head) 'lambda)))
		     (mapply head (cdr x)head))
		    (t (merror "illegal form - simplifya:~%~s" x))))
	     ((get head 'opers);; simplify forms in various ways
	      (let ((opers-list *opers-list)) (oper-apply x y)))
	     ((and (eq head 'mqapply)
		   (or (atom (cadr x))
		       (and (eq substp 'mqapply)
			    (or (eq (car (cadr x)) 'lambda)
				(eq (caar (cadr x)) 'lambda)))))
	      (cond ((or (symbolp (cadr x)) (not (atom (cadr x))))
		     (simplifya (cons (cons (cadr x) headcdr) (cddr x)) y))
		    ((or (not (memq 'array headcdr)) (not $subnumsimp))
		     (merror "improper value in functional position:~%~m" x))
		    (t (cadr x))))
	     ;; if we know nothing about the function, just do this
	     (t (simpargs x y))))))))


(defun rulechk (x) (or (mget x 'oldrules) (get x 'rules)))

(defun resimplify (x) (let ((dosimp t)) (simplifya x nil)))

(defun ssimplifya (x) (let ((dosimp t)) (simplifya x nil)))  ; temporary

(defun simpargs (x y)
  (if (or (eq (get (caar x) 'dimension) 'dimension-infix)
	  (get (caar x) 'binary))
      (twoargcheck x))
  (if (and (memq 'array (cdar x)) (null (margs x)))
      (merror "subscripted variable found with no subscripts."))
  (eqtest (if y x (let ((flag (memq (caar x) '(mlist mequal))))
		    (cons (ncons (caar x))
			  (mapcar #'(lambda (u)
				     (if flag (simplifya u nil)
					      (simpcheck u nil)))
				  (cdr x)))))
	  x))


;; there are several special cases handled by simplus. (comments by rjf)

;; equations + scalars
;;  v +(a=b)  -> v+a = v+b
;; (a=b) + (c=d) =  a+c = d+b

;; matrices (or lists)  + scalars or matrices
;; v+matrix([a,b],[c,d]) -> matrix([a+v,b+v],[c+v,d+v])
;; matrix(..)+matrix(...)-> sum of matrices if sizes are the same.
;; v+ [a,b,c] -> [v+a,v+b,v+c]
;; if list or matrix size incompatibilities occur, error from "fullmap"

;; matrix (or lists)  + equation
;; (a=b)+matrix produces  (matrix+a)=(matrix+b)

;; mrat form (canonical rational form)
;; mrat(e) + anything -> rat (e+anything)

;;%sum form
;; sum(a[i],i,1,n)+sum(a[i],i,n+1,m) -> sum(a[i],i,1,m)
;; (this is done by sumpls, and the circumstances have to be just right
;; to make it happen)

;; various other combinations are also handled.

;; sum+(a=b) -> sum+a = sum+b
;; sum+matrix -> matrix of (entry[i,j]+sum).
;; mrat+sum -> mrat, with one "rational variable" being a sum.
;; mrat+matrix -> matrix of (entry[i,j]+ "disrep'd" mrat).

;;  this kind of treatment is really uncomfortable.  add a new
;; data type and you have to make n coercion decisions and/or
;;  break the program. 


;; examples: 
;; a+b*%i+c+d*%i -> %i*d+c+%i*b +a.
;; interval(a,b)+interval(a,b) -> 2*interval(a,b)
;; infinity+infinity-> 2*infinity;  (limit(2*infinity) -> infinity, though)

;; an apology for the structure of the program -- it even has goto and
;; labels!  it was originally written by knut korsvold, circa 1963
;; when lisp was 3 years old and people didn't know better. -- rjf

(defun simplus (x w z)  ; w must be 1
  (prog (res check eqnflag matrixflag sumflag (len 0) (terms nil) h)
	(if (endp (cdr x)) (return 0)) ;case of ((mplus) ) = 0
	(if (endp (cddr x)) (return (simplifya (cadr x) nil))) ; case of ((mplus) x)

	;; check is used to try to preserve existing structure. In
	;; case no changes are made by simplification, eqtest will
	;; return check rather than the new consed structure.
;	(setq check x)

	;;this is the start of a big loop over the ((mplus)...)
    ;;structure handed to simplus
      (if z nil (setf x (dolist (i (cdr x) (cons (car x) terms))
			(incf len)
			(setf h (simplifya i nil))
			(cond ((mplusp h) (dolist (k (cdr h))(push k terms)))
			      (t (push h terms))))))

	(setq res (make-hash-table :size 15 :test #'eq))

	;; for identification purposes in case we use hash tables for other
	;; things.  consider this:
	;;	(setf (gethash '*tabletype* res) 'mplus)
start

        (setq x (cdr x))  ;; x was ((mplus) a b c ..) now x=(a b c...)
	
	(if (endp x) (go end)) ;; no more entries

	;;set w to the first item in x  (e.g. a)
	(setq w (car x))
	;; examine w.
   st1  
   (cond
    ((atom w) (plshash w res)) ;;w could be a hash-table too
    (t (case 
	(caar w)
	((mtimes mexpt)   (plshash w res)) ;quick check for usual
	;; stick in ((mplus) stuff) but simplify it..
	;; won'thappen
;;	(mplus 
;;	 (format t "~% mplus of mplus ~s~%" (cdr w))
;;	 (mapc #'(lambda(z)(plshash z res)) (cdr w)))
	
	(mrat (cond ((or eqnflag matrixflag
		    (and sumflag (not (memq 'trunc (cdar w))))
		    (spsimpcases (cdr x) w))
		(setq w (ratdisrep w)) (go st1))
	       (t (return (ratf (cons '(mplus)
				      (nconc (mapcar #'simplify (cons w (cdr x)))
					     (cdr res))))))))
	(mequal
	 (setq eqnflag
	       (if (not eqnflag)
		   w
		   (list (car eqnflag)
			 (add2 (cadr eqnflag) (cadr w))
			 (add2 (caddr eqnflag) (caddr w)))))	)
	((mlist $matrix)
	 (setq matrixflag
	       (cond ((not matrixflag) w)
		     ((and (or $doallmxops $domxmxops $domxplus
			       (and (eq (caar w) 'mlist) ($listp matrixflag)))
			   (or (not (eq (caar w) 'mlist)) $listarith))
		      (addmx matrixflag w))
		     (t (plshash w res) matrixflag))))
	(%sum
	 (setq sumflag t res (sumpls w res))
	 (setq w (car res) res (cdr res)))
	(t (plshash w res))))) ; other cases are just added in

;; debug (format t "~% res = ~s~%" (hash2list+ res))
   (go start)
;; termination
;; re-think this.  res is a hashtable now
   end
 (setq res (normalhash res))
;; if ret2list is t, make result a list here.
   (if ret2list (setq res (hash2list+ res)) )
    (if matrixflag
	    (setq res (cond ((zerop1 res) matrixflag)
			    ((and (or ($listp matrixflag)
				      $doallmxops $doscmxplus $doscmxops)
				  (or (not ($listp matrixflag)) $listarith))
			     (mxplusc res matrixflag))
			    (t (testp (pls matrixflag (pls res nil)))))))
;;	(setq res (eqtest res check))
	(return (if eqnflag
		    (list (car eqnflag)
			  (add2 (cadr eqnflag) res)
			  (add2 (caddr eqnflag) res))
		    res))))

(defun *red1 (x) (*red (cadr x) (caddr x)))

(defun *red (n d)
  (if $float (/ n d 1.0)  (/ n d)))

(defun num1 (a) (if (rationalp a) (numerator a) a))

(defun denom1 (a) (if (rationalp a) (denominator a) 1))

;; big change to handling complex numbers !!$$!!  -- use common lisp
(setq $%i  (complex 0 1))

(defun timesk (x y)  ; x and y are assumed to be already reduced
 (cond ((equal x 1) y)
       ((equal y 1) x)
       ((and (numberp x) (numberp y)) (times x y))
       ((or ($bfloatp x) ($bfloatp y)) ($bfloat (list '(mtimes) x y)))
       ((floatp x) (times x (fpcofrat y)))
       ((floatp y) (times y (fpcofrat x)))
       (t (timeskl x y))))

(defun timeskl (x y)
 (prog (u v g)
       (setq u (/ (num1 x) (denom1 y)))
       (setq v (/ (num1 y) (denom1 x)))
       (setq g (cond ((or (equal u 0) (equal v 0)) 0)
		     ((equal v 1) u)
		     ((and (numberp u) (numberp v)) (* u v))
		     ))
       (return (cond ((numberp g) g)
		     ($float (fpcofrat g))
		     (t g)))))

;; old program was elaborate conversion of rational to float to avoid
;; overflows.

(defun fpcofrat (ratno) (/ (cadr ratno) (caddr ratno) 1.0))


(defun expta (x y) (cond ((equal y 1) x)
			 ((numberp x) (exptb x (num1 y)))
			 (($bfloatp x) ($bfloat (list '(mexpt) x y)))
			 ((minusp (num1 y))
			  (/ (exptb (caddr x) (minus (num1 y)))
				(exptb (cadr x) (minus (num1 y)))))
			 (t (/ (exptb (cadr x) (num1 y))
				  (exptb (caddr x) (num1 y))))))

(defun exptb (a b)
 (cond ((equal a %e-val) (exp b))
       ((or (floatp a) (not (minusp b))) (expt a b))
       (t (setq b (expt a (minus b))) (/ 1 b))))

(defmacro ulist(&rest l)(cond ((null l)nil)
			      (t `(ucons ,(car l) (ulist ,@(cdr l))))))

(defun hash2list+ (h)
  (let (arr (i 0) (size 0))
    (declare (fixnum i))
    (cond ((hash-table-p h)

	   ;; first remove zero entries  (optional, as written here)
;;	   (setq h (normalhash h))
	   ;;next, copy the hash-table into an array
	   (setq arr (make-array (setq size (hash-table-count h))))
	  ; (setq i 0)
	   (maphash #'(lambda(key val)
			;;each element in the array will be (key . val)
			(setf (svref arr i) (cons key val))
			(incf i)) h)

	   ;; now sort the elements of the array by keys (keys are distinct)

	   (sort arr #'(lambda(k1 k2)(great (car k1)(car k2))))

	   ;; now run through entries computing a simplified form.
	   ;; we are not unique-ifying the running sum or the products
	   ;; at the moment.  the keys will be unique, though.

;; we could we do this in place in the array and then coerce to a list.
	   (setq h
		 (do ((i 0 (1+ i))
		(key)
		(val)
		(res)) ; three locals, initialized to nil
		     ;; the termination test
		     ;; if we want to save the hash version, we can do
		     ;; something like (cons (list 'mplus 'simp h) res)  below
	       ((= i size) (ucons *mps* res))  

	       (setf key (car (svref arr i)))
	       (setf val (cdr (svref arr i)))

	       (cond
		((eq val 1)		; 1*key -> key
		 (setq res (ucons key res)))
		((eq val 0))		; 0*key -> 0
		((eq key 1)		;constant number *1 -> number
		  (setq res (ucons val res)))
		;; all the other cases produce a ((mtimes simp) <num> expr)
		((or (atom key)(not(eq(caar key) 'mtimes)))
		 (setq res (ucons (ulist *mts* val key) res)))
		;; case of ((mtimes simp) 3 ((mtimes simp) x y))
		;; --> ((mtimes simp) 3 x y)
		(t (setq res (ucons (ucons *mts* (ucons val (cdr key))) res)) 
		   ))))
	   ;; check for degenerate cases ((mplus)) --> 0, ((mplus) a)  --> a
	   (cond ((null (cdr h)) 0)
		 ((null (cddr h)) (cadr h))
		 (t h)))
	
	  (t h);; not a hashtable. just return it.
	  )))


;; changed to use cl  /

(defun simpquot (x y z)
 (twoargcheck x)
 (cond
       ((and (numberp (cadr x)) (numberp (caddr x)) (not (zerop (caddr x))))
	(/ (cadr x) (caddr x)))
       (t (setq y (simplifya (cadr x) z))
	  (setq x (simplifya (list '(mexpt) (caddr x) -1) z))
	  (if (equal y 1) x (simplifya (list '(mtimes) y x) t)))))


(defun plshash (x table)

;; there are 6 cases for x, which is already simplified.

;; x= atomic name e.g.. alpha

;; x= ((mtimes)  <number>  <otherstuff>)  e.g. ((mtimes) 3 ((mexpt) y 3))

;; x= ((mtimes) <otherstuff>)  e.g. ((mtimes) alpha beta)

;; x= ((mplus) ....)  ;; won't happen.

;; x= ((<otherop>) ...)  e.g. ((mexpt) y 3)  or (($foo) bar)

;; x= hashtable of type mplus

;; the first and 5th are treated the same.

  (cond((atom x)
	(typecase x
		  (number
		   ;; add all numbers into hashtable value for key=1
		   (let((v (+ x (gethash 1 table 0))))
		     (if (eql v 0) (remhash 1 table)
		       (setf (gethash 1 table) v))))
		  (hash-table
		   ;; maybe another plus hash-table !
		  (plshashhash x table) ) ;merge them
		  ;; add all other atoms into their own buckets
		  (t  
		   (if (mzerop (mincf (gethash x table 0)))(remhash x table)))))
       ;;at this point x is not atomic, so it has a caar.
       ;; (caar x) = mplus may have been taken care of by simplus? check anyway
     ;;  ((eq (caar x) 'mplus) (mapc #'(lambda(z)(plshash z table))(cdr x)))
       ((not(eq (caar x) 'mtimes)) 
	;; add using macsyma's "addk" to existing entry: a rational or integer.
	(if (mzerop (mincf (gethash x table 0)))(remhash x table)))

       ;; the only other case of interest is mtimes.
       (t				; (caar x) = 'mtimes
	(let ((s (cadr x)))		;s is the 2nd item in the 
					;list ((mtimes ) s ...)
	  (if (mnump s) 
	      (let ((r 
		     (if (endp (cdddr x))
				;; this is the case of  ((mtimes) 3 u)
				;; the key is u, the value is 3
			 (caddr x)
				;; this is the case of ((mtimes) 3 u v)
				;; the key is ((mtimes) u v) the value is 3
		       (ucons *mts* (cddr x))) ))
	  (if (mzerop (mincf  (gethash r table 0) s)) (remhash r table)))
	    ;; this is the case of x = ((mtimes) u v ..) [s not number]
	     (if (mzerop (mincf  (gethash s table 0) 1)) (remhash s table)))))))

(defun plshashhash(x table)
  ;; x is a hash-table too
  (maphash #'(lambda(key val)
	       (mzerop (mincf  (gethash key table 0) val)))
	   x)
  table) ;;clobbers table..

;; this tests for (), ((mplus) 0 x) ((mplus) 0 x y))
(defun testp (x) (cond ((atom x) 0)
		       ((null (cddr x)) (cadr x))
		       ((zerop1 (cadr x))
			(cond ((null (cdddr x)) (caddr x)) (t (rplacd x (cddr x)))))
		       (t x)))


;; in case alike1 is still used somewhere,
(defun alike1 (x y) (eq  x y))


;; maps alike1 down two lists.

(defun alike (x y)
;; maybe should worry about hash lists too.?
 (do ((x x (cdr x)) (y y (cdr y))) ((atom x) (equal x y))
     (cond ((or (atom y) (not (alike1 (car x) (car y))))
	    (return nil)))))

;; for comparison
(defun old-great (x y)
 ;;(declare (object y))
 (cond ((atom x)
	(cond ((atom y)
	       (cond 
		     ((numberp x)
		      (cond ((numberp y)
			     (setq y (- x y))
			     (cond ((zerop y) (floatp x)) (t (plusp y))))))
		     ((constant x)
		      (cond ((constant y) (alphalessp y x)) (t (numberp y))))
		     ((mget x '$scalar)
		      (cond ((mget y '$scalar) (alphalessp y x)) (t 
								  (maxima-constantp y)
								  )))
		     ((mget x '$mainvar)
		      (cond ((mget y '$mainvar) (alphalessp y x)) (t t)))
		     (t (or (maxima-constantp y) (mget y '$scalar) 
			    (and (not (mget y '$mainvar)) (alphalessp y x))))))
	      (t (not (ordfna y x)))))
       ((atom y) (ordfna x y))
       ((eq (caar x) 'rat)
	(cond ((eq (caar y) 'rat)
	       (greaterp (times (caddr y) (cadr x)) (times (caddr x) (cadr y))))))
       ((eq (caar y) 'rat))
       ((memq (caar x) '(mbox mlabox)) (old-great (cadr x) y))
       ((memq (caar y) '(mbox mlabox)) (old-great x (cadr y)))
       ((or (memq (caar x) '(mtimes mplus mexpt %del))
	    (memq (caar y) '(mtimes mplus mexpt %del)))
	(ordfn x y))
       ((and (eq (caar x) 'bigfloat) (eq (caar y) 'bigfloat)) (mgrp x y))
       (t (do ((x1 (margs x) (cdr x1)) (y1 (margs y) (cdr y1))) (())
	      (cond ((null x1)
		     (return (cond (y1 nil)
				   ((not (alike1 (mop x) (mop y)))
				    (old-great (mop x) (mop y)))
				   ((memq 'array (cdar x)) t))))
		    ((null y1) (return t))
		    ((not (alike1 (car x1) (car y1)))
		     (return (old-great (car x1) (car y1)))))))))

;; a faster version of  test(n):=sum(a[i],i,1,n);

(defun $test(n)(let ((tab
		      (make-hash-table :size n :test #'eq)))
		 
		(do ((i 1 (1+ i)))
		    ((= i n) (hash2list+ tab))
		    (plshash (ulist '($a array) i) tab))))

;; from asum.lisp

(defun assume (pat &aux (ret2list t)) ;; make simplus return list
  (if (and (not (atom pat))
	   (eq (caar pat) 'mnot)
	   (eq (caaadr pat) '$equal))
      (setq pat `(($notequal) ,@(cdadr pat))))
  (let ((dummy (let (patevalled $assume_pos) 
	             (mevalp1 pat)))) ;; 2nd arg to mevalp1?
       (cond ((eq dummy t) '$redundant)
	     ((null dummy) '$inconsistent)
	     ((atom dummy) '$meaningless)
	     (t (learn pat t)))))


(defun dosum (exp ind low hi sump)
  (if (not (symbolp ind))
      (merror "improper index to ~:m:~%~m" (if sump '$sum '$product) ind))
  (unwind-protect 
   (prog (*i u lind l*i *hl (ret2list t))
	 
	 #+ignore (when (if sump (not $sumhack) (not $prodhack))
	       (assume (list '(mgeqp) ind low))
	       (if (not (eq hi '$inf))
		   (assume (list '(mgeqp) hi ind))))
	 (setq lind (cons ind nil))
	 (cond ((not (fixnump (setq *hl (m- hi low))))
		(setq exp (mevalsumarg exp ind low hi sump))
		(return (cons (if sump '(%sum) '(%product))
			      (list exp ind low hi))))
	       ((equal *hl -1) (return (if sump 0 1)))
	       ((signp l *hl)
		(cond ((and sump $sumhack)
		       (return
			(m- (dosum exp ind (m+ 1 hi) (m- low 1) t))))
		      ((and (not sump) $prodhack)
		       (return
			(m// (dosum exp ind (m+ 1 hi) (m- low 1) nil)))))
		(merror "lower bound to ~:m: ~m~%is greater than the upper bound: ~m"
			(if sump '$sum '$product) low hi)))
	 (setq l*i (list low))

	 (setq u (if sump
			  ;; empty hashtable if sum
		     (make-hash-table :test #'eq :size (max 8 
							    (truncate
							     (+ 1 hi (- low)) 0.7)))
		   ;; aim at 70%
			1));; 1 if product
	 (setq *i low)
	 (loop ;;until return
;(format t "*i= ~s~%" *i)
	  (setq u
		   (if sump
		       (plshash 
			 (mbinding (lind l*i) (meval exp))
			u)
		     ;;regular handling for times
		     (simplifya (list '(mtimes)
				     (mbinding (lind l*i) (meval exp))u)
			       t)))
	  ;;termination
	     (if (equal *i hi) (return u))
	     (setq *i  (car (rplaca l*i (if (numberp *i)(1+ *i)(m+ *i 1)))))
    	     ); end loop

	 #+ignore
	 (cond ((if sump (not $sumhack) (not $prodhack))
		(let ((ret2list t))
;;;	   (format t "forgetting ~s < ~s < ~s~%" low ind hi)
		  (forget (list '(mgeqp) ind low))
		  (if (not (eq hi '$inf)) (forget (list '(mgeqp) hi ind))))))
	 (return u))))







;;;
(defun dispuniq()(maphash #'(lambda(key val)(displa (cons key val))) uniqtab))
(defun listuniq()(maphash #'(lambda(key val)(format t "~s~%"
						    (cons key val))) uniqtab))


;;**********************************************************************;;

;;**********************************************************************;;
;; shortcut to greatness

(defun great (x y)
 (cond ((atom x)
	(cond ((atom y)
	       (if (numberp x)
		   (if (numberp y) 
		       (if (or(complexp x)(complex y))
			   (complex-great x y)
			 (> x y))
		     nil) ;; order numbers by value; complexes by real then im.
		 ;; if x is a number but y isn't, return nil.
		 (fastgreat x y) ) ;; x and y are atoms but not numbers
	       )
	      (t (fastgreat x y))))
       ;; x is not atom
       ((numberp y) t)
       ;; neither x nor y is a number (or atom).
       (t  (fastgreat x y))))

(defun complex-great (x y)(cond
			   ((> (realpart x)(realpart y)) t)
			   ((eql (realpart x)(realpart y))
			    (> (imagpart x)(imagpart y)))
			   ;; either equal or <. 
			   (t nil)))

(defun fastgreat (a b)
  (> (the fixnum (etree-si a #'old-great-comp))
     (the fixnum (etree-si b #'old-great-comp))))
    

;; new tree routine wants cmp similiar to 'strcmp'

(defun old-great-comp (a b)
 (let ((res
	(cond ((eql a b) 0)
	((old-great a b) 1)
	(t -1))))
   ;set temporary position for "a" in hashtable
   (setf (gethash a *great-hash*)(+ res (gethash b *great-hash*)))
   res))

;;; more fixes from rjf
;; from nforma.lisp
;; this allows for complex numbers as well as rational numbers in
;; common lisp to be displayed in macsyma/ maxima/   a better
;; way might be to fix dim-atom etc and avoid consing up "rat" stuff
;; here.  -- rjf 14 nov. 1991

;;; also allows for display of mplus encoded in hash table

(defun nformat (form)
  (cond ((atom form)
	 (typecase form
	   
	   ;;
	   ;;	   ((integer * -1) ;negative integers
	   ;;		    (list '(mminus) (minus form)))
	   ;;		   (integer form) ;positive integers unchanged
	   ;; GCL doesn't understand the above.
	   (integer
	    (if (< form 0) `((mminus) ,(- form)) form))
	   
	   
		   ;; uncertain what to do here...
		   (complex 
		    (let ((r(realpart form))
			  (i(imagpart form)))
		      (setq i 
			    (if (equalp i 1) '$%i
			      `((mtimes) ,i $%i)))
		      (if (equalp r 0) (nformat i)
			(nformat`((mplus)  ,r ,i)))))
		   ((or (double-float * (0.0d0)) (single-float * (0.0e0))) ;neg
			      (list '(mminus) (minus form)))
		   ((or double-float single-float) form)
		   ((rational 0 *) ;positive ratios
	            (list '(rat) (numerator form)(denominator form)))
		   (ratio  ;negative ratios
		      (list '(mminus) (nformat(minus form))))
;;; if we use hash tables for mplus, un-do it here
		   (hash-table (nformat(hash2list+ form)))
		   (t form) ;a symbol, string, etc?
			    
		   ))
;; i don't understand this next line, so i'm commenting it out .. it
;; breaks the display of #c(3 -4) . it seems to do something about 
;; forms without simp flags not being formatted  --rjf
;;	((null (cdar form)) form)
	((atom (car form))
	 form)
	((atom (caar form))
	 (case (caar form)
	       ;; maybe obsolete if we get rid of all non-cl ((rat) a b)
	       (rat
		(cond ((minusp (cadr form))
		       (list '(mminus) (list '(rat) (minus (cadr form)) (caddr form))))
		      (t  form  ; was.. (cons '(rat) (cdr form))
)))
	       (mmacroexpanded (nformat (caddr form)))
	       (mplus (form-mplus form))
	       (mtimes (form-mtimes form))
	       (mexpt (form-mexpt form))
	       (mrat (form-mrat form))
	       (mpois (nformat ($outofpois form)))
	       (bigfloat
		(if (minusp (cadr form))
		    (list '(mminus) (list (car form) (minus (cadr form)) (caddr form)))
		    (cons (car form) (cdr form))))
	       (t form)))))

(defun form-mexpt (form &aux exp)
  (cond ((and $sqrtdispflag (eql 1/2 (caddr form)))
	 (list '(%sqrt) (cadr form)))
	((and $sqrtdispflag (eql -1/2 (caddr form)))
	 (list '(mquotient) 1 (list '(%sqrt) (cadr form))))
	((and (or (and $%edispflag (eq '$%e (cadr form)))
		  (and $exptdispflag (not (eq '$%e (cadr form)))))
	      (not (atom (setq exp (nformat (caddr form)))))
	      (eq 'mminus (caar exp)))
	 (list '(mquotient) 1 (if (equal 1 (cadr exp)) (cadr form)
				  (list '(mexpt) (cadr form) (cadr exp)))))
	(t (cons '(mexpt) (cdr form)))))

;; from compar.lisp

(defun sign (x)
  (cond ((mnump x) (setq sign (rgrp x 0) 
			 minus nil 
			 odds nil 
			 evens nil))
	((atom x) 
	 (if (hash-table-p x) 
	     (sign (hash2list+ x))
	   (if (eq x '$%i) 
	       (imag-err x)
	   (sign-any x)))) ;; check this 10/03/RJF
	((eq (caar x) 'mtimes) (sign-mtimes x))
	((eq (caar x) 'mplus) (sign-mplus x))
	((eq (caar x) 'mexpt) (sign-mexpt x))
	((eq (caar x) '%log) (compare (cadr x) 1))
	((eq (caar x) 'mabs) (sign-mabs x))
	((memq (caar x) '(%csc %csch))
	 (sign (inv* (cons (ncons (zl-get (caar x) 'recip)) (cdr x)))))
	((specrepp x) (sign (specdisrep x)))
	((kindp (caar x) '$posfun) (sign-posfun x))
	((or (memq (caar x) '(%signum %erf))
	     (and (kindp (caar x) '$oddfun) (kindp (caar x) '$increasing)))
	 (sign-oddinc x))
	(t (sign-any x))))

(defun constp (x)
   (cond ((floatp x) 'float)
	 ((numberp x) 'numer)
	 ((symbolp x) (if (memq x '($%pi $%e $%phi $%gamma)) 'symbol))
	 ((atom x) 
		       ;;fix up for mplus being a hash table 
		       ;; really inefficient in the worst case  -- rjf 11/22/91
	 (if (typep x 'hash-table)  (constp (hash2list+ x)) 'symbol))
	 ((eq (caar x) 'rat) 'numer)
	 ((eq (caar x) 'bigfloat) 'bigfloat)


	 
	 ((specrepp x) (constp (specdisrep x)))
	 (t (do ((l (cdr x) (cdr l))
		 (dum)
		 (ans 'numer))
		((null l) ans)
	      (setq dum (constp (car l)))
	      (cond ((eq dum 'float) (return 'float))
		      ((eq dum 'numer))
		      ((eq dum 'bigfloat) (setq ans 'bigfloat))
		      ((eq dum 'symbol)
		       (if (eq ans 'numer) (setq ans 'symbol)))
		      (t (return nil)))))))

;;; from simp

;;; this is just a minor hack so that (a+b)-(a+b) works fast
;;; we need to fix simptimes and friends so that a number times
;;; a hash-table is fast, and more generally, an expression
;;; can be multiplied by a hashtable, as well as two hash-tables
;;; by each other.  this latter operation is really only done
;;; by expand, though.

(defun simpmin (x vestigial z)
  (declare (ignore z vestigial))
 (oneargcheck x)
 (let ((z (cadr x)))
 (if (atom z)
     (cond ((numberp z) (minus z))
	   ((hash-table-p z)(minushash z))
	   (t (list *mts* -1 z)))
   (simplifya (list '(mtimes) -1 (simplifya (cadr x) z)) t))))

;; make a new hashtable and change all the signs

(defun minushash (z)
 (let ((res (make-hash-table :size (hash-table-count z) :test #'eq)))
  (maphash #'(lambda(key val) (setf (gethash key res)(- val))) z)
res))

;;; remove zero entries from the hash table.  not usually needed.

(defun hashcleanup(z)
  (maphash #'(lambda(key value)
	       (if (= value 0) (remhash key z))) z)z)


;; if the hash table has only 0 or 1 elements it is either zero, a constant
;; or a monomial that is not a sum.
;; normal changes such tables.

(defun normalhash(z)
 (if (and (hash-table-p z)(< (hash-table-count z)2))
     (hash2list+ z) ;; zero or number
   z))

;; from comm.lisp

(defvar *diff-hash* nil) ; it is nil if not within deriv
(defvar *diff-hash-var* nil)

(defun deriv (e)
       (prog (exp z count *diff-hash* )
	     (declare (special *diff-hash*))

	     (cond ((null e) (wna-err '$diff))
		   ((null (cdr e)) (return (stotaldiff (car e))))
		   ((null (cddr e)) (nconc e '(1))))
	     (setq exp (car e) z (setq e (copy-top-level e)))
	loop (if (or (null derivlist) (zl-member (cadr z) derivlist)) (go doit))
					; derivlist is set by $ev
	     (setq z (cdr z))
	loop2(cond ((cdr z) (go loop))
		   ((null (cdr e)) (return exp))
		   (t (go noun)))
	doit (cond ((nonvarcheck (cadr z) '$diff))
		   ((null (cddr z)) (wna-err '$diff))
		   ((not (eq (ml-typep (caddr z)) 'fixnum)) (go noun))
		   ((minusp (setq count (caddr z)))
		    (merror "improper count to diff:~%~m" count)))
	(cond(ret2list nil) ;;we don't do hashing if ret2list = t
	     ;;check if it is already set up
	     ((and *diff-hash* (eq *diff-hash-var* (cadr z))))
	     (t
;	      (format t "set up diff-hash-var = ~s~%" (cadr z))
	      ;; set it up
	      (setq *diff-hash-var* (cadr z))
	      (setq *diff-hash* (make-hash-table :test #'eq))))
	loop1(cond ((zerop count) (rplacd z (cdddr z)) (go loop2))
		   ((equal (setq exp (sdiff exp (cadr z))) 0) (return 0)))
	     (setq count (f1- count))
	     (go loop1)
	noun (return (diff%deriv (cons exp (cdr e))))))

;;;  use hash representation for results of sdiff

(defun sdiffhash(e x)
 (let ((res (make-hash-table :size (hash-table-count e) :test #'eq)))
   (maphash #'(lambda(key val)
		(setq key 
		      (hash2list+
		       (if (eql val 1)(sdiff key x)
			 (sdiff (m* val key) x))))
		(unless (eql key 0) (plshash key  res)))
	    e)
(setf *debughash*    (normalhash res))))

;;; ordinarily you'd think this could be done much more compactly
;;; and modularly by a data-directed approach.  maybe this is so, but
;;; it wasn't written that way to start, and there are niceties of this
;;; program that would be difficult to change back.  so we leave it,
;;; mostly.

(defun sdiff (e x)  ; the args to sdiff are assumed to be simplified.
	 (declare (special *diff-hash* *diff-hash-var*))

  (cond ((alike1 e x) 1)
	((mnump e) 0)
	((symbolp e)(if (depends e x) (sdiffgrad e x) 0))
	;; check to see if we have a hash table for this variable
	;; and return the answer if we've already computed it.
	((and *diff-hash* (eq *diff-hash-var* x)(gethash e *diff-hash*)))
	;; check if it is a sum in hash-table form.
	((hash-table-p e) (sdiffhash e x))
	(t
 (case (caar e)
   (mtimes(sdifftimeshash (cdr e) x))
   (mrat (ratdx e x))
   (mplus (addnhash (sdiffmap (cdr e) x) t))
   (mexpt
     (let ((res(diffexpt e x)))
	     (cond ((and *diff-hash* (eq *diff-hash-var* x))
		    (setf (gethash e *diff-hash*)res)))
	     res))
#+ignore
((%log %plog)
    ;; this seems to change diff(log(abs(u)),x) to diff(log(u),x)
    (sdiffgrad (cond ((and (not (atom (cadr e))) (eq (caaadr e) 'mabs))
		      (cons (car e) (cdadr e)))
		     (t e))
		    x))
   ((%sum %product) (diffsumprod e x))
   (mnctimes
    (let (($dotdistrib t))
	 (add2 (ncmuln (cons (sdiff (cadr e) x) (cddr e)) t)
	       (ncmul2 (cadr e) (sdiff (cons '(mnctimes) (cddr e)) x)))))
   (mncexpt (diffncexpt e x))
   (%derivative
	 (cond ((or (atom (cadr e)) (memq 'array (cdaadr e))) (chainrule e x))
	       ((freel (cddr e) x) (diff%deriv (cons (sdiff (cadr e) x) (cddr e))))
	       (t (diff%deriv (list e x 1)))))
   ((%binomial $beta)
	 (let ((efact ($makefact e)))
	   (mul2 (factor (sdiff efact x)) (div e efact))))
   (%integrate (diffint e x))
   (%laplace (difflaplace e x))
   (%at (diff-%at e x))
   ((%realpart %imagpart)
    (list (cons (caar e) nil) (sdiff (cadr e) x)))
 ;; we could try diff of "if" (mcond), etc. etc. etc. or
 ;; try it via sdiffgrad.
   (t(cond
      ((mbagp e) (cons (car e) (sdiffmap (cdr e) x)))
      ((and $vect_cross (eq (caar e) '|$~|))
       (add2* `((|$~|) ,(cadr e) ,(sdiff (caddr e) x))
	      `((|$~|) ,(sdiff (cadr e) x) ,(caddr e))))
      
;checked in sdiffgrad also?
;      ((not (depends e x)) 0)
      (t (let ((res(sdiffgrad e x)))
	     (cond ((and *diff-hash* (eq *diff-hash-var* x))
		    ;; should we unhash res?
		    ;;(setq res (hash2list+ res))
;		   (format t "remember diff of ~s by ~s~%" e x)
		    (setf (gethash e *diff-hash*)res)))
	     res))))))))


(defun sdiffgrad (e x)
 (let* ((fun (caar e)) (grad (oldget fun 'grad)) args)
      (cond ((and (eq fun 'mqapply) (oldget (caaadr e) 'grad))
	     (sdiffgrad (cons (cons (caaadr e) nil) (append (cdadr e) (cddr e)))
			x))
	    ((or (eq fun 'mqapply) (null grad))
	     (if (not (depends e x)) 0 (diff%deriv (list e x 1))))
	    ((not (= (length (cdr e)) (length (car grad))))
	     (merror "wrong number of arguments for ~:m" fun))
	    (t (setq args (sdiffmap (cdr e) x))
	       (uniq(addn (mapcar
		      #'mul2
		      (cdr (substitutel
			    (cdr e) (car grad)
			    (do ((l1 (cdr grad) (cdr l1))
				 (args args (cdr args)) (l2))
				((null l1) (cons '(mlist) (nreverse l2)))
				(setq l2 (cons (cond ((equal (car args) 0) 0)
						     (t (car l1)))
					       l2)))))
		      args)
		     t))))))


(defun sdifftimeshash (l x)
       (prog (term left (out(make-hash-table :test #'eq) ))
	loop (setq term (car l) l (cdr l))
	     (plshash (muln (cons (sdiff term x) (append left l)) t) out)
	     (if (null l) (return out))
	     (setq left (ucons term left))
	     (go loop)))

(defun addnhash(l f)
 (setq f (make-hash-table :test #'eq))
 (mapc #'(lambda(r)(plshash r f)) l)
 f)

;; debugging help
(defun $u() (hash-table-count uniqtab))

;;; more changes to programs in simp.lisp
;;
;; i think that we know  r1 and r2 are  constant.

(defun exptrl (r1 r2)
  (cond ((equal r2 1) r1)		; r1^1 = r1
	((equal r2 1.0) (cond ((mnump r1) (addk 0.0 r1)) (t r1)))
					;r1^1.0 ->float(r1)
	((equal r2 bigfloatone) (cond ((mnump r1) ($bfloat r1)) (t r1)))
					;r1^1.0b0-> bfloat(r1)
	((zerop1 r1)
	 (cond ((or (zerop1 r2) (mnegp r2))
		(cond ((not errorsw)
		       (merror "~m has been generated" (list '(mexpt) r1 r2)))
		      (t (throw 'errorsw t))))
					; 0^0 -> error if errorsw= nil
	       (t (zerores r1 r2))))    ; 0^0 -> 0     if errorsw=t  (hmm..)
	((or (zerop1 r2) (onep1 r1))
	 (cond ((or ($bfloatp r1) ($bfloatp r2)) bigfloatone) ; 1^x or x^0 -> 1
	       ((or (floatp r1) (floatp r2)) 1.0)
	       (t 1)))
	((or ($bfloatp r1) ($bfloatp r2)) ($bfloat (list '(mexpt) r1 r2)))
					;bigfloat^bigfloat -> compute it
	((and (numberp r1) (integerp r2)) (expt r1 r2)) ;use cl
	((and (numberp r1) (floatp r2) (equalp r2 (round r2)))
	 (expt  r1 (round r2)))
					;number ^ (integer or float)  cases
	((or $numer (and (floatp r2) (or (plusp (num1 r1)) $numer_pbranch)))
					;(positive ^float)  or $numer=t
	 (expt r1 r2)
	 ;; was a big mess in here
	 )

	((and (floatp r1) (rationalp r2))
	 (expt r1 r2))
	((and (typep r1 'ratio)
	      (typep r2 'ratio)
	      (< -1 r2 1) )
	      ;(n/d)^r2 -> (n^r2)/(d^r2)
	 (m* (exptrl (numerator r1) r2)(exptrl (denominator r1) (- r2))))
	((complexp r1)(if (floatp(realpart r1))(expt r1 r2)
			(ulist *mes* r1 r2)))
	((and (minusp r1) 
	      (typep r2 'ratio)
	      (< -1 r2 1))
	      ;(-n)^r2 -> (-1^r2)*(n^r2)  ;; questionable sometimes
	 (if (and(eql -1 r1)(eql r2 1/2)) '#c(0 1)  
	     (m* (exptrl -1 r2)(exptrl(- r1)  r2))))
	((and (integerp r1) ;; now it's positive
	      (typep r2 'ratio)
	      (< -1 r2 1))
	 ;;(k* n^m)^r2 -> k^r2*n^(m*r2)  , approximately
	 (exroot r1 r2))

	(exptrlsw (list *mes* r1 r2))
	(t
	 ;; we need to compute
	 ;;2^(5/2) -> 4* 2^(1/2)
	 ;; (9/4)^(1/2) -> 3/2   hmmmm.
	 ;; let a = 2^5*3^7*5^2.
	 ;; do we want a^(1/4) to be 6*sqrt(5)*54^(1/4) ?
	 (let ((exptrlsw t) res)
	   (multiple-value-bind 
	    (intpart fracpart)
	    (truncate r2)
	    ;; now r2 = intpart + fracpart, where 0< fracpart <1
	    (setq res (expt r1 intpart)) ;; get that out of way
	    (if (eql (denominator r1) 1)(m* res (exroot r1 fracpart))
	      (m* res (exroot (numerator r1) fracpart)
		  (power* (exroot(denominator r1) fracpart)-1))))))))

;; compute  v^r for v an integer, r rational -see if we can extract
;; a perfect root.
;; this is not as powerful as simpnrt which will get sqrt(50)-> 5 sqrt(2)
;; but it is heaps faster.  what to do, what to do...
;; macsyma philosophy would be to put in an <expletive deleted> flag.
;; -- rjf  12/8/91

#+ignore (defun exroot (v fracpart)
  (multiple-value-bind (root rem)
		       (irootv v (denominator fracpart))
		       (if (eql rem 0) (expt root (numerator fracpart))
			 (list *mes* v fracpart))))


;; changes  (27)^(1/2) to  3*sqrt(3)
;;  does not do  (25)^(1/4) to sqrt(5).
;;  macsyma does this too, though.

(defun exroot(v fracpart)
  (let* ((d (denominator fracpart))
	 (n (isnrt v d)))
    (cond ((eql n 1)(list *mes* v fracpart))
	  ((eql (expt n d) v)(expt n (numerator fracpart)))
	  (t (list '(mtimes) (expt n (numerator fracpart))
		   (list *mes* (/ v (expt n d)) fracpart))))))
	   
(defun isnrt(k r) ;find largest n such that n^r*m=k 
  ;; returns n
  (let(($factorflag t) (f  1))
    (do ((i (pairem (cfactor k)) (cdr i)))
	((null i)  f)
	(if (>= (cdar i) r)
	    (setq f (* f (expt (caar i)(truncate (cdar i) r)))))
	  )))

;;rearrange (a b c d)  to ((a.b) (c.d))

(defun pairem (m)(cond ((null m) nil)
		       (t(cons (cons (car m)(cadr m)) (pairem (cddr m))))))


;; same as iroot, but returns multiple values

(defun irootv (a n) ; computes a^(1/n)  see fitch, sigsam bull nov 74
  (cond ((f< (haulong a) n) (values 1 (sub1 a)))
	(t  ;assumes integer a>0 n>=2
	 (do ((x (expt 2 (f1+ (quotient (haulong a) n)))
		 (difference x (quotient (plus n1 bk) n)))
	      (n1 (f1- n)) (xn) (bk))
	     (nil)
	   (cond ((signp le (setq bk (*dif x (*quo a (setq xn (expt x n1))))))
		  (return (values x (difference a (times x xn))))))))))


(defun nrthk2 (in *n) (power* in (/ 1 *n)))  ; computes in^(1/*n)

(defun simpsqrt (x vestigial z)
   (declare (ignore vestigial)) ;ignored.       
 (oneargcheck x)
 ;; changed to c ^ 1/2
 (simplifya (list '(mexpt) (cadr x) 1/2) z))

;; from comm.lisp

(defun $float (e) 
  (cond ((numberp e)
	 (coerce e (if (complexp e)'(complex double-float) 'double-float)))
	((or (atom e) (memq 'array (cdar e))) e)
;	((eq (caar e) 'rat) (fpcofrat e))
	((eq (caar e) 'bigfloat) (fp2flo e))
	((memq (caar e) '(mexpt mncexpt))
	 (list (ncons (caar e)) ($float (cadr e)) (caddr e)))
	(t (recur-apply #'$float e))))


;;from rat3e.lisp
(defvar vlist nil)

(defun newvar1 (x) 
       (cond ((numberp x) nil)
	     ((memalike x varlist) nil)
	     ((memalike x vlist) nil)
	     ((atom x)
	      (cond ((symbolp x)(putonvlist x))
		    ((hash-table-p x)
		     (maphash #'(lambda(key val)(newvar1 key)) x))
		    (t ;what else could it be?
		  (putonvlist x))))
	     ((memq (caar x)
		    '(mplus mtimes rat mdifference
			    mquotient mminus bigfloat))
	      (mapc (function newvar1) (cdr x)))
	     ((eq (caar x) 'mexpt)
	      (newvarmexpt x (caddr x) nil))
	     ((eq (caar x) 'mrat)
	      (and *withinratf* (memq 'trunc (cdddar x)) (throw 'ratf '%%))
	      (cond ($ratfac (mapc 'newvar3 (caddar x)))
		    (t (mapc (function newvar1) (reverse (caddar x))))))
	     (t (cond (*fnewvarsw (setq x (littlefr1 x))
				  (mapc (function newvar1)
					(cdr x))
				  (or (memalike x vlist)
				      (memalike x varlist)
				      (putonvlist x)))
		      (t (putonvlist x))))))

(defun prep1 (x &aux temp)
       (cond ((floatp x)
	      (cond ($keepfloat (cons x 1.0)) ((prepfloat x))))
	     ((integerp x) (cons (cmod x) 1))
	     #+cl ((rationalp x)
		      (cond ((null modulus)(cons  (numerator x) (denominator x)))
			    (t (cquotient (numerator x) (denominator x)))))
	     #+cl ((integerp x) (cons (cmod x) 1))
	     #+cl ((complexp x)(mformat t "cre (rational form)
 of the complex number ~m is unlikely to work.~%
using its real part.~%" x)
		   (realpart x))
	     
	     ((atom x)
	      (cond((hash-table-p x)
		    (let ((res (cons 0 1) ))
		   
		      (maphash #'(lambda (key val)
				   (setq res
					 (ratplus res
						  (rattimes (prep1 key)(prep1 val) t)))) x)
		      res))
		   
		  ((assolike x genpairs))
		   (t (newsym x))))
	     ((and $ratfac (assolike x genpairs)))
	     ((eq (caar x) 'mplus)
	      (cond ($ratfac
		     (setq x (mapcar 'prep1 (cdr x)))
		     (cond ((andmapc 'frpoly? x)
			    (cons (mfacpplus (mapl #'(lambda (x)
						      (rplaca x (caar x)))
						  x)) 
				  1))
			   (t (do ((a (car x) (facrplus a (car l)))
				   (l (cdr x) (cdr l)))
				  ((null l) a)))))
		    (t (do ((a (prep1 (cadr x)) (ratplus a (prep1 (car l))))
			    (l (cddr x) (cdr l)))
			   ((null l) a)))))
	     ((eq (caar x) 'mtimes)
	      (do ((a (savefactors (prep1 (cadr x)))
		      (rattimes a (savefactors (prep1 (car l))) sw))
		   (l (cddr x) (cdr l))
		   (sw (not (and $norepeat (memq 'ratsimp (cdar x))))))
		  ((null l) a)))
	     ((eq (caar x) 'mexpt)
	      (newvarmexpt x (caddr x) t))
	     ((eq (caar x) 'mquotient)
	      (ratquotient (savefactors (prep1 (cadr x)))
			   (savefactors (prep1 (caddr x)))))
	     ((eq (caar x) 'mminus)
	      (ratminus (prep1 (cadr x))))
	     ((eq (caar x) 'rat)
	      (cond (modulus (cons (cquotient (cmod (cadr x)) (cmod (caddr x))) 1))
		    (t (cons (cadr x) (caddr x)))))
	     ((eq (caar x) 'bigfloat)(bigfloat2rat x))
	     ((eq (caar x) 'mrat)
	      (cond ((and *withinratf* (memq 'trunc (car x)))
		     (throw 'ratf nil))
		    ((catch 'compatvl
		       (progn (setq temp (compatvarl (caddar x)
						     varlist
						     (cadddr (car x))
						     genvar))
			      t))
		     (cond ((memq 'trunc (car x))
			    (cdr ($taytorat x)))
			   ((and (not $keepfloat)
				 (or (pfloatp (cadr x)) (pfloatp (cddr x))))
			    (cdr (ratrep* ($ratdisrep x))))
			   ((sublis temp (cdr x)))))
		    (t (cdr (ratrep* ($ratdisrep x))))))
	     ((assolike x genpairs))
	     (t (setq x (littlefr1 x))
		(cond ((assolike x genpairs))
		      (t (newsym x))))))

;;;;;;;;;;;;;;;;;;;;
;;(in-package "maxima")
;;;...... tree stuff from nt.lisp
;; maintain in synchrony a tree tr = *great-tree*
;;  and a  hash table ht = *great-hash*.
;; we have a comparison function that we use for insertion into the
;; tree, but to find out the ordering of two entries a b, we use.

;; (> (gethash a ht) (gethash b ht))

;; if we don't know if the items are in the hashtable, we do
;; (> (etree-si a #'comp) (etree-si b #'comp))
;;
;; note: this uses just binary search tree.

;; search: o(1) if already there.
;; insert: usual o(log n) time, worst-case o(n) time

(defvar *great-tree* nil)
(defvar *great-hash* nil)

;;**********************************************************************;;

(defstruct (etree (:print-function etree-print))
  (left  nil) (right nil)  number data)

;; print a etree in infix order, indenting sub-etrees

(defun etree-print (etree str ind)
  (unless (null etree)
	  (etree-print (etree-left etree) str (1+ ind))
	   ;; tab to ind column, then print this datum
	  (format str "~vt~s~%" ind (etree-data etree))
	  (etree-print (etree-right etree) str (1+ ind))))


;;**********************************************************************;;

(defun etree-si (data comp)
  ;search or insert
  (or (gethash data *great-hash*)(etree-i data comp)))


;; etree-i  insert
;;inserts item in etree and  returns the order number 
;; assigned to  that data in our ordering.  updates hash table.

(defun etree-i (data comp
		&aux found lo hi
		     (lowest most-negative-fixnum)
		     (highest most-positive-fixnum)
		     (root *great-etree*)
		     (ht *great-hash*)
		     )
  (declare (special *great-hash* *great-etree*))
  
;; the value returned from etree-i is an integer n
;;  such that n=(gethash data ht).

;;  the 2 local programs use the globals
;;  data  (stuff we are inserting)
;;  comp  (a function to be used for comparison on data. returns -1, 0, 1)
;;  lo, hi, span the lowest and highest data nodes nearby.
;;  found is the value to be returned on exit.

  (labels
      ((etree-i-help (etree)
	 ;;(format t "~%etree-i-help etree=~s"etree)
	 (if (null etree)(make-node)
	   (case (funcall comp data (etree-data etree))
		 (0
		  ;; we have been asked to search in the etree for
		  ;; an expression that is apparently equal to
		  ;; but not eq to the node in the etree. this
		  ;; is an error, so we warn
		  (format t 
			  "warning ~s was given to the ordering function 
and was not unique." data)
		  ;; and then we install a pointer from the new data
		  ;; to the old comparison order. (this may eventually
		  ;; break if we re-space the hash table)
		  (setf (gethash data ht)
			(setq found(gethash(etree-data etree) ht)))
		  (throw 'etree-i found))
		 (1
		   (setf lo (etree-data etree))
		   (setf (etree-right etree)
		     (etree-i-help (etree-right etree)))
		   (throw 'etree-i found))
		 (-1
		  (setf hi (etree-data etree))
		  (setf (etree-left etree)
		    (etree-i-help (etree-left etree)))
		  (throw 'etree-i found)))))
       ;;end of etree-i-help

       ;; make-node always returns a new node, but sets found to the 
       ;; order value
       (make-node nil
		   (let (new loval hival newval)
		     (setq new (make-etree :data data))
		     (setf loval (gethash lo ht lowest))
		     (setf hival (gethash hi ht highest))
		     (setf newval (truncate (+ loval hival) 2))
		     (cond((= newval loval)
			   ;(format t "warning: ran out of unique numbers in ordering ~s~%" data)
			   ;; reset global etree and hash table
			 (respace)
			 (setf found
			   (etree-i data comp))
			 (throw 'etree-i found)))
		     (setf found newval)
		     (setf (gethash data ht) newval)
		     new ;; return 1 level so that we can put node in etree.
		     )))
    ;; body of etree-i
  ;  (trace make-node etree-i-help)
  (catch 'etree-i (etree-i-help root))
  found))

;; this changes two global variables *great-etree* and *great-hash*

(defun respace ( &aux (etree *great-etree*) (table *great-hash*) )
 (let*
    ((size (hash-table-count table))
     (arr (etree2arr etree size))		;temp array
     (ht  (make-hash-table :size (* 2 size) :test #'eq)) ;new table
     (increment (truncate  (- most-positive-fixnum most-negative-fixnum) 
			   (+ size 1))))

   ;; set up the new values in the hash table
   (do ((j 0 (1+ j))
	(newval (+ increment most-negative-fixnum)
		(+ newval increment)))
       ((= size j) nil)
     (declare (fixnum j increment newval))
     ;(format t "hash[~s] set to ~s~%" j newval)
     (setf (gethash (svref arr j) ht) newval))
   ;; set up the new etree
   (setf *great-etree* (arr2etree arr))
   ;; point to the new hash-table
   (setf *great-hash* ht)
  'alldone
   ))

;;**********************************************************************;;
;;test stuff

(defun comp (a b)
  (cond ((= a b )0)((> a b) 1)(t -1)))

(defun gcomp(a b)
  (cond ((eq a b) 0)((great a b) 1)(t -1)))

;;initialize hash and etree

(progn nil (setq *great-etree* (make-etree :data 0))
       (setq *great-hash* (make-hash-table :test #'eq))
       (setf (gethash 0 *great-hash*)0))

(defun l (a) ;return the order number of a
  (etree-si a #'comp))

(defun ci() ;check it out
  (format t "~% hash=~%")
  (listhash *great-hash*)
)

;; look at hash table

(defun listhash(h)(maphash #'(lambda(key val)(format t "~s -> ~s~%"
						    key val)) h))



;;arr2etree convert a (presumably sorted) array a into a etree.

(defun arr2etree(a)
  (labels ((a2t (lo hi)
	      (cond ((> lo hi) nil)
		    ((= lo hi) (make-etree :data (svref a lo)))
		    (t (let((m (truncate (+ lo hi)2)))
			 (make-etree :data (svref a m) 
				    :left (a2t lo (1- m))
				    :right (a2t (1+ m) hi)))))))
	(a2t 0 (1- (length a)))))

;;etree2arr converts a (presumably sorted) etree of h entries
;; to an array a[0..(h-1)].

(defun etree2arr (tr size) ;; size = number of entries in the etree
  (let ((res (make-array size)) (count 0))
    (labels ((t2a(node)
		 (unless (null node)
			 (t2a (etree-left node))
			 (setf (svref res count) (etree-data node))
			 (incf count)
		 ;; (if (= count size) "etree was bigger than u thought")
			 (t2a (etree-right node)))))
	    (t2a tr)
	    res)))

;;(in-package "maxima")
;; revised for hash stuff
;; this breaks  (a=b)*c,  matrix*scalar, and rational function stuff

(defun simptimes (x w z)  ; w must be 1
  (prog (res)
	;check eqnflag matrixflag sumflag not used here
	(if (endp (cdr x)) (return 1)) ;case of ((mtimes) ) = 1
	(if (endp (cddr x))(return(cadr x)))
	
	;; for extra comments, compare to simplus
;	(setq check x)
	(setq res (make-hash-table :size 8 :test #'eq))
   start(setq x (cdr x))
   (cond
    ;; stuff about multiplying a matrix by zero will have to be
    ;; re-examined sometime.
        
    ((null x) (go end)))
   (setq w (if z (car x) (simplifya (car x) nil)))
   ;; unused label st1  
   (cond
    ((atom w) nil)
    ;; dealing with mrat also will have to be re-examined
    ;; dealing with mequal, mlist, $matrix, %sum  ditto
    )
   (if (eql w 0) (return 0))
   (setq res (tmshash w res))
   (go start)
end (return (hash2list* res))
	
;; there are 3 more conds that have to do
;; with expansion, and treatment of sums, equations, etc. omitted here.
	
))


(defun tmshash(x table) ; x is to be inserted into the times-hashtable table

;; there are these 5 cases for x, which is already simplified.

;; x= atomic name e.g.. alpha

;; x= ((mtimes)   <stuff>)  e.g. ((mtimes) 3 ((mexpt) y 3)) or
;;                               ((mtimes) u v)

;; x= ((mplus) ....)

;; x= ((mexpt) ...)  e.g. ((mexpt) y 3)  or ((mexpt) y ((mplus) a b))
  
;; x= (<otherop> {args})  e.g. (($foo) bar)


;; the first and 5th are treated the same.

  ;; the rule <anything>^0 =1   is applied ruthlessly and without
  ;; discretion, so infinity^0 or 0^0, if it gets here, comes out 1.
  

  (cond((atom x)
	(typecase x
		  (number
		   ;; mult all numbers into hashtable value for key=1
		   ;; numbers include rational, complex, float.
		   ;; that is, value^1 is part of the product
		   (let((v (* x (gethash 1 table 1))))
		     (if (eql v 1) (remhash 1 table) ;dont save 1^1 value
		       (setf (gethash 1 table) v))))
		  (hash-table
		   ;; this would be a plus hash-table !
		   ;; we don't expand, in this version of tmshash
		   (tmshash (hash2list+ x) table))
		  ;; mult other atoms into their own buckets by
		  ;; incrementing the value.  that is, x^2 is
		  ;; key=x, value =2,  multiply by x to get
		  ;; key=x, value =3.
		  ;; x^a is key=x, value=a, multiply by x to get
		  ;; key=x, value=a+1
		  (t  
		   (let((v (mfast+ 1 (gethash x table 0))))
		     (if (eql v 0) (remhash x table) ; x^0 is 1, remove it
		       (setf (gethash x table) v)))
		   ))
	     )
       ;;at this point x is not atomic, so it has a caar.
       ;; (caar x) = mtimes may have been taken care of by simptimes?
       ;; check anyway
       ((eq (caar x) 'mtimes) (mapc #'(lambda(z)(tmshash z table))(cdr x)))
       ((not(eq (caar x) 'mexpt)) ;; treat same as atom
	;; 
	(let((v (mfast+ 1 (gethash x table 0))))
	  (if (eql v 0) (remhash x table) ;<stuff>^0 is 1, remove it
	    (setf (gethash x table) v)))
	)			

       ;; the only other case of interest is mexpt
       (t				; (caar x) = 'mexpt
	(let ((b (cadr x))
	      (e (caddr x))             ;b is the base, e the exponent
					;list ((mexpt )b e)
		 v)
	  (cond 
		((and (consp e)
		      (eq (caar e) 'mtimes)
		      (numberp (cadr e)))
		 ;; this is the case of, e.g.  z^(3*stuff)
		 ;; change to (z^stuff)^3.
		 (setq b (m^ b (if (cdddr e)
				   ;;((mtimes simp) 3 a b) -> 
				   ;; ((mtimes simp) a b)
				   (ucons *mts*(cddr e))
				 ;; ((mtimes simp) 3 z) -> z
				 (caddr e))))
		 (setq v (mfast+ (cadr e) (gethash b table 0)))
		 (if (eql v 0) (remhash b table) ;b^0 is 1, remove it
		   (setf (gethash b table) v)))
		;; remaining cases  x^(a+b) or 2^x or y^z
		(t
		 ;; this is the case of, e.g. z^(a+b)
		 (setq v (mfast+ e (gethash b table 0)))
		 (if (eql v 0) (remhash b table) ; remove b^0
		   (setf (gethash b table) v))))
		 )))
  table)

(defun hash2list* (h)
  (let (arr (i 0) (size 0))
    (declare (fixnum i))
    (cond ((hash-table-p h)

	   ;; copy the hash-table into an array
	   (setq arr (make-array (setq size (hash-table-count h))))
	 ;;  (setq i 0)
	   (maphash #'(lambda(key val)
			;;each element in the array will be (key . val)
			(setf (svref arr i) (cons key val))
			(incf i)) h)

	   ;; now sort the elements of the array by keys (keys are distinct)

	   (sort arr #'(lambda(k1 k2)(great (car k1)(car k2))))

	   ;; now run through entries computing a simplified form.
	   ;; we are not unique-ifying the running sum or the products
	   ;; at the moment.  the keys will be unique, though.

;; we could do this in place in the array and then coerce to a list.
	   (setq h
		 (do ((i 0 (1+ i))
		(key)
		(val)
		(res)) ; three locals, initialized to nil
		     ;; the termination test
		     ;; if we want to save the hash version, we can do
		     ;; something like (cons (list 'mplus 'simp h) res)  below
	       ((= i size) (ucons *mts* res))  

	       (setf key (car (svref arr i)))
	       (setf val (cdr (svref arr i)))

	       (cond
		((eq val 1)		; key^1 -> key
		 (setq res (ucons key res)))
		((eq val 0))		; key^0 -> 1  should not be here
		((eq key 1)		;constant number ^1 -> number
		  (setq res (ucons val res)))
		;; all the other cases produce a ((mexpt simp) base expon)
		(t (setq res (ucons (ulist *mes*  key val) res)) 
		   ))))
	   ;; check for degenerate cases ((mtimes)) --> 1, ((mtimes) a)  --> a
	   (cond ((null (cdr h)) 1)
		 ((null (cddr h)) (cadr h))
		 (t h)))
	
	  (t h);; not a hashtable. just return it.
	  )))

(defun mfast+(x y)(if (and (numberp x)(numberp y)) (+ x y)(m+ x y)))



;; patch to hash.cl to provide hcons and hlist
;;
;;				-[mon jan 29 22:30:51 1990 by layer]-
;; 

;;; a pile of stuff for hcons specific to allegro CL 4.3 was in a file
;;; s:/newmaxima/hayden/proj/hcons.lisp.  Removed

(defun hcons (kcar kcdr *uniq-table*)
  (ucons kcar kcdr))

(defmacro uconsm(a b)`(hcons ,a ,b uniqtab))

;; (c) 1990, 1991, Richard J. Fateman
;; (c) 1994 Richard J. Fateman
;; re-edited  10/2003 RJF
;; non-standard hash table feature used below

#+allegro (defmacro eq-hash (object)
  "Gives us a hashing of an object such that (eq a b) implies
   (= (eq-hash a) (eq-hash b))"
  `(the fixnum (excl::pointer-to-positive-fixnum ,object)))

#+allegro (defun car-cdr-eq (key1 key2)
  ;; special hash-code function for the unique-cons table
  (declare (optimize (speed 3)(safety 0)(debug 0)))
    ;; this is the test to see if two conses have eq cars and eq cdrs
    (and (eq (eq-hash (car key1))(eq-hash(car key2)))
	 (eq (eq-hash (cdr key1))(eq-hash(cdr key2)))))

#+allegro (defvar *uniq-table* (make-hash-table :test 'car-cdr-eq))
#-allegro  ;often as good
(defvar *uniq-table* (make-hash-table :test 'eq))

(defvar *uniq-atom-table* (make-hash-table :test #'eql))
(defmfun $uniq(x)(uniq x))  
(defun uniq (x)
  "Return a canonical representation that is EQUAL to x,
  such that (equal x y) => (eq (uniq x) (uniq y))"
  (typecase x
    ((or fixnum symbol) x)
    (atom (or (gethash x *uniq-atom-table*)
              (setf (gethash x *uniq-atom-table*) x)))
    (cons (ucons (uniq (car x))   ; this could check in 
                                  ; *uniq-table* first...
                 (uniq (cdr x))))))
n
(defvar *fakecons* '(car . cdr))

(defun ucons (x y)
  "Unique cons: (eq (ucons x y) (ucons x y)) is always true."
  (declare (special *fakecons* *uniq-table*)
	   (optimize (speed 3)(safety 0)#+allegro (debug 0)))
  (let((temp *fakecons*)(tt *uniq-table*))
    (setf (car temp) x (cdr temp) y) ;don't allocate yet.
    (cond ((gethash temp  tt))
	  ;;If already there, great. 
	  (t (setf (gethash temp  tt) temp)
	     (setf *fakecons* (cons 'car 'cdr))
	     temp))))

(defun umapcar(f x)(cond((null x)nil)
			(t (ucons (funcall f (car x))(umapcar f (cdr x))))))


(defun uappend(r s)(cond ((null r)s)
			 (t (ucons (car r)(uappend (cdr r) s)))))



(setf *mts* (ulist 'mtimes 'simp))
(setf (get 'mtimes 'msimpind) *mts*)
(setf *mps* (ulist 'mplus 'simp))
(setf (get 'mplus 'msimpind) *mps*)
(setf *mes*  (ulist 'mexpt 'simp))
(setf (get 'mexpt 'msimpind)*mes*)
(setq bigfloatzero (uniq bigfloatzero))
(setq bigfloatone (uniq bigfloatone))

(defun maxima-constantp(x)(constantp x))  ;;; what should this really be??


(defun addk (xx yy)	;;xx and yy are assumed to be already reduced
  (setf xx(if (and (listp xx)(eq (caar xx) 'rat)) (/ (cadr xx)(caddr xx)) xx))
  (setf yy(if (and (listp yy)(eq (caar yy) 'rat)) (/ (cadr yy)(caddr yy)) yy))
  
 (cond ((equal xx 0) yy)
       ((equal yy 0) xx)
       ((and (numberp xx) (numberp yy)) (+ xx yy))
       ((or ($bfloatp xx) ($bfloatp yy)) ($bfloat (list '(mplus) xx yy)))
       (t (error "addk given non-constants to add"))))


#+ignore 
(defun addk (xx yy)	;;xx and yy are assumed to be already reduced
  (setf xx(if (and (listp xx)(eq (caar xx) 'rat)) (/ (cadr xx)(caddr xx)) xx))
  (setf yy(if (and (listp yy)(eq (caar yy) 'rat)) (/ (cadr yy)(caddr yy)) yy))
  
 (setf xx (cond ((equal xx 0) yy)
       ((equal yy 0) xx)
       ((and (numberp xx) (numberp yy)) (+ xx yy))
       ((or ($bfloatp xx) ($bfloatp yy)) ($bfloat (list '(mplus) xx yy)))
       (t (error "addk given non-constants to add"))))
 ;; old style
 (cond((rationalp xx)(if (= 1 (numerator xx)) xx (list '(rat)(numerator xx)(denominator xx))))
      (t xx)))


(defun redo()(load  (compile-file "c:/lisp/nnewsimp.cl")))