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)))
|
---|