mp2.mlw 289 KB
Newer Older
MARCHE Claude's avatar
MARCHE Claude committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

(** {2 An Efficient Library for Arbitrary-Precision Integer Library}

    This code implements arbitrary-precision integers in the style of GMP.
    Supported operations are addition, multiplication, division on non-negative integers,
    and also comparisons and shifts.

    For detailed info, see the paper ``How to Get an Efficient yet
Verified Arbitrary-Precision Integer Library'' by Raphaël Rieu-Helft,
Claude Marché and Guillaume Melquiond, published in VSTTE'2017
proceedings, Springer LNCS.

   Important note: this Why3 code currently requires the branch `mp' from the
   {h <a href="https://gforge.inria.fr/scm/?group_id=2990">Why3 development git repository</a>}

MARCHE Claude's avatar
MARCHE Claude committed
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

  {h <h3> Instructions for extraction to C and runnings tests</h2>}

  Execute the {h <code> make </code>} command to extract and build the C file {h <code> N.c </code>} in the {h <code> build </code>} subdirectory.

  C compilation requires some GMP headers (the {h <code> longlong.h </code>} and {h <code> config.h </code>} headers are used), so an installation of GMP with sources is required to compile it (see {h <a href="https://gmplib.org#DOWNLOAD">here</a>}).

  The following commands test the library on random inputs and compare the results to GMP for verification:

  {h <pre>
    export GMP_DIR=/path/to/your/gmp/install
    make tests
  </pre>}

  To execute the benchmarks, an installation of GMP with assembly
  disabled is required.  First install GMP with assembly disabled in a
  fresh directory:

  {h <pre>
    cd /path/to/gmp/sources
    ./configure --disable-assembly --prefix=/where/to/install/gmp
    make
    make install
  </pre>}

  You are now ready to benchmark the library (gnuplot required to view the plots). The process takes a few minutes.
  {h <pre>
    export GMP_DIR=/path/to/gmp/sources
    export GMP_LIB=/where/to/install/gmp/lib
    export LD_LIBRARY_PATH=$GMP_LIB
    make plots
    </pre>}

MARCHE Claude's avatar
MARCHE Claude committed
49 50
*)

MARCHE Claude's avatar
MARCHE Claude committed
51 52
(** {2 The main Why3 module for non-negative integers} *)

53 54
module N

55 56
  use import array.Array
  use import map.Map
57 58
  use map.MapEq
  use map.Const
59
  use import int.Int
60

MARCHE Claude's avatar
MARCHE Claude committed
61 62
  (** {3 complements to map standard library} *)

63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
  predicate map_eq_sub_shift (x y:map int 'a) (xi yi sz:int) =
    forall i. 0 <= i < sz -> x[xi+i] = y[yi+i]

  let lemma map_eq_shift (x y:map int 'a) (xi yi sz k:int)
    requires { map_eq_sub_shift x y xi yi sz }
    requires { 0 <= k < sz }
    ensures { x[xi+k] = y[yi+k] }
  = ()

  let rec lemma map_eq_shift_zero (x y: map int 'a) (n m: int)
    requires { map_eq_sub_shift x y n n (m-n) }
    variant { m - n }
    ensures { MapEq.map_eq_sub x y n m }
  =
    if n < m then
    begin
      assert { forall i. 0 <= i < m-n -> x[n+i] = y[n+i] };
      assert { forall i. n <= i < m ->
81 82
                 let j = i - n in 0 <= j < m-n ->
                     x[n+j] = y[n+j] -> x[i] = y[i]};
83 84
      map_eq_shift_zero x y (n+1) m;
    end
85
    else ()
86

87 88
  use import mach.int.Int32
  use import ref.Ref
89 90 91 92
  use import mach.int.UInt64 as Limb
  use import int.Int
  use import int.Power

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
93 94
  meta compute_max_steps 0x100000

MARCHE Claude's avatar
MARCHE Claude committed
95 96
  (** {3 Long integers as arrays of libs} *)

97 98 99 100
  type limb = uint64

  lemma limb_max_bound: 1 <= max_uint64

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
101
  function l2i (x:limb) : int = Limb.to_int x
102

103
  function p2i (i:int32) : int = int32'int i
104

Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
105
  exception Break
Raphaël Rieu-Helft's avatar
Cleanup  
Raphaël Rieu-Helft committed
106
  exception Return32 int32
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
107
  exception ReturnLimb limb
108 109 110 111 112 113

  let lemma prod_compat_strict_r (a b c:int)
    requires { 0 <= a < b }
    requires { 0 < c }
    ensures { c * a < c * b }
  = ()
114 115 116 117 118
  let lemma prod_compat_r (a b c:int)
    requires { 0 <= a <= b }
    requires { 0 <= c }
    ensures { c * a <= c * b }
  = ()
119

MARCHE Claude's avatar
MARCHE Claude committed
120
  (** {3 Integer value of a natural number} *)
Raphaël Rieu-Helft's avatar
Cleanup  
Raphaël Rieu-Helft committed
121

122 123 124 125 126 127 128 129 130
  (** [value_sub x n m] denotes the integer represented by
     the digits x[n..m-1] with lsb at index n *)
  let rec ghost function value_sub (x:map int limb) (n:int) (m:int) : int
     variant {m - n}
   =
     if n < m then
       l2i x[n] + radix * value_sub x (n+1) m
       else 0

131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
  let rec lemma value_sub_frame (x y:map int limb) (n m:int)
    requires { MapEq.map_eq_sub x y n m }
    variant  { m - n }
    ensures  { value_sub x n m = value_sub y n m }
  =
    if n < m then value_sub_frame x y (n+1) m else ()

  let rec lemma value_sub_frame_shift (x y:map int limb) (xi yi sz:int)
    requires { map_eq_sub_shift x y xi yi sz }
    variant { sz }
    ensures { value_sub x xi (xi+sz) = value_sub y yi (yi+sz) }
 =
    if sz>0
    then begin
      map_eq_shift x y xi yi sz 0;
      assert { forall i. 0 <= i < sz-1 ->
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
147
                 let j = 1+i in x[xi+j] = y[yi+j] };
148 149 150 151 152 153 154 155 156
      value_sub_frame_shift x y (xi+1) (yi+1) (sz-1)
      end
    else assert { 1+2 = 3 }

  let rec lemma value_sub_tail (x:map int limb) (n m:int)
    requires { n <= m }
    variant  { m - n }
    ensures  {
      value_sub x n (m+1) =
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
157
        value_sub x n m + (Map.get x m) * power radix (m-n) }
158
  = [@vc:sp] if n < m then value_sub_tail x (n+1) m else ()(*assert { 1+2=3 }*)
159 160 161 162 163 164

  let rec lemma value_sub_concat (x:map int limb) (n m l:int)
    requires { n <= m <= l}
    variant  { m - n }
    ensures  {
      value_sub x n l =
165
        value_sub x n m + value_sub x m l * power radix (m-n) }
166 167 168 169 170 171 172 173
  =
  if n < m then
     begin
     assert {n<m};
     value_sub_concat x (n+1) m l
     end
  else ()

174 175
  let lemma value_sub_head (x:map int limb) (n m:int)
    requires { n < m }
176
    ensures { value_sub x n m = x[n] + radix * value_sub x (n+1) m }
177 178
  = value_sub_concat x n (n+1) m

179 180 181 182
  let lemma value_sub_update (x:map int limb) (i n m:int) (v:limb)
    requires { n <= i < m }
    ensures {
      value_sub (Map.set x i v) n m =
183
      value_sub x n m + power radix (i - n) * (v -(Map.get x i))
184 185 186 187
    }
  = assert { MapEq.map_eq_sub x (Map.set x i v) n i };
    assert { MapEq.map_eq_sub x (Map.set x i v) (i+1) m };
    value_sub_concat x n i m;
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
188 189 190
    value_sub_concat (Map.set x i v) n i m;
    value_sub_head x i m;
    value_sub_head (Map.set x i v) i m
191 192

  let rec lemma value_zero (x:map int limb) (n m:int)
193
    requires { MapEq.map_eq_sub x (Const.const Limb.zero_unsigned) n m }
194 195 196 197
    variant  { m - n }
    ensures  { value_sub x n m = 0 }
  = if n < m then value_zero x (n+1) m else ()

198 199 200 201 202 203 204 205 206 207
  let lemma value_sub_update_no_change (x: map int limb) (i n m: int) (v:limb)
     requires { n <= m }
     requires { i < n \/ m <= i }
     ensures { value_sub x n m = value_sub (Map.set x i v) n m }
  = value_sub_frame x (Map.set x i v) n m

  let lemma value_sub_shift_no_change (x:map int limb) (ofs i sz:int) (v:limb)
     requires { i < 0 \/ sz <= i }
     requires { 0 <= sz }
     ensures { value_sub x ofs (ofs + sz) =
208
               value_sub (Map.set x (ofs+i) v) ofs (ofs+sz) }
209 210
  = value_sub_frame_shift x (Map.set x (ofs+i) v) ofs ofs sz

MARCHE Claude's avatar
MARCHE Claude committed
211
  (** {3 Comparisons} *)
212

213
  let rec lemma value_sub_lower_bound (x:map int limb) (x1 x2:int)
214 215 216 217
    variant  { x2 - x1 }
    ensures  { 0 <= value_sub x x1 x2 }
  = if x2 <= x1 then () else
      begin
218
        value_sub_head x x1 x2;
219
        value_sub_lower_bound x (x1+1) x2
220 221 222 223 224 225 226 227
      end

  let rec lemma value_sub_upper_bound (x:map int limb) (x1 x2:int)
    requires { x1 <= x2 }
    variant  { x2 - x1 }
    ensures  { value_sub x x1 x2 < power radix (x2 - x1) }
  = if x1 = x2 then () else
      begin
228
      value_sub_tail x x1 (x2-1);
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
229
      assert { value_sub x x1 x2
230
               <= value_sub x x1 (x2-1) + power radix (x2-x1-1) * (radix - 1) };
231 232 233 234 235
      value_sub_upper_bound x x1 (x2-1)
      end

  let lemma value_sub_lower_bound_tight (x:map int limb) (x1 x2:int)
    requires { x1 < x2 }
236
    ensures  { power radix (x2-x1-1) *  l2i (Map.get x (x2-1)) <= value_sub x x1 x2 }
237
  = assert   { value_sub x x1 x2 = value_sub x x1 (x2-1)
238
               + power radix (x2-x1-1) * l2i (Map.get x (x2-1)) }
239 240 241 242 243 244

  let lemma value_sub_upper_bound_tight (x:map int limb) (x1 x2:int)
    requires { x1 < x2 }
    ensures  { value_sub x x1 x2 < power radix (x2-x1-1) *  (l2i (Map.get x (x2-1)) + 1) }
  = value_sub_upper_bound x x1 (x2-1)

245 246
  use import mach.c.C
  type t = ptr limb
247

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
248
  function value (x:t) (sz:int) : int =
249
     value_sub (pelts x) x.offset (x.offset + sz)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
250

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
251 252 253 254 255 256 257
  let lemma value_tail (x:t) (sz:int32)
    requires { 0 <= sz }
    ensures  { value x (sz+1) = value x sz + (pelts x)[x.offset + sz] * power radix sz }
  = value_sub_tail (pelts x) x.offset (x.offset + p2i sz)

  meta remove_prop axiom value_tail

258 259 260 261 262 263 264 265 266 267
  let lemma value_concat (x:t) (n m:int32)
    requires { 0 <= n <= m }
    ensures  { value x m
             = value x n + power radix n
                            * value_sub (pelts x) (x.offset + n) (x.offset + m) }

  = value_sub_concat (pelts x) x.offset (x.offset + p2i n) (x.offset + p2i m)

  meta remove_prop axiom value_concat

268 269 270
  function compare_int (x y:int) : int =
    if x < y then -1 else if x=y then 0 else 1

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
271 272 273 274 275 276 277 278
  let copy (r x:t) (sz:int32) : unit
    requires { valid x sz }
    requires { valid r sz }
    ensures { map_eq_sub_shift (pelts r) (pelts x) r.offset x.offset sz }
  =
    let zero = Int32.of_int 0 in
    let one = Int32.of_int 1 in
    let i = ref zero in
279 280
    let xp = ref (C.incr x zero) in
    let rp = ref (C.incr r zero) in
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
281 282 283 284 285 286 287 288 289 290 291
    while (Int32.(<) !i sz) do
      variant { p2i sz - p2i !i }
      invariant { 0 <= !i <= sz }
      invariant { map_eq_sub_shift (pelts r) (pelts x) r.offset x.offset !i }
      invariant { pelts !xp = pelts x }
      invariant { pelts !rp = pelts r }
      invariant { plength !xp = plength x }
      invariant { plength !rp = plength r }
      invariant { !xp.offset = x.offset + !i }
      invariant { !rp.offset = r.offset + !i }
      C.set !rp (C.get !xp);
292 293
      rp.contents <- C.incr !rp one;
      xp.contents <- C.incr !xp one;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
294 295 296 297
      i := Int32.(+) !i one;
    done


298 299
  (** [compare_same_size] compares [x[0..sz-1]] and [y[0..sz-1]] as unsigned integers. It corresponds to [GMPN_CMP]. *)
  let compare_same_size (x y:t) (sz:int32) : int32
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
300 301
    requires { valid x sz }
    requires { valid y sz }
302
    ensures { result = compare_int (value x sz) (value y sz) }
303 304 305 306 307
  =
   let i = ref sz in
   try
     while Int32.(>=) !i (Int32.of_int 1) do
       variant { p2i !i }
308 309
       invariant { 0 <= !i <= sz }
       invariant { forall j. !i <= j < sz ->
310
                   (pelts x)[x.offset+j] = (pelts y)[y.offset+j] }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
311
       assert { forall j. 0 <= j < sz - !i ->
312 313
                let k = !i+j in
                !i <= k < sz ->
314
                (pelts x)[x.offset+k] = (pelts y)[y.offset+k] /\
315 316 317
                (pelts x)[!i+x.offset+j] = (pelts y)[!i+y.offset+j] };
       value_sub_frame_shift (pelts x) (pelts y) (p2i !i+x.offset)
                             (p2i !i+y.offset) ((p2i sz) - (p2i !i));
318 319 320
       let ghost k = p2i !i in
       i := Int32.(-) !i (Int32.of_int 1);

321
       assert { 0 <= !i < sz };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
322 323
       let lx = get_ofs x !i in
       let ly = get_ofs y !i in
324
       if (not (Limb.(=) lx ly))
325
       then begin
326 327
            value_sub_concat (pelts x) x.offset (x.offset+k) (x.offset+p2i sz);
            value_sub_concat (pelts y) y.offset (y.offset+k) (y.offset+p2i sz);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
328 329 330
            assert { compare_int (value x sz)
                       (value y sz)
                   = compare_int (value x k) (value y k) };
331 332
            value_sub_tail (pelts x) x.offset (x.offset+k-1);
            value_sub_tail (pelts y) y.offset (y.offset+k-1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
333
            if Limb.(>) lx ly
334 335 336
            then begin
             value_sub_upper_bound (pelts y) y.offset (y.offset+k-1);
             value_sub_lower_bound (pelts x) x.offset (x.offset+k-1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
337
             assert { value x k - value y k =
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
338
                      (l2i lx - ly) * (power radix (k-1))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
339
                    - ((value y (k-1)) - (value x (k-1)))
340
                       };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
341
             assert { (lx - ly) * (power radix (k-1))
342
                      >= power radix (k-1)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
343
                      > ((value y (k-1)) - (value x (k-1)))
344
                       };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
345
             raise Return32 (Int32.of_int 1)
346 347
            end
            else begin
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
348
             assert { ly > lx };
349 350
             value_sub_upper_bound (pelts x) x.offset (x.offset+k-1);
             value_sub_lower_bound (pelts y) y.offset (y.offset+k-1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
351
             assert { value y k - value x k =
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
352
                    (ly - lx) * (power radix (k-1))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
353
                    - ((value x (k-1)) - (value y (k-1)))
354
                     };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
355
             assert { (ly - lx) * (power radix (k-1))
356
                      >= power radix (k-1)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
357
                      > ((value x (k-1)) - (value y (k-1)))
358
                     };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
359 360
            raise Return32 (Int32.(-_) (Int32.of_int 1))
            end
361
         end
362 363
       else ()
     done;
364
     value_sub_frame_shift (pelts x) (pelts y) x.offset y.offset (p2i sz);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
365
     Int32.of_int 0
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
366
   with Return32 r -> r
367 368 369 370
   end

   (* [is_zero] checks if [x[0..sz-1]] is zero. It corresponds to [mpn_zero_p]. *)
   let is_zero (x:t) (sz:int32) : int32
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
371
     requires { valid x sz }
372
     ensures { 0 <= Int32.to_int result <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
373
     ensures { Int32.to_int result = 1 <-> value x sz = 0 }
374 375
   =
     let i = ref sz in
376
     let uzero = Limb.of_int 0 in
377 378 379
     let lx = ref uzero in
     try
       while Int32.(>=) !i (Int32.of_int 1) do
380
         variant { p2i !i }
381 382
         invariant { 0 <= !i <= sz }
         invariant { value_sub (pelts x) (x.offset + !i) (x.offset + sz)=0 }
383 384
         let ghost k = p2i !i in
         i := Int32.(-) !i (Int32.of_int 1);
385
         assert { 0 <= !i < sz };
386
         lx := get_ofs x !i;
387
         if not (Limb.(=) !lx uzero)
388 389 390 391
         then begin
           value_sub_concat (pelts x) x.offset (x.offset+k) (x.offset + p2i sz);
           value_sub_lower_bound_tight (pelts x) x.offset (x.offset+k);
           value_sub_lower_bound (pelts x) (x.offset+k) (x.offset + p2i sz);
392
           raise Return32 (Int32.of_int 0);
393 394 395 396
         end
         else begin
           assert { 1+2=3 };
         end
397 398
       done;
       Int32.of_int 1
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
399
     with Return32 r -> r
400 401
     end

402 403
  (** [zero r sz] sets [(r,sz)] to zero. Corresponds to [mpn_zero]. *)
  let zero (r:t) (sz:int32) : unit
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
404
    requires { valid r sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
405
    ensures { value r sz = 0 }
406 407 408 409
  =
    let i = ref (Int32.of_int 0) in
    let lzero = Limb.of_int 0 in
    while Int32.(<) !i sz do
410 411
      invariant { 0 <= !i <= sz }
      variant { sz - !i }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
412
      invariant { value r !i = 0 }
413 414 415 416 417
      set_ofs r !i lzero;
      value_sub_tail (pelts r) r.offset (r.offset + p2i !i);
      i := Int32.(+) !i (Int32.of_int 1);
    done

MARCHE Claude's avatar
MARCHE Claude committed
418
  (** {3 Addition} *)
419

MARCHE Claude's avatar
MARCHE Claude committed
420
  (** [add_limb r x y sz] adds to [x] the value of the limb [y],
421
      writes the result in [r] and returns the carry. [r] and [x]
422
      have size [sz]. This corresponds to the function [mpn_add_1] *)
MARCHE Claude's avatar
MARCHE Claude committed
423
  (* r and x must be separated. This is enforced by Why3 regions in typing *)
424
  let add_limb (r x:t) (y:limb) (sz:int32) : limb
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
425 426
    requires { valid x sz }
    requires { valid r sz }
427
    requires { sz > 0 } (* ? GMP does the same for 0 and 1*)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
428 429
    ensures { value r sz + (power radix sz) * result =
              value x sz + y }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
430
    ensures { 0 <= result <= 1 }
431
    writes { r.data.elts }
432
  =
433
    let limb_zero = Limb.of_int 0 in
434 435 436
    let c = ref y in
    let lx = ref limb_zero in
    let i = ref (Int32.of_int 0) in
437
    while Int32.(<) !i sz && (not (Limb.(=) !c limb_zero)) do
438 439
      invariant { 0 <= !i <= sz }
      invariant { !i > 0 -> 0 <= !c <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
440 441
      invariant { value r !i + (power radix !i) * !c =
                  value x !i + y}
442
      variant { sz - !i }
443 444 445 446
      label StartLoop in
      lx := get_ofs x !i;
      let (res, carry) = add_with_carry !lx !c limb_zero in
      set_ofs r !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
447 448
      assert { value r !i + (power radix !i) * !c =
                  value x !i + y };
449
      c := carry;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
450 451 452 453 454
      value_tail r !i;
      value_tail x !i;
      assert { value r (!i+1) + (power radix (!i+1)) * !c
             = value x (!i+1) + y
             (* by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
455 456
             value r !i + (power radix !i) * !c
             = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
457
                               + (power radix !i) * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
458
             = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
459
                               + (power radix k) * radix * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
460 461
             = value r k + (power radix k) * (res + radix * !c)
             = value r k +
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
462
               (power radix k) * (!lx + (!c at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
463
             = value r k + (power radix k) * (!c at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
464
                               + (power radix k) * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
465
             = value x k + y + (power radix k) * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
466 467
             = value x !i + y*) };
      i := Int32.(+) !i (Int32.of_int 1);
468
    done;
469
    if Int32.(=) !i sz then !c
470 471
    else begin
    while Int32.(<) !i sz do
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
472
      invariant { !c  = 0 }
473
      invariant { 0 <= !i <= sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
474 475
      invariant { value r !i + (power radix !i) * !c =
                  value x !i + y}
476
      variant { sz - !i }
477 478
      lx := get_ofs x !i;
      set_ofs r !i !lx;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
479 480
      assert { value r !i + (power radix !i) * !c =
                  value x !i + y };
481 482
      let ghost k = p2i !i in
      i := Int32.(+) !i (Int32.of_int 1);
483 484
      value_sub_tail (pelts r) r.offset (r.offset + k);
      value_sub_tail (pelts x) x.offset (x.offset + k);
485 486 487 488 489
    done;
    !c
    end


490 491
  (** [add_limbs r x y sz] adds [x[0..sz-1]] and [y[0..sz-1]] and writes the result in [r].
      Returns the carry, either [0] or [1]. Corresponds to the function [mpn_add_n]. *)
492

493
  let add_limbs (r x y:t) (sz:int32) : limb
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
494 495 496
    requires { valid x sz }
    requires { valid y sz }
    requires { valid r sz }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
497
    ensures { 0 <= result <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
498 499
    ensures { value r sz + (power radix sz) * result =
            value x sz + value y sz }
500
    writes { r.data.elts }
501
    =
502
    let limb_zero = Limb.of_int 0 in
503 504 505 506 507
    let lx = ref limb_zero in
    let ly = ref limb_zero in
    let c = ref limb_zero in
    let i = ref (Int32.of_int 0) in
    while Int32.(<) !i sz do
508 509
      variant { sz - !i }
      invariant { 0 <= !i <= sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
510 511
      invariant { value r !i + (power radix !i) * !c =
                value x !i + value y !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
512
      invariant { 0 <= !c <= 1 }
513 514 515 516 517
      label StartLoop in
      lx := get_ofs x !i;
      ly := get_ofs y !i;
      let res, carry = add_with_carry !lx !ly !c in
      set_ofs r !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
518
      assert { value r !i + (power radix !i) * !c =
519 520
                value x !i + value y !i
               by value r !i = (value r !i at StartLoop) };
521
      c := carry;
522 523 524 525 526 527
      value_tail r !i;
      value_tail x !i;
      value_tail y !i;
      assert { value r (!i+1) + (power radix (!i+1)) * !c =
                value x (!i+1) + value y (!i+1)
              (*by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
528 529
              value r !i + (power radix !i) * !c
              = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
530
                   + (power radix !i) * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
531
              = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
532
                   + (power radix k) * radix * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
533 534
              = value r k + (power radix k) * (res + radix * !c)
              = value r k +
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
535
                  (power radix k) * (!lx + !ly + (!c at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
536
              = value r k + (power radix k) * (!c at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
537
                 + (power radix k) * (!lx + !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
538
              = value x k + value y k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
539
                 + (power radix k) * (!lx + !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
540 541 542 543 544 545
              = value x k + (power radix k) * !lx
                 + value y k + (power radix k) * !ly
              = value x !i
                 + value y k + (power radix k) * !ly
              = value x !i
                 + (value y k + (power radix k) * !ly)
546 547
              = value x !i + value y !i*) };
      i := Int32.(+) !i (Int32.of_int 1);
548 549 550
    done;
    !c

Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
551 552 553
  (** [add r x y sx sy] adds [(x, sx)] to [(y,sy)] and writes the
      result in [(r, sx)].  [sx] must be greater than or equal to
      [sy]. Returns carry, either 0 or 1. Corresponds to [mpn_add]. *)
554
  let add (r x y:t) (sx sy:int32) : limb
555
    requires { 0 <= sy <= sx }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
556 557 558
    requires { valid x sx }
    requires { valid y sy }
    requires { valid r sx }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
559 560
    ensures { value r sx + (power radix sx) * result =
              value x sx + value y sy }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
561
    ensures { 0 <= result <= 1 }
562
    writes { r.data.elts }
563
 =
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
564
    let limb_zero = Limb.of_int 0 in
565 566 567 568 569
    let lx = ref limb_zero in
    let ly = ref limb_zero in
    let c = ref limb_zero in
    let i = ref (Int32.of_int 0) in
    while Int32.(<) !i sy do
570 571
      variant { sy - !i }
      invariant { 0 <= !i <= sy }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
572 573
      invariant { value r !i + (power radix !i) * !c =
                value x !i + value y !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
574
      invariant { 0 <= !c <= 1 }
575 576 577 578 579
      label StartLoop in
      lx := get_ofs x !i;
      ly := get_ofs y !i;
      let res, carry = add_with_carry !lx !ly !c in
      set_ofs r !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
580 581
      assert { value r !i + (power radix !i) * !c =
                value x !i + value y !i };
582
      c := carry;
583 584 585 586 587 588
      value_tail r !i;
      value_tail x !i;
      value_tail y !i;
      assert { value r (!i+1) + (power radix (!i+1)) * !c =
                value x (!i+1) + value y (!i+1)
              (*by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
589 590
              value r !i + (power radix !i) * !c
              = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
591
                   + (power radix !i) * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
592
              = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
593
                   + (power radix k) * radix * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
594 595
              = value r k + (power radix k) * (res + radix * !c)
              = value r k +
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
596
                  (power radix k) * (!lx + !ly + (!c at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
597
              = value r k + (power radix k) * (!c at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
598
                 + (power radix k) * (!lx + !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
599
              = value x k + value y k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
600
                 + (power radix k) * (!lx + !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
601 602 603 604 605 606
              = value x k + (power radix k) * !lx
                 + value y k + (power radix k) * !ly
              = value x !i
                 + value y k + (power radix k) * !ly
              = value x !i
                 + (value y k + (power radix k) * !ly)
607 608
              = value x !i + value y !i*) };
      i := Int32.(+) !i (Int32.of_int 1);
609
    done;
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
610 611
    try
    begin while Int32.(<) !i sx do
612 613
      variant { sx - !i }
      invariant { sy <= !i <= sx }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
614 615
      invariant { value r !i + (power radix !i) * !c =
                value x !i + value y sy }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
616
      invariant { 0 <= !c <= 1 }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
617
      (if (Limb.(=) !c (Limb.of_int 0)) then raise Break);
618 619 620 621
      label StartLoop2 in
      lx := get_ofs x !i;
      let res, carry = add_with_carry !lx limb_zero !c in
      set_ofs r !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
622 623
      assert { value r !i + (power radix !i) * !c =
                value x !i + value y sy };
624
      c := carry;
625 626 627 628 629
      value_tail r !i;
      value_tail x !i;
      assert { value r (!i+1) + (power radix (!i+1)) * !c =
                value x (!i+1) + value y sy
              (*by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
630 631
              value r !i + (power radix !i) * !c
              = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
632
                   + (power radix !i) * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
633
              = value r k + (power radix k) * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
634
                   + (power radix k) * radix * !c
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
635 636
              = value r k + (power radix k) * (res + radix * !c)
              = value r k +
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
637
                  (power radix k) * (!lx + 0 + (!c at StartLoop2))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
638
              = value r k + (power radix k) * (!c at StartLoop2)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
639
                 + (power radix k) * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
640
              = value x k + value y sy
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
641
                 + (power radix k) * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
642
              = value x !i
643 644
                 + value y sy*) };
      i := Int32.(+) !i (Int32.of_int 1);
645
    done;
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
646 647 648 649 650 651 652 653
    assert { !i = sx }
    end
    with Break -> assert { !c = 0 }
    end;
    while Int32.(<) !i sx do
      variant { sx - !i }
      invariant { sy <= !i <= sx }
      invariant { !i = sx \/ !c = 0 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
654 655
      invariant { value r !i + power radix !i * !c =
                value x !i + value y sy }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
656 657 658
      assert { !c = 0 by !i < sx };
      lx := get_ofs x !i;
      set_ofs r !i !lx;
659 660 661 662 663 664
      value_tail r !i;
      value_tail x !i;
      assert { value r !i = value x !i + value y sy }; (* true with this, should not be needed *)
      assert { value r (!i+1) + power radix (!i+1) * !c
               = value x (!i+1) + value y sy
               (*
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
665
               by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
666 667 668 669 670 671 672
               value r !i + power radix !i * !c
                 = value r !i
                 = value r k + power radix k * !lx
               so value x !i
                  = value x k + power radix k * !lx
               so value r k
                  = value r k + power radix k * !c
673 674
                  = value x k + value y sy*) };
      i := Int32.(+) !i (Int32.of_int 1);
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
675
    done;
676 677
    !c

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706
  let add_in_place (x y:t) (sx sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid x sx }
    requires { valid y sy }
    ensures  { value x sx + (power radix sx) * result
               = value (old x) sx + value y sy }
    ensures  { 0 <= result <= 1 }
    ensures { forall j. j < x.offset \/ x.offset + sx <= j ->
              (pelts x)[j] = (pelts (old x))[j] }
    writes   { x.data.elts }
  =
    let ghost ox = { x } in
    let limb_zero = Limb.of_int 0 in
    let lx = ref limb_zero in
    let ly = ref limb_zero in
    let c = ref limb_zero in
    let i = ref (Int32.of_int 0) in
    while Int32.(<) !i sy do
      variant   { sy - !i }
      invariant { 0 <= !i <= sy }
      invariant { value x !i + (power radix !i) * !c =
                  value ox !i + value y !i }
      invariant { 0 <= !c <= 1 }
      invariant { forall j. !i <= j < sx ->
                  (pelts x)[x.offset + j] = (pelts ox)[x.offset + j] }
      invariant { forall j. j < x.offset \/ x.offset + sx <= j ->
                  (pelts x)[j] = (pelts (old x))[j] }
      label StartLoop in
      lx := get_ofs x !i;
707
      assert { !lx = (pelts ox)[ox.offset + !i] };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
708 709 710
      ly := get_ofs y !i;
      let res, carry = add_with_carry !lx !ly !c in
      set_ofs x !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
711 712 713 714 715 716
      assert { forall j. !i < j < sx ->
                 (pelts x)[x.offset + j]
                 = (pelts ox)[x.offset + j]
                 by (pelts x)[x.offset + j]
                 = (pelts (x at StartLoop))[x.offset + j]
                 = (pelts ox)[x.offset + j]};
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
717 718
      assert { value x !i + (power radix !i) * !c = value ox !i + value y !i };
      c := carry;
719 720 721 722 723 724
      value_tail x !i;
      value_tail ox !i;
      value_tail y !i;
      assert { value x (!i+1) + (power radix (!i+1)) * !c =
                value ox (!i+1) + value y (!i+1)
              (*by value ox k + (power radix k) * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743
                 = value ox !i
              so value x !i + (power radix !i) * !c
              = value x k + (power radix k) * res
                   + (power radix !i) * !c
              = value x k + (power radix k) * res
                   + (power radix k) * radix * !c
              = value x k + (power radix k) * (res + radix * !c)
              = value x k +
                  (power radix k) * (!lx + !ly + (!c at StartLoop))
              = value x k + (power radix k) * (!c at StartLoop)
                 + (power radix k) * (!lx + !ly)
              = value ox k + value y k
                 + (power radix k) * (!lx + !ly)
              = (value ox k + (power radix k) * !lx)
                 + (value y k + (power radix k) * !ly)
              = value ox !i
                 + (value y k + (power radix k) * !ly)
              = value ox !i
                 + (value y k + (power radix k) * !ly)
744 745
              = value ox !i + value y !i*) };
      i := Int32.(+) !i (Int32.of_int 1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
746 747 748 749 750 751 752 753 754 755 756 757 758 759 760
    done;
    try
    while Int32.(<) !i sx do
      variant   { sx - !i }
      invariant { sy <= !i <= sx }
      invariant { value x !i + (power radix !i) * !c =
                  value ox !i + value y sy }
      invariant { 0 <= !c <= 1 }
      invariant { forall j. !i <= j < sx ->
                  (pelts x)[x.offset + j] = (pelts ox) [x.offset + j] }
      invariant { forall j. j < x.offset \/ x.offset + sx <= j ->
                  (pelts x)[j] = (pelts (old x))[j] }
      (if (Limb.(=) !c limb_zero) then raise ReturnLimb limb_zero);
      label StartLoop2 in
      lx := get_ofs x !i;
761
      assert { !lx = (pelts ox)[ox.offset + !i] };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
762 763 764 765 766 767 768
      let res, carry = add_with_carry !lx limb_zero !c in
      value_sub_update_no_change (pelts x) (x.offset + p2i !i)
                                           (x.offset + p2i !i + 1)
                                           (x.offset + p2i sx) res;
      set_ofs x !i res;
      assert { value x !i + (power radix !i) * !c = value ox !i + value y sy };
      c := carry;
769
      assert { forall j. !i < j < sx ->
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
770
                  (pelts x)[x.offset + j] = (pelts ox) [x.offset + j] };
771 772 773 774 775
      value_tail ox !i;
      value_tail x !i;
      assert { value x (!i+1) + (power radix (!i+1)) * !c =
                value ox (!i+1) + value y sy
              (*by value ox k + (power radix k) * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
776 777 778 779 780 781 782 783 784 785 786 787 788 789 790
                 = value ox !i
              so
              value x !i + (power radix !i) * !c
              = value x k + (power radix k) * res
                   + (power radix !i) * !c
              = value x k + (power radix k) * res
                   + (power radix k) * radix * !c
              = value x k + (power radix k) * (res + radix * !c)
              = value x k +
                  (power radix k) * (!lx + 0 + (!c at StartLoop2))
              = value x k + (power radix k) * (!c at StartLoop2)
                 + (power radix k) * !lx
              = value ox k + value y sy
                 + (power radix k) * !lx
              = value ox !i
791 792
                 + value y sy*) };
      i := Int32.(+) !i (Int32.of_int 1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
793 794 795 796 797 798 799 800 801 802 803 804 805
    done;
    assert { !i = sx };
    !c
    with ReturnLimb n -> begin
      assert { n = 0 = !c };
      assert { forall j. x.offset + !i <= j < x.offset + sx
               -> (pelts x)[j] = (pelts ox)[j]
               by !i <= j - x.offset < sx
               so (pelts x)[x.offset + (j - x.offset)]
                  = (pelts ox)[x.offset + (j - x.offset)] };
      value_sub_frame (pelts x) (pelts ox) (x.offset + p2i !i) (x.offset + p2i sx);
      value_sub_concat (pelts x) x.offset (x.offset + p2i !i) (x.offset + p2i sx);
      value_sub_concat (pelts ox) x.offset (x.offset + p2i !i) (x.offset + p2i sx);
806
      assert { value x sx = value (old x) sx + value y sy };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
807 808 809 810
      n
      end
    end

811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870
  (** [incr x y sz] adds to [x] the value of the limb [y] in place.
      [x] has size [sz]. The addition must not overflow. This corresponds
      to [mpn_incr] *)
  let incr (x:t) (y:limb) (sz:int32) : unit
    requires { valid x sz }
    requires { sz > 0 }
    requires { value x sz + y < power radix sz }
    ensures  { value x sz = value (old x) sz + y }
    ensures { forall j. j < x.offset \/ x.offset + sz <= j ->
              (pelts x)[j] = (pelts (old x))[j] }
    writes   { x.data.elts }
  =
    let ghost ox = { x } in
    let c = ref y in
    let lx : ref limb = ref 0 in
    let i : ref int32 = ref 0 in
    while not (Limb.(=) !c 0) do
      invariant { 0 <= !i <= sz }
      invariant { !i = sz -> !c = 0 }
      invariant { !i > 0 -> 0 <= !c <= 1 }
      invariant { value x !i + (power radix !i) * !c
                  = value ox !i + y }
      invariant { forall j. !i <= j < sz ->
                  (pelts x)[x.offset + j] = (pelts ox)[x.offset + j] }
      invariant { forall j. j < x.offset \/ x.offset + sz <= j ->
                  (pelts x)[j] = (pelts ox)[j] }
      variant   { sz - !i }
      label StartLoop in
      lx := get_ofs x !i;
      assert { !lx = (pelts ox)[ox.offset + !i] };
      let (res, carry) = add_with_carry !lx !c 0 in (*TODO*)
      assert { res + radix * carry = !lx + !c }; (* TODO remove this *)
      value_sub_update_no_change (pelts x) (x.offset + p2i !i)
                                           (x.offset + p2i !i + 1)
                                           (x.offset + p2i sz) res;
      set_ofs x !i res;
      assert { forall j. !i < j < sz ->
                 (pelts x)[x.offset + j]
                 = (pelts ox)[x.offset + j] };
      assert { value x !i + (power radix !i) * !c = value ox !i + y };
      c := carry;
      value_tail x !i;
      value_tail ox !i;
      assert { value x (!i+1) + power radix (!i+1) * !c =
               value ox (!i+1) + y };
      i := Int32.(+) !i 1;
      assert { !i = sz -> !c = 0
               by value x sz + power radix sz * !c = value ox sz + y
                  so value ox sz + y < power radix sz
                  so 0 <= !c <= 1};
    done;
    value_concat x !i sz;
    value_concat ox !i sz;
    assert { forall j. x.offset + !i <= j < x.offset + sz ->
             (pelts x)[j] = (pelts ox)[j]
             by let k = j - x.offset in
                !i <= k < sz
                so (pelts x)[x.offset + k] = (pelts ox)[x.offset + k]};
    value_sub_frame (pelts x) (pelts ox) (x.offset + p2i !i) (x.offset + p2i sz)

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
871

Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
872 873 874 875
  (** [sub_limb r x y sz] substracts [y] from [(x, sz)] and writes
      the result to [(r, sz)]. Returns borrow, either 0 or
      1. Corresponds to [mpn_sub_1]. *)
  let sub_limb (r x:t) (y:limb) (sz:int32) : limb
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
876 877
    requires { valid x sz }
    requires { valid r sz }
878
    requires { 0 < sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
879 880
    ensures { value r sz - power radix sz * result
              = value x sz - y }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
881
    ensures { 0 <= result <= 1 }
882
    writes { r.data.elts }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
883 884 885 886 887
  =
    let limb_zero = Limb.of_int 0 in
    let b = ref y in
    let lx = ref limb_zero in
    let i = ref (Int32.of_int 0) in
888
    while Int32.(<) !i sz && (not (Limb.(=) !b limb_zero)) do
889 890
      invariant { 0 <= !i <= sz }
      invariant { !i > 0 -> 0 <= !b <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
891 892
      invariant { value r !i - power radix !i * !b
                  = value x !i - y }
893
      variant { sz - !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
894 895 896 897
      label StartLoop in
      lx := get_ofs x !i;
      let (res, borrow) = sub_with_borrow !lx !b limb_zero in
      set_ofs r !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
898 899
      assert { value r !i - power radix !i * !b =
                 value x !i - y };
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
900
      b := borrow;
901 902 903 904 905
      value_tail r !i;
      value_tail x !i;
      assert { value r (!i+1) - power radix (!i+1) * !b
                  = value x (!i+1) - y
             (*by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
906 907
               value r !i - power radix !i * !b
             = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
908
                 - power radix !i * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
909
             = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
910
                 - power radix k * radix * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
911 912
             = value r k + power radix k * (res - radix * !b)
             = value r k +
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
913
                 (power radix k) * (!lx - (!b at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
914
             = value r k - power radix k * (!b at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
915
                 + power radix k * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
916
             = value x k - y + power radix k * !lx
917
             = value x !i - y*)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
918
      };
919
      i := Int32.(+) !i (Int32.of_int 1);
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
920
    done;
921
    if Int32.(=) !i sz then !b
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
922 923
    else begin
    while Int32.(<) !i sz do
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
924
      invariant { !b = 0 }
925
      invariant { 0 <= !i <= sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
926 927
      invariant { value r !i - power radix !i * !b
                   = value x !i - y }
928
      variant { sz - !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
929 930
      lx := get_ofs x !i;
      set_ofs r !i !lx;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
931 932
      assert { value r !i - power radix !i * !b
                   = value x !i - y };
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
933 934 935 936 937 938 939 940 941 942 943 944
      let ghost k = p2i !i in
      i := Int32.(+) !i (Int32.of_int 1);
      value_sub_tail (pelts r) r.offset (r.offset + k);
      value_sub_tail (pelts x) x.offset (x.offset + k);
    done;
    !b
  end

  (** [sub_limbs r x y sz] substracts [(y, sz)] from [(x, sz)] and
      writes the result to [(r, sz)]. Returns borrow, either 0 or
      1. Corresponds to [mpn_sub_n]. *)
  let sub_limbs (r x y:t) (sz:int32) : limb
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
945 946 947
    requires { valid x sz }
    requires { valid y sz }
    requires { valid r sz }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
948
    ensures { 0 <= result <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
949 950
    ensures { value r sz - power radix sz * result
              = value x sz - value y sz }
951
    writes { r.data.elts }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
952
  =
953
    let limb_zero = Limb.of_int 0 in
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
954 955 956 957 958
    let lx = ref limb_zero in
    let ly = ref limb_zero in
    let b = ref limb_zero in
    let i = ref (Int32.of_int 0) in
    while Int32.(<) !i sz do
959 960
      variant { sz - !i }
      invariant { 0 <= !i <= sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
961 962
      invariant { value r !i - (power radix !i) * !b
                  = value x !i - value y !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
963
      invariant { 0 <= !b <= 1 }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
964 965 966 967 968
      label StartLoop in
      lx := get_ofs x !i;
      ly := get_ofs y !i;
      let res, borrow = sub_with_borrow !lx !ly !b in
      set_ofs r !i res;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
969 970
      assert { value r !i - (power radix !i) * !b =
      value x !i - value y !i };
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
971
      b := borrow;
972 973 974 975 976 977
      value_tail r !i;
      value_tail x !i;
      value_tail y !i;
      assert { value r (!i+1) - (power radix (!i+1)) * !b
                  = value x (!i+1) - value y (!i+1)
               (*by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
978 979
                 value r !i - power radix !i * !b
               = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
980
                   - power radix !i * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
981
               = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
982
                   - power radix k * radix * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
983
               = value r k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
984
                 + power radix k * (res - radix * !b)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
985
               = value r k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
986
                 + power radix k * (!lx - !ly - (!b at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
987
               = value r k - power radix k * (!b at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
988
                 + power radix k * (!lx - !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
989
               = value x k - value y k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
990
                 + power radix k * (!lx - !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
991
               = value x k - value y k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
992
                 + power radix k * !lx - power radix k * !ly
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
993 994 995 996
               = value x k + power radix k * !lx
                 - (value y k + power radix k * !ly)
               = value x !i
                 - (value y k + power radix k * !ly)
997
               = value x !i - value y !i*)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
998
        };
999
      i := Int32.(+) !i (Int32.of_int 1);
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
1000 1001 1002
      done;
      !b

1003
  (** [sub r x y sx sy] substracts [(y,sy)] from [(x, sx)] and writes the
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
1004 1005 1006
      result in [(r, sx)]. [sx] must be greater than or equal to
      [sy]. Returns borrow, either 0 or 1. Corresponds to [mpn_sub]. *)
  let sub (r x y:t) (sx sy:int32) : limb
1007
    requires { 0 <= sy <= sx }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
1008 1009 1010
    requires { valid x sx }
    requires { valid y sy }
    requires { valid r sx }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
1011 1012
    ensures { value r sx  - power radix sx * result
              = value x sx - value y sy }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
1013
    ensures { 0 <= result <= 1 }
1014
    writes { r.data.elts }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
1015 1016 1017 1018 1019 1020 1021 1022
  =
    let limb_zero = Limb.of_int 0 in
    let lx = ref limb_zero in
    let ly = ref limb_zero in
    let b = ref limb_zero in
    let i = ref (Int32.of_int 0) in
    let one = Int32.of_int 1 in
    while Int32.(<) !i sy do
1023 1024
      variant { sy - !i }
      invariant { 0 <= !i <= sy }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
1025 1026
      invariant { value r !i - power radix !i * !b =
                  value x !i - value y !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
1027
      invariant { 0 <= !b <= 1 }