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 DKFLars 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)}
} 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}]
} 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.tcl.tk/pub/tcl/mirror/ftp.procplace.com/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 [L4 ], 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
I think below is a better implementation, and that there are probably optimizations that could be done.
factor_num is pulled from https://wiki.tcl-lang.org/951
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]
}proc lcm { {num ""} {f 2} {i ""} } {
global cf
set num1 ""
set sum 0
set cf(fact) $i
foreach {s} $num { incr sum $s }
if { $sum > [llength $num] } {
foreach {n} $num {
if { [expr int($n/$f)-($n*1.0/$f)] == 0.0 } {
lappend num1 [expr ($n/$f)]
} else {
lappend num1 $n
}
}
if { $num == $num1 } {
if { $f == 2 } {
incr f
} else {
incr f 2
}
lappend cf(fact) 1
} else {
lappend cf(fact) $f
}
set num $num1
mcm $num $f $cf(fact)
}
set mcm 1
foreach {p} $cf(fact) { set mcm [expr $mcm*$p] }
return $mcm
}% lcm "4 8 5 3" 120