;; new, 7.2010 RJF
;; multiply by Kronecker trick...


(defun mul-kron(p q) ;; p, q are polynomials as #(expons) . #(coefs)
  (let ((v (padpower2 p q)))
    (frombig (mulbig (tobig p v)(tobig q v)) v)))

;; maximum abs value in an array

(defun arraymax (A)(let ((neg 0)(pos 0)(v 0))
		     (loop for i fixnum from 0 below (length A) do
			   (setf v (aref A i))
			   (if (< v 0)(if (< v neg)(setf neg v))
			     (if (> v pos)(setf pos v))))
		     (max pos (- neg))))

(defun padpower2(p1 p2)
 (ceiling (+ (integer-length (min (aref (car p1) (1-(length (car p1))))
				  (aref (car p2) (1-(length (car p2))))
			 )) ;; the max exponents
     (integer-length (arraymax (cdr p1)))
     (integer-length (arraymax (cdr p2))))))

(defmethod lisp2gmpz ((x fixnum))
  (let* ((a (alloc-gmpz))
	 (r (gmpz-z a)))
    (mpz_set_si r x)
    a))

(defun tobigG(p logv)
  (let* ((g (create_mpz_zero))
	 (r (gmpz-z g)))
    (do ((i (1- (length p)) (1- i)))
        ((< i 0) g)
      (mpz_mul_2exp r r logv);; res <- res*2^logv
      (mpz_add r r (gmpz-z(lisp2gmpz (aref p i)))))))

;;;;;;;;;;;;; above, works.  need  frombig, sortof adapted to produce pair of arrays
;;; and we are done.


    

(defun polymultg(p q);; use GMP bignum stuff to multiply polys with gmp coefs
  (let((logv (padpower2g p q)))
    (frombigG2PA (mp* (tobigG q logv)(tobigG p logv)) logv)))




 ;;;(in-package :gmp)
;;{\section{Appendix 2:  polynomials mapped to numbers}
;;\begin{verbatim}


;;; source for polysbyGMP.lisp  September 2003
;;; Using arrays of gmps.

;;; there is an issue here of using gmpz structures or just the raw gmp arrays.


(defun tobigG(p logv)
  (let ((res (create_mpz_zero)))
    (do ((i (1- (length p)) (1- i)))
        ((< i 0) res)
      (mpz_mul_2exp res res logv);; res <- res*2^logv
      (mpz_add res res (togmp (aref p i))))))

(defun frombigG(i logv)
  (let* ((ans nil)
         (q (create_mpz_zero))
         (r (create_mpz_zero))
         (one (create_mpz_zero))        ; will be one
         (vby2 (create_mpz_zero))       ; will be v/2
         (v (create_mpz_zero))          ; will be v
         (signum (>= (signum (aref i 1)) 0))); t if i non-neg
    (mpz_set_si one 1)                  ; set to one
    (mpz_mul_2exp vby2 one (1- logv))   ; set to v/2
    (mpz_mul_2exp v one logv)           ; set to v
    (if signum nil (mpz_neg i i))       ; make i positive
    (while (> (aref i 1) 0)             ; will be 0 when i is zero
      (setf r (create_mpz_zero))
      (mpz_fdiv_r_2exp r i logv)        ; set r to remainder
      (mpz_fdiv_q_2exp q i logv)        ; set q to quotient
      (cond ((> (mpz_cmp r vby2) 0)     ; it is a negative coef.
             (if signum (mpz_sub r r v) (mpz_sub r v r))
             (push r ans)
             (mpz_add i q one) )
            (t
             (if signum nil (mpz_neg r r))
             (push r ans)
             (setf i q))))
    (coerce (nreverse ans)array)))

(defun padpower2g(p1 p2);; p1 and p2 are arrays of gmps
  (+ (integer-length (min (length p1) (length p2)))
     (mpz_sizeinbase(maxcoefg p1) 2)
     (mpz_sizeinbase (maxcoefg p2) 2)))


(defun maxcoefg(p1)  
  ;; assume coefs are gmp numbers, use gmp arith
  ;; accumulate pos and negs separately to avoid computing/storing abs.
  ;; compare at end and make positive.
  (let ((pans (create_mpz_zero))(mans (create_mpz_zero)) temp)
    (do ((i (1- (length p1)) (1- i)))
        ((< i 0)(mpz_neg mans mans)     ; make minus ans into pos
         (if (> (mpz_cmp mans pans) 0) mans pans)); return larger
      (setf temp (aref p1 i))
      (print temp)
      ;; attempt to check the sign, NG
      (if (< (aref temp 1) 0)
          (if (< (mpz_cmp mans temp)0) nil (setf mans temp))
        (if (> (mpz_cmp pans temp) 0)nil (setf pans temp))))))

;; attempt to fix this..  programs changed underneath us??
;; buggy??
(defun maxcoefg(p1)  
  ;; assume coefs are gmp numbers, use gmp arith
  ;; accumulate pos and negs separately to avoid computing/storing abs.
  ;; compare at end and make positive.
  (let ((z (create_mpz_zero))(ans (create_mpz_zero)) (temp (create_mpz_zero)))
    (loop for i fixnum from 0 below (length p1) 
	finally (return ans) do
	  (setf temp (aref p1 i))
	  (format t "~% temp=~s, ans=~s" temp ans)
	  (if (< (mpz_cmp temp z) 0)
	      (mpz_neg temp temp))		;set temp to |temp|
	  (if (< (mpz_cmp temp ans) 0)
	      (setf ans temp)))))


(defun ppowerg(p n) (if (= n 1) p (polymultg p (ppowerg p (1- n)))))


#|
{\section{Appendix 3:  Multivariate polynomials mapped to univariate}
\begin{verbatim}
;;; Simple code for multiplying multivariate polynomials by mapping to univariate

(defun max-expons-multivar(p)
  ;; p is a polynomial in the form below. Return a list of the max exponents
  ;; (#((5 4 3) (9 7 0)) .   ;;exponents
  ;;  #(30 21)  )            ;;coefs
  ;;  30*x^5*y^4*z^3+  21*x^9*y^7*z^0
  (let ((h (map 'array #'(lambda(z)0) (aref (car p) 0)))) ;zero out an array
    (map nil #'(lambda(r) (setf h (map 'array #'max h r)))(car p))
    h))
(defun bound-expons-product(p1 p2)
  ;; returns a list of the maximum exponents that will occur in p1*p2, polynomials
  (map 'list #'+ (max-expons-multivar p1)(max-expons-multivar p2)))

(defun sumlist(h n)(cond((null h) nil) ;; (sumlist '(a b c) 0) returns list (a a+b a+b+c).
                        (t (let((z(+ (car h) n)))
                             (cons  z (sumlist (cdr h) z))))))

(defun mul-multivar-direct (p q)
  ;; use hash tables, multiply directly WITHOUT packing exponents
  (let* ((ans  (make-hash-table :test 'equal ))
         (ep (car p))    (cp (cdr p))
         (eq (car q))    (cq (cdr q))
         (cpi 0)         (epi 0))
    
    (dotimes (i (length ep) 
               ;; this answer is not in the right form.
               ;; needs to be sorted and converted to pair of arrays.
                (l2pam
                 (sort
                  (loop for x being the hash-key in ans 
                        and y being the hash-value in ans unless (= y 0)
                        collect (cons x y))
                  #'lexgp :key #'car)
                 (hash-table-count ans)))
      (setf epi (aref ep i) cpi (aref cp i))
      (dotimes (j (length eq))
        (incf (gethash (map 'list #'+ epi(aref eq j)) ans 0) 
              (* cpi(aref cq j)))))))

(defun l2pam  (z &optional (n (length z)))
    ;; list to pair of arrays, multivar
  (let ((anse (make-array n))
        (ansc (make-array n))
        (exponsize (length (caar z)))
        (j (1- n)))
    
    (declare (fixnum j n exponsize))
    (loop (if (< j 0) (return (cons anse ansc)))
      (setf (aref anse j) (make-array exponsize :initial-contents (caar z)))
      (setf (aref ansc j) (cdar z))
      (decf j)
      (pop z))))

(defun lexgp(a b)
  (cond((null b) nil)
       ((> (car a)(car b)) t)
       (t (and (= (car a)(car b))(lexgp (cdr a)(cdr b))))))

;; this version encodes multivariate into univariate, then multiplies

(defun mul-multivar(r s &optional (mulprog #'mul-naive)) 
  ;; multiply r and s; mulprog is the univariate mult prog to use
  (multiple-value-bind (in out)(make-coders r s) ; set up the encoder and decoder
    (let((prod 
          (funcall mulprog 
                   (cons (map 'array in (car r))(cdr r))
                   (cons (map 'array in (car s))(cdr s)))))
      (cons ;; the re-formed exponents
       (map 'array  out (car prod)) (cdr prod)))))

(defun make-coders(p1 p2);; returns 2 programs for encoding and decoding
  (let* ((z (nreverse(bound-expons-product p1 p2)))
         (fields (mapcar #'integer-length z))
         (sumfields (cons 0(sumlist fields 0)))
         (maskfields  (mapcar #'(lambda(r)(1- (ash 1 r))) fields)))
    ;; encoder  vector to integer
    (values
     #'(lambda(v)
         (reduce #'+ (map 'list  #'ash (reverse v)  sumfields)))
     ;; decoder  integer to vector
     #'(lambda(i)(let ((q (make-array 5 :adjustable t :fill-pointer 0)))
                   (do ((u fields (cdr u))
                        (masks maskfields (cdr masks)))
                       ((null u) (nreverse q))
                     (vector-push-extend (boole boole-and (car masks) i) q)
                     (setf i (ash i (- (car u))))) )))))

\end{verbatim}}

\section{Appendix 4  -  Radix Tree multiplication}
{\begin{verbatim}

s;;; A multi-way tree for storing sparse indexed items. Although I just
;;; made this up, it looks like it is basically around in the
;;; literature, sometimes called a radix tree or a crit bit tree.  It
;;; is an alternative to a hash table or an array, and is, in some
;;; sense, a data structure that can be made more array-like or more
;;; tree-like. The index or key is always an integer, perhaps a
;;; bignum.

;;; The size of the nodes is parameterized: each internal node is a
;;; small array, say of length 8 =2^logsize Then <key, data> is
;;; stored in a node by decomposing the key into 3 =logsize bits at
;;; a time per layer in tree.  thus 30-bit keys will have depth 10
;;; max.  The initial sub-keys are the rightmost bits, and the remaining
;;; keys are computed by shifting the keys by logsize at each level.

;;; The key is never in fact stored in the tree, but computed by position.
;;; There is one subtlety worth noting. 

;;; Note that with logsize 3, we can store items with keys
;;; [0,1,...,7] in the topmost node.  If we need to store an item with
;;; key 9, it has the same trailing 3 bits as 1. To accomodate, we
;;; move the 1 down exactly one level.  That is, the topmost node
;;; looks like [0,p,...,7] where p=[1,9,nil,nil...nil].

;;; so an key like 1  be at the top level (or maybe 1 down).

;;; CON: Keys must be integers or mapped to integers.  Traversing the
;;; tree "in order" is tricky. Numerically adjacent keys will not be
;;; adjacent in the tree, except in the trivial case where all keys
;;; fit in one node.  Nodes may be largely empty.  A key with N bits
;;; will take about (N/logsize) probes.

;;; PRO: we need not know how long the longest key is, which would be
;;; the case if we were using a string-type decomposition starting
;;; from the left.  Short keys will be located near the top of the
;;; tree. The keys are not stored, saving space. By increasing logsize
;;; we can make the performance closer to that of an array.  We never
;;; compare keys; we only extract sub-keys and use them as indexes
;;; into arrays.

;;; author RJF June 13, 2008
;;; made faster, Oct 22, 2008, RJF


;;;; THIS SHOULD BE REPLACED ENTIRELY WITH CODE FROM Nov 21, 2008

(eval-when (:compile-toplevel :load-toplevel :execute)
  (declaim (optimize (speed 3)(safety 1)))
 ;;   (defconstant logsize 5)           ; intermediate node will have 2^logsize slots
  (defconstant logsize 3)               ; best for debugging, not bad for running some tests.
)

(defmacro init-rtree()`(make-array #.(expt 2 logsize) :initial-element 0))
(defmacro leafp(x) `(numberp ,x))
(defmacro nodep(x) `(arrayp ,x))

(defun update-rtree(key val tree &aux node)
  ;; we use this for multiplication or addition of polynomials
  ;; but we could do anything with the update operation below
  ;; Declarations for speed.
  (declare (optimize (speed 3)(safety 0))
           (type (simple-array t (#.(expt 2 logsize))) tree))
  
  (cond  ((fixnump key)(update-rtree-fix key val tree))
         ;;above, optimized for fixnum key
         ((< key #.(expt 2 logsize))
          ;; we have found the level, or level-1 for the key.
          (cond ((nodep (setf node (aref tree key)))
                 ;; if a node here, must descend one level further in tree
                 (update-rtree-fix 0 val node)); and update that node's location 0.
                (t(setf (aref tree key) (+ node 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 #.(1- (expt 2 logsize)))))
              ;; mask off the relevant bits for subkey
              (declare (fixnum ind))
              (setf node (aref tree ind))
              (cond ((arrayp node)
                     (update-rtree (ash key #.(- logsize)) val node))
                    ;;descend into tree with rest of key.
                    (t (setf (aref tree ind)
                         (update-rtree 
                          (ash key #.(- logsize)) val 
                          (if (leafp node)
                              (update-rtree-fix 0 node (init-rtree))
                            (init-rtree)))))))))
  tree)

(defun query-rtree(key tree)            ;works. though we don't use it..
  (declare (optimize (speed 3)(safety 0))
           (type (simple-array t (#.(expt 2 logsize))) tree))
  (cond  ((< key #.(expt 2 logsize)) 
          (if (arrayp (aref tree key)) (query-rtree 0 (aref tree key))
            (aref tree key)))
         (t (let* ((ind(logand key #.(1- (expt 2 logsize))))
                   (h (aref tree ind)))
              (if h (query-rtree (ash key #.(- logsize)) h) nil)))))

(defun maptree (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"
  (declare (optimize (speed 3)(safety 0)))
  (labels 
      ((tt1 (path displace mt);; path = path to the key's location
            (cond ((null mt) nil)
                  ((leafp mt) (funcall fn path mt));; function applied to key  and val.
                  (t;; must be an array
                   (do ((i 0 (1+ i)))
                       ((= i #.(expt 2 logsize)))
                     (declare (fixnum i))
                     (tt1 (+ (ash i displace) path) (+ displace logsize)(aref mt i)))))))
    (tt1 0 0  mt)))

#| example
(setf mt (init-rtree))
(dotimes (i (expt 8 3))(update-rtree i i mt))
mt  ;;; a fully populated b-rtree with value n at location n, 0<= n < 8^3
(ordertree-nz mt)
|#

(defun mul-multivar-rtree(r s) 
  (declare(optimize (speed 3)(safety 0)))
  (multiple-value-bind (in out)(make-coders r s)
    (labels ((incode (z)(mapcar #'(lambda(h)(cons (car h)(funcall in (cdr h)))) z)))
        (setf r (ordertree-nz (mul-rtree (incode r) (incode s))))
        (dolist (h r r)(setf (car h)(funcall out (car h)))))))

(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)))))
            
    
  
(defun mul-rtree (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 (init-rtree)) )
    (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))
        (update-rtree (+ eri(aref es j)) (* cri(aref cs j)) ans)))))

(defun mul-rtree-nosort (r s);multiply two polynomials, in arrays, result in tree
  (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 (init-rtree)) )
    (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                 ;(L2PA (ordertree-nz ans))
                )
      (declare (fixnum i))
      (setf cri (aref cr i) eri(aref er i))
      (dotimes (j les)
        (declare (fixnum j))
        (update-rtree (+ eri(aref es j)) (* cri(aref cs j)) ans)
        ))))

;;; The program update-rt1 below is a more general update function.
;;; First we show how we could use it for multiplication or addition.

(defun update-rtree+ (key val tree) ;; same as update-rtree, using rt1
  (update-rt1 key val tree #'(lambda(val oldval)(+ val (or oldval 0)))))

(defun update-rt1(key val tree fn)
  (labels((update-rtree
           (key val tree)
           ;; we use this for multiplication or addition of polynomials
           ;; but we could do anything with the update operation below
           ;; carefully, perhaps overly, declared for speed.
           (declare (optimize (speed 3)(safety 0))
                    (type (simple-array t (#.(expt 2 logsize))) tree))
  
           (cond  ((fixnump key)(update-rt1-fix key val tree fn))
                  ((< key #.(expt 2 logsize)) ; we have found the level, or level-1 for the key.
                   (cond ((nodep (aref tree (the fixnum key))); if a node here, must descend one level further in tree
                          (update-rtree 0 val (aref tree (the fixnum key)))); and update that node's location 0.
                         (t
                          ;; we have a leaf node to update.
                          (setf (aref tree (the fixnum key)) 
                                (funcall fn val (aref tree (the fixnum key)))))))
                  ;;The key is still too big. Compute the subkey and the rest of the key by masking and shifting.
                  (t (let* ((ind(logand key #.(1- (expt 2 logsize)))); mask off the relevant bits for subkey
                            (h (aref tree ind)))
                       (declare (type (simple-array t (#.(expt 2 logsize))) tree)(fixnum ind))
                       (cond ((arrayp h)
                              (update-rtree (ash key #.(- logsize)) val h));descend into tree with rest of key.
                             ((leafp h) ; leaf, not link. We make a subtree and move the leaf down one level.
                              (setf (aref tree ind)
                                    (update-rtree (ash key #.(- logsize)) val (update-rtree 0 h (init-rtree)))))
                             ;; h is nil means the node was empty. Just make a new subtree and insert.
                             (t (setf (aref tree ind)(update-rtree (ash key #.(- logsize)) val (init-rtree))))))))
           tree))
    (update-rtree key val tree)))


(defun ordertree-nz(mt);;   order-tree the tree removing zeros.  sort after the whole tree is traversed
  (declare (optimize (speed 3)(safety 0))
           (type (simple-array t (#.(expt 2 logsize))) tree))
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
     (declare   (type (array t (*)) ans))
  (labels 
      ((tt1 (path displace mt);; path = path to the key's location
            (cond ((null mt) nil)
                  ;; line below changed to return nil if mt is zero
                  ((leafp mt)(unless (zerop mt) 
               (vector-push-extend (cons path mt) ans)))
                  (t;; mt must be an array
                   (do ((i #.(1- (expt 2 logsize)) (1- i)))
                       ((< i 0) ans)
                      (declare (fixnum i) (type (array t (*)) mt))
                     (tt1 (+ (ash i displace) path) (+ displace logsize)(aref mt i)))))))
    (sort (tt1 0 0 mt) #'< :key #'car))))


(defun update-rtree-fix(key val tree &aux node)  ;;  optimization for  key being fixnum; minor hacks
  (declare (optimize (speed 3)(safety 0))
           (type (simple-array t (#.(expt 2 logsize))) tree)
           (fixnum key))
  
  (cond  ((< key #.(expt 2 logsize))    ; we have found the level, or level-1 for the key.
          (cond ((nodep (setf node (aref tree key))); if a node here, must descend one level further in tree
                 (update-rtree-fix 0 val node)); and update that node's location 0.
                (t(setf (aref tree key) (+ node val)))));; put it here.
         
          (t (let ((ind(logand key #.(1- (expt 2 logsize))))); mask off the relevant bits for subkey
              (declare (fixnum ind))
              (setf node (aref tree ind))
              (cond ((arrayp node)
                     (update-rtree-fix (ash key #.(- logsize)) val node));descend into tree with rest of key.
                    (t (setf (aref tree ind)
                             (update-rtree-fix (ash key #.(- logsize)) val 
                                           (if (leafp node)
                                               (update-rtree-fix 0 node (init-rtree))
                                             (init-rtree)))))))))  
  tree)

(defun update-rt1-fix(key val tree fn)
  (labels((update-rtree;;internally, use fixnum keys
           (key val tree)
           (declare (optimize (speed 3)(safety 0))
                    (type (simple-array t (#.(expt 2 logsize))) tree)
                    (fixnum key))
  
           (cond  ((< key #.(expt 2 logsize)); we have found the level, or level-1 for the key.
                   (cond ((nodep (aref tree (the fixnum key))); if a node here, must descend one level further in tree
                          (update-rtree 0 val (aref tree (the fixnum key)))); and update that node's location 0.
                         (t
                          (setf (aref tree (the fixnum key)) 
                                (funcall fn val (aref tree (the fixnum key)))))))
                  
                  
                  (t (let* ((ind(logand key #.(1- (expt 2 logsize)))); mask off the relevant bits for subkey
                            (h (aref tree ind)))
                       (declare (type (simple-array t (#.(expt 2 logsize))) tree)(fixnum ind))
                       (cond ((arrayp h)
                              (update-rtree (ash key #.(- logsize)) val h));descend into tree with rest of key.
                             ((leafp h) ; leaf, not link. We make a subtree and move the leaf down one level.
                              (setf (aref tree ind)
                                    (update-rtree (ash key #.(- logsize)) val (update-rtree 0 h (init-rtree)))))
                             (t (setf (aref tree ind)(update-rtree (ash key #.(- logsize)) val (init-rtree))))))))
           tree))
    (update-rtree key val tree)))




\end{verbatim}}

\section{Appendix 4: Chinese Remainder Algorithm}
{How does this work?  Let $N = n_0n_1 ...n_k$.  Then the Chinese
Remainder Theorem tells us that we can represent any number x in the
range $-(N-1)/2$ to $+(N-1)/2$ by its residues modulo $n_0,~ n_1, n_2,~ ...,~ n_k$.

Conversion to Normal Form (Garner's Alg.) 
Converting to normal rep. takes $k^2$ steps.

Beforehand, compute
inverse of $n_1$ mod $n_0$, inverse of $n_2$ mod $n_0n_1$, and also the products $n_0n_1$, etc.

Aside: how to compute  these inverses:  Extended Euclidean Algorithm.

Input: x as a list of residues $\{u_i\}$: $u_i = x$ mod $n_i$

Output: x as an integer in $[-(N-1)/2,(N-1)/2]$.  (Other possibilities
include x in another range, also x as a rational fraction!)

Consider the mixed radix representation

$$x = v_0 + v_1*n0 + v-2*(n_0*n_1)+ ... +v_k*(n_0*...n_{k-1}) \eqno(*)$$

if we find the $v_i$, we are done, after $k$ more mults.
$v_0 = u_0,$  each being $x$ mod $n_0$

Next, computing remainder mod $n_1$ of each side of (*)

$u_1 =  v_0 +v_1*n_0   $mod $n_1$,   with everything else dropping out

so $v_0 = (u_1-v_0)* n_0^{-1} $    all mod $n_1$

in general,

  $v_k = ( u_k - [v_0+v_1*n_0 +...+v_{k-1}* (n_0...n_{k-2}) ]) * (n_0*...*n_{k-1})^{-1}$ mod $n_k.$

Note that all the $v_k$ are small numbers and the products of $\{n_i\}$  are precomputed. 
 
Cost: If we find all these values in sequence, we have $k^2$
multiplication operations, where $k$ is the number of primes needed.  In
practice we pick the $k$ largest primes we can, subject to the other constraints
(e.g. we would like the arithmetic to be single-precision, and fast). If
we can do 31-bit signed arithmetic fast, and $L$ is the largest number
in the answer,  $k$ is
about $2\log (L/2^{31}$)}
\section{Appendix 5: Heap management}
Code is provide in the shortprog.tex file.  Needs cleaning up. Also testing data.

\begin{thebibliography}{99}
{\bibitem{atlas}  ATLAS, Automatically Tuned Linear Algebra Software,
{\tt http://math-atlas.sourceforge.net/}}
{\bibitem{brown} Brown, W.S. ``A language and system for symbolic algebra on a digital computer'', 
      in Kuo, F. F. and J. F. Kaiser (eds.), {\em System Analysis by Digital Computer}, John Wiley, New York, 1966 349--369 }
\bibitem{fateman} Richard Fateman. ``Comparing the speed of programs for sparse polynomial multiplication,''
 SIGSAM Bulletin 37, 1 March 2003.
\bibitem{fredman} Fredman, M. L.  `` How good is the information
theory bound in sorting?''  {\em Theoret. Comput. Sci.}, 1:355-361,
1976.  

\bibitem{scj} Johnson, S.C., ``Sparse Polynomial Arithmetic'',
{\em SIGSAM Bull}., vol 8 issue 3 (August 1974) 63---71.

\bibitem{moenck}  Moenck, R. ``Practical Fast Polynomial Multiplication,'' ISSAC 1976 (SYMSAC 1976). ACM.
\bibitem{monagan} Michael Monagan, Roman Pearce.
{ ``Polynomial Division Using Dynamic Arrays, Heaps, and Packed Exponent Vectors''
in {\em Computer Algebra in Scientific Computing},
Lecture Notes in Computer Science Volume 4770/2007, 
2007
295--315}
\bibitem {yozo} Yozo Hida, 
``Data Structures and Cache Behavior of Sparse Polynomial Multiplication,'' Class project CS282, UC Berkeley, May, 2002. {\verb|http://www.cs.berkeley.edu/~fateman/papers/yozocache.pdf|}
\bibitem{zippel}
Zippel, R.
 ``{Interpolating polynomials from their values},''
 {\em J. Symb. Comput. 9}, 3, 375--403, 1990
% volume = {9},
% number = {3},
% year = {1990},
% issn = {0747-7171},
% pages = {375--403},
% doi = {http://dx.doi.org/10.1016/S0747-7171(08)80018-1},
% publisher = {Academic Press, Inc.},
% address = {Duluth, MN, USA},

\end{thebibliography}

\end{document}

(defun pdumps(r varlist)        ;list to string output, various hacks included
    (cond ((numberp r) r)
          (t
           (fix+-  ;; fix up a+-43*b to a-43*b.
            (subseq
             (apply #'concatenate 'string
                  (map 'list  #'(lambda(k)
                         (format nil "+~{~s~}" 
                           ;; do something later for negative coeffs, also -1.
                                 (let ((h (apply #'append
                                              (map 'list 
                                                #'(lambda(r s)
                                                    (cond ((= s 0) nil)
                                                          ((= s 1) (list r))
                                                          (t (list '* r '^ s))
                                                          ))
                                                varlist (cdr k)))))
                                   (cond ((= (car k) 1) h)
                                         ((= (car k) -1)(print h)  (cons '- (if (eq (car h) '*)(cdr h) h)) )
                                         (t(cons (car k) h))))))
                         
                (sort (copy-list r) #'lex> :key #'cdr))) 1))
           )))

;; if you don't like s= "x+-45",  do (fix+- s)

(defun lex> (a b)(cond ((null b) nil)((null a) t)
                       ((= (car a)(car b))(lex> (cdr a)(cdr b)))
                       ((> (car a)(car b)) t)
                       (t nil)))

;; (setf k '((630 14 74 3) (900 10 8 33) (441 18 77 0) (630 14 11 30)))
;; (pdumps k '(x y z))

;;(setf r '((1 0 1)(1 1 0)))            ;x + y
;;(setf s '((1 0 1)(-1 1 0)))           ;x - y
;;(pdumps (ptimes-multivar r r)'(x y)) ==> "x^2+2xy+y^2" .. could put "*" in there too..

(defun ptimes-univar-faster(r s)  ;; store (coef) rather than coef in hash-table. Faster updating we hope
  (declare(optimize (speed 3)(safety 0)))
  (let ((ans  (make-hash-table)))
    (maphash #'(lambda(ex co)           ; exponent,  coefficient in a list.
                 (maphash #'(lambda (ex2 co2)
                              (let* ((newex (+ ex ex2)))
                                (incf (car (or (gethash newex ans nil)(setf (gethash newex ans) (list 0))) )
                                      (* (car co)(car co2)))))r))s)
    ans))
                                
(defun list2hash-faster (k)  ;; k looks like ((30 . (5 4 3))(21 . (9 7 0)))  30*x^5*y^4*z^3+ ..
  (let ((ans (make-hash-table :test 'equal)))
    (dolist (i k ans)(setf (gethash (cdr i) ans) (list (car i))))))

(defun ptimes-multivar-faster(r s)
  (multiple-value-bind (in out)(make-coders r s)
    (labels ((incode (z)(mapcar #'(lambda(h)(cons (car h)(funcall in (cdr h)))) z)))
        (mapcar #'(lambda(h)(cons (car h)(funcall out (cdr h))))
                (hash2list-faster (ptimes-univar-faster (list2hash-faster (incode r))(list2hash-faster (incode s)))))
        )))
(defun hash2list-faster(h)
  (let ((ans nil))
    (maphash #'(lambda(key val)(unless (= 0 (car val)) (push (cons (car val) key) ans))) h)
    ans))

(defun ppower(x n)(if (= n 1) x (ptimes-multivar-faster x (ppower x (1- n)))))


(defun fix+-(s) ;; take a string like "a b+c+-d+e+-ef" and change it to "a b+c-d+e-ef"
  (do ((h (position #\+ s)(position #\+ s :start h)))
      ((null h) s)
     (if (char= (elt s (incf h)) #\-) ;look for a following  #\- ;string must not end with +
         (setf s (delete #\+ s :start (1- h) :end h)))))
    
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; sortx+y.lisp
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; sort X+Y various ways.
(defun make-rl(n d)
  ;;make random list of length n with spacing k, with 0<k<=d between elements
  (let ((res nil) (v (random d)))
    (dotimes (i n (nreverse res))
      (push v res)
      (setf v (+ 1 v (random d))))))

(defun makeX+Y4(x y n m)                ; n is length x, m is length y
  (declare(fixnum n m))
  (let* ((a (make-array (list n m)))
         (c1 0)
         (c2 0))
    (declare (fixnum c1 c2)(optimize (speed 3)(safety 0)))
    (dolist (i x a)
      (setf c2 0)
      (dolist (j y)
        (setf (aref a c1 c2) (+ i j))
        (incf c2))
      (incf c1))))

(defun sort-array(a)
  (let ((v (make-array (apply '* (array-dimensions a)) :displaced-to a)))
    (sort  v #'<)))

(defun makeX+Ylist(x y)
  (declare(fixnum n m)  (optimize (speed 3)(safety 0)))
    (loop for i in x collect 
          (loop for j in y collect (+ i j))))

(defun makeX+Ylistc (x y)
  (declare(fixnum n m)  (optimize (speed 3)(safety 0)))
    (loop for i in x append  ;;; one long list, with coordinates
          (loop for j in y collect (list (+ i j) (cons i j)))))

(defun mergelist(r) ;; r is a list of lists. merge them together.
  (declare (optimize (speed 3)(safety 0)))
  (if (null (cdr r)) (car r)
    (merge 'list (car r)(mergelist (cdr r)) #'<)))

;; this is faster

(defun mergelist2(r) ;; r is a list of lists. merge them together.
  (declare (optimize (speed 3)(safety 0)))
  (if (null (cdr r)) (car r)
    (let ((k (length r)))
      (merge 'list 
             (mergelist2 (nthcdr  (ash k -1) r))
             (mergelist2 (nthcdr  (- k (ash k -1)) (reverse r)) )
             #'<))))

;; to time, 
;; (setf n 1000 m 500)                                  ;(setf n 6 m 5) 

(defun test(n m)
 (progn (setf x (make-rl n 5)) 
 (setf y (make-rl m 5))
;; (format t "~%array version:")
;; (time (sort-array (makeX+Y4 x y n m)))      ;; 3.34 sec
;; (format t "~%~%list version:")
;; (time (mergelist (makeX+Ylist x y )))       ;; 8.30 sec
;; (format t "~%~%binary list merge version:")
;; (time (mergelist2 (makeX+Ylist x y )))     ;; 0.375 sec
;; (format t "~%~%binary list merge version with coordinates:")
;; (time (compresslist (mergelist2v (makeX+Ylistv x y )))) ;; 0.594+0.485GC = 1.079 sec
 (format t "~%~%nconc list merge version with coordinates:")
 (time (compresslist(car (mergelistvnconc (let ((r (makeX+Ylistv x y ))) (setf r (nconc r r)))))))
 ;; 0.688 sec,  0.359 with mymerge.
 (format t "~%~%hashing with coordinates")
 (time (sorthash (makeX+Yhash x y))) ;; 0.55 sec,
  (format t "~%~%heap management with coordinates, full loading")
 (time  (loadheap x y)) ;; first attempt, 2.7 secs. If we do our own vector-pop etc, 0.83 sec.
))

;;;; variations

(defun makeX+Ylistv(x y)
  (declare(fixnum n m)  (optimize (speed 3)(safety 0)))
    (loop for i in x collect 
          (loop for j in y collect (list (+ i j) (cons i j)) )))

(defun mergelist2v(r) ;; r is a list of lists. merge them together.
  (declare (optimize (speed 3)(safety 0)))
  (if (null (cdr r)) (car r)
    (let ((k (length r)))
      (merge 'list 
             (mergelist2v (nthcdr  (ash k -1) r))
             (mergelist2v (nthcdr  (- k (ash k -1)) (reverse r)) ) ;; don;t like this..
             #'<
             :key #'car))))

(defun compresslist(r)   (cl1 r)  r) 
;; changes
;;((4 (0 . 4)) (6 (0 . 6)) (7 (3 . 4)) (7 (0 . 7)) ..
;; to
;;((4 (0 . 4)) (6 (0 . 6)) (7 (0 . 7) (3 . 4)) ...

(defun cl1(r)                   ;destructive..
              (declare(optimize (speed 3)(safety 0)))
              (cond ((null (cdr r)) r)
                    ((= (caar r)(caadr r)); matching values
                     (setf (cdar r)(cons (cadadr r) (cdar r)))
                     (setf (cdr r)(cddr r))
                     (cl1 r))
                    (t (cl1 (cdr r)))) )

(defun mergelistvnconc(r) ;; r is a looped list of lists. merge them together.
  (declare (optimize (speed 3)(safety 0)))
  (cond ((eq (cdr r) r) r) ;; only one element, return it
        (t 
         (setf (car r)(merge 'list (car r)(cadr r) #'< :key #'car)) ; merge first 2 sublists
         (setf (cdr r)(cddr r))         ;remove 2nd sublist
         (mergelistvnconc (cdr r))              ;keep going.
         )))

;; try something like   (time (cl1(car (mergelistvnconc (let ((r (makeX+Ylistv x y ))) (setf r (nconc r r)))))))

(defun mymerge(r s) ;; an iterative mymerge
   (declare (optimize (speed 1)(safety 3)(debug 3)))
  (let((ans nil) (pt nil))
    ;;set initial link
    (cond((null r)(return-from mymerge  s))
         ((null s)(return-from mymerge  r))
         ((< (caar r)(caar s))(setf ans r) (setf r (cdr r))) ;; key is car
         (t (setf ans s)(setf s (cdr s))))
    (setf pt ans)
    
    (loop
      (cond((null r)(setf (cdr pt) s) (return t))
           ((null s)(setf (cdr pt) r) (return t))
           ((< (caar r)(caar s)) ;;key is car
            (setf (cdr pt) r)
            (setf r (cdr r)))
           (t (setf (cdr pt) s)
              (setf s (cdr s))))
;;      (format t "~% r=~s s=~s pt=~s" r s pt)
      (setf pt (cdr pt)))
    ans))
                    
(defun mergelistvnconc(r) ;; r is a looped list of lists. merge them together.
  (declare (optimize (speed 3)(safety 0)))
  (cond ((eq (cdr r) r) r) ;; only one element, return it
        (t 
         (setf (car r)(mymerge (car r)(cadr r))) ; merge first 2 sublists
         (setf (cdr r)(cddr r))         ;remove 2nd sublist
         (mergelistvnconc (cdr r))              ;keep going.
         )))

(defun makeX+Yhash(x y)
  (declare(optimize(speed 3)(safety 0)))
  (let* ((a (make-hash-table)))
    (dolist (i x a)
      (dolist (j y)
        (push (cons i j) (gethash  (+ i j) a nil) )))))

(defun sorthash(h)
  (let ((ans nil))
    (maphash #'(lambda (i v)(push (cons i v) ans)) h)
    (sort ans #'< :key #'car)))

(defun heap2list (h)   (loop for  i from 1 to (heap-count h) collect (heap-remove-fast h)))
(defun list2heap (content) (create-heap :initial-contents content))


;(setf x '(4 5 6 10 15))
;(setf y '(3 4 9 12 14 19))

;; major change so we can go left from any index
;; works.
(defun loadheap (x y)
  (let*  ((x1 (car x))
          (m (maplist #'(lambda (r)(list (+ (car r) x1) (cons r x)))
                      (cons (car y) y)));; add one element to heap initialization this way..
          (h (create-heap-ordered m))
          (ans (list (list -1)))        ;removed later
          (r nil)
          (alist nil)
          (nexti nil)
          (nextj nil))
    (loop
     (if (heap-empty-p h)(return ans))
     (setf r (heap-remove-fast h))
      ;;   (format t "~%removing r=~s~% pushing ~s" r (list (car r) (cons (caaadr r)(cadadr r))))
      ;; remove the smallest entry and see if it is the same as previous one?

     (if (/= (car r)(caar ans))
         (push (list (car r) (cons (caaadr r)(cadadr r))) ans)
       (setf (cdar ans)                 ;combine it with previous item
             (cons (cons (caaadr r)(cadadr r))(cdar ans))))
      ;; in fact we could just compute a1*b1+a2*b2+...
      ;; here we insert into the heap the location to the right if any
     (setf alist (cadr r))
     (cond (alist  
            (setf nexti (car alist))
            (setf nextj (cddr alist))
            (cond ((and nexti nextj)
                   ;;(format t "~%nexti =~s nextj=~s index=~s" nexti nextj  (+ (car nexti)(car nextj)))
                   (heap-insert h (list (+ (car nexti)(car nextj))
                                        (cons nexti nextj)))))   
            )))
    (cdr (nreverse ans))
    ))


;;; 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 (r s)
  `(let ((cr (car ,r))(cs (car ,s)))
     (declare (integer cr cs)) ;gheap
  (cond ((< cr cs) t)
        ((= cr cs) (< (the integer (caaadr ,r))(the integer (caaadr ,s))))
        )))

(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)
      (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))

(defun heap-count (heap)
   (declare (ignore heap))
  (1- fill-pointer))

(defun 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 (HEAP-A HEAP))."
  (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 ---

;;; -*- Lisp -*-
;;; msortx.lisp multiplication of polyns based on copying coefs into arrays.
;;; this uses double-float coefficients

(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; sort X+Y various ways.

(defun make-rac(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 :element-type 'fixnum))
        (coefs(make-array n :element-type 'double-float))
        ;(coefs(make-array n :element-type 'fixnum))
        
        (v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
        ;(+ 1 (random d))
        (+ 1.0d0 (random d))
        )
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))

(defun test(n m)
  (setf x (make-rac n 5))
  (setf y (make-rac m 5))
;  (format t "~% mulhashb time-ar")
;    (time (mulhashb x y))
  (format t "~% mulhashe time-ar")
    (time (mulhashe x y)))

;; multiply two polynomials via a hash-table. For comparison
(defparameter *zcoefs-lim* 50)
(defparameter *zcoefs* (make-array *zcoefs-lim* :element-type 'double-float))

(defun sorthash(h) ;; utility for mulhasha
  (let ((ans nil))
    (maphash #'(lambda (i v)(push (cons i v) ans)) h)
    (sort ans #'< :key #'car)))

(defun sorthashb(h) ;; utility for mulhashb
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)(setf (cdr v)(car v))
                       (setf (car v) i)
                       (unless (= (cdr v) 0.0d0) ;; normalize! remove 0 coefs
                         (vector-push-extend v ans)))
             h)
    (sort ans #'< :key #'car)))

(defun sorthashc(h ar) ;; utility for mulhashc
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)
                 (setf v (aref ar v)); get the double-float, by indirection
                       (unless (= v  0.0d0) ;; normalize! remove 0 coefs
                         (vector-push-extend (cons i v) ans)))
             h)
    (sort ans #'< :key #'car)))

#|(test 1000 500)

 hash time-ar
; cpu time (non-gc) 109 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  125 msec user, 0 msec system
; real time  141 msec
; space allocation:
;  4,519 cons cells, 12,496,136 other bytes, 0 static bytes
  ...)
  using fixnums for coefficients, the time is somewhat lower, but
  dramatically less space is used.

 cl-user(163):  (time (mulhashb x y))
; cpu time (non-gc) 78 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,569 cons cells, 282,152 other bytes, 0 static bytes

                                        
The real problem with coefficients here is that lisp doesn't have 
a convenient unit operation  a:=a+b*c which destroys a in place.

HERE
we can make that happen by doing  (incf (aref a i) (* (aref b j)(aref c k)))
for suitable a,b,c declared arrays.

Using double floats this way,

 (time (mulhashe x y)) ;; for double-floats
; cpu time (non-gc) 62 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,513 cons cells, 388,952 other bytes, 0 static bytes

;; still about 50% of the time is taken in hash-code stuff, giving
;; us an opportunity to do something else.


BUT USING POTENTIAL BIGNUMS...
(setf xx (make-racf 1000 5)) ;; this is important. (make-rac ..) produces simple-arrays of doubles for coefs
(setf yy (make-racf 500 5))

[2] cl-user(391): (time (mulhashf xx yy))
; cpu time (non-gc) 78 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,513 cons cells, 536,208 other bytes, 0 static bytes
#((3 . 20) (6 . 20) (8 . 10) (9 . 12) (11 . 10) (13 . 25) (14 . 10)



|#

(defun mulhashe (x y);;  this hash-table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table :test 'eq)) ;; eq is good enough for fixnums for ACL
        ;;(ans (make-hash-table :test 'eq)) ;; more general, if exponents are arbitrary bignums
        
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefs*)
        (zlim *zcoefs-lim*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
             (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
             (fixnum xei xy zcounter zlim)
             (double-float xco)
             (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
        (cond (loc ;; there is an entry in hash-table for this exponent. update it
               ;; this happens about n*m times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) ;; if we reach this, we need more working space.
                     
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize :element-type 'double-float)))
                       (declare (fixnum newsize)
                                (type (simple-array double-float(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       
                       (setf zcoefs newarray 
                             *zcoefs* newarray 
                             *zcoefs-lim* newsize 
                             zlim newsize))
                                     
                     )))
              )))))



;; attempt to do same as mulhashe but with avl tree.

;; hey, this is really pretty good!
;; return to coefs  being any number type
(defun mulhashf (x y);;  this hash-table does not store double-floats, but indexes into an array of df.
  (let ((ans (make-hash-table :test 'eql)) ;; eq is good enough for fixnums for ACL
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs (make-array 10))
        (zlim 10)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
             (type (simple-array t (*)) xcoefs ycoefs zcoefs)
             (fixnum xei xy zcounter zlim)
            ; (double-float xco)
             (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
        (cond (loc ;; there is an entry in hash-table for this exponent. update it
               ;; this happens about n*m times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) ;; if we reach this, we need more working space.
                           
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize)))
                       (declare (fixnum newsize)
                                (type (simple-array t(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       (setf zcoefs newarray    zlim newsize))
                                     
                     )))
              )))))

(defun make-racf(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 :element-type 'fixnum))

        (coefs(make-array n))
        ;(coefs(make-array n :element-type 'fixnum))
        
        (v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
        ;(+ 1 (random d))
        (+ 1 (random d))
        )
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))



(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))
        ;(coefs(make-array n :element-type 'fixnum))
        
        (v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
        ;(+ 1 (random d))
        (+ 1 (random d))
        )
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))


(defparameter *zcoefs-limg* 50)
(defparameter *zcoefsg* (make-array *zcoefs-limg*))

(defun mulhashg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefsg*)
        (zlim *zcoefs-limg*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
             (type (simple-array t (*)) xcoefs ycoefs zcoefs) ; any coefs
             (fixnum zcounter zlim)
             (integer xei  xy)
            ; (double-float xco)
             (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
        (declare (fixnum j))
        ;; look in the hash-table at location for this exponent
        (setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
        (cond (loc 
               ;; There is an entry in hash-table for this
               ;; exponent. Update corresponding array location. Note:
               ;; hash-table is NOT updated.  this happens about n*m
               ;; times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
               ;; Hash-Table is updated.
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) 
                     ;; if we reach this, we need more working space.
                     ;; This is rarely used code. If global array is big enough, it is never used.
                     ;; If global array is too small, it is made bigger.
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize)))
                       (declare (fixnum newsize)
                                (type (simple-array t(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       (setf zcoefs newarray  zlim newsize
                             *zcoefs-limg* zlim 
                             *zcoefsg* zcoefs))))))))))

;; will this be slower??? Yes.
;; put numbers IN the hash-table. Problem is the frequent re-hashing at line ***

;;; try skiplists
;;;no, this is really slow
(defun mulhashi (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (sl-init)) ;skip list
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (finger nil))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
             (type (simple-array t (*)) xcoefs ycoefs) ; any coefs
             (integer xei )
             (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)ans ) ;the skiplist
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i))        ;coefficient corresponding to that exponent
      (setf finger ans)
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf finger(SLupdate-finger  ans finger (+ xei (aref yexps j))  (* xco (aref ycoefs j))))

        ))))

(defun addar(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (cond 
   ((= (length(car x)) 0) y)
   ((= (length(car y)) 0) x)
   (t
    (let* 
        ((xexps (car x))
         (xi 0)
         (xlim  (1- (length xexps)))
         (xcoefs (cdr x))
         (yexps (car y))
         (yi 0)
         (ylim  (1- (length yexps)))
         (ycoefs (cdr y))
         (zexps (make-array (+ xlim ylim 1) :adjustable t :fill-pointer 0))
         (zcoefs (make-array (+ xlim ylim 1) :adjustable t :fill-pointer 0)))
      (declare (fixnum xi yi ylim xlim)(optimize (speed 3)(safety 0)))
      (loop 
       (cond 
        ;; check first if we have run out of x or y,
        ((> xi xlim)
         (if (> yi ylim) (return (cons zexps zcoefs))
           (let()
             (loop for i from yi to ylim do 
                   (vector-push (aref yexps i) zexps)
                   (vector-push (aref ycoefs i) zcoefs))
             (setf yi (1+ ylim)))))
        ((> yi ylim)
         (let()
           (loop for i from xi to xlim do 
                 (vector-push (aref xexps i) zexps)
                 (vector-push (aref xcoefs i) zcoefs))
           (setf xi (1+ xlim))))
                
        ((< (aref xexps xi)(aref yexps yi))
         (vector-push (aref xexps xi)  zexps)
         (vector-push (aref xcoefs xi) zcoefs)
         (incf xi))
        ((= (aref xexps xi)(aref yexps yi))
         (vector-push (aref xexps xi) zexps)
         (vector-push (+ (aref ycoefs yi) (aref xcoefs xi)) zcoefs)
         (incf xi) (incf yi))
        (t;;(> (aref xexps xi)(aref yexps yi))
         (vector-push (aref yexps yi)  zexps)
         (vector-push (aref ycoefs yi) zcoefs)
         (incf yi))))))))

(defun splitpol(z) ;; take a polynomial and return its two halves
  (let* ((h (length (car z)))
         (half (ash h -1))
         (rest (- h half)))
    (declare (fixnum h half rest)(optimize (speed 3)(safety 0)))
    (values
         (cons (make-array half :displaced-to (car z))
                     (make-array half :displaced-to (cdr z)))
         (cons (make-array rest :displaced-to (car z) :displaced-index-offset half)
                     (make-array rest :displaced-to (cdr z) :displaced-index-offset half)))))

(defun mular(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (let* ((xexps (car x))
         (xe 0)
         (xc 0)
         (xlen  (length xexps)))
    (declare (fixnum  xlim)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
           '(#() . #())) ; empty polynomial
           ((= 1 xlen)                  ; a monomial times a polynomial!
           (setf xe (aref xexps 0))
           (setf xc (aref (cdr x) 0))
           (cons (map 'vector #'(lambda (r)(+ r xe)) (car y))
                 (map 'vector #'(lambda (r)(* r xc)) (cdr y))))
           (t (multiple-value-bind (low hi)
                  (splitpol x)
                (addar (mular low y)(mular hi y)))))))

;; this works. not particularly compact. or fast. mulhashg is much faster.
(defun addar(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  (cond 
   ((= (length(car x)) 0) y)
   ((= (length(car y)) 0) x)
   (t
    (let* 
        ((xexps (car x))
         (xi 0)
         (xlim  (1- (length xexps)))
         (xcoefs (cdr x))
         (yexps (car y))
         (yi 0)
         (ylim  (1- (length yexps)))
         (ycoefs (cdr y))
         (m (+ xlim ylim)) ; maximum possible size.
         (zexps (make-array m :fill-pointer 0))
         (zcoefs (make-array m :fill-pointer 0)))
      (declare (fixnum xi yi ylim xlim)(optimize (speed 3)(safety 0)))
      (loop 
       (cond 
        ;; check first if we have run out of x or y,
        ((> xi xlim)
         (if (> yi ylim) (return (cons zexps zcoefs))
           (let()
             (loop for i from yi to ylim do 
                   (vector-push (aref yexps i) zexps)
                   (vector-push (aref ycoefs i) zcoefs))
             (setf yi (1+ ylim)))))
        ((> yi ylim)
         (let()
           (loop for i from xi to xlim do 
                 (vector-push (aref xexps i) zexps)
                 (vector-push (aref xcoefs i) zcoefs))
           (setf xi (1+ xlim))))
                
        ((< (aref xexps xi)(aref yexps yi))
         (vector-push (aref xexps xi)  zexps)
         (vector-push (aref xcoefs xi) zcoefs)
         (incf xi))
        ((= (aref xexps xi)(aref yexps yi))
         (vector-push (aref xexps xi) zexps)
         (vector-push (+ (aref ycoefs yi) (aref xcoefs xi)) zcoefs)
         (incf xi) (incf yi))
        (t;;(> (aref xexps xi)(aref yexps yi))
         (vector-push (aref yexps yi)  zexps)
         (vector-push (aref ycoefs yi) zcoefs)
         (incf yi))))))))


(defun mulhashg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefsg*)
        (zlim *zcoefs-limg*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
             (type (simple-array t (*)) xcoefs ycoefs zcoefs) ; any coefs
             (fixnum zcounter zlim)
             (integer xei  xy)
            ; (double-float xco)
             (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashg ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
        (declare (fixnum j))
        ;; look in the hash-table at location for this exponent
        (setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
        (cond (loc 
               ;; There is an entry in hash-table for this
               ;; exponent. Update corresponding array location. Note:
               ;; hash-table is NOT updated.  this happens about n*m
               ;; times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
               ;; Hash-Table is updated.
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) 
                     ;; if we reach this, we need more working space.
                     ;; This is rarely used code. If global array is big enough, it is never used.
                     ;; If global array is too small, it is made bigger.
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize)))
                       (declare (fixnum newsize)
                                (type (simple-array t(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       (setf zcoefs newarray  zlim newsize
                             *zcoefs-limg* zlim 
                             *zcoefsg* zcoefs))))))))))
(defun sorthashg(h ar) ;; utility for mulhashgi
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)
                 (setf v (aref ar v)); get the double-float, by indirection
                       (unless (= v  0.0d0) ;; normalize! remove 0 coefs
                         (vector-push-extend (cons i v) ans)))
             h)
    (breakup (sort ans #'< :key #'car))
    ))

(defun breakup (ar) 
  ;; ar is an array of (exp . coef).
  ;; separate them into a pair: one array of exponents and one of coefficients.
  (let* ((ll (length ar))
        (exps (make-array ll))
         (coefs (make-array ll))
         (tt nil))
    (dotimes (i ll (cons exps coefs))
      (setf tt (aref ar i))
      (setf (aref exps i) (car tt))
      (setf (aref coefs i)(cdr tt)))))


(defun PA2L (PA)  
  ;;  Pair of Arrays to List, length
  ;;  returns L a LIST of (exp . coef) of length n.
  ;; also returns n.

  (let* ((exps (car PA))
         (coefs (cdr PA))
         (n (length exps)))
    (declare (fixnum n))
    (do ((i (1- n) (1- i))
         (ans nil (cons (cons (aref exps i)(aref coefs i)) ans)))
        ((< i 0) (values ans n))
      (declare (fixnum i)))))



(defun L2PA (L n) 
       ;;List to Pair of Arrays
       ;; L is a LIST of (exp . coef) of length n.
       ;; create a pair: one array of exponents and one of coefficients.
  (let* ((exps (make-array n))
         (coefs (make-array n))
         (count 0))
    (declare (fixnum count))
    (dolist (i L (cons exps coefs))
      (setf (aref exps count) (car i))
      (setf (aref coefs count)(cdr i))
      (incf count))))
              
(defun addarL(r s);; r and s are sparse polys as Lists of (exp . coef).
  ;; return a similar list of their sum. destructive.
  
   (declare (optimize (speed 1)(safety 3)(debug 3)))
  (let((ans nil) (pt nil))
    ;;set initial link
    (cond((null r)(return-from addarL s))
         ((null s)(return-from addarL r))
         ((< (caar r)(caar s))(setf ans r) (setf r (cdr r))) ;; which arg do we destroy?
         (t (setf ans s)(setf s (cdr s))))
    (setf pt ans)
    
    (loop
      (cond((null r)(setf (cdr pt) s) (return ans))
           ((null s)(setf (cdr pt) r) (return ans))
           ((< (caar r)(caar s)) ;;key is car
            (setf (cdr pt) r)
            (setf r (cdr r)))
           ((= (caar r)(caar s)) 
;            (format t "~% r=~s s=~s pt=~s" r s pt)
            
            (setf (cdadr pt)(+ (cdar r)(cdar s)))
            (setf (cdr pt) r)
            (setf s (cdr s))
            (setf r (cdr r)))
           
           (t (setf (cdr pt) s)
              (setf s (cdr s))))
;;      (format t "~% r=~s s=~s pt=~s" r s pt)
      (setf pt (cdr pt)))))


(defun mularL(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  ;; this works but uses a lot of cons cells and is about 5X slower than mulhashg.
  (let* ((xexps (car x))
         (xcoefs (cdr x))
         (yexps (car y))
         (ycoefs (cdr y))
         (xe 0)
         (xc 0)
         (xlen  (length xexps))
         (ans nil)
         (i (1- (length yexps))))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
           '(#() . #())) ; empty polynomial
           ((= 1 xlen)                  ; a monomial times a polynomial!
            (setf xe (aref xexps 0))
            (setf xc (aref xcoefs 0))
            ;; probably just as fast to do this next line rather than loop
            ;;(map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) (car y)(cdr y))
            (loop
              (if (< i 0)(return t))
              (push (cons (+ xe (aref yexps i)) (* xc (aref ycoefs i))) ans)
              (decf i))
            ans)

           (t (multiple-value-bind (low hi)
                  (splitpol x)
                (addarL (mularL low y)(mularL hi y)))))))


(defun mularL(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  ;; this works but uses a lot of cons cells and is about 5X slower than mulhashg.
  (let* ((xexps (car x))
         (xcoefs (cdr x))
         (yexps (car y))
         (ycoefs (cdr y))
         (xe 0)
         (xc 0)
         (xlen  (length xexps))
         (ans nil)
         (i (1- (length yexps))))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
           '(#() . #())) ; empty polynomial
           ((= 1 xlen)                  ; a monomial times a polynomial!
            (setf xe (aref xexps 0))
            (setf xc (aref xcoefs 0))
            ;; probably just as fast to do this next line rather than loop
            ;;(map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) (car y)(cdr y))
            (loop
              (if (< i 0)(return t))
              (push (cons (+ xe (aref yexps i)) (* xc (aref ycoefs i))) ans)
              (decf i))
            ans)

           (t (multiple-value-bind (low hi)
                  (splitpol x)
                (addarL (mularL low y)(mularL hi y)))))))

;; subdivide both polys but only to some limited depth.
;; 
(defvar *cachelim* 10  )
(defun dosmallmult(x y) (ptimes x y)) ;; make sure result is right form.
  
(defun mularLd(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; but do the intermediate adds using lists.
  
  (let* ((xexps (car x))
         (xcoefs (cdr x))
         (yexps (car y))
         (ycoefs (cdr y))
         (xe 0)
         (xc 0)
         (xlen  (length xexps))
         (ans nil))
    (declare (fixnum  i xlen)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
           nil) ; empty polynomial
           ((= 1 xlen)                  ; a monomial times a polynomial!
            (setf xe (aref xexps 0))
            (setf xc (aref xcoefs 0))
            (map 'list #'(lambda (r s)(cons (+ r xe) (* s xc))) yexps ycoefs))
           ((< xlen *cachelim*)(if (< (length yexps) *cachelim*) 
                                   ;; both x and y are small enough
                                   (mularg x y)
                                 ;; x is small enough, y is not. split it.
                                 (multiple-value-bind (low hi)
                                     (splitpol y)
                                   (addarL (mularL low x)(mularL  hi x)))))
           ;; x is not small enough.
           (t (multiple-value-bind (low hi)
                  (splitpol x)
                (addarL (mularL low y )(mularL hi y)))))))

(defun testsplit(n z)
  (dotimes (i n)(splitpol z)))

(defun mularq(x y);; x and y are sparse polys in arrays as produced by (make-racg n m)
  ;; divide only down to cachelim
  
  (let* ((xexps (car x))
         (yexps (car y))
         (xe 0)
         (xc 0)
         (xlen  (length xexps)))
    (declare (fixnum  xlim)(optimize (speed 3)(safety 0)))
    (cond ((= 0 xlen)
           '(#() . #())) ; empty polynomial
           ((= 1 xlen)                  ; a monomial times a polynomial!
           (setf xe (aref xexps 0))
           (setf xc (aref (cdr x) 0))
           (cons (map 'vector #'(lambda (r)(+ r xe)) (car y))
                 (map 'vector #'(lambda (r)(* r xc)) (cdr y))))
           ((< xlen *cachelim*)(if (< (length yexps) *cachelim*) 
                                   ;; both x and y are small enough
                                   (mulhashgg x y);; or something good for cache
                                 ;; x is small enough, y is not. split it.

                                      (multiple-value-bind (low hi)
                                          (splitpol y)
                                        ;;(format t "~%low=~s~%~%high=~s" low hi)
                                        (addar (mularq low x)(mularq hi x)))))
           (t (multiple-value-bind (low hi)
                  (splitpol x)
                (addar (mularq low y)(mularq hi y)))))))



;; same as mulhashg but simple-array --> array.  Slighly slower.

(defun mulhashgg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefsg*)
        (zlim *zcoefs-limg*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (array t (*)) xexps yexps) ; any exponents
             (type (array t (*)) xcoefs ycoefs zcoefs) ; any coefs
             (fixnum zcounter zlim)
             (integer xei  xy)
            ; (double-float xco)
             (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashg ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
        (declare (fixnum j))
        ;; look in the hash-table at location for this exponent
        (setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
        (cond (loc 
               ;; There is an entry in hash-table for this
               ;; exponent. Update corresponding array location. Note:
               ;; hash-table is NOT updated.  this happens about n*m
               ;; times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
               ;; Hash-Table is updated.
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) 
                     ;; if we reach this, we need more working space.
                     ;; This is rarely used code. If global array is big enough, it is never used.
                     ;; If global array is too small, it is made bigger.
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize)))
                       (declare (fixnum newsize)
                                (type (array t(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       (setf zcoefs newarray  zlim newsize
                             *zcoefs-limg* zlim 
                             *zcoefsg* zcoefs))))))))))


;;end of listing

                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                

;;;; yet another heap-merge-multiplier, using a hacked together top-down heap by rjf
;;;
;;; multiplication of polyns based on sortx+y.lisp
;;;
;;; -*- Lisp -*-
;;; msort5.lisp multiplication of polyns based on
;;; sortx+y.lisp msort2-5 uses heaps implemented as lisp lists.  seems
;;; like a loser, 4X slower than hash, somewhat slower than arrays.
;;; msort5 uses lisp lists to construct a heap, each internal node n
;;; is, in this design (here . (left . right)) where here is a leaf
;;; node, and left/right may either be nil or an internal node or a
;;; leaf node.  A leaf node looks like a dotted pair consisting of
;;; (number . number).  They represent an exponent . coefficient pair.
;;; An empty heap looks like this: (nil)


(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 test(n m)
  (setf x (make-rlc n 5))
  (setf y (make-rlc m 5))
  (format t "~% heap time")
  (time (mulheap2 x y))  
  (format t "~% hash time")
    (time (mulhash x y)))

(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))


;; multiply two polynomials via a hash-table. For comparison
(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)))))))

(defun sorthash(h) ;; utility for mulhash
  (let ((ans nil))
    (maphash #'(lambda (i v)(push (cons i v) ans)) h)
    (sort ans #'< :key #'car)))

;;; HEAP stuff from heap-lisp2.
;;; see heap.lisp for docs.

(defmacro hlessfun (u v) `(< (car ,u) (car ,v)))

;;; implementation of heap data structure in lists.
;;; heapl is heap using lists, not arrays as traditional
;;; this version has leaf nodes as (number . stuff)

(defmacro leafp(x) `(numberp (car ,x)))

(defun create-heapl (ic)                ;initial contents.  ic  should be ordered.
  (labels 
      ((split(z);; used to balance the initial contents.
             (cond ((null z)(values nil nil))
                   (t  (let ((even z)
                             (odd (cdr z)))
                         (do ((i even (cdr i))
                              (j odd (cdr j)))
                             ((null j)(values even odd))
                           (setf (cdr i) (cddr i))
                           (setf (cdr j) (cddr j)))))))
           
       (ch1 (ic)
         (cond ((null ic) nil)
               ((null(cdr ic)) (car ic)) ;final leaf node
               (t (multiple-value-bind (a b)(split (cdr ic))
                    (cons (car ic)(cons (ch1 a)(ch1 b))))))))
    (cond ((null ic)(list nil))
          ((null (cdr ic)) ic)
          (t (ch1 ic)))))

;; see yn.lisp for discussion of this hackery to randomly choose
;; left/right subtree insertion, since we are inserting from the top
;; down. yn for yes or no.

(let* ((ybc 0)                          ;yes-no-bit-array-count
       (yl  #.(expt 2 16))              ;yes-no-bit-array-length
       (yb (make-array #.(expt 2 16) :element-type '(mod 2)));yes-no-bit-array
       (truenil (make-array 2 :initial-contents '(nil t))))
  (declare (fixnum ybc yl)
           (type (simple-array (mod 2) (*)) yb)
           (type (simple-array t (2)) truenil))
  ;; set up the lexically closed array of bits
  (loop for i from 0 to (1-(length yb)) do(setf (aref yb i)(random 2)))
  ;; this next function provides what is needed:  a nil or a t each time it is called
  ;; randomly, and much faster than  (defun nextyn()(= 1 (random 2)))
  (defun nextyn()
    (if (= ybc yl)(setf ybc -1))        ;loop around when you exhaust the bits
    (aref truenil (aref yb (incf ybc)))))

(defun heap-insertl (heap x)
    (declare (optimize (speed 3)(safety 0)))
    ;; a heap could be (nil) or  (here left .  right)
    (cond
     ((null (car heap)) (setf (car heap) x))
     ((leafp heap) (setf heap
                     (if (hlessfun x heap)
                         (if (nextyn)
                             `(,x ,heap)`(,x nil .,heap))
                       (if (nextyn)
                           `(,heap ,x)`(,heap nil .,x)))))

     ;;(here (car heap))
     ;;(left (car (cadr heap)))
     ;;(right (car (cddr heap)))
     ((hlessfun (car heap) x)  ;;;  if left or right is nil, put it there...
      (cond ((null (cadr heap))(setf (cadr heap) x))
            ((null (cddr heap))(setf (cddr heap) x))
            ;;otherwise move down to left or right...
            ((nextyn)
                 (setf (cddr heap) (heap-insertl (cddr heap) x)))
            (t   (setf (cadr heap) (heap-insertl (cadr heap) x)))) )
     (t;; x is <= here
      (cond ((null (cadr heap))(setf (cadr heap) (car heap)))
            ((null (cddr heap))(setf (cddr heap) (car heap)))
            ((nextyn)
                 (setf (cddr heap) (heap-insertl (cddr heap) (car heap))))
            (t   (setf (cadr heap) (heap-insertl (cadr heap) (car heap)))))
      (setf (car heap) x)))
    heap)


(defmacro heap-empty-p(x)`(null(car ,x)))

(defun heap-removel(heap)
  ;;(format t "~% remove from ~%~s " heap)
  (if (null (cdr heap)) (prog1 (car heap) (setf (car heap) nil))
    (let ((here (car heap))
          (newheap (merge-heapsl (cadr heap)(cddr heap))))
      (if (leafp newheap)(setf newheap (list newheap nil)))
      (setf (car heap)(car newheap))
      (setf (cdr heap)(cdr newheap))
      here)))

(defun merge-heapsl(a b) ;destructive
  (cond ((null a) b)
        ((null b) a)
        ((leafp a) (heap-insertl b a))
        ((leafp b) (heap-insertl a b))
        ((hlessfun (car a)(car b))
         (setf (cadr a)(merge-heapsl (cadr a)(cddr a)))
         (setf (cddr a) b)
         a)
        (t (setf (cadr b)(merge-heapsl (cadr b)(cddr b)))
           (setf (cddr b) a)
           b)))
  
;;;;;;;;;;
;;testing stuff    

(defun maxdepth (h c)                   ; initially heap 0
  (cond ((null h) c)
        ((leafp h) c)
        (t (max (maxdepth (cadr h) (1+ c))
                (maxdepth (cddr h) (1+ c))))))

(defun avdepth (h c)                    ; initially heap 0
  (cond ((null h) c)
        ((leafp h) c)
        (t (* 0.5d0(+ (avdepth (cadr h) (1+ c))
                      (avdepth (cddr h) (1+ c)))))))

;; this is OK. I think. but not so fast as mulhash
(defun mulheap2 (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)))
                      y       ))
          (h (create-heapl 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-removel h))
;;       (format t "~%removing r=~s" r)
      ;; (format t "~% (car r)=~s (caar  ans)=~s" (car r) (caar ans))
;;       (format t "~% newterm= ~s"   (* (cdr(caadr 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
;;           (format t "~% multiplying ~s by ~s incrementing ~s"  (cdaar(cdr r))(cdaddr r) (cdadr ans))
            (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-insertl h (cons (+ (caar nexti)(caar nextj))
                                         alist))))))) 
    (cdr (car ans))))





#|improvements, msort4  --> msort5

(time (mulheap2 x y))  in msort4.  x is 1000 terms, y is 500 terms
; cpu time (non-gc) 1,000 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  1,000 msec user, 0 msec system
; real time  1,047 msec
; space allocation:
;  2,005,505 cons cells, 0 other bytes, 0 static bytes
......
 (time (mulheap2 x y)) in msort5
; cpu time (non-gc) 937 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  953 msec user, 0 msec system
; real time  969 msec
; space allocation:
;  1,753,843 cons cells, 32 other bytes, 0 static bytes

.... still, faster and less space is hashing:

(time (mulhash x y))
; 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,501 cons cells, 240,752 other bytes, 0 static bytes

... and using an array-based heap is better than list heap, not as fast as hash.

(time (mulheap x y))
; cpu time (non-gc) 656 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  672 msec user, 0 msec system
; real time  688 msec
; space allocation:
;  505,515 cons cells, 2,184 other bytes, 0 static bytes

;;; inserted here... If the coefficients are doublefloats stored in a separate array, we
;;; have times like this:

;[2] cl-user(367):  (time (mulhashe x y))
; cpu time (non-gc) 62 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  62 msec user, 0 msec system
; real time  63 msec
; space allocation:
;  4,513 cons cells, 388,856 other bytes, 0 static bytes
#((3 . 2.0d0) (4 . 5.0d0) (5 . 3.0d0) (6 . 2.0d0) (7 . 12.0d0)


|#


;;;;msortxd2.  FASTEST HASH MULTIPLIER, set up for arrays of double-floats..

;;; -*- Lisp -*-
;;; msortx.lisp multiplication of polyns based on copying coefs into arrays.
;;; this uses double-float coefficients

(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; sort X+Y various ways.

(defun make-rac(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 :element-type 'fixnum))
        (coefs(make-array n :element-type 'double-float))
        ;(coefs(make-array n :element-type 'fixnum))
        
        (v (random d)))
    (dotimes (i n (cons exps coefs))
      (setf (aref coefs i)
        ;(+ 1 (random d))
        (+ 1.0d0 (random d))
        )
      (setf (aref exps i)v)
      (setf v (+ 1 v (random d))))))


(defun test(n m)
  (setf x (make-rac n 5))
  (setf y (make-rac m 5))
  (format t "~% mulhashb time-ar")
    (time (mulhashb x y))
  (format t "~% mulhashe time-ar")
    (time (mulhashe x y)))

;; multiply two polynomials via a hash-table. For comparison
(defparameter *zcoefs-lim* 50)
(defparameter *zcoefs* (make-array *zcoefs-lim* :element-type 'double-float))

#+ignore ;; bug because it relies on *zcoefs* being large enough. see mulhashe
(defun mulhashc (x y);;  this hash-table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table))
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefs*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
             (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
             ;;(type (simple-array fixnum(*)) xcoefs ycoefs)
             (fixnum xei xy zcounter)
             (double-float xco)
             (optimize (speed 3)(safety 0))
;;           (:explain :types)
             )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
        (cond((null  loc) 
              ;; if loc is nil, allocate a space in the answer array
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)) ; should check for overflow or something, here.
             (t ;; there is an entry in hash-table for this exponent
              (incf (aref zcoefs loc) ;; no number consing here either??
                    (*  xco (aref ycoefs j)))))))))

(defun sorthash(h) ;; utility for mulhasha
  (let ((ans nil))
    (maphash #'(lambda (i v)(push (cons i v) ans)) h)
    (sort ans #'< :key #'car)))

(defun sorthashb(h) ;; utility for mulhashb
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)(setf (cdr v)(car v))
                       (setf (car v) i)
                       (unless (= (cdr v) 0.0d0) ;; normalize! remove 0 coefs
                         (vector-push-extend v ans)))
             h)
    (sort ans #'< :key #'car)))

(defun sorthashc(h ar) ;; utility for mulhashc
  (let ((ans (make-array 10 :adjustable t :fill-pointer 0)))
    (maphash #'(lambda (i v)
                 (setf v (aref ar v)); get the double-float, by indirection
                       (unless (= v  0.0d0) ;; normalize! remove 0 coefs
                         (vector-push-extend (cons i v) ans)))
             h)
    (sort ans #'< :key #'car)))

#|(test 1000 500)

 hash time-ar
; cpu time (non-gc) 109 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  125 msec user, 0 msec system
; real time  141 msec
; space allocation:
;  4,519 cons cells, 12,496,136 other bytes, 0 static bytes
  ...)
  using fixnums for coefficients, the time is somewhat lower, but
  dramatically less space is used.

 cl-user(163):  (time (mulhashb x y))
; cpu time (non-gc) 78 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,569 cons cells, 282,152 other bytes, 0 static bytes

                                        
The real problem with coefficients here is that lisp doesn't have 
a convenient unit operation  a:=a+b*c which destroys a in place.

HERE
we can make that happen by doing  (incf (aref a i) (* (aref b j)(aref c k)))
for suitable a,b,c declared arrays.

Using double floats this way,

 (time (mulhashe x y)) ;; for double-floats
; cpu time (non-gc) 62 msec user, 0 msec system
; cpu time (gc)     16 msec user, 0 msec system
; cpu time (total)  78 msec user, 0 msec system
; real time  78 msec
; space allocation:
;  4,513 cons cells, 388,952 other bytes, 0 static bytes

;; still about 50% of the time is taken in hash-code stuff, giving
;; us an opportunity to do something else.

|#




;;; cleaned up version that doesn't run out of array space..

;; because mulhashd cannot declare zcoefs to be simple-array, this
;; is much slower than hashc.  See hashd that reallocates "by hand"
#+ignore
(defun mulhashd (x y);;  this hash-table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table))
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs  (make-array 5 :element-type 'double-float :adjustable t :fill-pointer 0))
        (loc nil)
        (xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
             (type (simple-array double-float(*)) xcoefs ycoefs) ; not zcoefs
             (fixnum xei xy zcounter)
             (double-float xco)
             (optimize (speed 3)(safety 0))
             )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))

      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
        (cond((null  loc) 
              ;; if loc is nil, allocate a space in the answer array
              ;; first, point to that spot from the hash-table.
              (setf (gethash xy ans) (fill-pointer zcoefs))
              ;; then store the result in zcoefs array.
              (vector-push-extend (* xco (aref ycoefs j)) zcoefs) ; no number consing??
              )                         ; prevents overflow by extending zcoefs.
             (t ;; there is an entry in hash-table for this exponent
              (incf (aref zcoefs loc) ;; no number consing here either??
                    (*  xco (aref ycoefs j)))))))))



;; try to do this right...
(defun mulhashe (x y);;  this hash-table does not store double-floats, but aindexes into an array of df.
  (let ((ans (make-hash-table :test 'eq)) ;; eq is good enough for fixnums for ACL
        ;;(ans (make-hash-table :test 'eq)) ;; more general, if exponents are arbitrary bignums
        
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefs*)
        (zlim *zcoefs-lim*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
             (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
             (fixnum xei xy zcounter zlim)
             (double-float xco)
             (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf loc (gethash (setf xy (+ xei(aref yexps j))) ans))
        (cond (loc ;; there is an entry in hash-table for this exponent. update it
               ;; this happens about n*m times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) ;; if we reach this, we need more working space.
                     
                     #+ignore  ;; if we do this, the array is not a simple-array and we lose efficiency
                     (let*((newsize (ash zlim 1))
                           (newarray (adjust-array zcoefs newsize :element-type 'double-float)))
                        (setf zcoefs newarray 
                             *zcoefs* newarray 
                             *zcoefs-lim* newsize 
                             zlim newsize))
                        
                       
                           
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize :element-type 'double-float)))
                       (declare (fixnum newsize)
                                (type (simple-array double-float(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       
                       (setf zcoefs newarray 
                             *zcoefs* newarray 
                             *zcoefs-lim* newsize 
                             zlim newsize))
                                     
                     )))
              )))))



;; attempt to do same as mulhashe but with avl tree.

(defun multree (x y)
  ;;  this version uses tree not hash-table. it does not store double-floats, but indexes into an array of df.
  (let ((ans (tree-init 0 0)) ;whatever kind of tree.
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefs*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array fixnum (*)) xexps yexps)
             (type (simple-array double-float(*)) xcoefs ycoefs zcoefs)
             (fixnum xei xy zcounter)
             (double-float xco)
             (optimize (speed 3)(safety 0))   )
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i))
      (setf xco (aref xcoefs i))
      (dotimes (j (length yexps))
        (declare (fixnum j))
        (setf loc (treesearch ans (setf xy (+ xei(aref yexps j))))) ;look for entry at index xy
        (cond((null loc) 
              ;; if loc is nil, allocate a space in the answer array
              (treeinsert ans xy zcounter) ; insert into ans, at index xy, the record zcounter
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter *zcoefs-lim*) ;; we need more working space.
                     (let* ((newsize (ash *zcoefs-lim* 1)) ; double the size
                            (newarray (make-array newsize :element-type 'double-float)))
                       (declare (fixnum newsize))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)(setf (aref newarray m)(aref zcoefs m)))
                       (setf zcoefs newarray *zcoefs* newarray *zcoefs-lim* newsize)))))
             
             (t ;; there is an entry in tree for this exponent
              (incf (aref zcoefs loc) ;; no number consing here either??
                    (*  xco (aref ycoefs j)))))))))






;;;;;;;;;;;;;;;;;;;;;;;;
;;; -*- Lisp -*-
;;; mulstream.lisp multiplication of polyns based on streams
;;; general, works, about 2.3X slower than mulhashg, but uses less storage.
;;;


(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; sort X+Y various ways.
(defvar x nil)
(defvar y nil)

(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))
  (format t "~% mulstr6 time")
  (time (mulstr6 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.

;; a heap is just an array.

(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))))


;; construct length(x) streams, a_0*y,  a_1*y, .... 

(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
                            )))))))))




;; set up a context for heap array fill pointer..

(let ((fill-pointer 0)) 

 #+ignore ;; maybe clearer what is going on
  (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."
    (declare (fixnum hole fp)
             (type (simple-array t (*)) aa)(optimize (speed 3)(safety 0)))
        
    (do ((child (lesser-child aa hole fp) (lesser-child aa hole fp)))
        ((or (>= (the fixnum child) fp) (hlessfun xx (aref aa child)))
         hole)
      (declare (fixnum child))
      (setf (aref aa hole) (aref aa child)
            hole child)))
  
  
;; #+ignore;; works also, maybe a teensy bit faster.
  (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 aaval)
               (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))

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


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

;; even better, not adjustable arrays, making sure 2nd arg is longer than first
;; notice storage is way down.
(time (mulstr3 x y))
; cpu time (non-gc) 235 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  235 msec user, 0 msec system
; real time  235 msec
; space allocation:
;  1,037 cons cells, 114,344 other bytes, 0 static bytes
(#(1 3 6 7 8 11 12 13 14 15 ...) . #(20 4 25 16 5 36 25 4 1 8 ...))



|#


;;

(defun mulstr3(x y)  ;;works
  
  (let* ((xl (length (car x)))
         (yl (length (car y)))
         (h (if (> (the fixnum xl)(the fixnum yl))(mularZ2 y x)(mularZ2 x y)));; initial heap
         (ml (+  (the fixnum xl)(the fixnum yl)))
         (r nil)
         (anse (make-array ml :adjustable t :fill-pointer 0))
         (ansc (make-array ml :adjustable t :fill-pointer 0))
         (preve 0)
         (prevc 0))
    (declare (integer preve prevc ) (fixnum xl yl ml))
  ;  (declare (special fill-pointer))
    (setf fill-pointer  (length h))

    (loop while (not (heap-empty-p)) 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 fill-pointer ))
                  )))                   ;re-insert in the heap
    (vector-push-extend prevc ansc)             ;put in last results before exit
    (vector-push-extend preve anse)
    (cons anse ansc)))


;; well, this part works if update-heap works.

  (defun mulstr6(x y);; hack to make heap management cuter.
    (declare (optimize(speed 3)(safety 0)))
    (let* ((hh (if (< (length (car x))(length (car y)))(mularZ2 x y)(mularZ2 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-heap 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)))
      

  ;; WORKS
  (defun update-heap (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))
           (left (aref h 2))
           (right (aref h 3))
           (lefti (car left))
           (righti (car right))
           (rprog (cddr top))
           (ta 0)(tb 0)
           ;(t1 nil)
           (t2 nil))
                                        ; (format t "~%top= ~s ~%left=~s ~%right=~s"top left right)
      (multiple-value-bind (newe newc)
          (funcall (the compiled-function rprog))
             ;; 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) 
                                        ; (format t "#")
              (setf ta (car top) tb (cadr top))
              (heap-remove h)
              (values ta tb))

             ((and                      ;(setf t1 )
               (<= newe lefti)
               (setf t2 (<= newe righti)))
                    ;; insert this node at the top of the heap
                            ;; and return it.
                                    ;; fill-pointer remains the same.
                                        ;                     (format t "^")
              (setf ta (car top)tb (cadr top))
              (setf (car top) newe (cadr top) newc)
              (values ta tb))
                                    ))
  ))

;; the idea here is to burp out one exponent-coefficient pair per call.
;; it works.

  (defun mulstr5(x y) ;; try to make a stream output 
  
  (let* ((xl (length (car x)))
         (yl (length (car y)))
         (h (if (> (the fixnum xl)(the fixnum yl))(mularZ2 y x)(mularZ2 x y)));; initial heap
         (r nil)
         (anse 0)
         (ansc 0)
         (preve 0)
         (prevc 0))
    (setf fill-pointer  (length h))
   ; (format t "~%0h=~s" h)
    #'(lambda() (prog nil
                 again
                  (loop
                    (cond ((heap-empty-p)(setf anse preve ansc prevc)
                                         (setf preve nil prevc nil)
                                         (return t))) ;exit from loop
                    ;; collect all terms with exponent preve
                    ;;  peek at next item
                    (cond ((not (equal preve (car (aref h 1))))
                           (setf anse preve ansc prevc)
          
                           (setf r (heap-remove h))
                           (setf preve (car r) prevc (cadr r))
                           (bubble-up h r)
                           (return t))
                          (t            ; preve= prevc
                           (setf r (heap-remove h))
                           (incf prevc (cadr r))
                           (bubble-up h r)))) ;keep looping
                  (if (equal ansc 0)(go again) ; don't return terms with 0 coef
                    (return (values anse ansc)))))))         

  (defun bubble-up(hh r)
    (declare (optimize (speed 3)(safety 0)))
        (multiple-value-bind (e1 c1);; advance the stream taken from heap
            (funcall (the compiled-function (cddr r)))
          
          (cond (e1                     ; check: is there more on this stream--
                 (setf (car r) e1 (cadr r) c1) ; update the stream
                 (heap-insert hh r fill-pointer))
                )))
   
)     ;; end of let


(defun muls (x y)(let ((k (mulstr5 x y))
                        (anse nil)
                        (ansc nil))
                    (loop (multiple-value-bind (e c) (funcall k)
                      (cond (e (push e anse)(push c ansc))
                            (t (return (cons (nreverse anse)(nreverse ansc)))))))))
;; about as fast as muls, but returns answer as pair of arrays.
(defun mula (x y)(let ((k (mulstr5 x y))
                       (anse (make-array 10 :adjustable t :fill-pointer 0))
                       (ansc (make-array 10 :adjustable t :fill-pointer 0)))
                    (loop (multiple-value-bind (e c) (funcall k)
                            (cond (e (vector-push-extend e anse)
                                     (vector-push-extend c ansc))
                                  (t (return (cons anse ansc))))))))




;; for comparison


(defparameter *zcoefs-limg* 50)
(defparameter *zcoefsg* (make-array *zcoefs-limg*))

(defun mulhashg (x y);;   exponents and coefs may be arbitrary here.
  (let ((ans (make-hash-table :test 'eql)) ;; eql is good enough for bignums too
        (xexps (car x))
        (xcoefs (cdr x))
        (yexps (car y))
        (ycoefs (cdr y))
        (xei 0)
        (xco 0)
        (zcoefs *zcoefsg*)
        (zlim *zcoefs-limg*)
        (zcounter 0)
        (loc nil)
        (xy 0))
    (declare (type (simple-array t (*)) xexps yexps) ; any exponents
             (type (simple-array t (*)) xcoefs ycoefs zcoefs) ; any coefs
             (fixnum zcounter zlim)
             (integer xei  xy)
            ; (double-float xco)
             (optimize (speed 3)(safety 0)))
    (dotimes (i (length xexps)(sorthashc ans zcoefs))
      (declare (fixnum i))
      (setf xei (aref xexps i)) ;exponent for x
      (setf xco (aref xcoefs i));coefficient corresponding to that exponent
      (dotimes (j (length yexps))
        (declare (fixnum j))
        ;; look in the hash-table at location for this exponent
        (setf loc (gethash (setf xy (+ xei (aref yexps j))) ans))
        (cond (loc 
               ;; There is an entry in hash-table for this
               ;; exponent. Update corresponding array location. Note:
               ;; hash-table is NOT updated.  this happens about n*m
               ;; times, where n=size(x), m=size(y).
               (incf (aref zcoefs loc) ;; no number consing here either??
                     (*  xco (aref ycoefs j))))
              (t
               ;; if loc is nil, allocate a space in the answer array.
               ;; This happens relatively infrequently, about k times, where k=size(x*y).
               ;; Hash-Table is updated.
              (setf (gethash xy ans) zcounter)
              ;; store the result in zcoefs array.
              (setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing??
              (incf zcounter)
              (cond ((= zcounter zlim) 
                     ;; if we reach this, we need more working space.
                     ;; This is rarely used code. If global array is big enough, it is never used.
                     ;; If global array is too small, it is made bigger.
                     (let* ((newsize (ash zlim 1)) ; double the size
                            (newarray (make-array newsize)))
                       (declare (fixnum newsize)
                                (type (simple-array t(*)) newarray))
                    ;;   (format t "~%zcoefs now size ~s" newsize)
                       (dotimes (m newsize)
                         (declare (fixnum m))
                         (setf (aref newarray m)(aref zcoefs m)))
                       (setf zcoefs newarray  zlim newsize
                             *zcoefs-limg* zlim 
                             *zcoefsg* zcoefs))))))))))

;;;;;;;;chaining version

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





;;;;;;;;;;;;;;;;;;;; dense, naive algorithm, fixnum..


(defun make-fixnum-array (A) (if (typep A '(simple-array fixnum))
                                 A
                               (make-array (length A)
                                             :element-type 'fixnum
                                             :initial-contents A)))

;; this next program is an attempt to beat the FFT on its home base,
;; but by using the most straight-forward algorithm taking n*m multiplies
;; using a storage scheme that is compact if the answer is dense.

(defun mul-dense-fix (x y)
  ;; fast and simple for small degree dense
  ;; exponents and coefficients are all fixnums, in the answer too.
  ;; if they are not, this answer will be wrong! 
  ;; we want y to be shorter, less copying..
  (if (< (length (car x))(length (car y))) ; reverse args if x is shorter
      (mul-dense-fix y x)

  (let* (
         (xe (car x))                   ;array of exponents of x
         (ye (car y))                   ;array of exponents of y
         (xc (cdr x))                   ;array of coefficients of x
          ;array of coefficients of y, put in typed array so we can iterate faster
         (yc (make-fixnum-array (cdr y)))
         (xei 0)(xci 0)
         (yl (length ye))
         (ans (make-array
               (+ 1 (aref xe (1- (length xe)))
                  (aref ye (1- yl)))    ; size of ans
               :element-type 'fixnum
               :initial-element 0)))
    (declare (fixnum yl xei xci)
             (type (simple-array t (*)) xe ye xc)
             (type (simple-array fixnum (*)) ans yc))
    (dotimes (i (length (car x)) (DA2PA ans))
      (declare (fixnum i))
      (setf xci (aref xc i))
      (unless (= xci 0)
        (setf xei (aref xe i))
        (dotimes(j yl)
          ;; multiply terms 
          (declare (fixnum j))
          (incf (aref ans (the fixnum (+ xei(the fixnum (aref ye j)))))
                (the fixnum(* xci (the fixnum (aref yc j)))))))))))


(defun mul-dense (x y)
  ;;  simple for small degree dense ;; faster than mul-naive if dense
  ;; exponents and coefficients are integers. compare to mul-dense-fix
  ;; if they are not, mul-dense-fix's  answer will be wrong! 
  ;; the idea here is to copy one of the polys to a (dense) array.
  ;; we want y to be shorter, less copying..
  (if (< (length (car x))(length (car y))) ; reverse args if x is shorter
      (mul-dense y x)
  
  (let* (
         (xe (car x))                   ;array of exponents of x
         (ye (car y))                   ;array of exponents of y
         (xc (cdr x))                   ;array of coefficients of x
          ;array of coefficients of y, put in typed array so we can iterate faster
         (yc (make-simple-array (cdr y)))
         (xei 0)(xci 0)
         (yl (length ye))
         (ans (make-array
               (+ 1 (aref xe (1- (length xe)))
                  (aref ye (1- yl)))    ; size of ans
               :element-type 'integer
               :initial-element 0)))
    ;; we could do a reality check here: if the degree is computed above is huge
    ;; compared to the number of terms, maybe this is a bad idea.
    (declare (fixnum yl)
             (type (simple-array integer (*)) xe ye xc)
             (type (simple-array integer (*)) ans yc))
    (dotimes (i (length (car x)) (DA2PA ans))
      (declare (fixnum i))
      (setf xci (aref xc i))
      (unless (= xci 0)
        (setf xei (aref xe i))
        (dotimes(j yl)
          ;; multiply terms 
          (declare (fixnum j))
          (incf (aref ans (+ xei (aref ye j)))
                (* xci (aref yc j)))))))))

)
(defun DA2PA (A) 
    ;;DENSE single array to Pair of Arrays
    ;; A is an ARRAY of (coef coef coef ...), implicit exponents (0 1 2 ...)
    ;; create a pair: one array of exponents and one of coefficients.
  (let* ((n  (count-if-not #'(lambda(r)(eq r 0))
                           A))
         (exps  (make-array n))
         (coefs (make-array n))
         (c 0))
     (declare (fixnum n)(type (simple-array t (*)) A exps coefs)
              (type (array t (*)) exps coefs))
     (do ((i 0  (1+ i))
          (j 0))
         ((> i n) (cons exps coefs))
       (declare (fixnum i j))
            (unless (= 0 (setf c (aref A i)))
              (setf (aref exps j) i (aref coefs j)c)
              (incf j)))))
|#