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'' [http://aleph0.clarku.edu/~djoyce/java/elements/elements.html] (Book 7, Propositions 1 and 2 [http://aleph0.clarku.edu/~djoyce/java/elements/bookVII/propVII2.html]). 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[http://theory.lcs.mit.edu/~clr/].) (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)} } ====== [glennj] ====== namespace path {::tcl::mathop ::tcl::mathfunc} proc gcdI {p q} { while {$q != 0} {lassign "$q [% $p $q]" p q} abs $p } ====== [DKF]: The fastest I've found so far is this one (which minimizes assignments and data movement): ====== proc gcdJ {p q} { while 1 { if {![set p [expr {$p % $q}]]} {return [expr {$q>0?$q:-$q}]} if {![set q [expr {$q % $p}]]} {return [expr {$p>0?$p:-$p}]} } } ====== Pure bytecode, minimized data movement, no memory allocation in loop. ---- ''[Cameron Laird] wrote in [comp.lang.tcl]:'' 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. [glennj] referring to the perl implementation at [http://www.rosettacode.org/wiki/Greatest_common_divisor#Iterative_binary_algorithm_3], this implementation is about twice as fast on my machine: ====== package require Tcl 8.5 namespace path {::tcl::mathop ::tcl::mathfunc} proc gcd_bin {p q} { if {$p == $q} {return [abs $p]} set p [abs $p] if {$q == 0} {return $p} set q [abs $q] if {$p < $q} {lassign [list $q $p] p q} set k 1 while {($p & 1) == 0 && ($q & 1) == 0} { set p [>> $p 1] set q [>> $q 1] set k [<< $k 1] } set t [expr {$p & 1 ? -$q : $p}] while {$t} { while {$t & 1 == 0} {set t [>> $t 1]} if {$t > 0} {set p $t} {set q [- $t]} set t [- $p $q] } return [* $p $k] } ====== Timing: set u 40902 set v 24140 time {gcd_bin $v $u} 100000 ;# ==> 9.54744 microseconds per iteration time {binary-gcd $v $u} 100000 ;# ==> 20.07662 microseconds per iteration ---- 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] } ====== [jnxn] The gcd_lcm function appears to have a bug when dealing with a list longer then 2 numbers. % gcd_lcm 4 8 5 3 1 480 % lcm 4 8 5 3 120 % ===== proc factor_num {num} { for {set i 2} {$i <= $num} {} { if {![expr {$num % $i}]} { lappend factors $i set num [expr {$num / $i}] } else { incr i } } return $factors } proc lcm {args} { set args_factor_list {} set unique_factor_list {} set unique_factor_list_count {} foreach arg_value $args { set factor_list [factor_num $arg_value] lappend args_factor_list $factor_list foreach factor $factor_list { if {[lsearch -all -exact -inline $unique_factor_list $factor] ne $factor} { lappend unique_factor_list $factor lappend unique_factor_list_count 0 } set factor_count [llength [lsearch -all -exact $factor_list $factor]] set unique_factor_index [lsearch -all -exact $unique_factor_list $factor] if {[lindex $unique_factor_list_count $unique_factor_index] <= $factor_count} { lset unique_factor_list_count $unique_factor_index $factor_count } } } set lcmexpr "int(" foreach unique_factor $unique_factor_list { set unique_factor_index [lsearch -all -exact $unique_factor_list $unique_factor] if {$unique_factor_index == [expr [llength $unique_factor_list] -1]} { append lcmexpr "pow($unique_factor,[lindex $unique_factor_list_count $unique_factor_index]))" } else { append lcmexpr "pow($unique_factor,[lindex $unique_factor_list_count $unique_factor_index])* " } } return [expr $lcmexpr] } ==== <> Mathematics | Algorithm | Arts and crafts of Tcl-Tk programming