;; A Lisp program to multiply univariate polynomials.
;; relying on Lisp bignum multiplication. Good if bignum mult is fast.

;;; 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)))
    (frombig (* (tobig p v)
		(tobig q v)) 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

#+ignore
(defun tobig(p v)
  (declare (optimize(speed 3)(safety 0)))
  (reduce #'+ (map 'array #'(lambda(r s)(ash s (* r v))) (car p)(cdr p))))

;; faster version
(defun tobig(p v)(declare (optimize(speed 3)(safety 0)))
  (let ((ans 0)
	(r (car p))			;exponents
	(s (cdr p)))			;coefs
    (loop for i fixnum from (1- (length r)) downto 1 do
	  (setf ans (ash (+ ans(aref s i))
			 (* v (- (aref r i)(aref r (1- i)))))))
    (ash (+ ans (aref s 0))(* v (aref r 0)))))
									

;; convert a big integer into a polynomial, v bits between coefs.
#+ignore
(defun frombig (i v) 
  (let* ((ansc (make-array 10 :adjustable t :fill-pointer 0))
	 (anse (make-array 10 :adjustable t :fill-pointer 0))
	 (expon 0)
	 (hv (ash 1 (1- v)))		;half of 2^v
	 (twov (ash hv 1))		;2^v
	 (signum (>= i 0)))		;keep track of sign
    (if signum nil (setf i (- i)))
    (loop (if (= i 0) (return (cons anse ansc)))
	  (multiple-value-bind (q r)
	      (truncate i twov)	  ;  could be done (faster) by logand and shift
	    (cond ((= r 0) (setf i q))	; don't store zero coefs.
		  ((> r hv)		; r is a negative coefficient
		   (setf r (if signum (- r twov)(- twov r))) ; compute the coeff
		   (vector-push-extend r ansc)   (vector-push-extend expon anse) (setf i (+ q 1)))
		  (t (vector-push-extend r ansc) (vector-push-extend expon anse) (setf i q)))
	    (incf expon)))))


(defun frombig (i v) 
  (let* ((ansc (make-array 10 :adjustable t :fill-pointer 0))
	 (anse (make-array 10 :adjustable t :fill-pointer 0))
	 (expon 0)
	 (hv (ash 1 (1- v)))		;half of 2^v
	 (twov (ash hv 1))		;2^v
	 (mask (1- twov))
	 (signum (>= i 0)))		;keep track of sign
    (if signum nil (setf i (- i)))
    (loop (if (= i 0) (return (cons anse ansc)))
      (let((q (ash i (- v)))
	   (r (logand i mask)))
	    (cond ((= r 0) (setf i q))	; don't store zero coefs.
		  ((> r hv)		; r is a negative coefficient
		   (setf r (if signum (- r twov)(- twov r))) ; compute the coeff
		   (vector-push-extend r ansc)   (vector-push-extend expon anse) (setf i (+ q 1)))
		  (t (vector-push-extend r ansc) (vector-push-extend expon anse) (setf i q)))
	    (incf expon)))))

;; that's all.