;;; 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.
;; throw out all the experimental stuff and just leave in the working stuff.

(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
(defvar *karatlim* 1000) ;; for running.

;; 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)
	;;  (type (simple-array integer (*)) carp cdrp)
	     )
    (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))))))

(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 (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 '(simple-array fixnum (*)))(coerce coefs '(simple-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 
		  '(simple-array fixnum (*)) ) 
	  (coerce coefs '(simple-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 '(simple-array fixnum (*)) )
	  (coerce coefs '(simple-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)))

|#