;;; -*- Lisp -*-
;;; msort5.lisp multiplication of polyns based on
;;; sortx+y.lisp msort2-5 uses heaps implemented as lisp lists.  seems
;;; like a loser, 4X slower than hash, somewhat slower than arrays.
;;; msort5 uses lisp lists to construct a heap, each internal node n
;;; is, in this design (here . (left . right)) where here is a leaf
;;; node, and left/right may either be nil or an internal node or a
;;; leaf node.  A leaf node looks like a dotted pair consisting of
;;; (number . number).  They represent an exponent . coefficient pair.
;;; An empty heap looks like this: (nil)


(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 test(n m)
  (setf x (make-rlc n 5))
  (setf y (make-rlc m 5))
  (format t "~% heap time")
  (time (mulheap2 x y))  
  (format t "~% hash time")
    (time (mulhash x y)))

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


;; multiply two polynomials via a hash table. For comparison
(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)))))))

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

;;; HEAP stuff from heap-lisp2.
;;; see heap.lisp for docs.

(defmacro hlessfun (u v) `(< (car ,u) (car ,v)))

;;; implementation of heap data structure in lists.
;;; heapl is heap using lists, not arrays as traditional
;;; this version has leaf nodes as (number . stuff)

(defmacro leafp(x) `(numberp (car ,x)))

(defun create-heapl (ic)		;initial contents.  ic  should be ordered.
  (labels 
      ((split(z);; used to balance the initial contents.
	     (cond ((null z)(values nil nil))
		   (t  (let ((even z)
			     (odd (cdr z)))
			 (do ((i even (cdr i))
			      (j odd (cdr j)))
			     ((null j)(values even odd))
			   (setf (cdr i) (cddr i))
			   (setf (cdr j) (cddr j)))))))
	   
       (ch1 (ic)
	 (cond ((null ic) nil)
	       ((null(cdr ic)) (car ic)) ;final leaf node
	       (t (multiple-value-bind (a b)(split (cdr ic))
		    (cons (car ic)(cons (ch1 a)(ch1 b))))))))
    (cond ((null ic)(list nil nil))
	  ((null (cdr ic)) ic)
	  (t (ch1 ic)))))

;; see yn.lisp for discussion of this hackery to randomly choose
;; left/right subtree insertion, since we are inserting from the top
;; down. yn for yes or no.

(let* ((ybc 0)				;yes-no-bit-array-count
       (yl  #.(expt 2 16))		;yes-no-bit-array-length
       (yb (make-array #.(expt 2 16) :element-type '(mod 2)));yes-no-bit-array
       (truenil (make-array 2 :initial-contents '(nil t))))
  (declare (fixnum ybc yl)
	   (type (simple-array (mod 2) (*)) yb)
	   (type (simple-array t (2)) truenil))
  ;; set up the lexically closed array of bits
  (loop for i from 0 to (1-(length yb)) do(setf (aref yb i)(random 2)))
  ;; this next function provides what is needed:  a nil or a t each time it is called
  ;; randomly, and much faster than  (defun nextyn()(= 1 (random 2)))
  (defun nextyn()
    (if (= ybc yl)(setf ybc -1))	;loop around when you exhaust the bits
    (aref truenil (aref yb (incf ybc)))))

(defun heap-insertl (heap x)
    (declare (optimize (speed 3)(safety 0)))
    ;; a heap could be (nil) or  (here left .  right)
    (cond
    
     ((leafp heap) (setf heap
		     (if (hlessfun x heap)
			 (if (nextyn)
			     `(,x ,heap)`(,x nil .,heap))
		       (if (nextyn)
			   `(,heap ,x)`(,heap nil .,x)))))
      ((null (car heap))	(setf (car heap) x))

     ;;(here (car heap))
     ;;(left (car (cadr heap)))
     ;;(right (car (cddr heap)))
     ((hlessfun (car heap) x)  ;;;  if left or right is nil, put it there...
      (cond ((null (cadr heap))(setf (cadr heap) x))
	    ((null (cddr heap))(setf (cddr heap) x))
	    ;;otherwise move down to left or right...
	    ((nextyn)
		 (setf (cddr heap) (heap-insertl (cddr heap) x)))
	    (t   (setf (cadr heap) (heap-insertl (cadr heap) x)))) )
     (t;; x is <= here
      (cond ((null (cadr heap))(setf (cadr heap) (car heap)))
	    ((null (cddr heap))(setf (cddr heap) (car heap)))
	    ((nextyn)
		 (setf (cddr heap) (heap-insertl (cddr heap) (car heap))))
	    (t   (setf (cadr heap) (heap-insertl (cadr heap) (car heap)))))
      (setf (car heap) x)))
    heap)


(defmacro heap-empty-p(x)`(null(car ,x)))

(defun heap-removel(heap)
  ;;(format t "~% remove from ~%~s " heap)
  (if (null (cdr heap)) (prog1 (car heap) (setf (car heap) nil))
    (let ((here (car heap))
	  (newheap (merge-heapsl (cadr heap)(cddr heap))))
      (if (or (null newheap)(leafp newheap))(setf newheap (list newheap nil)))
      (setf (car heap)(car newheap))
      (setf (cdr heap)(cdr newheap))
      here)))


(defun merge-heapsl(a b) ;destructive
  (cond ((null a) b)
	((null b) a)
	((leafp a) (heap-insertl b a))
	((leafp b) (heap-insertl a b))
	((hlessfun (car a)(car b))
	 (setf (cadr a)(merge-heapsl (cadr a)(cddr a)))
	 (setf (cddr a) b)
	 a)
	(t (setf (cadr b)(merge-heapsl (cadr b)(cddr b)))
	   (setf (cddr b) a)
	   b)))
  
;;;;;;;;;;
;;; testing

;; lists

(defun testl-in(heap data)(dolist (i data)(heap-insertl heap i)))
(defun testl-out(heap count)(dotimes (i count)(heap-removel heap)))

;; arrays

;;(defun testa-in(heap data)(dolist (i data)(heap-insert heap i)))
;;(defun testa-out(heap count)(dotimes (i count)(heap-remove heap)))
;; timing results
#|
;; data is 1000 integers between 0 and 1000
(defparameter data (loop for i from 1 to 1000 collect (list (random 1000))))
(defparameter data (loop for i from 1 to 100000 collect (list (random 100000))))


|#


(defun list2heapl(a)
  (let ((h (create-heapl (list a))))
    (loop for i in (cdr a) do (heap-insertl h i))
    h))

(defun heapl2list(h)  (loop while (car h)  collect (heap-removel h)))

	   

(defun maxdepth (h c)			; initially heap 0
  (cond ((null h) c)
	((leafp h) c)
	(t (max (maxdepth (cadr h) (1+ c))
		(maxdepth (cddr h) (1+ c))))))

(defun avdepth (h c)			; initially heap 0
  (cond ((null h) c)
	((leafp h) c)
	(t (* 0.5d0(+ (avdepth (cadr h) (1+ c))
		      (avdepth (cddr h) (1+ c)))))))

;; this is  ok I think. but not so fast as mulhash
(defun mulheap2 (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)))
		      y	      ))
	  (h (create-heapl 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-removel h))
;;       (format t "~%removing r=~s" r)
      ;; (format t "~% (car r)=~s (caar  ans)=~s" (car r) (caar ans))
;;       (format t "~% newterm= ~s"   (* (cdr(caadr 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
;;	     (format t "~% multiplying ~s by ~s incrementing ~s"  (cdaar(cdr r))(cdaddr r) (cdadr ans))
	    (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-insertl h (cons (+ (caar nexti)(caar nextj))
					 alist))))))) 
    (cdr (car ans))))





#|improvements, msort4  --> msort5

(time (mulheap2 x y))  in msort4.  x is 1000 terms, y is 500 terms


; cpu time (non-gc) 982 msec user, 0 msec system
; cpu time (gc)     111 msec user, 15 msec system
; cpu time (total)  1,093 msec user, 15 msec system
; real time  1,125 msec
; space allocation:
;  2,005,460 cons cells, 64 other bytes, 0 static bytes
((4 . 12) (5 . 26) (6 . 10) (8 . 9) (9 . 23) (10 . 4) (11 . 9)
	  (12 . 18) (13 . 11) (14 . 12) ...)
......
 (time (mulheap2 x y)) in msort5
; cpu time (non-gc) 968 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  984 msec user, 0 msec system
; real time  985 msec
; space allocation:
;  1,754,760 cons cells, 0 other bytes, 0 static bytes   <<<<<< note less memory
((4 . 12) (5 . 26) (6 . 10) (8 . 9) (9 . 23) (10 . 4) (11 . 9)
 (12 . 18) (13 . 11) (14 . 12) ...)
|#