;;; -*- 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
    ;; 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)
    (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 (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 :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) 
	       (< 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. 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)))

      (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 next-term-peek(rec)   	; record is a stream 5-tuple  return next expon 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  (+ (st-e rec)(aref *YYEXPS* i))))))
 
 
 
 
 )