source: CGBLisp/src/rat.lisp@ 8

#|
        $Id: rat.lisp,v 1.4 2009/01/22 04:07:12 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.         |
  *--------------------------------------------------------------------------*
|#
;;
;; Arithmetic on rational functions in n variables
;; They are represented as conses of two polynomials, represented 
;; as in the monom package, ordered lexicographically (the default
;; of the monom package).
;; This package is trivial in the sense that it does not produce
;; results in which numerator and denominator are relatively prime.
;; However, zero testing works correctly. This is all that is
;; needed for the resultant package
;;

(defpackage "RAT"
  (:export num denom rat+ rat- rat* rat/ scalar-times-rat scalar-div-rat
           rat-zerop rat-uminus rat-expt rat-constant rat-to-poly
           rat-simplify)
  (:use "ORDER" "POLY-GCD" "DIVISION" "POLY" "COMMON-LISP"))
 
(in-package "RAT")

(proclaim '(optimize (speed 0) (debug 3)))

(defun num (p) (car p))                 ;numerator
(defun denom (p) (cdr p))               ;denominator

(defun rat-simplify-2 (num denom)
  (let ((gcd (poly-gcd num denom)))
    (cons (poly-exact-divide num gcd)
          (poly-exact-divide denom gcd))))

(defun rat-simplify (p)
  (rat-simplify-2 (car p) (cdr p)))

(defun rat+ (p q)
  ;; p1/p2+q1/q2=(p1*q2+p2*q1)/(p2*q2)
  (rat-simplify-2 
   (poly+ (poly* (num p) (denom q))
                (poly* (denom p) (num q)))
   (poly* (denom p) (denom q))))

(defun rat- (p q)
  ;; p1/p2-q1/q2=(p1*q2+p2-q1)/(p2*q2)
  (rat-simplify-2
   (poly- (poly* (num p) (denom q))
          (poly* (denom p) (num q)))
   (poly* (denom p) (denom q))))

;; This could be more efficient if it simplified before multiplication
(defun rat* (p q)
  (rat-simplify-2
   (poly* (num p) (num q))
   (poly* (denom p) (denom q))))

(defun rat/ (p q)
  (rat-simplify-2
   (poly* (num p) (denom q))
   (poly* (denom p) (num q))))

(defun scalar-times-rat (scalar p)
  (cons (scalar-times-poly scalar (num p)) (denom p)))

;; scalar/rational-fun
(defun scalar-div-rat (scalar p)
  (cons (scalar-times-poly scalar (denom p)) (num p)))

(defun rat-zerop (p)
  (endp (car p)))

(defun rat-uminus (p)
  (cons (minus-poly (car p)) (cdr p)))

(defun rat-expt (p n)
  (cons (poly-expt (car p) n) (poly-expt (cdr p) n)))

(defun rat-constant (c n)
  "Make a constant rational function equal to c with n variables"
  (cons (list (cons (make-list n :initial-element 0) c))
        (list (cons (make-list n :initial-element 0) 1))))

(defun rat-to-poly (p)
  "Attempt to convert a rational function to a polynomial by
dividing numerator by denominator. Error if not divisible"
  (poly-exact-divide (car p) (cdr p)))