In this class, we'll look at an extended example of using generative
recursion to implement an algorithm (Gaussian elimination).

Consider a system of equalities:

  2x + 2y + 3z  = 10
  2x + 5y + 12z = 31
  4x + 1y + -2z = 1

In high school, you learned techniques for findings values of x, y,
and z that satisfy these equations.  One of the most common techniques 
is called Gaussian elimination.  Gaussian elimination proceeds
according to the following steps:

1. Convert the top row into an equivalent equation such that the
   coefficient of the first variable is 1.

2. Add a multiple of the first row to each other row in turn, such
   that the resulting row has a 0 coefficient on the first variable.

3. Ignore the first row and repeat these two steps on the remaning
   rows, this time eliminating the next variable in the equation.

So, we would transform the above equation as follows:

  2x + 2y + 3z  = 10
  2x + 5y + 12z = 31
  4x + 1y + -2z = 1

  1x + 1y + 3/2z  = 5
  2x + 5y + 12z = 31
  4x + 1y + -2z = 1

  1x + 1y + 3/2z  = 5
  0x + 3y + 9z = 21
  0x + -3y + -8z = -19

  1x + 1y + 3/2z  = 5
  0x + 1y + 3z = 7
  0x + -3y + -8z = -19

  1x + 1y + 3/2z  = 5
  0x + 1y + 3z = 7
  0x + 0y + 1z = 2

Now, finding values for x, y, and z is easy.  The last equation tells
us that z must be 2.  If we substitute 2 for z in the second equation,
then y must be 1.  If we substitute these values for z and y into the
first equation, we find that x must be 1.  Thus x=1,y=1,z=2 is the
solution to this system of equations.  Note that this set of equations
has a special form: it has a triangle in the lower left corner in
which all coefficients are zero.  These equations are in triangular form.

We want to write a program to perform Gaussian elimination.  First, we 
need a representation of equations.

;; An eqn is a non-empty list of numbers where the last number represents the
;; constant on the right hand side of the equation and the remaining numbers
;; represent coefficients of variables

;; Sample input based on previous example
;;
;;(list (list 2 2 3 10)
;;      (list 2 5 12 31)
;;      (list 4 1 -2 1)))

Now, we need to write two main programs:

1. triangulate: consumes a list of equations and returns a list of
equations.  The returned equations have the same solutions as the
original equations, except they are in triangular form.

2. solve: consumes a list of equations in triangular form and returns
a list of numbers.  The returned numbers are the solutions for the
variables in the equations.

We spent the rest of the class working on these two functions.  Here
are the solutions.  The text also discusses this problem in detail in
section 27.5 (page 385, section called "Gaussian elimination" for
those using the on-line text).

;; subtract-pairwise : eqn eqn -> eqn
;; subtracts the second equation from the first
(define (subtract-pairwise lon1 lon2)
  (cond [(empty? lon1) empty]
        [else (cons (- (first lon1) (first lon2))
                    (subtract-pairwise (rest lon1) (rest lon2)))]))

;; convert-leading-one : eqn -> eqn
;; creates an equivalent equation in which the first coefficient is 1
(define (convert-leading-one an-eqn)
  (local [(define divisor (first an-eqn))]
    (map (lambda (n) (/ n divisor)) an-eqn)))

;; zero-first-column : (listof eqn) -> (listof eqn)
;; creates an equivalent set of equations with zero coefficiants for the 
;; first variable in all but the first equation.
(define (zero-first-column loe)
  (local [(define row1 (first loe))
          ;; mult-eqn-by-factor : eqn num -> eqn
          (define (mult-eqn-by-factor an-eqn factor)
            (map (lambda (n) (* n factor)) an-eqn))]
    (cons row1
          (map (lambda (eqn)
                 (subtract-pairwise 
                  (mult-eqn-by-factor row1 (first eqn))
                  eqn))
               (rest loe)))))

;; triangulate : (listof eqn) -> (listof eqn)
;; convert a list of equations to triangular form
(define (triangulate loe)
  (cond [(empty? loe) empty]
        [else (local [(define zeroed-first
                        (zero-first-column
                         (cons (convert-leading-one (first loe))
                               (rest loe))))]
                (cons (first zeroed-first)
                      (map (lambda (eqn) (cons 0 eqn))
                           (triangulate (map rest (rest zeroed-first))))))]))

;; solve : eqns -> (listof number)
;; input must be triangular set of equations.  Returned list contains values for
;; the coefficients that solves the given system of equations.
(define (solve loe)
  (cond [(empty? loe) empty]
        [else
         (local [(define rest-sols (solve (map rest (rest loe))))
                 (define const-eqn
                   (subst rest-sols (rest (first loe))))]
           (cons (/ (eval const-eqn)
                    (first (first loe)))
                 rest-sols))]))
           
;; eval : (nelistof number) -> number
;; solves sum of constants where last number is on right-hand side and
;;  all other numbers are on left-hand side.  Returned number is on right side.
;; Note: can't use eqn in the contract because this list does not represent 
;;  coefficients any longer
(define (eval const-eqn)
  (cond [(empty? (rest const-eqn)) (first const-eqn)]
        [else (+ (- (first const-eqn))
                 (eval (rest const-eqn)))]))

;; subst : (listof number) eqn -> (listof number)
;; substitutes values for variables in an equantion.
;; the first argument contains values for the variables in the
;;  equation.  The first list is therefore one element shorter than the second
(define (subst sols an-eqn)
  (cond [(empty? sols) an-eqn]
        [else (cons (* (first sols) (first an-eqn))
                    (subst (rest sols) (rest an-eqn)))]))

This solution doesn't work in all cases.  If the first equation in the
list has a zero coefficient on the first variable, the above code will
generate an error (where?).  Thus, the actual algorithm must check for
this condition, and move a row with a non-zero coefficient to the top
of the list before triangulating.  Here's the code for performing the
swap.  We leave modifying the above code to perform the swap
appropriately as an exercise.
  
;; nonzero-row-to-top : (listof eqn) -> (listof eqn)
;; moves some row with non-zero coefficient in first variable to 
;;   first position in the list
(define (nonzero-row-to-top loe)
  (local [(define (find-non-zero-first loe)
            (foldr (lambda (an-eqn rest-eqn)
                     (cond [(not (zero? (first an-eqn))) an-eqn]
                           [else rest-eqn]))
                   false loe))
          (define swap-row (find-non-zero-first loe))]
    (cond [(boolean? swap-row) (error "All rows have zero in first column")]
          [else
           (cons swap-row (remove swap-row loe))])))