;;; -*- Lisp -*-
(eval-when (compile load) (declaim (optimize (speed 3)(safety 0))))
;;; Multiplication POLYALGORITHM for segmentation

;;;;;;;;;;;;; yet more hacks.  suppose the exponents are bignums.

#|strategy: decompose the input.  Let k = ( most-positive-fixnum-1) /2

P= p_0 + ... + p_n*x_n  into P_0 + x^k*P_1+ x^(k*r)*P_r

where P_0, ... P_r  have only fixnum exponents.
Same for Q.

Do the mults.

Put the pieces together. |#

(defstruct segpoly pieces shiftvals)

(defun segbyexp(p &optional(k #.(ash most-positive-fixnum -1)))
  ;; segment a polynomial by exponent range.
  ;;; find the next lowest exponent, say E
  ;;; encode all terms from E to E+k-1 as
  ;;; x^E* (poly-shifted by E)
  
  (let* ((n (aref (car p) (1- (length (car p)))))	; largest exponent
	 (pieces (make-array 2 :adjustable t :fill-pointer 0) )
	 (shiftvals (make-array 2 :adjustable t :fill-pointer 0) ) ; unknown number of pieces
	 )
    (cond ((< n k)			; no need to segment at all 
	   (vector-push-extend  p pieces)
	   (vector-push-extend  0 shiftvals)
	   (make-segpoly :pieces pieces :shiftvals shiftvals))
	  (t 
	   (let* ((lp (length (car p)))
		 (lim 0)(low 0) (expon 0)
		  (thepart nil)
		  (i 0))
	   ;; put all pieces with exponent less than i*n into ith cell.
	     ;; after reducing by k^i.

	     (setf i 0)	     
	     (loop while (< i lp) do
	     (setf low (aref (car p) i))
	     (setf lim (+ low k)) ; next piece has stuff up to lim
	     (setf thepart nil)
	     (loop while (and (< i lp)
			      (< (setf expon (aref (car p) i)) lim))
		 do
		   (push (cons (- expon low) (aref (cdr p) i)) thepart)
		   (incf i))
	     ;; check if there really were any parts in this interval... since the next
	     ;; exponent  is too big
	     (cond (thepart
		    (vector-push-extend (L2PA (nreverse thepart)) pieces)
		    (vector-push-extend low shiftvals))))
	     ;; at this point we have gone through the whole original polynomial
	     (make-segpoly :pieces pieces :shiftvals shiftvals))))))

;; this next program now works.

#+ignore
(defun segbygap(p g)
  ;; segment a polynomial by gap.
  ;;; any time there is a gap>g, make a new polynomial segment
  ;;; encode all terms from E to E+k-1 as
  ;;; x^E* (poly-shifted by E)
  (declare(optimize (speed 3)(safety 0)))
  (let* ((lp (length (car p)))
	 (lpm1 (1- lp))
	 (n (aref (car p) lpm1))	; largest exponent in polynomial p
	 )
    ;; no need to segment at all
    ;; though we might still benefit from shifting...
    (if (< n g)
	(return-from segbygap
	  (make-segpoly :pieces (vector p) :shiftvals(vector 0))))
	  
    (let* ((pieces (make-array 2 :adjustable t :fill-pointer 0) )
	   (shiftvals (make-array 2 :adjustable t :fill-pointer 0) ); unknown number of pieces
	   (expon 0)
	   (carp (car p));; the exponent array
	   (cdrp (cdr p));; the coefficient array
	   (low (aref carp 0))
	   (exps (make-array 2 :adjustable t :fill-pointer 0))
	   (coefs (make-array 2 :adjustable t :fill-pointer 0)) )
      (vector-push-extend 0 exps) ;; the first exponent is shifted to zero
      (vector-push-extend (aref cdrp 0) coefs) ;; the first coefficient

      (loop for i from 1 to lpm1 ;; for each exponent, we have to put it somewhere!
	   do
	   ;;(format t "~%low=~s, gap=~s" low   (- (setf expon (aref carp i))  (aref carp (1- i))))
	   (cond ((< (- (setf expon (aref carp i)) ;the gap is not too bid
			(aref carp (1- i)))g)
		  
		   (vector-push-extend  (- expon low) exps) ;extend this segment's exponents
		   (vector-push-extend  (aref cdrp i) coefs) ;extend coefs too
					;   (format t "~%extending =~s, ~s ~%i=~s" exps coefs i)
		   )
		  (t			; the gap is too big. 
		   ;; clear out the previous polynomial
		   (vector-push-extend (cons exps coefs) pieces)
		   (vector-push-extend low shiftvals)
		   ;; set up for the new one
		   (setf exps (make-array 2 :adjustable t :fill-pointer 0))
		   (setf coefs (make-array 2 :adjustable t :fill-pointer 0))
		   (vector-push-extend 0 exps) ;; the first exponent is shifted to zero
		   (vector-push-extend (aref cdrp i) coefs)
		   (setf low (aref carp i)))))
      
      ;; final exponent/coef
;      (vector-push-extend  (- (aref carp lpm1) low) exps) ;extend this segment's exponents
;      (vector-push-extend  (aref cdrp lpm1) coefs) ;extend coefs too
      (vector-push-extend (cons exps coefs) pieces)
      (vector-push-extend low shiftvals)

      (make-segpoly :pieces pieces :shiftvals shiftvals))))


(defun segbygap(p g)  
  ;; segment a polynomial by gap.
  ;;; any time there is a gap>g, make a new polynomial segment
  ;;; encode all terms from E to E+k-1 as
  ;;; x^E* (poly-shifted by E)
  (declare(optimize (speed 3)(safety 0)) )
  (let* ((lp (length (car p)))
	 (lpm1 (1- lp)) )
    (let* ((pieces (make-array 10 :adjustable t :fill-pointer 0) )
	   (shiftvals (make-array 10 :adjustable t :fill-pointer 0) ); unknown number of pieces
	   (expon 0)
	   (carp (car p));; the exponent array
	   (cdrp (cdr p));; the coefficient array
	   (low (aref carp 0))
	   (exps (make-array 10 :adjustable t :fill-pointer 0))
	   (coefs (make-array 10 :adjustable t :fill-pointer 0)) )
      (declare (fixnum lp lpm1) (type (array t (*) ) carp cdrp exps coefs) )
      (vector-push-extend 0 exps) ;; the first exponent is shifted to zero
      (vector-push-extend (aref cdrp 0) coefs) ;; the first coefficient

      (loop for i fixnum from 1 to lpm1 ;; for each exponent, we have to put it somewhere!
	   do
	   ;;(format t "~%low=~s, gap=~s" low   (- (setf expon (aref carp i))  (aref carp (1- i))))
	   (cond ((< (- (setf expon (aref carp i)) ;the gap is not too bid
			(aref carp (1- i)))g)
		  
		   (vector-push-extend  (- expon low) exps) ;extend this segment's exponents
		   (vector-push-extend  (aref cdrp i) coefs) ;extend coefs too
					;   (format t "~%extending =~s, ~s ~%i=~s" exps coefs i)
		   )
		  (t			; the gap is too big. 
		   ;; clear out the previous polynomial and its shift
		   (vector-push-extend (cons exps coefs) pieces)
		   (vector-push-extend low shiftvals)
		   ;; set up for the new segment by allocating new arrays
		   (setf exps (make-array 2 :adjustable t :fill-pointer 0))
		   (setf coefs (make-array 2 :adjustable t :fill-pointer 0))
		   (vector-push-extend 0 exps) ;; the first exponent is shifted to zero
		   (vector-push-extend (aref cdrp i) coefs)
		   (setf low (aref carp i)))))
      ;; clean up with last exponent/coef after ending the loop      
      (vector-push-extend (cons exps coefs) pieces)
      (vector-push-extend low shiftvals)
      
      (make-segpoly :pieces pieces :shiftvals shiftvals))))

(defun segbygap-count(p g)  
  ;; segment a polynomial by gap.  Just count them up!
  ;;; any time there is a gap>g, make a new polynomial segment
  ;;; encode all terms from E to E+k-1 as
  ;;; x^E* (poly-shifted by E)
  (declare(optimize (speed 3)(safety 0)) )
  (let* ((lp (length (car p)))
	 (lpm1 (1- lp)) )
    (let* ((pieces (make-array 10 :adjustable t :fill-pointer 0) )
	   (shiftvals (make-array 10 :adjustable t :fill-pointer 0) ); unknown number of pieces
	   (carp (car p));; the input exponent array
	   (low (aref carp 0)) ;; the smallest exponenet
	   (expcount 1) ;; initialize with first exponent
	   )
      (declare (fixnum lp lpm1 expcount) (type (array t (*) ) carp cdrp exps coefs) )
      (loop for i fixnum from 1 to lpm1 ;; for each exponent, we have to put it somewhere!
	   do
	   ;;(format t "~%low=~s, gap=~s" low   (- (setf expon (aref carp i))  (aref carp (1- i))))
	   (cond ((< (- (aref carp i) ;the gap is not too bid
			(aref carp (1- i)))g)
		  (incf expcount)
		  )
		  (t			; the gap is too big. 
		   ;; clear out the previous polynomial and its shift
		   (vector-push-extend expcount pieces) ; how many exponents in this segment
		   (vector-push-extend low shiftvals)   ; the shift value
		   ;; set up for the new segment
		   (setf expcount 1)
		   (setf low (aref carp i))
		   )))
      ;; clean up with last exponent/coef after ending the loop      
      (vector-push-extend (incf expcount) pieces)
      (vector-push-extend low shiftvals)
      
      (make-segpoly :pieces pieces :shiftvals shiftvals))))


;; This program returns a list of N triples, where N is the number of
;; segments in the polynomial p, when decomposed by gap size g.  Each
;; triple consists of E the number of elements in the segment, S the
;; shift value, and D the degree in that segment.

;; the polynomial segment is t^S*(some poly of degree d with E terms).
;; the polynomial will be of degree 0 if it is a single term.
;; 
(defun segbygap-countL(p g)   ;; use lists for storing data
  ;; segment a polynomial by gap.  Just count them up!
  ;;; any time there is a gap>g, make a new polynomial segment
  ;;; encode all terms from E to E+k-1 as
  ;;; x^E* (poly-shifted by E)
  (declare(optimize (speed 3)(safety 0)) )
  (let* ((lp (length (car p)))
	 (lpm1 (1- lp)) )
    (let* ((pieces nil )
	   (shiftvals nil ); unknown number of pieces
	   (carp (car p));; the input exponent array
	   (low (aref carp 0)) ;; the smallest exponenet
	   (expcount 1) ;; initialize with first exponent
	   (maxdeg 0)
	   )
      (declare (fixnum lp lpm1 expcount) (type (array t (*) ) carp cdrp exps coefs) )
      (loop for i fixnum from 1 to lpm1 ;; for each exponent, we have to put it somewhere!
	   do
	   (cond ((< (- (aref carp i) ;the gap is not too bid
			(aref carp (1- i)))g)
		  (incf expcount) )
		  (t			; the gap is too big. 
		   ;; clear out the previous polynomial and its shift
		   (push expcount pieces) ; how many exponents in this segment
		   (push (- (aref carp (1- i)) low) maxdeg) ;; max degree in this segment
		   (push low shiftvals)   ; the shift value
		   ;; set up for the new segment
		   (setf expcount 1)
		   (setf low (aref carp i))
		   )))
      ;; clean up with last exponent/coef after ending the loop      
      (push expcount pieces)
      (push low shiftvals)
      (push (- (aref carp lpm1)low)  maxdeg)
      ;; the length of each segment
      ;; the shift multiplier for each segment
      ;; the max degree for each segment
      (mapcar #'list  pieces shiftvals maxdeg))))

		     

(defun unseg (a) ;; a is a segmented polynomial.  return it as a single polynomial.
  (reduce #'add-poly
	  (let ((p (segpoly-pieces a))
		(s (segpoly-shiftvals a)))
	    (loop for i from 0 to (1- (length s)) collect
		  (cons  (map 'array #'(lambda(q)(+ q (aref s i))) (car (aref p i)))
			 (cdr (aref p i)))))))
  






;;
#|To  multiply segmented polynomials, one could consider a karatsuba-like scheme, or just NXM.
This is just NXM.
|#


(defun mul-seg(p q  &key (mulfun #'mul-naive) )

  (let ((result (cons (vector)(vector))) ;or whatever is needed..
	(pseg (segpoly-pieces p))
	(psv(segpoly-shiftvals p))
	(qseg (segpoly-pieces q))
	(qsv (segpoly-shiftvals q))
	(psvi 0)
)
    
    
    ;; if we had L  processors, we could do L of these NxM multiplies
    ;; at the same time and do an L-1 way merge. Here we just use linear merge, bad.
    (dotimes (i (length pseg) result)
      (setf psvi (aref psv i))
      (dotimes(j (length qseg))
	(setf result (add-poly result ; not necessarily fixnum exponents
			       (shiftexp
				(+ psvi (aref qsv j))
				(funcall mulfun  ;; fixnum exponents
					 (aref pseg i)(aref qseg j)))))))))

(defun shiftexp(n p) ;shift the exponents of polynomial p by multiplying by x^n, ie add n to exp
  (if (= n 1) p 
    (cons  ;; the "functional" way
	#+ignore  (map 'array #'(lambda(r)(+ n r))	(the (simple-array t (*))(car p)))
	  ;; the modify-in-place way
	  (loop for i from 0 to (1- (length r)) finally (return r) do (incf (aref r i)n))
	  
	(cdr p))))

(defun add-poly (p q) ;;works
  (let* ((anse (make-array 10 :adjustable t :fill-pointer 0)); used by seg mult
	 (ansc (make-array 10 :adjustable t :fill-pointer 0))
	 (i 0)(j 0)
	 (ep (car p))			;exponents
	 (eq (car q))
	 (cp (cdr p))			;coefs
	 (cq (cdr q))
	 (epi 0)(eqj 0)
	 (lp (1-(length ep)))
	 (lq (1-(length eq))))
    (declare (fixnum lp lq i j)
	     (type (simple-array t (*)) ep eq cp cq)
	     (type (array t (*)) anse ansc)
	     (integer epi eqj))
    (loop
   ;;   (format t "~%i=~s j=~s anse=~s" i j anse)
     (cond ((> i lp)
	    (do ((jj j (1+ jj)))
		((> jj lq) (return-from add-poly (cons 
						   (make-simple-array anse)
						   (make-simple-array ansc))))
	      (declare (fixnum jj))
	      (vector-push-extend  (aref eq jj) anse)
	      (vector-push-extend  (aref cq jj) ansc)))
	  
	   ((> j lq)
	    (do ((jj i (1+ jj)))
		((> jj lp) (return-from add-poly  (cons 
						   (make-simple-array anse)
						   (make-simple-array ansc))))
	      (declare (fixnum jj))
	      (vector-push-extend  (aref ep jj) anse)
	      (vector-push-extend  (aref cp jj) ansc)))

	   ((< (setf epi (aref ep i))(setf eqj (aref eq j)))
	    (vector-push-extend  epi anse)
	    (vector-push-extend  (aref cp i) ansc)
	    (incf i))
	   
	   ((> epi eqj)
	    (vector-push-extend  eqj anse)
	    (vector-push-extend  (aref cq j) ansc)
	    (incf j))
	   
	   (t ;; (= epi epj)
	    (vector-push-extend eqj anse) ;exponents are equal
	    (vector-push-extend (+ (aref cp i)(aref cq j)) ansc)
	    (incf i)(incf j))) )))

	   
	   
(defun test-segexp(u v n)
  (mul-seg (segbyexp u n)(segbyexp v n) n))

(defun test-seggap(u v n)
  (mul-seg (segbygap u n)(segbygap v n) n))
 
    
#|    
	
;; is fixnum faster than signed-byte 32?
(defun 
    f1(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 32) (*)) z))
  (incf (aref z 0)(aref z 1))  ;; 146 bytes!!
  
  z)

(defun 
    f12(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 32) (*)) z))
  (incf (the (signed-byte 32)(aref z 0))(aref z 1))  ;; only  15 bytes.
  z)

(defun 
    f2(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array fixnum (*)) z))
  (incf (aref z 0)(aref z 1))  ;; 15 bytes also
  z)

	

(defun 
    f3(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 24) (*)) z))
  (incf (aref z 0)(aref z 1))
  z)

(defun 
    f4(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 31) (*)) z)) ;; also 15 bytes
  (incf (aref z 0)(aref z 1)) ;; adding is same; multiplying of fixnum needs shift of 2  bits
  z)

(defun 
    f11(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 32) (*)) z))
  (setf (aref z 0)(the (signed-byte 32) (* (aref z 0)(aref z 1)))) ;; 35 bytes
  z)


(defun 
    g4(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 31) (*)) z)) ;; winner, same as fixnum, though.
  (setf (aref z 0)(*(aref z 0)(aref z 1)))
  z)

(defun 
    g5(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array fixnum (*)) z)) ;; winner, same as fixnum, though.
  (setf (aref z 0)(*(aref z 0)(aref z 1)))
  z)

(defun 
    g6(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array fixnum (*)) z)) ;; winner, same as fixnum, though.
  (setf (aref z 0)(the fixnum (*(aref z 0)(aref z 1))))
  z)

(defun 
    g7(z)
  (declare (optimize (speed 3)(safety 0))
	   (type (simple-array (signed-byte 31) (*)) z)) ;; winner, same as fixnum, though.
  (setf (aref z 0)(the (signed-byte 31) (*(aref z 0)(aref z 1))))
  z)
|#

;; test case for (1+x+y+z+t)^20 , squared... or something like it
#|
(setf f (cons #(1 40  #.(* 40 40) #.(* 40 40 40) #.(* 40 40 40 40))
	      #(1 1 1 1 1)))

(defun sq(x)(mul-heap x x))  ;; in mulheap5.lisp for example
(setf f4 (sq (sq f)))
(setf f16 (sq f4))
(setf f20 (mul-heap f4 f16))
|#

(defun segbygap-L(p g)   ;; use lists for storing data
  ;; segment a polynomial by gap.  Count AND create segments
  ;;; any time there is a gap>g, make a new polynomial segment
  ;;; encode all terms from E to E+k-1 as
  ;;; x^E* (poly-shifted by E)
  (declare(optimize (speed 3)(safety 0)) )
  (let* ((lp (length (car p)))
	 (lpm1 (1- lp)) )
    (let* ((pieces nil )
	   (shiftvals nil ); unknown number of pieces
	   (carp (car p));; the input exponent array
	   (cdrp (cdr p));; the input exponent array
	   (low (aref carp 0)) ;; the smallest exponenet
	   (expcount 1) ;; initialize with first exponent
	   (low-i 0)
	   (exps nil)
	   )
      (declare (fixnum lp lpm1 expcount) (type (array t (*) ) carp cdrp exps coefs) )
      (loop for i fixnum from 1 to lpm1 ;; for each exponent, we have to put it somewhere!
	   do
	   ;;(format t "~%low=~s, gap=~s" low   (- (setf expon (aref carp i))  (aref carp (1- i))))
	   (cond ((< (- (aref carp i) ;the gap is not too bid
			(aref carp (1- i)))g)
		  (incf expcount))
		  (t			; the gap is too big. 
		   ;; clear out the previous polynomial and its shift
		  ; (push expcount pieces) ; how many exponents in this segment
		   (setf exps (make-array expcount :element-type 'fixnum))
		   (loop for i from 0 to (1- expcount) do
			 (setf (aref exps i)(- (aref carp (+ i low-i)) low)))
		   (push (cons exps  (subseq cdrp low-i (+ low-i expcount)))
			 pieces)
		   
		   (push low shiftvals)   ; the shift value
		   ;; set up for the new segment
		   (setf expcount 1)
		   (setf low (aref carp i))
		   (setf low-i i)
		   )))
      ;; clean up with last exponent/coef after ending the loop      
         (setf exps (make-array expcount :element-type 'fixnum))
	 (loop for i from 0 to (1- expcount) do
	       (setf (aref exps i)(- (aref carp (+ i low-i)) low)))
	 (push (cons exps  (subseq cdrp low-i (+ low-i expcount)))
			 pieces)
	 (push low shiftvals)
      ;; reverse to make it look same for debugging.  reversal not necessary really
	(make-segpoly :pieces (nreverse pieces) :shiftvals (nreverse shiftvals))
      ;(make-segpoly :pieces pieces :shiftvals shiftvals)
      )))