[1] | 1 | head 1.4;
|
---|
| 2 | access;
|
---|
| 3 | symbols;
|
---|
| 4 | locks; strict;
|
---|
| 5 | comment @;;; @;
|
---|
| 6 |
|
---|
| 7 |
|
---|
| 8 | 1.4
|
---|
| 9 | date 2009.01.22.04.07.33; author marek; state Exp;
|
---|
| 10 | branches;
|
---|
| 11 | next 1.3;
|
---|
| 12 |
|
---|
| 13 | 1.3
|
---|
| 14 | date 2009.01.19.09.30.41; author marek; state Exp;
|
---|
| 15 | branches;
|
---|
| 16 | next 1.2;
|
---|
| 17 |
|
---|
| 18 | 1.2
|
---|
| 19 | date 2009.01.19.07.52.41; author marek; state Exp;
|
---|
| 20 | branches;
|
---|
| 21 | next 1.1;
|
---|
| 22 |
|
---|
| 23 | 1.1
|
---|
| 24 | date 2009.01.19.06.46.07; author marek; state Exp;
|
---|
| 25 | branches;
|
---|
| 26 | next ;
|
---|
| 27 |
|
---|
| 28 |
|
---|
| 29 | desc
|
---|
| 30 | @@
|
---|
| 31 |
|
---|
| 32 |
|
---|
| 33 | 1.4
|
---|
| 34 | log
|
---|
| 35 | @*** empty log message ***
|
---|
| 36 | @
|
---|
| 37 | text
|
---|
| 38 | @#|
|
---|
| 39 | $Id$
|
---|
| 40 | *--------------------------------------------------------------------------*
|
---|
| 41 | | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@@math.arizona.edu) |
|
---|
| 42 | | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
|
---|
| 43 | | |
|
---|
| 44 | | Everyone is permitted to copy, distribute and modify the code in this |
|
---|
| 45 | | directory, as long as this copyright note is preserved verbatim. |
|
---|
| 46 | *--------------------------------------------------------------------------*
|
---|
| 47 | |#
|
---|
| 48 | ;; Operations in the ring k(x2,x3,...,xn)[x1]
|
---|
| 49 | ;; The representation is (n1 a1 n2 a2 ... nk ak)
|
---|
| 50 | ;; where ak are rational functions in an arbitrary number of variables
|
---|
| 51 | ;; as in the "rat" package; this representation is referred to as RATPOLY
|
---|
| 52 | ;;
|
---|
| 53 | ;; Another representation is POLY1 and this corresponds to reppresenting
|
---|
| 54 | ;; polynomials as elements of k[x2,x3,...,xn][x1]
|
---|
| 55 | ;;
|
---|
| 56 | ;; In order to facilitate input, a function called
|
---|
| 57 | ;; poly-to-ratpoly is provided; it converts a polynomial in alist form
|
---|
| 58 | ;; to a form like above (plist form) with ai being rational
|
---|
| 59 | ;; Output is done by function ratpoly-print
|
---|
| 60 | ;;
|
---|
| 61 |
|
---|
| 62 | (defpackage "RATPOLY"
|
---|
| 63 | (:export ratpoly+ ratpoly- ratpoly* ratpoly-uminus scalar-times-ratpoly
|
---|
| 64 | rat-times-ratpoly ratpoly-divide ratpoly-remainder
|
---|
| 65 | ratpoly-gcd ratpoly-diff ratpoly-square-free
|
---|
| 66 | ratpoly-normalize ratpoly-resultant
|
---|
| 67 | deg lead ratpoly-discriminant ratpoly-print
|
---|
| 68 | poly-to-ratpoly ratpoly-to-poly poly-resultant)
|
---|
| 69 | (:use "RAT" "PRINTER" "MONOM" "DIVISION" "TERM" "POLY" "COMMON-LISP"))
|
---|
| 70 |
|
---|
| 71 | (in-package "RATPOLY")
|
---|
| 72 |
|
---|
| 73 | #+debug(proclaim '(optimize (speed 0) (debug 3)))
|
---|
| 74 | #-debug(proclaim '(optimize (speed 3) (safety 0)))
|
---|
| 75 |
|
---|
| 76 | ;; Arithmetic on polynomials in one variable
|
---|
| 77 | ;; The polynomial a1 * x^n1 + a2 * x^n2 + ... + ak * x^nk
|
---|
| 78 | ;; where ai are nonzero numbers and n1>n2>...>nk is represented as list
|
---|
| 79 | ;; (n1 a1 n2 a2 ... nk ak)
|
---|
| 80 | ;; This is the sparse recursive representation
|
---|
| 81 |
|
---|
| 82 | (defun ratpoly+ (p q)
|
---|
| 83 | "Add polynomials P and Q."
|
---|
| 84 | (cond ((endp p) q)
|
---|
| 85 | ((endp q) p)
|
---|
| 86 | ((= (car p) (car q))
|
---|
| 87 | (let ((s (rat+ (cadr p) (cadr q))))
|
---|
| 88 | (if (rat-zerop s) ;check for cancellation
|
---|
| 89 | (ratpoly+ (cddr p) (cddr q))
|
---|
| 90 | (append (list (car p) s)
|
---|
| 91 | (ratpoly+ (cddr p) (cddr q))))))
|
---|
| 92 | ((> (car p) (car q))
|
---|
| 93 | (append (list (car p) (cadr p))
|
---|
| 94 | (ratpoly+ (cddr p) q)))
|
---|
| 95 | (t (append (list (car q) (cadr q))
|
---|
| 96 | (ratpoly+ p (cddr q))))))
|
---|
| 97 |
|
---|
| 98 | (defun ratpoly- (p q)
|
---|
| 99 | ;;"Subtract polynomial Q from P."
|
---|
| 100 | (cond ((endp p) (ratpoly-uminus q))
|
---|
| 101 | ((endp q) p)
|
---|
| 102 | ((= (car p) (car q))
|
---|
| 103 | (let ((s (rat- (cadr p) (cadr q))))
|
---|
| 104 | (if (rat-zerop s) ;check for cancellation
|
---|
| 105 | (ratpoly- (cddr p) (cddr q))
|
---|
| 106 | (append (list (car p) s)
|
---|
| 107 | (ratpoly- (cddr p) (cddr q))))))
|
---|
| 108 | ((> (car p) (car q))
|
---|
| 109 | (append (list (car p) (cadr p))
|
---|
| 110 | (ratpoly- (cddr p) q)))
|
---|
| 111 | (t (append (list (car q) (rat-uminus (cadr q)))
|
---|
| 112 | (ratpoly- p (cddr q))))))
|
---|
| 113 |
|
---|
| 114 | (defun ratpoly-uminus (p)
|
---|
| 115 | (cond
|
---|
| 116 | ((endp p) nil)
|
---|
| 117 | (t (cons (car p) (cons (rat-uminus (cadr p)) (ratpoly-uminus (cddr p)))))))
|
---|
| 118 |
|
---|
| 119 | (defun ratpoly* (p q)
|
---|
| 120 | "Multiply polynomials P and Q."
|
---|
| 121 | (cond ((or (endp p) (endp q)) nil) ;p or q is 0 (represented by NIL)
|
---|
| 122 | ;; If p=p0+p1 and q=q0+q1 then pq=p0q0+p0q1+p1q
|
---|
| 123 | (t (append (list (+ (car p) (car q)) (rat* (cadr p) (cadr q)))
|
---|
| 124 | (ratpoly+ (ratpoly* (list (car p) (cadr p)) (cddr q))
|
---|
| 125 | (ratpoly* (cddr p) q))))))
|
---|
| 126 |
|
---|
| 127 | (defun scalar-times-ratpoly (scalar p)
|
---|
| 128 | "Multiply scalar SCALAR by a polynomial P."
|
---|
| 129 | (cond ((endp p) nil)
|
---|
| 130 | (t (cons (car p)
|
---|
| 131 | (cons (scalar-times-rat scalar (cadr p))
|
---|
| 132 | (scalar-times-ratpoly scalar (cddr p)))))))
|
---|
| 133 |
|
---|
| 134 | (defun rat-times-ratpoly (scalar p)
|
---|
| 135 | "Multiply rational function SCALAR by a polynomial P."
|
---|
| 136 | (cond ((endp p) nil)
|
---|
| 137 | (t (cons (car p)
|
---|
| 138 | (cons (rat* scalar (cadr p))
|
---|
| 139 | (rat-times-ratpoly scalar (cddr p)))))))
|
---|
| 140 |
|
---|
| 141 | (defun ratpoly-divide (f g)
|
---|
| 142 | "Divide polynomial F by G. Return quotient and remainder as multiple values."
|
---|
| 143 | (do* ((lp nil (- (car r) (car g)))
|
---|
| 144 | (lc nil (rat/ (cadr r) (cadr g)))
|
---|
| 145 | (q nil (nconc q (list lp lc)))
|
---|
| 146 | (r f (ratpoly- (cddr r) (ratpoly* (list lp lc) (cddr g)))))
|
---|
| 147 | ((or (endp r) (< (car r) (car g)))
|
---|
| 148 | (values q r))))
|
---|
| 149 |
|
---|
| 150 | (defun ratpoly-remainder (f g)
|
---|
| 151 | "The remainder of the division of a polynomial F by G."
|
---|
| 152 | (second (multiple-value-list (ratpoly-divide f g))))
|
---|
| 153 |
|
---|
| 154 | (defun ratpoly-gcd (f g)
|
---|
| 155 | "Return GCD of polynomials F and G."
|
---|
| 156 | (do ((h f s)
|
---|
| 157 | (s g (ratpoly-remainder h s)))
|
---|
| 158 | ((endp s) h)))
|
---|
| 159 |
|
---|
| 160 | (defun ratpoly-diff (f)
|
---|
| 161 | "Differentiate a polynomial."
|
---|
| 162 | (cond ((endp f) nil)
|
---|
| 163 | ((zerop (car f)) nil) ;degree 0
|
---|
| 164 | (t (append (list (1- (car f)) (scalar-times-rat (car f) (cadr f)))
|
---|
| 165 | (ratpoly-diff (cddr f))))))
|
---|
| 166 |
|
---|
| 167 | (defun ratpoly-square-free (f)
|
---|
| 168 | "Return the square-free part of a polynomial F."
|
---|
| 169 | (ratpoly-divide f (ratpoly-gcd f (ratpoly-diff f))))
|
---|
| 170 |
|
---|
| 171 | (defun ratpoly-normalize (f)
|
---|
| 172 | "Divide a non-zero polynomial by the coefficient at the highest power."
|
---|
| 173 | (rat-times-ratpoly (cons (denom (cadr f)) (num (cadr f))) f))
|
---|
| 174 |
|
---|
| 175 | ;; A modification of Euclidean algorithm
|
---|
| 176 | (defun ratpoly-resultant (f g)
|
---|
| 177 | "Return the resultant of polynomials F and G."
|
---|
| 178 | (do* ((r nil (ratpoly-remainder h s))
|
---|
| 179 | (res (list 0 (rat-constant 1 (length (car (caaadr f)))))
|
---|
| 180 | (rat-times-ratpoly
|
---|
| 181 | (scalar-times-rat
|
---|
| 182 | (expt -1 (* (deg h) (deg s)))
|
---|
| 183 | (rat-expt (lead s) (- (deg h) (deg r))))
|
---|
| 184 | res))
|
---|
| 185 | (h f s)
|
---|
| 186 | (s g r))
|
---|
| 187 | ((<= (deg s) 0)
|
---|
| 188 | (cond
|
---|
| 189 | ((or (endp s) (endp h)) nil)
|
---|
| 190 | ((plusp (deg h))
|
---|
| 191 | (rat-times-ratpoly (rat-expt (lead s) (deg h))
|
---|
| 192 | res))))))
|
---|
| 193 |
|
---|
| 194 | (defun deg (s)
|
---|
| 195 | (cond
|
---|
| 196 | ((endp s)
|
---|
| 197 | #+debugging(warn "ratpoly::deg: Calculating degree of 0 polynomial")
|
---|
| 198 | -1)
|
---|
| 199 | (t (car s))))
|
---|
| 200 |
|
---|
| 201 | (defun lead (s) (cadr s))
|
---|
| 202 |
|
---|
| 203 | (defun ratpoly-discriminant (p &aux (l (deg p)))
|
---|
| 204 | "The discriminant of a polynomial P."
|
---|
| 205 | (rat-times-ratpoly
|
---|
| 206 | (scalar-div-rat (expt -1 (mod (* l (1- l)) 2)) (lead p))
|
---|
| 207 | (ratpoly-resultant p (ratpoly-diff p))))
|
---|
| 208 |
|
---|
| 209 | (defun ratpoly-print (p vars &optional (stream t) (beg t) (p-orig p))
|
---|
| 210 | (when (endp p)
|
---|
| 211 | (when beg (format stream "0"))
|
---|
| 212 | (return-from ratpoly-print p-orig))
|
---|
| 213 | (if beg (format stream "((")
|
---|
| 214 | (format stream " + (("))
|
---|
| 215 | (poly-print (caadr p) (rest vars) stream)
|
---|
| 216 | (format stream ") / (")
|
---|
| 217 | (poly-print (cdadr p) (rest vars) stream)
|
---|
| 218 | (case (car p)
|
---|
| 219 | (1 (format stream ")) * ~a" (car vars)))
|
---|
| 220 | (0 (format stream "))" (car vars)))
|
---|
| 221 | (otherwise
|
---|
| 222 | (format stream ")) * ~a^~d" (car vars) (car p))))
|
---|
| 223 | (ratpoly-print (cddr p) vars stream nil p-orig))
|
---|
| 224 |
|
---|
| 225 | (defun poly-to-ratpoly (p)
|
---|
| 226 | (poly1-to-ratpoly (poly-to-poly1 p)))
|
---|
| 227 |
|
---|
| 228 | ;; Convert a polynomial to a polynomial in k[x2,...,xn][x1]
|
---|
| 229 | (defun poly-to-poly1 (p &aux (htab (make-hash-table)) q)
|
---|
| 230 | (dolist (term p)
|
---|
| 231 | (push (cons (cdar term) (cdr term))
|
---|
| 232 | (gethash (caar term) htab nil)))
|
---|
| 233 | (maphash #'(lambda (key val) (push (cons key val) q)) htab)
|
---|
| 234 | (mapcan #'(lambda (x &aux (deg (car x)) (coef (cdr x)))
|
---|
| 235 | (list deg
|
---|
| 236 | (sort-poly coef)))
|
---|
| 237 | (sort q #'> :key #'car)))
|
---|
| 238 |
|
---|
| 239 |
|
---|
| 240 | ;; Convert poly1 to ratpoly, i.e. add denominators=1
|
---|
| 241 | (defun poly1-to-ratpoly (p)
|
---|
| 242 | (unless (endp p)
|
---|
| 243 | (cons (car p)
|
---|
| 244 | (cons
|
---|
| 245 | (cons (cadr p)
|
---|
| 246 | (list
|
---|
| 247 | (cons (make-list (length (caaadr p)) :initial-element 0) 1)))
|
---|
| 248 | (poly1-to-ratpoly (cddr p))))))
|
---|
| 249 |
|
---|
| 250 | (defun ratpoly-to-poly1 (p)
|
---|
| 251 | "Convert every coefficient of ratpoly to polynomial if possible"
|
---|
| 252 | (cond
|
---|
| 253 | ((endp p) nil)
|
---|
| 254 | (t (cons (car p)
|
---|
| 255 | (cons (rat-to-poly (cadr p))
|
---|
| 256 | (ratpoly-to-poly1 (cddr p)))))))
|
---|
| 257 |
|
---|
| 258 | (defun poly1-to-poly (p)
|
---|
| 259 | "Convert a ratpoly, whose coeffs have been converted to poly,
|
---|
| 260 | into a poly structure, i.e. tack in powers of first variable."
|
---|
| 261 | (cond
|
---|
| 262 | ((endp p) nil)
|
---|
| 263 | (t
|
---|
| 264 | (append
|
---|
| 265 | (mapcar #'(lambda (x) (cons (cons (car p) (car x)) (cdr x))) (cadr p))
|
---|
| 266 | (poly1-to-poly (cddr p))))))
|
---|
| 267 |
|
---|
| 268 | (defun ratpoly-to-poly (p)
|
---|
| 269 | (poly1-to-poly (ratpoly-to-poly1 p)))
|
---|
| 270 |
|
---|
| 271 |
|
---|
| 272 | (defun poly-resultant (f g)
|
---|
| 273 | "Calculate resultant of F and G given in poly i.e. alist representation."
|
---|
| 274 | (ratpoly-to-poly
|
---|
| 275 | (ratpoly-resultant (poly-to-ratpoly f)
|
---|
| 276 | (poly-to-ratpoly g))))
|
---|
| 277 |
|
---|
| 278 | #|
|
---|
| 279 | ;;----------------------------------------------------------------
|
---|
| 280 | ;; Multi-variable GCD algorithm
|
---|
| 281 | ;; Roughly p. 134 of Davenport, Siret, Tournier
|
---|
| 282 | ;;----------------------------------------------------------------
|
---|
| 283 | (defun poly-gcd (A B)
|
---|
| 284 | (multiple-value-bind (Ap Ac)
|
---|
| 285 | (primitive-part (poly-to-poly1 A))
|
---|
| 286 | (multiple-value-bind (Bp Bc)
|
---|
| 287 | (primitive-part (poly-to-poly1 B))
|
---|
| 288 | (poly* (poly1-to-poly (primitive-part (euclid Ap Bp)))
|
---|
| 289 | (poly1-to-poly
|
---|
| 290 | (list 0
|
---|
| 291 | (cond
|
---|
| 292 | ((numberp Ac) (list (cons nil (gcd Ac Bc))))
|
---|
| 293 | (t (poly-gcd Ac Bc)))))))))
|
---|
| 294 |
|
---|
| 295 | ;; This operates on poly1
|
---|
| 296 | (defun content1 (A)
|
---|
| 297 | (cond
|
---|
| 298 | ((endp A) (error "content1: Content of 0 is not defined."))
|
---|
| 299 | ((endp (caaadr A)) ;1-variable poly
|
---|
| 300 | (content1-aux A)) ;a number
|
---|
| 301 | ((endp (cddr A)) (cadr A))
|
---|
| 302 | (t (poly-gcd (cadr A) (content1 (cddr A))))))
|
---|
| 303 |
|
---|
| 304 | ;; 1-variable case
|
---|
| 305 | (defun content1-aux (A)
|
---|
| 306 | (cond ((endp A) (error "content1-aux: Content of 0 is not defined."))
|
---|
| 307 | ((endp (cddr A)) (cdaadr A))
|
---|
| 308 | (t (gcd (cdaadr A) (content1-aux (cddr A))))))
|
---|
| 309 |
|
---|
| 310 | ;; operates on poly
|
---|
| 311 | (defun content (A &aux (A1 (poly-to-poly1 A)))
|
---|
| 312 | (content1 A1))
|
---|
| 313 |
|
---|
| 314 | ;; Operates on poly1
|
---|
| 315 | (defun primitive-part (A)
|
---|
| 316 | (let ((Ac (content1 A)))
|
---|
| 317 | (values (divide-coeffs A Ac) Ac)))
|
---|
| 318 |
|
---|
| 319 | ;; Operates on A in poly1 form and A in poly form of 1 variable less
|
---|
| 320 | (defun divide-coeffs (A Ac)
|
---|
| 321 | (cond
|
---|
| 322 | ((endp A) nil)
|
---|
| 323 | ((numberp Ac) ;ground case
|
---|
| 324 | (cons (car A)
|
---|
| 325 | (cons (list (cons nil (floor (cdaadr A) Ac)))
|
---|
| 326 | (divide-coeffs (cddr A) Ac))))
|
---|
| 327 | (t (cons (car A)
|
---|
| 328 | (cons
|
---|
| 329 | (multiple-value-bind (q r)
|
---|
| 330 | (divide (cadr A) (list Ac))
|
---|
| 331 | (unless (endp r) (error "divide-coeffs: not divisible."))
|
---|
| 332 | (car q))
|
---|
| 333 | (divide-coeffs (cddr A) Ac))))))
|
---|
| 334 |
|
---|
| 335 | ;; Euclid operates on poly1
|
---|
| 336 | (defun euclid (A B)
|
---|
| 337 | (ratpoly-to-poly1 (ratpoly-gcd (poly1-to-ratpoly A)
|
---|
| 338 | (poly1-to-ratpoly B))))
|
---|
| 339 |
|
---|
| 340 | |#
|
---|
| 341 |
|
---|
| 342 | @
|
---|
| 343 |
|
---|
| 344 |
|
---|
| 345 | 1.3
|
---|
| 346 | log
|
---|
| 347 | @*** empty log message ***
|
---|
| 348 | @
|
---|
| 349 | text
|
---|
| 350 | @d36 2
|
---|
| 351 | a37 2
|
---|
| 352 | ;;(proclaim '(optimize (speed 0) (debug 3)))
|
---|
| 353 | (proclaim '(optimize (speed 3) (safety 0)))
|
---|
| 354 | @
|
---|
| 355 |
|
---|
| 356 |
|
---|
| 357 | 1.2
|
---|
| 358 | log
|
---|
| 359 | @*** empty log message ***
|
---|
| 360 | @
|
---|
| 361 | text
|
---|
| 362 | @d36 2
|
---|
| 363 | a37 1
|
---|
| 364 | (proclaim '(optimize (speed 0) (debug 3)))
|
---|
| 365 | @
|
---|
| 366 |
|
---|
| 367 |
|
---|
| 368 | 1.1
|
---|
| 369 | log
|
---|
| 370 | @Initial revision
|
---|
| 371 | @
|
---|
| 372 | text
|
---|
| 373 | @d2 1
|
---|
| 374 | a2 1
|
---|
| 375 | $Id: ratpoly.lisp,v 1.14 1997/12/13 06:50:42 marek Exp $
|
---|
| 376 | d36 2
|
---|
| 377 | @
|
---|