mp2.mlw 286 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

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
103
  function p2i (i:int32) : int = Int32.to_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 159 160 161 162 163 164
  = "vc:sp" if n < m then value_sub_tail x (n+1) m else ()(*assert { 1+2=3 }*)

  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 665 666
      (*assert { value r !i + (power radix !i) * !c =
                value x !i + value y sy };*) (* false without this, cannotreduce with this *)
      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
667
               by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
668 669 670 671 672 673 674
               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
675 676
                  = value x k + value y sy*) };
      i := Int32.(+) !i (Int32.of_int 1);
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
677
    done;
678 679
    !c

Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
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 707 708
  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;
709
      assert { !lx = (pelts ox)[ox.offset + !i] };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
710 711 712
      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
713 714 715 716 717 718
      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
719 720
      assert { value x !i + (power radix !i) * !c = value ox !i + value y !i };
      c := carry;
721 722 723 724 725 726
      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
727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745
                 = 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)
746 747
              = value ox !i + value y !i*) };
      i := Int32.(+) !i (Int32.of_int 1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
    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;
763
      assert { !lx = (pelts ox)[ox.offset + !i] };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
764 765 766 767 768 769 770
      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;
771
      assert { forall j. !i < j < sx ->
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
772
                  (pelts x)[x.offset + j] = (pelts ox) [x.offset + j] };
773 774 775 776 777
      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
778 779 780 781 782 783 784 785 786 787 788 789 790 791 792
                 = 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
793 794
                 + value y sy*) };
      i := Int32.(+) !i (Int32.of_int 1);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
795 796 797 798 799 800 801 802 803 804 805 806 807
    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);
808
      assert { value x sx = value (old x) sx + value y sy };
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
809 810 811 812 813
      n
      end
    end


Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
814 815 816 817
  (** [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
818 819
    requires { valid x sz }
    requires { valid r sz }
820
    requires { 0 < sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
821 822
    ensures { value r sz - power radix sz * result
              = value x sz - y }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
823
    ensures { 0 <= result <= 1 }
824
    writes { r.data.elts }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
825 826 827 828 829
  =
    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
830
    while Int32.(<) !i sz && (not (Limb.(=) !b limb_zero)) do
831 832
      invariant { 0 <= !i <= sz }
      invariant { !i > 0 -> 0 <= !b <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
833 834
      invariant { value r !i - power radix !i * !b
                  = value x !i - y }
835
      variant { sz - !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
836 837 838 839
      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
840 841
      assert { value r !i - power radix !i * !b =
                 value x !i - y };
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
842
      b := borrow;
843 844 845 846 847
      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
848 849
               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
850
                 - power radix !i * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
851
             = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
852
                 - power radix k * radix * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
853 854
             = value r k + power radix k * (res - radix * !b)
             = value r k +
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
855
                 (power radix k) * (!lx - (!b at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
856
             = value r k - power radix k * (!b at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
857
                 + power radix k * !lx
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
858
             = value x k - y + power radix k * !lx
859
             = value x !i - y*)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
860
      };
861
      i := Int32.(+) !i (Int32.of_int 1);
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
862
    done;
863
    if Int32.(=) !i sz then !b
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
864 865
    else begin
    while Int32.(<) !i sz do
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
866
      invariant { !b = 0 }
867
      invariant { 0 <= !i <= sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
868 869
      invariant { value r !i - power radix !i * !b
                   = value x !i - y }
870
      variant { sz - !i }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
871 872
      lx := get_ofs x !i;
      set_ofs r !i !lx;
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
873 874
      assert { value r !i - power radix !i * !b
                   = value x !i - y };
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
875 876 877 878 879 880 881 882 883 884 885 886
      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
887 888 889
    requires { valid x sz }
    requires { valid y sz }
    requires { valid r sz }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
890
    ensures { 0 <= result <= 1 }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
891 892
    ensures { value r sz - power radix sz * result
              = value x sz - value y sz }
893
    writes { r.data.elts }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
894
  =
895
    let limb_zero = Limb.of_int 0 in
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
896 897 898 899 900
    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
901 902
      variant { sz - !i }
      invariant { 0 <= !i <= sz }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
903 904
      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
905
      invariant { 0 <= !b <= 1 }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
906 907 908 909 910
      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
911 912
      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
913 914 915 916 917 918
      b := borrow;
      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);
      value_sub_tail (pelts y) y.offset (y.offset + k);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
919 920
      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
921
               by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
922 923
                 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
924
                   - power radix !i * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
925
               = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
926
                   - power radix k * radix * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
927
               = value r k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
928
                 + power radix k * (res - radix * !b)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
929
               = value r k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
930
                 + power radix k * (!lx - !ly - (!b at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
931
               = value r k - power radix k * (!b at StartLoop)
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
932
                 + power radix k * (!lx - !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
933
               = value x k - value y k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
934
                 + power radix k * (!lx - !ly)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
935
               = value x k - value y k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
936
                 + power radix k * !lx - power radix k * !ly
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
937 938 939 940 941
               = 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 !i
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
942 943 944 945
        };
      done;
      !b

946
  (** [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
947 948 949
      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
950
    requires { 0 <= sy <= sx }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
951 952 953
    requires { valid x sx }
    requires { valid y sy }
    requires { valid r sx }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
954 955
    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
956
    ensures { 0 <= result <= 1 }
957
    writes { r.data.elts }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
958 959 960 961 962 963 964 965
  =
    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
966 967
      variant { sy - !i }
      invariant { 0 <= !i <= sy }
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
968 969
      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
970
      invariant { 0 <= !b <= 1 }
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
971 972 973 974 975
      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
976 977
      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
978 979 980 981 982 983
      b := borrow;
      let ghost k = p2i !i in
      i := Int32.(+) !i one;
      value_sub_tail (pelts r) r.offset (r.offset + k);
      value_sub_tail (pelts x) x.offset (x.offset + k);
      value_sub_tail (pelts y) y.offset (y.offset + k);
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
984 985
      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
986
              by
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
987 988
              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
989
                - power radix !i * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
990
              = value r k + power radix k * res
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
991
                - (power radix k) * radix * !b
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
992
              = value r k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
993
                + power radix k * (res - radix * !b)
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
994
              = value r k
Raphaël Rieu-Helft's avatar
Raphaël Rieu-Helft committed
995
                + power radix k * (!lx - !ly - (!b at StartLoop))
Raphael Rieu-Helft's avatar
Raphael Rieu-Helft committed
996
              = value r k - (power radix k) * (!b at StartLoop)