[1] | 1 | #|
|
---|
| 2 | $Id: rat.lisp,v 1.4 2009/01/22 04:07:12 marek Exp $
|
---|
| 3 | *--------------------------------------------------------------------------*
|
---|
| 4 | | Copyright (C) 1994, Marek Rychlik (e-mail: rychlik@math.arizona.edu) |
|
---|
| 5 | | Department of Mathematics, University of Arizona, Tucson, AZ 85721 |
|
---|
| 6 | | |
|
---|
| 7 | | Everyone is permitted to copy, distribute and modify the code in this |
|
---|
| 8 | | directory, as long as this copyright note is preserved verbatim. |
|
---|
| 9 | *--------------------------------------------------------------------------*
|
---|
| 10 | |#
|
---|
| 11 | ;;
|
---|
| 12 | ;; Arithmetic on rational functions in n variables
|
---|
| 13 | ;; They are represented as conses of two polynomials, represented
|
---|
| 14 | ;; as in the monom package, ordered lexicographically (the default
|
---|
| 15 | ;; of the monom package).
|
---|
| 16 | ;; This package is trivial in the sense that it does not produce
|
---|
| 17 | ;; results in which numerator and denominator are relatively prime.
|
---|
| 18 | ;; However, zero testing works correctly. This is all that is
|
---|
| 19 | ;; needed for the resultant package
|
---|
| 20 | ;;
|
---|
| 21 |
|
---|
| 22 | (defpackage "RAT"
|
---|
| 23 | (:export num denom rat+ rat- rat* rat/ scalar-times-rat scalar-div-rat
|
---|
| 24 | rat-zerop rat-uminus rat-expt rat-constant rat-to-poly
|
---|
| 25 | rat-simplify)
|
---|
| 26 | (:use "ORDER" "POLY-GCD" "DIVISION" "POLY" "COMMON-LISP"))
|
---|
| 27 |
|
---|
| 28 | (in-package "RAT")
|
---|
| 29 |
|
---|
| 30 | #+debug(proclaim '(optimize (speed 0) (debug 3)))
|
---|
| 31 | #-debug(proclaim '(optimize (speed 3) (safety 0)))
|
---|
| 32 |
|
---|
| 33 | (defun num (p) (car p)) ;numerator
|
---|
| 34 | (defun denom (p) (cdr p)) ;denominator
|
---|
| 35 |
|
---|
| 36 | (defun rat-simplify-2 (num denom)
|
---|
| 37 | (let ((gcd (poly-gcd num denom)))
|
---|
| 38 | (cons (poly-exact-divide num gcd)
|
---|
| 39 | (poly-exact-divide denom gcd))))
|
---|
| 40 |
|
---|
| 41 | (defun rat-simplify (p)
|
---|
| 42 | (rat-simplify-2 (car p) (cdr p)))
|
---|
| 43 |
|
---|
| 44 | (defun rat+ (p q)
|
---|
| 45 | ;; p1/p2+q1/q2=(p1*q2+p2*q1)/(p2*q2)
|
---|
| 46 | (rat-simplify-2
|
---|
| 47 | (poly+ (poly* (num p) (denom q))
|
---|
| 48 | (poly* (denom p) (num q)))
|
---|
| 49 | (poly* (denom p) (denom q))))
|
---|
| 50 |
|
---|
| 51 | (defun rat- (p q)
|
---|
| 52 | ;; p1/p2-q1/q2=(p1*q2+p2-q1)/(p2*q2)
|
---|
| 53 | (rat-simplify-2
|
---|
| 54 | (poly- (poly* (num p) (denom q))
|
---|
| 55 | (poly* (denom p) (num q)))
|
---|
| 56 | (poly* (denom p) (denom q))))
|
---|
| 57 |
|
---|
| 58 | ;; This could be more efficient if it simplified before multiplication
|
---|
| 59 | (defun rat* (p q)
|
---|
| 60 | (rat-simplify-2
|
---|
| 61 | (poly* (num p) (num q))
|
---|
| 62 | (poly* (denom p) (denom q))))
|
---|
| 63 |
|
---|
| 64 | (defun rat/ (p q)
|
---|
| 65 | (rat-simplify-2
|
---|
| 66 | (poly* (num p) (denom q))
|
---|
| 67 | (poly* (denom p) (num q))))
|
---|
| 68 |
|
---|
| 69 | (defun scalar-times-rat (scalar p)
|
---|
| 70 | (cons (scalar-times-poly scalar (num p)) (denom p)))
|
---|
| 71 |
|
---|
| 72 | ;; scalar/rational-fun
|
---|
| 73 | (defun scalar-div-rat (scalar p)
|
---|
| 74 | (cons (scalar-times-poly scalar (denom p)) (num p)))
|
---|
| 75 |
|
---|
| 76 | (defun rat-zerop (p)
|
---|
| 77 | (endp (car p)))
|
---|
| 78 |
|
---|
| 79 | (defun rat-uminus (p)
|
---|
| 80 | (cons (minus-poly (car p)) (cdr p)))
|
---|
| 81 |
|
---|
| 82 | (defun rat-expt (p n)
|
---|
| 83 | (cons (poly-expt (car p) n) (poly-expt (cdr p) n)))
|
---|
| 84 |
|
---|
| 85 | (defun rat-constant (c n)
|
---|
| 86 | "Make a constant rational function equal to c with n variables"
|
---|
| 87 | (cons (list (cons (make-list n :initial-element 0) c))
|
---|
| 88 | (list (cons (make-list n :initial-element 0) 1))))
|
---|
| 89 |
|
---|
| 90 | (defun rat-to-poly (p)
|
---|
| 91 | "Attempt to convert a rational function to a polynomial by
|
---|
| 92 | dividing numerator by denominator. Error if not divisible"
|
---|
| 93 | (poly-exact-divide (car p) (cdr p)))
|
---|