;;; -*- Lisp -*- ;;; various ways of multiplying sparse polynomials. ;;; using lists, hashtables, trees, streams, heaps, and just hacking the naive method. ;;; ;;; TRYING GMP numbers. ;;; :ld gmp.dll ;;; :ld simplegmp.fasl (eval-when (compile load) (declaim (optimize (speed 3)(safety 0)))) ;; set up a lexically available context for handling several ;; variables which would otherwise be global, and require more ;; effort to access from compiled code. ;; (let ((*YYEXPS* nil) (*YYCOEFS* nil) (*YYLEN* 0) (*zcoefs-lim* 50) (*zcoefs* (make-array 50)) (fill-pointer 0)) (declare (fixnum *YYLEN* fill-pointer) (type (simple-array t (*)) *YYEXPS* *YYCOEFS*)) ;; input format: a univariate polynomial is represented by a Lisp cons ;; of 2 arrays. The first array is exponents, in ascending order. The ;; second array is corresponding coefficients. We expect that the ;; usual way of representing numbers in Lisp is used, whereby a small ;; number (in magnitude less than, say 2^29) is represented "in place" ;; shifted right a few bits so it is tagged as a small number; larger ;; numbers, of unlimited length, are pointed to by the entry in the ;; same field. Thus "short" arithmetic in Lisp is slightly ;; handicapped by the need to do some tag checking and shifting. ;;This makes the transition to arbitrary precision, as well as to ;; float, double-float, and even complex numbers programmatically ;; smoother. ;; If the Lisp compiler is informed (via declarations) that it has an array of small ;; numbers, it may dispense with the tags and generate code to use ordinary 32-bit (or ;; 64-bit) arithmetic. ;; example of the input format: 45x^10+25x^2+ 9 [note: 9 = 9x^0] ;; ( #(0 2 10) . (9 25 45)) ;; Now to multiplicatio of polynomials. ;; First, the really obvious way, using a double-nested loop, ;; inserting each of the new exponent-coefficient pairs in ;; the right place in a linked list. (defun mul-naiveG (x y);; naive but with a finger into the answer use GMP (let* ((ans (list (cons 0 0))) ; makes iteration neater (finger ans) (yl (length (the simple-array (car y)))) (xl (length (the simple-array (car x)))) (xe (car x)) (ye (car y)) (xc (cdr x)) (yc (cdr y)) ; (gtemp (gmp::create_mpz_zero)) ; (g-inside (gmp::gmpz-z gtemp)) ) (declare(fixnum yl xl) (type (simple-array t (*)) xe ye xc yc)) (dotimes (i xl (L2PA (cdr ans))) (dotimes(j yl) ;; multiply terms ; (gmp::mpz_mul g-inside (gmp::gmpz-z )(gmp::gmpz-z(aref yc j))) (setf finger (addtermG(+ (aref xe i)(aref ye j)) (aref xc i) (aref yc j) finger))) (setf finger ans)))) ;; We slightly improve the addterm program that returns an ;; appropriate finger into the list for insertion. ;; Instead of the recursive version we use an iterative format. (defun gmpmul(c1 c2) (let ((r (gmp::create_mpz_zero))) (gmp::mpz_mul (gmp::gmpz-z r)(gmp::gmpz-z c1)(gmp::gmpz-z c2)) r)) (defun addtermG(e c1 c2 L) ;iterative version with do ; 258 bytes (do ((ahead (caadr L)(caadr L))) ((null ahead) (setf (cdr L)(list (cons e (gmpmul c1 c2))))) (cond ((> e ahead) (setf L (cdr L))) ((= e ahead) (gmp::mpz_addmul(gmp::gmpz-z (cdadr L)) (gmp::gmpz-z c1)(gmp::gmpz-z c2)) ;; like (incf (cdadr L) (* c1 c2)) (return (cdr L))) (t (setf (cdr L)(cons (cons e (gmpmul c1 c2)) (cdr L))) (return (cdr L)))))) ;; OK, mul-naive works. How much faster is it than mul-naive? ;; On one example of a 500 X 1000 term multiplication, mul-naive takes 10.500 seconds. ;; The improved one takes 0.047 seconds, and is therefore 223 times faster. ;; Next step is to make a testing program. ;; First, we specify a program to create some random data. (defun make-racgG(n d c) ;; GMP version; change for extra arg c ;; make Random Array of increasing exponents of length n with ;; spacing k, with 0<k<=d between elements and with Coefficients that ;; are random integers between 1 and c. c may be, in General, arbitrary precision. (let ((exps (make-array n)) (coefs(make-array n)) (v (random d))) (dotimes (i n (cons exps coefs)) (setf (aref coefs i) (gmp::lisp2gmpz(+ 1 (random c)))) (setf (aref exps i) v) (setf v (+ 1 v (random d)))))) (defvar *x* nil);; these will be global polynomials for testing (defvar *y* nil) (defun test0(n m &optional (d 5) (c 5)) ;try these parameters for now. (setf *x* (make-racgG n d c)) (setf *y* (make-racgG m d c)) (format t "~% mul-naiveG time - just use insertion into single list, using GMP") (time (mul-naiveG *x* *y*)) ;; We will add more lines here later (format t "~% mul-hash1G time - hash1 using GMP") (time (mul-hash1G *x* *y*)) (format t "~% mul-strG time - heapmerge using GMP") (time (mul-strG *x* *y*)) ) (defun mul-hash1G (x y);; exponents may be arbitrary integers. coefficients numbers. (let ((ans (make-hash-table :test 'eql));; eql is good enough for bignums too (xexps (car x)) (xcoefs (cdr x)) (yexps (car y)) (ycoefs (cdr y)) (xei 0) (xco 0) (zcoefs *zcoefs*) (zlim *zcoefs-lim*) (zcounter 0) (loc nil) (xy 0) (gg nil)) (declare (type (simple-array t (*)) xexps yexps); any exponents (type (simple-array t (*)) xcoefs ycoefs zcoefs); any coefs (fixnum zcounter zlim) (integer xei xy) (optimize (speed 3)(safety 0))) (dotimes (i (length xexps)(sorthash1 ans zcoefs)) (declare (fixnum i)) (setf xei (aref xexps i)) ;exponent for x (setf xco (gmp::gmpz-z(aref xcoefs i))) ;coefficient corresponding to that exponent (dotimes (j (length yexps)) (declare (fixnum j)) ;; look in the hash table at location for this exponent (setf loc (gethash (setf xy (+ xei (aref yexps j))) ans)) (cond (loc ;; There is an entry in hash table for this ;; exponent. Update corresponding array location. Note: ;; hashtable is NOT updated. this happens about n*m ;; times, where n=size(x), m=size(y). ;;(incf (aref zcoefs loc);; no number consing here either?? ;; (* xco (aref ycoefs j))) ;; gmp question... we could just store INsides; save memory, screw up GC? (setf gg(aref zcoefs loc)) ; (gmp::mpz_addmul (gmp::gmpz-z gg) xco (gmp::gmpz-z (aref ycoefs j))) ) (t ;; if loc is nil, allocate a space in the answer array. ;; This happens relatively infrequently, about k times, where k=size(x*y). ;; Hashtable is updated. (setf (gethash xy ans) zcounter) ;; store the result in zcoefs array. ;;(setf (aref zcoefs zcounter) (* xco (aref ycoefs j))) ; no number consing?? (setf gg (gmp::create_mpz_zero)) (gmp::mpz_mul (gmp::gmpz-z gg) xco (gmp::gmpz-z (aref ycoefs j))) (setf (aref zcoefs zcounter) gg) (incf zcounter) (cond ((= zcounter zlim) ;;if we reach this, we need more working space. ;; This is rarely used code. If global array is ;; big enough, it is never used. If global array ;; is too small, it is made bigger. (let* ((newsize (ash zlim 1)); double the size ;; (newarray (make-array newsize)) (newarray (make-array newsize :initial-element 0))) (declare (fixnum newsize) (type (simple-array t(*)) newarray)) ;;(format t "~% sizes: newarray ~s zcoefs ~s newsize ~s" ;; (length newarray)(length zcoefs) newsize) (dotimes ;(m newsize) ;wrong (m zlim) (declare (fixnum m)) (setf (aref newarray m)(aref zcoefs m))) (setf zcoefs newarray zlim newsize *zcoefs-lim* zlim *zcoefs* zcoefs)))))))))) ;;; heap stuff here.. modified for GMP (defstruct (st) ;; A kind of stream. each structure contains 5 items ;; which by reference to arrays of XX and YY , can burp out ;; successive exponents and coefficients of the product. (exp 0 :type integer) (coef 0 :type integer) (i 0 :type fixnum) (v 0 :type integer) (e 0 :type integer)) (defun make-simple-array (A) (if (typep A 'simple-array) A (make-array (length A) :element-type (array-element-type A) :initial-contents A))) (defun mularZ4(x y) (declare (optimize (safety 3)(speed 0))) (let* ((xexps (car x)) (xcoefs (cdr x)) (yexps (make-simple-array (car y))) (ycoefs(make-simple-array (cdr y))) (yexp0 (aref yexps 0)) (ycoef0 (aref ycoefs 0)) ;should be gmpz (ans (make-array (1+(length xexps)))) (ylen (length (car y))) (xe 0) (v 0)) ; (declare (fixnum i) (type (simple-array t (*)) ans xexps xcoefs yexps ycoefs)) (setf *YYEXPS* yexps *YYCOEFS* ycoefs *YYLEN* ylen) (setf (aref ans 0) (make-st :exp 0 :coef (gmp::create_mpz_zero) :i ylen :v (gmp::create_mpz_zero) :e 0)) (dotimes (j (length xexps) ans) (setf v (gmp::lisp2gmpz(aref xcoefs j))) (setf xe (aref xexps j)) (setf (aref ans (1+ j)) (make-st :exp (+ xe yexp0) :coef (gmpmul v ycoef0) :i 1 :v v :e xe) )))) (defun update-strheap (h) ;;return top of heap h, which is necessarily not empty. ;; take successor in stream of top of heap, and insert it in heap. (declare(type (simple-array t(*)) h)(optimize (speed 3)(safety 0))) (let* ((top (aref h 1)) (ta (st-exp top))(tb (st-coef top)) (newe (updatestrfun top))) ;; in case newe <= lefti and newe <= righti. put it at the top. ;; avoiding all heap programs ;; top of heap is location 1 (ignore 0) ;; left branch is location 2 ;; right branch is location 3 (cond ;if end of stream from this subproduct, remove top and make heap smaller! ((null newe) (heap-remove h)) ;; we've modified the top node! leave it there! ((and(<= newe (st-exp (aref h 2))) (<= newe (st-exp (aref h 3))))) ;; it's not the top node. remove it and put the (revised) node back in. (t (heap-insert h (heap-remove h) fill-pointer))); the new item doesn't fit there (values ta tb))) (defun mul-strG(x y);; hacked to insert at top of heap if possible (declare (optimize(speed 3)(safety 0))) (let* ((hh (if (< (length (car x))(length (car y)))(mularZ4 x y)(mularZ4 y x))) (alen 10) ; heuristic(* (length (car x))(length (car y)))) ;; initial heap (anse (make-array alen :adjustable t :element-type 'integer :fill-pointer 0)) (ansc (make-array alen :adjustable t :element-type 'integer :fill-pointer 0)) (preve 0) (prevc (gmp::lisp2gmpz 0)) (prevci (gmp::gmpz-z prevc))) (declare (integer preve prevc) (type (array t(*)) anse ansc)) (setf fill-pointer (length hh)) (loop while (not (heap-empty-p)) do (multiple-value-bind (newe newc) (update-strheap hh) (declare (integer newe) ) ;; (format t "~% preve ~s prevc ~s prevci ~s newe ~s newc~s" preve prevc prevci newe newc) ; take next term off heap and insert its successor (cond ((= preve newe) ; is this same exponent as previous? ;; (incf prevc newc) ; if so, add coefficient. (gmp::mpz_add prevci prevci (gmp::gmpz-z newc)) preve) (t (unless (= 0 (aref prevci 1));; test for zero value ; to normalize: remove zero coeffs (vector-push-extend prevc ansc) (vector-push-extend preve anse)) (setf prevc newc) ;a gmp number (setf prevci (gmp::gmpz-z prevc)) (setf preve newe))); initialize new coefficient )) (vector-push-extend prevc ansc) ;put in last results before exit (vector-push-extend preve anse) (cons anse ansc))) (defun updatestrfun(rec) ; record is a stream 5-tuple (declare(type (simple-array t(*)) *YYCOEFS* *YYEXPS*) (optimize (speed 3)(safety 0)) (type st rec)) (let ((i (st-i rec))) (cond ((>= i *YYLEN*) nil);; no more items on stream of products (t (incf (st-i rec)) (gmp::mpz_mul(gmp::gmpz-z (st-coef rec)) (gmp::gmpz-z (st-v rec)) (gmp::gmpz-z (aref *YYCOEFS* i))) ;; (setf (st-coef rec) (* (the integer (st-v rec))(the integer (aref *YYCOEFS* i)))) (setf (st-exp rec) (+ (the integer (st-e rec))(the integer (aref *YYEXPS* i))) ; returns this val ))))) (defmacro heap-empty-p () "Returns non-NIL if & only if the heap contains no items." `(<= fill-pointer 1)) (defmacro hlessfun (a b) `(< (st-exp ,a) (st-exp ,b))) (defmacro lesser-child (a parent fp) "Return the index of the lesser child. If there's one child, return its index. If there are no children, return (FILL-POINTER (HEAP-A HEAP))." `(let* ((left;;(+ (the fixnum ,parent)(the fixnum ,parent)) (ash ,parent 1)) (right (1+ (the fixnum left)))) (declare (fixnum left right ,fp ,parent ) (optimize(speed 3)(safety 0)) (type (simple-array t (*)) ,a)) (cond ((>= left ,fp) ,fp) ((= right ,fp) left) ((hlessfun (aref ,a left) (aref ,a right)) left) (t right)))) (defun percolate-down (aa hole xx fp) "Private. Move the HOLE down until it's in a location suitable for X. Return the new index of the hole." (let ((xxval (st-exp xx))) (declare (fixnum hole fp)(integer xxval) (type (simple-array t (*)) aa)(optimize (speed 3)(safety 0))) (do* ((child (lesser-child aa hole fp) (lesser-child aa hole fp)) (aaval (aref aa child)(aref aa child))) ((or (>= (the fixnum child) fp) (< xxval (the fixnum(st-exp aaval)))) hole) (declare (fixnum child)) (setf (aref aa hole) aaval hole child)))) (defun percolate-up (a hole xx) "Private. Moves the HOLE until it's in a location suitable for holding X. Does not actually bind X to the HOLE. Returns the new index of the HOLE. The hole itself percolates down; it's the X that percolates up." (declare (type (simple-array t (*)) a)(fixnum hole) (optimize (speed 3)(safety 0))) (setf (aref a 0) xx) (do ((i hole parent) (parent (ash hole -1);;(floor (/ hole 2)) (ash parent -1);;(floor (/ parent 2)) )) ;; potential to speed up line below by declaration if a, x are fixnum, ((not (hlessfun xx (aref a parent))) i) (declare (fixnum i parent)) (setf (aref a i) (aref a parent)))) (defun heap-insert (a xx fp) "Insert a new element into the heap. Returns the heap.";; rjf (declare (type (simple-array t (*)) a)(fixnum fill-pointer) (optimize (speed 3)(safety 0))) ;; (format t ".") (setf (aref a fp) nil) (incf fill-pointer) ; the global one ;; Move the hole from the end towards the front of the ;; queue until it is in the right position for the new ;; element. (setf (aref a (percolate-up a fp xx)) xx)) (defun heap-remove(a) ;; assumes non-empty!! if empty next answer is bogus. ;; answer after that is an error (fill pointer can't pop) (declare(type (simple-array t(*)) a)(optimize (speed 3)(safety 0))) (let* ((xx (aref a 1)) (fp (decf fill-pointer)) ;faster, fp is local now (last-object (aref a fp))) (declare (fixnum fp)) (setf (aref a (the fixnum (percolate-down a 1 last-object fp))) last-object) xx)) ) (defun sorthash1(h ar) ;; utility for mul-hash1 (let* ((htc (hash-table-count h)) (anse (make-array htc :fill-pointer 0)) (ansc (make-array htc :fill-pointer 0))) (maphash #'(lambda (i v) (setf v (aref ar v)); get the coefficient, by indirection (vector-push i anse) (vector-push v ansc)) h) (qs anse ansc 0 (1- (length anse))) ;; qs is a quick sort we wrote, with a twist. (cons anse ansc) )) (defun qs(a b lo hi) ;; a re-implementation of quicksort, where entries in a are used as keys ; ;for sorting b too. (let ((pivotloc 0)) (declare(type (array t(*)) a b) (fixnum lo hi pivotloc)) (cond ((< lo hi) (setf pivotloc (partition a b lo hi)) (qs a b lo (1- pivotloc)) (qs a b (1+ pivotloc) hi) nil)))) (defmacro swapab(i j) `(let((temp& (aref a ,i))) (setf (aref a ,i) (aref a ,j)) (setf (aref a ,j) temp&) (setf temp& (aref b ,i)) ;also update b. this is why we wrote qs (setf (aref b ,i) (aref b ,j)) (setf (aref b ,j) temp&) )) (defun partition(a b lo hi) ;; the partition step in quicksort (declare(type (array t(*)) a b ) (fixnum lo hi)) (let ((left lo)(right hi)(pivot (aref a lo)) (bpivot (aref b lo))) ;pick a[lo] for pivot (declare(fixnum left right)) (loop while (< left right) do (while (< pivot (aref a right)) (decf right)) (while (and (<= (aref a left) pivot)(< left hi)) (incf left)) ;; at this point a[left]>pivot and a[right]<=pivot. ;; or left=hi (if (< left right) (swapab left right)) ) (setf (aref a lo) (aref a right)) (setf (aref b lo) (aref b right)) (setf (aref a right) pivot) (setf (aref b right) bpivot) ;also update b array. This is why we wrote qs right)) (defun setupG() (let ((k1 nil)(k2 nil) (gone (gmp::lisp2gmpz 1)) (gtwo (gmp::lisp2gmpz 2))) (setf *r* (cons #(0 1 100 10000 1000000) (vector gone gone gone gone gone))) ; 1+x+y+z+t (setf *s2* (power *r* 20)) (setf k1 (make-array 10626 :initial-contents (car *s2*))) (setf k2 (make-array 10626 :initial-contents (cdr *s2*))) (setf (aref k2 0) gtwo) (setf *s2p1* (cons k1 k2)))) (defun power(p n)(if (eql n 1) p(mul-naiveG p (power p (1- n)))))