Integer square root on machine integers
An efficient computation of the square root of machine integers, using bitwise operatons.
Authors: Claude Marché / Sylvain Dailler
Topics: Bitwise operations
Tools: 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) 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 endVersion 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 endVersion 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
download ZIP archive