;;; -*- Lisp -*-
;;; mulstream.lisp multiplication of polyns based on streams
;;; now exponents are fixnums and coeffs.  NOT really justifiable...
;;; not faster.

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

;; sort X+Y various ways.

(defun make-rac(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 :element-type 'fixnum))
	(coefs(make-array n :element-type 'double-float))
	;(coefs(make-array n :element-type 'fixnum))
	
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	;(+ 1 (random d))
	(+ 1.0d0 (random d))
	)
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))


(defun test(n m)
  (setf x (make-rac n 5))
  (setf y (make-rac m 5))
;  (format t "~% mulhashb time-ar")
;    (time (mulhashb x y))
  (format t "~% mulhashe time-ar")
    (time (mulhashe x y)))

;; multiply two polynomials via a hash table. For comparison
(defparameter *zcoefs-lim* 50)
(defparameter *zcoefs* (make-array *zcoefs-lim* :element-type 'double-float))

(defun make-racf(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 :element-type 'fixnum))
	;(coefs(make-array n))
	(coefs(make-array n :element-type 'fixnum))
	
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	;(+ 1 (random d))
	(+ 1 (random d))
	)
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))

(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))
	;(coefs(make-array n :element-type 'fixnum))
	
	(v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
	;(+ 1 (random d))
	(+ 1 (random d))
	)
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))


;;; see heap.lisp for docs.

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

;; heap is just an array.

;;(defstruct heap a  max-count)

(defmacro hlessfun (a b) `(< (the fixnum (car ,a))(the fixnum (car ,b))))
;(defmacro hlessfun (a b) `(< (car ,a) (car ,b)))

(defmacro heap-empty-p (heap)
  "Returns non-NIL if & only if the heap contains no items."
  `(<= (the fixnum 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))."
  (declare (optimize(speed 3)(safety 0))(type (simple-array t (*)) a) (fixnum parent fp))
  `(let* ((left (ash ,parent 1))
         (right (1+ left)))
    (declare (fixnum left right))
    (cond ((>= left ,fp) ,fp)
          ((= right ,fp) left)
          ((hlessfun (aref ,a left) (aref ,a right)) left)
          (t right))))

(defun percolate-down (aa hole x fp)
  "Private. Move the HOLE down until it's in a location suitable for X.
Return the new index of the hole."
    (declare (fixnum hole fp)(type (simple-array t (*)) aa)(optimize (speed 3)(safety 0)))
	
  (do ((child (lesser-child aa hole fp) (lesser-child aa hole fp)))
      ((or (>= (the fixnum child) fp) (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
  
    ;; Append a hole for the new element.
    ;; (vector-push-extend nil a) 
    ;; assume enough room...
    ;; (vector-push nil a)
   
   ;; (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)
 ;;   (format t "-")
  (declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
  (let*  ((x (aref a 1))
	 (last-object	
	  (progn (decf fill-pointer)(aref a fill-pointer) )))
    ;; could declare array a element-type if we knew it if was fixnum
    ;; here and elsewhere.

    (setf (aref a (percolate-down a 1 last-object fill-pointer)) last-object)
	     x))

(defun heap-peek (heap)
  "Return the first element in the heap, but don't remove it.  It'll
 be Erroneous if the heap is empty.   (Should that be an error?)"
  (aref heap  1))



(provide "heap")


(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))(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()
			      (declare (fixnum index xco expi)) ;;;****
			      (incf index)
			    (if (>= index ylim) nil
			    (values (+ xexpi (the fixnum(aref yexps index))) ; the exponent
				    (* xco (the fixnum (aref ycoefs index)))) ; the coef
			    )))))))))


(defun mulstr2(x y);; buggy . spends 28% of time in percolate-down, 17% in lesser-child
  (let* ((ans nil)
	 (h (mularZ2 x y));; initial heap
	 (r nil)
	 (fill-pointer  (length h))
	 )
	 
    (declare (special fill-pointer))
;;	  (format t "~% h=~s" h)
    (push (cons 0 0) ans);; initial answer
    (loop while (not (heap-empty-p h)) do
	  (setf r (heap-remove-fast h))	; next term
;	  (format t "~% ans=~s, r=~s" ans r)
	  (if (= (car r)(caar ans))
	      (incf (cdar ans)(cadr r)) ; add coefficient.
	    (push (cons (car r)(cadr r)) ans)); put in new term.
	  (multiple-value-bind (e1 c1)
	      (funcall (cddr r))
	    (cond (e1 
		   (setf (car r) e1 (cadr r) c1)
	;	   (format t "~% new r is ~s" r)
		   (heap-insert h r)))))
    ;;update r and re-insert in heap
    ;; could put in right order by using tconc.
   ;(cdr  (nreverse ans))
    ;;ans
    ))