;;; -*- Lisp -*-
;;; msortx.lisp multiplication of polyns based on copying coefs into arrays.
;;; this uses double-float coefficients

(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; sort X+Y various ways.

(defun make-rac(n d) 
  ;;make random list of length n with spacing k, with 0<k<=d between elements
  ;; with coefficients that are random integers between 1 and d
  (let ((exps (make-array n :element-type 'fixnum))
	(coefs(make-array n :element-type 'double-float))
	;(coefs(make-array n :element-type 'fixnum))
	
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	;(+ 1 (random d))
	(+ 1.0d0 (random d))
	)
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))


(defun test(n m)
  (setf x (make-rac n 5))
  (setf y (make-rac m 5))
;  (format t "~% mulhashb time-ar")
;    (time (mulhashb x y))
  (format t "~% mulhashe time-ar")
    (time (mulhashe x y)))

;; multiply two polynomials via a hash table. For comparison
(defparameter *zcoefs-lim* 50)
(defparameter *zcoefs* (make-array *zcoefs-lim* :element-type 'double-float))

#+ignore ;; bug because it relies on *zcoefs* being large enough. see mulhashe
(defun mulhashc (x y);;  this hash table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table))
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs *zcoefs*)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
	     (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
	     ;;(type (simple-array fixnum(*)) xcoefs ycoefs)
	     (fixnum xei xy zcounter)
	     (double-float xco)
	     (optimize (speed 3)(safety 0))
;;	     (:explain :types)
	     )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
	(cond((null  loc) 
	      ;; if loc is nil, allocate a space in the answer array
	      (setf (gethash xy ans) zcounter)
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)) ; should check for overflow or something, here.
	     (t ;; there is an entry in hash table for this exponent
	      (incf (aref zcoefs loc) ;; no number consing here either??
		    (*  xco (aref ycoefs j)))))))))

(defun sorthash(h) ;; utility for mulhasha
  (let ((ans nil))
    (maphash #'(lambda (i v)(push (cons i v) ans)) h)
    (sort ans #'< :key #'car)))

(defun sorthashb(h) ;; utility for mulhashb
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)(setf (cdr v)(car v))
		       (setf (car v) i)
		       (unless (= (cdr v) 0.0d0) ;; normalize! remove 0 coefs
			 (vector-push-extend v ans)))
	     h)
    (sort ans #'< :key #'car)))

(defun sorthashc(h ar) ;; utility for mulhashc
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)
		 (setf v (aref ar v)); get the double-float, by indirection
		       (unless (= v  0.0d0) ;; normalize! remove 0 coefs
			 (vector-push-extend (cons i v) ans)))
	     h)
    (sort ans #'< :key #'car)))

#|(test 1000 500)

 hash time-ar
; cpu time (non-gc) 109 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  125 msec user, 0 msec system
; real time  141 msec
; space allocation:
;  4,519 cons cells, 12,496,136 other bytes, 0 static bytes
  ...)
  using fixnums for coefficients, the time is somewhat lower, but
  dramatically less space is used.

 cl-user(163):  (time (mulhashb x y))
; cpu time (non-gc) 78 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,569 cons cells, 282,152 other bytes, 0 static bytes

					
The real problem with coefficients here is that lisp doesn't have 
a convenient unit operation  a:=a+b*c which destroys a in place.

HERE
we can make that happen by doing  (incf (aref a i) (* (aref b j)(aref c k)))
for suitable a,b,c declared arrays.

Using double floats this way,

 (time (mulhashe x y)) ;; for double-floats
; cpu time (non-gc) 62 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
					;  4,513 cons cells, 388,952 other bytes, 0 static bytes

					

;; still about 50% of the time is taken in hash-code stuff, giving
;; us an opportunity to do something else.


BUT USING POTENTIAL BIGNUMS...
(setf xx (make-racf 1000 5)) ;; this is important. (make-rac ..) produces simple-arrays of doubles for coefs
(setf yy (make-racf 500 5))

[2] cl-user(391): (time (mulhashf xx yy))
; cpu time (non-gc) 78 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,513 cons cells, 536,208 other bytes, 0 static bytes
#((3 . 20) (6 . 20) (8 . 10) (9 . 12) (11 . 10) (13 . 25) (14 . 10)



|#




;;; cleaned up version that doesn't run out of array space..

;; because mulhashd cannot declare zcoefs to be simple-array, this
;; is much slower than hashc.  See hashd that reallocates "by hand"
#+ignore
(defun mulhashd (x y);;  this hash table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table))
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs  (make-array 5 :element-type 'double-float :adjustable t :fill-pointer 0))
	(loc nil)
	(xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
	     (type (simple-array double-float(*)) xcoefs ycoefs) ; not zcoefs
	     (fixnum xei xy zcounter)
	     (double-float xco)
	     (optimize (speed 3)(safety 0))
	     )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))

      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
	(cond((null  loc) 
	      ;; if loc is nil, allocate a space in the answer array
	      ;; first, point to that spot from the hash table.
	      (setf (gethash xy ans) (fill-pointer zcoefs))
	      ;; then store the result in zcoefs array.
	      (vector-push-extend (* xco (aref ycoefs j)) zcoefs) ; no number consing??
	      )				; prevents overflow by extending zcoefs.
	     (t ;; there is an entry in hash table for this exponent
	      (incf (aref zcoefs loc) ;; no number consing here either??
		    (*  xco (aref ycoefs j)))))))))



;; try to do this right...
(defun mulhashe (x y);;  this hash table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table :test 'eq)) ;; eq is good enough for fixnums for ACL
	;;(ans (make-hash-table :test 'eq)) ;; more general, if exponents are arbitrary bignums
	
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs *zcoefs*)
	(zlim *zcoefs-lim*)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
	     (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
	     (fixnum xei xy zcounter zlim)
	     (double-float xco)
	     (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
	(cond (loc ;; there is an entry in hash table for this exponent. update it
	       ;; this happens about n*m times, where n=size(x), m=size(y).
	       (incf (aref zcoefs loc) ;; no number consing here either??
		     (*  xco (aref ycoefs j))))
	      (t
	       ;; if loc is nil, allocate a space in the answer array.
	       ;; This happens relatively infrequently, about k times, where k=size(x*y).
	      (setf (gethash xy ans) zcounter)
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)
	      (cond ((= zcounter zlim) ;; if we reach this, we need more working space.
		     
		     #+ignore  ;; if we do this, the array is not a simple-array and we lose efficiency
		     (let*((newsize (ash zlim 1))
			   (newarray (adjust-array zcoefs newsize :element-type 'double-float)))
		        (setf zcoefs newarray 
			     *zcoefs* newarray 
			     *zcoefs-lim* newsize 
			     zlim newsize))
			
		       
			   
		     (let* ((newsize (ash zlim 1)) ; double the size
			    (newarray (make-array newsize :element-type 'double-float)))
		       (declare (fixnum newsize)
				(type (simple-array double-float(*)) newarray))
		    ;;   (format t "~%zcoefs now size ~s" newsize)
		       (dotimes (m newsize)
			 (declare (fixnum m))
			 (setf (aref newarray m)(aref zcoefs m)))
		       
		       (setf zcoefs newarray 
			     *zcoefs* newarray 
			     *zcoefs-lim* newsize 
			     zlim newsize))
		     		     
		     )))
	      )))))



;; attempt to do same as mulhashe but with avl tree.

#+ignore
(defun multree (x y)
  ;;  this version uses tree not hashtable. it does not store double-floats, but indexes into an array of df.
  (let ((ans (tree-init 0 0)) ;whatever kind of tree.
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs *zcoefs*)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
	     (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
	     (fixnum xei xy zcounter)
	     (double-float xco)
	     (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf loc (treesearch ans (setf xy (+ xei(aref yexps j))))) ;look for entry at index xy
	(cond((null loc) 
	      ;; if loc is nil, allocate a space in the answer array
	      (treeinsert ans xy zcounter) ; insert into ans, at index xy, the record zcounter
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)
	      (cond ((= zcounter *zcoefs-lim*) ;; we need more working space.
		     (let* ((newsize (ash *zcoefs-lim* 1)) ; double the size
			    (newarray (make-array newsize :element-type 'double-float)))
		       (declare (fixnum newsize))
		    ;;   (format t "~%zcoefs now size ~s" newsize)
		       (dotimes (m newsize)(setf (aref newarray m)(aref zcoefs m)))
		       (setf zcoefs newarray *zcoefs* newarray *zcoefs-lim* newsize)))))
	     
	     (t ;; there is an entry in tree for this exponent
	      (incf (aref zcoefs loc) ;; no number consing here either??
		    (*  xco (aref ycoefs j)))))))))


;; hey, this is really pretty good!
;; return to coefs  being any number type
(defun mulhashf (x y);;  this hash table does not store double-floats, but indexes into an array of df.
  (let ((ans (make-hash-table :test 'eql)) ;; eq is good enough for fixnums for ACL
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs (make-array 10))
	(zlim 10)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
	     (type (simple-array t (*)) xcoefs ycoefs zcoefs)
	     (fixnum xei xy zcounter zlim)
	    ; (double-float xco)
	     (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
	(cond (loc ;; there is an entry in hash table for this exponent. update it
	       ;; this happens about n*m times, where n=size(x), m=size(y).
	       (incf (aref zcoefs loc) ;; no number consing here either??
		     (*  xco (aref ycoefs j))))
	      (t
	       ;; if loc is nil, allocate a space in the answer array.
	       ;; This happens relatively infrequently, about k times, where k=size(x*y).
	      (setf (gethash xy ans) zcounter)
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)
	      (cond ((= zcounter zlim) ;; if we reach this, we need more working space.
			   
		     (let* ((newsize (ash zlim 1)) ; double the size
			    (newarray (make-array newsize)))
		       (declare (fixnum newsize)
				(type (simple-array t(*)) newarray))
		    ;;   (format t "~%zcoefs now size ~s" newsize)
		       (dotimes (m newsize)
			 (declare (fixnum m))
			 (setf (aref newarray m)(aref zcoefs m)))
		       (setf zcoefs newarray    zlim newsize))
		     		     
		     )))
	      )))))

(defun make-racf(n d) 
  ;;make random list of length n with spacing k, with 0<k<=d between elements
  ;; with coefficients that are random integers between 1 and d
  (let ((exps (make-array n :element-type 'fixnum))
	(coefs(make-array n))
	;(coefs(make-array n :element-type 'fixnum))
	
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	;(+ 1 (random d))
	(+ 1 (random d))
	)
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))



(defun make-racg(n d) 
  ;;make random list of length n with spacing k, with 0<k<=d between elements
  ;; with coefficients that are random integers between 1 and d
  (let ((exps (make-array n )) ;; not guaranteed to be fixnums any more
	(coefs(make-array n))
	;(coefs(make-array n :element-type 'fixnum))
	
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	;(+ 1 (random d))
	(+ 1 (random d))
	)
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))


(defparameter *zcoefs-limg* 50)
(defparameter *zcoefsg* (make-array *zcoefs-limg*))

(defun mulhashg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs *zcoefsg*)
	(zlim *zcoefs-limg*)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
	     (type (simple-array t (*)) xcoefs ycoefs zcoefs) ; any coefs
	     (fixnum zcounter zlim)
	     (integer xei  xy)
	    ; (double-float xco)
	     (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
	(declare (fixnum j))
	;; look in the hash table at location for this exponent
	(setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
	(cond (loc 
	       ;; There is an entry in hash table for this
	       ;; exponent. Update corresponding array location. Note:
	       ;; hashtable is NOT updated.  this happens about n*m
	       ;; times, where n=size(x), m=size(y).
	       (incf (aref zcoefs loc) ;; no number consing here either??
		     (*  xco (aref ycoefs j))))
	      (t
	       ;; if loc is nil, allocate a space in the answer array.
	       ;; This happens relatively infrequently, about k times, where k=size(x*y).
	       ;; Hashtable is updated.
	      (setf (gethash xy ans) zcounter)
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)
	      (cond ((= zcounter zlim) 
		     ;; if we reach this, we need more working space.
		     ;; This is rarely used code. If global array is big enough, it is never used.
		     ;; If global array is too small, it is made bigger.
		     (let* ((newsize (ash zlim 1)) ; double the size
			    (newarray (make-array newsize)))
		       (declare (fixnum newsize)
				(type (simple-array t(*)) newarray))
		    ;;   (format t "~%zcoefs now size ~s" newsize)
		       (dotimes (m newsize)
			 (declare (fixnum m))
			 (setf (aref newarray m)(aref zcoefs m)))
		       (setf zcoefs newarray  zlim newsize
			     *zcoefs-limg* zlim 
			     *zcoefsg* zcoefs))))))))))

;; will this be slower??? Yes.
;; put numbers IN the hashtable. Problem is the frequent re-hashing at line ***

#+ignore
(defun mulhashh (x y);;   exponents may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eq is good enough for bignums too
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array t  (*)) xexps yexps) ; not fixnums guaranteed.
	     (type (simple-array t (*)) xcoefs ycoefs)
	     (fixnum zcounter zlim)
	    ; (double-float xco)
	     (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps) (sorthash ans))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
	(cond (loc ;; there is an entry in hash table for this exponent. update it
	       ;; this happens about n*m times, where n=size(x), m=size(y).
	       (incf (gethash xy ans) ;;***
		     (*  xco (aref ycoefs j))))
	      (t	      (setf (gethash xy ans)  (*  xco (aref ycoefs j)))))))))


;;; try skiplists
;;;no, this is really slow
(defun mulhashi (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (sl-init)) ;skip list
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(finger nil))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
	     (type (simple-array t (*)) xcoefs ycoefs) ; any coefs
	     (integer xei )
	     (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)ans ) ;the skiplist
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i))	;coefficient corresponding to that exponent
      (setf finger ans)
      (dotimes (j (length yexps))
	(declare (fixnum j))
	(setf finger(SLupdate-finger  ans finger (+ xei (aref yexps j))  (* xco (aref ycoefs j))))

	))))

(defun addar(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (cond 
   ((= (length(car x)) 0) y)
   ((= (length(car y)) 0) x)
   (t
    (let* 
	((xexps (car x))
	 (xi 0)
	 (xlim  (1- (length xexps)))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (yi 0)
	 (ylim  (1- (length yexps)))
	 (ycoefs (cdr y))
	 (zexps (make-array (+ xlim ylim 1) :adjustable t :fill-pointer 0))
	 (zcoefs (make-array (+ xlim ylim 1) :adjustable t :fill-pointer 0)))
      (declare (fixnum xi yi ylim xlim)(optimize (speed 3)(safety 0)))
      (loop 
       (cond 
	;; check first if we have run out of x or y,
	((> xi xlim)
	 (if (> yi ylim) (return (cons zexps zcoefs))
	   (let()
	     (loop for i from yi to ylim do 
		   (vector-push (aref yexps i) zexps)
		   (vector-push (aref ycoefs i) zcoefs))
	     (setf yi (1+ ylim)))))
	((> yi ylim)
	 (let()
	   (loop for i from xi to xlim do 
		 (vector-push (aref xexps i) zexps)
		 (vector-push (aref xcoefs i) zcoefs))
	   (setf xi (1+ xlim))))
		
	((< (aref xexps xi)(aref yexps yi))
	 (vector-push (aref xexps xi)  zexps)
	 (vector-push (aref xcoefs xi) zcoefs)
	 (incf xi))
	((= (aref xexps xi)(aref yexps yi))
	 (vector-push (aref xexps xi) zexps)
	 (vector-push (+ (aref ycoefs yi) (aref xcoefs xi)) zcoefs)
	 (incf xi) (incf yi))
	(t;;(> (aref xexps xi)(aref yexps yi))
	 (vector-push (aref yexps yi)  zexps)
	 (vector-push (aref ycoefs yi) zcoefs)
	 (incf yi))))))))

(defun splitpol(z) ;; take a polynomial and return its two halves
  (let* ((h (length (car z)))
	 (half (ash h -1))
	 (rest (- h half)))
    (declare (fixnum h half rest)(optimize (speed 3)(safety 0)))
    (values
	 (cons (make-array half :displaced-to (car z))
		     (make-array half :displaced-to (cdr z)))
	 (cons (make-array rest :displaced-to (car z) :displaced-index-offset half)
		     (make-array rest :displaced-to (cdr z) :displaced-index-offset half)))))

(defun mular(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (let* ((xexps (car x))
	 (xe 0)
	 (xc 0)
	 (xlen  (length xexps)))
    (declare (fixnum  xlim)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
	   '(#() . #())) ; empty polynomial
	   ((= 1 xlen)			; a monomial times a polynomial!
	   (setf xe (aref xexps 0))
	   (setf xc (aref (cdr x) 0))
	   (cons (map 'vector #'(lambda (r)(+ r xe)) (car y))
		 (map 'vector #'(lambda (r)(* r xc)) (cdr y))))
	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		(addar (mular low y)(mular hi y)))))))
;; BROKEN
#+ignore

(defun addar(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;;; version with no fill pointers..
  (cond 
   ((= (length(car x)) 0) y)
   ((= (length(car y)) 0) x)
   (t
    (let* 
	((xexps (car x))
	 (xi 0)
	 (xlim  (1- (length xexps)))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (yi 0)
	 (ylim  (1- (length yexps)))
	 (ycoefs (cdr y))
	 (m (+ xlim ylim)) ; maximum possible size.
	 (zexps (make-array m))
	 (zi 0)
	 (zcoefs (make-array m)))
	 
      (declare (fixnum xi yi ylim xlim zi)(optimize (speed 3)(safety 0)))
      (loop 
       (cond 
	;; check first if we have run out of x or y,
	((null (aref xexps xi)) ;run out of x
	 (if (null (aref yexps yi)) (return (cons zexps zcoefs))
	   (let()
	     (loop for i from yi to ylim do 
		   (if (null (aref yexps i))(return-from addar (cons zexps zcoefs)))
		   (setf (aref zexps zi) (aref yexps i))
		   (setf (aref zcoefs zi) (aref ycoefs i))
		   (incf zi))
	     (setf yi (1+ ylim)))))
	((> yi ylim)
	 (let()
	   (loop for i from xi to xlim do 
		 (if (null (aref xexps i))(return-from addar (cons zexps zcoefs)))
		    (setf (aref zexps zi) (aref xexps i))
		    (setf (aref zcoefs zi) (aref xcoefs i))
		    (incf zi))
	   (setf xi (1+ xlim))))
		
	((< (aref xexps xi)(aref yexps yi))
	 (setf (aref zexps zi)(aref xexps xi))
	 (setf (aref zcoefs zi)(aref xcoefs xi))
	 (incf zi)
	 (incf xi))
	((= (aref xexps xi)(aref yexps yi))
	 (setf (aref zexps zi) (aref xexps xi))
	 (setf (aref zcoefs zi) (+ (aref ycoefs yi) (aref xcoefs xi)))
	 (incf xi) (incf yi) (incf zi))
	(t;;(> (aref xexps xi)(aref yexps yi))
	 (setf (aref zexps zi)(aref yexps xi))
	 (setf (aref zcoefs zi)(aref ycoefs xi))
	 (incf zi)
	 (incf yi))))))))


;; this works. not particularly compact. or fast. mulhashg is much faster.
(defun addar(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (cond 
   ((= (length(car x)) 0) y)
   ((= (length(car y)) 0) x)
   (t
    (let* 
	((xexps (car x))
	 (xi 0)
	 (xlim  (1- (length xexps)))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (yi 0)
	 (ylim  (1- (length yexps)))
	 (ycoefs (cdr y))
	 (m (+ xlim ylim)) ; maximum possible size.
	 (zexps (make-array m :fill-pointer 0))
	 (zcoefs (make-array m :fill-pointer 0)))
      (declare (fixnum xi yi ylim xlim)(optimize (speed 3)(safety 0)))
      (loop 
       (cond 
	;; check first if we have run out of x or y,
	((> xi xlim)
	 (if (> yi ylim) (return (cons zexps zcoefs))
	   (let()
	     (loop for i from yi to ylim do 
		   (vector-push (aref yexps i) zexps)
		   (vector-push (aref ycoefs i) zcoefs))
	     (setf yi (1+ ylim)))))
	((> yi ylim)
	 (let()
	   (loop for i from xi to xlim do 
		 (vector-push (aref xexps i) zexps)
		 (vector-push (aref xcoefs i) zcoefs))
	   (setf xi (1+ xlim))))
		
	((< (aref xexps xi)(aref yexps yi))
	 (vector-push (aref xexps xi)  zexps)
	 (vector-push (aref xcoefs xi) zcoefs)
	 (incf xi))
	((= (aref xexps xi)(aref yexps yi))
	 (vector-push (aref xexps xi) zexps)
	 (vector-push (+ (aref ycoefs yi) (aref xcoefs xi)) zcoefs)
	 (incf xi) (incf yi))
	(t;;(> (aref xexps xi)(aref yexps yi))
	 (vector-push (aref yexps yi)  zexps)
	 (vector-push (aref ycoefs yi) zcoefs)
	 (incf yi))))))))


(defun mulhashg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs *zcoefsg*)
	(zlim *zcoefs-limg*)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
	     (type (simple-array t (*)) xcoefs ycoefs zcoefs) ; any coefs
	     (fixnum zcounter zlim)
	     (integer xei  xy)
	    ; (double-float xco)
	     (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashg ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
	(declare (fixnum j))
	;; look in the hash table at location for this exponent
	(setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
	(cond (loc 
	       ;; There is an entry in hash table for this
	       ;; exponent. Update corresponding array location. Note:
	       ;; hashtable is NOT updated.  this happens about n*m
	       ;; times, where n=size(x), m=size(y).
	       (incf (aref zcoefs loc) ;; no number consing here either??
		     (*  xco (aref ycoefs j))))
	      (t
	       ;; if loc is nil, allocate a space in the answer array.
	       ;; This happens relatively infrequently, about k times, where k=size(x*y).
	       ;; Hashtable is updated.
	      (setf (gethash xy ans) zcounter)
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)
	      (cond ((= zcounter zlim) 
		     ;; if we reach this, we need more working space.
		     ;; This is rarely used code. If global array is big enough, it is never used.
		     ;; If global array is too small, it is made bigger.
		     (let* ((newsize (ash zlim 1)) ; double the size
			    (newarray (make-array newsize)))
		       (declare (fixnum newsize)
				(type (simple-array t(*)) newarray))
		    ;;   (format t "~%zcoefs now size ~s" newsize)
		       (dotimes (m newsize)
			 (declare (fixnum m))
			 (setf (aref newarray m)(aref zcoefs m)))
		       (setf zcoefs newarray  zlim newsize
			     *zcoefs-limg* zlim 
			     *zcoefsg* zcoefs))))))))))
(defun sorthashg(h ar) ;; utility for mulhashgi
  (let ((ans (make-array (hash-table-count h) :fill-pointer 0)))
    (maphash #'(lambda (i v)
		 (setf v (aref ar v)); get the double-float, by indirection
		       (unless (= v  0.0d0) ;; normalize! remove 0 coefs
			 (vector-push (cons i v) ans)))
	     h)
    (breakup (sort ans #'< :key #'car))
    ))

(defun breakup (ar) 
  ;; ar is an array of (exp . coef).
  ;; separate them into a pair: one array of exponents and one of coefficients.
  (let* ((ll (length ar))
	(exps (make-array ll))
	 (coefs (make-array ll))
	 (tt nil))
    (dotimes (i ll (cons exps coefs))
      (setf tt (aref ar i))
      (setf (aref exps i) (car tt))
      (setf (aref coefs i)(cdr tt)))))


(defun PA2L (PA)  
  ;;  Pair of Arrays to List, length
  ;;  returns L a LIST of (exp . coef) of length n.
  ;; also returns n.

  (let* ((exps (car PA))
	 (coefs (cdr PA))
	 (n (length exps)))
    (declare (fixnum n))
    (do ((i (1- n) (1- i))
	 (ans nil (cons (cons (aref exps i)(aref coefs i)) ans)))
	((< i 0) (values ans n))
      (declare (fixnum i)))))



(defun L2PA (L n) 
       ;;List to Pair of Arrays
       ;; L is a LIST of (exp . coef) of length n.
       ;; create a pair: one array of exponents and one of coefficients.
  (let* ((exps (make-array n))
	 (coefs (make-array n))
	 (count 0))
    (declare (fixnum count))
    (dolist (i L (cons exps coefs))
      (setf (aref exps count) (car i))
      (setf (aref coefs count)(cdr i))
      (incf count))))
	      
(defun addarL(r s);; r and s are sparse polys as Lists of (exp . coef).
  ;; return a similar list of their sum. destructive.
  
   (declare (optimize (speed 1)(safety 3)(debug 3)))
  (let((ans nil) (pt nil))
    ;;set initial link
    (cond((null r)(return-from addarL s))
	 ((null s)(return-from addarL r))
	 ((< (caar r)(caar s))(setf ans r) (setf r (cdr r))) ;; which arg do we destroy?
	 (t (setf ans s)(setf s (cdr s))))
    (setf pt ans)
    
    (loop
      (cond((null r)(setf (cdr pt) s) (return ans))
	   ((null s)(setf (cdr pt) r) (return ans))
	   ((< (caar r)(caar s)) ;;key is car
	    (setf (cdr pt) r)
	    (setf r (cdr r)))
	   ((= (caar r)(caar s)) 
;	     (format t "~% r=~s s=~s pt=~s" r s pt)
	    
	    (setf (cdadr pt)(+ (cdar r)(cdar s)))
	    (setf (cdr pt) r)
	    (setf s (cdr s))
	    (setf r (cdr r)))
	   
	   (t (setf (cdr pt) s)
	      (setf s (cdr s))))
;;      (format t "~% r=~s s=~s pt=~s" r s pt)
      (setf pt (cdr pt)))))


(defun mularL(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  ;; this works but uses a lot of cons cells and is about 5X slower than mulhashg.
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (xe 0)
	 (xc 0)
	 (xlen  (length xexps))
	 (ans nil)
	 (i (1- (length yexps))))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
	   '(#() . #())) ; empty polynomial
	   ((= 1 xlen)			; a monomial times a polynomial!
	    (setf xe (aref xexps 0))
	    (setf xc (aref xcoefs 0))
	    ;; probably just as fast to do this next line rather than loop
	    ;;(map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) (car y)(cdr y))
	    (loop
	      (if (< i 0)(return t))
	      (push (cons (+ xe (aref yexps i)) (* xc (aref ycoefs i))) ans)
	      (decf i))
	    ans)

	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		(addarL (mularL low y)(mularL hi y)))))))


(defun mularL(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  ;; this works but uses a lot of cons cells and is about 5X slower than mulhashg.
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (xe 0)
	 (xc 0)
	 (xlen  (length xexps))
	 (ans nil)
	 (i (1- (length yexps))))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
	   '(#() . #())) ; empty polynomial
	   ((= 1 xlen)			; a monomial times a polynomial!
	    (setf xe (aref xexps 0))
	    (setf xc (aref xcoefs 0))
	    ;; probably just as fast to do this next line rather than loop
	    ;;(map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) (car y)(cdr y))
	    (loop
	      (if (< i 0)(return t))
	      (push (cons (+ xe (aref yexps i)) (* xc (aref ycoefs i))) ans)
	      (decf i))
	    ans)

	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		(addarL (mularL low y)(mularL hi y)))))))

#+ignore
;; subdivide BOTH polys!
(defun mularL(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  ;; this works but uses a lot of cons cells and is about 5X slower than mulhashg.
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (xe 0)
	 (xc 0)
	 (xlen  (length xexps))
	 (ans nil)
	 (i (1- (length yexps))))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
	   '(#() . #())) ; empty polynomial
	   ((= 1 xlen)			; a monomial times a polynomial!
	    (setf xe (aref xexps 0))
	    (setf xc (aref xcoefs 0))
	    ;; probably just as fast to do this next line rather than loop
	    ;;(map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) (car y)(cdr y))
	    (loop
	      (if (< i 0)(return t))
	      (push (cons (+ xe (aref yexps i)) (* xc (aref ycoefs i))) ans)
	      (decf i))
	    ans)

	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		;; this 2-way subdivision is a loser. answer is ok, but slow.
		(addarL (mularL y low )(mularL y hi)))))))




;; subdivide both polys but only to some limited depth.
;; 
(defvar *cachelim* 10  )
(defun dosmallmult(x y) (ptimes x y)) ;; make sure result is right form.
  
(defun mularLd(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (xe 0)
	 (xc 0)
	 (xlen  (length xexps))
	 (ans nil))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
	   nil) ; empty polynomial
	   ((= 1 xlen)			; a monomial times a polynomial!
	    (setf xe (aref xexps 0))
	    (setf xc (aref xcoefs 0))
	    (map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) yexps ycoefs))
	   ((< xlen *cachelim*)(if (< (length yexps) *cachelim*) 
				   ;; both x and y are small enough
				   (mular x y)
				 ;; x is small enough, y is not. split it.
				 (multiple-value-bind (low hi)
				     (splitpol y)
				   (addarL (mularL low x)(mularL  hi x)))))
	   ;; x is not small enough.
	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		(addarL (mularL low y )(mularL hi y)))))))

(defun testsplit(n z)
  (dotimes (i n)(splitpol z)))

(defun mularq(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; divide only down to cachelim
  
  (let* ((xexps (car x))
	 (yexps (car y))
	 (xe 0)
	 (xc 0)
	 (xlen  (length xexps)))
    (declare (fixnum  xlim)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
	   '(#() . #())) ; empty polynomial
	   ((= 1 xlen)			; a monomial times a polynomial!
	   (setf xe (aref xexps 0))
	   (setf xc (aref (cdr x) 0))
	   (cons (map 'vector #'(lambda (r)(+ r xe)) (car y))
		 (map 'vector #'(lambda (r)(* r xc)) (cdr y))))
	   ((< xlen *cachelim*)(if (< (length yexps) *cachelim*) 
				   ;; both x and y are small enough
				   (mulhashgg x y);; or something good for cache
				 ;; x is small enough, y is not. split it.

				      (multiple-value-bind (low hi)
					  (splitpol y)
					;;(format t "~%low=~s~%~%high=~s" low hi)
					(addar (mularq low x)(mularq hi x)))))
	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		(addar (mularq low y)(mularq hi y)))))))



;; same as mulhashg but simple-array --> array.  Slighly slower.

(defun mulhashgg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
	(xexps (car x))
	(xcoefs (cdr x))
	(yexps (car y))
	(ycoefs (cdr y))
	(xei 0)
	(xco 0)
	(zcoefs *zcoefsg*)
	(zlim *zcoefs-limg*)
	(zcounter 0)
	(loc nil)
	(xy 0))
    (declare (type (array t (*)) xexps yexps) ; any exponents
	     (type (array t (*)) xcoefs ycoefs zcoefs) ; any coefs
	     (fixnum zcounter zlim)
	     (integer xei  xy)
	    ; (double-float xco)
	     (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashg ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
	(declare (fixnum j))
	;; look in the hash table at location for this exponent
	(setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
	(cond (loc 
	       ;; There is an entry in hash table for this
	       ;; exponent. Update corresponding array location. Note:
	       ;; hashtable is NOT updated.  this happens about n*m
	       ;; times, where n=size(x), m=size(y).
	       (incf (aref zcoefs loc) ;; no number consing here either??
		     (*  xco (aref ycoefs j))))
	      (t
	       ;; if loc is nil, allocate a space in the answer array.
	       ;; This happens relatively infrequently, about k times, where k=size(x*y).
	       ;; Hashtable is updated.
	      (setf (gethash xy ans) zcounter)
	      ;; store the result in zcoefs array.
	      (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
	      (incf zcounter)
	      (cond ((= zcounter zlim) 
		     ;; if we reach this, we need more working space.
		     ;; This is rarely used code. If global array is big enough, it is never used.
		     ;; If global array is too small, it is made bigger.
		     (let* ((newsize (ash zlim 1)) ; double the size
			    (newarray (make-array newsize)))
		       (declare (fixnum newsize)
				(type (array t(*)) newarray))
		    ;;   (format t "~%zcoefs now size ~s" newsize)
		       (dotimes (m newsize)
			 (declare (fixnum m))
			 (setf (aref newarray m)(aref zcoefs m)))
		       (setf zcoefs newarray  zlim newsize
			     *zcoefs-limg* zlim 
			     *zcoefsg* zcoefs))))))))))



;;; make this prettier for paper

#+ignore
(defun mul-sac(x y)
  ;; mul-sparse-array-cache(x y)
  ;; x and y are sparse polys in arrays as produced by (make-racg n m)
  (let* ((xexps (car x))  ;; an array of exponents for x
	 (yexps (car y)))  ;; an array of exponents for y
    (cond  ((< (length xexps) *cachelim*)
	    (if (< (length yexps) *cachelim*) 
		;; both x and y are small enough
		(mul-sparse-inside-cache x y);; something good for cache
	      ;; x is small enough, y is not. split it.
	      (multiple-value-bind (low hi)
		  (splitpol y)
		(add-sparse-array (mul low x)(mul-sac hi x)))))
	   (t (multiple-value-bind (low hi)
		  (splitpol x)
		(add-sparse-array (mul-sac low y)(mul-sac hi y)))))))



;;; try to queue up subproblems so we can re-use scratch space..
(defvar *subprobs* nil)
(defun mularX(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (let* ((xexps (car x))
	 (xlen  (length xexps)))
    (declare (fixnum  xlim)(optimize (speed 3)(safety 0)))
    (cond ((= 1 xlen)			; a monomial times a polynomial!
	   (push (list (aref (car x) 0) (aref (cdr x) 0) y) *subprobs*))
	  (t (multiple-value-bind (low hi)
		 (splitpol x)
	       (mularX low y)(mularX hi y))))))

#+ignore
(defun cleanup (sp)
;; clean up queue of subproblems
  (let ((ans nil))			;an array??
    
    ;; pop a guy off *subprob*, execute the multiply,
    ;; push on a list to add them.  re-use arrays.
 (dolist (i *subprobs* ans)
	   (let ((xe (car i))
		 (xc (cadr i)))
	    ;; do this next part in fixed locations.
		 
	   (cons (map 'vector #'(lambda (r)(+ r xe)) (caaddr i))
		 (map 'vector #'(lambda (r)(* r xc)) (cdaddr i)))
	   ;;hmm, maybe should do this in pairs? or with heap?	   
	   (push (addar);; etc
		 )))))

   
#+ignore
(defun mulartop(x y)
  (let ((*subprobs* nil)
	(*ans*))
    (mularY x y)
    ;;hm.  one way to handle this is to form N streams
    ;; each stream has an exponent/coefficient pair on the top.
    ;; collect all these into a heap and add up coeffs at the top of each stream.
    ;; similar to what we tried with lists.
    (collect *subprobs*)
    ))


#+ignore
(defun mularY(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (let* ((xexps (car x))
	 (yexp0 (aref (car y) 0))
	 (ycoefs (cdr y))
	 (h (create-heap))
	 (anse (make-array 100 :adjustable t :fill-pointer 0))
	 (ansc (make-array 100 :adjustable t :fill-pointer 0))
	 r rexp ri xex xco)
	 
    (declare (optimize (speed 3)(safety 0)))
    (dotimes (i (1- (length xexps)))
      (insert-heap h  ;; key is CAR of list of 3 items
		   (list (+ yexp0(aref xexps i)) ;; lead x exponent + lead y exponent
			 i		; index counter into y array
			 (aref (car x) i) ; the x exponent
			 (aref (cdr x) i) ; the x coefficient
			 )))
    (vector-push 0 anse)  add 0*x^0 to the answer.
    (vector-push 0 ansc)
    (setf ansetop 0) 
    
    (loop (if (empty-heap h) (return (cons anse ansc)))
      (setf r (remove-heap-fast h))
      (setf rexp (car r))
      (setf ri   (cadr r))
      (setf xex  (caddr r))
      (setf xco  (cadddr r))
      (cond ((= rexp ansetop)		; same exponent. Add into ansc
	     (incf (ansc (1- (fill-pointer ansc)))
		   (* xco (aref ycoefs ri)))
	     (insert-heap h (list (+ xex (aref yexps ri))
				  (1+ ri)
				  xex
				  ??
				  )))))))


;;; hm. try creating a bunch of streams.



(defun mularZ(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (yexp0 (aref yexps 0))
	 (ycoef0 (aref ycoefs 0))
	 (ylim (length (car y)))
	 (ans nil))
	 
    (declare (optimize (speed 3)(safety 0)))
	 (dotimes (i (length xexps) ans)
	   (push
	    (list (+ yexp0 (aref xexps i)) ;; the lowest exponent
		   (* (aref xcoefs i)ycoef0) ; its coefficient
		  (let ((index -1)	; index into y
			(xexpi (aref xexps i))
			(xco (aref xcoefs i)))
			#'(lambda()
			    (incf index)
			    (if (>= index ylim) nil
			    (values (+ xexpi (aref yexps index)) ; the exponent
				    (* xco (aref ycoefs index))) ; the coef
			    ))))	    
	    ans)   )))

(defparameter *ht* (make-hash-table))
      
(defun update-st(w) ;; w looks like ( enumber  cnumber  function)
  (declare (special *ht*))
  (multiple-value-bind (e1 c1)
      (funcall (caddr w))
    ;; update the hash table, as a debugging tool..
    (cond (e1 (incf (gethash e1 *ht* 0) c1) ;; omitting this line makes it much faster
	      (setf (car w) e1 (cadr w) c1)
	      )
	  (t nil))))

(defun domuls (wlist)
  (dolist (i wlist *ht*)
    (loop while (update-st i))))

(defun mulstr(x y)  ;; oddly enough, this is a pretty good program!! 
  (let ((*ht* (make-hash-table)))
    (domuls (mularZ x y))
    (sorthashs *ht*)))


(defun sorthashs(h) ;; utility for mulstr
  (let ((ans (make-array (hash-table-count h) :fill-pointer 0)))
    (maphash #'(lambda (i v)
		 (vector-push (cons i v) ans))		 
	     h)
    (breakup (sort ans #'< :key #'car))    ))



(defun mularZ1(x y);; x and y are sparse polys in arrays as produced by (make-racg n m). return array
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (yexp0 (aref yexps 0))
	 (ycoef0 (aref ycoefs 0))
	 (ylim (length (car y)))
	 (ans (make-array (length xexps))))
	 
    (declare (optimize (speed 3)(safety 0))(type (simple-array t (*)) ans xexps xcoefs yexps ycoefs))
	 (dotimes (i (length xexps) ans)
	    (setf (aref ans i)
	      (list (+ yexp0 (aref xexps i)) ;; the lowest exponent
		   (* (aref xcoefs i)ycoef0) ; its coefficient
		  (let ((index -1)	; index into y
			(xexpi (aref xexps i))
			(xco (aref xcoefs i)))
			#'(lambda()
			    (incf index)
			    (if (>= index ylim) nil
			    (values (+ xexpi (aref yexps index)) ; the exponent
				    (* xco (aref ycoefs index))) ; the coef
			    ))))))))





;; -*- Lisp -*-
;;; see heap.lisp for docs.

(eval-when (compile load) (declaim (special fill-pointer)
				   (optimize (speed 3)(safety 0))))

;; heap is just an array.

;;(defstruct heap a  max-count)

(defmacro hlessfun (a b) `(< (car ,a) (car ,b)))

(defun percolate-down (aa hole x)
  "Private. Move the HOLE down until it's in a location suitable for X.
Return the new index of the hole."
    (declare (fixnum hole)(type (simple-array t (*)) aa)(optimize (speed 3)(safety 0)))
	
  (do ((child (lesser-child aa hole) (lesser-child aa hole)))
      ((or (>= child fill-pointer) (hlessfun x (aref aa child)))
       hole)
      (setf (aref aa hole) (aref aa child)
	    hole child)))

(defun percolate-up (a hole x)
  "Private.  Moves the HOLE until it's in a location suitable for holding
X.  Does not actually bind X to the HOLE.  Returns the new
index of the HOLE.  The hole itself percolates down; it's the X
that percolates up."
    (declare (type (simple-array t (*)) a)(fixnum hole)
	     (optimize (speed 3)(safety 0)))
    (setf (aref a 0) x)
    (do ((i hole parent)
	 (parent (ash hole -1);;(floor (/ hole 2))
		 (ash parent -1);;(floor (/ parent 2))
		 ))
	
	
	;; potential to speed up line below by declaration if a, x are fixnum,
	((not (hlessfun x (aref a parent))) i)
      (declare (fixnum i parent))

      (setf (aref a i) (aref a parent))))


  ;; important: initial-contents must have the first element duplicated.
  ;; e.g.  (3 5 8)  should be (3 3 5 8)

(defun heap-empty-p (heap)
  "Returns non-NIL if & only if the heap contains no items."
  (declare (ignore heap))
  (<= fill-pointer 1))

(defun heap-insert (a x)
  "Insert a new element into the heap.  Returns the heap." ;; rjf
  
    ;; Append a hole for the new element.
    ;; (vector-push-extend nil a) 
    ;; assume enough room...
    ;; (vector-push nil a)
   
   ;; (format t ".")
    (setf (aref a fill-pointer) nil)
    (incf fill-pointer)
    ;; Move the hole from the end towards the front of the
    ;; queue until it is in the right position for the new
    ;; element.
    (setf (aref a (percolate-up a (1- fill-pointer) x)) x))

(defun heap-remove-fast(a) 
  ;; assumes non-empty!! if empty next answer is bogus.
  ;; answer after that is an error (fill pointer can't pop)
 ;;   (format t "-")

  (let*  ((x (aref a 1))
	 (last-object			;(prog1 (aref a fill-pointer) (decf fill-pointer))
	  (progn (decf fill-pointer)(aref a fill-pointer) )))
    ;; could declare array a element-type if we knew it if was fixnum
    ;; here and elsewhere.
    (declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
    (setf (aref a (percolate-down a 1 last-object)) last-object)
	     x))

(defun heap-peek (heap)
  "Return the first element in the heap, but don't remove it.  It'll
 be Erroneous if the heap is empty.   (Should that be an error?)"
  (aref heap  1))

(defun lesser-child (a parent)
  "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))."
  (declare (optimize(speed 3)(safety 0))(type (simple-array t (*)) a) (fixnum parent))
  (let* ((left (ash parent 1)) ;;(* parent 2 )
         (right (1+ left))
         (fp fill-pointer))
    (declare (fixnum left fp right))
    (cond ((>= left fp) fp)
          ((= right fp) left)
          ((hlessfun (aref a left) (aref a right)) left)
          (t right))))

(provide "heap")

;;; --- end of file ---















(defun mularZ2(x y);; x and y are sparse polys in arrays as produced by (make-racg n m). return array
  ;; hack for heap management with a starter entry in ans
  (let* ((xexps (car x))
	 (xcoefs (cdr x))
	 (yexps (car y))
	 (ycoefs (cdr y))
	 (yexp0 (aref yexps 0))
	 (ycoef0 (aref ycoefs 0))
	 (ylim (length (car y)))
	 (ans (make-array (1+ (length xexps)))))
	 
    (declare (optimize (speed 3)(safety 0))(type (simple-array t (*)) ans xexps xcoefs yexps ycoefs))
    (setf (aref ans 0)(cons 0 (cons 0 #'(lambda() nil)))) ; blocker for heap management
    (dotimes (i (length xexps) ans)
	    (setf (aref ans (1+ i))
	      (cons (+ yexp0 (aref xexps i)) ;; the lowest exponent
		  (cons (* (aref xcoefs i)ycoef0) ; its coefficient
		  (let ((index -1)	; index into y
			(xexpi (aref xexps i))
			(xco (aref xcoefs i)))
			#'(lambda()
			    (incf index)
			    (if (>= index ylim) nil
			    (values (+ xexpi (aref yexps index)) ; the exponent
				    (* xco (aref ycoefs index))) ; the coef
			    )))))))))


(defun mulstr2(x y);; buggy . spends 28% of time in percolate-down, 17% in lesser-child
  (let* ((ans nil)
	 (h (mularZ2 x y));; initial heap
	 (r nil)
	 (k nil)
	 (fill-pointer (1- (length h))))
    (declare (special fill-pointer))
;;	  (format t "~% h=~s" h)
    ;; 
    (push (cons 0 0) ans);; initial answer
    (loop while (not (heap-empty-p h)) do
	  (setf r (heap-remove-fast h))	; next term
;;	  (format t "~% r=~s" r)
	  (if (= (car r)(caar ans))
	      (incf (cdar ans)(cadr r)) ; add coefficient.
	    (push (cons (car r)(cadr r)) ans)); put in new term.
	  (if (setf k (update-st2 r))		
	      (heap-insert h k)		;update r and re-insert in heap
	    ))
    ;;(sorthash ans #'< :key #'car)
    (nreverse ans)))
		

(defun update-st2(w) ;; w looks like ( enumber  cnumber  function)
  (declare (special *ht*))
  (multiple-value-bind (e1 c1)
      (funcall (cddr w))

    ;; update the hash table, as a debugging tool..
    (cond (e1 ;(incf (gethash e1 *ht* 0) c1) ;; omitting this line makes it much faster
	      (setf (car w) e1 (cadr w) c1) w )
	  (t nil))))