diff --git a/math32.tal b/math32.tal
index 67f4037..9789b25 100644
--- a/math32.tal
+++ b/math32.tal
@@ -56,8 +56,9 @@
 %DEBUG { #ff #0e DEO }
 %RTN  { JMP2r }
 %TOR { ROT ROT } ( a b c -> c a b )
-( %TOR2 { ROT2 ROT2 } )
 %COMPLEMENT32 { SWP2 #ffff EOR2 SWP2 #ffff EOR2 }
+%DUP4 { OVR2 OVR2 }
+%POP4 { POP2 POP2 }
 
 ( bitcount: number of bits needed to represent number )
 ( equivalent to floor[log2[x]] + 1 )
@@ -411,3 +412,23 @@
     ,&cur0 LDR2 ,&cur1 LDR2 #01 ;rshift32 JSR2 ,&cur1 STR2 ,&cur0 STR2 ( cur >>= 1 )
     ,&cur0 LDR2 ,&cur1 LDR2 ;non-zero32 JSR2 ,&loop JCN ( if cur>0, loop. else we're done )
     RTN
+
+( greatest common divisor - euclidean algorithm )
+@gcd32 ( x** y** -> z** )
+    &loop ( x y )
+    DUP4 ( x y y )
+    ;is-zero32 JSR2 ( x y y=0? )
+    ,&done JCN ( x y )
+    DUP4 ( x y y )
+    STH2 STH2 ( x y [y] )
+    ;mod32 JSR2 ( r=x%y [y] )
+    STH2r ( rhi rlo yhi [ylo] )
+    ROT2 ( rlo yhi rhi [ylo] )
+    ROT2 ( yhi rhi rlo [ylo] )
+    STH2r ( yhi rhi rlo ylo )
+    ROT2 ( yhi rlo ylo rhi )
+    ROT2 ( yhi ylo rhi rlo )
+    ,&loop JMP
+    &done
+    POP4 ( x )
+    RTN
diff --git a/test-math32.tal b/test-math32.tal
index 8a1d952..8bdd3d4 100644
--- a/test-math32.tal
+++ b/test-math32.tal
@@ -53,6 +53,7 @@ RTN
     ;buf LDA LIT '- EQU ;test-sub32 JCN2
     ;buf LDA LIT '/ EQU ;test-div32 JCN2
     ;buf LDA LIT '% EQU ;test-mod32 JCN2
+    ;buf LDA LIT 'G EQU ;test-gcd32 JCN2
     ;buf LDA LIT 'L EQU ;test-lshift32 JCN2
     ;buf LDA LIT 'R EQU ;test-rshift32 JCN2
     ;buf LDA LIT 'B EQU ;test-bitcount32 JCN2
@@ -125,6 +126,7 @@ RTN
 @test-sub32 ;sub32 ;binary-32-test JMP2
 @test-div32 ;div32 ;binary-32-test JMP2
 @test-mod32 ;mod32 ;binary-32-test JMP2
+@test-gcd32 ;gcd32 ;binary-32-test JMP2
 @test-lshift32 ;lshift32 ;binary-32-8-32-test JMP2
 @test-rshift32 ;rshift32 ;binary-32-8-32-test JMP2
 @test-bitcount32 ;bitcount32 ;unary-32-8-test JMP2
diff --git a/tester.py b/tester.py
index e0bcaa6..ff116e3 100644
--- a/tester.py
+++ b/tester.py
@@ -54,6 +54,9 @@ def pipe():
 
 def bitcount(x):
     return floor(log(x, 2)) + 1
+
+def gcd(a, b):
+    return a if b == 0 else gcd(b, a % b)
         
 def main():
     trials = 100
@@ -64,6 +67,7 @@ def main():
     test(p, trials, b'*', [('x', u32), ('y', u32)], u32, lambda x, y: x * y)
     test(p, trials, b'/', [('x', u32), ('y', u32)], u32, lambda x, y: x // y)
     test(p, trials, b'%', [('x', u32), ('y', u32)], u32, lambda x, y: x % y)
+    test(p, trials, b'G', [('x', u32), ('y', u32)], u32, gcd)
     test(p, trials, b'L', [('x', u32), ('y', u5)], u32, lambda x, y: x << y)
     test(p, trials, b'R', [('x', u32), ('y', u5)], u32, lambda x, y: x >> y)
     test(p, trials, b'B', [('x', u32)], u8, bitcount)