;; more hacking, 11/16/08 RJF
;; change short nodes to 2-element arrays
;; update elements by (setf (aref ..) (insert-tree...))
;; instead of operating with destructive operations alone. Makes updating possible.
;; also, allow for different size nodes.
;; see r4tree.lisp for comments
;; try hacking more, declaring a few fixnums, perhaps overoptimistically
;; also, we have a program to pick out the smallest index, useful if
;; we want to program a DIVISION!!
(eval-when (:compile-toplevel :load-toplevel :execute)
(declaim (optimize (speed 3)(safety 0)))
(defparameter *mlogsizes* (list -8 -5 -3))
;; (defparameter *mlogsizes* (list -8 -8 -8)) ;alternative, shallower trees
(defparameter *sizes* (mapcar #'(lambda(r)(expt 2 (- r))) *mlogsizes*))
(defparameter *sizesm1* (mapcar #'(lambda(r)(1- (expt 2 (- r)))) *mlogsizes*))
(nconc *sizes* (cddr *sizes*))
(nconc *sizesm1* (cddr *sizesm1*))
(nconc *mlogsizes* (cddr *mlogsizes*)) )
;; temporary, to make tracing possible
#|
(defparameter *mlogsizes* '(-8 -5 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3))
(defparameter *sizes* '(256 32 8 8 8 8 8 8 8 8 8 8))
(defparameter *sizesm1* '(255 31 7 7 7 7 7 7 7 7 7 7))
|#
(defmacro leafp(x) `(numberp ,x))
(defmacro nodep(x) `(arrayp ,x))
(defmacro makeshortnode(key val)`(vector ,key ,val))
(defmacro makefullnode (s)`(make-array ,s :initial-element 0))
(defmacro shortkey-key (h)`(aref ,h 0))
(defmacro shortkey-val (h)`(aref ,h 1))
(defmacro shortnodep(h)`(eq 2 (length ,h)))
(defmacro longnodep(h)`(< 2 (the fixnum (length ,h))))
(defparameter *debug* nil)
(defun db()(setf *debug* (null *debug*)))
(defun irt(key val tree)
(insert-rtree key val tree *sizes* *sizesm1* *mlogsizes*))
#+ignore ;; modify version to debug..
(defun insert-rtree(key val tree sizes sizesm1 mlogsizes &aux node)
(declare(optimize(speed 3)(safety 0))
(type (simple-array t (*)) tree) (fixnum key))
(cond
((shortnodep tree) ; it is a single-entry node
(cond ((eq key (shortkey-key tree))
(setf (shortkey-val tree) val));; case where key matches exactly. replace value
(t
(let ((newtree (insert-rtree
(shortkey-key tree) (shortkey-val tree)
(makefullnode (car sizes))
(cdr sizes) (cdr sizesm1)(cdr mlogsizes))))
;;; insert A NEW FULL-SIZE NODE HERE
(setf tree (insert-rtree
key val newtree (cdr sizes) (cdr sizesm1)(cdr mlogsizes)))))))
;; the node (i.e. tree) is array >2 in length
;; next case is that the key is small enough to fit in this array
((< (the fixnum key)(the fixnum (length tree)))
;; we have found the level, or level-1 for the key.
(cond ((nodep (setf node (aref tree key))); not a leaf
;; if a node here, must descend one level further in tree; key now is 0
(setf (aref tree key)
(insert-rtree 0 val node(cdr sizes) (cdr sizesm1)(cdr mlogsizes)))); and update that node's location 0.
(t
(setf (aref tree key) val) )));; put it here.
;;The key has too many bits to insert at this level.
;;Compute the subkey and separate the rest of the key by masking and shifting.
(t (let ((ind (logand key (the fixnum(car sizesm1))));rightmost logsize bits
(newkey (ash key (the fixnum(car mlogsizes))))); reduced key
(declare (fixnum ind newkey))
(setf node (aref tree ind))
(cond ((not (nodep node));; a single item, must be moved
(cond ((eq (the fixnum node) 0)
(setf (aref tree ind) (makeshortnode newkey val)))
(t
(insert-rtree newkey val
(setf (aref tree ind) (makefullnode (car sizes)))
;; (setf (aref tree ind) (makefullnode (cadr sizes)))
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;;sizes sizesm1 mlogsizes
)
(insert-rtree 0;;ind
node (aref tree ind)
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;; sizes sizesm1 mlogsizes
))))
(t
(setf (aref tree ind)
(insert-rtree
newkey val node
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;; sizes sizesm1 mlogsizes
)))))))
tree)
;; debug this
(defun insert-rtree(key val tree sizes sizesm1 mlogsizes &aux node)
(declare(optimize(speed 3)(safety 0))
(type (simple-array t (*)) tree) (fixnum key))
(cond
((shortnodep tree) ; it is a single-entry node
(cond ((eq key (shortkey-key tree))
(setf (shortkey-val tree) val));; case where key matches exactly. replace value
(t
(let ((newtree (insert-rtree
(shortkey-key tree) (shortkey-val tree)
(makefullnode (car sizes)) ;;***
sizes sizesm1 mlogsizes
;; (cdr sizes) (cdr sizesm1)(cdr mlogsizes)
)))
;;; insert A NEW FULL-SIZE NODE HERE
(setf tree (insert-rtree
key val newtree
;;(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
sizes sizesm1 mlogsizes
))))))
;; the node (i.e. tree) is array >2 in length
;; next case is that the key is small enough to fit in this array
((< (the fixnum key)(the fixnum (length tree)))
;; (assert (= (length tree) (car sizes)))
;; we have found the level, or level-1 for the key.
(cond ((nodep (setf node (aref tree key))); not a leaf
;; if a node here, must descend one level further in tree; key now is 0
(setf (aref tree key)
(insert-rtree 0 val node(cdr sizes) (cdr sizesm1)(cdr mlogsizes)))); and update that node's location 0.
(t
(setf (aref tree key) val) )));; put it here.
;;The key has too many bits to insert at this level.
;;Compute the subkey and separate the rest of the key by masking and shifting.
(t (let ((ind (logand key (the fixnum(car sizesm1))));rightmost logsize bits
(newkey (ash key (the fixnum(car mlogsizes))))); reduced key
(declare (fixnum ind newkey))
(setf node (aref tree ind))
(cond ((not (nodep node));; a single item, must be moved
(cond ((eq (the fixnum node) 0)
(setf (aref tree ind) (makeshortnode newkey val)))
(t
(insert-rtree newkey val
(setf (aref tree ind) (makefullnode (car sizes)))
;; (setf (aref tree ind) (makefullnode (cadr sizes)))
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;;sizes sizesm1 mlogsizes
)
(insert-rtree 0;;ind
node (aref tree ind)
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;; sizes sizesm1 mlogsizes
))))
(t
(setf (aref tree ind)
(insert-rtree
newkey val node
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;; sizes sizesm1 mlogsizes
)))))))
tree)
(defun urt(key val tree)
(update-rtree key val tree *sizes* *sizesm1* *mlogsizes*))
(defun update-rtree(key val tree sizes sizesm1 mlogsizes &aux node)
(declare(type (simple-array t (*)) tree) (fixnum key))
(cond
((shortnodep tree) ; it is a single-entry node
(cond ((eq key (shortkey-key tree))
(incf (shortkey-val tree) val));; case where key matches exactly. replace value
(t
(let ((newtree (insert-rtree
(shortkey-key tree) (shortkey-val tree)
(makefullnode (car sizes))
sizes sizesm1 mlogsizes
;;(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
)))
;;; insert A NEW FULL-SIZE NODE HERE
(setf tree (update-rtree
key val newtree
;;(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
sizes sizesm1 mlogsizes
))))))
;; the node (i.e. tree) is array >2 in length
;; next case is that the key is small enough to fit in this array
((< (the fixnum key)(the fixnum (length tree)))
;; we have found the level, or level-1 for the key.
(cond ((nodep (setf node (aref tree key))); not a leaf
;; if a node here, must descend one level further in tree; key now is 0
(setf (aref tree key)
;; (insert-rtree 0 val node(cdr sizes) (cdr sizesm1)(cdr mlogsizes))
(update-rtree 0 val node(cdr sizes) (cdr sizesm1)(cdr mlogsizes))
))
; and update that node's location 0.
(t
(incf (aref tree key) val) )));; put it here.
;;The key has too many bits to insert at this level.
;;Compute the subkey and separate the rest of the key by masking and shifting.
(t (let ((ind (logand key (the fixnum (car sizesm1))));rightmost logsize bits
(newkey (ash key (the fixnum (car mlogsizes))))); reduced key
(declare (fixnum ind newkey))
(setf node (aref tree ind))
(cond ((not (nodep node));; a single item, must be moved
(cond ((eq node 0)
(setf (aref tree ind) (makeshortnode newkey val)))
(t
(update-rtree newkey val
;; (setf (aref tree ind) (makefullnode (car sizes)))
(setf (aref tree ind) (makefullnode (cadr sizes)))
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;;sizes sizesm1 mlogsizes
)
(update-rtree 0;;ind
node (aref tree ind)
(cdr sizes)(cdr sizesm1)(cdr mlogsizes)
;;sizes sizesm1 mlogsizes
))))
(t
(setf (aref tree ind)
(update-rtree
newkey val node
(cdr sizes) (cdr sizesm1)(cdr mlogsizes)
;;sizes sizesm1 mlogsizes
)))))))
tree)
(defun qrt(key tree) ;; get value for key in rtree tree. 0 means "no entry"
(query-rtree key tree *sizes* *sizesm1* *mlogsizes*))
(defun query-rtree(key tree sizes sizesm1 mlogsizes) ;; returns 0 if no entry in the tree for that key
(declare (fixnum key))
(cond
((leafp tree) 0)
((shortnodep tree)
(if (eq (shortkey-key tree) (ash key (car mlogsizes)))(shortkey-val tree) 0))
;; case 1: the key should be in this node.
((< key (the fixnum (car sizes)))
(let ((z (aref tree key)))
(if (nodep z);; if the node is an array
(query-rtree 0 z (cdr sizes)(cdr sizesm1)(cdr mlogsizes));; look in the zero location
z)));; the node is NOT an array: this is the value.
(t (let* ((ind(logand key (car sizesm1)));;otherwise compute the offset
(h (aref tree ind))) ; and the subtree to use for further search
(if h (query-rtree (ash key (car mlogsizes)) h (cdr sizes)(cdr sizesm1)(cdr mlogsizes)) 0)))))
(defun map-rtree (fn mt);; order appears jumbled, actually it is in a reversed-radix order
;; apply fn, a function of 2 arguments, to key and val, for each entry in tree.
;; order is "radix reversed"
(labels
((tt1 (path displace mt mlogsizes);; path = path to the key's location
(declare (fixnum path displace))
(cond ((leafp mt)
(unless (= 0 mt)(funcall fn path mt)));; function applied to key and val.
((shortnodep mt)
(funcall fn (+ path (ash (the fixnum (shortkey-key mt))displace ))
(shortkey-val mt)))
(t;; must be an array
(let ((m (length mt)))
(declare (fixnum m))
(do ((i 0 (1+ i)))
((= i m ))
(declare (fixnum i))
(tt1 (+ (the fixnum (ash i displace)) path)
(- displace (car mlogsizes)) ;um, fix this?
(aref mt i) (cdr mlogsizes))))))))
(tt1 0 0 mt *mlogsizes* )))
(defun ordertree-nz(mt)
(declare (type (simple-array t (*)) tree))
(let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
(declare (type (array t (*)) ans))
(map-rtree #'(lambda(key val)
(vector-push-extend (cons key val) ans))
mt)
(sort ans #'(lambda(a b)(declare (fixnum a b)) (< a b)) ;faster because of declarations?
;;#'<
:key #'car)))
(defun mul-rtree4 (r s) ;multiply two polynomials, in arrays
(declare (optimize (speed 3)(safety 0)))
(let* ((er (car r)) (es (car s))
(cr (cdr r)) (cs (cdr s))
(cri 0) (eri 0)
(ler (length er))
(les (length es))
(ans (makefullnode (car *sizes*))))
(declare (type (simple-array integer (*)) cr cs er es)
(fixnum ler les)) ;maybe some others are fixnums too..
;; do the NXM multiplies into the answer tree
(dotimes (i ler (A2PA (ordertree-nz ans)))
(declare (fixnum i))
(setf cri (aref cr i) eri(aref er i))
(dotimes (j les)
(declare (fixnum j))
(urt (+ eri(aref es j)) (* cri(aref cs j)) ans
)))))
;; we could save space/time by doing a:=a+b*c in one step if we are doing GMP arith.
(defun mul-rtree4-nosort(r s) ;multiply two polynomials, in arrays
(declare (optimize (speed 3)(safety 0)))
(let* ((er (car r)) (es (car s))
(cr (cdr r)) (cs (cdr s))
(cri 0) (eri 0)
(ler (length er))
(les (length es))
(ans (makefullnode (car *sizes*))))
(declare (type (simple-array integer (*)) cr cs er es)
(fixnum ler les)) ;maybe some others are fixnums too..
;; do the NXM multiplies into the answer tree
(dotimes (i ler ans ) ;;; (A2PA (ordertree-nz ans))) ;;dont sort
(declare (fixnum i))
(setf cri (aref cr i) eri(aref er i))
(dotimes (j les)
(declare (fixnum j))
(urt (+ eri(aref es j)) (* cri(aref cs j)) ans)))))
(defun A2PA(A) ;; array of (exp . coef) to 2 arrays.
(let* ((LA (length A))
(V nil)
(anse (make-array LA))
(ansc (make-array LA)))
(dotimes (i LA (cons anse ansc))
(setf V (aref A i))
(setf (aref anse i) (car V))
(setf (aref ansc i) (cdr V)))))
#|
;; example loop iteration..
;; (loop for i across z when (> i 0) do (print i))
;; (loop for i across z maximizing i)
;; this picks the item with the smallest index and returns a pair:
;; (index . val)
;; useful for division!
(defun map-rtreexx (fn mt);; order appears jumbled, actually reversed-radix order
;; apply fn, a function of 2 arguments, to key and val, for each entry in tree.
;; order is "radix reversed"
(labels
((tt1 (path displace mt);; path = path to the key's location
(declare (fixnum path displace) (special *sm*))
(cond ((leafp mt) (unless (= 0 mt)(funcall fn path mt)
));; function applied to key and val.
((shortnodep(aref mt 0))
(funcall fn (+ (ash (the fixnum (aref mt 2)) displace) path) (aref mt 1)))
(t;; must be an array
(do* ((i 0 (1+ i))
(path2 path (+ (ash i displace) path)))
((= i #.(expt 2 logsize)))
(declare (fixnum i))
;; we cut this off if
;; (aref mt i) is a leaf and index is too big
(if (and (leafp (aref mt i))
(> path2 *sm*)) (return) ;from do*
(tt1 path2 (+ displace logsize)(aref mt i))))))))
(tt1 0 0 mt)))
'(defun smallest(mt)
(let ((*sm* most-positive-fixnum)
(val 0))
(declare (special *sm*))
(map-rtreexx #'(lambda(k v)
; (format t "~%testing k=~s"k)
(if (< k *sm*)
(setf *sm* k val v)))
mt)
(cons *sm* val)))
|#
;;; testing, compare simply insertion of random numbers into a structure
#|
(setf mt (makefullnode 256))
(irt 3 3 mt)
(irt 259 259 mt)
(irt 16777216 16777216 mt)
|# ;testing, comparing to hash table
(defun t1(N r)
(let ((k 0)
(mt (makefullnode 256)))
(declare (fixnum k))
(dotimes (i N mt)
(setf k (random r))
(irt k k mt))))
(defun t2(N r);; use hash table
(let ((k 0)
(ht (make-hash-table)))
(declare (fixnum k))
(dotimes (i N ht)
(setf k (random r))
(setf (gethash k ht) k))))
(defun t3(N r);; just do random
(let ((k 0)
(ht (make-hash-table)))
(declare (fixnum k))
(dotimes (i N ht)
(setf k (random r))
(fff k) ;; calling a dummy function so compiler does not remove previous line.
)))
(defun fff(k) k)
(defun t4(N r);; use array. sometimes faster than t3!
(let ((k 0)
(ar (make-array r :initial-element 0)))
(declare (fixnum k) (type (simple-array t (*)) ar))
(dotimes (i N ar)
(setf k (random r))
(setf (aref ar k) k))))
;; (t1 1000000 2000000)is 1031ms
;; (t2 1000000 2000000)is 781ms; more memory
;; (t3 1000000 2000000)is 140ms [just calling (random), then dummy fn.
;; (t4 1000000 2000000)is 203ms
;; (t1 1000000 20000)is 656ms, redo 344
;; (t2 1000000 20000)is 516ms; more memory , redo 281
;; (t3 1000000 20000)is 140ms , redo 188
;; (t4 1000000 20000)is 156ms , redo 125
;; (t1 1000000 2000)is 250ms
;; (t2 1000000 2000)is 344ms; more memory , redo 250
;; (t3 1000000 2000)is 140ms, redo 125
;; t4 is 156 msg ;
;; (t1 1000000 255)is 188ms, redo 203
;; (t2 1000000 255)is 288ms; much more memory, redo 266
;; (t3 1000000 255)is 140ms, redo 187
;; t4 is 172; redo 109 (!) [cost of array update is less than dummy function call]
;; subtract off the time for the loop and the generation of random numbers..
;; conclusion: hashing is competitive for insertion, maybe 35-40% faster often
;; highlight: rtree is faster by 3X if all elements are in the first level node
;; for items of size up to 20,000 compare 516 vs 376, favoring hashing
;; for items of size up to 2,000 compare 110 vs 104, favoring rtree
;; for items of size up to 255 compare 48 vs 148, favoring rtree
;; for items of size up to
;; 2,000,000 compare 891 vs 641, favoring hashing
;; BUT
;; even if we have to allocate an array of size 2,000,000, t4 using arrays
;; is fastest. with times of 43ms, 16ms, 16ms, 32ms. [we redid these timings
;; with a longer (X10) loop. They are in the right ballpark.]
;; putting this into the mix, using an array is between 14X faster and 3X faster
;; for insertion.
;; even if we have to allocate an array of size 20,000,000, t4 using arrays
;; is fastest. with times of 43ms, 16ms, 16ms, 32ms. [we redid these timings
;; with a longer (X10) loop. They are in the right ballpark.]
;; putting this into the mix, using an array is between 14X faster and 3X faster
;; for insertion.
;; what about an array of size 200,000,000 ? My lisp system refused to allocate it.
;; Error: An allocation request for 800000016 bytes caused a need for [too much memory from OS]
;; however, the other guys still worked.
;; (time (prog nil (t1 1000000 200000000) nil))
;; took 1703 ms. (radix)
;; t2 1015ms (hash)
;; t3 172 ms (empty loop)
;; t4 can't do it.
;; t1 used ; 484,367 cons cells, 56,991,112 other bytes, 0 static bytes
;; t2 used ; ; 22 cons cells, 63,469,160 other bytes, 0 static bytes
;; conclude: simple array is a winner unless it can't fit. otherwise hashing is good.
#| more timing.
cl-user(32): (setf q (make-racg 10000 50 5)r (make-racg 10000 50 5))
(#(45 58 82 130 152 202 212 241 285 305 ...)
. #(5 2 1 4 5 1 5 1 3 4 ...))
cl-user(33): (time (mul-rtree4 q r))
; cpu time (non-gc) 42,718 msec user, 0 msec system
; cpu time (gc) 954 msec user, 15 msec system
; cpu time (total) 43,672 msec user, 15 msec system
; real time 44,344 msec
; space allocation:
; 584,645 cons cells, 16,295,368 other bytes, 0 static bytes
(#(94 107 128 131 141 165 178 179 191 201 ...)
. #(5 2 20 1 8 4 20 4 8 5 ...))
cl-user(38): (time (setf ans3 (mul-heap q r)))
; cpu time (non-gc) 64,171 msec (00:01:04.171) user, 0 msec system
; cpu time (gc) 47 msec user, 0 msec system
; cpu time (total) 64,218 msec (00:01:04.218) user, 0 msec system
; real time 65,938 msec (00:01:05.938)
; space allocation:
; 280 cons cells, 10,928,736 other bytes, 0 static bytes
(#(94 107 128 131 141 165 178 179 191 201 ...)
. #(5 2 20 1 8 4 20 4 8 5 ...))
[1] cl-user(39): (length (car ans3))
510519
[1] cl-user(40): (aref (car ans3) (1- *))
511991
[1] cl-user(41): (time (mul-naivef q r))
; cpu time (non-gc) 49,844 msec user, 0 msec system
; cpu time (gc) 422 msec user, 16 msec system
; cpu time (total) 50,266 msec user, 16 msec system
; real time 51,266 msec
; space allocation:
; 1,022,065 cons cells, 4,090,064 other bytes, 8,192 static bytes
(#(94 107 128 131 141 165 178 179 191 201 ...)
. #(5 2 20 1 8 4 20 4 8 5 ...))
;; try this (setf q (make-racg 5000 50 5)r (make-racg 1000 50 5))
(time (mul-dense-fix q r)) is 0.078 sec
(time (mul-dense q r)) is 0.219 sec
(time (mul-fft q r)) is 1.41 sec
(time (mul-rtree4 q r)) is 1.625 sec
(time (mul-hash1 q r)) is 1.687 sec
(time (mul-heap5 q r)) is 2.53 sec ;; in newmulr2s.fasl
(time (mul-naive q r)) is 4.03 sec
(setf q (make-racg 5000 500 5) r (make-racg 1000 500 5))
(time (mul-dense-fix q r)) 0.203 sec
(time (mul-dense q r)) 0.375 sec
(time (mul-heap5 q r)) 3.37 sec
(time (mul-hash1 q r)) 8.1 sec
(time (mul-rtree4 q r)) 8.4 sec
(time (mul-naive q r)) 129.19 sec
(time (mul-fft q r)) ;stopped. out of storage.
(setf q (make-racg 5000 1 5)r (make-racg 1000 1 5))
(time (mul-fft q r)) 0.016 sec
(time (mul-dense-fix q r)) 0.031 sec
(time (mul-dense q r)) 0.188 sec
(time (mul-toom q r)) 0.218 sec
(time (mul-naive q r)) 0.30 sec
(time (mul-karat q r)) 0.500 sec
(time (mul-rtree4 q r)) 0.67 sec
(time (mul-hash1 q r)) 0.81 sec
(time (mul-heap5 q r)) 2.17 sec ; from newmulr2s
(time (mul-heap5 q r)) 2.76 sec ; from newmulr2
(setf q (make-racg 5000 10000 5)r (make-racg 1000 10000 5)) ;; really sparse
(time (mul-heap5 q r)) 4.703 sec
(time (mul-hash1 q r)) 22.875 sec 329 megabytes
(time (mul-rtree4 q r)) 22.875 sec
(time (mul-naive q r)) 320.812 sec much less space... 74 megabytes
|#
;;; try to put the heap stuff from newmulr2s in this file too
;;; -*- Lisp -*-
;;;;; excerpt from newmul with Roman's version of chaining, sort of.
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))
;;; this tries to do f:=f+a*b instead of c:=a*b; f:=f+c. could be faster with GMP.
;;; attempt to shorten because the "insert at top" heuristic doesn't really work often enough.
(defmacro heap-empty-p ()
"Returns non-NIL if & only if the heap contains no items."
`(<= fill-pointer 1))
(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)
;;; use the line below if exponents are fixnums!!
((< (the fixnum (st-exp (aref ,a left)))(the fixnum (st-exp (aref ,a right)))) left)
;; ((< (st-exp (aref ,a left))(st-exp (aref ,a right))) left) ;more general
(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 6 items, not 5 as in previous version
;; which by reference to arrays of XX and YY , can burp out
;; successive exponents and coefficients of the product.
(exp 0 :type integer)
(coef1 0 :type integer) ;; could be doubles, or gmp numbers or ..
(coef2 1 :type integer) ;; the second item, to be multiplied by coef1 and added..
(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) ;; set up initial heap.
(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 (type (simple-array t (*))
ans xexps xcoefs yexps ycoefs))
(setf *YYEXPS* yexps *YYCOEFS* ycoefs *YYLEN* ylen)
(setf (aref ans 0) (make-st :exp 0 :coef1 0 :coef2 1 :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) :coef1 (* 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)
(< (the fixnum xxval) (the fixnum(st-exp aaval))) ;faster, less general
;; (< xxval (the fixnum(st-exp aaval))) ;faster, less general
;; (< xxval (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 a. 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 fill-pointer)
(optimize (speed 3)(safety 0)))
(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))
(setf c (aref a parent))
(cond
((= (st-exp xx) (st-exp c));; aha, we can insert this item by chaining here.
(incf (st-coef c) (* (st-coef1 xx) (st-coef2 xx)))
;; Line above is one op-code in GMP, using no persistent storage outside c...
;(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, eliminated: peeking at top
(defun mul-heap5(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)(newc1 0) (newc2 1) )
(declare (integer preve prevc newe newc1 newc2)
(type (array t(*)) anse ansc hh))
(setf fill-pointer (length hh))
(loop while (not (heap-empty-p))
do
(setf top (heap-remove4 hh))
(setf newe (st-exp top) newc1 (st-coef1 top) newc2 (st-coef2 top))
(cond ((= preve newe) ; is this same exponent as previous?
(incf prevc (* newc1 newc2))) ; 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 (* newc1 newc2)) ; 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-coef1 rec)(st-v rec))
(setf (st-coef2 rec)(aref *YYCOEFS* i)) ;; don't multiply yet, save for later
;; later we can do c:=c+a*b in one fell swoop.
(setf (st-exp rec)
(+ (st-e rec)(aref *YYEXPS* i)))
rec))))
)
#|
(defun h(z) (setf u (make-racg z 10000 5)v (make-racg z 10000 5))
(setf ansr (mul-rtree4 u v) ansn (mul-naive u v))
(pol= ansr ansn))
;; test for same answer..
|#