Initial revision

diff --git a/numbertheory.txt b/numbertheory.txt
new file mode 100644 (file)
index 0000000..23db6de
--- /dev/null
+++ b/numbertheory.txt
@@ -0,0 +1,99 @@
+// Elementary Algorithms in Number Theory
+
+val N: ℕ;
+val M: ℕ;
+val K: ℕ;
+
+type Base   = ℕ[N];
+type Exp    = ℕ[M];
+type Result = ℕ[N^M];
+type nat    = ℕ[K];
+
+////////////////////////////////////////////////////////////////////
+// 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
+   ensures result = x^n;
+{
+   var res:Result = 1; 
+   var locn:Exp = n;
+   var locx:ℕ[N^(2*M)] = x;
+   while locn > 0 do
+     invariant n ≥ locn ∧ locn ≥ 0 ∧ x^n = res*(locx^locn);
+     decreases locn;
+   {
+     if locn%2 = 1 then
+     {
+       res = locx*res;
+     }
+     locx = locx^2;
+     locn = locn/2;
+   }
+   return res;
+}
+
+////////////////////////////////////////////////////////////////////
+// Algorithm: Calcluate x^n by the classic Algorithm
+////////////////////////////////////////////////////////////////////
+proc xpowern(x:Base,n:Exp): Result
+  ensures result = x^n;
+{
+   var res:Result = 1;
+   var i:Nat[M+1] = 0;
+   while i < n do
+     invariant 0 ≤ i ∧ i ≤ n ∧ res = x^i;
+     decreases n-i+1;
+   {
+     res = res * x;
+     i = i+1;
+   }
+   return res;
+}
+
+proc main(): ()
+{
+   var s:Base = 0;
+   var t:Exp  = 0;
+   
+   while s < N do
+   {
+         s=s+1;
+         while t<N do
+         {
+            t = t+1;
+            print s,t,xpowern(s,t);
+         }
+   }
+   
+}
+
+
+
+
+
+
+
+
+