2013-02-20

SICP Exercise 2.93: Rational Functions

Modify the rational-arithmetic package to use generic operations, but change make-rat so that it does not attempt to reduce fractions to lowest terms. Test your system by calling make-rational on two polynomials to produce a rational function
(define p1 (make-polynomial 'x '((2 1)(0 1))))
(define p2 (make-polynomial 'x '((3 1)(0 1))))
(define rf (make-rational p2 p1))
Now add rf to itself, using add. You will observe that this addition procedure does not reduce fractions to lowest terms.
Let's start by looking at our current rational package implementation:
(define (install-rational-package)
  ;; internal procedures
  (define (numer x) (car x))
  (define (denom x) (cdr x))
  (define (make-rat n d)
    (if (and (integer? n) (integer? d))
        (let ((g (gcd n d)))
          (cons (/ n g) (/ d g)))
        (error "non-integer numerator or denominator"
               (list n d))))
  (define (add-rat x y)
    (make-rat (+ (* (numer x) (denom y))
                 (* (numer y) (denom x)))
              (* (denom x) (denom y))))
  (define (sub-rat x y)
    (make-rat (- (* (numer x) (denom y))
                 (* (numer y) (denom x)))
              (* (denom x) (denom y))))
  (define (mul-rat x y)
    (make-rat (* (numer x) (numer y))
              (* (denom x) (denom y))))
  (define (div-rat x y)
    (make-rat (* (numer x) (denom y))
              (* (denom x) (numer y))))
  (define (equ? x y)
    (and (= (numer x) (numer y))
         (= (denom x) (denom y))))
  (define (=zero? x) (= (numer x) 0))
  (define (addd x y z)
    (make-rat (+ (* (numer x) (denom y) (denom z))
                 (* (denom x) (numer y) (denom z))
                 (* (denom x) (denom y) (numer z)))
              (* (denom x) (denom y) (denom z))))
  (define (rational->real r) (make-real (/ (numer r) (denom r))))
  (define (rational->integer r) (make-integer (round (/ (numer r) (denom r)))))
  (define (sqrt-rat x)
    (let ((root (sqrt (/ (numer x) (denom x)))))
      (make-complex-from-real-imag (make-real (real-part root))
                                   (make-real (imag-part root)))))
  (define (square-rat x)
    (mul-rat x x))
  (define (arctan-rat x y)
    (atan (/ (numer x) (denom x))
          (/ (numer y) (denom y))))
  (define (cosine-rat x)
    (cos (/ (numer x) (denom x))))
  (define (sine-rat x)
    (sin (/ (numer x) (denom x))))
  (define (negate-rat x)
    (make-rat (- (numer x)) (denom x)))
  ;; interface to rest of the system
  (define (tag x) (attach-tag 'rational x))
  (put 'add '(rational rational)
       (lambda (x y) (tag (add-rat x y))))
  (put 'sub '(rational rational)
       (lambda (x y) (tag (sub-rat x y))))
  (put 'mul '(rational rational)
       (lambda (x y) (tag (mul-rat x y))))
  (put 'div '(rational rational)
       (lambda (x y) (tag (div-rat x y))))
  (put 'equ? '(rational rational) equ?)
  (put '=zero? '(rational) =zero?)
  (put 'addd '(rational rational rational)
       (lambda (x y z) (tag (addd x y z))))
  (put 'sqrt '(rational)
       (lambda (x) (make-real (sqrt-rat x))))
  (put 'square '(rational)
       (lambda (x) (tag (square-rat x))))
  (put 'arctan '(rational rational)
       (lambda (x y) (make-real (arctan-rat x y))))
  (put 'cosine '(rational)
       (lambda (x) (make-real (cosine-rat x))))
  (put 'sine '(rational)
       (lambda (x) (make-real (sine-rat x))))
  (put 'negate '(rational)
       (lambda (x) (tag (negate-rat x))))
  (put 'make 'rational
       (lambda (n d) (tag (make-rat n d))))
  (put-coercion 'rational 'real rational->real)
  (put-coercion 'rational 'integer rational->integer)
  'done)

(define (make-rational n d)
  ((get 'make 'rational) n d))
We've added a lot in there and made a lot of changes since the rational package was first introduced. As a result it's not simply a case of doing a find-and-replace, mapping the primitive operations of +, -, * and / to their generic equivalents of add, sub, mul and div, although that does form part of the modifications.
Let's go through the package and identify the changes we need to make.

Rational Representation

Changing the values we store for the numerator and denominator from primitive integers to tagged types does not affect the overarching representation; we can still use a pair to hold these. As a result numer and denom remain unchanged. However, we do need to modify the creation of the rational representation provided by make-rat:
  • In order to produce mathematically correct rational numbers, make-rat currently restricts us to constructing rational representations using primitive integer numerator and denominator values. Extending the rational package so that it can also be used to represent rational functions requires us to allow polynomials as well. We're doing this by modifying the package so that it uses tagged types for its numerator and denominator. However, while we'll now accept tagged types, in order to ensure mathematical correctness we should ensure that the numerator and denominator are restricted to integer or polynomial types.
  • The integer and real packages include procedures for performing type conversions to the rational type (i.e. integer->rational and real->rational respectively). Both invoke make-rational as part of their type conversion, and pass in primitive integer values. If we're going to modify the rational package to use tagged types and generic operations then our changes need to account for these type conversion procedures. One approach would be to modify the integer->rational and real->rational procedures so that they pass tagged integer values to make-rational. However, we can limit the scope of our changes to the rational package alone if, in addition to integer and polynomial tagged types, we also accept primitive integers as the numerator and/or denominator for make-rat and have make-rat automatically convert them to tagged integer types. To do this we simply use a cond in make-rat, with the predicates in the first two clauses testing respectively whether the numerator or the denominator is a primitive integer, creating a tagged integer and recursively calling with the tagged value if the predicate matches.
These changes look as follows:
(define (numer x) (car x))
(define (denom x) (cdr x))
(define (valid-component? c)
  (memq (type-tag c) '(integer polynomial)))
(define (make-rat n d)
  (cond ((integer? n) (make-rat (make-integer n) d))
        ((integer? d) (make-rat n (make-integer d)))
        ((and (valid-component? n) (valid-component? d)) (cons n d))
        (else (error
               "numerator and denominator must both be integer or polynomial types"
               (list n d)))))
We've introduced the procedure valid-component? here which checks whether a value is tagged as an integer or polynomial. This is used by make-rat to perform the final validation of the types of the passed arguments once it's checked that any primitive integers have been wrapped.

Basic Arithmetic Operations

For the four basic operations of addition, subtraction, multiplication and division a find-and-replace will suffice. We replace + with add, - with sub and * with mul (the primitive division operation, /, is unused here):
(define (add-rat x y)
  (make-rat (add (mul (numer x) (denom y))
                 (mul (numer y) (denom x)))
            (mul (denom x) (denom y))))
(define (sub-rat x y)
  (make-rat (sub (mul (numer x) (denom y))
                 (mul (numer y) (denom x)))
            (mul (denom x) (denom y))))
(define (mul-rat x y)
  (make-rat (mul (numer x) (numer y))
            (mul (denom x) (denom y))))
(define (div-rat x y)
  (make-rat (mul (numer x) (denom y))
            (mul (denom x) (numer y))))
The procedure square-rat simply calls through to mul-rat without inspecting its argument, so it can remain unchanged:
(define (square-rat x)
  (mul-rat x x))
There are a few other basic arithmetic operations, equ?, =zero?, addd and negate-rat, which need only minor tweaks:
  • We can replace = with a call to the generic equ? operation in the equ? procedure. Similarly, we can replace the comparison of the numerator with 0 in =zero? with a call to the generic =zero? on the numerator. However, as these rational operations have identical names to the generic operations, the globally scoped generic operations are obscured by the rational package's implementations. As a result, the interpreter will not invoke the global generic operations with these names but will instead recursively invoke the rational package's implementations, which will then fail as the passed values will not be of the correct type. To resolve this we'll rename the procedures within the rational package to equ-rat? and =zero-rat? respectively and change the corresponding put invocations which install them appropriately.
  • Substituting add and mul for + and * respectively into the procedure addd, which we built as part of exercise 2.82, will not work. The primitive operation * is used as a ternary operation in the current implementation of addd. However, our generic arithmetic operation add is a binary operation. We can address this simply by converting the implementation of addd to use add-rat twice in succession: adding the first two arguments together first, which produces a rational result, and then adding the third argument to the result of that.
  • Similarly, simply substituting sub for - into the procedure negate-rat will not work. In this case - is currently used as a unary operation with the semantics that this negates the primitive integer value passed to it. However, our generic arithmetic operation sub is a binary operation. Fixing this is simpler than for addd: we make the substitution, but we pass 0 as the first argument to sub, and x as the second (as 0 - n = -n).
Here are the updated operations:
(define (equ-rat? x y)
  (and (equ? (numer x) (numer y))
       (equ? (denom x) (denom y))))
(define (=zero-rat? x) (=zero? (numer x)))
(define (addd x y z)
  (add-rat (add-rat x y) z))
(define (negate-rat x)
  (make-rat (sub zero (numer x)) (denom x)))
…
(put 'equ? '(rational rational) equ-rat?)
(put '=zero? '(rational) =zero-rat?)

Type Conversion, Square Root and Trigonometric Operations

The remaining operations are:
  • The type conversion operations: rational->real and rational->integer.
  • The square root operation: sqrt-rat.
  • The trigonometric operations: arctan-rat, cosine-rat and sine-rat.
These all share a couple of somewhat problematic properties:
  1. Applying these operations to polynomials is either complicated, or does not make sense. We can make our lives easier in this exercise if we state that we'll only support these operations for rational numbers; rational functions are not supported.
  2. When applied to integer-tagged values these operations require access to the underlying primitive integer values in order that we can compute their values easily. This breaks the encapsulation of the tagged types but is necessary as all of the operations need to calculate the result of dividing the numerator by the denominator in order to apply a primitive operation to the result. Note that we can't simply use div to calculate this as applying div to the two tagged integer values representing the numerator and denominator will just produce the original tagged rational value we started with! We'll deal with resolving this problem below...

Avoiding the Issue: Using real Operations

Let's deal with sqrt-rat and the trigonometric operations first. Note that, while implementing these directly for rational numbers involves breaking encapsulation, no such problem arises if we first raise a rational number to a real number, and then delegate the calculation to that package! As a result, each of these operations can be achieved by:
  1. Validating that the parameter(s) to the operation are rational numbers and not rational functions (i.e. both the numerator and denominator are integers).
  2. Raising the rational numbers to real numbers. Of course we'll need to re-tag the rational numbers as such before attempting this, otherwise raise won't know what type of number it's dealing with! Note that, assuming we fix our type conversion operations, this raise will be guaranteed to succeed and result in a real number.
  3. Applying the corresponding generic operation (e.g. square-root for sqrt-rat, and so on). The generic operations will then apply the implementations defined for real numbers and so avoid us having to implement them directly!
We can encapsulate the first step of this process in a procedure that checks whether or not all rationals in a list are rational numbers, with integer numerators and denomniators, as follows:
(define (all-rational-numbers? rs)
  (cond ((null? rs) true)
        ((and (eq? 'integer (type-tag (numer (car rs))))
              (eq? 'integer (type-tag (denom (car rs)))))
         (all-rational-numbers? (cdr rs)))
        (else false)))
We can then capture the full process in a procedure which takes a generic operation that is to be applied and a list of rationals to apply it to, tests that they are all rational numbers, tags and raises the rationals and finally applies the generic operation using apply as follows:
(define (apply-as-real-if-all-rational-numbers f rs)
  (if (all-rational-numbers? rs)
      (apply f (map raise (map tag rs)))
      (error (list "numerator and denominator must both be integer to apply " f))))
The final steps are to update sqrt-rat and the trigonometric operations to utilize this procedure, and to change the installation of these procedures via put. The latter is necessary as the results of apply-as-real-if-all-rational-numbers will be a valid tagged value if it succeeds, so there's no need to turn the results into a real number anymore. This allows us to remove the wrapping λ-functions here:
(define (sqrt-rat x)
  (apply-as-real-if-all-rational-numbers square-root (list x)))
(define (arctan-rat x y)
  (apply-as-real-if-all-rational-numbers arctan (list x y)))
(define (cosine-rat x)
  (apply-as-real-if-all-rational-numbers cosine (list x)))
(define (sine-rat x)
  (apply-as-real-if-all-rational-numbers sine (list x)))
…
(put 'sqrt '(rational) sqrt-rat)
(put 'arctan '(rational rational) arctan-rat)
(put 'cosine '(rational) cosine-rat)
(put 'sine '(rational) sine-rat)

Type Conversion and Breaking Encapsulation

All that's left is to deal with the type conversion operations rational->real and rational->integer. While we're only going to support rational numbers here, it's pretty much impossible to avoid the stumbling block that in order to perform these conversions we need to understand the internal representation of the integer type from within the rational package in order that we can perform the primitive division necessary to support the type conversion. I.e. we end up breaking the encapsulation.
In order to contain the scope of the damage I extracted out a procedure, to-primitive, which converts a rational representing a rational number into its equivalent Scheme primitive numerical value. This restricts the breakage to this procedure alone. The two type conversion operations invoke this, blissfully unaware of the issue, and do not need to access the internals of the numerator and denominator in order to perform the division themselves. Note that to-primitive returns 0 if either the numerator or denominator is a polynomial. This is to allow the existing drop behaviour to continue function correctly - if it's not possible to perform a conversion then the value remains unchanged.
Here's the code:
(define (to-primitive x)
  (let ((n (numer x))
        (d (denom x)))
    (if (and (eq? 'integer (type-tag n)) (eq? 'integer (type-tag d)))
        (/ (contents n) (contents d))
        0)))
(define (rational->real r)
      (make-real (to-primitive r)))
(define (rational->integer r)
      (make-integer (round (to-primitive r))))
Yuck! I do apologise!
Of course, if anyone's got a better solution here then please shout!

Rational Functions In Action

Okay, having completed our implementation, let's see it in action... Here's the tests from the book:
> (define p1 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                zero
                                                (make-integer 1))))
> (define p2 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                zero
                                                zero
                                                (make-integer 1))))
> (define rf (make-rational p2 p1))
> (add rf rf)
(rational (polynomial x sparse-terms (term 5 (integer . 2))
                                     (term 3 (integer . 2))
                                     (term 2 (integer . 2))
                                     (term 0 (integer . 2)))
           polynomial x sparse-terms (term 4 (integer . 1))
                                     (term 2 (integer . 2))
                                     (term 0 (integer . 1)))
As you can see this adds the rational function to itself correctly but, as noted in the exercise statement, fails to reduce the result to lowest terms.
Now let's exercise the package with some rational fractions and check we haven't broken anything obvious:
> (define r1 (make-rational 10 3))
> (define r2 (make-rational 3 4))
> (add r1 r2)
(rational (integer . 49) integer . 12)
> (sub r1 r2)
(rational (integer . 31) integer . 12)
> (mul r1 r2)
(rational (integer . 30) integer . 12)
> (div r1 r2)
(rational (integer . 40) integer . 9)
> (equ? r1 r2)
#f
> (equ? r1 r1)
#t
> (=zero? r1)
#f
> (=zero? (make-rational zero 4))
#t
> (addd r1 r1 r1)
(rational (integer . 270) integer . 27)
> (square r2)
(rational (integer . 9) integer . 16)
> (square-root (square r2))
(rational (integer . 3) integer . 4)
> (negate r1)
(rational (integer . -10) integer . 3)
For the trigonometric operations we'll compare the results from our package with the results of Scheme's built-in trigonometric operations. These should be identical as the former ultimately delegates to the latter...
> (arctan r1 r2)
(rational (integer . 3038763055969609) integer . 2251799813685248)
> (exact->inexact (/ 3038763055969609 2251799813685248))
1.34948188444711
> (atan 10/3 3/4)
1.34948188444711
> (cosine r2)
(rational (integer . 6590467434422559) integer . 9007199254740992)
> (exact->inexact (/ 6590467434422559 9007199254740992))
0.731688868873821
> (cos 3/4)
0.731688868873821
> (sine r2)
(rational (integer . 6139656131284749) integer . 9007199254740992)
> (exact->inexact (/ 6139656131284749 9007199254740992))
0.681638760023334
> (sin 3/4)
0.681638760023334
All looks good from our brief test... Now let's move on to writing GCD for polynomials in exercise 2.94!

No comments:

Post a Comment