;;; -*- Lisp -*-
;;;;; excerpt from newmul with Roman's version of chaining, sort of.
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))
;; set up a lexically available context for handling several
;; variables which would otherwise be global, and require more
;; effort to access from compiled code.
(defmacro heap-empty-p ()
"Returns non-NIL if & only if the heap contains no items."
`(<= fill-pointer 1))
(defmacro hlessfun (a b) `(< (st-exp ,a) (st-exp ,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))))
;; put together in a lexical context to make access faster to global vars
(let ((*YYEXPS* nil)
(*YYCOEFS* nil)
(*YYLEN* 0)
(fill-pointer 0))
(declare (fixnum *YYLEN* fill-pointer)
(type (simple-array t (*)) *YYEXPS* *YYCOEFS*))
(defstruct (st)
;; A kind of stream. each structure contains 5 items
;; which by reference to arrays of XX and YY , can burp out
;; successive exponents and coefficients of the product.
(exp 0 :type integer)
(coef 0 :type integer)
(i 0 :type fixnum)
(v 0 :type integer) (e 0 :type integer)); v, e are read-only
(defun make-simple-array (A) (if (typep A 'simple-array) A
(make-array (length A)
:element-type (array-element-type A)
:initial-contents A)))
(defun mularZ4(x y)
(declare (optimize (safety 3)(speed 0)))
(let* ((xexps (car x))
(xcoefs (cdr x))
(yexps (make-simple-array (car y)))
(ycoefs(make-simple-array (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) (make-st :exp 0 :coef 0 :i ylen :v 0 :e 0))
(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 percolate-down4 (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 (st-exp 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(st-exp aaval))))
hole)
(declare (fixnum child) (type st aaval))
(setf (aref aa hole) aaval hole child))))
;; important: initial-contents of heap must have an element 0 that
;; is not really used. for convenience we could have it equal to the first element.
;; e.g. (3 5 8) should be (3 3 5 8)
(defun heap-insert4 (a xx);; xx is new element
"Insert a new element into the heap. Heap is changed";; rjf
(declare (type (simple-array t (*)) a)(fixnum fill-pointer)
(integer newe)
(optimize (speed 3)(safety 0)))
(let ((c nil) ; potential chain point
(fp (incf fill-pointer)) )
(declare (type (simple-array t (*)) a)(fixnum hole fp)
(optimize (speed 3)(safety 0)))
;; (format t "~%in insert4, fill-pointer is ~s" fill-pointer)
;; (setf (aref a (decf fp)) xx)
(decf fp)
(setf (aref a fp) xx)
(do ((i fp parent)
(parent (ash fp -1);; initialize to parent of fp
(ash parent -1);; and its parent, next time
))
((>= (st-exp xx) (st-exp(aref a parent)));; is xx>= parent? if so, exit
a);; inserted the item at location i.
(declare (fixnum i parent))
;; (format t "~%insert4: c=~s parent=~s xx=~s"c parent xx)
(setf c (aref a parent))
(cond
;; #+ignore
;; program is right if this clause is ignored. But not as fast.
((= (st-exp xx) (st-exp c));; aha, we can insert this item by chaining here.
;(format t "~% insert4 updated ~s " c)
; (format t "~% heap is ~s" a)
(incf (st-coef c) (st-coef xx))
;(format t "~% to ~s " c)
; (format t "~% heap is ~s" a)
(setf (aref a fill-pointer) nil)
(decf fill-pointer);; we don't need the extra heap spot for this item
(if (setf xx (next-term xx)) ;try inserting the successor to this item then.
(heap-insert4 a xx))
(return a)) ;exit from do loop
(t (setf (aref a i) c
(aref a parent) xx) a)))))
(defun heap-remove4(a)
;; returns nil if the heap is empty, otherwise the removed item is returned
(declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0)))
(if (heap-empty-p) nil
(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 fp) nil) ;; optional. makes tracing neater
(setf (aref a (percolate-down4 a 1 last-object fp)) last-object)
xx)))
;; this is the multiplication program
(defun mul-heap4(x y)
(declare (optimize(speed 3)(safety 0)))
(let* ((hh (if (< (length (car x))(length (car y)))(mularZ4 x y)(mularZ4 y x)))
(alen 10)
;; initial heap
(anse (make-array alen :adjustable t :element-type 'integer :fill-pointer 0))
(ansc (make-array alen :adjustable t :element-type 'integer :fill-pointer 0))
(preve 0) (prevc 0)
(top nil)
(newe 0)(newc 0) )
(declare (integer preve prevc newe newc)
(type (array t(*)) anse ansc hh))
(setf fill-pointer (length hh))
(loop while (setf top (heap-remove4 hh)) do
(setf newe (st-exp top) newc (st-coef top))
;; remove top item from heap
(cond
((= preve newe) ; is this same exponent as previous?
(incf prevc newc)) ; if so, add coefficient, and loop
(t (unless (= prevc 0) ; normalize: remove zero coeffs
(vector-push-extend prevc ansc)
(vector-push-extend preve anse))
(setf prevc newc) ; initialize new coefficient
(setf preve newe)))
(if (setf top (next-term top)) ;; is there a successor?
(heap-insert4 hh top))) ;; if so, insert it in heap
;; end of loop, push out the final monomial before exit from routine.
(vector-push-extend prevc ansc)
(vector-push-extend preve anse)
(cons anse ansc)))
(defun next-term(rec) ; record is a stream 5-tuple return succesor or nil
(declare(type (simple-array t(*)) *YYCOEFS* *YYEXPS*)
(optimize (speed 3)(safety 0))
(type st rec))
(let ((i (st-i rec)))
(cond ((>= i *YYLEN*) nil);; no more items on stream of products
(t (incf (the integer (st-i rec)))
(setf (st-coef rec)
(* (st-v rec)(aref *YYCOEFS* i)))
(setf (st-exp rec)
(+ (st-e rec)(aref *YYEXPS* i)))
rec)))))