examples_in_progress/gmp_square_root.mlw

   1
   2 (*
   3   An efficient algorithm proposed by P. Zimmermann in 1999 to compute
   4   square roots of arbitrarily large integers
   5
   6   Following the Coq proof described in
   7
   8     A Proof of GMP Square Root
   9     Yves Bertot, Nicolas Magaud, and Paul Zimmermann
  10     Journal of Automated Reasoning 29: 225–252, 2002.
  11
  12   and available at
  13
  14     http://www-sop.inria.fr/lemme/AOC/SQRT/index.html
  15 *)
  16
  17 module GmpModel
  18
  19   use int.Int
  20   use array.Array
  21
  22   constant beta' : int
  23   axiom one_lt_beta' : 1 < beta'
  24   constant beta : int = 2 * beta'
  25
  26   type memory = array int
  27
  28   val m: memory
  29     (** invariant: each cell contains an integer in [\[0..beta\[] *)
  30
  31   type pointer = int
  32   type size = int
  33
  34   predicate valid (m: memory) (p: pointer) (s: size) =
  35     0 <= p /\ p + s <= length m
  36
  37   predicate no_overlap (p1: pointer) (s1: size) (p2: pointer) (s2: size) =
  38     p1 + s1 <= p2 \/ p2 + s2 <= p1
  39
  40   function i (m: memory) (p: pointer) (s: size) : int
  41     (** the GMP integer in [m\[p..p+s\[] *)
  42
  43   predicate modifies (m1 m2: memory) (p: pointer) (s: size) =
  44     forall q: pointer. (q < p \/ p + s <= q) -> m2[q] = m1[q]
  45     (** nothing modified apart from [\[p..p+s\[] *)
  46
  47   axiom frame:
  48     forall m1 m2: memory, p: pointer, s: size.
  49     modifies m1 m2 p s ->
  50     forall a: pointer, l: size. no_overlap p s a l -> i m2 a l = i m1 a l
  51
  52   predicate modifies2 (m1 m2: memory)
  53                       (p1: pointer) (s1: size) (p2: pointer) (s2: size) =
  54     forall q: pointer. (   (q < p1 \/ p1 + s1 <= q)
  55                         /\ (q < p2 \/ p2 + s2 <= q)) -> m2[q] = m1[q]
  56     (** nothing modified apart from [\[p1..p1+s1\[] and [\[p2..p2+s2\[] *)
  57
  58   axiom frame2:
  59     forall m1 m2: memory, p1 p2: pointer, s1 s2: size.
  60     modifies2 m1 m2 p1 s1 p2 s2 ->
  61     forall a: pointer, l: size. no_overlap p1 s1 a l -> no_overlap p2 s2 a l ->
  62     i m2 a l = i m1 a l
  63
  64   predicate is_normalized (n h: int) = n < h <= 4 * n
  65
  66   function sqr (z: int) : int = z*z
  67
  68   let int_of_bool (b: bool) : int = if b then 1 else 0
  69
  70   function mult_bool (c: bool) (n: int) : int = if c then n else 0
  71
  72 end
  73
  74 module GmpAuxiliaryfunctions
  75
  76   use export GmpModel
  77   use int.Int
  78   use int.ComputerDivision
  79   use int.Power
  80   use ref.Ref
  81   use array.Array
  82
  83   val rb: ref bool
  84
  85   val mpn_sub_n (r a b: pointer) (l: size) : unit
  86     requires { valid m r l }
  87     requires { valid m a l }
  88     requires { valid m b l }
  89     requires { no_overlap r l a l }
  90     requires { no_overlap r l b l }
  91     requires { no_overlap a l b l }
  92     writes   { m, rb }
  93     ensures  { i m r l = (mult_bool !rb (power beta l) + i (old m) a l)
  94                          - i (old m) b l }
  95     ensures  { modifies (old m) m r l }
  96
  97   val mpn_divrem (q a: pointer) (la: size) (b: pointer) (lb: size) : unit
  98     requires { valid m a la }
  99     requires { valid m b lb }
 100     requires { valid m q (la - lb) }
 101     requires { lb <= la }
 102     requires { power beta lb <= 2 * i m b lb }
 103     requires { no_overlap a la b lb }
 104     requires { no_overlap q (la - lb) a la }
 105     requires { no_overlap q (la - lb) b lb }
 106     writes   { m, rb }
 107     ensures  { i (old m) a la =
 108                  (mult_bool !rb (power beta (la - lb)) + i m q (la - lb)) *
 109                  (i (old m) b lb) + i m a lb }
 110     ensures  { i m a lb < i (old m) b lb }
 111     ensures  { modifies2 (old m) m q (la - lb) a lb }
 112
 113   val mpn_sqr_n (r a: pointer) (l: size) : unit
 114     requires { valid m a l }
 115     requires { valid m r (2*l) }
 116     requires { no_overlap r (2*l) a l }
 117     writes   { m }
 118     ensures  { i m r (2*l) = sqr (i (old m) a l) }
 119     ensures  { modifies (old m) m r (2*l) }
 120
 121   val mpn_sqrtrem2 (s r n: pointer) : unit
 122     requires { valid m n 2 }
 123     requires { valid m s 1 }
 124     requires { valid m r 1 }
 125     requires { no_overlap r 1 s 1 }
 126     requires { is_normalized (i m n 2) (power beta 2) }
 127     writes   { m, rb }
 128     ensures  { let s = i m s 1 in
 129                let r = i m r 1 + mult_bool !rb beta in
 130                i (old m) n 2 = s*s + r /\ 0 <= r <= 2*s }
 131     ensures  { modifies (old m) m n 2 }
 132
 133   val rz: ref int
 134
 135   val single_and_1 (a: pointer) : unit
 136     requires { valid m a 1 }
 137     writes   { rz }
 138     ensures  { !rz = mod m[a] 2 }
 139
 140   val mpn_rshift2 (a b: pointer) (l: size) : unit
 141     requires { valid m a l }
 142     requires { valid m b l }
 143     requires { b <= a } (* ? *)
 144     requires { no_overlap a l b l }
 145     writes   { m }
 146     ensures  { i (old m) b l = 2 * i m a l + mod (old m)[b] 2 }
 147     ensures  { modifies (old m) m a l }
 148
 149
 150 (*
 151   let square_s_and_sub (np sp nn l h q c b: int)
 152   = mpn_sqr_n (np + nn) sp l;
 153     mpn_sub_n np np (np + nn) (2 * l);
 154     ----
 155     b := !q + (bool2Z !rb);
 156   if (nat_eq_bool !l !h) then
 157      c := !c - !b
 158   else
 159      begin
 160        (mpn_sub_1 (plus np (mult (S (S O)) !l))
 161                   (plus np (mult (S (S O)) !l))
 162                   (S O) !b);
 163        c := !c - (bool2Z !rb)
 164      end
 165 *)
 166
 167 end
 168
 169 module GmpSquareRoot
 170
 171   use GmpAuxiliaryfunctions
 172   use ref.Ref
 173
 174 (**
 175   let division_step (np sp: pointer) (nn l h: size) (c q: ref int) : unit
 176   = if !q <> 0 then mpn_sub_n (np + 2*l) (np + 2*l) (sp + l) h;
 177     mpn_divrem sp (np + l) nn (sp + l + h);
 178     q := !q + int_of_bool !rb;
 179     single_and_1 sp;
 180     let c = !rz in
 181     mpn_rshift2 sp sp l;
 182     m[sp + l-1] <- l_or (s + l-1) (shift_1 !q);
 183     q := shift_2 !q;
 184     if !c <> 0 then begin
 185       mpn_add_n (np + l) (np + l) (sp + l) h;
 186       c := int_of_bool !rb
 187     end
 188
 189   let rec mpn_dq_sqrtrem (sp np: pointer) (nn: size) : unit
 190     requires { 0 < nn }
 191     variant  { nn }
 192   = if nn = 1 then
 193       mpn_sqrtrem2 sp np np
 194     else begin
 195       let l = div nn 2 in
 196       let h = nn - l in
 197       mpn_dq_sqrtrem (sp + l) (np + 2*l) h;
 198       let q = if !rb then 1 else 0 in
 199       let c = ref 0 in
 200       division_step np sp nn l h c q;
 201       square_s_and_sub np sp nn l h q c b;
 202       mpn_add_1 (sp + l) (sp + l) h q;
 203       let q = int_of_bool !rb in
 204       correct_result np sp nn l h q c b;
 205       rb := !c > 0
 206     end
 207 **)
 208
 209 end