;; 4/4/08 RJF
;; I can't find the hashtable polynomial multiplication, so I'll rewrite it.

#| assume a polynomial in 1 variable (implicit), is represented by a hashtable whose keys are fixnums and coefficients are integers.  Thus 5*x^15+4*x^90 is built by

(setf r (make-hash-table :test 'eql))
(setf (gethash 5 r) 30)			;30*x^5
(setf (gethash 9 r) 21)			; 21*x^9 

;; with some minor hacks, coefficients could, themselves, be hash
;; tables recursively.  some issues that usually need to be resolved,
;; especially in the multivariate case, is the consistent
;; representation of constants, and especially 0, as well as the
;; ordering of variables and storage of their names.

|#

(defun pdump(r)(maphash #'(lambda(key val)(format t "~%~a*x^~a" val key)) r))

(defun ptimes(r s)
  (declare  (hash-table r s) (optimize (speed 3)(safety 0)))
  (let ((ans  (make-hash-table :test 'eq :size 
					(* 10 (max(hash-table-count r) (hash-table-count s)))
			       )))
    (maphash #'(lambda(ex co)		; exponent, coefficient
		 (declare (fixnum ex) (integer co))
		 (maphash #'(lambda (ex2 co2)
			      (declare (fixnum ex2)(integer co2))
			      (incf (gethash (+ ex ex2) ans 0) (* co co2)))r))
	     s)
    ans))

(defun ppower(r n)
  (cond((= n 0) 1)
       ((= n 1) r)
       (t (ptimes r (ppower r (1- n))))))

(defun list2hash(k)  ;; k looks like ((30 . 5)(21 . 9)). Make a hashtable.
  (let ((ans (make-hash-table :test '=)))
    (dolist (i k ans)(setf (gethash (cdr i) ans) (car i)))))

(defun pplus-into (r s)  ;; r+s --> r
    (maphash #'(lambda(ex co) ; exponent, coefficient
			      (incf (gethash ex r 0) co))
     s)   r)


;; an encoding of (1+x+y+z)  with %,y=x^41,z=x^(41^2);

(setf w (list2hash '(( 1 . 1681) (1 . 41) (1 . 1) (1 . 0)))) 

;; try (setf f (ppower w 20)) (time (ptimes f f)) ;; benchmark test
;; discussed by rjf and sci.math.symblic recently.


;; ringing in some changes

(defun pdumps(r)			; like pdump but sorted; 0 coeffs removed
  (let ((ans nil))
    (maphash #'(lambda(key val)(unless (= 0 val)(push (cons key val) ans))) r)
    (map nil #'(lambda(k)(format t "~%~a*x^~a" (cdr k)(car k)))
	 (sort ans #'> :key #'car))
    (values)))

;; a different approach to the program

;; accumulate the answer by using an array in the middle.
;; actually faster than hashing ptimes above in our example of f*f;
;; but potentially a loser in the case of even sparser data.

(setf w2q  '(( 1 . 1681) (1 . 41) (1 . 1) (1 . 0))) ;; sample polynomial

(defun qtimes(r s)
  (declare   (optimize (speed 3)(safety 0)))
  (let ((ans (make-array (+ (cdar r)(cdar s) 1) :initial-element 0)))
    (map nil #'(lambda(pp)		; (exponent .  coefficient)
		 (map nil #'(lambda (qq)
			      (incf (aref ans(+ (cdr pp)(cdr qq))) (* (car pp)(car qq)))
					   ) r)) s)
    (tolist ans)))

(defun tolist(a)
  (let ((ans nil) (term 0))
    (dotimes (i (length a) ans)
     (setf term (aref a i)) 
      (unless (= 0 term)(push (cons  term i) ans)))))

(defun qpower(r n)
  (cond((= n 0) '( (0 . 1)))
       ((= n 1) r)
       (t (qtimes r (qpower r (1- n))))))

;;; yet another approach.  here we use exclusively lists.  we compute
;;; each of the products of monomial * polynomial2 and merge the
;;; results into a master list.  loses in 3 ways: 
;;; 1. excess consing,
;;; 2. merging is (eventually) not so fast as the master list grows,
;;; 3. bignum arithmetic in lisp is not as fast as gmp and does
;;; number-consing.  This last problem is same for anything in regular
;;; lisp

(defun rtimes(r s)
  (declare   (optimize (speed 3)(safety 0)))
  (let ((ans nil))
    (map nil #'(lambda(pp)		;; coefficient . exponent
		 (let ((row nil))
		   (map nil #'(lambda (qq)
				(push (cons (* (car pp)(car qq))
					    (+ (cdr pp)(cdr qq)))
				      row))r)
		 (setf ans (mergerow ans row)))) s)
    (nreverse ans)))

(defun mergerow (a b)(declare (optimize (speed 3)(safety 0)))
  (cond ((null a) b)((null b)a)
	((= (cdar a)(cdar b))
	 (rplaca (car a)  (+ (caar a)(caar b)))
	 (cons (car a)
	       (mergerow (cdr a)(cdr b))))
	((< (cdar a)(cdar b))
	 (cons (car a)(mergerow (cdr a) b)))
	(t (cons (car b)(mergerow a (cdr b))))))


;;; stimes  Another chance. 
;;; This time make the short row an array or two.


(defun stimes(r s) ;; assume r is no shorter than s, or make it happen
  (declare   (optimize (speed 3)(safety 0)))
  (let* ((ans nil)
	 (rlength (length r))
	 (row (make-array rlength)) ; one array, reused
	 (spot nil)
	 (qq r))
    (dotimes (i rlength)(setf (aref row i) (cons 0 0)))
    (map nil #'(lambda(pp)		;; coefficient . exponent
		 (dotimes (i rlength)
		   (setf spot (aref row i))
		   ;;(format t "~% cp=~s cq=~s"  pp(car qq))
		   (rplaca spot  (* (car pp)(caar qq)))
		   (rplacd spot  (+ (cdr pp)(cdar qq))))
		 (setf ans  (mergerowarray ans row 0 rlength))
		 (pop qq))
	 s)
	 (nreverse ans)))

;; this program is buggy.

(defun mergerowarray (a v j m)		; destructive of input a
  (declare (fixnum j m)
	   (optimize (speed 3)(safety 0)))
  (cond ((= j m) a)			; nothing more to do
	(t (let ((z (aref v j)))	; top item
	     
	     (cond((null a)
		   (push z a) 
		   (mergerowarray a v j m)
		   a)
		  ((= (cdar a)(cdr z))
		   (rplaca (car a) (+ (caar a)(car z)))
		   (mergerowarray (cdr a)v (1+ j) m)
		   a)
		  ((< (cdar a)(cdr z))
		   (mergerowarray (cdr a) v j m)
		   a)
		  (t (cons (cons (car z)(cdr z))  (mergerowarray a v (1+ j) m))))))))



;;; yet another version
;;; broken  4/5/08
(defun ttimes(r s) ;; assume r is no shorter than s, or make it happen
  (declare   (optimize (speed 3)(safety 0)))
  (let* ((ans nil);;  empty tree
	 (rlength (length r))
	 (row (make-array rlength)) ; one array, reused
	 (spot nil)
	 (qq r))
    (dotimes (i rlength)(setf (aref row i) (cons 0 0)))
    (map nil #'(lambda(pp)		;; coefficient . exponent
		 (dotimes (i rlength)
		   (setf spot (aref row i))
		   ;;(format t "~% cp=~s cq=~s"  pp(car qq))
		   (rplaca spot  (* (car pp)(caar qq)))
		   (rplacd spot  (+ (cdr pp)(cdar qq))))
		 (setf ans  (mrt ans row 0 rlength))
		 (pop qq))
	 s)
    (nreverse ans)))

;; tree looks like
;; nil empty
;; (( coe . exp) ptr ptr) ;; internal node
;; (coe .exp) leaf node.

(defun mrt (tt row j rlength) ;merge row with tree
  
  (if (= j rlength) tt
    (let ((spot (aref row j)))
      (cond ((null tt) 
	     (setf tt (list (cons (car spot)(cdr spot)) nil nil)) ;copy spot
	     (mrt tt row (1+ j) rlength)
		       tt) ;; wasteful?
	    ((= (cdar tt) (cdr spot))	; equal exponents
	     (rplaca (car tt)(+ (car spot)(caar tt)))
	     (mrt tt row (1+ j) rlength)
	     tt
	     )
	    ((> (cdar tt)(cdr spot))
	     (mrt(cadr tt) row j rlength) ;to left
	     tt)
	    (t
	     (mrt (caddr tt) row j rlength) ;to right
	    tt)))))


;;; variant on ptimes
;;; here the values in the hash table at location exp  are (coef . exp)
;;; this means that once you get the value out, you can smash the coef.

;;; um, works 4/4/08

(defun ptimesx(r s)
  (declare  (hash-table r s) (optimize (speed 3)(safety 0)))
  (let ((ans  (make-hash-table 
	       :test 'eq )))
    (maphash #'(lambda(ex val)		; index is exponent, val is (coef. exponent)
		 (declare (fixnum ex))
		 (maphash #'(lambda (ex2 val2)
			      (declare (fixnum ex2))
			      (let* ((expon (+ ex ex2))
				     (spot (gethash expon ans)))
				(cond (spot 
				       (incf (car spot) 
					     (* (car val)(car val2))))
				      (t (setf (gethash expon ans)
					   (cons  (* (car val)(car val2)) expon))))))r))
	     s)
    ans))

(defun list2hashx(k)  ;; k looks like ((30 . 5)(21 . 9)). Make a hashtable.
  (let ((ans (make-hash-table :test '=)))
    (dolist (i k ans)(setf (gethash (cdr i) ans) i))))

(defun pdumpsx(r)			; like pdump but sorted; 0 coeffs removed
  (let ((ans nil))
    (maphash #'(lambda(key val)(declare (ignore key))
		      (unless (= 0 (car val))(push  val ans))) r)
    (map nil #'(lambda(k)(format t "~%~a*x^~a" (car k)(cdr k)))
	 (sort ans #'> :key #'cdr))
    (values)))

(defun ppowerx(r n)
  (cond((= n 0) 1)
       ((= n 1) r)
       (t (ptimesx r (ppowerx r (1- n))))))


;;; another version, using gmp numbers
;;; works.
;;; not especially fast though.  2.8 seconds on benchmark,

(defun ptimesg(r s)
  (declare  (hash-table r s) (optimize (speed 3)(safety 0)))
  (let ((ans  (make-hash-table 
	       :test 'eq )))
    (maphash #'(lambda(ex val)		; index is exponent, val is (coef. exponent)
		 (declare (fixnum ex))
		 (maphash #'(lambda (ex2 val2)
			      (declare (fixnum ex2))
			      (let* ((expon (+ ex ex2))
				     (spot (gethash expon ans)))
				(declare (fixnum expon))
				(cond (spot 
				       (mpz_addmul (car spot) (car val)(car val2)))
				      (t (setf spot
					   (cons  
					    (create_mpz_zero) expon))
					 (mpz_mul (car spot)(car val)(car val2))			 
					 (setf (gethash expon ans) spot)
					   ))))r))
	     s)
    ans))

(defun list2hashg(k)  ;; k looks like ((30 . 5)(21 . 9)). Make a hashtable.
  (let ((ans (make-hash-table :test '=)))
    (dolist (i k ans)(setf (gethash (cdr i) ans) (cons (lisp2gmp (car i))(cdr i))))))

(defun pdumpsg(r)			; like pdump but sorted; 0 coeffs removed
  (let ((ans nil))
    (maphash #'(lambda(key val)(declare (ignore key))
		      ;;(unless (= 0 (car val))(push  val ans))
		      (push  val ans)
		      ) r)
    (map nil #'(lambda(k)(format t "~%~a*x^~a" (cm(car k))(cdr k)))
	 (sort ans #'> :key #'cdr))
    (values)))

(defun ppowerg(r n)
  (cond((= n 0) 1)
       ((= n 1) r)
       (t (ptimesg r (ppowerg r (1- n))))))


;; shorter. slower? yes. creates too much stuff.. 6.3 sec

#+ignore
(defun ptimesg(r s)
  (declare  (hash-table r s) (optimize (speed 3)(safety 0)))
  (let ((ans  (make-hash-table 
	       :test 'eq )))
    (maphash #'(lambda(ex val)		; index is exponent, val is (coef. exponent)
		 (declare (fixnum ex))
		 (maphash #'(lambda (ex2 val2)
			      (declare (fixnum ex2))
			      (let* ((expon (+ ex ex2))
				     (spot (gethash expon ans (cons (create_mpz_zero) expon))))
				(mpz_addmul (car spot) (car val)(car val2))
				))r))
	     s)
    ans))


;; 2.7 sec
(defun ptimesg(r s)
  (declare  (hash-table r s) (optimize (speed 3)(safety 0)))
  (let ((ans  (make-hash-table 
	       :test 'eq ))) ;faster than eql and much faster than =
    (maphash #'(lambda(ex val)		; index is exponent, val is (coef. exponent)
		 (declare (fixnum ex))
		 (maphash #'(lambda (ex2 val2)
			      (declare (fixnum ex2))
			      (let* ((expon (+ ex ex2))
				     (spot (gethash expon ans)))
				(cond ((null spot)
				       (setf (gethash expon ans)
					 (setf spot (cons (create_mpz_zero) expon)))))
				(mpz_addmul (car spot) (car val)(car val2))
				))r))
	     s)
    ans))


;; memq hacking

(defun m(a  b)(declare(optimize (speed 3)(safety 0)))
       (or (and b (eq a (car b)))
	   (and (cdr b)(eq a (cadr b)))
	   (and (cddr b)(eq a (caddr b)))
	   (member a (cdddr b) :test #'eq)))