Integer Square Root: Hacker's version

We present certified version of programs for computing the square root of machine integers. This implementation is a very efficient one making use of low-level bitwise operations. The algorithm originated from Warren's book "Hacker's Delight", Figure 11.4.

The reference C code is as follows

int isqrt4(unsigned x) {
     unsigned m, y, b;
     m = 0x40000000;
     y = 0;
     while(m != 0) {              // Do 16 times.
        b = y | m;
        y = y >> 1;
        if (x >= b) {
           x = x - b;
           y = y | m;
        }
        m = m >> 2;
     }
     return y;
  }

We present 3 different implementations respectively for 16, 32 and 64
bits. The proofs become significantly more difficult when increasing
the number of bits, so a need for more manual steps shows up in the 64
bits version.
Version 16 bits
module VonNeumann16

  use ref.Ref
  use 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_add2: forall x y.
    sqr (add x y) = add (sqr x) (mul y (add (mul (2:t) x) y))

  let isqrt16 (x:t) : t
    ensures { ule (sqr result) x }
    ensures { ule x (pred (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 not (eq !bits (0:t)) do
      variant { !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) is the square of (res / 2^m) *)
      invariant { sub x !num = sqr !res_g }
      invariant {ule (add !res_g (pow2 !m)) (0x100:t)}
      (* x < (res_g + 2^m)^2 *)
      invariant { ule x (pred (sqr (add !res_g (pow2 !m)))) }

      (* m is not null, let's subtract 1 *)
      assert { !m <> (0:t) };
      label L1 in
      ghost m := pred !m;
      assert { (!m at 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;
        label L in
        res := bw_or !res !bits;
        assert { !res = add (!res at 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;
    !res

end


Version 32 bits
module VonNeumann32

  use ref.Ref
  use bv.BV32

  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_add2: forall x y.
    sqr (add x y) = add (sqr x) (mul y (add (mul (2:t) x) y))

  let isqrt32 (x:t) : t
    ensures { ule (sqr result) x }
    ensures { ule x (pred (sqr (succ result)))  }
  = let num = ref x in
    let bits = ref (0x4000_0000:t) in
    let res = ref (0:t) in
    let ghost m = ref (16:t) in
    let ghost bits_g = ref (0x8000:t) in   (* 2^{m-1} *)
    let ghost res_g = ref (0:t) in     (* res / 2^m *)
    while not (eq !bits (0:t)) do
      variant { !bits }
      (* 0 <= m <= 16 *)
      invariant { ule !m (16: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^16 *)
      invariant { ult !res_g (0x1_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 }
      invariant {ule (add !res_g (pow2 !m)) (0x1_0000:t)}
      (* x < (res_g + 2^m)^2 *)
      invariant { ule x (pred (sqr (add !res_g (pow2 !m)))) }

      (* m is not null, let's subtract 1 *)
      assert { !m <> (0:t) };
      label L1 in
      ghost m := pred !m;
      assert { (!m at 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;
        label L in
        res := bw_or !res !bits;
        assert { !res = add (!res at 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;
    !res

end


Version 64 bits
module VonNeumann64

  use ref.Ref
  use 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_add2: forall x y.
    sqr (add x y) = add (sqr x) (mul y (add (mul (2:t) x) y))

  let isqrt64 (x:t) : t
    ensures { ule (sqr result) x }
    ensures { ule x (pred (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 not (eq !bits (0:t)) do
      variant { !bits }
      (* 0 <= m <= 32 *)
      invariant { ule !m (32:t) }
      (* bits_g = 2^{m-1} or 0 if 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 of 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) is the square of (res / 2^m) *)
      invariant { sub x !num = sqr !res_g }
      invariant {ule (add !res_g (pow2 !m)) (0x1_0000_0000:t)}
      (* x < (res_g + 2^m)^2 *)
      invariant { ule x (pred (sqr (add !res_g (pow2 !m)))) }

      (* m is not null, let's subtract 1 *)
      assert { !m <> (0:t) };
      label L1 in
      ghost m := pred !m;
      assert { (!m at 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;
        label L in
        res := bw_or !res !bits;
        assert { !res = add (!res at 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