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