; File: BinaryGCD.scm
; Author: Keith Schwarz (htiek@cs.stanford.edu)
;
; An algorithm for computing gcd(a, b) using the binary representations of a
; and b.  This algorithm has the same asymptotic complexity as Euclid's
; algorithm, but it can be implemented without division or modulus operations
; by instead relying on shifts.
;
; The algorithm is based on the following recursive structure, based on whether
; the input numbers are odd or even:
;
;   gcd(a, 0)           = a
;   gcd(2n, 2m)         = 2 gcd(n, m)
;   gcd(2n, 2m + 1)     = gcd(n, 2m + 1)
;   gcd(2n + 1, 2m)     = gcd(2n + 1, m)
;   gcd(2n + 1, 2m + 1) = gcd(2n - 2m, 2m + 1)
;
; We prove the correctness of each of these cases individually.  To do so, we
; rely on the fact that gcd(m, n) can be found by writing out each number as a
; product of primes.  The primes common to the representations of both m and n
; are then the gcd of the two numbers.  As an exception to the rule, if either
; number is 0, then the gcd is the other number.
;
; Given this, let us prove the first four cases correct.  First, gcd(a, 0) = a
; is true by the special case of the above rule.  Next, we can see that
; gcd(2n, 2m) = 2 gcd(n, m) must be true because if 2 is a common prime factor
; of the quantities, then the resulting gcd must be the gcd of the numbers with
; 2 removed as a factor from each times 2.  Finally, the two cases
;
;   gcd(2n, 2m + 1) = gcd(n, 2m + 1)
;   gcd(2n + 1, 2m) = gcd(2n + 1, m)
;
; are symmetric to one another.  Since in this case one of the numbers has two
; as a factor and the other does not, then one number has 2 in its prime
; factorization with multiplicity at least one, and the other number does not
; have 2 as a prime factor.  Consequently, 2 will not appear anywhere in the
; two numbers' gcd.  If we remove the 2 from the number that contains it as a
; prime factor, we therefore do not change the gcd.
;
; Finally, we must prove that
;
;   gcd(2n + 1, 2m + 1) = gcd(2n - 2m, 2m + 1)
;
; To see this, consider any divisor of 2n - 2m and 2m + 1.  We claim that this
; number must also be a divisor of 2n + 1 and 2m + 1.  To see this, call that
; number q and let
;
;   2m + 1  = aq
;   2n - 2m = bq
;
; then we have that
;
;   2n + 1 = 2n - 2m + 2m + 1 = aq + bq = (a + b)q
;
; so the number divides 2n + 1.  Similarly, any divisor of 2n + 1 and 2m + 1
; must be a divisor of 2n - 2m, since if
;
;   2n + 1 = cq
;   2m + 1 = dq
;
; then
;
;   2n - 2m = (2n + 1) - (2m + 1) = cq - dq = (c - d)q
;
; Thus gcd(2n + 1, 2m + 1) = gcd(2n - 2m, 2m + 1) as required.
;
; The binary GCD algorithm is quite fast; it always terminates after O(lg ab) 
; steps, where a and b are the numbers whose GCD we want to compute.  To see 
; this, we claim that after any two consecutive applications of the above 
; rules, either a or b will have dropped by at least a factor of two.  
; Consequently, we can only apply O(lg ab) rules before one of the terms
; drops below one and we can apply the first rule as our base case.  To prove
; this claim, consider all the rules in turn.
;
;  - gcd(2n, 2m) = 2 gcd(n, m).  Both n and m drop by a factor of two.
;  - gcd(2n, 2m + 1) = gcd(n, 2m + 1).  n drops by a factor of two.
;  - gcd(2n + 1, 2m) = gcd(2n + 1, m).  m drops by a factor of two.
;  - gcd(2n + 1, 2m + 1) = gcd(2n - 2m, 2m + 1).  This does not necessarily
;    decrease the first term by a factor of two directly.  However, at this
;    point the first term is even and the second odd, so we will apply the
;    second of the above rules to the term to drop it by a factor of two on
;    the next iteration.
;
; Consequently, we need only O(lg ab) steps to compute gcd(a, b), which is on
; par with Euclid's algorithm.

(define binary-gcd 
  (lambda (a b)
    (cond ((< a b) (binary-gcd b a))       ; For simplicity, assume a >= b.
          ((zero? b) a)                    ; gcd(a, 0) = a
          ((and (even? a) (even? b))       ; gcd(2n, 2m) = 2 gcd(n, m)
           (* 2 (binary-gcd (/ a 2) (/ b 2))))
          ((and (even? a) (odd? b))        ; gcd(2n, 2m + 1) = gcd(n, 2m + 1)
           (binary-gcd (/ a 2) b))
          ((and (odd? a) (even? b))        ; gcd(2n + 1, 2m) = gcd(2n + 1, m)
           (binary-gcd a (/ b 2)))
          (else (binary-gcd (- a b) b))))) ; gcd(2n + 1, 2m + 1) = 
                                           ;    gcd(2n - 2m, 2m + 1)