;;; experimentation in encoding polynomials for faster multiplication.
;;; karatsuba but using a sparse encoding!!
;;; July, 2010
;;; Richard Fateman

;; try to get it all to work with exponent array fixnums.


(eval-when (compile load)
  (declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
  (declaim (inline svref = eql make-array)))

;;An implementation of Karatsuba multiplication

#|Assume p is of degree 2n-1,  i.e.  p= a*x^n+ b  where a is degree n-1 and b is degree n-1.
 assume q is also of degree 2n-1,  i.e.  q= c*x^n+d 

compute U= (a+b)*(c+d)  =  a*c+a*d+b*c+b*d,  a polynomial of degree k=2*(n-1) 
compute S= a*c   of degree k
compute T= b*d   of degree k
compute U = U-S-T = a*d+b*c of degree k

Then V=p*q = S*x^(2n)+U*x^n + T.

from first principles,  V =p*q should be of degree 2*(2n-1) = 4n-2. 
from the form above, k+2n= 4n-2.
note that U*x^n is of degree 2*(n-1)+n =  3n-1 so that it OVERLAPS the S*x^(2n) component.
note that T     is of degree 2*(n-1)   =  2n-2 so that it OVERLAPS the U*x^(n)   component.


  (the trick is that, used recursively to multiply a*c etc,
    this saves coefficient mults, instead of nm, max(n,m)^1.5)

if p, q, are of degree (say) 10 or so, use conventional arithmetic. Just a guess.
What about this 2n-1 business? Find the max degree of p, q.  If it is odd, we
are set.  If it is even, 2n, do this:
  p=a*x^n+ b  where b is of degree n-1, but a is of degree n. Same deal works. I think.
|#

;;(defvar *karatlim* 12)
;;(defparameter *karatlim* 40)  ;; an educated guess -- don't use karatsuba for small probs
(defvar *karatlim* 4) ;; for testing

;; representation is (cons #(e0 e1 ....)  #(c0 c1 ...))  ; array of expons and coefs, ascending degree.

(defun degr(p)	  (aref (car p)(1-  (length (car p)))))


(defun truncatepoly(p v)		;  break up polynomial p into two pieces at degree v. 
  ;; (0 ...) (v ....).. actually, the 2nd value will be the first exponent >= v, not necessarily v.
  (declare (fixnum v)(optimize (speed 3)(safety 0)))
  (let ((spot(position v (car p) :test #'<=))
	(carp (car p))(cdrp (cdr p)))
    (declare (fixnum spot))
    (if (null spot)(values p '( #(). #()))
    (values (cons (make-array spot :displaced-to carp :element-type 'fixnum)
		  (make-array spot :displaced-to cdrp))
	    (cons (make-array (- (length carp) spot ) :displaced-to carp 
			      :displaced-index-offset spot :element-type 'fixnum)
		  (make-array (- (length carp) spot) :displaced-to cdrp :displaced-index-offset spot))))))

#+ignore
(defun truncatepoly(p v)		;  break up polynomial p into two pieces at degree v. 
  ;; (0 ...) (v ....).. actually, the 2nd value will be the first exponent >= v, not necessarily v.
  (declare (fixnum v)(optimize (speed 3)(safety 0)))
  (let ((spot(position v (car p) :test #'<=))
	(carp (car p))(cdrp (cdr p)))
    (declare (fixnum spot))
    (if (null spot)(values p '( #(). #()))
    (values (cons (make-array spot :displaced-to carp)
		  (make-array spot :displaced-to cdrp))
	    (cons (make-array (- (length carp) spot ) :displaced-to carp 
			      :displaced-index-offset spot 
			      ;:element-type 'fixnum
			      )
		  (make-array (- (length carp) spot) :displaced-to cdrp 
			      :displaced-index-offset spot))))))



;;#+ignore
;;; THIS IS IT

(defun mul-ks(pp qq) ;; by Karatsuba.,  use displaced arrays ;; sparse rep.
  (let* ((plim (length (car pp)))	; number of terms in p
	 (qlim (length (car qq))))	; number of terms in q
    (declare (fixnum qlim plim))
    (if (or (< plim *karatlim*)
	    (< qlim *karatlim*))
;;	(mul-dense pp qq) 

	#+ignore
      (mul-dense (cons(make-simple-array (car pp))
			(make-simple-array (cdr pp)))
		   (cons(make-simple-array (car qq))
			(make-simple-array (cdr qq))))
      
      (mul-dense-ns2 pp qq)
      
      ;; need to have a version of mul-dense that does not require simple arrays.
      


      ;; main program
      (let ((v (+ 1 (ash (max (degr pp)(degr qq)) -1)))) ;half of the  max degree of inputs, +1
	(multiple-value-bind
	    (A B)	  (truncatepoly pp v)
	  (setf B (kpolyshiftleft B (- v)))
	  (multiple-value-bind
	      (C D)    (truncatepoly qq v)
	    (setf D (kpolyshiftleft D  (- v)))
	    (let ((U  (mul-ks (kpolyadd A B)  (kpolyadd C D)))
		  (TT (mul-ks A C))
		  (S  (mul-ks B D))
		  (aa nil))
					;	    (format t "~%U=~s~%TT=~s~%S=~s" U TT S)
	      (setf U (kpolysub U S))
	      (setf U (kpolysub U TT))
					;	    (format t "~%U=~s" U)
	      (setf aa (kpolyaddshift TT U v))
	      (kpolyaddshift aa S  (+ v v))
	      )))))))

#+ignore ;; can't use this..
(defun kpolyshiftleft(p k)		; like multiplying by x^k, destructive
  (let ((A (car p)))
    (loop for i fixnum from 0 below (length A)
	finally (return p) do
	  (incf (aref A i) k)  )))

(defun kpolyshiftleft(p k)		; like multiplying by x^k, NON-destructive
  (let ((A 
	 (make-array (length (car p)):element-type 'fixnum :initial-contents (car p))
	 ))
    (loop for i fixnum from 0 below (length A)
	finally (return (cons A (cdr p))) do
	  (incf (aref A i) k)  )))

(defun kpolyadd(p q)
  (let* ((exps nil)(coefs nil)
	 (carp (car p))			;exponents of p ascending
	 (cdrp (cdr p))
	 (carq (car q))
	 (cdrq (cdr q))
	 (i (1- (length (car p))))
	 (j (1- (length (car q))))
	 (cc 0)
	 (carpi (aref carp i))
	 (carqj (aref carq j)))

    (loop 
     (cond ((< j 0)			; clean up
	    (loop for k fixnum from i downto 0 do 
		  (push (aref carp k) exps)
		  (push (aref cdrp k) coefs))
	    (return))

	   ((< i 0)			; clean up
	    (loop for k fixnum from j downto 0 do 
		  (push (aref carq k) exps)
		  (push (aref cdrq k) coefs))
	    (return)))
    ;  (format t "~% i=~s j=~s" i j)
     (setf carpi (aref carp i) carqj (aref carq j))
      (cond ((= carpi carqj) 
	     (cond ((= 0 (setf cc (+ (aref cdrp i)(aref cdrq j))))nil)
		    (t (push carpi exps)  (push cc coefs)))
	    (decf i)(decf j))
	   ((> carpi carqj) (push carpi exps)
	    (push (aref cdrp i) coefs)
	    (decf i))
	   (t  (push carqj exps)
	       (push (aref cdrq j) coefs)
	       (decf j))))
    (cons (coerce exps '(array fixnum (*)))(coerce coefs '(array integer (*))))))


(defun kpolyaddshift(p q n)  ;;p+ x^n*q
  (let* ((exps nil)(coefs nil)
	 (carp (car p))			;exponents of p ascending
	 (cdrp (cdr p))
	 (carq (car q))
	 (cdrq (cdr q))
	 (i (1- (length (car p))))
	 (j (1- (length (car q))))
	 (cc 0)
	 (carpi (aref carp i))
	 (carqj (aref carq j)))

    (loop 
     (cond ((< j 0)			; clean up
	    (loop for k fixnum from i downto 0 do 
		  (push (aref carp k) exps)
		  (push (aref cdrp k) coefs))
	    (return))

	   ((< i 0)			; clean up
	    (loop for k fixnum from j downto 0 do 
		  (push (+ n (aref carq k)) exps)
		  (push (aref cdrq k) coefs))
	    (return)))
    ;  (format t "~% i=~s j=~s" i j)
     (setf carpi (aref carp i) carqj (+ n (aref carq j)))
      (cond ((= carpi carqj) 
	     (cond ((= 0 (setf cc (+ (aref cdrp i)(aref cdrq j))))nil)
		    (t (push carpi exps)  (push cc coefs)))
	    (decf i)(decf j))
	   ((> carpi carqj) (push carpi exps)
	    (push (aref cdrp i) coefs)
	    (decf i))
	   (t  (push carqj exps)
	       (push (aref cdrq j) coefs)
	       (decf j))))
    (cons (coerce exps 
		  '(array fixnum (*))		  ) 
	  (coerce coefs '(array integer (*))))))

(defun kpolysub(p q)
  (let* ((exps nil)(coefs nil)
	 (carp (car p))			;exponents of p ascending
	 (cdrp (cdr p))
	 (carq (car q))
	 (cdrq (cdr q))
	 (i (1- (length (car p))))
	 (j (1- (length (car q))))
	 (cc 0)
	 (carpi (aref carp i))
	 (carqj (aref carq j)))

    (loop 
     (cond ((< j 0)			; clean up
	    (loop for k fixnum from i downto 0 do 
		  (push (aref carp k) exps)
		  (push (aref cdrp k) coefs))
	    (return))

	   ((< i 0)			; clean up
	    (loop for k fixnum from j downto 0 do 
		  (push (aref carq k) exps)
		  (push (- (aref cdrq k)) coefs))
	    (return)))
    ;  (format t "~% i=~s j=~s" i j)
     (setf carpi (aref carp i) carqj (aref carq j))
      (cond ((= carpi carqj)
	     (cond ((= 0 (setf cc (- (aref cdrp i)(aref cdrq j))))nil)
		   (t (push carpi exps)  (push cc coefs)))
	    (decf i)(decf j))
	   ((> carpi carqj) (push carpi exps)
	    (push (aref cdrp i) coefs)
	    (decf i))
	   (t  (push carqj exps)
	       (push (aref cdrq j) coefs)
	       (decf j))))
    (cons (coerce exps '(array fixnum (*))		  )
	  (coerce coefs '(array integer (*))))))

  
;; benchmarks
#|
(setf ones (cons (coerce (loop for i from 0 to 999 collect i) '(array fixnum (*)))
		 (coerce (loop for i from 1 to 1000 collect 1) 'array)))

;;; don't declare exponents to be fixnum....

(setf ones (cons (coerce (loop for i from 0 to 999 collect i) 'simple-array )
		 (coerce (loop for i from 1 to 1000 collect 1) 'simple-array)))


(setf lones (cons (coerce (loop for i from 0 to 9 collect i) '(array fixnum (*)))
		  (coerce (loop for i from 1 to 10 collect 1) 'array)))
(setf jones (cons (coerce (loop for i from 0 to 9 collect i) '(array fixnum (*)))
		 (coerce (loop for i from 1 to 10 collect 1) 'array)))

(setf ss (cons (coerce (loop for i from 0 to 2 collect i) '(array fixnum (*)))
		 (coerce (loop for i from 1 to 3 collect 1) 'array)))


|#


;; copied from muldense.lisp
;; and changed..


#+ignore
(defun mul-dense (x y)
;;  (format t "~%x is ~s, .  ~s  y is ~s . ~s "	  (type-of (car x))(type-of (cdr x))(type-of (car y))(type-of (cdr y)))
  
  (if (typep (car x) '(simple-array fixnum)) nil     (setf x (cons(coerce (car x) '(array fixnum (*))) (cdr x))))
;  (if (typep (cdr x) 'simple-array) nil     (setf x (cons(car x) (coerce (cdr x) '(array integer (*))))))

  (if (typep (car y) '(simple-array fixnum)) nil     (setf y (cons(coerce (car y) '(array fixnum (*))) (cdr y))))
;  (if (typep (cdr y) 'simple-array) nil     (setf y (cons(car y) (coerce (cdr y) '(array integer (*))))))

;;    (format t "~%after conversion ~%xp is ~s, .  ~s  y is ~s . ~s "	  (type-of (car x))(type-of (cdr x))(type-of (car y))(type-of (cdr y)))
  ;;  simple for small degree dense 
  ;; can be much faster than mul-naive, esp. if answer is dense
  ;; exponents and coefficients are integers. compare to mul-dense-fix
  ;; if they are not, mul-dense-fix's  answer will be wrong! 
  (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)))
	   (yc (cdr 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 (make-array
		 dense-check
		 :element-type 'integer ;; assumes coeffs are integers only
		 ;; note: an integer can be a bignum.
		 :initial-element 0)))
      (declare (fixnum yl)
	       
	       (type (array integer (*)) xc yc) 
	       (type (array fixnum (*)) xe ye )
	       (type (array integer (*)) xe ye)
	       (type (array integer (*)) ans) 
	       (optimize (speed 3)(safety 0))
	       )
#+ignore      (cond ((= 0 dense-check) (error "mul-dense is not be a good choice for polynomial mult of ~s by ~s terms, since it needs to allocate an array of size ~s" numtermsX numtermsY ASize)
	     (error "leaving mult program, need to fill in some alternative")))

      (dotimes (i (length (car x)) 
		 (DA2PA ans)
		 )

	(declare (fixnum i))
;		(print 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)))))))))




;;; another version with start/stop in arrays, rather than displaced arrays which slow up



(defmacro expons(x) `(car ,x))

(defmacro coefs(x) `(cdr ,x))

(defvar mulbuf (make-array 2000 :element-type 'integer :initial-element 0))

(defun mul-denses (x y)  ; simplest case of multiplyng full polynomials
  (mul-densess x y 0 (1-(length (car x))) 0 (1- (length (car y)))))

(defun mul-densess (x y xbot xtop ybot ytop)  ; multiply sections of polynomials
  ;;bot and top are indexes into polys x and y
  ;; expons(x) is the array of exponents.
  ;; coefs(x) is the array of coefficients, same length as coefs(x).
  (declare (fixnum  xbot xtop ybot ytop)
;;	   (:explain :types)
	   (optimize (speed 3)(safety 0)))
  (let ((xlen (- xtop xbot)) ; how many terms from x are we using?
	(ylen (- ytop ybot)))
  (if (< xlen ylen) ; reverse args if x is shorter
      (mul-densess y x ybot ytop xbot xtop)
  
    (let* ((xe (expons x))			;array of exponents of x
	   (ye (expons y))			;array of exponents of y
	   (xc (coefs x))			;array of coefficients of x
	   (yc (coefs y))
	   (xei 0)(xci 0)
	   (maxXe (aref xe xtop)) ;maximum degree in x
	   (maxYe (aref ye ytop))
	   (minXe (aref xe xbot)) ;minimum degree in x
	   (minYe (aref ye ybot))
	   ;; allocate an array for all terms 
	   ;; with exponenets from minXe+minYe through maxXe+maxYe
	   (ASize (+ 1 (- maxXe minXe)(- maxYe minYe)))
	   ;; we could allocate this array from some buffer
	   (ans mulbuf)
)
      (declare 	(fixnum ASize maxXe maxYe minXe minYe)
		(type (array integer (*)) xc yc ans) 
		(type (array fixnum (*)) xe ye ))
      (if (< (length mulbuf) ASize) 
	  (setf ans (setf mulbuf (make-array ASize :element-type 'integer))))
       (loop for i fixnum from 0 below ASize      do (setf (aref ans i) 0))
      
      (loop for i fixnum from xbot 
	  to xtop 
	  finally (return (DA2PAs ans ASize (+ minXe minYe)))
	  do 
	    (setf xci (aref xc i))
	    ;; xci will not be zero by hypothesis
	    (setf xei (aref xe i))
	    (loop for j fixnum from ybot
		to ytop do
		  (incf (aref ans (+ xei (aref ye j)))
			(* xci (aref yc j)))))))))
				    

(defun DA2PAs (A size offset) 
    ;;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.
  ;; offset the exponenets by offset.
  (declare (type (array integer (*)) A)
	    (optimize (speed 3)(safety 0)))
  (let* (
	 (n
	  (count-if-not #'zerop A :start 0 :end size))
	 ;; on one example this line above takes 19% of the time!
	 ;(n (1- (length A))) faster but not simplest.
         (exps  (make-array n :element-type 'fixnum))
         (coefs (make-array n :element-type 'integer))
         (c 0))
    (declare (fixnum n)
	     (type (array integer (*)) coefs)
	     (type (array fixnum (*)) exps)
	     (optimize (speed 3)(safety 0)))
     (if (= 0 offset)
	 (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)))
       ;; non-zero offset
     (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 offset) (aref coefs j)c)
              (incf j))))))



;;; we need programs that add polys from bot to top,  shifted.
;;; yeech.
(defun foldpoly(x)(foldpolyat x (+ 1 (ash  (degr x) -1))))


(defun foldpolyat(x h)  ;; add the top and the bottom half of a polynomial.
  ;; for example,  A+x^n*B  where A,B are of degree < n --> A+B  
    (add-densess x x 0 h (1- h)  (1-(length (expons x))) (1- h)))


(defun add-densess (x y xbot xtop ybot ytop offset) ; add sections of polynomials
  ;;bot and top are indexes into polys x and y
  ;; expons(x) is the array of exponents.
  ;; coefs(x) is the array of coefficients, same length as coefs(x).
  (declare (fixnum  xbot xtop ybot ytop)
	   ;;	   (:explain :types)
	   (optimize (speed 3)(safety 0)))

  (let ((xe (expons x))			;array of exponents of x
	(ye (expons y))			;array of exponents of y
	(xc (coefs x))			;array of coefficients of x
	(yc (coefs y)))
    (cond ((>= xbot (length xe))
	   (return-from add-densess (polyleftshift y offset)))
	  ((>= ybot (length ye)) (return-from add-densess x))
	  ((>= xtop (length xe))(setf xtop (1- (length xe))))
	  ((>= ytop (length ye))(setf ytop (1- (length ye))))
	  ((> offset ybot) (error "illegal offset ~s in add-densess" offset)))
    (let*(
	  (xlen (- xtop xbot))	 ; how many terms from x are we using?
	  (ylen (- ytop ybot))
	  (maxXe (aref xe xtop))	;maximum degree in x
	  (maxYe (aref ye ytop))
	  (minXe (aref xe xbot))	;minimum degree in x
	  (minYe (aref ye ybot))
	  ;; allocate an array for all terms 
	  ;; with exponenets from minXe+minYe through maxXe+maxYe
	  (ASize (1+ (- (max maxXe maxYe)(min minXe minYe))))
	  ;; we could allocate this array from some buffer
	  (ans mulbuf))
      (declare 	(fixnum ASize maxXe maxYe minXe minYe)
		(type (array integer (*)) xc yc ans) 
		(type (array fixnum (*)) xe ye ))
					;      (print ASize)
      (if (< (length mulbuf) ASize) 
	  (setf ans (setf mulbuf (make-array ASize :element-type 'integer))))
      (loop for i fixnum from 0 below ASize  do (setf (aref ans i) 0))
      (loop for i fixnum from xbot below xtop  do (setf (aref ans i) (aref xc i)))
      (loop for j fixnum from ybot below ytop  do (incf (aref ans (- j offset)) (aref yc j)))
      (DA2PAs ans ASize 0)))
  )


#+ignore ;; not working
(defun mul-ks(pp qq) ;; by Karatsuba.,  use displaced arrays ;; sparse rep.
  (let* ((plim (length (car pp)))	; number of terms in p
	 (qlim (length (car qq))))	; number of terms in q
    (declare (fixnum qlim plim))
    (if (or (< plim *karatlim*)
	    (< qlim *karatlim*))
	(mul-dense pp qq) 

      ;; main program
      (let ((v (+ 1 (ash (max (degr pp)(degr qq)) -1)))) ;half of the  max degree of inputs, +1
	(multiple-value-bind
	    (A B)	  (truncatepoly pp v)
	  (setf B (kpolyshiftleft B (- v)))
	  (multiple-value-bind
	      (C D)    (truncatepoly qq v)
	    (setf D (kpolyshiftleft D  (- v)))
	    (let ((U  (mul-ks (foldpolyat pp v)  (foldpolyat qq v)))
		  (TT (mul-ks A C))
		  (S  (mul-ks B D))
		  (aa nil))
					;	    (format t "~%U=~s~%TT=~s~%S=~s" U TT S)
	      (setf U (kpolysub U S))
	      (setf U (kpolysub U TT))
					;	    (format t "~%U=~s" U)
	      (setf aa (kpolyaddshift TT U v))
	      (kpolyaddshift aa S  (+ v v))
	      )))))))


(defun tri(pp v)	(multiple-value-bind
	    (A B)	  (truncatepoly pp v)
	    (setf B (kpolyshiftleft B (- v)))
	    (kpolyadd A B)))

;; (try tt 5) should be same as (foldpolyat tt 5) ??

(defun tr2(b n)(declare (type (array integer (*)) b) (fixnum n)
		      (optimize (speed 3)(safety 0)))
       (loop for i fixnum from 0 below n
	   do (setf (aref b i) 0))
       b)

(defun tr3(b n)(declare (fixnum n) ;(type (array integer (*)) b)
		      (optimize (speed 3)(safety 0)))
       (make-array n :element-type 'integer :initial-element 0))

(defun tr4(b n)(declare (type (array t (*)) b) (fixnum n)
			(optimize (speed 3)(safety 0)))
  (loop for i fixnum from 0 below n
	do (setf (aref b i) 0))
  b)

(defun tr5(b n)(declare (simple-array b) (fixnum n)
			(optimize (speed 3)(safety 0)))
  (loop for i fixnum from 0 below n
	do (setf (aref b i) 0))
  b)

(defun tr6(b n)(declare ;(simple-array b)
			(fixnum n)
			(optimize (speed 3)(safety 0)))
  (loop for i fixnum from 0 below n
	do (setf (aref b i) 0))
  b)

(defun tr7(b n)(declare (simple-array  b)
			(fixnum n)
			(optimize (speed 3)(safety 0)))
  (loop for i fixnum from 0 below n
	do (setf (svref b i) 0))
  b)


(defun tr2b(b n)(declare (type (array integer (*)) b) (fixnum n)
			 (optimize (speed 3)(safety 0)))
       (loop for i fixnum from 0 below n
	   do (setf (svref b i) 0))  ;;; the winner!!!
       b) ;; 35 bytes


(defun tr2(b n)(declare (type (array integer (*)) b) (fixnum n)
		      (optimize (speed 3)(safety 0)))
       (loop for i fixnum from 0 below n
	   do (setf (aref b i) 0))
       b);; 59 bytes



(defun tr2c(b n)(declare (type (simple-array fixnum (*)) b) (fixnum n)
;;			 (:explain :types :inlining)
		      (optimize (speed 3)(safety 0)))
       (loop for i fixnum from 0 below n
	   do (setf (aref b i) 0))
       b);; 35 bytes


(defun tr2d(b n)(declare (type (simple-array t (*)) b) (fixnum n)
			 ;(:explain :types :inlining)
		      (optimize (speed 3)(safety 0)))
       (loop for i fixnum from 0 below n
	   do (setf (aref b i) 0))
       b)

(defun tr2e(b n)(declare (type (simple-array integer (*)) b) (fixnum n)
			 ;(:explain :types :inlining)
		      (optimize (speed 3)(safety 0)))
       (loop for i fixnum from 0 below n
	   do (setf (aref b i) 0))
       b)


;(defvar mulbuf (make-array 2000 :element-type 'integer :initial-element 0))
;; simple-array
(setf mulbufsa (make-array 2000  :initial-element 0)) ;; simple array not 'integer

;(print 'bye)


;;  (time (dotimes (i 1000000)(tr3 mulbuf 1000)))
;;  (time (dotimes (i 1000000)(tr4 mulbufsa 1000)))
;;  (time (dotimes (i 1000000) (tr2 mulbuf 1000)))