isqrt.mlw 1.97 KB
Newer Older
1

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

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

6
  use int.Int
7

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
  predicate isqrt_spec (x res:int) =
    res >= 0 /\ sqr res <= x < sqr (res + 1)
end

(** {2 Simple algorithm} *)

module Simple

26
  use int.Int
27
  use ref.Refint
28
  use Square
29

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

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

module NewtonMethod

55 56 57 58
  use int.Int
  use mach.int.Int
  use ref.Ref
  use 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
    if x <= 3 then 1 else
65 66 67 68 69 70 71 72 73 74 75
    let ref y = x in
    let ref z = (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 <- (x / z + z) / 2;
76
      (* A few hints to prove preservation of the last invariant *)
77 78 79 80 81 82 83 84 85 86
      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)
             = sqr (2 * z + 2) - 4 * x
             >= sqr (a + y + 1) - 4 * x
             > sqr (a + y + 1) - 4 * (a * y + y)
             = sqr (a + 1 - y)
87
             >= 0 }
88
    done;
89 90 91
    assert { y * y <= div x y * y
             by y <= div x y };
    y
92 93

end