Added trial division algorithm for primes

diff --git a/numbertheory.txt b/numbertheory.txt
index 9cd4d16..f83a239 100644 (file)
--- a/numbertheory.txt
+++ b/numbertheory.txt
@@ -12,6 +12,56 @@ type nat    = ℕ[K];
 // Elementary division
 pred divides(m:nat,n:nat) ⇔ ∃p:nat. m*p = n;
 
+////////////////////////////////////////////////////////////////////
+// Algorithm: Calcluate Integer Root of a
+// See http://www.inf.fh-flensburg.de/lang/se/veri/aufg/wurzel.htm
+////////////////////////////////////////////////////////////////////
+proc integerroot(a:nat): nat
+   ensures result^2 ≤ a ∧ a < (result+1)^2;
+{
+   var x:nat = 0;
+   var y:ℕ[(K+1)^2] = 1;
+   var z:nat = 1;
+
+   while y ≤ a do
+      invariant x^2 ≤ a ∧ y = x^2+z ∧ z = 2*x+1;
+   {
+      x=x+1;
+      z=z+2;
+      y=y+z;
+   }
+
+   return x;
+}
+
+////////////////////////////////////////////////////////////////////
+// prime numbers
+// trial division from 3 ... sqrt(p) + 1
+////////////////////////////////////////////////////////////////////
+proc isprime(p:nat): Bool
+{
+   var res:Bool = true;
+   var i:nat = 2;
+
+   if p ≥ 2 then
+   {
+      while i < integerroot(p)+1 ∧ res = true do
+         invariant 2 ≤ i ∧ i ≤ p;
+         decreases p-i;
+      {
+         res = ¬ divides(i,p);
+         i = i+1;
+      }
+   }
+   else
+   {
+      res = false;
+   }
+
+   return res;
+}
+
+
 // Quadratic residue
 // a ∈ ℤ/pℤ is a quadratic residue ⇔
 // x^2 - a has a zero in ℤ/pℤ
@@ -38,28 +88,6 @@ proc isquadraticresidue(a:nat,p:nat): Bool
 
 
 ////////////////////////////////////////////////////////////////////
-// Algorithm: Calcluate Integer Root of a
-// See http://www.inf.fh-flensburg.de/lang/se/veri/aufg/wurzel.htm
-////////////////////////////////////////////////////////////////////
-proc integerroot(a:nat): nat
-   ensures result^2 ≤ a ∧ a < (result+1)^2;
-{
-   var x:nat = 0;
-   var y:Nat[(K+1)^2] = 1;
-   var z:nat = 1;
-   
-   while y ≤ a do
-      invariant x^2 ≤ a ∧ y = x^2+z ∧ z = 2*x+1;
-   {
-      x=x+1;
-      z=z+2;
-      y=y+z;
-   }
-   
-   return x;
-}
-
-////////////////////////////////////////////////////////////////////
 // Algorithm: Calcluate x^n by Right-To-Left Exponentation
 ////////////////////////////////////////////////////////////////////
 proc r2lexponentation(x:Base,n:Exp): Result