;; A Lisp program to multiply univariate polynomials.
;;  relying on GMP multiplication.
;;; use "Kronecker substitution" to multiply polynomials p q
;;; by internally encoding the coefficients of a poly into a single bignum.

(defun mul(p q) ;; p, q are polynomials as #(expons) . #(coefs)
  (let ((v (padpower2 p q)))
    (setf p (tobig p v) q (tobig q v))
    (mpz_mul (gmpz-z p)(gmpz-z p)(gmpz-z q))
    (frombig p v)))

;; some utility programs
;; maximum abs value in an array

(defun arraymax(A)
  (reduce  #'(lambda(r s)(max r (abs s))) A :initial-value 0))

(defun padpower2(p1 p2) ; number of bits in largest coefficient of p1*p2
  (+ 1(integer-length 
	    (* (min (aref (car p1) (1-(length (car p1)))) ; max degree of p1
		    (aref (car p2) (1-(length (car p2)))))
	       (arraymax (cdr p1)) ;max coef of p1
	       (arraymax (cdr p2))))))

;;; some examples
;;;(setf b '(#(1 3 4) . #(10 40 100000)))  is 10*x+40*x^3+100000*x^4
;;;(setf c '(#(1 3 40) . #(10 -40 100000))) is 10*x-40*x^3+100000*x^4
;;;(mul b c) should be
;;; (#(2 5 6 7 41 43 44) . #(100 1000000 -1600 -4000000 1000000 4000000 10000000000))

;; convert a polynomial to a big integer

(defun tobig (p v)(declare (optimize(speed 3)(safety 0)))
       (let* ((ansG (create_mpz_zero))
	      (ans (gmpz-z ansG))
	      (c 0)
	      (r (car p))		;exponents
	      (s (cdr p)))		;coefs
	 (loop for i fixnum from (1- (length r)) downto 1 do
	       ;; checking for fixnum replacing next line with more code,
	       ;;(mpz_add ans ans (gmpz-z (lisp2gmpz (aref s i))))
	       (if (fixnump (setf c (aref s i)))
		   (if (>= c 0) (mpz_add_ui ans ans c)(mpz_sub_ui ans ans (- c)))
		 (mpz_add ans ans (gmpz-z (lisp2gmpz c ))))
	       
	       (mpz_mul_2exp ans ans (* v (- (aref r i)(aref r (1- i))))))
	 (mpz_add ans ans (gmpz-z(lisp2gmpz (aref s 0))))
	 (mpz_mul_2exp ans ans (* v (aref r 0)))
       ansG))

;; convert a big integer into a polynomial, v bits between coefs.

(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))
	 (rz(create_mpz_zero))
	 (r (gmpz-z rz))
	 (q (gmpz-z in))
	(vby2 (gmpz-z(create_mpz_zero))) ; will be v/2
         (v (gmpz-z(create_mpz_zero)))	; will be v
	 (expon 0)
	 (one (gmpz-z(lisp2gmpz 1)))
         (signum (>= (signum (aref q 1)) 0))) ; t if i non-neg
    (declare (fixnum expon))
    (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 q q))       ; make i positive
    (while (> (aref q 1) 0)		; will be 0 when i is zero
      
      (mpz_fdiv_r_2exp r q logv)        ; set r to remainder by power of 2 (mask)
      (mpz_fdiv_q_2exp q q 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_ui q q 1))
            (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)))
            ))
      (incf expon))
    (decf expon)
    (cons anse (nreverse ansc)))))

(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))
         (rz(create_mpz_zero))
	 (r (gmpz-z rz))
	 (q (gmpz-z in))
	(vby2 (gmpz-z(create_mpz_zero))) ; will be v/2
         (v (gmpz-z(create_mpz_zero)))	; will be v
	 (expon 0)
	 	 (one (gmpz-z(lisp2gmpz 1)))
         (signum (>= (signum (aref q 1)) 0))) ; t if i non-neg
    (declare (fixnum expon))
    (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 q q))       ; make i positive
    (while (> (aref q 1) 0)		; will be 0 when i is zero
      
      (mpz_fdiv_r_2exp r q logv)        ; set r to remainder by power of 2 (mask)
      (mpz_fdiv_q_2exp q q 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_ui q q 1))
            (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)))    ))
      (incf expon))
    (decf expon)
    (cons anse (nreverse ansc))))