;;; -*- Lisp -*-
;;;;; excerpt from newmul with Roman's version of chaining, sort of.
;;; BUGGY?

(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.
;; 
(defvar *x* nil)
(defvar *y* nil)
  
(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))))


(let ((*YYEXPS* nil)
      (*YYCOEFS* nil)
      (*YYLEN* 0)
      (fill-pointer 0))
  (declare (fixnum *YYLEN* fill-pointer)
	   (type (simple-array t (*)) *YYEXPS* *YYCOEFS*))

  
  (defun make-racg(n d c) 
    ;; make Random Array of increasing exponents of length n with
    ;; spacing k, with 0<k<=d between elements and with Coefficients that
    ;; are random integers between 1 and c. d may be, in General, arbitrary precision.
    (let ((exps (make-array n))
	  (coefs(make-array n))
	  (v (random d)))
      (dotimes (i n (cons exps coefs))
	(setf (aref coefs i)
	      (+ 1 (random c)))
	(setf (aref exps i)v)
	(setf v (+ 1 v (random d))))))


  (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-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 (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 must have the first element duplicated.
  ;; e.g.  (3 5 8)  should be (3 3 5 8)
  
  
  ;;;THIS IS WRONG?
  
  ;; try again.
  
  (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)
      (do ((i fp parent)
	   (parent (ash fp -1);; initialize to parent of fp
		   (ash parent -1);; and its parent, next time
		   ))
	  ((or (< parent 1) (> (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.
	 ((= (st-exp xx) (st-exp c));; aha, we can insert this item by chaining here.
	  (incf (st-coef c) (st-coef xx))
;;	  (format t "~% insert4 updated c to ~s " c)
	   (decf fill-pointer) ;; we don't need the extra spot for this item ?????
	  (if (setf xx (next-term xx))	;try inserting the successor to this item then.
	      (heap-insert4   a xx ));; is this right?
	  (return a)	  )
					;exit from do ; the new heap.
	 (t;; swap locations i and parent. Should hack this to be faster.
	  (setf (aref a i) c
		(aref a parent) xx
		))
	 ))))
	 
	
  
  (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)))
  ;;  (format t "~% before remove, heap=~s, ~%fillpointer ~s" a fill-pointer)
    (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-down a 1 last-object fp)) last-object)
	;;  	  (setf (aref a (percolate-down a 1 last-object fp)) last-object)
;;	  (format t "~% after remove, heap=~s, ~%fillpointer ~s" a fill-pointer)
	  xx)))


 (defun mul-heap4(x y);; look for "chain" opportunity.  ;; clean up program.
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ4 x y)(mularZ4 y x)))
	   (alen  10)			; heuristic(* (length (car x))(length (car y))))
	   ;; 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
;;	    (format t "~% newe= ~s newc= ~s heap=~s" newe newc hh)
;;	    (format t "~% anse= ~s ansc= ~s " anse ansc)
	;;    (format t "~% preve= ~s prevc= ~s " preve prevc)
	
	      (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?
		;(setf top (next-term (copy-st top))) ;; is there a successor?
			(heap-insert4 hh top)) ;; if so, insert it in heap
	    
	   ;; (heap-insert4 hh top)
	    )
      ;; 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)
		     (* (the integer (st-v rec))(the integer (aref *YYCOEFS* i))))
	       (setf (st-exp rec) 
		     (+ (the integer (st-e rec))(the integer (aref *YYEXPS* i)))
		     )
	       rec))))
 

 (defun testu()(mul-heap4 u u)))

(defun test4()(time (mul-heap4 *x* *y*)))

(defvar *r* nil)
(defvar *s2* nil)
(defvar *s2p1* nil)

(defun setup()
  (let ((k1 nil)(k2 nil))
    (setf *r* (cons #(0 1 100 10000 1000000)#(1 1 1 1 1))); 1+x+y+z+t
    (setf *s2* (power *r* 20))
    (setf k1 (make-array 10626 :initial-contents (car *s2*)))
    (setf k2 (make-array 10626 :initial-contents (cdr *s2*)))
    (setf (aref k2 0) 2)
    (setf *s2p1* (cons k1 k2))))




(defun setup2() ;; smaller problem
  (let ((k1 nil)(k2 nil)(r2 nil)(len 0))
    (setf r2 (cons #(0 1 100 10000 )#(1 1 1  1))) ; 1+x+y+z   ;; not +t
    (setf *s2* (power r2 20))
    (setf len (length(car *s2*)))
    (setf k1 (make-array len :initial-contents (car *s2*)))
    (setf k2 (make-array len :initial-contents (cdr *s2*)))
    (setf (aref k2 0) 2)
    (setf *s2p1* (cons k1 k2))))

(defun power(p n)(if (eql n 1) p(mul-naive p (power p (1- n)))))
(defun power2(p n)(if (eql n 1) p(mul-heap p (power2 p (1- n)))))

(defun arr=(r s)
  (and (= (length r)(length s))
       (dotimes (i (length r) t)
	 (if (/= (aref r i)(aref s i)) (return nil)))))

(defun pol=(r s)(and (arr= (car r)(car s))
		     (arr= (cdr r)(cdr s))))

(defun diff1(r s)
  (dotimes (i (min (length s)(length r))  'nodiff)
    (cond ((= (aref r i)(aref s i))nil)
	  (t (return (format nil "differs at ~s: ~s ~s" i (aref r i)(aref s i)))))))