;;; A Lisp program to multiply univariate polynomials.
;;  relying on GMP multiplication. Specific to Allegro CL, 32bit.
;;  though a small change should work for 64 bit
;;; use "Kronecker substitution" to multiply polynomials p q
;;; by internally encoding the coefficients of a poly into a single bignum.

(eval-when (compile load)
  (declaim (optimize (speed 3) (safety 0) (space 0) (compilation-speed 0)))
  (declaim (inline svref = eql make-array round)))

;; some utility programs
;; maximum abs value in an array
;;(defun arraymax(A)  (reduce  #'(lambda(r s)(max r (abs s))) A :initial-value 0))
;; this version is about 50% faster, reducing many calls to abs.

(defun arraymax (A)
  (declare (simple-array A) (optimize (speed 3)(safety 0)))
  (let ((big 0)(small 0)(this 0))
    (declare (integer big small this))
    (loop for i from 0 below (length A) do
	  (setf this (aref A i))
	  (if (> this big)(setf big this)(if (< this small)(setf small this))))
    (setf this (- small))
    (if (> this big) this big)))

(defun degr(p)(aref (car p) (1-(length (car p)))))


(defun padpower2(p1 p2)			; number of bits in largest coefficient of p1*p2
   (integer-length  
	  (* (+ 1 (max 1 (min (degr p1)	; max degree of p1
			     (degr 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 (fixnum v)(optimize(speed 3)(safety 0)))
       (let* ((ansG (create_mpz_zero))
	      (ans (gmpz-z ansG))
	      (c 0)
	      (r (car p))		;exponents
	      (s (cdr p)))		;coefs
	 (declare (simple-array r s))
	 (loop for i fixnum from (1- (length r)) downto 1 do
	       (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 (L2G c)))) ; c is a bignum
	       (mpz_mul_2exp ans ans (the fixnum(* v (the fixnum (- (the fixnum (aref r i))
								    (the fixnum (aref r (1- i)))))))))
	 (mpz_add ans ans (gmpz-z (L2G (aref s 0)))) ;clean up last coef
	 (mpz_mul_2exp ans ans (* v (aref r 0)))
	 ansG))

;;; convert gmp number to a lisp number
;;; An inefficient hack:
;;(defmethod gmpz2lisp((x gmpz))(parse-integer (create_string_from_mpz x)))
;;; Here is a faster way

(defun G2L (x) 
  (declare (optimize (speed 3)(safety 0)))
  (let* ((addr (aref (gmpz-z x) 2))
	 (sum 0))
    (declare (integer sum)(simple-array addr))
    (loop for i fixnum  from  (ash (1-(abs(aref (gmpz-z x) 1))) 2) downto 0 by 4 do
	  (setf sum (+ (ash sum 32)
		        ;; next line tested only on 32-bit lisp and 32-bit GMP
		       (sys:memref-int  addr 0  i :unsigned-natural))))
    (if (< (aref (gmpz-z x)1) 0) (- sum) sum))) ; change to negative sign if necessary.

#+ignore
(defun G2L(x) ;; fake one that really only copies a number. Not faster though.
  (declare (optimize (speed 3)(safety 0)))
  (let* ((ans (alloc-gmpz2 (ash (abs(aref (gmpz-z x) 1)) 5)))
	 (in (gmpz-z ans)))
    (mpz_add_ui in (gmpz-z x) 0)
    ans))


;; utility to show limbs in a GMP number.

(defun limbs(x)
    (let ((addr (aref (gmpz-z x) 2)))
      (loop for i fixnum  from  0 below (abs(aref (gmpz-z x) 1))
	  collect  (sys:memref-int  addr 0  (ash i 2) :unsigned-natural))))

;; the main program:

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

(eval-when (compile load)(require :llstructs))

(defmacro nth-bigitm (x n) `(sys:memref ,x -14 ,n :unsigned-word))

(defun L2G (x)				;  x is lisp bignum, answer is gmp number
  (declare (optimize (speed 3)(safety 0)))
  (let* ((out (create_mpz_zero))	;the wrapper
	 (r (gmpz-z out))		;the inside
	 (negsign (if (< x 0) t)))
    (if negsign (setf x (- x)))
    (if (fixnump x)(mpz_add_ui r r x)
      (do ((i (1-(excl::bm_size x))(1- i)))
	  ((< i 0) r)
	(declare (fixnum i))
	(mpz_mul_2exp r r 16)
	(mpz_add_ui  r r  (nth-bigitm x (+ i i)))	))
    (if negsign (mpz_neg r r))
    out))

(defvar *gmpone* (L2G 1))


(defun frombig2 (in logv deg) ;; change a bignum in GMP to a polynomial.
  (declare (fixnum logv deg)(optimize (speed 3)(safety 0)))
  (let* ((ansc (make-array (1+ deg))) ;; assume NO COEFS ARE ZERO. arrays are now simple-arrays.
	 (anse (make-array (1+ deg))) ;; exponent array in this case carries no info..
	 (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 *gmpone*))
	 (signum (>= (signum (aref q 1)) 0))) ; t if i non-neg
    (declare (fixnum expon)(simple-array ansc anse))
    (loop for i fixnum from 0 to deg do (setf (aref anse i) i)) ;assume dense
	     
    (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 input positive
    (while (> (aref q 1) 0)  ; will be 0 when input is reduced to 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) 1)	; it is a negative coef.
	     (if signum (mpz_sub r r v) (mpz_sub r v r))
	     (setf (aref ansc expon) (G2L rz))
	     (mpz_add_ui q q 1))
	    (t (if signum nil (mpz_neg r r))
	       (setf (aref ansc expon) (G2L rz))  ))
      (incf expon))
    (cons anse ansc)))


(defun padpowern(p1 n)	; number of bits in largest coefficient of p1^n
   (* n (integer-length  
	 (max 1  (degr p1)))))

(defun pow (p n)			;n integer >=1
  (let ((v (padpowern p n))
	(ansdeg (* n (degr p))))
    (setf p (tobig p v))
    (mpz_pow_ui (gmpz-z p)(gmpz-z p) n)
    (frombig2 p v ansdeg)))