From: Christoph Fuerst <ch.fuerst@gmx.at>
Date: Tue, 21 Mar 2017 18:14:01 +0000 (+0100)
Subject: Initial revision
X-Git-Url: http://git.risc.jku.at/gitweb/?a=commitdiff_plain;h=9629c0039208a791d1cd11955e4b4fee8075e6bd;p=cfuerst%2Fformal-numbers.git

Initial revision
---

9629c0039208a791d1cd11955e4b4fee8075e6bd
diff --git a/numbertheory.txt b/numbertheory.txt
new file mode 100644
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);
+   	  }
+   }
+   
+}
+
+
+
+
+
+
+
+
+