;;; -*- Lisp -*-

;;; naive multiplication of polys using bit vector
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))

;; use this only if the highest degree is not extravagant
;; may win on a kind of middling-sparse domain, but maybe not.
;; better than naive, sometimes.  Better than dense??? Where??
;; better than the best of dense or heap5? no..

;; mul-naive-b seems to be the best of the programs in this file.

 ;; works just fine. try to speed up in -bb
(defun mul-naive-b (x y);; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (cond 
     
   #+ignore  ((> hideg 1000000)
      (format t  "~%mul-naive-b can't be used: too high degree ~s, using heap" hideg)
      (mul-heap5 x y))
   #+ignore  ((> (* 20 hideg) (*  xl yl))
      (format t  "~%mul-naive-b shouldn't be used: too low degree ~s, using dense" hideg)
      (mul-dense x y))
   
     (t(let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	    (xc (cdr x))  (yc (cdr y))
	    (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	    (termcount 0))
	 (declare(fixnum xei j ee termcount)
		 (type (simple-array t (*)) xc yc ansy ans)
		 (type (simple-array bit (*)) ansbit))
	 ;; compute exponent "skeleton"
	 ;; keeping track of maximum number of non-zero coefs possible
	 ;; in termcount.

	 (dotimes (i xl)
	   (declare (fixnum i))
	   (setf xei (aref xe i))
	   (dotimes(j yl)
	     (declare (fixnum j))
	     (setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
	     (unless (eq 1 (sbit ansbit ee))
	       (setf (sbit ansbit ee) 1)
	       (incf termcount))) )
	   
	 ;; now we know exactly how many terms in answer, allocate storage.
	 ;; we could do this:  if the number of terms is nearly the same as hideg,
	 ;; just use mul-dense.
	 (setf ae (make-array termcount))
	 (setf ac (make-array termcount :initial-element 0))
	 ;; initialize exponents. j=0 to start
	 (dotimes (i (1+ hideg) ae)
	   (declare(fixnum i))
	   (unless (zerop (aref ansbit i))
	     (setf (aref ae j) i)
	     (incf j)))
      
	 (dotimes (i xl (cons ae ac))
	   (declare (fixnum i))
	   (setf xei (aref xe i) xci (aref xc i))
	   (dotimes(j yl)
	     ;; multiply terms for real, this time.
	     (declare (fixnum j))

	     (binsearch-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			       (* (aref yc j) xci)))))))))

;; works. faster than mul-naive-b by 30% or more
(defun mul-naive-bb (x y);; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (cond 
     
   #+ignore  ((> hideg 1000000)
      (format t  "~%mul-naive-b can't be used: too high degree ~s, using heap" hideg)
      (mul-heap5 x y))
   #+ignore  ((> (* 20 hideg) (*  xl yl))
      (format t  "~%mul-naive-b shouldn't be used: too low degree ~s, using dense" hideg)
      (mul-dense x y))
   
     (t(let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	    (xc (cdr x))  (yc (cdr y)) (low 0)
	    (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	    (termcount 0))
	 (declare(fixnum xei j ee termcount low)
		 (type (simple-array t (*)) xc yc ansy ans)
		 (type (simple-array bit (*)) ansbit))
	 ;; compute exponent "skeleton"
	 ;; keeping track of maximum number of non-zero coefs possible
	 ;; in termcount.

	 (dotimes (i xl)
	   (declare (fixnum i))
	   (setf xei (aref xe i))
	   (dotimes(j yl)
	     (declare (fixnum j))
	     (setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
	     (unless (eq 1 (sbit ansbit ee))
	       (setf (sbit ansbit ee) 1)
	       (incf termcount))) )
	   
	 ;; now we know exactly how many terms in answer, allocate storage.
	 ;; we could do this:  if the number of terms is nearly the same as hideg,
	 ;; just use mul-dense.
	 (setf ae (make-array termcount))
	 (setf ac (make-array termcount :initial-element 0))
	 ;; initialize exponents. j=0 to start
	 (dotimes (i (1+ hideg) ae)
	   (declare(fixnum i))
	   (unless (zerop (aref ansbit i))
	     (setf (aref ae j) i)
	     (incf j)))
      
	 (dotimes (i xl (cons ae ac))
	   (declare (fixnum i))
	   (setf xei (aref xe i) xci (aref xc i))
	   (setf low -1)
	   (dotimes(j yl)
	     ;; multiply terms for real, this time.
	     (declare (fixnum j))
	     (setf low
	       (binsearch-updateL ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
				  (* (aref yc j) xci) (1+ low) termcount)))))))))


(defun mul-naive-bbb (x y);; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (yemax  (aref ye (1- yl)))
	 (hideg (+ (aref xe (1- xl))yemax ));degree of answer
	 )
    (declare(fixnum yl xl hideg yemax)
	    (type (simple-array t (*)) xe ye))
    (cond 
     
   #+ignore  ((> hideg 1000000)
      (format t  "~%mul-naive-b can't be used: too high degree ~s, using heap" hideg)
      (mul-heap5 x y))
   #+ignore  ((> (* 20 hideg) (*  xl yl))
      (format t  "~%mul-naive-b shouldn't be used: too low degree ~s, using dense" hideg)
      (mul-dense x y))
   
     (t(let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	    (xc (cdr x))  (yc (cdr y)) (low 0)(high 0)
	    (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	    (termcount 0))
	 (declare(fixnum xei j ee termcount low high yemax)
		 (type (simple-array t (*)) xc yc ansy ans)
		 (type (simple-array bit (*)) ansbit))
	 ;; compute exponent "skeleton"
	 ;; keeping track of maximum number of non-zero coefs possible
	 ;; in termcount.

	 (dotimes (i xl)
	   (declare (fixnum i))
	   (setf xei (aref xe i))
	   (dotimes(j yl)
	     (declare (fixnum j))
	     (setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
	     (unless (eq 1 (sbit ansbit ee))
	       (setf (sbit ansbit ee) 1)
	       (incf termcount))) )
	   
	 ;; now we know exactly how many terms in answer, allocate storage.
	 ;; we could do this:  if the number of terms is nearly the same as hideg,
	 ;; just use mul-dense.
	 (setf ae (make-array termcount))
	 (setf ac (make-array termcount :initial-element 0))
	 ;; initialize exponents. j=0 to start
	 (dotimes (i (1+ hideg) ae)
	   (declare(fixnum i))
	   (unless (zerop (aref ansbit i))
	     (setf (aref ae j) i)
	     (incf j)))
	 ;(format t "~%termcount=~s" termcount)
	 (dotimes (i xl (cons ae ac))
	   (declare (fixnum i))
	   (setf xei (aref xe i) xci (aref xc i))
	   (setf low -1)
	   (setf high  (1+ (binsearch-updateL ae ac (+ ;(aref xe (1- xl))
						  xei
						  yemax) 0 0  termcount)))
	   (dotimes (j yl) 
	     (declare (fixnum j))
	     (setf low
	       (binsearch-updateL ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
				  (* (aref yc j) xci) (1+ low) high)))))))))




(defun binsearch-update (indexA valA  ind val)
  ;;  binary search and update. All exponents are in array indexA, all values in array valA
  ;; find j such that  ind=indexA[j]. Then add val to element valA[j]   
  ;; by construction there is such a j.  Also indexA is sorted.
  (declare (fixnum ind)
	   (type (simple-array t (*)) indexA valA))

  (let ((lo 0)
	(hi (length indexA))
	(p 0)
	(midpt 0))
    (declare (fixnum lo hi ind p midpt count))
    (loop
      (setf p (ash (- hi lo) -1))
      (if (eq p 0)
	  ;; all values are in the array, so we found it.
	  (return(incf (aref valA 0) val)))  ;; 0? or lo
      (setf p (+ lo p))
      (setf midpt (the fixnum (aref indexA p)))	; midpoint, about.
      (cond ((eq midpt ind)
	     (return(incf (aref valA p) val))) ;; the only other way out
	    ((< ind midpt) 
	     (setf hi p))
	    (t (setf lo p))))))



(defun binsearch-updateL (indexA valA  ind val lo hi)
  ;;  binary search and update. All exponents are in array indexA, all values in array valA
  ;; find j such that  ind=indexA[j]. Then add val to element valA[j]   
  ;; by construction there is such a j.  Also indexA is sorted.
  (declare (fixnum ind lo hi)
	   (type (simple-array t (*)) indexA valA))
  (let ((p 0)
	(midpt 0))
    (declare (fixnum p midpt))
   ; (setf lo 0)
   ; (setf hi (length indexA))
    (loop
      (setf p (ash (- hi lo) -1))
      (cond((eq p 0) ;; either lo=hi or lo=hi-1
	    ;; all values are in the array, so we found it.
	    (incf (aref valA lo) val)
	 ;   (prin1 ".")
					;  (return lo)
	    (return 0)))
      (setf p (+ lo p))
      (setf midpt (the fixnum (aref indexA p)))	; midpoint, about.
      (cond ((eq midpt ind)
	     (incf (aref valA p) val)
	     (return p)) ;; the only other way out
	    ((< ind midpt)
	     (setf hi p))
	    (t (setf lo p))))))
	  


;; use logical operations to count bits.
;; this does not work if the mask integers get too large
;; does only O(n)+O(m) operations to find out where all the exponents are.

(defun mul-naive-c (x y) ;; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
    (let* ((xe (car x))
	   (ye (car y))
	   (yl (length (the simple-array ye)))
	   (xl (length (the simple-array xe)))
	   (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	   (xc (cdr x))
	   (yc (cdr y))
	   (j 0)
	   (xei 0)
	   (xci 0)
	   (ae #())
	   (ac #())
	   (maskx 0)
	   (maskxy 0)
	   (termcount 0))
      (declare(fixnum yl xl xei j termcount)
	      (type (simple-array t (*)) xe ye xc yc ))
      (if (> hideg 1000000) (error "mul-naive-c can't be used: too high degree ~s" hideg))

      (dotimes (i xl)
	(setf maskx (logior maskx (ash 1 (aref xe i))))) ;; place exponents into mask
	
      (dotimes(j yl)
	(setf maskxy (logior maskxy (ash maskx (the fixnum (aref ye j))))))
      (setf termcount (logcount maskxy))
      ;; now we know exactly how many terms in answer, allocate storage.
      ;; we could do this:  if the number of terms is nearly the same as hideg,
      ;; just use mul-dense.
      
      (setf ae (make-array termcount))
      (setf ac (make-array termcount :initial-element 0))
      (setf j 0)
      ;; initialize exponents
      (dotimes (i (1+ hideg))
	(cond ((logbitp i maskxy)
	       (setf (aref ae j) i)
	       (incf j))))
   ;   (format t "~%hideg=~s maskxy=~s ~% ~%ae=~s~%ac=~s" hideg maskxy ae ac)
      (dotimes (i xl (cons ae ac))
		  (setf xei (aref xe i) xci (aref xc i))

	(dotimes(j yl)
	  ;; multiply terms for real, this time.
	  (binsearch-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			    (* (aref yc j) xci))))))


(defun mul-naive-d (x y) ;; SAME AS -c  except interp-update
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
    (let* ((xe (car x))
	   (ye (car y))
	   (yl (length (the simple-array ye)))
	   (xl (length (the simple-array xe)))
	   (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	   (xc (cdr x))
	   (yc (cdr y))
	   (j 0)
	   (xei 0)
	   (xci 0)
	   (ae #())
	   (ac #())
	   (maskx 0)
	   (maskxy 0)
	   (termcount 0))
      (declare(fixnum yl xl xei j termcount)
	      (type (simple-array t (*)) xe ye xc yc ))
      (if (> hideg 1000000) (error "mul-naive-c can't be used: too high degree ~s" hideg))

      (dotimes (i xl)
	(setf maskx (logior maskx (ash 1 (aref xe i))))) ;; place exponents into mask
	
      (dotimes(j yl)
	(setf maskxy (logior maskxy (ash maskx (the fixnum (aref ye j))))))
      (setf termcount (logcount maskxy))
      ;; now we know exactly how many terms in answer, allocate storage.
      ;; we could do this:  if the number of terms is nearly the same as hideg,
      ;; just use mul-dense.
      
      (setf ae (make-array termcount))
      (setf ac (make-array termcount :initial-element 0))
      (setf j 0)
      ;; initialize exponents
      (dotimes (i (1+ hideg))
	(cond ((logbitp i maskxy)
	       (setf (aref ae j) i)
	       (incf j))))
   ;   (format t "~%hideg=~s maskxy=~s ~% ~%ae=~s~%ac=~s" hideg maskxy ae ac)
      (dotimes (i xl (cons ae ac))
		  (setf xei (aref xe i) xci (aref xc i))
	(dotimes(j yl)
	  ;; multiply terms for real, this time.

	  (interp-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			    (* (aref yc j) xci))))))


#|
 (setf q (make-racg 5000 100 5)r (make-racg 1000 100 5))
[3] cl-user(195): (start-profiler)
:sampling
[3] cl-user(196): (dotimes (i 5)(mul-naive-c q r))
nil
[3] cl-user(197): (show-flat-profile)
Sampling stopped after 6694 samples taken for the time profiler.

Sample represents 13.1 seconds of processor time (out of a total of 13.1)

Times below 1.0% will be suppressed.

  %     %     self  total            self   total  Function
 Time  Cum.   secs   secs    calls ms/call ms/call   name
 61.9  61.9    8.1    9.0                          ... binsearch-update
  9.3  71.3    1.2    1.2                          "logical_or"
  7.5  78.7    1.0    1.0                          "shift_left"
  5.6  84.3    0.7   12.7                          ... mul-naive-c
  5.4  89.7    0.7    0.7                          excl::.inv-aref_1d
  3.1  92.9    0.4    0.4                          excl::aref_1d
  1.9  94.8    0.2    0.2                          "integer_multiply"
  1.9  96.6    0.2    0.2                          excl::*_2op
  1.0  97.7    0.1    0.1                          "where_am_i_2"
[3] cl-user(198): (start-profiler)
:sampling
[3] cl-user(199): (dotimes (i 5)(mul-naive-b q r))
nil
[3] cl-user(200): (show-flat-profile)
Sampling stopped after 5962 samples taken for the time profiler.

Sample represents 11.6 seconds of processor time (out of a total of 11.6)

Times below 1.0% will be suppressed.

  %     %     self  total            self   total  Function
 Time  Cum.   secs   secs    calls ms/call ms/call   name
 76.7  76.7    8.9    9.8                          ... binsearch-update
  7.6  84.4    0.9   11.2                          ... mul-naive-b
  4.7  89.0    0.5    0.5                          excl::.inv-aref_1d
  3.3  92.3    0.4    0.4                          excl::aref_1d
  2.2  94.6    0.3    0.3                          "integer_multiply"
  2.1  96.7    0.2    0.2                          excl::*_2op
  1.1  97.8    0.1    0.1                          length
  1.1  98.9    0.1    0.1                          "where_am_i_2"
  [3] cl-user(201): 


  cl-user(7):   (setf q (make-racg 5000 100 5)r (make-racg 1000 100 5))
(#(72 137 189 281 284 357 390 403 452 475 ...)
 . #(5 3 4 3 5 3 2 3 4 3 ...))
cl-user(8): (time (mul-dense q r))
; cpu time (non-gc) 250 msec user, 0 msec system
; cpu time (gc)     47 msec user, 0 msec system
; cpu time (total)  297 msec user, 0 msec system
; real time  297 msec
; space allocation:
;  5 cons cells, 3,604,792 other bytes, 0 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(9): (time (mul-heap5 q r))
; cpu time (non-gc) 6,233 msec user, 0 msec system
; cpu time (gc)     63 msec user, 0 msec system
; cpu time (total)  6,296 msec user, 0 msec system
; real time  6,453 msec
; space allocation:
;  223 cons cells, 5,288,464 other bytes, 0 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(10): (time (mul-rtree4 q r))
; cpu time (non-gc) 2,829 msec user, 31 msec system
; cpu time (gc)     265 msec user, 16 msec system
; cpu time (total)  3,094 msec user, 47 msec system
; real time  3,266 msec
; space allocation:
;  372,816 cons cells, 10,019,784 other bytes, 0 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(11): (time (mul-naive-c q r))
; cpu time (non-gc) 2,924 msec user, 16 msec system
; cpu time (gc)     264 msec user, 0 msec system
; cpu time (total)  3,188 msec user, 16 msec system
; real time  3,468 msec
; space allocation:
;  1,043 cons cells, 230,349,328 other bytes, 0 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(12): (time (mul-naive-b q r))
; cpu time (non-gc) 2,422 msec user, 0 msec system
; cpu time (gc)     0 msec user, 0 msec system
; cpu time (total)  2,422 msec user, 0 msec system
; real time  2,438 msec
; space allocation:
;  7 cons cells, 2,427,784 other bytes, 0 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(13): (time (mul-naive q r))
; cpu time (non-gc) 16,093 msec user, 0 msec system
; cpu time (gc)     157 msec user, 15 msec system
; cpu time (total)  16,250 msec user, 15 msec system
; real time  16,718 msec
; space allocation:
;  597,461 cons cells, 2,389,888 other bytes, 0 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(14): (time (mul-hash1 q r))
; cpu time (non-gc) 5,392 msec user, 16 msec system
; cpu time (gc)     296 msec user, 0 msec system
; cpu time (total)  5,688 msec user, 16 msec system
; real time  5,828 msec
; space allocation:
;  21 cons cells, 21,513,776 other bytes, 8,192 static bytes
(#(127 157 192 211 222 244 274 276 281 308 ...)
 . #(20 15 12 5 9 16 12 3 15 25 ...))
cl-user(15): 

  |#


#|  public int interpolationSearch(int[] sortedArray, int toFind){
  // Returns index of toFind in sortedArray, or -1 if not found
  int low = 0;
  int high = sortedArray.length - 1;
  int mid;
 
  while (sortedArray[low] < toFind && sortedArray[high] >= toFind) {
   mid = low + ((toFind - sortedArray[low]) * (high - low)) / (sortedArray[high] - sortedArray[low]);
 
   if (sortedArray[mid] < toFind)
    low = mid + 1;
   else if (sortedArray[mid] > toFind)      //Repetition of the comparison code is forced by syntax limitations.
    high = mid - 1;
   else
    return mid;
  }
 
  if (sortedArray[low] == toFind)
   return low;
  else
   return -1; // Not found
 }
|#


(defun iS    (sortedArray  toFind)
    ;;interpolationSearch
  (let ((low 0)
	(high (1-(length sortedArray))) 
	(mid 0)
	(sal 0)
	(sah 0)
	(sam 0)
	)
    (declare (fixnum sal sah mid high low sam))
    
    (while  (and (< (setf sal(aref sortedArray low)) toFind)
		 (>= (setf sah (aref sortedArray high)) toFind))
      
      (setf mid (+ low (truncate(* (- toFind sal)(- high low))
				(- sah sal))))
      (if (< (setf sam(aref sortedArray mid)) toFind) (setf low (+ mid 1))
	(if (> sam toFind) (setf high (- mid 1))
	  (return-from iS mid))))
    (if (= sal toFind)  (return-from iS low) (error "not found"))))

;; declared stuff..
;; THIS DOES NOT WORK.  See next version later
#+ignore
(defun interp-update    (sortedArray valA  toFind val)
    ;;interpolationSearch and update.
  (let ((low 0)
	(high (1-(length sortedArray))) 
	(mid 0)
	(sal 0)
	(sah 0)
	(sam 0)
	)
    (declare (fixnum sal sah mid high low sam toFind)
	      (type (simple-array t (*)) sortedArray valA))
    
    (while  (and (< (setf sal(aref sortedArray low)) toFind)
		 (>= (setf sah (aref sortedArray high)) toFind))
      
      (setf mid (+ low (the fixnum 
			 (truncate(the fixnum (* (the fixnum(- toFind sal))(the fixnum(- high low))))
				  (the fixnum (- sah sal))))))
      
      (cond ((< (setf sam(aref sortedArray mid)) toFind) (setf low (+ mid 1)))
	    ((> sam toFind)(setf high (- mid 1)))
	    (t (incf (aref valA mid) val)
	       (return-from interp-update mid)))
      )
    ;; it is in there, by construction, so if we get here, we have found it.
     (incf (aref valA low) val)
     (return-from interp-update low)))

;; another version, works

(defun interp-update    (sortedArray valA  toFind val)
  ;;interpolationSearch and update.
  (declare (fixnum toFind)
	    (type (simple-array t (*)) sortedArray valA))
  (let* ((low 0)
	(high (1-(length sortedArray))) 
	(mid 0)
	(sal (aref sortedArray low))
	(sah (aref sortedArray high))
	(sam 0))
    (declare (fixnum sal sah mid high low sam toFind)
	  ;;   (:explain :types :calls)
	    )
    
    (while  (and (< sal toFind)
		 (>= sah toFind))

    (setf mid (+ low (truncate (* (- toFind sal)
				    (- high low))
					 (- sah sal))))
 
 #+ignore ;; fails to work!
 (setf mid (+ low (truncate (the fixnum (* (- toFind sal)
				    (- high low)))
				 (the fixnum (- sah sal)))))
      
      

      
      (setf sam (aref sortedArray mid))
      (cond ((<  sam toFind) 
	     (setf low (+ mid 1))
	     (setf sal (aref sortedArray low)))
	    ((> sam toFind)
	     (setf high (- mid 1))
	     (setf sah (aref sortedArray high)))
	    (t (incf (aref valA mid) val)
	       (return-from interp-update mid))))
    ;; it is in there, by construction, so if we get here, we have found it.
    (incf (aref valA low) val)
    low))


(defun mul-naive-c2 (x y) ;; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
    (let* ((xe (car x))
	   (ye (car y))
	   (yl (length (the simple-array ye)))
	   (xl (length (the simple-array xe)))
	   (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	   (xc (cdr x))
	   (yc (cdr y))
	   (j 0)
	   (xei 0)
	   (xci 0)
	   (ae #())
	   (ac #())
	   (maskx 0)
	   (maskxy 0)
	   (termcount 0)
	   (bytefield (byte 1 0)))
      (declare(fixnum yl xl xei j termcount)
	      (dynamic-extent bytefield maskx maskxy)

	      (type (simple-array t (*)) xe ye xc yc ))
      (if (> hideg 1000000) (error "mul-naive-c can't be used: too high degree ~s" hideg))

      (dotimes (i xl);;
		;;(setf (byte-position bytefield) (aref xe i))
		(setf (cdr bytefield) (aref xe i))
		(setf (ldb bytefield maskx) 1)
		)      ;; place exponents into mask

      #+ignore 
      ;; this is slower! I thought building the integer maskx first would help. 
      ;; it even takes more memory.
      
   (do ((i (1- xl) (1- i)))
	  ((< i 0) nil)
	(declare (fixnum i))
	(setf (cdr bytefield) (aref xe i)) ;byte-position
	(setf (ldb bytefield maskx) 1))

	
      (dotimes(j yl)
	(setf maskxy (logior maskxy (ash maskx (the fixnum (aref ye j))))))
      (setf termcount (logcount maskxy))
      ;; now we know exactly how many terms in answer, allocate storage.
      ;; we could do this:  if the number of terms is nearly the same as hideg,
      ;; just use mul-dense.
      
      (setf ae (make-array termcount))
      (setf ac (make-array termcount :initial-element 0))
      (setf j 0)
      ;; initialize exponents
      (dotimes (i (1+ hideg))
	(cond ((logbitp i maskxy)
	       (setf (aref ae j) i)
	       (incf j))))
   ;   (format t "~%hideg=~s maskxy=~s ~% ~%ae=~s~%ac=~s" hideg maskxy ae ac)
      (dotimes (i xl (cons ae ac))
		  (setf xei (aref xe i) xci (aref xc i))

	(dotimes(j yl)
	  ;; multiply terms for real, this time.
	  (binsearch-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			    (* (aref yc j) xci))))))


;; slower because interp-update does too much work.

(defun mul-naive-c2i (x y) ;; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
    (let* ((xe (car x))
	   (ye (car y))
	   (yl (length (the simple-array ye)))
	   (xl (length (the simple-array xe)))
	   (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	   (xc (cdr x))
	   (yc (cdr y))
	   (j 0)
	   (xei 0)
	   (xci 0)
	   (ae #())
	   (ac #())
	   (maskx 0)
	   (maskxy 0)
	   (termcount 0)
	   (bytefield (cons 1 2)))
      (declare(fixnum yl xl xei j termcount)
	      (dynamic-extent bytefield maskx maskxy)
	      (type (simple-array t (*)) xe ye xc yc ))
      (if (> hideg 1000000) (error "mul-naive-c can't be used: too high degree ~s" hideg))

      (dotimes (i xl);;
	;;(setf maskx (logior maskx (ash 1 (aref xe i))))
	(setf (cdr bytefield) (aref xe i))
	(setf (ldb bytefield maskx) 1)
	) ;; place exponents into mask
	
      (dotimes(j yl)
	(setf maskxy (logior maskxy (ash maskx (the fixnum (aref ye j)))))
		   )
      (setf termcount (logcount maskxy))
      ;; now we know exactly how many terms in answer, allocate storage.
      ;; we could do this:  if the number of terms is nearly the same as hideg,
      ;; just use mul-dense.
      
      (setf ae (make-array termcount))
      (setf ac (make-array termcount :initial-element 0))
      (setf j 0)
      ;; initialize exponents
      (dotimes (i (1+ hideg))
	(cond ((logbitp i maskxy)
	       (setf (aref ae j) i)
	       (incf j))))
   ;   (format t "~%hideg=~s maskxy=~s ~% ~%ae=~s~%ac=~s" hideg maskxy ae ac)
      (dotimes (i xl (cons ae ac))
		  (setf xei (aref xe i) xci (aref xc i))
	(dotimes(j yl)
	  ;; multiply terms for real, this time.
	  (interp-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			 (* (aref yc j) xci))))))



;; compiler testing

#+ignore
(defun mnb (x y);; arrange the data structure by computing a skeleton of answer exponents
  (declare (optimize (speed 3)(safety 0)))
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	 (xc (cdr x))  (yc (cdr y))
	 (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	 (termcount 0))
      (declare(fixnum xei j ee te>rmcount)
	      (type (simple-array t (*)) xc yc ansy ans)
	      (type (simple-array bit (*)) ansbit))
      ;; compute exponent "skeleton"
      ;; keeping track of maximum number of non-zero coefs possible
      ;; in termcount.

      (dotimes (i xl termcount)
	(declare (fixnum i))
	(setf xei (aref xe i))
	(dotimes(j yl)(declare (fixnum j))
		(setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
		(unless (eq 1 (aref ansbit ee))
		  (setf (aref ansbit ee) 1)
		  (incf termcount))) ))))


#+ignore
(defun mnb2 (x y);; arrange the data structure by computing a skeleton of answer exponents
  (declare (optimize (speed 3)(safety 0)))
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	 (xc (cdr x))  (yc (cdr y))
	 (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	 (termcount 0))
      (declare(fixnum xei j ee te>rmcount)
	      (type (simple-array t (*)) xc yc ansy ans)
	      (type (simple-array bit (*)) ansbit))
      ;; compute exponent "skeleton"
      ;; keeping track of maximum number of non-zero coefs possible
      ;; in termcount.

      (dotimes (i xl termcount)
	(declare (fixnum i))
	(setf xei (sbit xe i))
	(dotimes(j yl)(declare (fixnum j))
		(setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
		(unless (eq 1 (sbit ansbit ee))
		  (setf (sbit ansbit ee) 1)
		  (incf termcount))) ))))

;; this works but is way too slow . assumes toFind is a double-float
(defun interp-updatef    (sortedArray valA  toFind val)
  ;;interpolationSearch and update. using floats.
  (declare (double-float toFind)
	    (type (simple-array t (*)) sortedArray valA))
  (let* ((low 0)
	 (high (1-(length sortedArray))) 
	 (mid 0)
	 (sal (aref sortedArray low)); double
	 (sah  (aref sortedArray high)) ;double
	 (sam 0.0d0)
	)

    (declare (fixnum mid high low)
	     (double-float sal sah sam toFind) )
    
    (while  (and (< sal toFind)
		 (>= sah toFind))

      (setf mid (+ low			;fixnum
		   (truncate (* (- toFind sal) ;float
				(- high low)) ;fix
			     (- sah sal) ;float
			     )))
      
#+ignore ;; does not work
  (setf mid (+ low	(* (- high low)	;fixnum
		   (truncate (- toFind sal) ;float
			     (- sah sal) ;float
			     ))))
 
      (setf sam (coerce (aref sortedArray mid) 'double-float))
      (cond ((<  sam toFind) 
	     (setf low (+ mid 1))
	     (setf sal (coerce (aref sortedArray low) 'double-float)))
	    ((> sam toFind)
	     (setf high (- mid 1))
	     (setf sah (coerce (aref sortedArray high) 'double-float)))
	    (t (incf (aref valA mid) val)
	       (return-from interp-updatef mid))))
    ;; it is in there, by construction, so if we get here, we have found it.
    (incf (aref valA low) val)
    low))
(defun mul-naive-bi (x y);; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (cond 
   #+ignore  ((> hideg 1000000)
      (format t  "~%mul-naive-b can't be used: too high degree ~s, using heap" hideg)
      (mul-heap5 x y))
   #+ignore  ((> (* 20 hideg) (*  xl yl))
      (format t  "~%mul-naive-b shouldn't be used: too low degree ~s, using dense" hideg)
      (mul-dense x y))
     (t(let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	    (xc (cdr x))  (yc (cdr y))
	    (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	    (termcount 0))
	 (declare(fixnum xei j ee termcount)
		 (type (simple-array t (*)) xc yc ansy ans)
		 (type (simple-array bit (*)) ansbit))
	 ;; compute exponent "skeleton"
	 ;; keeping track of maximum number of non-zero coefs possible
	 ;; in termcount.

	 (dotimes (i xl)
	   (declare (fixnum i))
	   (setf xei (aref xe i))
	   (dotimes(j yl)
	     (declare (fixnum j))
	     (setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
	     (unless (eq 1 (sbit ansbit ee))
	       (setf (sbit ansbit ee) 1)
	       (incf termcount))) )
	   
	 ;; now we know exactly how many terms in answer, allocate storage.
	 ;; we could do this:  if the number of terms is nearly the same as hideg,
	 ;; just use mul-dense.
	 (setf ae (make-array termcount))
	 (setf ac (make-array termcount :initial-element 0))
	 ;; initialize exponents. j=0 to start
	 (dotimes (i (1+ hideg) ae)
	   (declare(fixnum i))
	   (unless (zerop (aref ansbit i))
	     (setf (aref ae j) i)
	     (incf j)))
      
	 (dotimes (i xl (cons ae ac))
	   (declare (fixnum i))
	   (setf xei (aref xe i) xci (aref xc i))
	   (dotimes(j yl)
	     ;; multiply terms for real, this time.
	     (declare (fixnum j))

	     (interp-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			     (* (aref yc j) xci)))))))))


;; use interp, float. not a winner, but try double-float exponents...
(defun mul-naive-bif (x y);; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (cond 
   #+ignore  ((> hideg 1000000)
      (format t  "~%mul-naive-b can't be used: too high degree ~s, using heap" hideg)
      (mul-heap5 x y))
   #+ignore  ((> (* 20 hideg) (*  xl yl))
      (format t  "~%mul-naive-b shouldn't be used: too low degree ~s, using dense" hideg)
      (mul-dense x y))
     (t(let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	    (xc (cdr x))  (yc (cdr y))
	    (j 0)  (xei 0)(fei 0.0d0)  (xci 0)  (ee 0) (ae #())  (ac #())
	    (termcount 0))
	 (declare(fixnum  j ee termcount xei)(double-float fei)
		 (type (simple-array t (*)) xc yc ansy ans)
		 (type (simple-array bit (*)) ansbit))
	 ;; compute exponent "skeleton"
	 ;; keeping track of maximum number of non-zero coefs possible
	 ;; in termcount.

	 (dotimes (i xl)
	   (declare (fixnum i))
	   (setf xei (aref xe i))
	   (dotimes(j yl)
	     (declare (fixnum j))
	     (setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
	     (unless (eq 1 (sbit ansbit ee))
	       (setf (sbit ansbit ee) 1)
	       (incf termcount))) )
	   
	 ;; now we know exactly how many terms in answer, allocate storage.
	 ;; we could do this:  if the number of terms is nearly the same as hideg,
	 ;; just use mul-dense.
	 (setf ae (make-array termcount))
	 (setf ac (make-array termcount :initial-element 0))
	 ;; initialize exponents. j=0 to start
	 (dotimes (i (1+ hideg) ae)
	   (declare(fixnum i))
	   (unless (zerop (aref ansbit i))
	     (setf (aref ae j) (coerce i 'double-float))
	     (incf j)))
      
	 (dotimes (i xl (cons ae ac))
	   (declare (fixnum i))
	   (setf fei (coerce (aref xe i) 'double-float) xci (aref xc i))
	   (dotimes(j yl)
	     ;; multiply terms for real, this time.
	     (declare (fixnum j))

	     (interp-updatef ae ac  (+ fei (aref ye j)) ; a double float
			       (* (aref yc j) xci)))))))))





(defun linsearch-update (indexA valA  ind val guess) ; guess could be (ash (length indexA) -1)
  ;;  linear search and update. All exponents are in array indexA, all values in array valA
  ;; find j such that  ind=indexA[j]. Then add val to element valA[j]   
  ;; by construction there is such a j.  Also indexA is sorted.
  (declare (fixnum ind guess)
	   (type (simple-array t (*)) indexA valA))
  (let ((p guess)
	(thisind 0))
    (declare (fixnum p thisind))
    (loop
      (setf thisind (aref indexA p))
      (cond ((<=  ind thisind)
	     (cond ((eq thisind ind)
		    (incf (aref valA p) val)
		    (return p))
		   (t(decf p))))
	    (t (incf p))))))


;; use linear search.  Really a loser.
(defun mul-naive-bl (x y);; arrange the data structure by computing a skeleton of answer exponents
  
  ;; if hideg (computed below) is huge, we may not want to do this at all.
  (let* ((xe (car x))
	 (ye (car y))
	 (yl (length (the simple-array ye)))
	 (xl (length (the simple-array xe)))
	 (hideg (+ (aref xe (1- xl))(aref ye (1- yl))));degree of answer
	 )
    (declare(fixnum yl xl hideg)
	    (type (simple-array t (*)) xe ye))
    (cond 
     
   #+ignore  ((> hideg 1000000)
      (format t  "~%mul-naive-b can't be used: too high degree ~s, using heap" hideg)
      (mul-heap5 x y))
   #+ignore  ((> (* 20 hideg) (*  xl yl))
      (format t  "~%mul-naive-b shouldn't be used: too low degree ~s, using dense" hideg)
      (mul-dense x y))
   
     (t(let((ansbit (make-array (1+ hideg) :element-type 'bit :initial-element 0))
	    (xc (cdr x))  (yc (cdr y)) (guess 0)
	    (j 0)  (xei 0)  (xci 0)  (ee 0) (ae #())  (ac #())
	    (termcount 0))
	 (declare(fixnum xei j ee termcount guess)
		 (type (simple-array t (*)) xc yc ansy ans)
		 (type (simple-array bit (*)) ansbit))
	 ;; compute exponent "skeleton"
	 ;; keeping track of maximum number of non-zero coefs possible
	 ;; in termcount.

	 (dotimes (i xl)
	   (declare (fixnum i))
	   (setf xei (aref xe i))
	   (dotimes(j yl)
	     (declare (fixnum j))
	     (setf ee (the fixnum(+ xei(the fixnum (aref ye j)))))
	     (unless (eq 1 (sbit ansbit ee))
	       (setf (sbit ansbit ee) 1)
	       (incf termcount))) )
	   
	 ;; now we know exactly how many terms in answer, allocate storage.
	 ;; we could do this:  if the number of terms is nearly the same as hideg,
	 ;; just use mul-dense.
	 (setf ae (make-array termcount))
	 (setf ac (make-array termcount :initial-element 0))
	 ;; initialize exponents. j=0 to start
	 (dotimes (i (1+ hideg) ae)
	   (declare(fixnum i))
	   (unless (zerop (aref ansbit i))
	     (setf (aref ae j) i)
	     (incf j)))
      
	 (dotimes (i xl (cons ae ac))
	   (declare (fixnum i))
	   (setf xei (aref xe i) xci (aref xc i))
	   (setf guess -1)
	   (dotimes(j yl)
	     ;; multiply terms for real, this time.
	     (declare (fixnum j))

	     (linsearch-update ae ac  (the fixnum(+ xei(the fixnum (aref ye j))))
			       (* (aref yc j) xci) (1+ guess)))))))))



#+ignore

(defun linsearch-update (indexA valA  ind val guess) ; guess could be (ash (length indexA) -1)
  ;;  linear search and update. All exponents are in array indexA, all values in array valA
  ;; find j such that  ind=indexA[j]. Then add val to element valA[j]   
  ;; by construction there is such a j.  Also indexA is sorted.
  (declare (fixnum ind guess)
	   (type (simple-array t (*)) indexA valA))
  (let ((p guess)
	(thisind 0))
    (declare (fixnum p thisind))
    (loop
      (setf thisind (aref indexA p))
      (cond ((<=  ind thisind)
	     (cond ((eq thisind ind)
		    (incf (aref valA p) val)
		    (return p))
		   (t(decf p))))
	    (t (incf p))))))

;; recursive binsearch

;; hm. doesn't work

(defun bsur(indexA valA  ind val)	;binary search and update, recursive
  (declare (fixnum ind)
	   (type (simple-array t (*)) indexA valA))
  (let ((p 0)(midpt 0))
    (declare (fixnum p midpt))
  (labels 
      ((bsu (lo hi)
	 ;;  binary search and update. All exponents are in array indexA, all values in array valA
	 ;; find j such that  ind=indexA[j]. Then add val to element valA[j]   
	 ;; by construction there is such a j.
	    (declare (fixnum lo hi ind p midpt))
	      (setf p (+ lo (ash (- hi lo) -1))
		    midpt (aref indexA p)) ; midpoint, about.
	      (cond ((or (eq midpt ind)(eq p 0))(incf (aref valA p) val));; the only way out
		    ((< ind midpt)(bsu lo p))
		    (t(bsu p hi)))))
    (bsu 0 (1- (length indexA))))))







#|
(setf q (make-racg 5000 50 5) r (make-racg 5000 50 5))

cl-user(67):  (time (mul-heap5 q r))
; cpu time (non-gc) 26,984 msec user, 0 msec system  5 megabytes

cl-user(65):   (time (mul-naive-b q r))
; cpu time (non-gc) 9,172 msec user, 0 msec system

cl-user(64):  (time (mul-naive-bb q r))
; cpu time (non-gc) 5,922 msec user, 0 msec system   2 megabytes

cl-user(68): (time (mul-dense q r))
; cpu time (non-gc) 1,031 msec user, 0 msec system   3 megabytes

 cl-user(70): (time (mul-hash1 q r))
; cpu time (non-gc) 8,656 msec user, 16 msec system
[3] cl-user(76):  (time (mul-fft q r))
; cpu time (non-gc) 1,343 msec user, 0 msec system  10X as much storage, 34 megabytes


|#