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

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

;; sort X+Y various ways.

(defun test(n m)
  (setf x (make-racg n 5))
  (setf y (make-racg m 5))
  (format t "~% mulhashg time")
  (time (mulhashg x y))
  (format t "~% mulstr3 time")
  (time (mulstr3 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.

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

;; a heap is just an array.


(defmacro heap-empty-p (heap)
  "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)
  "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))."
  (declare (type (simple-array t (*)) a) (fixnum parent fill-pointer))
  `(let* ((left (ash ,parent 1))
         (right (1+ left)))
    (declare (fixnum left right))
    (cond ((>= left fill-pointer) fill-pointer)
          ((= right fill-pointer) left)
          ((hlessfun (aref ,a left) (aref ,a right)) left)
          (t right))))

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

(let ((fill-pointer 0))

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

;;; 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*  ((x (aref a 1))
	   (last-object (aref a (decf fill-pointer))))
    (setf (aref a (the fixnum (percolate-down a 1 last-object))) last-object)
    x))

(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
			    )))))))))



#|
[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 ...))
|#

(defun mulstr3(x y)		;;works.  answer in array;;
  (let* ((h (mularZ2 x y));; initial heap
	 (r nil)
	 (anse (make-array 10 :adjustable t :fill-pointer 0))
	 (ansc (make-array 10 :adjustable t :fill-pointer 0))
	 (preve 0)
	 (prevc 0))
;;    (declare (special fill-pointer))
    (setf fill-pointer  (length h))

    (loop while (not (heap-empty-p h)) 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))
		  )))			;re-insert in the heap
    (vector-push-extend prevc ansc) ;put in last results before exit
    (vector-push-extend preve anse)
    (cons anse ansc)))


;; the idea here is to burp out one exponent-coefficient pair per call.
;; buggy.
(defun mulstream(x y)
  (let* ((h (mularZ2 x y));; initial heap
	 (r (list 0 0))
	 (preve 0)
	 (prevc 0)
	 (done nil))
 ;;    (declare (special fill-pointer))
     (setf fill-pointer  (length h))
   ; (format t "~%heap= ~s~% fill-pointer= ~s" h fill-pointer)
    #'(lambda ()
		;(format t "~%fill-pointer=~s" fill-pointer)
    (loop while (not (heap-empty-p h)) do
	  (setf r (heap-remove h))	; take next term off heap
	  ; (format t "~%r= ~s" r)
	 
	  (cond ((= preve (car r))	; is this same exponent as previous?
		 (incf prevc (cadr r))); if so, add coefficient.
		(t 
		 (setf prevc (cadr r))
		 (setf preve (car r))

	  (multiple-value-bind (e1 c1) ;; advance the stream taken from heap

			      (funcall (cddr r))
			    (cond (e1	; if 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))
				  ))
			  
			  (unless (= prevc 0) ; normalize: remove zero coeffs
			    (return  (values preve prevc)))))
					; initialize new coefficient
)			;re-insert in the heap

    (if (not done)(progn (setf done t) (values preve prevc))
      nil))))



)