source: CGBLisp/src/ratpoly.lisp@ 1

#|
        $Id: ratpoly.lisp,v 1.4 2009/01/22 04:07:33 marek Exp $ 
  *--------------------------------------------------------------------------*
  |  Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu)    |
  |    Department of Mathematics, University of Arizona, Tucson, AZ 85721    |
  |                                                                          |
  | Everyone is permitted to copy, distribute and modify the code in this    |
  | directory, as long as this copyright note is preserved verbatim.         |
  *--------------------------------------------------------------------------*
|#
;; Operations in the ring k(x2,x3,...,xn)[x1]
;; The representation is (n1 a1 n2 a2 ... nk ak)
;; where ak are rational functions in an arbitrary number of variables 
;; as in the "rat" package; this representation is referred to as RATPOLY
;;
;; Another representation is POLY1 and this corresponds to reppresenting
;; polynomials as elements of k[x2,x3,...,xn][x1]
;;
;; In order to facilitate input, a function called 
;; poly-to-ratpoly is provided; it converts a polynomial in alist form
;; to a form like above (plist form) with ai being rational
;; Output is done by function ratpoly-print
;;
 
(defpackage "RATPOLY"
  (:export ratpoly+ ratpoly- ratpoly* ratpoly-uminus scalar-times-ratpoly
           rat-times-ratpoly ratpoly-divide ratpoly-remainder
           ratpoly-gcd ratpoly-diff ratpoly-square-free
           ratpoly-normalize ratpoly-resultant
           deg lead ratpoly-discriminant ratpoly-print
           poly-to-ratpoly ratpoly-to-poly poly-resultant)
  (:use "RAT" "PRINTER" "MONOM"  "DIVISION" "TERM" "POLY" "COMMON-LISP"))

(in-package "RATPOLY")

#+debug(proclaim '(optimize (speed 0) (debug 3)))
#-debug(proclaim '(optimize (speed 3) (safety 0)))

;; Arithmetic on polynomials in one variable
;; The polynomial a1 * x^n1 + a2 * x^n2 + ... + ak * x^nk
;; where ai are nonzero numbers and n1>n2>...>nk is represented as list
;;      (n1 a1 n2 a2 ... nk ak)
;; This is the sparse recursive representation

(defun ratpoly+ (p q)
  "Add polynomials P and Q."
  (cond ((endp p) q)
        ((endp q) p)
        ((= (car p) (car q))
         (let ((s (rat+ (cadr p) (cadr q))))
           (if (rat-zerop s)            ;check for cancellation
               (ratpoly+ (cddr p) (cddr q))
             (append (list (car p) s)
                     (ratpoly+ (cddr p) (cddr q))))))
        ((> (car p) (car q))
         (append (list (car p) (cadr p))
                 (ratpoly+ (cddr p) q)))
        (t (append (list (car q) (cadr q))
                   (ratpoly+ p (cddr q))))))

(defun ratpoly- (p q)
  ;;"Subtract polynomial Q from P."
  (cond ((endp p) (ratpoly-uminus q))
        ((endp q) p)
        ((= (car p) (car q))
         (let ((s (rat- (cadr p) (cadr q))))
           (if (rat-zerop s)            ;check for cancellation
               (ratpoly- (cddr p) (cddr q))
             (append (list (car p) s)
                     (ratpoly- (cddr p) (cddr q))))))
        ((> (car p) (car q))
         (append (list (car p) (cadr p))
                 (ratpoly- (cddr p) q)))
        (t (append (list (car q) (rat-uminus (cadr q)))
                   (ratpoly- p (cddr q))))))

(defun ratpoly-uminus (p)
  (cond
   ((endp p) nil)
   (t (cons (car p) (cons (rat-uminus (cadr p)) (ratpoly-uminus (cddr p)))))))

(defun ratpoly* (p q)
  "Multiply polynomials P and Q."
  (cond ((or (endp p) (endp q)) nil)    ;p or q is 0 (represented by NIL)
        ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
        (t (append (list  (+ (car p) (car q)) (rat* (cadr p) (cadr q)))
                   (ratpoly+ (ratpoly* (list (car p) (cadr p)) (cddr q))
                             (ratpoly* (cddr p) q))))))

(defun scalar-times-ratpoly (scalar p)
  "Multiply scalar SCALAR by a polynomial P."
  (cond ((endp p) nil)
        (t (cons (car p)
                 (cons (scalar-times-rat scalar (cadr p))
                       (scalar-times-ratpoly scalar (cddr p)))))))

(defun rat-times-ratpoly (scalar p)
  "Multiply rational function SCALAR by a polynomial P."
  (cond ((endp p) nil)
        (t (cons (car p)
                 (cons (rat* scalar (cadr p))
                       (rat-times-ratpoly scalar (cddr p)))))))

(defun ratpoly-divide (f g)
  "Divide polynomial F by G. Return quotient and remainder as multiple values."
  (do* ((lp nil (- (car r) (car g)))
        (lc nil (rat/ (cadr r) (cadr g)))
        (q nil (nconc q (list lp lc)))
        (r f (ratpoly- (cddr r) (ratpoly* (list lp lc) (cddr g)))))
      ((or (endp r) (< (car r) (car g)))
       (values q r))))

(defun ratpoly-remainder (f g)
  "The remainder of the division of a polynomial F by G."
  (second (multiple-value-list (ratpoly-divide f g))))

(defun ratpoly-gcd (f g)
  "Return GCD of polynomials F and G."
  (do ((h f s)
       (s g (ratpoly-remainder h s)))
      ((endp s) h)))

(defun ratpoly-diff (f)
  "Differentiate a polynomial."
  (cond ((endp f) nil)
        ((zerop (car f)) nil)           ;degree 0
        (t (append (list (1- (car f)) (scalar-times-rat (car f) (cadr f)))
                   (ratpoly-diff (cddr f))))))

(defun ratpoly-square-free (f)
  "Return the square-free part of a polynomial F."
  (ratpoly-divide f (ratpoly-gcd f (ratpoly-diff f))))

(defun ratpoly-normalize (f)
  "Divide a non-zero polynomial by the coefficient at the highest power."
  (rat-times-ratpoly (cons (denom (cadr f)) (num (cadr f))) f))
  
;; A modification of Euclidean algorithm 
(defun ratpoly-resultant (f g)
  "Return the resultant of polynomials F and G."
  (do* ((r nil (ratpoly-remainder h s))
        (res (list 0 (rat-constant 1 (length (car (caaadr f)))))
             (rat-times-ratpoly
              (scalar-times-rat 
               (expt -1 (* (deg h) (deg s)))
               (rat-expt (lead s) (- (deg h) (deg r))))
              res))
        (h f s)
        (s g r))
      ((<= (deg s) 0)
       (cond 
        ((or (endp s) (endp h)) nil)
        ((plusp (deg h)) 
         (rat-times-ratpoly (rat-expt (lead s) (deg h))
                            res))))))

(defun deg (s) 
  (cond 
   ((endp s) 
    #+debugging(warn "ratpoly::deg: Calculating degree of 0 polynomial")
    -1)
   (t (car s))))

(defun lead (s) (cadr s))

(defun ratpoly-discriminant (p &aux (l (deg p)))
  "The discriminant of a polynomial P."
  (rat-times-ratpoly 
   (scalar-div-rat (expt -1 (mod (* l (1- l)) 2)) (lead p))
   (ratpoly-resultant p (ratpoly-diff p))))

(defun ratpoly-print (p vars &optional (stream t) (beg t) (p-orig p))
  (when (endp p)
    (when beg (format stream "0")) 
    (return-from ratpoly-print p-orig))
  (if beg (format stream "((")
    (format stream " + (("))
  (poly-print (caadr p) (rest vars) stream)
  (format stream ") / (")
  (poly-print (cdadr p) (rest vars) stream)
  (case (car p)
    (1  (format stream ")) * ~a" (car vars)))
    (0  (format stream "))" (car vars)))
    (otherwise  
     (format stream ")) * ~a^~d" (car vars) (car p))))
  (ratpoly-print (cddr p) vars stream nil p-orig))

(defun poly-to-ratpoly (p)
  (poly1-to-ratpoly (poly-to-poly1 p)))

;; Convert a polynomial to a polynomial in k[x2,...,xn][x1]
(defun poly-to-poly1 (p &aux (htab (make-hash-table)) q) 
  (dolist (term p)
    (push (cons (cdar term) (cdr term))
          (gethash (caar term) htab nil)))
  (maphash #'(lambda (key val) (push (cons key val) q)) htab)
  (mapcan #'(lambda (x &aux (deg (car x)) (coef (cdr x)))
                      (list deg
                            (sort-poly coef)))
                  (sort q #'> :key #'car)))


;; Convert poly1 to ratpoly, i.e. add denominators=1
(defun poly1-to-ratpoly (p)
  (unless (endp p)
    (cons (car p) 
          (cons 
           (cons (cadr p)
                 (list
                  (cons (make-list (length (caaadr p)) :initial-element 0) 1)))
           (poly1-to-ratpoly (cddr p))))))

(defun ratpoly-to-poly1 (p)
  "Convert every coefficient of ratpoly to polynomial if possible"
  (cond 
   ((endp p) nil) 
   (t (cons (car p)
            (cons (rat-to-poly (cadr p))
                  (ratpoly-to-poly1 (cddr p)))))))
  
(defun poly1-to-poly (p)
  "Convert a ratpoly, whose coeffs have been converted to poly,
into a poly structure, i.e. tack in powers of first variable."
  (cond 
   ((endp p) nil)
   (t
    (append
     (mapcar #'(lambda (x) (cons (cons (car p) (car x)) (cdr x))) (cadr p))
     (poly1-to-poly (cddr p))))))

(defun ratpoly-to-poly (p)
  (poly1-to-poly (ratpoly-to-poly1 p)))
  

(defun poly-resultant (f g)
  "Calculate resultant of F and G given in poly i.e. alist representation." 
  (ratpoly-to-poly
   (ratpoly-resultant (poly-to-ratpoly f)
                      (poly-to-ratpoly g))))

#|
;;----------------------------------------------------------------
;; Multi-variable GCD algorithm
;; Roughly p. 134 of Davenport, Siret, Tournier
;;----------------------------------------------------------------
(defun poly-gcd (A B) 
  (multiple-value-bind (Ap Ac)
      (primitive-part (poly-to-poly1 A))
    (multiple-value-bind (Bp Bc)
        (primitive-part (poly-to-poly1 B))
      (poly* (poly1-to-poly (primitive-part (euclid Ap Bp)))
             (poly1-to-poly 
              (list 0 
                    (cond
                     ((numberp Ac) (list (cons nil (gcd Ac Bc))))
                     (t (poly-gcd Ac Bc)))))))))

;; This operates on poly1
(defun content1 (A)
  (cond 
    ((endp A) (error "content1: Content of 0 is not defined."))
    ((endp (caaadr A))                  ;1-variable poly
     (content1-aux A))                  ;a number
    ((endp (cddr A)) (cadr A))
    (t (poly-gcd (cadr A) (content1 (cddr A))))))

;; 1-variable case
(defun content1-aux (A)
  (cond ((endp A) (error "content1-aux: Content of 0 is not defined."))
        ((endp (cddr A)) (cdaadr A))
        (t (gcd (cdaadr A) (content1-aux (cddr A))))))

;; operates on poly
(defun content (A &aux (A1 (poly-to-poly1 A)))
  (content1 A1))
       
;; Operates on poly1
(defun primitive-part (A)
  (let ((Ac (content1 A)))
    (values (divide-coeffs A Ac) Ac)))

;; Operates on A in poly1 form and A in poly form of 1 variable less
(defun divide-coeffs (A Ac)
  (cond 
    ((endp A) nil)
    ((numberp Ac)               ;ground case
     (cons (car A) 
             (cons (list (cons nil (floor (cdaadr A) Ac)))
                   (divide-coeffs (cddr A) Ac))))
    (t (cons (car A) 
             (cons 
              (multiple-value-bind (q r)
                  (divide (cadr A) (list Ac))
                (unless (endp r) (error "divide-coeffs: not divisible."))
                (car q))
              (divide-coeffs (cddr A) Ac))))))

;; Euclid operates on poly1
(defun euclid (A B)
  (ratpoly-to-poly1 (ratpoly-gcd (poly1-to-ratpoly A)
                                 (poly1-to-ratpoly B))))

|#