From: Christoph Fuerst <ch.fuerst@gmx.at>
Date: Tue, 11 Apr 2017 19:59:06 +0000 (+0200)
Subject: Added sieve of erathostenes
X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=072cd08324c3348ae465e8c5f35cb5d9ca64e05a;p=cfuerst%2Fformal-numbers.git

Added sieve of erathostenes
---

diff --git a/sieveEratosthenes.txt b/sieveEratosthenes.txt
index 77b5d46..4790088 100644
--- a/sieveEratosthenes.txt
+++ b/sieveEratosthenes.txt
@@ -1,33 +1,53 @@
-
 val MAX: ℕ;
 
 type nat = ℕ[MAX];
 
 pred divides(m:nat,n:nat) ⇔ ∃p:nat. m⋅p = n;
-pred isPrime(m:nat) ⇔ 1 < m ∧ ∀n:nat. divides(n,m) ⇒ (n≤1 ∨ n=m);
+pred isPrime(m:nat) ⇔ (m = 0) ∨ (m=1) ∨
+                      (1 < m ∧ ∀n:nat. divides(n,m) ⇒ (n≤1 ∨ n=m));
 
-proc SieveOfErathostenes(m:ℕ[MAX]): Array[MAX,nat]
-   requires m>1;
-   ensures ∀p:ℕ[MAX]. result[p]=1 ⇒ isPrime(p+1);
+
+proc IntegerRoot(a:nat): nat
+   ensures result^2 ≤ a ∧ a < (result+1)^2;
 {
-   var ret:Array[MAX,nat] := Array[MAX,nat](0);
-   var i:ℕ[MAX];
-   var j:ℕ[2*MAX];   
-   
-   ret[0] = 1;
-   ret[1] = 1;
-   
-//   for i=2;i<N;i=i+1 do
+   var x:nat = 0;
+   var y:ℕ[(MAX+1)^2] = 1;
+   var z:ℕ[(MAX+1)^2] = 1;
+
+   while y ≤ a do
+      invariant x^2 ≤ a ∧ y = x^2+z ∧ z = 2*x+1;
    {
-//      if ret[i-1] = 1 then
-      {
-         ret[i-1] = 1;
-         for j=(i-1)^2;j<MAX;j=j+i do
-            ret[j] = 1;
-      }
-      print ret;
+      x=x+1;
+      z=z+2;
+      y=y+z;
    }
 
+   return x;
+}
+
+
+proc SieveOfErathostenes(m:nat): Array[MAX,Bool]
+   requires m>1;
+   ensures ∀p:nat. ((p < m ∧ result[p]=true) ⇒ isPrime(p));
+{
+   var ret:Array[MAX,Bool] := Array[MAX,Bool](true);
+   var i:nat;
+   var k:nat = 0;
+   var j:nat;
+
+   for i=2; i≤IntegerRoot(m); i=i+1 do
+      decreases IntegerRoot(m)+1-i;
+   {
+	  if ret[i] = true then
+	  {
+             j = i^2;
+	     for k=0;i^2+k*i<MAX;k=k+1 do
+	     {
+                j=i^2+k*i;
+	        ret[j] = false;
+	     }
+	  }
+   }
    return ret;
 }
 
@@ -36,11 +56,11 @@ proc SetMin(s:Set[ℕ[MAX]]): ℕ[MAX]
    ensures ∀e ∈ s. result ≤ e;
 {
    var ret:ℕ[MAX] := MAX;
-   
+
    for x in s do
      if x < ret then
         ret = x;
-   
+
    return ret;
 }
 
@@ -53,17 +73,12 @@ proc SieveOfErathostenesSet(m:ℕ[MAX]): Set[ℕ[MAX]]
    var locmi:ℕ[MAX];
 
    while start ≠ ∅[ℕ[MAX]] do
+      invariant ∀e ∈ ret. isPrime(e);
    {
       locmi = SetMin(start);
 	  ret = ret ∪ {locmi};
 	  start = start \ ({locmi} ∪ { if divides(locmi,x) then x else locmi | x ∈ start});
-   }  
+   }
    return ret;
 }
 
-proc main(): ()
-{
-   SieveOfErathostenesSet(MAX);
-}
-
-