(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))
(defstruct (st (:type list))
(exp 0 :type integer)(coef 0 :type integer) (i 0 :type fixnum)
(v 0 :type integer) (e 0 :type integer) ;read-only
)
;; we could change the type from list to something else, but then hlesspfun needs to change
(let ((*YYEXPS* nil)( *YYCOEFS* nil)( *YYLEN* 0)(fill-pointer 0))
(declare (fixnum *YYLEN* fill-pointer))
;; similar to mularZ2, but the each stream consists of 4 items:
;; current exponent, current coefficient, multiplier for this
;; stream, and index of the next array element to extract from the
;; second polynomial. Example:
;; Let X = a[0]*x^e[0]+ a[1]*x^e[1]+ ....
;; Let Y = b[0]*x^f[0]+ b[1]*x^f[1]+ ....
;; the first stream will represent a[0]* Y
;; the second stream will represent a[1]* Y ...
;; the first stream will start out as
;; { e[0]+f[0] , a[0]*b[0], a[0], 1 a[0] e[0]} = {exp,coef,v,i, e}
;; the second { e[1]+f[1], a[1]*b[0], 2, a[1], a[1], e[1]}
;; updating the first stream will do two things:
;; 1. return exp and coef
;; 2. update the stream to contain { e[0]+f[1] , a[0]*b[1], 2, a[0], e[0]}
;; if a stream runs off the end of Y, exp is returned as NIL rather than a number.
;; 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
(defun mularZ4(x y)
(let* ((xexps (car x))
(xcoefs (cdr x))
(yexps (car y))
(ycoefs (cdr y))
(yexp0 (aref yexps 0))
(ycoef0 (aref ycoefs 0))
(ans (make-array (1+(length xexps))))
(ylen (length (car y)))
(xe 0)
(v 0)) ;
(declare (fixnum i) (type (simple-array t (*))
ans xexps xcoefs yexps ycoefs))
(setf *YYEXPS* yexps *YYCOEFS* ycoefs *YYLEN* ylen)
(setf (aref ans 0) (list 0 0 0 0 ylen))
(dotimes (j (length xexps) ans)
(setf v (aref xcoefs j))
(setf xe (aref xexps j))
(setf (aref ans (1+ j))
(make-st :exp (+ xe yexp0) :coef (* v ycoef0) :i 1 :v v :e xe) ))))
(defun update-simheap (h)
;;return top of heap h, which is necessarily not empty.
;; take successor in stream of top of heap, and insert it in heap.
(declare(type (simple-array t(*)) h)(optimize (speed 3)(safety 0)))
(let* ((top (aref h 1))
(ta (car top))(tb (cadr top))
(newe (updatesimfun top)))
;; in case newe <= lefti and newe <= righti. put it at the top.
;; avoiding all heap programs
;; top of heap is location 1 (ignore 0)
;; left branch is location 2
;; right branch is location 3
(cond((null newe) (heap-remove h));end of stream from this subproduct
;; we've modified the top node! leave it there!
((and(<= newe (car (aref h 2))) (<= newe (car (aref h 3)))))
;; it's not the top node. remove it and put the (revised) version back in
(t (heap-insert h (heap-remove h) fill-pointer))); the new item doesn't fit there
(values ta tb)))
(defun mulsim(x y);; hacked to insert at top of heap if possible
(declare (optimize(speed 3)(safety 0)))
(let* ((hh (if (< (length (car x))(length (car y)))(mularZ4 x y)(mularZ4 y x)))
;; initial heap
(anse (make-array 10 :adjustable t :fill-pointer 0))
(ansc (make-array 10 :adjustable t :fill-pointer 0))
(preve 0)
(prevc 0))
(declare (integer preve prevc))
(setf fill-pointer (length hh))
(loop while (not (heap-empty-p)) do
(multiple-value-bind
(newe newc)
(update-simheap hh)
(declare (integer newe newc))
; take next term off heap and insert its successor
(cond ((= preve newe) ; is this same exponent as previous?
(incf prevc newc)) ; if so, add coefficient.
(t (unless (= prevc 0); normalize: remove zero coeffs
(vector-push-extend prevc ansc)
(vector-push-extend preve anse))
(setf prevc newc)
(setf preve newe))); initialize new coefficient
))
(vector-push-extend prevc ansc) ;put in last results before exit
(vector-push-extend preve anse)
(cons anse ansc)))
#+ignore;; somehow this doesn't compile right.
(defun updatesimfun(rec) ; record is a stream 5-tuple
(declare(type (simple-array t(*)) *YYCOEFS* *YYEXPRS*) )
(let ((i (caddr rec)));; (st-i rec fails?)
(cond ((>= i *YYLEN*) nil);; no more items on stream of products
(t (incf (st-i rec))
(setf
(st-coef rec) (* (st-v rec) (aref *YYCOEFS* i)))
(setf
(st-exp rec) (+ (st-e rec) (aref *YYEXPS* i)); returns this val
)))))
;; this works
(defun updatesimfun(rec) ; record is a stream 5-tuple
(declare(type (simple-array t(*)) *YYCOEFS* *YYEXPRS*) )
(let ((i (caddr rec))
(v (cadddr rec))
(e (car(cddddr rec))))
(cond ((>= i *YYLEN*) nil);; no more items on stream of products
(t (incf (caddr rec))
(setf (cadr rec) (* v (aref *YYCOEFS* i)))
(setf (car rec) (+ e (aref *YYEXPS* i)); returns this val
)))))
(defmacro heap-empty-p ()
"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 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))."
`(let* ((left;;(+ (the fixnum ,parent)(the fixnum ,parent))
(ash ,parent 1))
(right (1+ (the fixnum left))))
(declare (fixnum left right ,fp ,parent )
(optimize(speed 3)(safety 0))
(type (simple-array t (*)) ,a))
(cond ((>= left ,fp) ,fp)
((= right ,fp) left)
((hlessfun (aref ,a left) (aref ,a right)) left)
(t right))))
;; We want to get fast access to a sort-of global variable "fill-pointer"
;; without incurring the overhead for a full global access. We can have
;; the best of both worlds by setting up a lexical environment in which
;; fill-pointer is bound, and defining programs in that environment.
;; The names of the programs (like percolate-down) are available in
;; any context including the global context.
;; (heap management parts of this were inspired by code that is
;; Copyright (c) 2002, 2003 Gene Michael Stover., GPL.
;; The code was found ;; at http://cybertiggyr.com/lisp-heap/
(defun percolate-down (aa hole xx fp)
"Private. Move the HOLE down until it's in a location suitable for X.
Return the new index of the hole."
(let ((xxval (car xx)))
(declare (fixnum hole fp)(integer xxval)
(type (simple-array t (*)) aa)(optimize (speed 3)(safety 0)))
(do* ((child (lesser-child aa hole fp) (lesser-child aa hole fp))
(aaval (aref aa child)(aref aa child)))
((or (>= (the fixnum child) fp) (< xxval (the fixnum(car aaval))))
hole)
(declare (fixnum child))
(setf (aref aa hole) aaval hole child))))
(defun percolate-up (a hole xx)
"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) xx)
(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 xx (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 xx fp)
"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 fp) nil)
(incf fill-pointer) ; the global one
;; 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 fp xx)) xx))
;;; Later, 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)
(declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
(let* ((xx (aref a 1))
(fp (decf fill-pointer)) ;faster, fp is local now
(last-object (aref a fp)))
(declare (fixnum fp))
(setf (aref a (the fixnum (percolate-down a 1 last-object fp))) last-object)
xx))
;;; end of the heap stuff. Note that this has been hacked to provide
;;; less generality and some consequent advantages in time or
;;; space. The original package on which this is based is far more flexible.
;;;
)
;; end of lex environment for *YY..*