Version 33 of Greatest common denominator

Updated 2007-03-08 20:32:23

The most common and usually most efficient algorithm to compute the GCD is called Euclid's algorithm and is based on finding the remainder of the smaller number repeatedly subtracted from the larger. It appears in Euclid's Elements [L1 ] (Book 7, Propositions 1 and 2 [L2 ]). It may have been derived from an earlier algorithm due to Eudoxus around 375 BCE. Euclid's algorithm may hold the distinction of being the oldest nontrivial algorithm (the Russian peasant algorithm used for multiplication may be older). Binary GCD below gives a non Euclid algorithm for computing GCD. KPV

The running time for Euclid's algorithm is: O(log(q/gcd(p,q))) log is based phi. (For details see Introduction to Algorithms by Cormen, Leiserson and Rivest[L3 ].)

(moved out of the Bag of algorithms) - RS started this with:

  proc gcdA {p q} {
     while {1} {
        if {$p==$q} {
            # Termination case
            return $p
        } elseif {$p>$q} {
            # Swap p and q
            set t $p
            set p $q
            set q $t
        }
        # Work out q%p ('cos we're looping...)
        incr q [expr {-$p}]
     }
  } ;#RS, with recursion removed by DKF

Lars H: It is simpler and faster to use % (remainder) instead of decrementing. (Consider what happens when the above tries to compute the GCD of 1000 and 10. OTOH, Euclid indeed stated everything in terms of repeated subtractions, so that has some historical accuracy credit.) Furthermore the above control structures can be greatly simplified:

  proc gcdB {p q} {
     set p [expr {abs($p)}]
     set q [expr {abs($q)}]
     while {$q != 0} {
        set r [expr {$p % $q}]
        set p $q
        set q $r
     }
     set p
  }

There is however an even shorter way of doing it:

  proc gcdC {p q} {
     set p [expr {abs($p)}]
     set q [expr {abs($q)}]
     while {$q != 0} {set q [expr {$p % [set p $q]}]}
     set p
  }

CL wants to flatten that even more:

   proc gcdE {p q} {
      while {$q != 0} {set q [expr {abs($p) % [set p [expr abs($q)]]}]}
      set p
   } 

KPV how about

   proc gcdH {p q} {
     while {$q != 0} {set q [expr {$p % [set p $q]}]}
     expr {abs($p)}
   }

Cameron Laird wrote in the comp.lang.tcl newsgroup: Simpler:

  proc gcdD {a b} {
     if {$b==0} {
        return $a
     }
     gcdD $b [expr {$a % $b}]
  }

RS has this expr variation:

 proc gcdE {a b} {expr {$b==0? $a: [gcdE $b [expr {$a%$b}]]}}

You can also do it iteratively, with arbitrary-precision integers, ...

Tom Poindexter added: ... e.g., by using Mpexpr, it's a one-liner:

 % package require Mpexpr
 1.0
 % mpformat %r [mpexpr 60/100.0]
 3/5

Mpexpr: ftp://ftp.procplace.com/pub/tcl/sorted/math/Mpexpr-1.0/


Extended Euclid's Algorithm KPV

It is easy to rewrite the algorithm to include additional useful information. Specifically, the integer coefficients x and y such that

   d = gcd(p,q) = px + qy

 proc Xgcd {p q} {
    if {$q == 0} {
        return [list $p 1 0]
    }
    foreach {d x y} [Xgcd $q [expr {$p % $q}]] break
    return [list $d $y [expr {$x - $y * int($p/$q)}]]
 }

On a theoretical note, the analysis of Euclid's algorithm brings up a surprising connection to Fibonacci numbers. It turns out the worst case input for the algorithm are two consecutive Fibonacci numbers.


Binary GCD - a non Euclid algorithm GCD KPV

 GCD(p,q) =  if p and q are both even           => gcd(p,q) = 2gcd(p/2,q/2)
             if p is odd and q is even          => gcd(p,q) = gcd(p,q/2)
             if p is even and q is odd          => gcd(p,q) = gcd(p/2,q)
             if p and q are both odd            => gcd(p,q) = gcd((p-q)/2,q)

Here's a iterative implementation of the above description:

 proc binary-gcd {a b} {
    if {$a == $b} { return $a }

    set doubles 1
    while {$a != $b} {
        if {$a < $b} { foreach {a b} [list $b $a] break }
        if {$a & 1} {                           
            if {$b & 1} {                       
                set a [expr {($a + $b) / 2}]    ;# A & B odd
                continue
            }
            set b [expr {$b/2}]                 ;# A odd, B even
            continue
        }
        if {$b & 1} {                           ;# A even, B odd
            set a [expr {$a/2}]
            continue
        }
        set doubles [expr {$doubles * 2}]       ;# A and b even
        set a [expr {$a/2}]
        set b [expr {$b/2}]
    }
    return [expr {$a * $doubles}]
 }

This method avoids the modulus computation found in the other methods, and only uses parity tests, subtractions and shifts. On some computers, this may be faster. In tcl, binary-gcd was 10 times slower than gcdC.


Closely related to the Extended Euclid Algorithm are the various procedures used to compute convergents to a regular continued fraction. Fraction Math has one that computes the nearest ratio of two integers to an arbitrary floating-point number, and Notes on continued fractions has more, including the continued-fraction representations of quadratic surds.


Performance - I (BBH) changed the names of the above to be unique and ran this little performance comparison:

 set nums {10 11 15 17 35 95 160 1000 1111 1113}
 puts "  each block contains the time for \[gcdX \$p \$q]/\[gcdX \$q \$p] for 5000 iterations"
 puts " |-----------------------------------------------------------|"
 puts " |   P & Q   |    gcdA   |    gcdB   |    gcdC   |    gcdD   |"
 set i 0
 foreach p $nums {
     incr i
     foreach q [lrange $nums $i end] {
         puts " |-----------+-----------+-----------+-----------+-----------|"
         puts -nonewline  [format " | %4d %4d |" $p $q]
         foreach f {A B C D} {
             # time for p q AND for Q P to see if it makes any difference
             puts -nonewline [format " %4d/%4d |" [lindex [time {gcd$f $p $q} 5000] 0] [lindex [time {gcd$f $q $p} 5000] 0]]
         }
         puts ""
         update
     }
 }
 puts " |-----------------------------------------------------------|"

and the results on my PIII running Win2K tcl 8.3.4 on 23 Jan 2002 were

 |-----------------------------------------------------------|
 |   P & Q   |    gcdA   |    gcdB   |    gcdC   |    gcdD   |
 |-----------+-----------+-----------+-----------+-----------|
 |   10   11 |   52/  56 |   28/  22 |   24/  20 |   28/  22 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10   15 |   26/  28 |   26/  24 |   26/  20 |   28/  24 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10   17 |   44/  46 |   32/  28 |   28/  26 |   36/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10   35 |   34/  36 |   26/  22 |   24/  20 |   28/  24 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10   95 |   54/  54 |   26/  22 |   24/  20 |   28/  22 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10  160 |   68/  68 |   22/  20 |   20/  20 |   22/  20 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10 1000 |  341/ 350 |   24/  20 |   22/  18 |   22/  20 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10 1111 |  415/ 437 |   24/  24 |   22/  20 |   28/  22 |
 |-----------+-----------+-----------+-----------+-----------|
 |   10 1113 |  401/ 403 |   28/  28 |   24/  22 |   32/  28 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11   15 |   44/  46 |   30/  30 |   28/  24 |   36/  30 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11   17 |   46/  50 |   32/  28 |   28/  24 |   36/  34 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11   35 |   54/  52 |   30/  26 |   24/  24 |   30/  28 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11   95 |   70/  70 |   34/  32 |   28/  28 |   40/  36 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11  160 |   90/  92 |   32/  28 |   26/  24 |   34/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11 1000 |  348/ 351 |   28/  26 |   24/  22 |   32/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11 1111 |  364/ 365 |   22/  20 |   20/  18 |   24/  18 |
 |-----------+-----------+-----------+-----------+-----------|
 |   11 1113 |  372/ 373 |   30/  24 |   28/  22 |   36/  30 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15   17 |   54/  58 |   28/  28 |   24/  22 |   34/  28 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15   35 |   34/  34 |   26/  24 |   22/  24 |   34/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15   95 |   54/  48 |   26/  24 |   22/  20 |   28/  24 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15  160 |   62/  62 |   30/  26 |   24/  22 |   32/  28 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15 1000 |  246/ 244 |   30/  24 |   26/  22 |   32/  26 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15 1111 |  306/ 310 |   26/  24 |   24/  20 |   26/  24 |
 |-----------+-----------+-----------+-----------+-----------|
 |   15 1113 |  300/ 294 |   30/  22 |   22/  22 |   26/  24 |
 |-----------+-----------+-----------+-----------+-----------|
 |   17   35 |   80/  82 |   24/  24 |   22/  20 |   28/  22 |
 |-----------+-----------+-----------+-----------+-----------|
 |   17   95 |   64/  64 |   38/  32 |   30/  28 |   40/  36 |
 |-----------+-----------+-----------+-----------+-----------|
 |   17  160 |   72/  76 |   32/  28 |   28/  24 |   36/  30 |
 |-----------+-----------+-----------+-----------+-----------|
 |   17 1000 |  238/ 238 |   34/  32 |   28/  28 |   40/  36 |
 |-----------+-----------+-----------+-----------+-----------|
 |   17 1111 |  258/ 260 |   32/  28 |   26/  26 |   36/  30 |
 |-----------+-----------+-----------+-----------+-----------|
 |   17 1113 |  262/ 266 |   28/  26 |   24/  24 |   32/  26 |
 |-----------+-----------+-----------+-----------+-----------|
 |   35   95 |   44/  44 |   32/  30 |   26/  26 |   36/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   35  160 |   54/  68 |   36/  30 |   26/  26 |   36/  30 |
 |-----------+-----------+-----------+-----------+-----------|
 |   35 1000 |  128/ 130 |   32/  28 |   26/  26 |   34/  30 |
 |-----------+-----------+-----------+-----------+-----------|
 |   35 1111 |  162/ 166 |   34/  32 |   28/  26 |   40/  36 |
 |-----------+-----------+-----------+-----------+-----------|
 |   35 1113 |  136/ 138 |   30/  24 |   26/  22 |   32/  26 |
 |-----------+-----------+-----------+-----------+-----------|
 |   95  160 |   54/  54 |   32/  30 |   26/  24 |   36/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   95 1000 |   92/  94 |   32/  28 |   26/  24 |   36/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 |   95 1111 |  100/ 102 |   42/  42 |   34/  32 |   54/  50 |
 |-----------+-----------+-----------+-----------+-----------|
 |   95 1113 |  128/ 128 |   38/  36 |   34/  30 |   44/  40 |
 |-----------+-----------+-----------+-----------+-----------|
 |  160 1000 |   50/  52 |   26/  22 |   24/  20 |   28/  22 |
 |-----------+-----------+-----------+-----------+-----------|
 |  160 1111 |  120/ 120 |   38/  34 |   30/  30 |   44/  40 |
 |-----------+-----------+-----------+-----------+-----------|
 |  160 1113 |  144/ 148 |   36/  34 |   28/  28 |   42/  36 |
 |-----------+-----------+-----------+-----------+-----------|
 | 1000 1111 |  475/ 475 |   28/  36 |   44/  24 |   32/  32 |
 |-----------+-----------+-----------+-----------+-----------|
 | 1000 1113 |  104/ 104 |   44/  42 |   36/  34 |   54/  48 |
 |-----------+-----------+-----------+-----------+-----------|
 | 1111 1113 | 2303/2143 |   32/  26 |   26/  28 |   34/  26 |
 |-----------------------------------------------------------|

KBK is curious how they compare on 'bad' inputs such as 2178309 and 1346269 - the smallest pair that require 32 trips through the loop, by the way.

See also: Least common multiple


The title of this page should have been greatest common divisor. There is no such thing as a greatest common denominator. There is a least common denominator, which is the least common multiple of the denominators.


AMG, being silly, presents GCD's cousin, Least Common Divisor. It only works right for positive arguments, but oh well.

 proc lcd {args} {return 1}

RVB GCD, LCM for list of 1 or more integers

 proc gcd_lcm {args} {
   set arg1 [lindex $args 0]
   set args [lrange $args 1 end]
   set gcd $arg1
   set lcm 1
   foreach arg $args {
     set lcm [expr {$lcm * $arg}]
     while {$arg != 0} {
       set t $arg
       set arg [expr {$gcd % $arg}]
       set gcd $t
     }
   }
   set lcm [expr ($arg1/$gcd)*$lcm]
   return [list $gcd $lcm]
 }


Category Mathematics - Category Algorithm - Arts and crafts of Tcl-Tk programming