;;; multiplication of polyns based on sortx+y.lisp
;;;
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; sort X+Y various ways.
(defun make-rlc(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 ((res nil) (v (random d)))
    (dotimes (i n (nreverse res))
      (push (cons v (1+ (random d))) res)
      (setf v (+ 1 v (random d))))))

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

(defun test(n m)
  (setf x (make-rlc n 5))
  (setf y (make-rlc m 5))
  (format t "~% heap time")
  (time (mulheap x y))  
  (format t "~% hash time")
    (time (mulhash x y)))

;; works.
(defun make-tconc(r)(let ((h (list r)))(cons h h)))

(defmacro tconc(tc new) ;; changes tconc structure.
  `(let()(setf (cdr ,tc) (setf (cddr ,tc)(list ,new)))
  ,tc))

(defun mulheap (x y);; will not be normalized [could have 0 coeffs]
  (let*  ((x1 (caar x))
	  (m (maplist #'(lambda (r)(cons (+ (caar r) x1) (cons r x)))
		      (cons (car y) y)));; add one element to heap initialization this way..
	  (h (create-heap-ordered m))
	  (ans (make-tconc (list -1)))	;stopper removed later
	  (r nil)
	  (alist nil)
	  (nexti nil)
	  (nextj nil))
    (declare (optimize (speed 3)(safety 0)))
    (loop
     (if (heap-empty-p h)(return ans))
     (setf r (heap-remove-fast h))
      ;;   (format t "~%removing r=~s" r)
      ;;   (format t "~% (car r)=~s (caar  ans)=~s" (car r) (caar ans))
   ;;  (format t "~% newterm= ~s"   (* (cdaar(cdr r))(cdaddr r)))
     (setf alist (cdr r))
     (cond ((/= (car r)(caadr ans));; same as previous exponent, if not,
	    ;; put in a new exponent-coefficient pair on ans.
	    ;;    (format t "~% multiplying ~s by ~s to set into term ~s"  (cdaar(cdr r))(cdaddr r) r)
	    (setf (cdr r)(* (cdaar(cdr r))(cdaddr r)))
	    ;; otherwise put this one on the list for the answer
	    ;; we should check here if previous item was coefficient 0, since we
	    ;; have all contributors added in. If 0, remove it.
	    (tconc ans r))		;whew
	   ;; just compute a1*b1+a2*b2+...
	   (t (incf (cdadr ans) (* (cdaar(cdr r))(cdaddr r))))) 
      ;;     (format t "~%ans=~s" ans)
      ;;     (format t "~%alist=~s" alist)
     (cond (alist  
	    (setf nexti (car alist))
	    (setf nextj (cddr alist))
	    (setf (cdr alist)(cddr alist))
	  
	    ;; (format t "~% nexti= ~s ~%nextj= ~s" nexti nextj)
	    ;; (format t "~% nextexp=~s"  (+ (caar nexti)(caar nextj)))
	    (cond ((and nexti nextj)
		   (heap-insert h (cons (+ (caar nexti)(caar nextj))
					;(cons nexti nextj)
					alist))))	 
	    ))) (cdr (car ans))))	; in right order, without reversing


(defun mulhash (x y)  ;; will not be normalized [could have 0 coeffs]
  (let ((ans (make-hash-table)))
    (dolist (i x (sorthash ans))
      (dolist(j y)
	(incf (gethash (+ (car i)(car j)) ans 0)
	      (* (cdr i)(cdr j)))))))
	  


;;; HEAP stuff
;;; REMOVING assumptions about data being fixnums. Just fixnum array indexes
;;; see heap.lisp for docs.
;;; here we make array big enough so adjustable is not needed.

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


(defmacro hlessfun (u v)
  `(let* ((r ,u)(s ,v)
	  (cr (car r))(cs (car s)))
     (declare (integer cr cs)
					;(fixnum cr cs)
	      )	
  (cond ((< cr cs) t)
	;;((= cr cs) (< (the integer (caaadr r))(the integer (caaadr s))))
	)))

#+ignore 
(defmacro hlessfun (u v)		;slower.
  `(let ((r ,u)(s ,v))
     (case (signum (- (the integer (car r))(the integer(car s))))
       (-1 t)
       (0 (< (the integer (caaadr r))(the integer (caaadr s))))
       (1 nil))))

(defun percolate-down (a hole x)
  "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 (*)) a)(optimize (speed 3)(safety 0)))
	
  (do ((child (lesser-child a hole) (lesser-child a hole)))
      ((or (>= child fill-pointer) (hlessfun x (aref a child)))
       hole)
    (declare (fixnum child fill-pointer))
      (setf (aref a hole) (aref a 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 (fixnum hole)  (type (simple-array t (*)) a)
	   (optimize (speed 3)(safety 0)))
    (setf (aref a 0) x)
    (do ((i hole parent)
	 (parent (ash hole -1) (ash parent -1) ))
	((not (hlessfun x (aref a parent))) i)
      (declare (fixnum i parent))
      (setf (aref a i) (aref a parent))))

(defun create-heap-ordered
  (initial-contents)
  ;; used only if initial-contents is already sorted list.
  ;; important: initial-contents must have the first element duplicated.
  ;; e.g.  (3 5 8)  should be (3 3 5 8)
 ;; (format t "*")
  
  "Initialize the indicated heap.  If INITIAL-CONTENTS is a non-empty
list, the heap's contents are initialized to the values in that
list; they are assumed already ordered according to hlessfun.  INITIAL-CONTENTS must
be a list or NIL."
  (let* ((n (length initial-contents))
	 (heap
	  (make-array  n
		       :initial-contents initial-contents
		       :element-type t)))
    (setf fill-pointer n)		;global value for this heap
    heap))

(defmacro heap-empty-p (heap)
   (declare (ignore heap))
  "Returns non-NIL if & only if the heap contains no items."
  `(= 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)

    (declare (type (simple-array t (*)) a)(fixnum fill-pointer)
	     (optimize (speed 3)(safety 0)))
         ;;(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)
 
  (declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
     ;;(format t "-")
  (let* ((x (aref a 1))
	 (last-object		
	  (progn (decf fill-pointer)(aref a fill-pointer) ) ))
    (setf (aref a (percolate-down a 1 last-object)) last-object)
    x))

(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."
  (declare (type (simple-array t (*)) a)
	   (optimize(speed 3)(safety 0)) (fixnum parent) )
  (let* ((left (ash parent 1))
         (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))))
;;; --- end of heap stuff ---

#|  test results
 heap time
; cpu time (non-gc) 859 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  859 msec user, 0 msec system
; real time  859 msec
; space allocation:
;  506,481 cons cells, 4,024 other bytes, 0 static bytes
 hash time
; cpu time (non-gc) 172 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  172 msec user, 0 msec system
; real time  172 msec
; space allocation:
;  31,287 cons cells, 240,752 other bytes, 0 static bytes
((8 . 8) (11 . 12) (12 . 14) (13 . 16) (14 . 4) (15 . 3) (16 . 7)
 (17 . 10) (18 . 20) (19 . 8) ...)

|#