;; just enough program support to do muldense fast.

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


(defvar mulbuf (make-array 8  :element-type 'integer))

(defun reset-mulbuf(&optional (size 8))
  (setf mulbuf (make-array size :element-type 'integer)))


;; july 12, 2010, help from Duane.. 
;; works even if exp and coefs not simple arrays...

(defun mul-dense (x y) ;; pre-allocated array for workspace
  
  ;;  simple for small degree dense 
  ;; can be much faster than mul-naive, esp. if answer is dense
  ;; exponents and coefficients are integers. 
  (declare    (optimize (speed 3)(safety 0)))
  (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
	   (xc (cdr x))	;array of coefficients of x
	   ;;array of exps and coefficients of y,
	   ;; put in TYPED array so we can iterate faster
	   (yc (make-simple-array (cdr y)))
	   (ye (make-simple-array (car y)))
	   (xei 0)(xci 0)
	   (yl (length ye))
	   (numtermsX (length xe))
	   (numtermsY (length ye))
	   (maxXe (aref xe (1- numtermsX)))
	   (maxYe (aref ye (1- numtermsY)))
	   (ASize (+ 1 maxXe maxYe))
	   ;; (dense-check (dense-check (* numtermsX numtermsY) ASize))
	   (ans mulbuf))
      (declare (fixnum yl) (integer xei xci)
	       (type (simple-array integer (*))  ye yc)
	       (type (array integer (*))  xe xc) ;; we don't know if these are simple
	       (type (simple-array integer (*)) ans) )
      (if (< (length mulbuf) ASize) 
	  (setf ans (reset-mulbuf ASize)))
      (loop for i fixnum from 0 below ASize 
	    do (setf (svref ans i) 0))
      (dotimes (i (length (car x)) 
		  (DA2PAx ans ASize))	;return dense-array to pair of arrays
	(declare (fixnum i))
	(setf xci (aref xc i))
	;;(unless (= xci 0) by hypothesis, input has no zeros
	(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 DA2PAx (A size)
    ;;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.
  (declare (type (simple-array integer (*)) A)
	   (fixnum size)
	    (optimize (speed 3)(safety 0)))
  (let* ((n(count-if-not #'zerop A :start 0 :end size))
         (exps  (make-array n :element-type 'fixnum))
         (coefs (make-array n :element-type 'integer))
         (c 0))
    (declare (fixnum n)
	     (type (simple-array integer (*)) coefs)
	     (type (simple-array fixnum (*)) exps))
     (do ((i 0  (1+ i))
          (j 0))
	 ((>= j n) (cons exps coefs))
       (declare (fixnum i j))
       (unless (eq 0 (setf c (aref A i)))
              (setf (aref exps j) i (aref coefs j)c)
              (incf j)))))


#|  some tests.

(setf expons(let ((ex 0))(coerce (loop for i from 1 to 2000 collect 
				       (prog1 ex (incf ex 
					(1+(random 60))
						       )))
		    '(simple-array fixnum (*)))))

(setf coefs (coerce  (loop for i from 1 to 2000 collect (random 4)) '(simple-array fixnum (*))))

(setf p (cons expons coefs))

(setf q (cons (let ((ex 0))(coerce (loop for i from 1 to 8192 collect 
					 (prog1 ex (incf ex 5)))
				   '(simple-array fixnum (*))))
	      coefs))

(setf f0 (cons #(0 1 41 1681 68921) #(1 1 1 1 1)))

(setf g0 (cons #(0 1 5) #(1 1 1)))

(setf h0 (cons #(0 1) #(10 10)))
(setf xv200 (cons (coerce  (loop for i from 0 to 199 collect i) 'vector)
		  (coerce  (loop for i from 0 to 199 collect 1) 'vector)))

(setf xvs200 (cons (coerce  (loop for i from 0 to 199 collect (* 1000 i)) 'vector)
		 (coerce  (loop for i from 0 to 199 collect 1) 'vector)))
|#