Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

a neat use of bignums: multiplying polynomials by "*"

63 views
Skip to first unread message

Richard Fateman

unread,
Jul 3, 2010, 3:24:03 PM7/3/10
to
A program that some might find amusing;
Richard Fateman (7/2010)


;; A Lisp program to multiply univariate polynomials over the integers.
;; 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
(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)(reduce #'+ (map 'array #'(lambda(r s)(* s (ash 1 (* r
v)))) (car p)(cdr p))))

;; convert a big integer into a polynomial, v bits between coefs.
(defun frombig (i v)
(let* ((ansc nil)
(anse nil)
(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 (coerce (nreverse anse) 'array)
(coerce (nreverse ansc) 'array))))
(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
(push r ansc) (push expon anse) (setf i (+ q 1)))
(t (push r ansc) (push expon anse) (setf i q)))
(incf expon)))))

;; that's all.

0 new messages