;; new, 7.2010 RJF
;; multiply by Kronecker trick...
(load "gmp-gen.fasl")

(defun mul-kron(p q) ;; p, q are polynomials as #(expons) . #(coefs)
  (let* ((v (padpower2 p q))
	 (pr (tobig  (car p)(cdr p) v))
	 (qr (tobig  (car q)(cdr q) v))	 )
    (mpz_mul (gmpz-z qr) (gmpz-z qr)(gmpz-z pr))
    (frombig  qr v) ))

;; maximum abs value in an array
(defun arraymax (A)(let ((neg 0)(pos 0)(v 0))
		     (declare (optimize (speed 3)(safety 0)))
		     (loop for i fixnum from 0 below (length A) do
			   (setf v (aref A i))
			   (if (< v 0)(if (< v neg)(setf neg v))
			     (if (> v pos)(setf pos v))))
		     (max pos (- neg))))

(defun padpower2(p1 p2)
  (declare (optimize (speed 3)(safety 0)))
  (ceiling (+ 
	    (integer-length 
	       (min (aref (car p1) (1-(length (car p1))))
		    (aref (car p2) (1-(length (car p2))))  )    ) 
	      ;; the max exponents
	    (integer-length (arraymax (cdr p1)))
	    (integer-length (arraymax (cdr p2))) )
	   ))
#+ignore  ;; should be faster, but is broken  ...

(defun tobig (expons coefs logv)	;; change a polynomial into a bignum in GMP
  (declare (fixnum logv)(optimize (speed 3)(safety 0)))
  (let* ((g (create_mpz_zero))
	 (r (gmpz-z g)) ;the inside of object g
	 (j 0)
	 (c 0))
    
    (declare (fixnum j))
    (dotimes (i (length expons) g) ;return g when done
      (declare (fixnum i))
      (mpz_mul_2exp r r (* logv (-(aref expons i) j)));; res <- res*2^(logv*delta)
      (if (fixnump (setf c (aref coefs 0)))
	(if (> c 0) (mpz_add_ui r r c) (mpz_sub_ui r r (- c))) ;adding in fixnums
	(mpz_add r r (gmpz-z(lisp2gmpz c))))
 ;; converting bignums here
      (setf j (aref expons i)))))

;;; works but may be slower..

(defun tobig (expons coefs logv)	;; change a polynomial into a bignum in GMP
  (declare (fixnum logv)(optimize (speed 3)(safety 0)))
  (let* ((g (create_mpz_zero))
	 (r (gmpz-z g)) ;the inside of object g
	 (c 0)
	 (sr (create_mpz_zero))
	 (s (gmpz-z sr)))
    (declare (fixnum j))
    (dotimes (i (length expons) g) ;return g when done
      (declare (fixnum i))
      (cond ((fixnump (setf c (aref coefs i))) ;convert coef to gmp
	     (mpz_set_si s c))
	    (t (setf s (gmpz-z(lisp2gmpz c)))))
      (mpz_mul_2exp s s (* logv (aref expons i))) ;; scale the coefficient by the exponent
      (mpz_add r r s) ;; add it in.  ;; we wanted to add then shift... what did we do wrong?
    )))

;; works
(defun frombig (in logv) ;; change a bignum in GMP to a polynomial.
 (declare (fixnum logv)(optimize (speed 3)(safety 0)))
  (if (< logv 53)(frombig-small-logv in logv)
  (let* ((anse (make-array 100 :adjustable t :fill-pointer 0))
	 (ansc (make-array 100 :adjustable t :fill-pointer 0))
         (q (gmpz-z(create_mpz_zero)))
         (rz(create_mpz_zero))
	 (r (gmpz-z rz))
	 (i (gmpz-z in))
         (one (gmpz-z(create_mpz_zero)))        ; will be one
	(vby2 (gmpz-z(create_mpz_zero))) ; will be v/2
         (v (gmpz-z(create_mpz_zero)))	; will be v
	 (expon 0)
         (signum (>= (signum (aref i 1)) 0))) ; t if i non-neg
    (declare (fixnum expon))
    (mpz_set_si one 1)                  ; set to one
    (mpz_mul_2exp vby2 one (1- logv))   ; set to v/2
    (mpz_mul_2exp v one logv)           ; set to v
    (if signum nil (mpz_neg i i))       ; make i positive
    (while (> (aref i 1) 0)		; will be 0 when i is zero
      
      (mpz_fdiv_r_2exp r i logv)        ; set r to remainder by power of 2 (mask)
      (mpz_fdiv_q_2exp q i logv)	; set q to quotient  by power of 2 (shift)
      
      (cond ((> (mpz_cmp r vby2) 0)     ; it is a negative coef.
             (if signum (mpz_sub r r v) (mpz_sub r v r))
	     (cond ((= 0 (aref r 1)) (break "~%r=0")))

             (vector-push-extend (gmpz2lisp rz) ansc)
	     (vector-push-extend expon anse)
             (mpz_add i q one))
            (t
	     (cond ((= 0 (aref r 1)	;(mpz_cmp r zero)
		       ;(mpz_get_d r)
		       )nil)		; leave out zeros.
		   (t (if signum nil (mpz_neg r r))
		      (vector-push-extend 
		       ;(mpz_get_d r)  
		       ;; we could use that simple line above if double-float is ok.
		      ;; (round (mpz_get_d r))
		       ;;we could use that simple line above if logv is less than 53,
		       ;; and we want integers.
		       ;; or we can use the little song and dance below, generally
		       (let ((k (mpz_get_d r)))
			 (declare (double-float k))
			 (if (< #.(float(-(expt 2 53)) 1.0d0)
				    k
				    #.(float(-(expt 2 53)) 1.0d0))
				 (round k)
			       (gmpz2lisp rz)))
		       ansc)
		      (vector-push-extend expon anse)))
             (setf i q)))
      (incf expon))
    (decf expon)
    (cons anse	  ansc	  ) )))

;;; this is the version that is actually run in tests with smallish coefficients (< 2^53)..

(defun frombig-small-logv (in logv)
    (declare (fixnum logv)(optimize (speed 3)(safety 0)))
  (let* ((anse (make-array 100 :adjustable t :fill-pointer 0))
	 (ansc (make-array 100 :adjustable t :fill-pointer 0))
         (q (gmpz-z(create_mpz_zero)))
         (rz(create_mpz_zero))
	 (r (gmpz-z rz))
	 (i (gmpz-z in))
         (one (gmpz-z(create_mpz_zero)))        ; will be one
	
	(vby2 (gmpz-z(create_mpz_zero))) ; will be v/2
         (v (gmpz-z(create_mpz_zero)))	; will be v
	 (expon 0)
         (signum (>= (signum (aref i 1)) 0))) ; t if i non-neg
    (declare (fixnum expon))
    (mpz_set_si one 1)                  ; set to one
    (mpz_mul_2exp vby2 one (1- logv))   ; set to v/2
    (mpz_mul_2exp v one logv)           ; set to v
    (if signum nil (mpz_neg i i))       ; make i positive
    (while (> (aref i 1) 0)		; will be 0 when i is zero
      
      (mpz_fdiv_r_2exp r i logv)        ; set r to remainder by power of 2 (mask)
      (mpz_fdiv_q_2exp q i logv)	; set q to quotient  by power of 2 (shift)
      
      (cond ((> (mpz_cmp r vby2) 0)     ; it is a negative coef.
             (if signum (mpz_sub r r v) (mpz_sub r v r))

             (vector-push-extend (gmpz2lisp rz) ansc)
	     (vector-push-extend expon anse)
             (mpz_add i q one))
            (t
	     (cond ;;#+ignore
		   ((= 0 (aref r 1)	;(mpz_cmp r zero) or (mpz_get_d r)
		       )nil)		; leave out zeros.
		   
		   (t (if signum nil (mpz_neg r r))
		      (vector-push-extend 
		       ;;(mpz_get_d r)  
		       ;; we could use that simple line above if double-float is ok.
		       (round (mpz_get_d r))
		       ;;we could use that simple line above if logv is less than 53,
		       ;; and we want integers.
		       ansc)
		      (vector-push-extend expon anse)))
             (setf i q)))
      (incf expon))
    (decf expon)
    (cons anse ansc)))

;;(setf b '(#(1 3 4) . #(10 40 100000)))
;; (mul-kron b b)   
;; (mul-dense b b)
;; (setf c '(#(1 3 4) . #(10 40 100)))  ; also


;; do the whole thing without GMP, in DECIMAL.. This works.


(defun mul-kron10(p q) ;; p, q are polynomials as #(expons) . #(coefs)
  (let* ((v (padpower10 p q))
	 (pr (tobig10  (car p)(cdr p) v))
	 (qr (tobig10  (car q)(cdr q) v)))
    (frombig10 (* qr pr) v)))

(defun padpower10(p1 p2)
  (ceiling (* 0.30103d0 ;; log 2 base 10
	      (+ (integer-length 
		      (min (aref (car p1) (1-(length (car p1))))
			   (aref (car p2) (1-(length (car p2)))))) 
		     ;; the max exponents
		     (integer-length (arraymax (cdr p1)))
		     (integer-length (arraymax (cdr p2)))))))
  
(defun tobig10(a b v)(reduce #'+ (map 'array #'(lambda(r s)(* s (expt 10 (* r v)))) a b)))

(defun frombig10 (i v) 
  (let ((ansc nil)
	(anse nil)
	(expon 0)
	(signum (>= i 10))
	(expt10v (expt 10 v))
	(hexpt10v (expt 10 (1- v))))
    (if signum nil (setf i (- i)))
    (loop (if (= i 0) (return (cons (coerce (nreverse anse) 'array)
				    (coerce (nreverse ansc) 'array))))
      (multiple-value-bind (q r)
	  (truncate i expt10v)
	(cond ((= r 0) (setf i q))
	      ((> r hexpt10v)
	       (setf r (if signum (- r expt10v)(- exprt10v r)))
	       (push r ansc)
	       (push expon anse)
	       (format t "~%i=~s" i)
	       (setf i (+ q 9)))
	      (t (push r ansc)
		 (push expon anse)
		 (setf i q)))
	(incf expon)))))


;; use lisp, but in binary, ok for negative numbers.

(defun mul-kronL(p q) ;; p, q are polynomials as #(expons) . #(coefs)
  (let* ((v (padpower2 p q))
	 (pr (tobigL  (car p)(cdr p) v))
	 (qr (tobigL  (car q)(cdr q) v)))
    (frombigL (* qr pr) v)))

(defun tobigL(a b v)(reduce #'+ (map 'array #'(lambda(r s)(* s (expt 2 (* r v)))) a b)))


(defun frombigL (i v) 
  (let ((ansc nil)
	(anse nil)
	(expon 0)
	(twov (ash 1 v))
	(hv (ash 1  (1- v)))
	(signum (>= i 0)))
    (if signum nil (setf i (- i)))
    (loop (if (= i 0) (return (cons (coerce (nreverse anse)'array)
				    (coerce (nreverse ansc) 'array))))
	  (multiple-value-bind (q r)
	      (truncate i twov)
	    (cond ((= r 0) (setf i q))
		  ((> r hv) 
		   (setf r (if signum (- r twov)(- twov r))) ; compute the coeff
		   (push r ansc)
		   (push expon anse)
		   (setf i (+ q 1)))
		  (t (push r ansc)
		     (push expon anse)
		     (setf i q)))
	    (incf expon)
	    ))))