## Integer square root on machine integers

An efficient computation of the square root of machine integers, using bitwise operatons.

Auteurs: Claude Marché / Sylvain Dailler

Catégories: Bitwise operations

Outils: Why3

see also the index (by topic, by tool, by reference, by year)

# 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)

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) };
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)

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) };
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)

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) };
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

```