From bb690892f5d30fbd124b3e49e22c9fda4186a3cd Mon Sep 17 00:00:00 2001
From: Christoph Fuerst <ch.fuerst@gmx.at>
Date: Fri, 14 Apr 2017 19:23:50 +0200
Subject: [PATCH] initial revision of Steins algorithm

---
 stein.txt | 71 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 71 insertions(+)
 create mode 100644 stein.txt

diff --git a/stein.txt b/stein.txt
new file mode 100644
index 0000000..b7bda62
--- /dev/null
+++ b/stein.txt
@@ -0,0 +1,71 @@
+
+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;
+}
+
-- 
2.1.4