examples/isqrt.mlw

   1
   2 (** {1 Integer square root} *)
   3
   4 module Square
   5
   6   use int.Int
   7
   8   function sqr (x:int) : int = x * x
   9
  10   lemma sqr_non_neg: forall x:int. sqr x >= 0
  11
  12   lemma sqr_increasing:
  13     forall x y:int. 0 <= x <= y -> sqr x <= sqr y
  14
  15   lemma sqr_sum :
  16     forall x y : int. sqr(x+y) = sqr x + 2*x*y + sqr y
  17
  18   predicate isqrt_spec (x res:int) =
  19     res >= 0 /\ sqr res <= x < sqr (res + 1)
  20 end
  21
  22 (** {2 Simple algorithm} *)
  23
  24 module Simple
  25
  26   use int.Int
  27   use ref.Refint
  28   use Square
  29
  30   let isqrt (x:int) : int
  31     requires { x >= 0 }
  32     ensures  { isqrt_spec x result }
  33   = let ref count = 0 in
  34     let ref sum = 1 in
  35     while sum <= x do
  36       invariant { count >= 0 }
  37       invariant { x >= sqr count }
  38       invariant { sum = sqr (count+1) }
  39       variant   { x - count }
  40       count += 1;
  41       sum += 2 * count + 1
  42     done;
  43     count
  44
  45   let main ()
  46     ensures { result = 4 }
  47   = isqrt 17
  48
  49 end
  50
  51 (** {2 Another algorithm, in the style of Newton-Raphson} *)
  52
  53 module NewtonMethod
  54
  55   use int.Int
  56   use mach.int.Int
  57   use ref.Ref
  58   use Square
  59
  60   let sqrt (x : int) : int
  61     requires { x >= 0 }
  62     ensures  { isqrt_spec x result }
  63   = if x = 0 then 0 else
  64     if x <= 3 then 1 else
  65     let ref y = x in
  66     let ref z = (1 + x) / 2 in
  67     while z < y do
  68       variant { y }
  69       invariant { z > 0 }
  70       invariant { y > 0 }
  71       invariant { z = div (div x y + y) 2 }
  72       invariant { x < sqr (y + 1) }
  73       invariant { x < sqr (z + 1) }
  74       y <- z;
  75       z <- (x / z + z) / 2;
  76       (* A few hints to prove preservation of the last invariant *)
  77       assert { x < sqr (z + 1)
  78         by let a = div x y in
  79               x < a * y + y
  80            so a + y <= 2 * z + 1
  81            so sqr (a + y + 1) <= sqr (2 * z + 2)
  82            so 4 * (sqr (z + 1) - x)
  83              = sqr (2 * z + 2) - 4 * x
  84              >= sqr (a + y + 1) - 4 * x
  85              > sqr (a + y + 1) - 4 * (a * y + y)
  86              = sqr (a + 1 - y)
  87              >= 0 }
  88     done;
  89     assert { y * y <= div x y * y
  90              by y <= div x y };
  91     y
  92
  93 end