;; more hacking, 11/16/08 RJF
;; change short nodes to 2-element arrays
;; update elements by (setf (aref ..)  (insert-tree...))
;; instead of operating with destructive operations alone. Makes updating possible.

;; also, allow for different size nodes.

;;  see r4tree.lisp for comments   

;; try hacking more, declaring a few fixnums, perhaps overoptimistically
;; also, we have a program to pick out the smallest index, useful if
;; we want to program a DIVISION!!

(eval-when (:compile-toplevel :load-toplevel :execute)
  (declaim (optimize (speed 3)(safety 0)))
    
  (defparameter *mlogsizes* (list -8 -5 -3))
  ;; (defparameter *mlogsizes* (list -8 -8 -8)) ;alternative, shallower trees
  (defparameter *sizes* (mapcar #'(lambda(r)(expt 2 (- r))) *mlogsizes*))
  (defparameter *sizesm1* (mapcar #'(lambda(r)(1- (expt 2 (- r)))) *mlogsizes*))
  (nconc *sizes* (cddr *sizes*))
  (nconc *sizesm1* (cddr *sizesm1*))
  (nconc *mlogsizes* (cddr *mlogsizes*)) )

;; temporary, to make tracing possible

#|
(defparameter *mlogsizes* '(-8  -5 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3))
(defparameter *sizes*     '(256 32 8 8 8 8 8 8 8 8 8 8))
(defparameter *sizesm1*   '(255 31 7 7 7 7 7 7 7 7 7 7))
|#



(defmacro leafp(x) `(numberp ,x))
(defmacro nodep(x) `(arrayp ,x))
(defmacro makeshortnode(key val)`(vector ,key ,val))
(defmacro makefullnode (s)`(make-array ,s :initial-element 0))
(defmacro shortkey-key (h)`(aref ,h 0))
(defmacro shortkey-val (h)`(aref ,h 1))
(defmacro shortnodep(h)`(eq 2 (length ,h)))
(defmacro longnodep(h)`(< 2 (the fixnum (length ,h))))

(defparameter *debug* nil)

(defun db()(setf *debug* (null *debug*)))

(defun irt(key val tree)
 (insert-rtree key val tree *sizes* *sizesm1* *mlogsizes*))

#+ignore  ;; modify version to debug..
(defun insert-rtree(key val tree sizes sizesm1 mlogsizes &aux node)
  (declare(optimize(speed 3)(safety 0))
	  (type (simple-array t (*)) tree) (fixnum key))
  (cond 
   ((shortnodep tree)			; it is a single-entry node
    (cond ((eq  key (shortkey-key tree))
	   (setf (shortkey-val tree) val));; case where key matches exactly. replace value
	  (t
	   (let ((newtree (insert-rtree 
			   (shortkey-key tree) (shortkey-val tree) 
			   (makefullnode (car sizes)) 
			   (cdr sizes) (cdr sizesm1)(cdr mlogsizes))))
		  ;;; insert A NEW FULL-SIZE NODE HERE
	     (setf tree (insert-rtree 
			 key val newtree (cdr sizes) (cdr sizesm1)(cdr mlogsizes)))))))
   ;; the node (i.e. tree) is array >2 in length
   ;; next case is that the key is small enough to fit in this array
   ((< (the fixnum key)(the fixnum (length tree)))
    ;; we have found the level, or level-1 for the key.
    (cond ((nodep (setf node (aref tree key))); not a leaf
	   ;; if a node here, must descend one level further in tree; key now is 0
	   (setf (aref tree key) 
		 (insert-rtree 0 val node(cdr sizes)  (cdr sizesm1)(cdr mlogsizes)))); and update that node's location 0.
	  (t 
	   (setf (aref tree key) val)	)));; put it here.
   ;;The key has too many bits to insert at this level. 
   ;;Compute the subkey and separate the rest of the key by masking and shifting.
   (t (let ((ind (logand key (the fixnum(car sizesm1))));rightmost logsize bits
	    (newkey (ash key (the fixnum(car mlogsizes))))); reduced key
	(declare (fixnum ind newkey))
	(setf node (aref tree ind))
	(cond ((not (nodep node));; a single item, must be moved
	       (cond ((eq (the fixnum node) 0)
		      (setf (aref tree ind) (makeshortnode newkey  val)))
		     (t
		      (insert-rtree newkey val
				    (setf (aref tree ind) (makefullnode (car sizes)))
				  ;;  (setf (aref tree ind) (makefullnode (cadr sizes)))
				    (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
				    ;;sizes sizesm1 mlogsizes	
				    )
		      (insert-rtree 0;;ind 
				    node (aref tree ind)
				    (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
				   ;; sizes sizesm1 mlogsizes				  
				    ))))
	      (t
	       (setf (aref tree ind)
		     (insert-rtree 
		      newkey val node 
		      (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
		    ;;  	    sizes sizesm1 mlogsizes				  
		      )))))))
  tree)


;; debug this

(defun insert-rtree(key val tree sizes sizesm1 mlogsizes &aux node)
  (declare(optimize(speed 3)(safety 0))
	  (type (simple-array t (*)) tree) (fixnum key))
  
  (cond 
   ((shortnodep tree)			; it is a single-entry node
    (cond ((eq  key (shortkey-key tree))
	   (setf (shortkey-val tree) val));; case where key matches exactly. replace value
	  (t
	   (let ((newtree (insert-rtree 
			   (shortkey-key tree) (shortkey-val tree) 
			   (makefullnode (car sizes))  ;;***
			   sizes sizesm1 mlogsizes
			 ;;  (cdr sizes) (cdr sizesm1)(cdr mlogsizes)
			   )))
		  ;;; insert A NEW FULL-SIZE NODE HERE
	     (setf tree (insert-rtree 
			 key val newtree 
			 ;;(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
			 sizes sizesm1 mlogsizes
			 ))))))
   ;; the node (i.e. tree) is array >2 in length
   ;; next case is that the key is small enough to fit in this array
   ((< (the fixnum key)(the fixnum (length tree)))
;;      (assert (= (length tree) (car sizes)))

    ;; we have found the level, or level-1 for the key.
    (cond ((nodep (setf node (aref tree key))); not a leaf
	   ;; if a node here, must descend one level further in tree; key now is 0
	   (setf (aref tree key) 
		 (insert-rtree 0 val node(cdr sizes)  (cdr sizesm1)(cdr mlogsizes)))); and update that node's location 0.
	  (t 
	   (setf (aref tree key) val)	)));; put it here.
   ;;The key has too many bits to insert at this level. 
   ;;Compute the subkey and separate the rest of the key by masking and shifting.
   (t (let ((ind (logand key (the fixnum(car sizesm1))));rightmost logsize bits
	    (newkey (ash key (the fixnum(car mlogsizes))))); reduced key
	(declare (fixnum ind newkey))
	(setf node (aref tree ind))
	(cond ((not (nodep node));; a single item, must be moved
	       (cond ((eq (the fixnum node) 0)
		      (setf (aref tree ind) (makeshortnode newkey  val)))
		     (t
		      (insert-rtree newkey val
				    (setf (aref tree ind) (makefullnode (car sizes)))
				  ;;  (setf (aref tree ind) (makefullnode (cadr sizes)))
				    (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
				    ;;sizes sizesm1 mlogsizes	
				    )
		      (insert-rtree 0;;ind 
				    node (aref tree ind)
				    (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
				   ;; sizes sizesm1 mlogsizes				  
				    ))))
	      (t
	       (setf (aref tree ind)
		     (insert-rtree 
		      newkey val node 
		      (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
		    ;;  	    sizes sizesm1 mlogsizes				  
		      )))))))
  tree)

(defun urt(key val tree)
  (update-rtree key val tree *sizes* *sizesm1* *mlogsizes*))

(defun update-rtree(key val tree sizes sizesm1 mlogsizes &aux node)
  (declare(type (simple-array t (*)) tree) (fixnum key))
  (cond 
   ((shortnodep tree)			; it is a single-entry node
    (cond ((eq  key (shortkey-key tree))
	   (incf (shortkey-val tree) val));; case where key matches exactly. replace value
	  (t
	   (let ((newtree (insert-rtree 
			   (shortkey-key tree) (shortkey-val tree) 
			   (makefullnode (car sizes)) 
			   sizes sizesm1 mlogsizes
			   ;;(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
			   )))
		  ;;; insert A NEW FULL-SIZE NODE HERE
	     (setf tree (update-rtree 
			 key val newtree
			 ;;(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
			  sizes sizesm1 mlogsizes
			 ))))))
   ;; the node (i.e. tree) is array >2 in length
   ;; next case is that the key is small enough to fit in this array
   ((< (the fixnum key)(the fixnum (length tree)))
    ;; we have found the level, or level-1 for the key.
    (cond ((nodep (setf node (aref tree key))); not a leaf
	   ;; if a node here, must descend one level further in tree; key now is 0
	   (setf (aref tree key) 
	     ;; (insert-rtree 0 val node(cdr sizes) (cdr sizesm1)(cdr mlogsizes))
	      (update-rtree 0 val node(cdr sizes) (cdr sizesm1)(cdr mlogsizes))
	     ))
					; and update that node's location 0.
	  (t 
	   (incf (aref tree key) val)	)));; put it here.
   ;;The key has too many bits to insert at this level. 
   ;;Compute the subkey and separate the rest of the key by masking and shifting.
   (t (let ((ind (logand key (the fixnum (car sizesm1))));rightmost logsize bits
	    (newkey (ash key (the fixnum (car mlogsizes))))); reduced key
	(declare (fixnum ind newkey))

	(setf node (aref tree ind))
	(cond ((not (nodep node));; a single item, must be moved
	       (cond ((eq node 0)
		      (setf (aref tree ind) (makeshortnode newkey  val)))
		     (t
		      (update-rtree newkey val
				  ;;  (setf (aref tree ind) (makefullnode (car sizes)))
				   (setf (aref tree ind) (makefullnode (cadr sizes)))
				    
				    (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
				    ;;sizes sizesm1 mlogsizes
				    )
		      (update-rtree 0;;ind 
				    node (aref tree ind)
				    (cdr sizes)(cdr sizesm1)(cdr mlogsizes)
				    ;;sizes sizesm1 mlogsizes
				    ))))
	      (t
	       (setf (aref tree ind)
		     (update-rtree 
		      newkey val node 
		     (cdr sizes)  (cdr sizesm1)(cdr mlogsizes)
		      ;;sizes sizesm1 mlogsizes
		      )))))))
  tree)


(defun qrt(key  tree)  ;; get value for key in rtree tree. 0 means "no entry"
    (query-rtree key  tree *sizes* *sizesm1* *mlogsizes*))
  
(defun query-rtree(key tree sizes sizesm1 mlogsizes) ;; returns 0 if no entry in the tree for that key
  (declare (fixnum key))
  (cond 
   ((leafp tree) 0)
   ((shortnodep tree)
    (if (eq (shortkey-key tree) (ash key (car mlogsizes)))(shortkey-val tree) 0))
   ;; case 1: the key should be in this node.
   ((< key (the fixnum (car sizes)))
    (let ((z (aref tree key)))
    (if (nodep z);; if the node is an array
	(query-rtree 0 z (cdr sizes)(cdr sizesm1)(cdr mlogsizes));; look in the zero location
     z)));; the node is NOT an array: this is the value.
   (t (let* ((ind(logand key (car sizesm1)));;otherwise compute the offset
	     (h (aref tree ind)))	; and the subtree to use for further search
	(if h (query-rtree (ash key (car mlogsizes)) h (cdr sizes)(cdr sizesm1)(cdr mlogsizes)) 0))))) 



(defun map-rtree (fn mt);; order appears jumbled, actually it is in a reversed-radix order
  ;; apply fn, a function of 2 arguments, to key and val, for each entry in tree.
  ;; order is "radix reversed"
   (labels 
       ((tt1 (path displace mt mlogsizes);; path = path to the key's location
	  (declare (fixnum path displace))
	  (cond ((leafp mt) 
		 (unless (= 0 mt)(funcall fn path mt)));; function applied to key  and val.
	       ((shortnodep mt)
		(funcall fn (+ path (ash (the fixnum (shortkey-key mt))displace )) 
			 (shortkey-val mt)))
	       (t;; must be an array
		(let ((m (length mt)))
		  (declare (fixnum m))
		(do ((i 0 (1+ i)))
		    ((= i m ))
		  (declare (fixnum i))
		  (tt1 (+ (the fixnum (ash i displace)) path) 
		       (- displace (car mlogsizes)) ;um, fix this?
			  (aref mt i) (cdr mlogsizes))))))))
     (tt1 0 0 mt *mlogsizes* )))

(defun ordertree-nz(mt)
  (declare (type (simple-array t (*)) tree))
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (declare  (type (array t (*)) ans))
    (map-rtree #'(lambda(key val)
		   (vector-push-extend (cons key val) ans))
	       mt)
    (sort ans #'(lambda(a b)(declare (fixnum a b)) (< a b)) ;faster because of declarations?
	  ;;#'< 
	  :key #'car)))

(defun mul-rtree4 (r s)			;multiply two polynomials, in arrays
  (declare (optimize (speed 3)(safety 0)))
  (let* ((er (car r)) (es (car s))
	 (cr (cdr r)) (cs (cdr s))
	 (cri 0) (eri 0)
	 (ler (length er))
	 (les (length es))
	 (ans (makefullnode (car *sizes*))))
    (declare (type (simple-array integer (*)) cr cs er es)
	     (fixnum ler les))		;maybe some others are fixnums too..
    ;; do the NXM multiplies into the answer tree
    (dotimes (i ler (A2PA (ordertree-nz ans)))
      (declare (fixnum i))
      (setf cri (aref cr i) eri(aref er i))
      (dotimes (j les)
	(declare (fixnum j))
	(urt (+ eri(aref es j)) (* cri(aref cs j)) ans
		      )))))

;; we could save space/time by doing a:=a+b*c  in one step if we are doing GMP arith.


(defun mul-rtree4-nosort(r s)			;multiply two polynomials, in arrays
  (declare (optimize (speed 3)(safety 0)))
  (let* ((er (car r)) (es (car s))
	 (cr (cdr r)) (cs (cdr s))
	 (cri 0) (eri 0)
	 (ler (length er))
	 (les (length es))
	 (ans (makefullnode (car *sizes*))))
    (declare (type (simple-array integer (*)) cr cs er es)
	     (fixnum ler les))		;maybe some others are fixnums too..
    ;; do the NXM multiplies into the answer tree
    (dotimes (i ler ans ) ;;; (A2PA (ordertree-nz ans))) ;;dont sort
      (declare (fixnum i))
      (setf cri (aref cr i) eri(aref er i))
      (dotimes (j les)
	(declare (fixnum j))
	(urt (+ eri(aref es j)) (* cri(aref cs j)) ans)))))



(defun A2PA(A) ;; array of (exp . coef) to 2 arrays.
  (let* ((LA (length A))
	 (V nil)
	 (anse (make-array LA))
	 (ansc (make-array LA)))
    (dotimes (i LA (cons anse ansc))
      (setf V (aref A i))
      (setf (aref anse i) (car V))
      (setf (aref ansc i) (cdr V)))))


#|

;; example loop iteration..
;; (loop for i across z when (> i 0) do (print i))
;; (loop for i across z maximizing i)
	  

;; this picks the item with the smallest index and returns a pair:
;; (index . val)
;; useful for division!


(defun map-rtreexx (fn mt);; order appears jumbled, actually reversed-radix order
  ;; apply fn, a function of 2 arguments, to key and val, for each entry in tree.
  ;; order is "radix reversed"
   (labels 
       ((tt1 (path displace mt);; path = path to the key's location
	  (declare (fixnum path displace) (special *sm*))
	  
	  (cond ((leafp mt) (unless (= 0 mt)(funcall fn path mt)
				    ));; function applied to key  and val.
	       ((shortnodep(aref mt 0))
		(funcall fn (+ (ash (the fixnum (aref mt 2)) displace) path) (aref mt 1)))
	       (t;; must be an array
		(do* ((i 0 (1+ i))
		     (path2 path (+ (ash i displace) path)))
		    ((= i #.(expt 2 logsize)))
		  (declare (fixnum i))
		  ;; we cut this off if
		  ;; (aref mt i) is a leaf and index is too big
		  (if (and (leafp (aref mt i))
			       (> path2 *sm*)) (return) ;from do*
		  (tt1  path2 (+ displace logsize)(aref mt i))))))))
    (tt1 0 0 mt)))


'(defun smallest(mt)
  (let ((*sm* most-positive-fixnum)
	(val 0))
    (declare (special *sm*))
    (map-rtreexx #'(lambda(k v)
		    ;		   (format t "~%testing k=~s"k)
		     (if (< k *sm*)
			 (setf *sm* k val v)))
		 mt)
   (cons *sm* val)))

|#

;;; testing, compare simply insertion of random numbers into a structure

#|
(setf mt (makefullnode 256))
(irt 3 3 mt)
(irt 259 259 mt)
(irt 16777216 16777216 mt)
|#					;testing, comparing to hash table

(defun t1(N r)
  (let ((k 0)
	(mt (makefullnode 256)))
    (declare (fixnum k))
    (dotimes (i N mt)
    (setf k (random r))
      (irt k k mt))))

(defun t2(N r);; use hash table
  (let ((k 0)
	(ht (make-hash-table)))
    (declare (fixnum k))
    (dotimes (i N ht)
    (setf k (random r))
      (setf (gethash k ht) k))))

(defun t3(N r);; just do random
  (let ((k 0)
	(ht (make-hash-table)))
    (declare (fixnum k))
    (dotimes (i N ht)
      (setf k (random r))
      (fff k) ;; calling a dummy function so compiler does not remove previous line.
   )))
(defun fff(k) k)

(defun t4(N r);; use array.  sometimes faster than t3!
  (let ((k 0)
	(ar (make-array r :initial-element 0)))
    (declare (fixnum k) (type (simple-array t (*)) ar))
    (dotimes (i N ar)
    (setf k (random r))
      (setf (aref ar k) k))))

;;  (t1 1000000 2000000)is 1031ms
;;  (t2 1000000 2000000)is  781ms; more memory
;;  (t3 1000000 2000000)is  140ms  [just calling (random), then dummy fn.
;;  (t4 1000000 2000000)is  203ms


;;  (t1 1000000 20000)is 656ms, redo 344
;;  (t2 1000000 20000)is 516ms; more memory  , redo 281
;;  (t3 1000000 20000)is 140ms , redo 188 
;;  (t4 1000000 20000)is 156ms , redo 125 


;;  (t1 1000000 2000)is 250ms
;;  (t2 1000000 2000)is 344ms; more memory , redo 250
;;  (t3 1000000 2000)is 140ms, redo 125
;;  t4 is 156 msg ; 

;;  (t1 1000000 255)is 188ms, redo 203
;;  (t2 1000000 255)is 288ms; much more memory, redo 266
;;  (t3 1000000 255)is 140ms, redo 187
;; t4 is 172; redo 109 (!)  [cost of array update is less than dummy function call]


;; subtract off the time for the loop and the generation of random numbers..

;;  conclusion: hashing is competitive for insertion, maybe 35-40% faster often
;;  highlight: rtree is faster by 3X if all elements are in the first level node

;;   for items of size up to 20,000 compare 516 vs 376, favoring hashing
;;   for items of size up to  2,000 compare 110 vs 104, favoring rtree
;;   for items of size up to    255 compare  48 vs 148, favoring rtree
;;   for items of size up to 
;;                        2,000,000 compare 891 vs 641, favoring hashing

;; BUT

;; even if we have to allocate an array of size 2,000,000, t4 using arrays
;; is fastest. with times of 43ms, 16ms, 16ms, 32ms.  [we redid these timings
;; with a longer (X10) loop. They are in the right ballpark.]

;; putting this into the mix, using an array is between 14X faster and 3X faster
;; for insertion.



;; even if we have to allocate an array of size 20,000,000, t4 using arrays
;; is fastest. with times of 43ms, 16ms, 16ms, 32ms.  [we redid these timings
;; with a longer (X10) loop. They are in the right ballpark.]

;; putting this into the mix, using an array is between 14X faster and 3X faster
;; for insertion.

;; what about an array of size 200,000,000 ?  My lisp system refused to allocate it.
;; Error: An allocation request for 800000016 bytes caused a need for [too much memory from OS]

;; however, the other guys still worked.
;;  (time (prog nil  (t1 1000000 200000000) nil))
;; took 1703 ms. (radix)
;; t2 1015ms (hash) 
;; t3 172 ms (empty loop)
;; t4 can't do it.


;; t1 used ;  484,367 cons cells, 56,991,112 other bytes, 0 static bytes
;; t2 used ; ;     22 cons cells, 63,469,160 other bytes, 0 static bytes

;; conclude: simple array is a winner unless it can't fit. otherwise hashing is good.



#| more timing.
 cl-user(32):  (setf q (make-racg 10000 50 5)r (make-racg 10000 50 5))
(#(45 58 82 130 152 202 212 241 285 305 ...)
 . #(5 2 1 4 5 1 5 1 3 4 ...))
 cl-user(33): (time (mul-rtree4 q r))
; cpu time (non-gc) 42,718 msec user, 0 msec system
; cpu time (gc)     954 msec user, 15 msec system
; cpu time (total)  43,672 msec user, 15 msec system
; real time  44,344 msec
; space allocation:
;  584,645 cons cells, 16,295,368 other bytes, 0 static bytes
(#(94 107 128 131 141 165 178 179 191 201 ...)
 . #(5 2 20 1 8 4 20 4 8 5 ...))

 cl-user(38): (time (setf ans3 (mul-heap q r)))
; cpu time (non-gc) 64,171 msec (00:01:04.171) user, 0 msec system
; cpu time (gc)     47 msec user, 0 msec system
; cpu time (total)  64,218 msec (00:01:04.218) user, 0 msec system
; real time  65,938 msec (00:01:05.938)
; space allocation:
;  280 cons cells, 10,928,736 other bytes, 0 static bytes
(#(94 107 128 131 141 165 178 179 191 201 ...)
 . #(5 2 20 1 8 4 20 4 8 5 ...))
[1] cl-user(39): (length (car ans3))
510519
[1] cl-user(40): (aref (car ans3) (1- *))
511991
[1] cl-user(41): (time (mul-naivef q r))
; cpu time (non-gc) 49,844 msec user, 0 msec system
; cpu time (gc)     422 msec user, 16 msec system
; cpu time (total)  50,266 msec user, 16 msec system
; real time  51,266 msec
; space allocation:
;  1,022,065 cons cells, 4,090,064 other bytes, 8,192 static bytes
(#(94 107 128 131 141 165 178 179 191 201 ...)
 . #(5 2 20 1 8 4 20 4 8 5 ...))


;; try this  (setf q (make-racg 5000 50 5)r (make-racg 1000 50 5))
(time (mul-dense-fix q r)) is 0.078 sec
(time (mul-dense q r))     is 0.219 sec
(time (mul-fft q r))       is 1.41 sec
(time (mul-rtree4 q r))    is 1.625 sec
(time (mul-hash1 q r))     is 1.687 sec
(time (mul-heap5 q r))     is 2.53 sec  ;; in newmulr2s.fasl
(time (mul-naive q r))     is 4.03 sec

(setf q (make-racg 5000 500 5) r (make-racg 1000 500 5))
(time (mul-dense-fix q r))    0.203 sec
(time (mul-dense q r))        0.375 sec
(time (mul-heap5 q r))        3.37 sec
(time (mul-hash1 q r))        8.1 sec
(time (mul-rtree4 q r))       8.4 sec
(time (mul-naive q r))      129.19 sec
(time (mul-fft q r))			;stopped. out of storage.

(setf q (make-racg 5000 1 5)r (make-racg 1000 1 5))
(time (mul-fft q r))          0.016 sec  
(time (mul-dense-fix q r))    0.031 sec
(time (mul-dense q r))        0.188 sec
(time (mul-toom q r))         0.218 sec
(time (mul-naive q r))        0.30  sec
(time (mul-karat q r))        0.500 sec
(time (mul-rtree4 q r))       0.67  sec
(time (mul-hash1 q r))        0.81  sec
(time (mul-heap5 q r))        2.17  sec	; from newmulr2s
(time (mul-heap5 q r))        2.76  sec  ; from newmulr2

(setf q (make-racg 5000 10000 5)r (make-racg 1000 10000 5)) ;; really sparse
(time (mul-heap5 q r))       4.703 sec
(time (mul-hash1 q r))      22.875 sec   329 megabytes
(time (mul-rtree4 q r))     22.875 sec 
(time (mul-naive q r))     320.812 sec  much less space...   74 megabytes

|#


;;; try to put the heap stuff from newmulr2s in this file too


;;; -*- Lisp -*-
;;;;; excerpt from newmul with Roman's version of chaining, sort of.
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))
;;; this tries to do f:=f+a*b  instead of c:=a*b; f:=f+c. could be faster with GMP.

;;; attempt to shorten because the "insert at top" heuristic doesn't really work often enough.
  
(defmacro heap-empty-p ()
    "Returns non-NIL if & only if the heap contains no items."
    `(<= fill-pointer 1))

(defmacro lesser-child (a parent fp)
    "Return the index of the lesser child.  If there's one child,
 return its index.  If there are no children, return 
 (FILL-POINTER (HEAP-A HEAP))."
    `(let* ((left;;(+ (the fixnum ,parent)(the fixnum  ,parent))
	     (ash ,parent 1))
	    (right (1+ (the fixnum left))))
       (declare (fixnum left right ,fp ,parent )
		(optimize(speed 3)(safety 0))
		(type (simple-array t (*)) ,a))
       (cond ((>= left ,fp) ,fp)
	     ((= right ,fp) left)
	     ;;; use the line below if exponents are fixnums!!
	     ((< (the fixnum (st-exp (aref ,a left)))(the fixnum (st-exp (aref ,a right)))) left)
	     ;; ((< (st-exp (aref ,a left))(st-exp (aref ,a right))) left) ;more general
	     (t right))))

;; put together in a lexical context to make access faster to global vars
(let ((*YYEXPS* nil)
      (*YYCOEFS* nil)
      (*YYLEN* 0)
      (fill-pointer 0))
  (declare (fixnum *YYLEN* fill-pointer)
	   (type (simple-array t (*)) *YYEXPS* *YYCOEFS*))

  (defstruct (st) 
    ;; A kind of stream.  each structure contains 6 items, not 5 as in previous version
    ;; which by reference to arrays of XX and YY , can burp out
    ;; successive exponents and coefficients of the product.
  
    (exp 0 :type integer)
    (coef1 0 :type integer)  ;; could be doubles, or gmp numbers or ..
    (coef2 1 :type integer)  ;; the second item,  to be multiplied by coef1 and added..
    (i 0 :type fixnum)  
    (v 0 :type integer) (e 0 :type integer)); v, e are read-only 

  (defun make-simple-array (A) (if (typep A 'simple-array) A 
				 (make-array (length A)
					     :element-type (array-element-type A)
					     :initial-contents A)))
  
  (defun mularZ4(x y)  ;; set up initial heap.
    (declare (optimize (safety 3)(speed 0)))
    (let* ((xexps (car x))
	   (xcoefs (cdr x))
	   (yexps (make-simple-array (car y)))
	   (ycoefs(make-simple-array (cdr y)))
	   (yexp0 (aref yexps 0))
	   (ycoef0 (aref ycoefs 0))
	   (ans (make-array (1+(length xexps))))
	   (ylen (length (car y)))
	   (xe 0)
	   (v 0))			;
      (declare (type (simple-array t (*))
		     ans xexps xcoefs yexps ycoefs))
      (setf *YYEXPS* yexps *YYCOEFS* ycoefs *YYLEN* ylen)
      (setf (aref ans 0) (make-st :exp 0 :coef1 0 :coef2 1 :i ylen :v 0 :e 0))
      (dotimes (j (length xexps) ans)
	(setf v (aref xcoefs j))
	(setf xe (aref xexps j))
	(setf (aref ans (1+ j)) 
	  (make-st :exp (+ xe yexp0) :coef1 (* v ycoef0) :i 1 :v v :e xe) ))))

  (defun percolate-down4 (aa hole xx fp)
    "Private. Move the HOLE down until it's in a location suitable for X.
Return the new index of the hole."

    (let ((xxval (st-exp xx)))
      (declare (fixnum hole fp)(integer xxval)
	       (type (simple-array t (*)) aa)(optimize (speed 3)(safety 0)))
      (do* ((child (lesser-child aa hole fp) (lesser-child aa hole fp))
	    (aaval (aref aa child)(aref aa child)))
	  ((or (>= (the fixnum child) fp) 
	       (< (the fixnum xxval) (the fixnum(st-exp aaval))) ;faster, less general
	       
	   ;;    (< xxval (the fixnum(st-exp aaval))) ;faster, less general
	       
	    ;;   (< xxval (st-exp aaval))
	       )
	   hole)
	(declare (fixnum child) (type st aaval))
	(setf (aref aa hole)  aaval  hole child))))
  
  ;; important: initial-contents of heap must have an element 0 that
  ;; is not really used. for convenience we could have it equal to the first element.
  ;; e.g.  (3 5 8)  should be (3 3 5 8)
    
  (defun heap-insert4 (a xx);; xx is new element
    "Insert a new element into the heap a. Heap is changed";; rjf
    (declare (type (simple-array t (*)) a)(fixnum fill-pointer)
	     (integer newe)
	     (optimize (speed 3)(safety 0)))
    
    (let ((c nil)			; potential chain point
	  (fp (incf fill-pointer)) )
      (declare (type (simple-array t (*)) a)(fixnum hole fp fill-pointer)
	       (optimize (speed 3)(safety 0)))
      (decf fp)
      (setf (aref a  fp) xx) 
      (do ((i fp parent)
	   (parent (ash fp -1);; initialize to parent of fp
		   (ash parent -1);; and its parent, next time
		   ))
	  ((>= (st-exp xx) (st-exp(aref a parent)));; is xx>= parent? if so, exit
	   a);; inserted the item at location i.
	(declare (fixnum i parent))
	(setf c (aref a parent))
	(cond
	 ((= (st-exp xx) (st-exp c));; aha, we can insert this item by chaining here.
	  (incf (st-coef c) (* (st-coef1 xx) (st-coef2 xx)))
	  ;; Line above is one op-code in GMP, using no persistent storage outside c...
					;(format t "~%              to ~s " c)
					;	  (format t "~% heap is ~s" a)
	  (setf (aref a fill-pointer) nil)
	  (decf fill-pointer);; we don't need the extra heap spot for this item 
	  (if (setf xx  (next-term xx))	;try inserting the successor to this item then.
	      (heap-insert4  a  xx))	
	  (return a)) ;exit from do loop
	 (t  (setf (aref a i) c 
		   (aref a parent) xx) a)))))
	
  (defun heap-remove4(a) 
    ;; returns nil if the heap is empty, otherwise the removed item is returned
    (declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
    (if (heap-empty-p) nil
        (let*  ((xx (aref a 1))
		(fp (decf fill-pointer)) ;faster, fp is local now
		(last-object (aref a fp)))
	  (declare (fixnum fp))
	  (setf (aref a fp) nil)  ;; optional. makes tracing neater
	  (setf (aref a (percolate-down4 a 1 last-object fp)) last-object)
	  xx)))
  
   ;; this is the multiplication program, eliminated: peeking at top
 (defun mul-heap5(x y)
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ4 x y)(mularZ4 y x)))
	   (alen  10)
	   ;; initial heap
	   (anse (make-array alen :adjustable t :element-type 'integer :fill-pointer 0))
	   (ansc (make-array alen :adjustable t :element-type 'integer :fill-pointer 0))
	   (preve 0) (prevc 0)
	   (top nil)
	   (newe 0)(newc1 0) (newc2 1) )
      (declare (integer preve prevc newe newc1 newc2) 
	       (type (array t(*)) anse ansc hh))
      
      (setf fill-pointer  (length hh))
      (loop while  (not (heap-empty-p))
	    do
	    (setf top (heap-remove4 hh))
	    (setf newe (st-exp top) newc1 (st-coef1 top) newc2 (st-coef2 top))
	    (cond   ((= preve newe)		; is this same exponent as previous?
		     (incf prevc (* newc1 newc2))) ; if so, add coefficient, and loop
		    (t (unless (= prevc 0) ; normalize: remove zero coeffs
			 (vector-push-extend prevc ansc)
			 (vector-push-extend preve anse))
		       (setf prevc (* newc1 newc2)) ; initialize new coefficient
		       (setf preve newe)))
	    
	    (if (setf top (next-term top));; is there a successor?
		(heap-insert4 hh top)));; if so, insert it in heap
      ;; end of loop, push out the final monomial before exit from routine.
      (vector-push-extend prevc ansc)
      (vector-push-extend preve anse)
      (cons anse ansc)))
 
 (defun next-term(rec)   	; record is a stream 5-tuple  return succesor or nil
    (declare(type (simple-array t(*)) *YYCOEFS* *YYEXPS*)
	    (optimize (speed 3)(safety 0)) 
	    (type st rec))
    (let ((i (st-i rec)))
      (cond ((>= i *YYLEN*) nil);; no more items on stream of products
	    (t (incf (the integer (st-i rec)))
	      ;;; (setf (st-coef rec)(* (st-v rec)(aref *YYCOEFS* i)))
	       
	       (setf (st-coef1 rec)(st-v rec))
	       (setf (st-coef2 rec)(aref *YYCOEFS* i)) ;; don't multiply yet, save for later
	       ;; later we can do  c:=c+a*b in one fell swoop.
	       	       
	       (setf (st-exp rec) 
		     (+ (st-e rec)(aref *YYEXPS* i)))
	       rec)))) 
 )




#|

 (defun h(z) (setf u (make-racg z 10000 5)v (make-racg z 10000 5)) 
			 (setf ansr (mul-rtree4 u v) ansn (mul-naive u v))
			 (pol= ansr ansn))
			 ;; test for same answer..
			 
|#