;;; 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)))
|#