--- /dev/null
+
+val N:ℕ;
+
+type int = ℤ[-N,N];
+type bigint = ℤ[-2*N,2*N];
+
+pred divides(m:int,n:int) ⇔ ∃p:int. m*p = n;
+
+fun gcd(m:int,n:int): int
+ requires m ≠ 0 ∨ n ≠ 0;
+= choose result:int with
+ divides(result,m) ∧ divides(result,n) ∧
+ ¬∃r:int. divides(r,m) ∧ divides(r,n) ∧ r > result;
+
+
+proc GCDStein(a:int,b:int): int
+ requires a≠0 ∨ b≠0;
+ ensures result = gcd(a,b);
+{
+ var res:int = 0;
+ var k:int = 0;
+ var loca:int = a;
+ var locb:int = b;
+ var t:bigint;
+
+ if loca=0 then
+ {
+ res = locb;
+ }
+ else
+ {
+ while loca%2 = 0 ∧ locb%2 = 0 do
+ {
+ loca = loca/2;
+ locb = locb/2;
+ k = k+1;
+ }
+ if loca%2 = 1 then
+ {
+ t = -b;
+ }
+ else
+ {
+ t = a;
+ }
+ while t ≠ 0 do
+ {
+ while t%2 = 0 do
+ {
+ t = t/2;
+ }
+ if t > 0 then
+ {
+ loca = t;
+ }
+ else
+ {
+ locb = -t;
+ }
+ t = loca-locb;
+ print t;
+ }
+ res = loca*2^k;
+ }
+
+ if res < 0 then
+ res = res*(-1);
+
+ return res;
+}
+