head	1.4;
access;
symbols;
locks; strict;
comment	@;;; @;


1.4
date	2009.01.22.04.07.33;	author marek;	state Exp;
branches;
next	1.3;

1.3
date	2009.01.19.09.30.41;	author marek;	state Exp;
branches;
next	1.2;

1.2
date	2009.01.19.07.52.41;	author marek;	state Exp;
branches;
next	1.1;

1.1
date	2009.01.19.06.46.07;	author marek;	state Exp;
branches;
next	;


desc
@@


1.4
log
@*** empty log message ***
@
text
@#|
	$Id$	
  *--------------------------------------------------------------------------*
  |  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))))

|#

@


1.3
log
@*** empty log message ***
@
text
@d36 2
a37 2
;;(proclaim '(optimize (speed 0) (debug 3)))
(proclaim '(optimize (speed 3) (safety 0)))
@


1.2
log
@*** empty log message ***
@
text
@d36 2
a37 1
(proclaim '(optimize (speed 0) (debug 3)))
@


1.1
log
@Initial revision
@
text
@d2 1
a2 1
	$Id: ratpoly.lisp,v 1.14 1997/12/13 06:50:42 marek Exp $	
d36 2
@