;;; -*- 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))))
)