2013-02-28

SICP Exercise 2.95: GCD and Non-Integer Values

Define P1, P2, and P3 to be the polynomials
P1: x2 - 2x + 1
P2: 11x2 + 7
P3: 13x + 5
Now define Q1 to be the product of P1 and P2 and Q2 to be the product of P1 and P3, and use greatest-common-divisor (exercise 2.94) to compute the GCD of Q1 and Q2. Note that the answer is not the same as P1. This example introduces noninteger operations into the computation, causing difficulties with the GCD algorithm. To understand what is happening, try tracing gcd-terms while computing the GCD or try performing the division by hand.
First let's define P1, P2 and P3 as p1, p2 and p3 respectively:
> (define p1 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 1)
                                                (make-integer -2)
                                                (make-integer 1))))
> (define p2 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 11)
                                                zero
                                                (make-integer 7))))
> (define p3 (make-polynomial-from-coeffs 'x
                                          (list (make-integer 13)
                                                (make-integer 5))))
We then use these to define Q1 and Q2 as q1 and q2 respectively:
> (define q1 (mul p1 p2))
> (define q2 (mul p1 p3))
> q1
(polynomial x dense-terms (integer . 11)
                          (integer . -22)
                          (integer . 18)
                          (integer . -14)
                          (integer . 7))
> q2
(polynomial x dense-terms (integer . 13)
                          (integer . -21)
                          (integer . 3)
                          (integer . 5))
Finally we calculate the GCD of Q1 and Q2:
> (greatest-common-divisor q1 q2)
(polynomial x dense-terms (rational (integer . 18954) integer . 2197)
                          (rational (integer . -37908) integer . 2197)
                          (rational (integer . 1458) integer . 169))
As the authors note in the exercise statement, the answer is not the same as P1, it does introduce non-integer operations into the computation and these then cause difficulties with the GCD algorithm. However, if you look closely at the numbers you might notice a correlation with P1. Note that 2197 / 13 = 169. We can use this to first express all the coefficients as multiples of 1/169:
  (18954/2197)x2 - (37908/2197)x + 1458/169
= ((18954/13)/(2197/13))x2 - ((37908/13)/(2197/13))x + 1458/169
= (1458/169)x2 - (2916/169)x + 1458/169
Now let's normalize this by multiplying all the coefficients by 169/1458:
  ((1458/169)x2 - (2916/169)x + 1458/169) * (169/1458)
= ((1458*169)/(169*1458))x2 - ((2916*169)/(169*1458))x + (1458*169)/(169*1458)
= ((1458*169)/(169*1458))x2 - ((2916*169)/(169*1458))x + (1458*169)/(169*1458)
= (246402/246402)x2 - (492804/246402)x + 246402/246402
= x2 - 2x + 1
= P1
As you can see, the calculated GCD's term coefficients have the same relative ratios as the coefficients of P1, i.e. 1:-2:1.
Let's trace the calculation. The Scheme interpreter I'm using has a trace procedure, which takes a procedure as an argument and turns on tracing of all calls to that passed procedure. However, I found this to be limited in use: when a trace output line is "too long" it arbitrarily collapses each argument passed to the procedure to "#". This is obviously no use for what we'd like to achieve here. As a result, I quickly hacked a trace directly into the procedures gcd-terms and div-terms so that I could get a better picture of what's going on.
With this in place my traced call produced the following output, with the point where non-integer values enter the calculation highlighted in yellow:
> (greatest-common-divisor q1 q2)
gcd-terms:
+-> (dense-terms (integer . 11)
|                (integer . -22)
|                (integer . 18)
|                (integer . -14)
|                (integer . 7))
+-> (dense-terms (integer . 13)
|                (integer . -21)
|                (integer . 3)
|                (integer . 5))
| div-terms:
| +-> (dense-terms (integer . 11)
| |                (integer . -22)
| |                (integer . 18)
| |                (integer . -14)
| |                (integer . 7))
| +-> (dense-terms (integer . 13)
| |                (integer . -21)
| |                (integer . 3)
| |                (integer . 5))
| | div-terms:
| | +-> (dense-terms (rational (integer . -55) integer . 13)
| | |                (rational (integer . 201) integer . 13)
| | |                (rational (integer . -237) integer . 13)
| | |                (integer . 7))
| | +-> (dense-terms (integer . 13)
| | |                (integer . -21)
| | |                (integer . 3)
| | |                (integer . 5))
| | | div-terms:
| | | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | +-> (dense-terms (integer . 13)
| | | |                (integer . -21)
| | | |                (integer . 3)
| | | |                (integer . 5))
| | | +-= ((sparse-terms)
| | |      (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | |                    (term 1 (rational (integer . -37908) integer . 2197))
| | |                    (term 0 (rational (integer . 1458) integer . 169))))
| | +-= ((sparse-terms (term 0 (rational (integer . -55) integer . 169)))
| |      (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| |                    (term 1 (rational (integer . -37908) integer . 2197))
| |                    (term 0 (rational (integer . 1458) integer . 169))))
| +-= ((sparse-terms (term 1 (rational (integer . 11) integer . 13))
|                    (term 0 (rational (integer . -55) integer . 169)))
|      (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
|                    (term 1 (rational (integer . -37908) integer . 2197))
|                    (term 0 (rational (integer . 1458) integer . 169))))
| gcd-terms:
| +-> (dense-terms (integer . 13)
| |                (integer . -21)
| |                (integer . 3)
| |                (integer . 5))
| +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| |                 (term 1 (rational (integer . -37908) integer . 2197))
| |                 (term 0 (rational (integer . 1458) integer . 169)))
| | div-terms:
| | +-> (dense-terms (integer . 13)
| | |                (integer . -21)
| | |                (integer . 3)
| | |                (integer . 5))
| | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | div-terms:
| | | +-> (dense-terms (rational (integer . 208209690) integer . 41641938)
| | | |                (rational (integer . -32032260) integer . 3203226)
| | | |                (integer . 5))
| | | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | | div-terms:
| | | | +-> (sparse-terms)
| | | | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | | | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | | | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | | | +-= ((sparse-terms)
| | | |      (sparse-terms))
| | | +-= ((sparse-terms (term 0 (rational (integer . 457436688930)
| | |                                       integer . 789281292852)))
| | |      (sparse-terms))
| | +-= ((sparse-terms (term 1 (rational (integer . 28561) integer . 18954))
| |                    (term 0 (rational (integer . 457436688930)
| |                                       integer . 789281292852)))
| |      (sparse-terms))
| | gcd-terms:
| | +-> (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| | |                 (term 1 (rational (integer . -37908) integer . 2197))
| | |                 (term 0 (rational (integer . 1458) integer . 169)))
| | +-> (sparse-terms)
| | +-= (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
| |                   (term 1 (rational (integer . -37908) integer . 2197))
| |                   (term 0 (rational (integer . 1458) integer . 169)))
| +-= (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
|                   (term 1 (rational (integer . -37908) integer . 2197))
|                   (term 0 (rational (integer . 1458) integer . 169)))
+-= (sparse-terms (term 2 (rational (integer . 18954) integer . 2197))
                  (term 1 (rational (integer . -37908) integer . 2197))
                  (term 0 (rational (integer . 1458) integer . 169)))
(polynomial x dense-terms (rational (integer . 18954) integer . 2197)
                          (rational (integer . -37908) integer . 2197)
                          (rational (integer . 1458) integer . 169))
As you can see the non-integer values enter the calculation when the first invocation of div-terms by gcd-terms (via remainder-terms) makes a recursive call to itself. The fraction first enters the calculations when div-terms calculates the coefficient to use for the first term of the result. This is calculated by dividing the coefficient of the highest-order term from the first term-list by the coefficient of the highest-order term from the second term-list. Provided the latter is a factor of the former we will avoid non-integer values entering the computation. However in this case the coefficients are 11 and 13 respectively which, being prime numbers, are relatively prime to each other and so share no factors. Not only does this mean the coefficient of the first term of the result will be non-integer, but it also affects the calculations of both the coefficients for the other terms in the result and the remainder.
You may also have noticed in the trace above that the denominators and numerators of the rational numbers that appear in our calculations get exponentially larger as the algorithm progresses. This is a known limitation of Euclidean division, which is the algorithm we're using to implement div-terms. In the next exercise we'll address this by first implementing a procedure which calculates pseudo-remainders, which guarantees that no non-integer coefficents will result from the calculation, and then incorporating this into our gcd-terms implementation.

No comments:

Post a Comment