isqrt.mlw 2.09 KB
Newer Older
1

2
(** {1 Integer square root} *)
3

4
module Square
MARCHE Claude's avatar
isqrt  
MARCHE Claude committed
5

6 7
  use import int.Int

Andrei Paskevich's avatar
Andrei Paskevich committed
8
  function sqr (x:int) : int = x * x
MARCHE Claude's avatar
isqrt  
MARCHE Claude committed
9

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

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

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

18 19 20 21 22 23 24 25 26 27 28 29
  predicate isqrt_spec (x res:int) =
    res >= 0 /\ sqr res <= x < sqr (res + 1)
end

(** {2 Simple algorithm} *)

module Simple

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

30 31
  let isqrt (x:int) : int
    requires { x >= 0 }
32
    ensures { isqrt_spec x result }
33
  = let count = ref 0 in
34 35
    let sum = ref 1 in
    while !sum <= x do
MARCHE Claude's avatar
MARCHE Claude committed
36 37 38
      invariant { !count >= 0 }
      invariant { x >= sqr !count }
      invariant { !sum = sqr (!count+1) }
39
      variant   { x - !count }
40 41 42 43 44
      count := !count + 1;
      sum := !sum + 2 * !count + 1
    done;
    !count

45 46 47
  let main ()
    ensures { result = 4 }
  = isqrt 17
MARCHE Claude's avatar
isqrt  
MARCHE Claude committed
48

49
end
MARCHE Claude's avatar
isqrt  
MARCHE Claude committed
50

51
(** {2 Another algorithm, in the style of Newton-Raphson} *)
52 53 54 55 56

module NewtonMethod

  use import int.Int
  use import int.ComputerDivision
57
  use import ref.Ref
58
  use import Square
59

60 61 62 63
  let sqrt (x : int) : int
    requires { x >= 0 }
    ensures  { isqrt_spec x result }
  = if x = 0 then 0 else
64 65
    if x <= 3 then 1 else
    let y = ref x in
66
    let z = ref (div (1 + x) 2) in
67 68
    while !z < !y do
      variant { !y }
69 70 71
      invariant { !z > 0 }
      invariant { !y > 0 }
      invariant { !z = div (div x !y + !y) 2 }
72 73
      invariant { x < sqr (!y + 1) }
      invariant { x < sqr (!z + 1) }
74
      y := !z;
75 76
      z := div (div x !z + !z) 2;
      (* A few hints to prove preservation of the last invariant *)
MARCHE Claude's avatar
MARCHE Claude committed
77 78 79 80 81 82
      assert { x < sqr (!z + 1)
        by let a = div x !y in
              x < a * !y + !y
           so a + !y <= 2 * !z + 1
           so sqr (a + !y + 1) <= sqr (2 * !z + 2)
           so 4 * (sqr (!z + 1) - x)
83 84 85 86 87
             = 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 }
88
    done;
MARCHE Claude's avatar
MARCHE Claude committed
89 90
    assert { !y * !y <= div x !y * !y
             by !y <= div x !y };
91 92 93
    !y

end