;;; -*- Lisp -*-
;;; mulstream.lisp multiplication of polyns based on streams
;;; general, works, about 2.3X slower than mulhashg, but uses less storage.
;;;


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

;; sort X+Y various ways.
(defvar x nil)
(defvar y nil)

(defun test(n m)
  (setf x (make-racg n 5))
  (setf y (make-racg m 5))
  (format t "~% mulhashg time - best hash table version")
  (time (mulhashg x y))
  (format t "~% mulstr3 time - use a heap in an array, insert from bottom")
  (time (mulstr3 x y))
  (format t "~% mulstr6 time - use a heap in an array, insert from top sometimes")
  (time (mulstr6 x y))  
  (format t "~% mul9 time - use streams, binary merging")
  (time (mul9 x y))  

   (format t "~% mul12 time - just use insertion into single list")
  (time (mul12 x y))  
)

(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))
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	(+ 1 (random d)))
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))

;;; see heap.lisp for docs.

;; a heap is just an array.

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

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

(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)
          ((hlessfun (aref ,a left) (aref ,a right)) left)
          (t right))))


;; construct length(x) streams, a_0*y,  a_1*y, .... 

(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))(fixnum j)(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 0)
			      (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
			    )))))))))




;; set up a context for heap array fill pointer..

(let ((fill-pointer 0)) 

 #+ignore ;; maybe clearer what is going on
  (defun percolate-down (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."
    (declare (fixnum hole fp)
	     (type (simple-array t (*)) aa)(optimize (speed 3)(safety 0)))
	
    (do ((child (lesser-child aa hole fp) (lesser-child aa hole fp)))
	((or (>= (the fixnum child) fp) (hlessfun xx (aref aa child)))
	 hole)
      (declare (fixnum child))
      (setf (aref aa hole) (aref aa child)
	    hole child)))
  
  
;; #+ignore;; works also, maybe a teensy bit faster.
  (defun percolate-down (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 (car xx)))

      (declare (fixnum hole fp)(integer xxval aaval)
	       (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) (< xxval (the fixnum(car aaval))))
	   hole)
	(declare (fixnum child))
	(setf (aref aa hole)  aaval  hole child))))

  (defun percolate-up (a hole xx)
    "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) xx)
    (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 xx (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-insert (a xx fp)
    "Insert a new element into the heap.  Returns the heap.";; rjf
    (declare (type (simple-array t (*)) a)(fixnum fill-pointer)
	     (optimize (speed 3)(safety 0)))
    ;; (format t ".")
    (setf (aref a fp) nil)
    (incf fill-pointer)			; the global one
    ;; 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 fp xx)) xx))

;;; here's a trick to try, if we believe Roman Pearce's comment:
;;; Instead of percolating down, see what is going to be inserted next,
;;; and try to insert it at the top. This saves a lot of time if it works.

  (defun heap-remove(a) 
    ;; assumes non-empty!! if empty next answer is bogus.
    ;; answer after that is an error (fill pointer can't pop)
    ;;   (format t "-")
    (declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
    (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 (the fixnum (percolate-down a 1 last-object fp))) last-object)
      xx))


  #|
[2c] cl-user(931): (time (mulhashg xxx yyy))
; cpu time (non-gc) 94 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  94 msec user, 0 msec system
; real time  93 msec
; space allocation:
;  4,514 cons cells, 511,264 other bytes, 0 static bytes
(#(4 6 7 9 10 11 12 13 14 15 ...) . #(9 6 9 10 2 19 2 3 10 5 ...))

;;; another... with array answers..

[5] cl-user(1002):  (time (mulstr3 xxx yyy))
; cpu time (non-gc) 219 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  219 msec user, 0 msec system
; real time  297 msec
; space allocation:
;  2,135 cons cells, 142,856 other bytes, 0 static bytes
(#(4 6 7 9 10 11 12 13 14 15 ...) . #(9 6 9 10 2 19 2 3 10 5 ...))

;; even better, not adjustable arrays, making sure 2nd arg is longer than first
;; notice storage is way down.
(time (mulstr3 x y))
; cpu time (non-gc) 235 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  235 msec user, 0 msec system
; real time  235 msec
; space allocation:
;  1,037 cons cells, 114,344 other bytes, 0 static bytes
(#(1 3 6 7 8 11 12 13 14 15 ...) . #(20 4 25 16 5 36 25 4 1 8 ...))



|#


;;

(defun mulstr3(x y)  ;;works
  
  (let* ((xl (length (car x)))
	 (yl (length (car y)))
	 (h (if (> (the fixnum xl)(the fixnum yl))(mularZ2 y x)(mularZ2 x y)));; initial heap
	 (ml (+  (the fixnum xl)(the fixnum yl)))
	 (r nil)
	 (anse (make-array ml :adjustable t :fill-pointer 0))
	 (ansc (make-array ml :adjustable t :fill-pointer 0))
	 (preve 0)
	 (prevc 0))
    (declare (integer preve prevc ) (fixnum xl yl ml))
  ;  (declare (special fill-pointer))
    (setf fill-pointer  (length h))

    (loop while (not (heap-empty-p)) do
	  (setf r (heap-remove h))	; take next term off heap
	  (cond ((= preve (car r))	; is this same exponent as previous?
		 (incf prevc (cadr r)))	; if so, add coefficient.
		(t (unless (= prevc 0)	; normalize: remove zero coeffs
		     (vector-push-extend prevc ansc)
		     (vector-push-extend preve anse))
		   (setf prevc (cadr r))
		   (setf preve (car r)))); initialize new coefficient
	  (multiple-value-bind (e1 c1);; advance the stream taken from heap
	      (funcall (cddr r))
	    (cond (e1			; check: is there more on this stream--
		   (setf (car r) e1 (cadr r) c1); update the stream
					;     	  (format t "~% inserting r=~s" r)
		   (heap-insert h r fill-pointer ))
		  )))			;re-insert in the heap
    (vector-push-extend prevc ansc)		;put in last results before exit
    (vector-push-extend preve anse)
    (cons anse ansc)))


;; well, this part works if update-heap works.

  (defun mulstr6(x y);; hack to make heap management cuter.
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ2 x y)(mularZ2 y x)))
        ;; initial heap
	   (anse (make-array 10 :adjustable t :fill-pointer 0))
	   (ansc (make-array 10 :adjustable t :fill-pointer 0))
	   (preve 0)
	   (prevc 0))
      (declare (integer preve prevc))
      (setf fill-pointer  (length hh))
      (loop while (not (heap-empty-p)) do
	    (multiple-value-bind
		(newe newc)
		(update-heap hh)
	      (declare (integer newe newc))
					; take next term off heap and insert its successor
	      (cond ((= preve newe)	; is this same exponent as previous?
		     (incf prevc newc))	; if so, add coefficient.
		    (t (unless (= prevc 0); normalize: remove zero coeffs
			 (vector-push-extend prevc ansc)
			 (vector-push-extend preve anse))
		       (setf prevc newc)
		       (setf preve newe))); initialize new coefficient
	      ))
      (vector-push-extend prevc ansc)	;put in last results before exit
      (vector-push-extend preve anse)
      (cons anse ansc)))
      
  #+ignore
  (defun update-heap (hh)
    ;;return top of heap h, which is necessarily not empty.
        ;; take successor in stream of top of heap, and insert it in heap.
    (declare (optimize (speed 3)(safety 0)))
    (let* ((r (heap-remove hh))
	   (newe (car r))
	   (newc (cadr r)))
      (bubble-up hh r)
      (values newe newc))  )

            ;; try next more elaborate version

	    	    ;; this next guy also WORKS

  ;; WORKS
  (defun update-heap (h)
    ;;return top of heap h, which is necessarily not empty.
        ;; take successor in stream of top of heap, and insert it in heap.
    (declare(type (simple-array t(*)) h)(optimize (speed 3)(safety 0)))
    (let* ((top (aref h 1))
	   (left (aref h 2))
	   (right (aref h 3))
	   (lefti (car left))
	   (righti (car right))
	   (rprog (cddr top))
	   (ta 0)(tb 0)
	   ;(t1 nil)
	   (t2 nil))
					; (format t "~%top= ~s ~%left=~s ~%right=~s"top left right)
      (multiple-value-bind (newe newc)
	  (funcall (the compiled-function rprog))
	     ;; case, newe <= lefti  and newe <= righti.  put it at the top.
	           ;; avoiding all heap programs
		        ;; top of heap is location 1  (ignore 0)
			   ;; left branch is location 2
			         ;; right branch is location 3

	(cond((null newe) 
					; (format t "#")
	      (setf ta (car top) tb (cadr top))
	      (heap-remove h)
	      (values ta tb))

	     ((and			;(setf t1 )
	       (<= newe lefti)
	       (setf t2 (<= newe righti)))
	      	    ;; insert this node at the top of the heap
		     	    ;; and return it.
			      	    ;; fill-pointer remains the same.
					;		      (format t "^")
	      (setf ta (car top)tb (cadr top))
	      (setf (car top) newe (cadr top) newc)
	      (values ta tb))
	         	   ;; hardly ever happens to be true, but
			     	   ;; this clause is WRONG, somehow.
	     #+ignore  
	     (t1 
	      ;; but  NOT (newe <= righti). 
	      ;; therefore  righti < newe <= lefti
	      (format t "~% ~s <~s<= ~s" righti newe lefti)
	      (setf ta (car top)tb (cadr top))
					;promote the right node to top
	      (setf (aref h 1) right) 
					;put newe on right
	      (setf (car top) newe (cadr top) newc)
	      (setf (aref h 3)top)
	      (format t "~%    ta=~s tb=~s" ta tb)
	      (values ta tb) )
	     (t2 
					;		   (format t ">")
	      (setf ta (car top)tb (cadr top))
					;promote the left node to top
	      (setf (aref h 1) left) 
					;put newe on left
	      (setf (car top) newe (cadr top) newc)
	      (setf (aref h 2)top)
	      (values ta tb) )
	     (t

	      	;; the new item doesn't fit in the top 3.
					;	   (format t "-")
	      (let ((r (heap-remove h)))
		(setf ta (car top)tb (cadr top))
		(setf (car r) newe (cadr r) newc)
		(heap-insert h r fill-pointer)
	       (values ta tb))
	      )))))


;; the idea here is to burp out one exponent-coefficient pair per call.
;; it works.

  (defun mulstr5(x y) ;; try to make a stream output 
  
  (let* ((xl (length (car x)))
	 (yl (length (car y)))
	 (h (if (> (the fixnum xl)(the fixnum yl))(mularZ2 y x)(mularZ2 x y)));; initial heap
	 (r nil)
	 (anse 0)
	 (ansc 0)
	 (preve 0)
	 (prevc 0))
    (setf fill-pointer  (length h))
   ; (format t "~%0h=~s" h)
    #'(lambda() (prog nil
		 again
		  (loop
		    (cond ((heap-empty-p)(setf anse preve ansc prevc)
					 (setf preve nil prevc nil)
					 (return t))) ;exit from loop
		    ;; collect all terms with exponent preve
		    ;;  peek at next item
		    (cond ((not (equal preve (car (aref h 1))))
			   (setf anse preve ansc prevc)
			   (setf r (heap-remove h))
			   (setf preve (car r) prevc (cadr r))
			   (bubble-up h r)
			   (return t))
			  (t		; preve= prevc
			   (setf r (heap-remove h))
			   (incf prevc (cadr r))
			   (bubble-up h r)))) ;keep looping
		  (if (equal ansc 0)(go again) ; don't return terms with 0 coef
		    (return (values anse ansc)))))))	     

  (defun bubble-up(hh r)
    (declare (optimize (speed 3)(safety 0)))
        (multiple-value-bind (e1 c1);; advance the stream taken from heap
	    (funcall (the compiled-function (cddr r)))
	  
	  (cond (e1			; check: is there more on this stream--
		 (setf (car r) e1 (cadr r) c1) ; update the stream
		 (heap-insert hh r fill-pointer))
		)))

  #+ignore
  (defun testheap(u)
    (setf h (mularZ2 u u))
    (setf fill-pointer (length h))
    (format t  "~%h=~s" h)
    (loop while (not (heap-empty-p)) do
	  (setf r (heap-remove h)) 
	  (format t "~%r=~s~%h=~s" r h)))
  
)     ;; end of let


(defun muls (x y)(let ((k (mulstr5 x y))
			(anse nil)
			(ansc nil))
		    (loop (multiple-value-bind (e c) (funcall k)
		      (cond (e (push e anse)(push c ansc))
			    (t (return (cons (nreverse anse)(nreverse ansc)))))))))
;; about as fast as muls, but returns answer as pair of arrays.
(defun mula (x y)(let ((k (mulstr5 x y))
		       (anse (make-array 10 :adjustable t :fill-pointer 0))
		       (ansc (make-array 10 :adjustable t :fill-pointer 0)))
		    (loop (multiple-value-bind (e c) (funcall k)
			    (cond (e (vector-push-extend e anse)
				     (vector-push-extend c ansc))
				  (t (return (cons anse ansc))))))))




;; for comparison


(defparameter *zcoefs-limg* 50)
(defparameter *zcoefsg* (make-array *zcoefs-limg*))
(defun sorthashc(h ar) ;; utility for mulhashg
  (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) ;; normalize! remove 0 coefs
			 (vector-push-extend (cons i v) ans)))
	     h)
    (sort ans #'< :key #'car)))
(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))))))))))

;;; this program, mul7 works too, but is 10X slower than mulhashg
;;; unless you are allowing the representation of the answer as
;;; a stream, namely that produced by mulstr7. In which case the
;;; answer is instantaneous, pretty much.


(defun mergestreams(r s)
  ;; r and s each look like (exp coef . prog)
  (let ((er (car r))
	(es (car s)))
    (cond ((null er) s)
	  ((null es) r)
	  ((< er es)(cons er (cons (cadr r)
				   #'(lambda()(mergestreams (funcall (cddr r)) s)))))
	  ((< es er)(cons es (cons (cadr s)
				   #'(lambda()(mergestreams (funcall (cddr s)) r)))))
	  (t ;; equal!
	   (cons es (cons (+ (cadr s)(cadr r))
			  #'(lambda()(mergestreams (funcall (cddr s))(funcall (cddr r))))))))))

;; less consing, but buggy
(defun mergestreams2(r s)
  ;; r and s each look like (exp coef . prog)
  (let ((er (car r))
	(es (car s))
	(p nil) (q nil))
    (cond ((null er) s)
	  ((null es) r)
	  ((< er es)(setf p (cddr r))
	   (setf (cddr r) #'(lambda()(mergestreams2 (funcall p) s)))
	   r)
	  ((< es er)(setf q (cddr s))
	   (setf (cddr s) #'(lambda()(mergestreams2 (funcall q) r))) s)
	  (t;; equal!
	   (setf (cadr r)(+(cadr r)(cadr s)))
	   (setf q (cddr s) p (cddr r))
	   (setf (cddr r) #'(lambda()(mergestreams2 (funcall p)(funcall q))))
	   r))))


(defun mulstr7(x y) ;just merge pieces.
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ3 x y)(mularZ3 y x)))
	   (hl (length hh)))
	   (do ((ans (aref hh 1)(mergestreams (aref hh i) ans))
		(i 2 (1+ i)))
	       ((= i hl) ans))))

(defun mulstr7x(x y) ;just merge pieces.
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ3 x y)(mularZ3 y x)))
	   (hl (length hh)))
	   (do ((ans (aref hh 1)(mergestreams2 (aref hh i) ans))
		(i 2 (1+ i)))
	     ((= i hl) ans))))


(defun mularZ3(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))(fixnum j)(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 0)
			      (xexpi (aref xexps i))
			      (xco (aref xcoefs i)))
			  (labels((theprog()
				    (incf index)
				    (if (>= index ylim) nil
				      (cons (+ xexpi (aref yexps index)) ; the exponent
					    (cons  (* xco (aref ycoefs index)) ;the coef
						   #'theprog)))))
			    
			    #'theprog))))))))

(defun mul7(r s)
  (let ((k (mulstr7 r s))
	(ans nil))
    (loop  while k do (push (cons (car k)(cadr k)) ans)
		      (setf k (funcall (cddr k))))
    (nreverse ans)))

(defun mul7x(r s)
  (let ((k (mulstr7x r s))
	(ans nil))
    (loop  while k do (push (cons (car k)(cadr k)) ans)
		      (setf k (funcall (cddr k))))
     (nreverse ans)))


(defun mulstr8(x y) ;just merge pieces, but in a tree form. 49% faster than mul7
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ3 x y)(mularZ3 y x)))
	   (hl (length hh))
	   (stack (list (aref hh 1))))
      	   (do ((i 2 (1+ i)))
	       ((= i hl)nil)
	     (push (aref hh i) stack))
	   (while (cdr stack)  ;; more than one
		  (setf stack (pairmerge stack)))
	   (car stack)))

(defun mulstr8x(x y) ;just merge pieces, but in a tree form. 49% faster than mul7
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ3 x y)(mularZ3 y x)))
	   (hl (length hh))
	   (stack (list (aref hh 1))))
      	   (do ((i 2 (1+ i)))
	       ((= i hl)nil)
	     (push (aref hh i) stack))
	   (while (cdr stack)  ;; more than one
		  (setf stack (pairmerge2 stack)))
	   (car stack)))

(defun pairmerge (s)
  (if (null (cdr s)) s 
	(cons (mergestreams (car s)(cadr s))
	      (pairmerge (cddr s)))))

(defun pairmerge2 (s)
  (if (null (cdr s)) s 
	(cons (mergestreams2 (car s)(cadr s))
	      (pairmerge2 (cddr s)))))

(defun mul8(r s)
  (let ((k (mulstr8 r s))
	(ans nil))
    (loop  while k do (push (cons (car k)(cadr k)) ans)
		      (setf k (funcall (cddr k))))
    (nreverse ans)))

(defun mul8x(r s) ;; less consing, faster than mul8 by 1/4. 2X faster than mul7
  (let ((k (mulstr8x r s))
	(ans nil))
    (loop  while k do (push (cons (car k)(cadr k)) ans)
		      (setf k (funcall (cddr k))))
     (nreverse ans)))



(defun mul9 (x y)(let ((k (mulstr8x x y)) ;; like mul8x with arrays
		       (anse (make-array 10 :adjustable t :fill-pointer 0))
		       (ansc (make-array 10 :adjustable t :fill-pointer 0)))
		    (loop 
		      (vector-push-extend (car k) anse)
		      (vector-push-extend (cadr k) ansc)
		      (if (not(setf k (funcall (cddr k))))
			  (return (cons anse ansc))))))


(defun mul10 (x y)
  (let* ((ans (list (cons 0  0)))
	 (ptr ans)
	 (e 0)(c 0)
	 (xl (length (car x)))
	 (yl (length (car y)))
	 (xe (car x))
	 (ye (car y))
	 (xc (cdr x))
	 (yc (cdr y)))
    (declare (optimize (speed 3)(safety 0))(fixnum xl yl)
	     (type (simple-array t (*)) xe ye xc yc))
	     
    (dotimes (i xl ans)
      (dotimes(j yl)
	;; multiply terms 
	(setf e (+ (aref xe i)(aref ye j))
	      c (* (aref xc i)(aref yc j)))
	(cond ((< e (caar ptr))		; insert from beginning of ans
	       (setf ans (addterm e c ans)))
	      (t (setf ptr (addterm e c ptr))))))))

(defun addterm(e c L)
  (let ((ahead (cdr L)))
    (cond ((null ahead)
	   (setf (cdr L)(list (cons e c))))
	  ((= e (caar ahead))
	   (incf (cdar ahead) c))
	  ((< e (caar ahead))
	   (setf (cdr L)(cons (cons e c) ahead)))
	  (t (addterm e c (cdr L))))
    L))

	  

(defun mul11 (x y) ;; like mul10 but with better fingers into the answer for insertion
  (let* ((ans (list (cons 0  0)))
	 (ptr0 ans)
	 (ptr1 ans)
	 (e 0)(c 0)
	 (xl (length (car x)))
	 (yl (length (car y)))
	 (xe (car x))
	 (ye (car y))
	 (xc (cdr x))
	 (yc (cdr y)))
    (declare (optimize (speed 3)(safety 0))(fixnum xl yl)
	     (type (simple-array t (*)) xe ye xc yc))
	     
    (dotimes (i xl (cdr ans))
      (dotimes(j yl)
	;; multiply terms 
	(setf e (+ (aref xe i)(aref ye j))
	      c (* (aref xc i)(aref yc j)))
	(cond ((< e (caar ptr1))
	       
	       (if (< e (caar ptr0))
		   (progn;		(format t "0")
		     (setf ptr1 (addterm11 e c ans)))
		 (progn	;		(format t "1")
		   (setf ptr1(addterm11 e c ptr0))))) ;used never
	      (t ;(format t "2")
		 (setf ptr1 (addterm11 e c ptr1)))))
      (setf ptr0 ptr1))))


(defun addterm11(e c L) 
  (let ((ahead (cdr L)))
    (cond ((null ahead)
	   (setf (cdr L)(list (cons e c))))
	  ((= e (caar ahead))
	   (incf (cdar ahead) c)
	   (cdr L))
	  ((< e (caar ahead))
	   (setf (cdr L)(cons (cons e c) ahead))
	   (cdr L))
	  (t (addterm11 e c (cdr L))))  ))


(defun mul12 (x y) ;; like mul10 but with better fingers into the answer for insertion
  (let* ((ans (list (cons 0 0)))
	 (ptr1 ans)
	 (e 0)(c 0)
	 (xl (length (car x)))
	 (yl (length (car y)))
	 (xe (car x))
	 (ye (car y))
	 (xc (cdr x))
	 (yc (cdr y)))
    (declare (optimize (speed 3)(safety 0))(fixnum xl yl)
	     (type (simple-array t (*)) xe ye xc yc))
	     
    (dotimes (i xl (breakupL (cdr ans)))
      (dotimes(j yl)
	;; multiply terms 
	(setf e (+ (aref xe i)(aref ye j))
	      c (* (aref xc i)(aref yc j)))
	(cond ((< e (caar ptr1))
	        (setf ptr1 (addterm11 e c ans)))
	      (t(setf ptr1 (addterm11 e c ptr1))))) )))

(defun breakupL (L) 
  ;; L is a list of (exp . coef).
  ;; separate them into a pair: one array of exponents and one of coefficients.
  (let* ((ll (length L))
	 (exps (make-array ll))
	 (coefs (make-array ll))
	 (count 0))
    (declare (optimize (speed 3)(safety 0))(fixnum ll count)
	     (type (simple-array t (*)) exps coefs))
    (dolist (i L (cons exps coefs))
      
      (setf (aref exps count) (car i))
      (setf (aref coefs count)(cdr i))
      (incf count))))