isqrt von neumann continued, fully proved in 16 or 32 bits

but not yet in 64
the specs are awful, need cleaning

examples/in_progress/isqrt_von_neumann.mlw

 
-module VonNeumann
-
 (*
 
 algorithm from Warren's "Hacker's Delight" Figure 11.4
@@ -23,6 +21,98 @@ int isqrt4(unsigned x) {
 
 *)
 
+
+
+module VonNeumann16
+
+  use import ref.Ref
+  use import bv.BV16
+
+  function sqr (x:t) : t = mul x x
+  function succ (x:t) : t = add x (1:t)
+  function pred (x:t) : t = sub x (1:t)
+  function pow2 (n:t) : t = lsl_bv (1:t) n
+  predicate is_pow2 (x:t) (n:t) = bw_and x (pred (pow2 n)) = (0:t)
+
+  lemma sqr_add: forall x y.
+    sqr (add x y) = add (sqr x) (add (sqr y) (mul (2:t) (mul x y)))
+
+  lemma sqr_aux: forall x. pred (sqr x) = mul (succ x) (pred x)
+  
+  let isqrt16 (x:t) : t
+    ensures { ule (sqr result) x }
+    ensures { result = (0xff:t) \/ ult x (sqr (succ result))  }
+  = let num = ref x in
+    let bits = ref (0x4000:t) in
+    let res = ref (0:t) in
+    let ghost m = ref (8:t) in
+    let ghost bits_g = ref (0x80:t) in   (* 2^{m-1} *)
+    let ghost res_g = ref (0:t) in     (* res / 2^m *)
+    while !bits <> (0:t) do
+      variant { t'int !bits }
+      (* 0 <= m <= 8 *)
+      invariant { ule !m (8:t) }
+      (* bits_g = 2^{m-1} ou 0 si m=0 *)
+      invariant { !bits_g = if !m=(0:t) then (0:t) else pow2 (pred !m) }
+      (* bits = bits_g^2 *)
+      invariant { !bits = sqr !bits_g }
+      (* res_g multiple de 2^m *)
+      invariant { is_pow2 !res_g !m }
+      (* res_g smaller than 2^8 *)
+      invariant { ult !res_g (0x100:t) }
+      (* res = res_g * 2^m *)
+      invariant { !res = mul !res_g (pow2 !m) }
+      (* num <= x *)
+      invariant { ule !num x }
+      (* (x - num) est le carré de (res / 2^m) *)
+      invariant { sub x !num = sqr !res_g }
+      (* x < (res_g + 2^m)^2 *)
+      invariant { ule x (mul (succ (add !res_g (pow2 !m))) (pred (add !res_g (pow2 !m)))) }
+
+      (* m is not null, let's subtract 1 *)
+      assert { !m <> (0:t) };
+      'L1:
+      ghost m := pred !m;
+      assert { (at !m 'L1) = succ !m };
+      assert { !res = mul !res_g (pow2 (succ !m)) };
+      assert { !bits_g = pow2 !m };
+      let b = bw_or !res !bits in
+      assert { b = add !res !bits };
+      assert { b = mul (add (mul (2:t) !res_g) !bits_g) (pow2 !m) };
+      res := lsr_bv !res (1:t);
+      assert { !res = mul !res_g (pow2 !m) };
+      if uge !num b then
+        begin
+        num := sub !num b;
+        'L:
+        res := bw_or !res !bits;
+        assert { !res = add (at !res 'L) !bits };
+        assert { !res = mul (add !res_g !bits_g) (pow2 !m) };
+        res_g := add !res_g !bits_g
+        end;
+      bits := lsr_bv !bits (2:t);
+      ghost bits_g := lsr_bv !bits_g (1:t)
+    done;
+    assert { !m = (0:t) };
+    assert { !bits = (0:t) };
+    assert { !res = !res_g };
+    assert { pred (sqr (add !res_g (pow2 !m))) = 
+             mul (succ (add !res_g (pow2 !m))) (pred (add !res_g (pow2 !m))) };
+    assert { ule x (pred (sqr (add !res_g (pow2 !m)))) };
+    assert { !res_g <> (0xff:t) -> sqr (add !res_g (pow2 !m)) <> (0:t) };
+    assert { !res_g <> (0xff:t) -> ult x (sqr (add !res_g (pow2 !m))) };
+    assert { sqr (succ !res) = add (sqr !res_g) (succ (mul (2:t) !res_g)) };
+    !res
+
+
+end
+
+
+
+
+
+module VonNeumann32
+
   use import ref.Ref
   use import bv.BV32
 
@@ -35,7 +125,7 @@ int isqrt4(unsigned x) {
   lemma sqr_add: forall x y.
     sqr (add x y) = add (sqr x) (add (sqr y) (mul (2:t) (mul x y)))
 
-  let isqrt4 (x:t) : t
+  let isqrt32 (x:t) : t
     ensures { ule (sqr result) x }
     (* ensures { ult x (sqr (succ result))  } *)
   = let num = ref x in
@@ -62,6 +152,8 @@ int isqrt4(unsigned x) {
       invariant { ule !num x }
       (* (x - num) est le carré de (res / 2^m) *)
       invariant { sub x !num = sqr !res_g }
+      (* x < (res_g + 2^m)^2 *)
+      (* invariant { ule x (mul (succ (add !res_g (pow2 !m))) (pred (add !res_g (pow2 !m)))) } *)
 
       (* m is not null, let's subtract 1 *)
       assert { !m <> (0:t) };
@@ -89,20 +181,95 @@ int isqrt4(unsigned x) {
     done;
     assert { !m = (0:t) };
     assert { !bits = (0:t) };
+    assert { !res = !res_g };
+    (* assert { sqr (succ !res) = add (sqr !res_g) (succ (mul (2:t) !res_g)) }; *)
     !res
 
 
 end
 
-(*
-(*
-let ghost i = ref (0:t) in
-    let ghost bits = ref (30:t) in
-*)
-invariant { !bits = sub (30:t) (add !i !i) }
-      invariant { !i = (16:t) /\ !m = (0:t)
-                  \/
-                  ult !i (16:t) /\ !m = lsl_bv (1:t) !bits }
-      invariant { let ly = lsr_bv !y !bits in
-                  ule (mul ly ly) x0 }
-*)
+
+
+
+
+
+
+
+
+
+
+
+
+
+module VonNeumann64
+
+  use import ref.Ref
+  use import bv.BV64
+
+  function sqr (x:t) : t = mul x x
+  function succ (x:t) : t = add x (1:t)
+  function pred (x:t) : t = sub x (1:t)
+  function pow2 (n:t) : t = lsl_bv (1:t) n
+  predicate is_pow2 (x:t) (n:t) = bw_and x (pred (pow2 n)) = (0:t)
+
+  lemma sqr_add: forall x y.
+    sqr (add x y) = add (sqr x) (add (sqr y) (mul (2:t) (mul x y)))
+
+  let isqrt64 (x:t) : t
+    ensures { ule (sqr result) x }
+    (* ensures { ult x (sqr (succ result))  } *)
+  = let num = ref x in
+    let bits = ref (0x4000_0000_0000_0000:t) in
+    let res = ref (0:t) in
+    let ghost m = ref (32:t) in
+    let ghost bits_g = ref (0x8000_0000:t) in   (* 2^{m-1} *)
+    let ghost res_g = ref (0:t) in     (* res / 2^m *)
+    while !bits <> (0:t) do
+      variant { t'int !bits }
+      (* 0 <= m <= 32 *)
+      invariant { ule !m (32:t) }
+      (* bits_g = 2^{m-1} ou 0 si m=0 *)
+      invariant { !bits_g = if !m=(0:t) then (0:t) else pow2 (pred !m) }
+      (* bits = bits_g^2 *)
+      invariant { !bits = sqr !bits_g }
+      (* res_g multiple de 2^m *)
+      invariant { is_pow2 !res_g !m }
+      (* res_g smaller than 2^32 *)
+      invariant { ult !res_g (0x1_0000_0000:t) }
+      (* res = res_g * 2^m *)
+      invariant { !res = mul !res_g (pow2 !m) }
+      (* num <= x *)
+      invariant { ule !num x }
+      (* (x - num) est le carré de (res / 2^m) *)
+      invariant { sub x !num = sqr !res_g }
+
+      (* m is not null, let's subtract 1 *)
+      assert { !m <> (0:t) };
+      'L1:
+      ghost m := pred !m;
+      assert { (at !m 'L1) = succ !m };
+      assert { !res = mul !res_g (pow2 (succ !m)) };
+      assert { !bits_g = pow2 !m };
+      let b = bw_or !res !bits in
+      assert { b = add !res !bits };
+      assert { b = mul (add (mul (2:t) !res_g) !bits_g) (pow2 !m) };
+      res := lsr_bv !res (1:t);
+      assert { !res = mul !res_g (pow2 !m) };
+      if uge !num b then
+        begin
+        num := sub !num b;
+        'L:
+        res := bw_or !res !bits;
+        assert { !res = add (at !res 'L) !bits };
+        assert { !res = mul (add !res_g !bits_g) (pow2 !m) };
+        res_g := add !res_g !bits_g
+        end;
+      bits := lsr_bv !bits (2:t);
+      ghost bits_g := lsr_bv !bits_g (1:t)
+    done;
+    assert { !m = (0:t) };
+    assert { !bits = (0:t) };
+    !res
+
+
+end