;;; -*- Lisp -*-
;;;;; excerpt from newmul with Roman's version of chaining, sort of.
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))
;; set up a lexically available context for handling several
;; variables which would otherwise be global, and require more
;; effort to access from compiled code.
  
(defmacro heap-empty-p ()
    "Returns non-NIL if & only if the heap contains no items."
    `(<= fill-pointer 1))

(defmacro hlessfun (a b) `(< (st-exp ,a) (st-exp ,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))))

;; 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 5 items
    ;; which by reference to arrays of XX and YY , can burp out
    ;; successive exponents and coefficients of the product.
  
    (exp 0 :type integer)
    (coef 0 :type integer) 
    (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) 
    (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 (fixnum i) (type (simple-array t (*))
				ans xexps xcoefs yexps ycoefs))
      (setf *YYEXPS* yexps *YYCOEFS* ycoefs *YYLEN* ylen)
      (setf (aref ans 0) (make-st :exp 0 :coef 0 :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) :coef (* 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) (< xxval (the fixnum(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. 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)
	       (optimize (speed 3)(safety 0)))

      ;;      (format t "~%in insert4, fill-pointer is ~s" fill-pointer)
      ;;   (setf (aref a (decf fp)) xx)
      (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))
	;;	    (format t "~%insert4: c=~s parent=~s xx=~s"c parent xx)
	(setf c (aref a parent))
	(cond
	 ;; #+ignore
	 ;; program is right if this clause is ignored. But not as fast.
	 ((= (st-exp xx) (st-exp c));; aha, we can insert this item by chaining here.
					;(format t "~% insert4 updated ~s  " c)
					; (format t "~% heap is ~s" a)
	  (incf (st-coef c) (st-coef xx))
					;(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
 (defun mul-heap4(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)(newc 0) )
      (declare (integer preve prevc newe newc) 
	       (type (array t(*)) anse ansc hh))
      
      (setf fill-pointer  (length hh))
      (loop while  (setf top (heap-remove4 hh))   do
	    (setf newe (st-exp top) newc (st-coef top))
	      ;; remove top item from heap
	      (cond 
	       ((= preve newe)		; is this same exponent as previous?
		(incf prevc newc))	; 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 newc)	; 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-exp rec) 
		     (+ (st-e rec)(aref *YYEXPS* i)))
	       rec)))))