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