Wiki Agenda Contact Version française

Integer square root

Simple implementation in O(sqrt(n)).


Authors: Claude Marché

Topics: Non-linear Arithmetic / Tricky termination

Tools: Why3

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



Integer square root

module Square

  use import int.Int

  function sqr (x:int) : int = x * x

  lemma sqr_non_neg: forall x:int. sqr x >= 0

  lemma sqr_increasing:
    forall x y:int. 0 <= x <= y -> sqr x <= sqr y

  lemma sqr_sum :
    forall x y : int. sqr(x+y) = sqr x + 2*x*y + sqr y

  predicate isqrt_spec (x res:int) =
    res >= 0 /\ sqr res <= x < sqr (res + 1)
end

Simple algorithm

module Simple

  use import int.Int
  use import ref.Ref
  use import Square

  let isqrt (x:int) : int
    requires { x >= 0 }
    ensures { isqrt_spec x result }
  = let count = ref 0 in
    let sum = ref 1 in
    while !sum <= x do
      invariant { !count >= 0 /\ x >= sqr !count /\ !sum = sqr (!count+1) }
      variant   { x - !count }
      count := !count + 1;
      sum := !sum + 2 * !count + 1
    done;
    !count

  let main ()
    ensures { result = 4 }
  = isqrt 17

end

Another algorithm, in the style of Newton-Raphson

module NewtonMethod

  use import int.Int
  use import int.ComputerDivision
  use import ref.Ref
  use import Square

  let sqrt (x : int) : int
    requires { x >= 0 }
    ensures  { isqrt_spec x result }
  = if x = 0 then 0 else
    if x <= 3 then 1 else
    let y = ref x in
    let z = ref (div (1 + x) 2) in
    while !z < !y do
      variant { !y }
      invariant { !z > 0 }
      invariant { !y > 0 }
      invariant { !z = div (div x !y + !y) 2 }
      invariant { x < sqr (!y + 1) }
      invariant { x < sqr (!z + 1) }
      y := !z;
      z := div (div x !z + !z) 2;
      (* A few hints to prove preservation of the last invariant *)
      let ghost a = div x !y in
      assert { x < a * !y + !y };
      assert { a + !y <= 2 * !z + 1 };
      assert { 0 <= a + !y + 1 };
      assert { sqr (a + !y + 1) <= sqr (2 * !z + 2) };
      assert { 4 * (sqr (!z + 1) - x)
             = sqr (2 * !z + 2) - 4 * x
             >= sqr (a + !y + 1) - 4 * x
             > sqr (a + !y + 1) - 4 * (a * !y + !y)
             = sqr (a + 1 - !y)
             >= 0 }
    done;
    assert { !y * !y <= div x !y * !y by !y <= div x !y };
    !y

end

download ZIP archive