examples/multiprecision/div.mlw

   1 module Div
   2
   3   use int.Int
   4   use mach.int.Int32
   5   use import mach.int.UInt64GMP as Limb
   6   use int.Power
   7   use ref.Ref
   8   use mach.c.C
   9   use array.Array
  10   use map.Map
  11   use types.Types
  12   use types.Int32Eq
  13   use types.UInt64Eq
  14   use lemmas.Lemmas
  15   use compare.Compare
  16   use util.UtilOld
  17   use add.AddOld
  18   use sub.SubOld
  19   use logical.LogicalUtil
  20   use logical.LogicalOld
  21   use int.EuclideanDivision
  22
  23   (** Based on Niels Möller and Torbjörn Granlund, “Improved
  24       division by invariant integers” 2010 *)
  25
  26   use int.MinMax as MM
  27
  28   predicate reciprocal (v d:limb) =
  29     v = (div (radix*radix - 1) (d)) - radix
  30
  31   let lemma fact_div (x y z:int)
  32     requires { y > 0 }
  33     ensures { div (x + y * z) y = (div x y) + z }
  34   =
  35     assert { div (x + y * z) y = (div x y) + z
  36              by x + y * z = y * (div (x + y * z) y) + mod (x + y * z) y
  37              so mod (x + y * z) y = mod (y * z + x) y = mod x y
  38              so x + y * z = y * (div (x + y * z) y) + mod x y
  39              so x = y * div x y + mod x y
  40              so x + y * z = y * div x y + mod x y + y * z
  41              so y * (div (x + y * z) y) + mod x y
  42                 = y * div x y + mod x y + y * z
  43              so y * (div (x + y * z) y) = y * div x y + y * z
  44                                         = y * ((div x y) + z)
  45              so y <> 0
  46              so div (x + y * z) y = div x y + z }
  47
  48   let invert_limb (d:limb) : limb
  49     requires { d >= div radix 2 }
  50     ensures { reciprocal result d }
  51   =
  52     let v = div2by1 Limb.uint64_max
  53                     (Limb.uint64_max - d)
  54                     d in
  55     fact_div (radix * radix - 1) (l2i d) (- radix);
  56     assert { v = (div (radix*radix - 1) (d)) - radix
  57              by
  58              radix - 1 + radix * (radix - 1 - d)
  59              = radix - 1 + radix * (radix - 1) - radix * d
  60              = radix - 1 + radix * radix - radix - radix * d
  61              = radix * radix - 1 - radix * d
  62              so
  63              radix - 1 + radix * (radix - 1 - d)
  64              = radix * radix - 1 - radix * d
  65              so
  66              v
  67              = div ((radix - 1) + radix * (radix - 1 - d)) (d)
  68              = div (radix * radix - 1 - radix * d) (d)
  69              = div (radix * radix - 1) (d) - radix
  70            };
  71    v
  72
  73   (** Divide a two-word integer by a one-word integer given the
  74   reciprocal of the divisor. *)
  75   let div2by1_inv (uh ul d v:limb) : (limb,limb)
  76     requires { d >= div radix 2 }
  77     requires { uh < d }
  78     requires { reciprocal v d }
  79     returns { q, r -> l2i q * d + l2i r = ul + radix * uh }
  80     returns { _q, r -> 0 <= l2i r < d }
  81   =
  82     let ghost k = radix * radix - (radix + l2i v) * l2i d in
  83     let ghost u = l2i ul + radix * l2i uh in
  84     assert { 1 <= k <= d };
  85     let l,h = mul_double v uh in
  86     let sl,c = add_with_carry l ul 0 in
  87     let (sh,ghost c') = add_with_carry uh h c in  (* <c',sh,sl> = <uh, ul> + <h,l> *)
  88     assert { sl + radix * sh + radix * radix * c'
  89              = l + radix * h + ul + radix * uh };
  90     assert { c' = 0
  91              by
  92              uh < d
  93              so v * uh <= v * d
  94              so k = radix * radix - (radix + v) * d
  95                   = radix * radix - radix * d - v * d
  96              so v * d = radix * radix - radix * d - k
  97                               = radix * (radix - d) - k
  98              so k > 0
  99              so v * d < radix * (radix - d)
 100              so v * uh < radix * (radix - d)
 101              so l + radix * h = v * uh
 102              so l + radix * h < radix * (radix - d)
 103              so uh <= d - 1
 104              so radix * uh <= radix * (d - 1) = radix * d - radix
 105              so l + radix * h + radix * uh
 106                 < radix * (radix - d) + radix * uh
 107                 <= radix * (radix - d) + radix * d - radix
 108                 <= radix * (radix - d + d) - radix = radix * radix - radix
 109              so ul < radix
 110              so l + radix * h + ul + radix * uh
 111                 = l + radix * h + radix * uh + ul
 112                 < radix * radix - radix + ul
 113                 < radix * radix - radix + radix = radix * radix
 114              so sl + radix * sh + radix * radix * c'
 115                 = l + radix * h + ul + radix * uh
 116                 < radix * radix
 117              so radix * radix * c' <= sl + radix * sh + radix * radix * c'
 118              so radix * radix * c' < radix * radix
 119      };
 120     assert { sl + radix * sh = l + radix * h + ul + radix * uh
 121              = v * uh + ul + radix * uh
 122              = ul + (radix + v) * uh };
 123     let qh = ref (sh:limb) in
 124     let ql = ref sl in
 125     let ghost q0 = l2i !ql in
 126     let ghost cq = l2i sh + 1 in (*candidate quotient*)
 127     let ghost cr = l2i ul - cq * l2i d + radix * l2i uh in (*candidate remainder*)
 128     assert { cq * d + cr = u};
 129     qh := add_mod !qh 1;
 130     assert { !qh = mod cq radix };
 131     let p = mul_mod !qh d in
 132     let r = ref (sub_mod ul p) in
 133     let ghost r' = !r in
 134     assert { r' = mod cr radix
 135              by
 136              let a = (- div (!qh * d) radix) in
 137              r' = !r
 138              = mod (ul - p) radix
 139              = mod (ul - mod (!qh * d) radix) radix
 140              = mod (radix * a
 141                    + ul - mod (!qh * d) radix) radix
 142              = mod (ul - mod (!qh * d) radix
 143                            - radix * div (!qh * d) radix) radix
 144              = mod (ul - !qh * d) radix
 145              = mod (ul - mod cq radix * d) radix
 146              = mod (radix * (- (div cq radix)) * d + ul - mod cq radix * d) radix
 147              = mod (ul - (mod cq radix + radix * div cq radix) * d) radix
 148              = mod (ul - cq * d) radix
 149              = mod (radix * uh + ul - cq * d) radix
 150              = mod (ul - cq * d + radix * uh) radix
 151              = mod cr radix };
 152     assert { radix * cr = uh * k + ul * (radix - d) + q0 * d - radix * d };
 153     prod_compat_strict_r (l2i ul) radix (radix - l2i d);
 154     prod_compat_strict_r (l2i d) radix (radix - q0);
 155     assert { (* Theorem 2 of Möller&Granlund 2010 *)
 156              (MM.max (radix - d) (q0 + 1)) - radix <= cr < MM.max (radix - d) q0
 157              by radix * cr = uh * k + ul * (radix - d) + q0 * d - radix * d
 158              so (uh * k + ul * (radix - d) >= 0
 159                 by uh >= 0 /\ k >= 0 /\ ul >= 0 /\ radix - d >= 0)
 160              so radix * cr >= q0 * d - radix * d
 161              so radix * cr >= - radix * d
 162              so cr >= - d
 163              so radix * cr >= q0 * d - radix * d = (q0 - radix) * d
 164              so radix > d
 165              so radix - q0 > 0
 166              so d * (radix-q0) < radix * (radix - q0)
 167              so (q0 - radix) * d > (q0 - radix) * radix
 168              so radix * cr > (q0 - radix) * radix
 169              so cr > q0 - radix
 170              so (let m = MM.max (radix - d) (q0 +1) in
 171                 cr >= m - radix
 172                 by (cr + radix >= - d + radix
 173                    /\ (cr + radix > q0 so cr + radix >= q0 + 1))
 174                    so cr + radix >= m)
 175              so 0 < k <= d so 0 <= uh < d
 176              so k * uh < k * d <= d * d
 177              so radix * cr < d * d + ul * (radix - d) + q0 * d - radix * d
 178              so ul * (radix - d) < radix * (radix - d)
 179              so radix * cr < d * d + radix * (radix - d) + q0 * d - radix * d
 180              so (radix * cr < (radix - d) * (radix - d) + q0 * d
 181                 by
 182                 d * d + radix * (radix - d) + q0 * d - radix * d
 183                         = radix * (radix - d) + d * d - radix * d + q0 * d
 184                         = radix * (radix - d) + (d - radix) * d +  q0 * d
 185                         = radix * (radix - d) - d * (radix - d) + q0 * d
 186                         = (radix - d) * (radix - d) + q0 * d)
 187              so let m = MM.max (radix - d) q0 in
 188                 radix - d <= m
 189              so (radix - d) * (radix - d) <= m* (radix - d)
 190              so (q0 * d <= m * d by 0 <= q0 <= m /\ 0 < d)
 191              so radix * cr < (radix - d) * (radix - d) + q0 * d
 192                            <= m* (radix - d) + q0 * d
 193                            <= m* (radix - d) + m * d
 194                            = m * radix
 195              so cr < m
 196                            };
 197     assert { cr >= 0 -> r' = cr };
 198     assert { cr < 0 ->
 199            ( r' = cr + radix
 200              by cr >= MM.max (radix - d) (q0 + 1) - radix
 201              so cr >= - d
 202              so cr + radix >= radix - d >= 0
 203              so 0 <= cr + radix < radix
 204              so mod (cr + radix) radix = mod cr radix
 205              so r' = mod (cr + radix) radix ) };
 206     assert { cr < 0 ->
 207                 ( !r > !ql
 208                 by MM.max (radix - d) (q0 + 1) >= q0 + 1 > q0
 209                 so cr >= (MM.max (radix - d) (q0 +1)) - radix > q0 - radix
 210                 so r' = cr + radix > q0 - radix + radix = q0 ) };
 211     assert { 1 <= cq <= radix };
 212     assert { (!qh = cq \/ (!qh = 0 /\ cq = radix)
 213               by (1 <= cq < radix -> !qh = mod cq radix = cq)
 214               so (cq = radix -> !qh = 0) ) };
 215     assert { cq = radix ->
 216              (cr < 0
 217                 by cq * d + cr = u
 218                 so uh <= d - 1
 219                 so 1 + uh <= d
 220                 so ul < radix
 221                 so u = ul + radix * uh
 222                      < radix + radix * uh
 223                      = radix * (1 + uh)
 224                      <= radix * d
 225                 so u < radix * d
 226                 so radix * d + cr = u
 227                 so radix * d + cr < radix * d
 228                 so cr < 0) };
 229     assert { 1 <= cq < radix -> !qh = cq /\ !qh * d + cr = u };
 230     if !r > !ql
 231     then
 232     begin
 233       qh := sub_mod !qh 1;
 234       r := add_mod !r d;
 235       assert { cr >= 0 ->
 236                   (!r = cr + d
 237                   by r' = cr
 238                   so r' < MM.max (radix - d) q0
 239                   so r' > q0
 240                   so 0 <= r' < radix - d
 241                   so d <= r' + d < radix
 242                   so !r = mod (r' + d) radix = r' + d) };
 243       assert { cr >= 0 ->
 244                   ( !r >= d
 245                   by r' = cr >= 0
 246                   so !r = r' + d >= d ) };
 247       assert { cr < 0 ->
 248                   ( !r = r' + d - radix
 249                   by r' = cr + radix < radix
 250                   so cr >= MM.max (radix - d) (q0 + 1) - radix
 251                         >= radix - d - radix = - d
 252                   so r' = cr + radix >= radix - d
 253                   so !r = mod (r' + d) radix
 254                   so radix + radix >= r' + d >= radix
 255                   so !r = mod (r' + d) radix = r' + d - radix ) };
 256       assert { cr < 0 ->
 257                   ( 0 <= !r < d
 258                   by r' = cr + radix < radix
 259                   so !r = mod (r' + d) radix = r' + d - radix
 260                   so !r >= 0
 261                   so !r = r' + d - radix < d ) };
 262       assert { cq = radix ->
 263                 ( !qh * d + !r = u
 264                 by cq * d + cr = u
 265                 so cr < 0
 266                 so r' = cr + radix
 267                 so u = radix * d + cr
 268                      = (radix - 1) * d + d + cr
 269                      = (radix - 1) * d + d + r' - radix
 270                 so r' = cr + radix >= MM.max (radix - d) (q0 + 1)
 271                    >= radix - d
 272                 so radix + radix >= d + r' >= radix
 273                 so !r = mod (d + r') radix = d + r' - radix
 274                 so (radix - 1) * d + !r = u
 275                 so !qh = mod ((mod cq radix) - 1) radix
 276                            = mod (-1) radix
 277                            = radix - 1
 278                 so !qh * d + !r = u
 279                ) };
 280       assert { !r = cr + d by [@case_split] cr >= 0 \/ cr < 0 };
 281       assert { 1 <= cq <= radix ->
 282                ( !qh * d + !r = u
 283                by cq * d + cr = u
 284                so !qh = cq - 1
 285                so !qh * d + cr + d = u
 286                so !r = cr + d ) };
 287     end
 288     else
 289     begin
 290        assert { cr >= 0 };
 291        assert { 1 <= cq < radix };
 292     end;
 293     assert { !qh * d + !r = ul + radix * uh
 294             by [@case_split] cq = radix \/ 1 <= cq < radix };
 295     if !r >= d
 296     then begin
 297       assert { cr >= 0 };
 298       assert { !qh < radix - 1
 299                by
 300                !qh * d = ul + radix * uh  - !r
 301                so uh <= d - 1
 302                so ul + radix * uh - !r
 303                   <= ul + radix * (d - 1) - !r
 304                   = ul + radix * d - radix - !r
 305                   = (ul - radix)  + radix * d - !r
 306                   <  radix * d - !r
 307                   <= radix * d - d
 308                   = (radix - 1) * d
 309                so !qh * d < (radix - 1) * d
 310                so d > 0
 311                so !qh < radix - 1 };
 312       qh := !qh + 1;
 313       r := !r - d;
 314       assert { 0 <= !r < d };
 315       assert { !qh * d + !r = ul + radix * uh };
 316     end;
 317     assert { 0 <= !r < d };
 318     assert { !qh * d + !r = ul + radix * uh };
 319     (!qh,!r)
 320
 321 (** `wmpn_divrem_1 q x sz y` divides `(x, sz)` by `y`, writes the quotient
 322     in `(q, sz)` and returns the remainder. Corresponds to
 323     `mpn_divrem_1`. *)
 324 (* TODO develop further decimal points (qxn) *)
 325 let wmpn_divrem_1 (q x:t) (sz:int32) (y:limb) : limb
 326     requires { valid x sz }
 327     requires { valid q sz }
 328     requires { writable q }
 329     requires { 0 < sz }
 330     requires { 0 < y }
 331     ensures { value x sz
 332               = value q sz * y + result }
 333     ensures { result < y }
 334   =
 335     let msb = sz - 1 in
 336     let lx = ref 0 in
 337     let i = ref msb in
 338     let r = ref 0 in
 339     (*normalize divisor*)
 340     let clz = count_leading_zeros y in
 341     if (clz > 0)
 342     then begin
 343       let ghost mult = power 2 (p2i clz) in
 344       let ry = lsl y (Limb.of_int32 clz) in
 345       assert { ry = mult * y };
 346       let ghost tlum = power 2 (Limb.length - p2i clz) in
 347       assert { tlum * mult = radix };
 348       let v = invert_limb ry in
 349       while (!i >= 0) do
 350         variant   { !i }
 351         invariant { -1 <= !i <= msb }
 352         invariant { !r < ry }
 353         invariant { mult * value_sub (pelts x) (x.offset + !i + 1) (x.offset + sz)
 354                     = value_sub (pelts q) (q.offset + !i + 1)
 355                                           (q.offset + sz)
 356                       * ry
 357                       + !r }
 358         invariant { !r <= radix - mult }
 359         invariant { mod (!r) mult = 0 }
 360         assert { !i >= 0 };
 361         label StartLoop in
 362         lx := C.get_ofs x !i;
 363         (*TODO lshift in place would simplify things*)
 364         let l,h = lsld_ext !lx (Limb.of_int32 clz) in
 365         mod_mult mult (l2i y) 0;
 366         assert { !r + h < ry
 367                  by
 368                  let drm = div (!r) mult in
 369                  let dym = div (ry) mult in
 370                  mod (!r) mult = 0
 371                  so !r = mult * drm
 372                  so mod (ry) mult
 373                     = mod (mult * (y) + 0) mult
 374                     = mod 0 mult
 375                     = 0
 376                  so ry = mult * dym
 377                  so !r < ry
 378                  so 0 < ry - !r
 379                       = mult * dym - mult * drm
 380                       = mult * (dym - drm)
 381                  so mult > 0
 382                  so dym - drm > 0
 383                  so dym >= drm + 1
 384                  so h < mult
 385                  so !r + h = mult * drm + h
 386                     < mult * drm + mult
 387                     = mult * (drm + 1)
 388                     <= mult * dym = l2i ry };
 389         assert { !r + h < radix by
 390                  !r + h < ry < radix };
 391         let (qu,rem) = div2by1_inv (!r + h) l ry v in
 392         mod_mult mult (l2i y * l2i qu) (l2i rem);
 393         mod_mult mult (tlum * (l2i !r + l2i h)) (l2i l);
 394         assert { mod (rem) mult = 0
 395                  by
 396                  ry * qu + rem
 397                  = (radix * (!r + h) + l)
 398                  so
 399                  mult * y * qu + rem
 400                  = (mult * tlum * (!r + h) + l)
 401                  so mod (mult * y * qu + rem) mult
 402                     = mod (mult * tlum * (!r + h) + l) mult
 403                  so mult > 0
 404                  so mod (mult * (y * qu) + rem) mult
 405                     = mod (rem) mult
 406                  so mod (mult * tlum * (!r + h) + l) mult
 407                     = mod (l) mult
 408                     = 0
 409                  };
 410         let ghost mer = div (l2i rem) mult in
 411         assert { rem <= radix - mult
 412                  by
 413                  mod (rem) mult = 0
 414                  so mult * mer = l2i rem < radix = mult * tlum
 415                  so mult > 0
 416                  so 0 < mult * tlum - mult * mer = mult * (tlum - mer)
 417                  so tlum - mer > 0
 418                  so mer < tlum
 419                  so rem = mult * mer <= mult * (tlum - 1) = radix - mult
 420                  };
 421         r:=rem;
 422         assert { qu * ry + !r = l + radix * h + radix * (!r at StartLoop) };
 423                (* coerced div2by1 postcondition *)
 424         value_sub_update_no_change (pelts q) (q.offset + p2i !i)
 425                                    (q.offset + 1 + p2i !i)
 426                                    (q.offset + p2i sz) qu;
 427         C.set_ofs q !i qu;
 428         assert { mult * value_sub (pelts x) (x.offset + !i + 1) (x.offset + sz)
 429                     = value_sub (pelts q) (q.offset + !i + 1)
 430                                           (q.offset + sz)
 431                       * ry
 432                       + (!r at StartLoop) }; (* previous invariant is still true *)
 433         value_sub_head (pelts x) (x.offset + int32'int !i) (x.offset + p2i sz);
 434         value_sub_head (pelts q) (q.offset + int32'int !i) (q.offset + p2i sz);
 435         assert { l + radix * h = mult * !lx }; (*lsld_ext postcondition *)
 436         assert { mult * value_sub (pelts x) (x.offset + !i)
 437                                             (x.offset + sz)
 438                  = mult * !lx
 439                    + radix * (mult * value_sub (pelts x) (x.offset + !i + 1)
 440                                               (x.offset + sz))
 441                  by (pelts x)[x.offset + !i] = !lx
 442                  so value_sub (pelts x) (x.offset + !i) (x.offset + sz)
 443                     = !lx + radix * value_sub (pelts x) (x.offset + !i + 1)
 444                                               (x.offset + sz) }; (*nonlinear*)
 445         assert { value_sub (pelts q) (q.offset + !i) (q.offset + sz) * ry
 446                  = qu * ry
 447                    + radix
 448                       * (value_sub (pelts q) (q.offset + !i + 1) (q.offset + sz)
 449                         * ry)
 450                  by (pelts q)[q.offset + !i] = qu
 451                  so value_sub (pelts q) (q.offset + !i) (q.offset + sz)
 452                     = qu + radix * value_sub (pelts q) (q.offset + !i + 1)
 453                                               (q.offset + sz) }; (*nonlinear*)
 454         assert { mult * value_sub (pelts x) (x.offset + !i)
 455                                             (x.offset + sz)
 456                 = value_sub (pelts q) (q.offset + !i) (q.offset + sz) * ry
 457                   + !r };
 458         i := !i - 1;
 459       done;
 460       let ghost res = lsr !r (Limb.of_int32 clz) in
 461       assert { value x sz = value q sz * y + res
 462                by !r = res * mult
 463                so mult * value x sz
 464                   = value q sz * ry + !r
 465                   = value q sz * y * mult + !r
 466                   = value q sz * y * mult + res * mult
 467                   = (value q sz * y + res) * mult };
 468       lsr !r (Limb.of_int32 clz) end
 469     else begin
 470       let v = invert_limb y in
 471       while (!i >= 0) do
 472         variant   { !i }
 473         invariant { -1 <= !i <= msb }
 474         invariant { !r < y }
 475         invariant { value_sub (pelts x) (x.offset + !i + 1) (x.offset + sz)
 476                     = value_sub (pelts q) (q.offset + !i + 1)
 477                                           (q.offset + sz)
 478                       * y
 479                       + !r }
 480         assert { !i >= 0 };
 481         label StartLoop in
 482         let ghost k = p2i !i in
 483         lx := C.get_ofs x !i;
 484         let (qu, rem) = div2by1_inv !r !lx y v in
 485         r := rem;
 486         value_sub_update_no_change (pelts q) (q.offset + p2i !i)
 487                                    (q.offset + 1 + p2i !i)
 488                                    (q.offset + p2i sz) qu;
 489         C.set_ofs q !i qu;
 490         i := !i - 1;
 491         value_sub_head (pelts x) (x.offset + k) (x.offset + p2i sz);
 492         value_sub_head (pelts q) (q.offset + k) (q.offset + p2i sz);
 493         assert { value_sub (pelts x) (x.offset + !i + 1) (x.offset + sz)
 494                     = value_sub (pelts q) (q.offset + !i + 1)
 495                                           (q.offset + sz)
 496                       * y
 497                       + !r
 498                  by (pelts q)[q.offset + k] = qu
 499                  so (pelts x)[x.offset + k] = !lx
 500                  so value_sub (pelts x) (x.offset + !i + 1) (x.offset + sz)
 501                     = value_sub (pelts x) (x.offset + k) (x.offset + sz)
 502                     = !lx + radix * value_sub (pelts x) (x.offset + k + 1)
 503                                                         (x.offset + sz)
 504                     = !lx + radix * (value_sub (pelts q) (q.offset + k + 1)
 505                                                          (q.offset + sz)
 506                                      * y + (!r at StartLoop))
 507                     = !lx + radix * (!r at StartLoop)
 508                           + radix * (value_sub (pelts q) (q.offset + k + 1)
 509                                                          (q.offset + sz)
 510                                      * y)
 511                     = qu * y + !r
 512                           + radix * (value_sub (pelts q) (q.offset + k + 1)
 513                                                          (q.offset + sz)
 514                                      * y)
 515                     = (qu + radix * value_sub (pelts q) (q.offset + k + 1)
 516                                                         (q.offset + sz))
 517                       * y
 518                       + !r
 519                     = value_sub (pelts q) (q.offset + !i + 1)
 520                                           (q.offset + sz)
 521                       * y
 522                       + !r };
 523       done;
 524       !r
 525     end
 526
 527
 528   predicate reciprocal_3by2 (v dh dl:limb) =
 529     v = div (radix*radix*radix -1) (dl + radix * dh) - radix
 530
 531   let div3by2_inv (uh um ul dh dl v: limb) : (limb,limb,limb)
 532     requires { dh >= div radix 2 }
 533     requires { reciprocal_3by2 v dh dl }
 534     requires { um + radix * uh < dl + radix * dh }
 535     returns { q, rl, rh -> uint64'int q * dl + radix * q * dh
 536                            + uint64'int rl + radix * uint64'int rh
 537                   = ul + radix * um + radix * radix * uh }
 538     returns { _q, rl, rh -> 0 <= uint64'int rl + radix * uint64'int rh < dl + radix * dh }
 539   = [@vc:do_not_keep_trace] (* traceability info breaks the proof *)
 540     let ghost d = l2i dl + radix * l2i dh in
 541     let ghost u = l2i ul + radix * (l2i um + radix * l2i uh) in
 542     let q1 = ref 0 in
 543     let r0 = ref 0 in
 544     let r1 = ref 0 in
 545     let l,h = mul_double v uh in
 546     let sl, c = add_with_carry um l 0 in
 547     let sh, ghost c' = add_with_carry uh h c in
 548     assert { sl + radix * sh + radix * radix * c'
 549              = um + radix * uh + v * uh };
 550     assert { c' = 0
 551              by
 552              um + radix * uh < d
 553              so radix * uh < d
 554              so radix * radix * radix - 1
 555                 = d * (div (radix * radix * radix - 1) d)
 556                   + mod (radix * radix * radix - 1) d
 557                 >= d * (div (radix * radix * radix - 1) d)
 558              so radix * (um + radix * uh + v * uh)
 559                 < radix * (d + v * uh)
 560                 = radix * d + v * radix * uh
 561                 <= radix * d + v * d
 562                 = (div (radix * radix * radix - 1) d) * d
 563                 <= radix * radix * radix - 1
 564                 < radix * radix * radix
 565              so um + radix * uh + v * uh < radix * radix
 566              so sl + radix * sh + radix * radix * c' < radix * radix
 567              so radix * radix * c' < radix * radix
 568      };
 569     q1 := sh;
 570     let ghost q0 = l2i sl in
 571     let ghost cq = l2i !q1 + 1 in (*candidate quotient*)
 572     q1 := add_mod !q1 1;
 573     assert { !q1 = mod cq radix };
 574     let p = mul_mod dh sh in
 575     r1 := sub_mod um p;
 576     label CQuot in
 577     let ghost a = div (l2i um - l2i dh * l2i sh) radix in
 578     let tl, th = mul_double sh dl in
 579     let il, b = sub_with_borrow ul tl 0 in
 580     let (ih, ghost b') = sub_with_borrow !r1 th b in
 581     assert { il + radix * ih - radix * radix * b'
 582              = ul + radix * !r1 - sh * dl };
 583     let bl,b2 = sub_with_borrow il dl 0 in
 584     let bh, ghost b2' = sub_with_borrow ih dh b2 in
 585     assert { bl + radix * bh - radix * radix * b2'
 586              = il + radix * ih - dl - radix * dh };
 587     mod_mult (radix * radix) (uint64'int b')
 588              (uint64'int ul + radix * uint64'int !r1
 589                 - uint64'int sh * uint64'int dl - uint64'int dl
 590                 - radix * uint64'int dh);
 591     assert { bl + radix * bh
 592              = mod (ul + radix * !r1
 593                 - sh * dl- dl
 594                 - radix * dh) (radix * radix)
 595              by
 596              mod (il + radix * ih - dl - radix * dh) (radix * radix)
 597              = mod (radix * radix * (-b2') + bl + radix * bh) (radix * radix)
 598              = mod (bl + radix * bh) (radix * radix)
 599              = bl + radix * bh
 600              so
 601              bl + radix * bh
 602              = mod (il + radix * ih
 603                 - dl - radix * dh) (radix * radix)
 604              so il + radix * ih
 605                 = radix * radix * b' + ul + radix * !r1
 606                   - sh * dl
 607              so mod (il + radix * ih
 608                 - dl - radix * dh) (radix * radix)
 609                 = mod (radix * radix * b' + ul + radix * !r1
 610                       - sh * dl - dl - radix * dh)
 611                   (radix * radix)
 612                 = mod (ul + radix * !r1
 613                 - sh * dl - dl
 614                 - radix * dh) (radix * radix) };
 615     r1 := bh;
 616     r0 := bl;
 617     let ghost r' = l2i !r0 + radix * l2i !r1 in
 618     let ghost cr = u - d * cq in
 619     assert {  r' = mod cr(radix * radix)
 620               by
 621               (!r1 at CQuot = mod (um - dh * sh) radix
 622                 by let a' = div (dh * sh) radix in
 623                    dh * sh = p + radix * a'
 624                 so !r1 at CQuot = mod (um - p) radix
 625                    = mod (radix * a' + um - dh * sh) radix
 626                    = mod (um - dh * sh) radix )
 627               so um - dh * sh = a * radix + !r1 at CQuot
 628               so !r0 + radix * !r1
 629                  = mod (ul + radix * (!r1 at CQuot)
 630                   - sh * dl - dl
 631                   - radix * dh) (radix * radix)
 632               so ul + radix * (!r1 at CQuot)
 633                    - sh * dl - dl - radix * dh
 634                  = ul + radix * (um - dh * sh - a * radix)
 635                    - sh * dl - dl - radix * dh
 636                  = ul + radix * um - radix * dh * sh
 637                    - radix * radix * a - sh * dl - dl
 638                    - radix * dh
 639                  = ul + radix * um - radix * dh * (sh + 1)
 640                    - radix * radix * a - sh * dl - dl
 641                  = ul + radix * um - radix * dh * (sh + 1)
 642                    - radix * radix * a - dl * (sh + 1)
 643                  = ul + radix * um
 644                    - (dl + radix * dh) * (sh + 1)
 645                    - radix * radix * a
 646                  = ul + radix * um - d * cq - radix * radix * a
 647                  = u - radix * radix * uh - d * cq - radix * radix * a
 648                  = cr + radix * radix * (- a - uh)
 649               so (*let y = - a - uh in*)
 650                  mod (ul + radix * (!r1 at CQuot)
 651                   - sh * dl - dl
 652                   - radix * dh) (radix * radix)
 653                  = mod (radix * radix * (-a - uh) + cr)
 654                        (radix * radix)
 655                  = mod cr (radix*radix)
 656                };
 657     let ghost m = MM.max (radix * radix - d) (q0 * radix) in
 658     assert { (* Theorem 3 of Moller&Granlund 2010 *)
 659              m - radix * radix <= cr < m
 660             by
 661             let k = radix * radix * radix - (radix + v) * d in
 662             reciprocal_3by2 v dh dl
 663             so let m3 = radix * radix * radix - 1 in
 664                (radix + v) * d = d * div m3 d = m3 - mod m3 d
 665             so (k = 1 + mod m3 d
 666                 by k = radix * radix * radix - (radix + v) * d
 667                      = m3 + 1 - (radix + v) * d
 668                      = m3 + 1 - m3 + mod m3 d
 669                      = 1 + mod m3 d)
 670             so 1 <=  k <= d
 671             so q0 + radix * sh = (radix + v) * uh + um
 672             so cq = sh + 1
 673             so radix * cq = radix * sh + radix
 674                = (radix + v) * uh + um - q0 + radix
 675             so (radix * cr = k * uh + (radix * radix - d) * um
 676                  + radix * ul + d * q0 - d * radix
 677                by radix * cr = radix * (u - cq * d)
 678                = radix * u
 679                  - ((radix + v) * uh + um - q0 + radix) * d
 680                = radix * u - d * (radix + v) * uh
 681                  - d * um + d * q0 - d * radix
 682                = radix * u - (radix * radix * radix - k) * uh
 683                  - d * um + d * q0 - d * radix
 684                = (radix * radix * radix * uh + radix * radix * um
 685                  + radix * ul) - (radix * radix * radix - k) * uh
 686                  - d * um + d * q0 - d * radix
 687                = k * uh + radix * radix * um + radix * ul
 688                  - d * um + d * q0 - d * radix
 689                = k * uh + (radix * radix - d) * um + radix * ul
 690                  + d * q0 - d * radix )
 691             so (cr >= m - radix * radix
 692                by (
 693                   k >= 0 so radix * radix - d >= 0
 694                so uh >= 0 so um >= 0 so ul >= 0
 695                so  k * uh + (radix * radix - d) * um + radix * ul
 696                    >= 0
 697                so radix * cr >= d * q0 - d * radix
 698                so q0 >= 0 so d >= 0
 699                so d * q0 >= 0
 700                so radix * cr >= - d * radix
 701                so cr >= -d = radix * radix - d - radix * radix
 702                so radix * cr >= d * (q0 - radix)
 703                so (
 704                     (radix - q0) * d < (radix - q0) * radix * radix
 705                     by let rq = radix - q0 in let r2 = radix * radix in
 706                     rq > 0 /\ d < r2
 707                     so rq * d < rq * r2
 708                   )
 709                so d * (q0 - radix) > radix * radix * (q0 - radix)
 710                so radix * cr > radix * radix * (q0 - radix)
 711                so cr > radix * (q0 - radix) = radix * q0 - radix * radix
 712                   ))
 713             so cr < m
 714             by (
 715                 let bbd = radix * radix - d in
 716                 bbd > 0 /\ bbd <= m /\ q0 * radix <= m
 717                 so (bbd * bbd <= bbd * m
 718                     by [@case_split]
 719                     (bbd = m \/ (bbd < m so bbd * bbd < bbd * m)))
 720                 so (d*(radix * q0) <= d * m
 721                     by [@case_split]
 722                     (radix * q0 = m \/ (radix * q0 < m so d > 0 so d * (radix * q0) < d * m)))
 723                 so if uh <= dh - 1
 724                 then
 725                   let dm = dh - 1 in
 726                   uh <= dm
 727                   so
 728                   k * uh <= k * dm
 729                   so (k * dm <= d * dm
 730                      by k <= d /\ 0 <= dm
 731                      so [@case_split] (k = d \/ dm = 0 \/
 732                                      (k < d /\ dm > 0 so k * dm < d * dm)))
 733                   so k * uh <= d * dm
 734                   so
 735                   bbd * um <= bbd * (radix - 1)
 736                   so
 737                   radix * cr
 738                   = k * uh + (radix * radix - d) * um
 739                     + radix * ul + d * q0 - radix * d
 740                   <= d * dm + bbd * um
 741                     + radix * ul + d * q0 - radix * d
 742                   <= d * dm + bbd * (radix - 1)
 743                     + radix * ul + d * q0 - radix * d
 744                   < d * dm + bbd * (radix - 1)
 745                     + radix * radix + d * q0 - radix * d
 746                   so radix * radix * cr
 747                   < radix * (d * dm + bbd * (radix - 1)
 748                     + radix * radix + d * q0 - radix * d)
 749                   = d * radix * (dh - 1) + bbd * radix * (radix - 1)
 750                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 751                   = d * radix * dh - d * radix + bbd * radix * (radix - 1)
 752                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 753                   = d * (d - dl) - d * radix + bbd * radix * (radix - 1)
 754                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 755                   = d * d - d * radix + bbd * radix * (radix - 1)
 756                     + radix * radix * radix + radix * d * q0 - radix * radix * d - d * dl
 757                   so (d * dl >= 0 by d >= 0 /\ dl >= 0)
 758                   so radix * radix * cr
 759                   < d * d - d * radix + bbd * radix * (radix - 1)
 760                     + radix * radix * radix + radix * d * q0 - radix * radix * d - d * dl
 761                   <= d * d - d * radix + bbd * radix * (radix - 1)
 762                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 763                   = d * d - d * radix + bbd * (radix * radix - radix)
 764                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 765                   = d * d - d * radix + bbd * radix * radix - (radix * radix - d) * radix
 766                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 767                   = d * d - d * radix + bbd * radix * radix
 768                     + radix * d - radix * radix * radix
 769                     + radix * radix * radix + radix * d * q0 - radix * radix * d
 770                   = d * d + bbd * radix * radix - radix * radix * d + radix * d * q0
 771                   = bbd * radix * radix - d * (radix * radix - d) + radix * d * q0
 772                   = bbd * radix * radix - d * bbd + radix * d * q0
 773                   = bbd * bbd + d * (radix * q0)
 774                   <= bbd * m + d * (radix * q0)
 775                   <= bbd * m + d * m
 776                   = radix * radix * m
 777                   so cr < m
 778                  else
 779                   uh = dh
 780                   so
 781                   (um <= dl - 1
 782                   by um + radix * uh < dl + radix * dh)
 783                   so (radix * radix - d) * um <= (radix * radix - d) * (dl - 1)
 784                   so
 785                   ( radix * radix * cr
 786                     < radix * radix * m
 787                     - (radix - dl) * (radix * radix * radix - d * (1+ radix))
 788                  by radix * cr
 789                   = k * dh + (radix * radix - d) * um
 790                     + radix * ul + d * q0 - radix * d
 791                   <= d * dh + (radix * radix - d) * um
 792                     + radix * ul + d * q0 - radix * d
 793                   <= d * dh + (radix * radix - d) * (dl - 1)
 794                     + radix * ul + d * q0 - radix * d
 795                   < d * dh + (radix * radix - d) * (dl - 1)
 796                     + radix * radix + d * q0 - radix * d
 797                   so radix * radix * cr
 798                   < radix * (d * dh + (radix * radix - d) * (dl - 1)
 799                     + radix * radix + d * q0 - radix * d)
 800                   = d * radix * dh
 801                     + (radix * radix - d) * (dl - 1) * radix
 802                     + radix * radix * radix + d * q0 * radix - radix * radix * d
 803                   = d * (d - dl)
 804                     + (radix * radix - d) * (radix * dl - radix)
 805                     + radix * radix * radix + d * q0 * radix - radix * radix * d
 806                   = d * d - d * dl + radix * radix * radix * dl
 807                     - d * radix * dl + d * radix - radix * radix * radix
 808                     + radix * radix * radix + d * q0 * radix - radix * radix * d
 809                   = d * d - d * dl + radix * radix * radix * dl
 810                     - d * radix * dl + d * radix + d * q0 * radix
 811                     - radix * radix * d
 812                   = d * d - radix * radix * d + d * radix + d * q0 * radix
 813                     + dl * (radix * radix * radix - d - d * radix)
 814                   = d * (d - radix * radix) + d * radix + d * q0 * radix
 815                     + dl * (radix * radix * radix - d - d * radix)
 816                   = bbd * (-d) + d * radix + d * q0 * radix
 817                     + dl * (radix * radix * radix - d - d * radix)
 818                   = bbd * (bbd - radix * radix) + d * radix + d * q0 * radix
 819                     + dl * (radix * radix * radix - d - d * radix)
 820                   = bbd * bbd + d * q0 * radix
 821                     - bbd * radix * radix + d * radix
 822                     + dl * (radix * radix * radix - d * (1 + radix))
 823                   = bbd * bbd + d * q0 * radix
 824                     - (radix * radix - d) * radix * radix + d * radix
 825                     + dl * (radix * radix * radix - d * (1 + radix))
 826                   = bbd * bbd + d * q0 * radix
 827                     - radix * ((radix * radix - d) * radix - d)
 828                     + dl * (radix * radix * radix - d * (1 + radix))
 829                   = bbd * bbd + d * q0 * radix
 830                     - radix *  (radix * radix * radix - d * radix - d)
 831                     + dl * (radix * radix * radix - d * (1 + radix))
 832                   = bbd * bbd + d * q0 * radix
 833                     - radix * (radix * radix * radix - d * (1+ radix))
 834                     + dl * (radix * radix * radix - d * (1 + radix))
 835                   = bbd * bbd + d * q0 * radix
 836                     - (radix - dl) * (radix * radix * radix - d * (1+ radix))
 837                   <= bbd * m + d * q0 * radix
 838                     - (radix - dl) * (radix * radix * radix - d * (1+ radix))
 839                   <= bbd * m + d * m
 840                     - (radix - dl) * (radix * radix * radix - d * (1+ radix))
 841                   = (bbd + d) * m
 842                     - (radix - dl) * (radix * radix * radix - d * (1+ radix))
 843                   = radix * radix * m
 844                     - (radix - dl) * (radix * radix * radix - d * (1+ radix))
 845                   )
 846                   so
 847                     (cr < m by
 848                     if d <= radix * (radix - 1)
 849                     then (radix + 1) * d <= radix * (radix - 1) * (radix + 1)
 850                          = radix * (radix * radix - 1)
 851                          = radix * radix * radix - radix
 852                          < radix * radix * radix
 853                          so (radix * radix * radix - d * (1+ radix)) > 0
 854                          so radix - dl > 0
 855                          so (radix - dl) * (radix * radix * radix
 856                                                - d * (1+ radix))
 857                             > 0
 858                          so
 859                          radix * radix * cr
 860                          < radix * radix * m
 861                            - (radix - dl) * (radix * radix * radix
 862                            - d * (1+ radix))
 863                          < radix * radix * m
 864                          so radix * radix * cr < radix * radix * m
 865                     else
 866                         dl + radix * dh = d > radix * (radix - 1)
 867                         so dl < radix
 868                         so dl + radix * dh < radix * (1 + dh)
 869                         so radix - 1 < 1 + dh
 870                         so dh > radix - 2
 871                         so dh = radix - 1
 872                         so uh = dh
 873                         so d >= radix * (radix - 1) +1
 874                         so d * (radix + 1)
 875                            >= (radix * (radix - 1) + 1) * (radix +1)
 876                             = radix * (radix * radix - 1) + radix + 1
 877                             = radix * radix * radix - radix + radix + 1
 878                             = radix * radix * radix + 1
 879                         so
 880                         (d * div (radix * radix * radix - 1) d
 881                              <= d * (radix + 1)
 882                           by d * div (radix * radix * radix - 1) d
 883                              <= radix * radix * radix - 1
 884                              < radix * radix * radix + 1
 885                              <= d * (radix + 1))
 886                         so (let a = div (radix * radix * radix - 1) d in
 887                             a < radix + 1
 888                             by d > 0
 889                             so (forall x y z. x * z < y * z /\ z > 0 -> x < y)
 890                             so (forall x y. x * d < y * d -> x < y)
 891                             so let r = radix + 1 in
 892                                a * d < r * d
 893                                so a < r)
 894                         so v = div (radix * radix * radix - 1) d - radix
 895                                  < radix + 1 - radix = 1
 896                         so v = 0
 897                         so sh = uh = radix - 1
 898                         so cq = sh + 1 = radix
 899                         so cr = u - cq * d
 900                               = u - radix * d
 901                               = ul + radix * (um + radix * dh)
 902                                        - radix * (dl + radix * dh)
 903                               = ul + radix * (um - dl)
 904                         so um <= dl - 1
 905                         so 1 + um - dl <= 0
 906                         so ul < radix
 907                         so cr = ul + radix * (um - dl)
 908                               < radix + radix * (um - dl)
 909                               = radix * (1 + um - dl) <= 0
 910                         so cr < 0 <= m
 911                    )
 912                 )
 913             };
 914     assert { cr >= 0 -> r' = cr };
 915     assert { cr < 0 -> r' = radix * radix + cr
 916              by
 917              m >= radix * radix - d
 918              so cr >= m - radix * radix >= -d
 919              so cr + radix * radix >= radix * radix - d >= 0
 920              so 0 <= cr + radix * radix < radix * radix
 921              so mod (radix * radix + cr) (radix*radix) = mod cr (radix*radix)
 922              so r' = mod (radix * radix + cr) (radix*radix) };
 923     assert { cr < 0 -> !r1 >= sl
 924              by m >= radix * q0
 925              so cr >= m - radix * radix >= radix * q0 - radix * radix
 926              so r' = radix * radix + cr >= radix * q0
 927              so r' = radix * !r1 + !r0 >= radix * q0
 928              so !r0 < radix
 929              so r' < radix * !r1 + radix = radix * (!r1 + 1)
 930              so radix * q0 < radix * (!r1 + 1)
 931              so sl = q0 < !r1 + 1 };
 932     assert { 1 <= cq <= radix };
 933     assert { 1 <= cq < radix -> !q1 = cq so !q1 * d + cr = u };
 934     assert { cq = radix ->
 935              (cr < 0
 936                 by cq * d + cr = u
 937                 so um + radix * uh <= d - 1
 938                 so radix * d + cr = ul
 939                                     + radix * (um + radix * uh)
 940                                   <= ul + radix * (d - 1)
 941                                   = ul - radix + radix * d
 942                                   < radix * d
 943              )
 944            };
 945     label PreCorrections in
 946     if !r1 >= sl
 947     then begin
 948       q1 := sub_mod !q1 1;
 949       assert { !q1 = cq - 1
 950                by
 951                if cq = radix
 952                then
 953                  (!q1 at PreCorrections)
 954                  = mod cq radix = mod radix radix= 0
 955                  so !q1 = mod (0 - 1) radix = radix - 1 = cq - 1
 956                else
 957                  0 <= cq - 1 < radix - 1
 958                  so (!q1 at PreCorrections) = cq
 959                  so !q1 = mod (cq - 1) radix = cq - 1
 960                  };
 961       let rl, c = add_with_carry !r0 dl 0 in
 962       let rh, ghost c' = add_with_carry !r1 dh c in
 963       assert { rl + radix * rh = mod (r' + d) (radix * radix)
 964                by radix * radix * c' + rl + radix * rh
 965                   = r' + d
 966                so mod (r' + d) (radix * radix)
 967                   = mod (radix * radix * c' + rl + radix * rh)
 968                       (radix * radix)
 969                   = mod (rl + radix * rh) (radix * radix)  };
 970       assert { rl + radix * rh = cr + d
 971                by
 972                if cr >= 0
 973                then r' = cr
 974                     so rl + radix * rh = mod (cr + d) (radix * radix)
 975                     so cr < MM.max (radix * radix - d) (q0*radix)
 976                     so (cr >= q0 * radix
 977                        by
 978                           r' = radix * !r1 + !r0
 979                              >= radix * !r1
 980                              >= radix * q0)
 981                     so cr < radix * radix - d
 982                     so cr + d < radix * radix
 983                     so (cr + d >= 0  by cr + d >= cr)
 984                     so mod (cr + d) (radix * radix) = cr + d
 985                else
 986                     r' = cr + radix * radix
 987                     so cr >= m - radix * radix
 988                     so r' >= m >= radix * radix - d
 989                     so r' + d >= radix * radix
 990                     so r' < radix * radix
 991                     so d < radix * radix
 992                     so r' + d < radix * radix + radix * radix
 993                     so mod (r' + d) (radix * radix)
 994                        = r' + d - radix * radix
 995                        = cr + d
 996              };
 997       r1 := rh;
 998       r0 := rl;
 999       assert { !q1 * d + !r0 + radix * !r1 = u
1000                by
1001                cq * d + cr = u
1002                so !q1 = cq - 1
1003                so !r0 + radix * !r1 = cr + d
1004                so !q1 * d + !r0 + radix * !r1
1005                   = (cq - 1) * d + cr + d
1006                   = cq * d - d + cr + d
1007                   = cq * d + cr };
1008     end
1009     else assert { !q1 * d + r' = u
1010                   by cr >= 0
1011                   so r' = cr
1012                   so 1 <= cq < radix
1013                   so !q1 * d + cr = u };
1014     assert { !q1 * d + !r0 + radix * !r1 = u };
1015     label PreRemAdjust in
1016     if [@extraction:unlikely] (!r1 > dh) || (!r1 = dh && !r0 >= dl)
1017     then begin
1018       let bl, b = sub_with_borrow !r0 dl 0 in
1019       let bh, ghost b'= sub_with_borrow !r1 dh b in
1020       assert { b' = 0 };
1021       assert { bl + radix * bh = !r0 + radix * !r1 - d };
1022       assert { !q1 < radix - 1
1023                by !q1 * d + !r0 + radix * !r1 = u
1024                so !r0 + radix * !r1 >= d
1025                so um + radix * uh <= d - 1
1026                so u = ul + radix * (um + radix * uh)
1027                     <= ul + radix * (d - 1)
1028                     < radix + radix * (d-1)
1029                     = radix * d
1030                so (!q1 * d < (radix - 1) * d
1031                   by
1032                   !q1 * d = u - (!r0 + radix * !r1)
1033                               <= u - d
1034                               < radix * d - d
1035                               = (radix - 1) * d )
1036                };
1037       q1 := add_mod !q1 1;
1038       assert { !q1 = (!q1 at PreRemAdjust) + 1 };
1039       r1 := bh;
1040       r0 := bl;
1041       assert { !q1 * d + !r0 + radix * !r1 = u
1042                by
1043                !q1 * d + !r0 + radix * !r1
1044                = ((!q1 at PreRemAdjust) + 1) * d
1045                  + (!r0 + radix * !r1 at PreRemAdjust) - d
1046                = (!q1 * d + !r0 + radix * !r1 at PreRemAdjust)
1047                 };
1048     end;
1049     assert { 0 <= !r0 + radix * !r1 < d };
1050     (!q1,!r0,!r1)
1051
1052   let lemma bounds_imply_rec3by2 (v dh dl:limb)
1053     requires { radix * radix * radix - (dl + radix * dh)
1054                <= (radix + v) * (dl + radix * dh)
1055                < radix * radix * radix }
1056     ensures { reciprocal_3by2 v dh dl }
1057   = ()
1058
1059   let reciprocal_word_3by2 (dh dl:limb) : limb
1060     requires { dh >= div radix 2 }
1061     ensures { reciprocal_3by2 result dh dl }
1062   =
1063     let ghost d = l2i dl + radix * l2i dh in
1064     let v = ref (invert_limb dh) in
1065     assert { radix * radix - dh
1066              <= (radix + !v) * dh
1067              < radix * radix
1068              by
1069              radix + !v = div (radix * radix - 1) (dh) };
1070     let p = ref (mul_mod dh !v) in
1071     assert { (radix + !v) * dh
1072              = radix * (radix-1)
1073                + !p
1074              by
1075              mod ((radix + !v) * dh) radix
1076              = mod (radix * dh + dh * !v) radix
1077              = mod (dh * !v) radix = l2i !p
1078              so
1079              div ((radix + !v) * dh) radix = radix - 1
1080              so
1081              (radix + !v) * dh
1082              = radix * div ((radix + !v) * dh) radix
1083                + mod (dh * !v) radix
1084              = radix * (radix - 1) + !p
1085              };
1086     label Estimate in
1087     p := add_mod !p dl;
1088     if !p < dl (* carry out *)
1089     then begin
1090       assert { (!p at Estimate) + dl >= radix };
1091       assert { (!p at Estimate) + dl = radix + !p };
1092       assert { !v >= 1
1093                by
1094                (!p at Estimate) + dl >= radix
1095                so (!p at Estimate) > 0
1096              };
1097       assert { (radix + !v) * dh + dl
1098                = radix * (radix - 1) + radix + !p };
1099       label Carry in
1100       if !p >= dh
1101       then begin
1102         v := !v - 1;
1103         p := !p - dh;
1104         assert { (radix + !v) * dh + dl
1105                  = radix * (radix - 1) + radix + !p
1106                };
1107       end;
1108       label Borrow in
1109       v := !v - 1;
1110       assert { !p < dh };
1111       p := sub_mod !p dh;
1112       assert { !p = radix + !p at Borrow - dh };
1113     end;
1114     assert { (radix + !v) * dh * radix + radix * dl
1115                = radix * radix * (radix - 1) + radix * !p
1116              by (radix + !v) * dh + dl
1117                  = radix * (radix - 1) + !p };
1118     assert { radix * radix - dh
1119              <= (radix + !v) * dh + dl
1120              < radix * radix };
1121     let tl, th = mul_double !v dl in
1122     label Adjust in
1123     p := add_mod !p th;
1124     if !p < th (* carry out *)
1125     then begin
1126       assert { (!p at Adjust) + th >= radix };
1127       assert { (!p at Adjust) + th = radix + !p
1128                by (!p at Adjust) + th < radix + radix
1129                so div ((!p at Adjust) + th) radix = 1
1130                so !p = mod ((!p at Adjust) + th) radix
1131                so (!p at Adjust) + th
1132                   = radix * div ((!p at Adjust) + th) radix
1133                     + mod ((!p at Adjust) + th) radix
1134                   = radix + !p
1135              };
1136       assert { !v >= 1
1137                by
1138                th <> 0
1139                so !v <> 0
1140              };
1141       if !p > dh || (!p = dh && tl >= dl)
1142       then begin
1143         assert { tl + radix * !p >= d };
1144         v := !v - 1;
1145         assert { (radix + !v) * dh * radix + radix * dl
1146                    + !v * dl
1147                  = radix * radix * radix
1148                    + radix * !p + tl - d
1149                  by
1150                  (radix + !v) * dh * radix + radix * dl
1151                    + !v * dl
1152                  = (radix + !v at Adjust - 1) * dh * radix
1153                    + radix * dl
1154                    + (!v at Adjust - 1) * dl
1155                  = (radix + !v at Adjust) * dh * radix
1156                    + radix *  dl
1157                    + (!v at Adjust) * dl - radix * dh
1158                    - dl
1159                  = radix * radix * (radix - 1) + radix * (!p at Adjust)
1160                    + (!v at Adjust) * dl - radix * dh
1161                    - dl
1162                  = radix * radix * (radix - 1) + radix * (!p at Adjust)
1163                    + radix * th + tl - d
1164                  = radix * radix * (radix - 1) + radix * (radix + !p)
1165                    + tl - d
1166                  = radix * radix * (radix - 1) + radix * radix + radix * !p
1167                    + tl - d
1168                  = radix * radix * radix + radix * !p + tl - d
1169             };
1170       end;
1171       assert { radix * radix * radix
1172                <= (radix + !v) * dh * radix + radix * dl
1173                    + !v * dl
1174                < radix * radix * radix + d };
1175       v := !v - 1;
1176     end;
1177     bounds_imply_rec3by2 !v dh dl;
1178     !v
1179
1180   (* `(x, sz)` is normalized if its most significant bit is set. *)
1181   predicate normalized (x:t) (sz:int32) =
1182     valid x sz
1183     /\ (pelts x)[x.offset + sz - 1] >= div radix 2
1184
1185   use mul.Mul
1186
1187   let div_sb_qr (q x:t) (sx:int32) (y:t) (sy:int32) : limb
1188     requires { 3 <= sy <= sx }
1189     requires { valid x sx }
1190     requires { valid y sy }
1191     requires { valid q (sx - sy) }
1192     requires { writable q }
1193     requires { writable x }
1194     requires { normalized y sy }
1195     ensures { value (old x) sx =
1196               (value q (sx - sy)
1197               + power radix (sx - sy) * result)
1198                 * value y sy
1199               + value x sy }
1200     ensures { value x sy < value y sy }
1201     ensures { 0 <= result <= 1 }
1202   = [@vc:do_not_keep_trace] (* traceability info breaks the proof *)
1203     let xp = ref (C.incr x (sx - 2)) in
1204     let qp = ref (C.incr q (sx - sy)) in
1205     let dh = C.get_ofs y (sy - 1) in
1206     assert { dh >= div radix 2 by normalized y sy };
1207     let dl = C.get_ofs y (sy - 2) in
1208     let v = reciprocal_word_3by2 dh dl in
1209     let i = ref (sx - sy) in
1210     let mdn = 2 - sy in
1211     let ql = ref 0 in
1212     let xd = C.incr !xp mdn in
1213     let ghost vy = value y (p2i sy) in
1214     let x1 = ref 0 in
1215     let x0 = ref 0 in
1216     let r = wmpn_cmp xd y sy in
1217     let qh = (*begin
1218                ensures { r >= 0 -> result = 1 }
1219                ensures { r < 0 -> result = 0 }*)
1220                if (r >= 0)
1221                then 1
1222                else 0
1223              (*end*) in
1224     label PreAdjust in
1225     begin
1226       ensures { value (old x) sx =
1227                   (value !qp (sx - sy - !i)
1228                    + qh * power radix (sx - sy - !i))
1229                   * vy * power radix !i
1230                   + value x (sy + !i - 1)
1231                   + power radix (sy + !i - 1) * !x1 }
1232       ensures { value_sub (pelts x) (!xp.offset + mdn)
1233                  (!xp.offset + mdn + sy - 1)
1234                  + power radix (sy - 1) * !x1
1235                  < vy }
1236       ensures { dl + radix * dh
1237                   >= (pelts x)[(!xp).offset] + radix * !x1 }
1238     let ghost ox = pelts x in
1239     begin [@vc:sp]
1240     if (not (qh = 0))
1241     then begin
1242          assert { qh = 1 };
1243          let ghost b = wmpn_sub_n_in_place xd y sy in
1244          begin
1245              ensures { value (x at PreAdjust) sx =
1246                   (value !qp (sx - sy - !i)
1247                    + qh * power radix (sx - sy - !i))
1248                   * vy * power radix !i
1249                   + value x sx }
1250              ensures { value_sub (pelts x) (!xp.offset + mdn)
1251                  (!xp.offset + mdn + sy)
1252                  < vy }
1253            value_sub_upper_bound (pelts x) xd.offset (xd.offset + p2i sy);
1254            assert { b = 0 };
1255            assert { value (xd at PreAdjust) sy
1256                     = value xd sy + vy };
1257            value_sub_concat (pelts x) x.offset xd.offset (xd.offset + p2i sy);
1258            value_sub_concat ox x.offset xd.offset (xd.offset + p2i sy);
1259            value_sub_frame (pelts x) ox x.offset xd.offset;
1260            assert { value (x at PreAdjust) sx
1261              = value x sx
1262              + power radix (sx - sy) * vy
1263              by
1264              value_sub (pelts (x at PreAdjust)) x.offset xd.offset
1265              = value_sub (pelts x) x.offset xd.offset
1266              so pelts (xd at PreAdjust) = pelts (x at PreAdjust)
1267              so value_sub (pelts (x at PreAdjust)) xd.offset (xd.offset + sy)
1268                 = value (xd at PreAdjust) sy
1269              so value (x at PreAdjust) sx
1270              = value_sub (pelts (x at PreAdjust)) x.offset xd.offset
1271                + power radix (sx - sy)
1272                * value_sub (pelts (x at PreAdjust)) xd.offset (xd.offset + sy)
1273              = value_sub (pelts x) x.offset xd.offset
1274                + power radix (sx - sy)
1275                * value (xd at PreAdjust) sy
1276              = value_sub (pelts x) x.offset xd.offset
1277                + power radix (sx - sy)
1278                * (value xd sy + vy)
1279              = value_sub (pelts x) x.offset xd.offset
1280                + power radix (sx - sy)
1281                * (value_sub (pelts x) (xd.offset) (xd.offset + sy) + vy)
1282              = value_sub (pelts x) x.offset xd.offset
1283                + power radix (sx - sy)
1284                * value_sub (pelts x) (xd.offset) (xd.offset + sy)
1285                + power radix (sx - sy) * vy
1286              = value x sx
1287                + power radix (sx - sy) * vy
1288            };
1289            value_sub_tail (pelts x) x.offset (x.offset + p2i sy + p2i !i - 1);
1290            assert { value (x at PreAdjust) sx =
1291                   (value !qp (sx - sy - !i)
1292                    + qh * power radix (sx - sy - !i))
1293                   * vy * power radix !i
1294                   + value x sx
1295                  by
1296                    !i = sx - sy
1297                  so power radix (sx - sy - !i) = 1
1298                  so value !qp (sx - sy - !i) = 0 };
1299            value_sub_lower_bound_tight (pelts y) y.offset (y.offset + p2i sy);
1300            assert { value_sub (pelts x) (!xp.offset + mdn)
1301                  (!xp.offset + mdn + sy)
1302                  < vy
1303                  by
1304                  value_sub (pelts x) (!xp.offset + mdn)
1305                            (!xp.offset + mdn + sy)
1306                  = value xd sy
1307                  = value (xd at PreAdjust) sy - vy
1308                  so value (xd at PreAdjust) sy
1309                     < power radix sy
1310                  so vy >= dh * power radix (sy - 1)
1311                  so 2 * vy >= 2 * dh * power radix (sy - 1)
1312                  so 2 * dh >= radix
1313                  so 2 * dh * power radix (sy - 1) >= radix * power radix (sy - 1)
1314                  so 2 * vy >= radix * power radix (sy - 1) = power radix sy
1315                  so value (xd at PreAdjust) sy < 2 * vy
1316                  so value (xd at PreAdjust) sy - vy < vy };
1317            end
1318          end
1319     else begin
1320          assert { value_sub (pelts x) (!xp.offset + mdn)
1321                  (!xp.offset + mdn + sy)
1322                  < vy
1323                  by r < 0 };
1324          assert { value (x at PreAdjust) sx =
1325                   (value !qp (sx - sy - !i)
1326                    + qh * power radix (sx - sy - !i))
1327                    * vy * power radix !i
1328                    + value x sx
1329                  by qh = 0
1330                  so sx - sy - !i = 0
1331                  so (value !qp (sx - sy - !i)
1332                    + qh * power radix (sx - sy - !i)) = 0 };
1333          end
1334     end;
1335     let ghost gx1 = (C.get_ofs !xp 1) in
1336     value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
1337     value_sub_upper_bound_tight (pelts y) y.offset (y.offset + p2i sy - 1);
1338     value_sub_tail (pelts x) (!xp.offset) (!xp.offset + p2i sy - 1);
1339     value_sub_lower_bound_tight (pelts x) (!xp.offset) (!xp.offset + p2i sy - 1);
1340     assert { dl + radix * dh
1341              >= (pelts x)[(!xp).offset] + radix * gx1
1342              by value_sub (pelts x) (!xp.offset + mdn)
1343                                     (!xp.offset + mdn + sy)
1344                 < vy
1345              so value y (sy - 1) < (dl + 1) * power radix (sy - 1 - 1)
1346              so vy = dh * power radix (sy - 1)
1347                      + value y (sy - 1)
1348                      < dh * power radix (sy - 1)
1349                      + (dl + 1) * power radix (sy - 1 - 1)
1350                    = power radix (sy - 2) * (dl+1 + radix * dh)
1351              so !xp.offset + mdn + sy - 1 = !xp.offset + 1
1352              so (pelts x)[!xp.offset + mdn + sy - 1]
1353                 = (pelts x)[!xp.offset + 1] = gx1
1354              so value_sub (pelts x) (!xp.offset + mdn) (!xp.offset + mdn + sy)
1355                 = gx1 * power radix (sy - 1)
1356                      + value_sub (pelts x) (!xp.offset + mdn)
1357                                            (!xp.offset + mdn + sy - 1)
1358                 >= gx1 * power radix (sy - 1)
1359                    + (pelts x)[!xp.offset] * power radix (sy - 1 - 1)
1360                 = power radix (sy - 2)
1361                   * ((pelts x) [!xp.offset] + radix * gx1)
1362              so power radix (sy - 2) * ((pelts x) [!xp.offset] + radix * gx1)
1363                 < power radix (sy - 2) * (dl+1 + radix * dh)
1364              so (pelts x) [!xp.offset] + radix * gx1
1365                 < dl + 1 + radix * dh
1366              };
1367     value_sub_tail (pelts x) (!xp.offset + p2i mdn)
1368                    (!xp.offset + p2i mdn + p2i sy - 1);
1369     value_sub_tail (pelts x) x.offset (x.offset + p2i sy + p2i !i - 1);
1370     x1 := (C.get_ofs !xp 1);
1371     assert { value_sub (pelts x) (!xp.offset + mdn)
1372                  (!xp.offset + mdn + sy - 1)
1373                  + power radix (sy - 1) * !x1
1374                  < vy
1375              by
1376              !xp.offset + mdn + sy - 1 = !xp.offset + 1
1377              so value_sub (pelts x) (!xp.offset + mdn)
1378                  (!xp.offset + mdn + sy - 1)
1379                  + power radix (sy - 1) * !x1
1380              = value_sub (pelts x) (!xp.offset + mdn)
1381                  (!xp.offset + mdn + sy - 1)
1382                  + power radix (sy - 1) * (pelts x)[!xp.offset + 1]
1383              = value_sub (pelts x) (!xp.offset + mdn)
1384                  (!xp.offset + mdn + sy - 1)
1385                  + power radix (sy - 1)
1386                    * (pelts x)[!xp.offset + mdn + sy - 1]
1387              = value_sub (pelts x) (!xp.offset + mdn)
1388                          (!xp.offset + mdn + sy)
1389              < vy };
1390     assert { value (x at PreAdjust) sx =
1391                   (value !qp (sx - sy - !i)
1392                    + qh * power radix (sx - sy - !i))
1393                   * vy * power radix !i
1394                   + value x (sy + !i - 1)
1395                   + power radix (sy + !i - 1) * !x1
1396                 by value (x at PreAdjust) sx =
1397                    (value !qp (sx - sy - !i)
1398                    + qh * power radix (sx - sy - !i))
1399                    * vy * power radix !i
1400                    + value x sx
1401                 so sx = sy + !i
1402                 so x.offset + sy + !i - 1 = !xp.offset + 1
1403                 so (pelts x)[x.offset + sy + !i - 1] =
1404                    (pelts x)[!xp.offset + 1]= !x1
1405                 so value x sx
1406                    = value x (sx - 1)
1407                      + power radix (sx -1) * (pelts x)[x.offset + sx - 1]
1408                    = value x (sy + !i - 1)
1409                      + power radix (sy + !i - 1) * (pelts x)[x.offset + sy + !i - 1]
1410                 so value x (sy + !i - 1)
1411                   + power radix (sy + !i - 1) * !x1
1412                   = value x sx
1413            };
1414     end;
1415     while (!i > 0) do
1416       variant { !i }
1417       invariant { 0 <= !i <= sx - sy }
1418       invariant { (!qp).offset = q.offset + !i }
1419       invariant { (!xp).offset = x.offset + sy + !i - 2 }
1420       invariant { plength !qp = plength q }
1421       invariant { !qp.min = q.min }
1422       invariant { !qp.max = q.max }
1423       invariant { plength !xp = plength x }
1424       invariant { !xp.min = x.min }
1425       invariant { !xp.max = x.max }
1426       invariant { pelts !qp = pelts q }
1427       invariant { pelts !xp = pelts x }
1428       invariant { writable !qp /\ writable !xp }
1429       invariant { value (old x) sx =
1430                   (value !qp (sx - sy - !i)
1431                    + qh * power radix (sx - sy - !i))
1432                   * vy * power radix !i
1433                   + value x (sy + !i - 1)
1434                   + power radix (sy + !i - 1) * !x1 }
1435       invariant { value_sub (pelts x) (!xp.offset + mdn)
1436                  (!xp.offset + mdn + sy - 1)
1437                  + power radix (sy - 1) * !x1
1438                  < vy }
1439       invariant { dl + radix * dh
1440                   >= (pelts x)[(!xp).offset] + radix * !x1 }
1441       label StartLoop in
1442       let ghost k = int32'int !i in
1443       i := !i - 1;
1444       let ghost s = int32'int sy + int32'int !i - 1 in
1445       xp.contents <- C.incr !xp (-1);
1446       let xd = C.incr !xp mdn in
1447       let nx0 = C.get_ofs !xp 1 in
1448       if [@extraction:unlikely] (!x1 = dh && nx0 = dl)
1449       then begin
1450         ql := Limb.uint64_max;
1451         value_sub_concat (pelts x) x.offset xd.offset (xd.offset + p2i sy);
1452         value_sub_tail (pelts xd) xd.offset (xd.offset + p2i sy - 1);
1453         let ghost vlx = value xd (p2i sy - 1) in
1454         assert { value xd sy
1455                  = vlx + power radix (sy - 1) * dl
1456                  by value xd sy
1457                     = vlx + power radix (sy - 1)
1458                             * (pelts xd)[xd.offset + sy - 1]
1459                  so xd.offset + sy - 1 = !xp.offset + mdn + sy - 1
1460                                            = !xp.offset + 1
1461                  so pelts xd = pelts !xp
1462                  so (pelts xd)[xd.offset + sy - 1]
1463                     = (pelts !xp)[!xp.offset + 1] = dl
1464                };
1465         value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
1466         value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 2);
1467         let ghost vly = value y (p2i sy - 2) in
1468         assert { vy = vly + power radix (sy - 2) * dl
1469                           + power radix (sy - 1) * dh
1470                  by (pelts y)[y.offset + sy - 1] = dh
1471                  so (pelts y)[y.offset + sy - 2] = dl
1472                  so
1473                     vy = value y (sy - 1)
1474                           + power radix (sy - 1) * dh
1475                     = vly + power radix (sy - 2) * dl
1476                           + power radix (sy - 1) * dh };
1477         begin
1478           ensures { value_sub (pelts xd) (xd.offset+1) (xd.offset + p2i sy - 1)
1479                         + power radix (sy - 2) * dl
1480                         + power radix (sy - 1) * dh
1481                      < vy }
1482           value_sub_tail (pelts xd) (xd.offset + 1) (xd.offset + p2i sy - 2);
1483           assert { value_sub (pelts x) (!xp.offset at StartLoop + mdn)
1484                              (!xp.offset at StartLoop + mdn + sy - 1)
1485                    = value_sub (pelts xd) (xd.offset+1) (xd.offset + p2i sy - 1)
1486                         + power radix (sy - 2) * dl
1487                    by
1488                    pelts x = pelts xd
1489                    so !xp.offset at StartLoop + mdn = xd.offset + 1
1490                    so !xp.offset at StartLoop + mdn + sy - 1 = xd.offset + sy
1491                    so xd.offset + sy - 1 = !xp.offset + 1
1492                    so pelts xd = pelts !xp
1493                    so (pelts xd)[xd.offset + sy - 1] = (pelts !xp)[!xp.offset+1] = dl
1494                    so value_sub (pelts x) (!xp.offset at StartLoop + mdn)
1495                              (!xp.offset at StartLoop + mdn + sy - 1)
1496                     = value_sub (pelts xd) (xd.offset+1) (xd.offset + sy)
1497                     = value_sub (pelts xd) (xd.offset+1) (xd.offset + p2i sy - 1)
1498                       + power radix (sy - 2)
1499                         * (pelts xd)[xd.offset + p2i sy - 1]
1500                     = value_sub (pelts xd) (xd.offset+1) (xd.offset + p2i sy - 1)
1501                       + power radix (sy - 2) * dl
1502                  };
1503          assert { !x1 = dh };
1504         end;
1505         label SubMax in
1506         let ghost xc = Array.copy (x.data) in
1507         value_sub_frame (pelts x) xc.elts x.offset (x.offset + p2i !i);
1508         let ghost b = wmpn_submul_1 xd y sy !ql in
1509         begin
1510           ensures { value x !i
1511                  = value (x at SubMax) !i }
1512           assert { forall j. x.offset <= j < x.offset + !i
1513                           -> (pelts x)[j] = xc.elts[j]
1514                    by
1515                    (pelts x)[j] = (pelts x at SubMax)[j]
1516                    so
1517                    ((pelts x at SubMax)[j] = xc.elts[j]
1518                    by
1519                    0 <= j /\ j < xc.Array.length
1520                    ) };
1521           value_sub_frame (pelts x) xc.elts x.offset (x.offset + p2i !i);
1522         end;
1523         value_sub_upper_bound (pelts xd) xd.offset (xd.offset + p2i sy);
1524         value_sub_lower_bound (pelts xd) xd.offset (xd.offset + p2i sy);
1525         value_sub_head (pelts xd) xd.offset (xd.offset + p2i sy - 1);
1526         assert { vlx < radix * vly
1527                  by
1528                    vlx = value_sub (pelts xd at SubMax) xd.offset
1529                                    (xd.offset + sy - 1)
1530                        = (pelts xd at SubMax)[xd.offset]
1531                          + radix * value_sub (pelts xd at SubMax)
1532                                              (xd.offset + 1)
1533                                              (xd.offset + sy - 1)
1534                    so value_sub (pelts xd at SubMax) (xd.offset + 1)
1535                                            (xd.offset + sy - 1)
1536                       + power radix (sy - 2) * dl
1537                       + power radix (sy - 1) * dh
1538                       < vy
1539                       = vly + power radix (sy - 2) * dl
1540                             + power radix (sy - 1) * dh
1541                    so value_sub (pelts xd at SubMax) (xd.offset + 1)
1542                                            (xd.offset + sy - 1)
1543                       < vly
1544                    so value_sub (pelts xd at SubMax) (xd.offset + 1)
1545                                            (xd.offset + sy - 1)
1546                       <= vly - 1
1547                    so vlx = (pelts xd at SubMax)[xd.offset]
1548                             + radix * value_sub (pelts xd at SubMax)
1549                                                 (xd.offset + 1)
1550                                                 (xd.offset + sy - 1)
1551                           <= (pelts xd at SubMax)[xd.offset]
1552                             + radix * (vly - 1)
1553                           < radix + radix * (vly - 1)
1554                           = radix * vly
1555                };
1556         assert { b = dh
1557                  by
1558                    value xd sy
1559                  = value (xd at SubMax) sy
1560                    - (!ql) * vy
1561                    + power radix sy * b
1562                  so !ql = radix - 1
1563                  so 0 <= value xd sy < power radix sy
1564                  so radix * power radix (sy - 2) = power radix (sy - 1)
1565                  so radix * power radix (sy - 1) = power radix sy
1566                  so value xd sy
1567                     = power radix (sy - 1) * dl + vlx
1568                       - (radix - 1) * vy
1569                       + power radix sy * b
1570                     = power radix (sy - 1) * dl + vlx
1571                       - radix * (vly + power radix (sy - 2) * dl
1572                                 + power radix (sy - 1) * dh)
1573                       + vy + power radix sy * b
1574                     = power radix (sy - 1) * dl + vlx
1575                       - radix * vly - radix * power radix (sy - 2) * dl
1576                       - radix * power radix (sy - 1) * dh
1577                       + vy + power radix sy * b
1578                     = power radix (sy - 1) * dl + vlx
1579                       - radix * vly - power radix (sy - 1) * dl
1580                       - power radix sy * dh
1581                       + vy + power radix sy * b
1582                     = power radix sy * (b - dh)
1583                       + vlx - radix * vly + vy
1584                  so vlx < radix * vly
1585                  so (0 <= vlx - radix * vly + vy < power radix sy
1586                       by
1587                       vy - radix * vly
1588                       = vly + power radix (sy - 2) * dl
1589                           + power radix (sy - 1) * dh
1590                           - radix * vly
1591                       = power radix (sy - 2) * (dl + radix * dh)
1592                         - vly * (radix - 1)
1593                       so let pr2 = power radix (sy - 2) in
1594                       0 <= vly < pr2
1595                       so 0 <= vly * (radix - 1) < pr2 * (radix - 1)
1596                       so vy - radix * vly
1597                          >= pr2 * (dl + radix * dh)
1598                            - pr2 * (radix - 1)
1599                          = pr2 * (dl + radix * dh - (radix - 1))
1600                       so dh + radix * dh - (radix - 1) >= 0
1601                       so pr2 >= 0
1602                       so vy - radix * vly
1603                          >= pr2 * (dl + radix * dh - (radix - 1)) >= 0
1604                       so vlx - radix * vly < 0
1605                       so vlx - radix * vly + vy < vy < power radix sy
1606                     )
1607                  so - (power radix sy)
1608                     < power radix sy * (b - dh)
1609                     < power radix sy
1610                  so - 1 < b - dh < 1
1611                };
1612         value_sub_concat (pelts x) x.offset xd.offset (x.offset + s);
1613         x1 := C.get_ofs !xp 1;
1614         qp.contents <- C.incr !qp (-1);
1615         value_sub_update_no_change (pelts q) (!qp).offset
1616                             ((!qp).offset + 1)
1617                             ((!qp).offset + p2i sx - p2i sy - p2i !i)
1618                             !ql;
1619         label QUp in
1620         C.set !qp !ql;
1621         assert { value_sub (pelts q) ((!qp).offset + 1)
1622                               ((!qp).offset + sx - sy - !i)
1623                     = value (!qp at StartLoop)
1624                               (sx - sy - k)
1625                 by value_sub (pelts q) ((!qp).offset + 1)
1626                               ((!qp).offset + sx - sy - !i)
1627                    = value_sub (pelts q at QUp) (!qp.offset + 1)
1628                               ((!qp).offset + sx - sy - !i)
1629                    = value (!qp at StartLoop) (sx - sy - k) };
1630         value_sub_head (pelts q) (!qp).offset
1631           ((!qp).offset + p2i sx - p2i sy - p2i !i);
1632         value_sub_tail (pelts x) x.offset (x.offset + p2i sy + p2i !i - 1);
1633         assert { value xd (sy - 1)
1634                  + power radix (sy - 1) * !x1
1635                  = value (xd at SubMax) sy
1636                    + power radix sy * (!x1 at StartLoop)
1637                    - !ql * vy
1638                  by
1639                  value xd sy
1640                  = value (xd at SubMax) sy
1641                    - (!ql) * vy
1642                    + power radix sy * b
1643                  so b = dh = !x1 at StartLoop
1644                  so pelts !xp = pelts x = pelts xd
1645                  so ((pelts xd)[xd.offset + sy - 1] = !x1
1646                     by
1647                     xd.offset = x.offset + !i
1648                     so (!xp).offset = x.offset + !i + sy - 2
1649                     so (!xp).offset + 1 = xd.offset + sy - 1
1650                     so (pelts xd)[xd.offset + sy - 1]
1651                        = (pelts !xp)[(!xp).offset + 1]
1652                        = !x1
1653                     )
1654                  so value xd sy
1655                     = value xd (sy - 1)
1656                       + power radix (sy - 1) * (pelts xd)[xd.offset + sy - 1]
1657                     = value xd (sy - 1)
1658                       + power radix (sy - 1) * !x1
1659                };
1660         (* refl *)
1661         assert { value (old x) sx =
1662                   (value !qp (sx - sy - !i)
1663                    + qh * power radix (sx - sy - !i))
1664                   * vy * power radix !i
1665                   + value x (sy + !i - 1)
1666                   + power radix (sy + !i - 1) * !x1
1667                  by
1668                  pelts !xp = pelts x = pelts xd
1669                  so
1670                  value xd sy
1671                  = value (xd at SubMax) sy
1672                    - (!ql) * vy
1673                    + power radix sy * b
1674                  = value (xd at SubMax) sy
1675                    - (!ql) * vy
1676                    + power radix sy * dh
1677                  so (value x s
1678                     = value x !i
1679                       + power radix !i
1680                         * value xd (sy - 1)
1681                     by
1682                       xd.offset = x.offset + !i
1683                       so x.offset + s = xd.offset + sy - 1
1684                       so value_sub (pelts x) (x.offset + !i) (x.offset + s)
1685                          = value xd (sy - 1)
1686                       so value x s
1687                     = value x !i
1688                       + power radix !i
1689                         * value_sub (pelts x) (x.offset + !i)
1690                                               (x.offset + s)
1691                     = value x !i
1692                       + power radix !i
1693                         * value xd (sy - 1))
1694                  so (power radix s
1695                     = power radix !i * power radix (sy - 1)
1696                     by
1697                     let n = !i in
1698                     let m = sy - 1 in
1699                     let x = radix in
1700                     power x s = power x (n + m)
1701                     so (power x (n + m) = power x n * power x m
1702                        by 0 <= n
1703                        so 0 <= m
1704                        so forall x:int, n:int, m:int.
1705                                  0 <= n -> 0 <= m ->
1706                                  power x (n + m) = (power x n * power x m)))
1707                  so (value x s + power radix s * !x1
1708                     = value x !i
1709                       + power radix !i * (value xd sy)
1710                     by
1711                     value x s + power radix s * !x1
1712                     = value x !i
1713                       + power radix !i
1714                         * value xd (sy - 1)
1715                       + power radix (!i + sy - 1) * !x1
1716                     = value x !i
1717                       + power radix !i *
1718                         (value xd (sy - 1)
1719                          + power radix (sy - 1) * !x1)
1720                     = value x !i
1721                       + power radix !i * (value xd sy)
1722                     )
1723                  so (value (x at StartLoop) (sy + k - 1)
1724                     = value (x at SubMax) !i
1725                       + power radix !i
1726                         * value (xd at SubMax) sy
1727                     by
1728                     pelts xd at SubMax = pelts x at SubMax
1729                     so x.offset at SubMax + !i = xd.offset at SubMax
1730                     so
1731                     value (x at StartLoop) (sy + k - 1)
1732                     = value_sub (pelts x at SubMax) (x at SubMax).offset
1733                                 (xd.offset at SubMax)
1734                       + power radix !i
1735                         * value_sub (pelts x at SubMax)
1736                            (xd.offset at SubMax)
1737                            (xd.offset at SubMax + sy)
1738                     so value_sub (pelts x at SubMax) (x at SubMax).offset
1739                                  (xd at SubMax).offset
1740                        = value (x at SubMax) !i
1741                     so value_sub (pelts x at SubMax) (xd.offset at SubMax)
1742                                                      (xd.offset at SubMax + sy)
1743                        = value (xd at SubMax) sy
1744                     )
1745                  so value x !i
1746                     = value (x at SubMax) !i
1747                  so value x s + power radix s * !x1
1748                     = value (x at StartLoop) (sy + k - 1)
1749                       + power radix !i
1750                         * (value xd sy
1751                            - value (xd at SubMax) sy)
1752                     = value (x at StartLoop) (sy + k - 1)
1753                       + power radix !i
1754                         * (- (!ql) * vy
1755                            + power radix sy * b)
1756                     = value (x at StartLoop) (sy + k - 1)
1757                       + power radix !i
1758                         * (- (!ql) * vy
1759                            + power radix sy * (!x1 at StartLoop))
1760                  so value !qp (sx - sy - !i)
1761                     = !ql + radix *
1762                           value_sub (pelts q) ((!qp).offset + 1)
1763                                     ((!qp).offset + sx - sy - !i)
1764                  so (value_sub (pelts q) ((!qp).offset + 1)
1765                               ((!qp).offset + sx - sy - !i)
1766                     = value (!qp at StartLoop)
1767                               (sx - sy - k)
1768                      by value (!qp at StartLoop) (sx - sy - k)
1769                       = value_sub (pelts q at StartLoop)
1770                                   (!qp.offset + 1) (!qp.offset + sx - sy - !i))
1771                  so value !qp (sx - sy - !i)
1772                     = !ql + radix * value (!qp at StartLoop)
1773                                     (sx - sy - k)
1774                  so power radix (sx - sy - !i)
1775                     = radix * power radix (sx - sy - k)
1776                  so radix * power radix !i = power radix k
1777                  so (power radix !i * power radix sy
1778                       = power radix (sy + k - 1)
1779                       by !i + sy = sy + k - 1
1780                       so power radix !i * power radix sy
1781                          = power radix (!i + sy))
1782                  so (value !qp (sx - sy - !i)
1783                       + qh * power radix (sx - sy - !i))
1784                       * vy * power radix !i
1785                       + value x (sy + !i - 1)
1786                       + power radix (sy + !i - 1) * !x1
1787                     = (!ql + radix * value (!qp at StartLoop)
1788                                     (sx - sy - k)
1789                       + qh * power radix (sx - sy - !i))
1790                       * vy * power radix !i
1791                       + value x (sy + !i - 1)
1792                       + power radix (sy + !i - 1) * !x1
1793                     = (!ql + radix * value (!qp at StartLoop)
1794                                     (sx - sy - k)
1795                       + radix * qh * power radix (sx - sy - k))
1796                       * vy * power radix !i
1797                       + value x (sy + !i - 1)
1798                       + power radix (sy + !i - 1) * !x1
1799                     = (!ql + radix * value (!qp at StartLoop)
1800                                     (sx - sy - k)
1801                       + radix * qh * power radix (sx - sy - k))
1802                       * vy * power radix !i
1803                       + value x s
1804                       + power radix s * !x1
1805                     = !ql * vy * power radix !i
1806                       + radix * (value (!qp at StartLoop)
1807                                     (sx - sy - k)
1808                                  + qh * power radix (sx - sy - k))
1809                         * vy * power radix !i
1810                       + value x s
1811                       + power radix s * !x1
1812                     = !ql * vy * power radix !i
1813                       + (value (!qp at StartLoop)
1814                                     (sx - sy - k)
1815                                  + qh * power radix (sx - sy - k))
1816                         * vy * radix * power radix !i
1817                       + value x s
1818                       + power radix s * !x1
1819                     = !ql * vy * power radix !i
1820                       + (value (!qp at StartLoop)
1821                                     (sx - sy - k)
1822                                  + qh * power radix (sx - sy - k))
1823                         * vy * power radix k
1824                       + value x s
1825                       + power radix s * !x1
1826                     = !ql * vy * power radix !i
1827                       + (value (!qp at StartLoop)
1828                                     (sx - sy - k)
1829                                  + qh * power radix (sx - sy - k))
1830                         * vy * power radix k
1831                       + value (x at StartLoop) (sy + k - 1)
1832                       + power radix !i
1833                         * (- (!ql) * vy
1834                            + power radix sy * (!x1 at StartLoop))
1835                     = (value (!qp at StartLoop)
1836                                     (sx - sy - k)
1837                                  + qh * power radix (sx - sy - k))
1838                         * vy * power radix k
1839                       + value (x at StartLoop) (sy + k - 1)
1840                       + power radix !i * power radix sy
1841                         * (!x1 at StartLoop)
1842                     = (value (!qp at StartLoop)
1843                                     (sx - sy - k)
1844                                  + qh * power radix (sx - sy - k))
1845                         * vy * power radix k
1846                       + value (x at StartLoop) (sy + k - 1)
1847                       + power radix (sy + k - 1) * (!x1 at StartLoop)
1848                     = value (old x) sx
1849                   };
1850         assert { value_sub (pelts x) (!xp.offset + mdn)
1851                    (!xp.offset + mdn + sy - 1)
1852                    + power radix (sy - 1) * !x1
1853                    < vy
1854                    by
1855                      pelts x = pelts xd
1856                    so xd.offset = !xp.offset + mdn
1857                    so !xp.offset + mdn + sy - 1 = xd.offset + sy - 1
1858                    so
1859                      value xd (sy - 1)
1860                      = value_sub (pelts xd) xd.offset (xd.offset + sy - 1)
1861                      = value_sub (pelts x) (!xp.offset + mdn)
1862                        (!xp.offset + mdn + sy - 1)
1863                    so value_sub (pelts x) (!xp.offset + mdn)
1864                                 (!xp.offset + mdn + sy - 1)
1865                       + power radix (sy - 1) * !x1
1866                       = value (xd at SubMax) sy
1867                         + power radix sy * (!x1 at StartLoop)
1868                         - !ql * vy
1869                    so value (xd at SubMax) sy =
1870                       vlx + power radix (sy - 1) * dl
1871                    so vlx < radix * vly
1872                    so (value (xd at SubMax) sy
1873                         + power radix sy * (!x1 at StartLoop)
1874                       < radix * vy
1875                       by
1876                       !x1 at StartLoop = dh
1877                       so power radix sy = radix * power radix (sy - 1)
1878                       so power radix (sy - 1) = radix * power radix (sy - 2)
1879                       so value (xd at SubMax) sy
1880                           + power radix sy * (!x1 at StartLoop)
1881                          = vlx + power radix (sy - 1) * dl
1882                                + power radix sy * dh
1883                          < radix * vly + power radix (sy - 1) * dl
1884                                        + power radix sy * dh
1885                          = radix * vly + radix * power radix (sy - 2) * dl
1886                                        + radix * power radix (sy - 1) * dh
1887                          = radix * (vly + power radix (sy - 2) * dl
1888                                         + power radix (sy - 1) * dh)
1889                          = radix * vy
1890                       )
1891                    so !ql = radix - 1
1892                    so value (xd at SubMax) sy
1893                         + power radix sy * (!x1 at StartLoop)
1894                         - !ql * vy
1895                       < radix * vy - (radix - 1) * vy
1896                       = vy
1897                };
1898           value_sub_tail (pelts x) (!xp.offset + p2i mdn) (!xp.offset);
1899           value_sub_upper_bound (pelts y) (y.offset) (y.offset + p2i sy - 2);
1900           value_sub_lower_bound (pelts x) (!xp.offset + p2i mdn) (!xp.offset);
1901           assert { dl + radix * dh
1902                    >= (pelts x)[(!xp).offset] + radix * !x1
1903                    by
1904                    vy = vly + power radix (sy - 2)
1905                               * (dl + radix * dh)
1906                    so value_sub (pelts x) (!xp.offset + mdn)
1907                    (!xp.offset + mdn + sy - 1)
1908                    + power radix (sy - 1) * !x1
1909                    < vy
1910                    so !xp.offset + mdn + sy - 1 = !xp.offset + 1
1911                    so power radix (sy - 1) = power radix (sy - 2) * radix
1912                    so - mdn = sy - 2
1913                    so vy
1914                       > value_sub (pelts x) (!xp.offset + mdn)
1915                                             (!xp.offset + mdn + sy - 1)
1916                         + power radix (sy - 1) * !x1
1917                       = value_sub (pelts x) (!xp.offset + mdn)
1918                                             (!xp.offset + 1)
1919                         + power radix (sy - 1) * !x1
1920                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
1921                         + power radix (- mdn) * (pelts x)[(!xp).offset]
1922                         + power radix (sy - 1) * !x1
1923                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
1924                         + power radix (sy - 2) * (pelts x)[(!xp).offset]
1925                         + power radix (sy - 1) * !x1
1926                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
1927                         + power radix (sy - 2) * (pelts x)[(!xp).offset]
1928                         + power radix (sy - 2) * radix * !x1
1929                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
1930                         + power radix (sy - 2)
1931                           * ((pelts x)[(!xp).offset] + radix * !x1)
1932                       >= power radix (sy - 2)
1933                           * ((pelts x)[(!xp).offset] + radix * !x1)
1934                    so vly < power radix (sy - 2)
1935                    so vy < power radix (sy - 2)
1936                            + power radix (sy - 2)
1937                               * (dl + radix * dh)
1938                            = power radix (sy - 2)
1939                              * (1 + dl + radix * dh)
1940                    so power radix (sy - 2)
1941                       * ((pelts x)[(!xp).offset] + radix * !x1)
1942                       < power radix (sy - 2) * (1 + dl + radix * dh)
1943                    so power radix (sy - 2)
1944                       * ((pelts x)[(!xp).offset] + radix * !x1
1945                         - (1 + dl + radix * dh))
1946                       < 0
1947                    so (pelts x)[(!xp).offset] + radix * !x1
1948                         - (1 + dl + radix * dh)
1949                       < 0
1950                  };
1951       end
1952       else begin
1953         assert { dl + radix * dh
1954                   > (pelts x)[(!xp).offset + 1] + radix * !x1
1955                   by
1956                   dl + radix * dh
1957                   >= (pelts x)[(!xp).offset + 1] + radix * !x1
1958                   so dh >= !x1
1959                   so [@case_split] dh <> !x1
1960                    \/ (dh = !x1
1961                       /\ dl <> (pelts x)[(!xp).offset + 1])
1962                   so
1963                    [@case_split] dh > !x1 \/
1964                    (dh = !x1 /\ dl > (pelts x)[(!xp).offset + 1])
1965         };
1966         label SmallDiv in
1967         let ghost vlx = value xd (p2i sy - 2) in
1968         let xp0 = C.get !xp in
1969         let xp1 = C.get_ofs !xp 1 in
1970         begin
1971           ensures { value xd sy =
1972                     vlx
1973                     + power radix (sy - 2) * (xp0 + radix * xp1) }
1974         value_sub_tail (pelts xd) xd.offset (xd.offset + p2i sy - 1);
1975         value_sub_tail (pelts xd) xd.offset (xd.offset + p2i sy - 2);
1976         value_sub_upper_bound (pelts xd) xd.offset (xd.offset + p2i sy - 2);
1977         assert { value xd sy
1978                  = vlx + power radix (sy - 2)
1979                          * (xp0 + radix * xp1)
1980                  by xd.offset + sy - 2 = !xp.offset
1981                  so (pelts xd)[xd.offset + sy - 1] = xp1
1982                  so (pelts xd)[xd.offset + sy - 2] = xp0
1983                  so pelts xd = pelts !xp
1984                  so value xd sy
1985                     = value xd (sy - 1)
1986                       + power radix (sy - 1)
1987                             * (pelts xd)[xd.offset + sy - 1]
1988                     = value xd (sy - 2)
1989                       + power radix (sy - 2)
1990                         * (pelts xd)[xd.offset + sy - 2]
1991                       + power radix (sy - 1)
1992                         * (pelts xd)[xd.offset + sy - 1]
1993                     = vlx
1994                       + power radix (sy - 2) * xp0
1995                       + power radix (sy - 1) * xp1
1996                     = value xd (sy - 2)
1997                       + power radix (sy - 2) * xp0
1998                       + power radix (sy - 2) * radix * xp1
1999                     = vlx + power radix (sy - 2)
2000                          * (xp0 + radix * xp1)
2001                };
2002         end;
2003         let qu, rl, rh =
2004             div3by2_inv !x1 xp1 xp0 dh dl v in
2005         ql := qu;
2006         x1 := rh;
2007         x0 := rl;
2008         label SubProd in
2009         value_sub_concat (pelts x) x.offset xd.offset
2010                            (x.offset + p2i sy + k - 1);
2011         let ghost xc = Array.copy (x.data) in
2012         value_sub_frame (pelts x) xc.elts x.offset (x.offset + p2i !i);
2013         let cy = wmpn_submul_1 xd y (sy - 2) !ql in
2014         label PostSub in
2015         begin
2016           ensures { value x !i
2017                  = value (x at SubProd) !i }
2018           assert { forall j. x.offset <= j < x.offset + !i
2019                           -> (pelts x)[j] = xc.elts[j]
2020                    by
2021                    (pelts x)[j] = (pelts x at SubProd)[j]
2022                    so
2023                    ((pelts x at SubProd)[j] = xc.elts[j]
2024                    by
2025                    0 <= j /\ j < xc.Array.length
2026                    ) };
2027           value_sub_frame (pelts x) xc.elts x.offset (x.offset + p2i !i);
2028         end;
2029         let (cy1:limb) = [@vc:sp] if (!x0 < cy) then 1 else 0 in
2030         x0 := sub_mod !x0 cy;
2031         let (cy2:limb) = [@vc:sp] if (!x1 < cy1) then 1 else 0 in
2032         x1 := sub_mod !x1 cy1;
2033         assert { 0 <= cy2 <= 1 };
2034         (* assert { cy2 = 1 -> rh = 0 }; (* and cy > rl *)*)
2035         value_sub_update (pelts x) (!xp).offset xd.offset
2036                                    (xd.offset + p2i sy - 1) !x0;
2037         value_sub_update_no_change (pelts x) (!xp).offset
2038                                    x.offset (x.offset + p2i !i) !x0;
2039         value_sub_update_no_change (pelts x) (!xp).offset
2040                                    xd.offset (xd.offset + p2i sy - 2) !x0;
2041         C.set !xp !x0;
2042         assert { value x !i
2043                  = value (x at SubProd) !i
2044                  by
2045                  value x !i
2046                  = value (x at PostSub) !i
2047                  = value (x at SubProd) !i };
2048         value_sub_tail (pelts x) xd.offset (xd.offset + p2i sy - 1);
2049         begin
2050           ensures { value xd (sy - 1)
2051                    + power radix (sy - 1) * !x1
2052                    - power radix sy * cy2
2053                  = value (xd at SubProd) sy
2054                    + power radix sy * (!x1 at StartLoop)
2055                    - !ql * vy }
2056           assert { value xd (sy - 2)
2057                    = value (xd at PostSub) (sy - 2) };
2058           value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
2059           value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 2);
2060           let ghost vly = value y (p2i sy - 2) in
2061           assert { vy = vly + power radix (sy - 2)
2062                               * (dl + radix * dh)
2063                    by (pelts y)[y.offset + sy - 1] = dh
2064                    so (pelts y)[y.offset + sy - 2] = dl
2065                    so
2066                       vy = value y (sy - 1)
2067                             + power radix (sy - 1) * dh
2068                          = vly + power radix (sy - 2) * dl
2069                             + power radix (sy - 1) * dh
2070                    so power radix (sy - 1)
2071                       = power radix (sy - 2) * radix };
2072           assert { value xd (sy - 2)
2073                    - power radix (sy - 2) * cy
2074                    = vlx - !ql * vly
2075                    by
2076                    value xd (sy - 2)
2077                      - power radix (sy - 2) * cy
2078                    = value (xd at PostSub) (sy - 2)
2079                      - power radix (sy - 2) * cy
2080                    = vlx - !ql * vly
2081                  };
2082           assert { power radix sy
2083                    = power radix (sy - 2) * radix * radix };
2084           assert { xp0 + radix * xp1
2085                    + radix * radix * !x1 at StartLoop
2086                    - !ql * (dl + radix * dh)
2087                   = rl + radix * rh };
2088           begin ensures { value (xd at SubProd) sy
2089                      + power radix sy * (!x1 at StartLoop)
2090                      - !ql * vy
2091                    = value xd (sy - 2)
2092                      - power radix (sy - 2) * cy
2093                      + power radix (sy - 2) *
2094                        (rl + radix * rh) }
2095             assert { value (xd at SubProd) sy
2096                      = vlx + power radix (sy - 2) * xp0
2097                            + power radix (sy - 1) * xp1 }; (*nonlinear*)
2098             assert { !ql * vy = !ql * vly
2099                                 + power radix (sy - 2)
2100                                   * (!ql * (dl + radix * dh)) }; (*nonlinear*)
2101            end;
2102            begin ensures { value xd (sy - 2)
2103                      - power radix (sy - 2) * cy
2104                      + power radix (sy - 2) * (rl + radix * rh)
2105                     =  value xd (sy - 1)
2106                     + power radix (sy - 1) * !x1
2107                     - power radix sy * cy2 }
2108            value_sub_tail (pelts xd) xd.offset (xd.offset + p2i sy - 2);
2109            assert { value xd (sy - 1)
2110                     = value xd (sy - 2)
2111                       + power radix (sy - 2) * !x0
2112                     by (pelts xd)[xd.offset + sy - 2] = !x0
2113                     so value xd (sy - 1)
2114                        = value_sub (pelts xd) xd.offset (xd.offset + sy - 1)
2115                        = value_sub (pelts xd) xd.offset (xd.offset + sy - 2)
2116                          + power radix (sy - 2) * !x0
2117                        = value xd (sy - 2)
2118                          + power radix (sy - 2) * !x0 };
2119            assert { rl + radix * rh - cy
2120                      = !x0 + radix * !x1 - power radix 2 * cy2
2121                      by
2122                         (!x0 - radix * cy1 = rl - cy
2123                          by
2124                          !x0 = mod (rl - cy) radix
2125                          so - radix < rl - cy < radix
2126                          so (if rl < cy
2127                              then cy1 = 1
2128                                   /\ (- radix < rl - cy < 0
2129                                      so
2130                                      div (rl - cy) radix = - 1
2131                                      so rl - cy
2132                                      = radix * div (rl - cy) radix
2133                                        + mod (rl - cy) radix
2134                                      = !x0 - radix
2135                                      = !x0 - radix * cy1)
2136                              else cy1 = 0 /\ rl - cy = l2i !x0)) }
2137                     (* nonlinear *)
2138            end;
2139         end;
2140         if [@extraction:unlikely] (not (cy2 = 0))
2141         then begin
2142           label Adjust in
2143           assert { cy2 = 1 };
2144           begin ensures { !ql > 0 }
2145             value_sub_lower_bound (pelts y) y.offset (y.offset + p2i sy - 1);
2146             value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
2147             value_sub_upper_bound (pelts xd) xd.offset (xd.offset + p2i sy - 1);
2148             assert { !ql > 0
2149                       by
2150                         (value xd (sy - 1)
2151                          + power radix (sy - 1) * !x1
2152                          - power radix sy * cy2
2153                          < 0
2154                          by
2155                          value xd (sy - 1) < power radix (sy - 1)
2156                          so !x1 <= radix - 1
2157                          so  value xd (sy - 1)
2158                              + power radix (sy - 1) * !x1
2159                            < power radix (sy - 1)
2160                              + power radix (sy - 1) * !x1
2161                            = power radix (sy - 1) * (1 + !x1)
2162                            <= power radix (sy - 1) * radix
2163                            = power radix sy
2164                          so  value xd (sy - 1)
2165                              + power radix (sy - 1) * !x1
2166                              - power radix sy * cy2
2167                            < power radix sy - power radix sy * cy2
2168                            = 0
2169                        )
2170                       so value (xd at SubProd) sy
2171                          + power radix sy * (!x1 at StartLoop)
2172                          - !ql * vy
2173                          < 0
2174                       so (value (xd at SubProd) sy
2175                          + power radix sy * (!x1 at StartLoop) >= 0
2176                          by value (xd at SubProd) sy >= 0
2177                          so !x1 at StartLoop >= 0
2178                          so power radix sy * (!x1 at StartLoop) >= 0
2179                          )
2180                       so !ql * vy > 0
2181                       so vy = value_sub (pelts y)
2182                                 y.offset (y.offset + sy - 1)
2183                               + power radix (sy - 1) * dh
2184                       so dh > 0
2185                       so vy > 0
2186                    };
2187           end;
2188           value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
2189           value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 2);
2190           let ghost vly = value y (p2i sy - 2) in
2191           assert { vy = vly + power radix (sy - 2) * dl
2192                             + power radix (sy - 1) * dh
2193                    by (pelts y)[y.offset + sy - 1] = dh
2194                    so (pelts y)[y.offset + sy - 2] = dl
2195                    so
2196                       vy = value y (sy - 1)
2197                             + power radix (sy - 1) * dh
2198                          = vly + power radix (sy - 2) * dl
2199                             + power radix (sy - 1) * dh };
2200           begin
2201             ensures { value xd (sy - 1)
2202                       + power radix (sy - 1) * !x1
2203                       >= power radix sy - vy }
2204             assert { value xd (sy - 1)
2205                       + power radix (sy - 1) * !x1
2206                      = power radix sy + value (xd at SubProd) sy
2207                        + power radix sy * (!x1 at StartLoop)
2208                        - !ql * vy };
2209             assert {  value (xd at SubProd) sy
2210                        + power radix sy * (!x1 at StartLoop)
2211                        - !ql * vy
2212                       >= - vy
2213                       by
2214                          value (xd at SubProd) sy
2215                          = vlx + power radix (sy - 2) * (xp0 + radix * xp1)
2216                       so xp0 + radix * xp1 + radix * radix * (!x1 at StartLoop)
2217                          = !ql * (dl + radix * dh) + rl + radix * rh
2218                       so power radix (sy - 1) = power radix (sy - 2) * radix
2219                       so vy = vly + power radix (sy - 2) * (dl + radix * dh)
2220                       so (!ql * vly < vy
2221                            by
2222                           vly <= power radix (sy - 2)
2223                           so !ql < radix
2224                           so !ql * vly <= !ql * power radix (sy - 2)
2225                                        < radix * power radix (sy - 2)
2226                                        = power radix (sy - 1)
2227                           so vy = vly + power radix (sy - 2) * (dl + radix * dh)
2228                           so dh >= div radix 2 > 1
2229                           so vly >= 0
2230                           so dl >= 0
2231                           so vy >= power radix (sy - 2) * radix * dh
2232                                 > power radix (sy - 2) * radix * 1
2233                                 = power radix (sy - 1)
2234                           )
2235                       so - !ql * vly > - vy
2236                       so vlx >= 0
2237                       so power radix sy = power radix (sy - 2) * radix * radix
2238                       so value (xd at SubProd) sy
2239                        + power radix sy * (!x1 at StartLoop)
2240                        - !ql * vy
2241                          = vlx + power radix (sy - 2) * (xp0 + radix * xp1)
2242                            + power radix sy * (!x1 at StartLoop)
2243                            - !ql * vy
2244                          = vlx + power radix (sy - 2) * (xp0 + radix * xp1)
2245                            + power radix (sy - 2)
2246                              * radix * radix * (!x1 at StartLoop)
2247                            - !ql * vy
2248                          = vlx + power radix (sy - 2)
2249                                  * (xp0 + radix * xp1
2250                                      + radix * radix * (!x1 at StartLoop))
2251                            - !ql * vy
2252                          = vlx + power radix (sy - 2) *
2253                                (!ql * (dl + radix * dh) + rl + radix * rh)
2254                            - !ql * vy
2255                          = vlx + power radix (sy - 2) *
2256                                (!ql * (dl + radix * dh) + rl + radix * rh)
2257                            - !ql * (vly
2258                                     + power radix (sy - 2) * (dl + radix * dh))
2259                          = vlx + power radix (sy - 2) * (rl + radix * rh)
2260                            - !ql * vly
2261                          >= power radix (sy - 2) * (rl + radix * rh)
2262                             - !ql * vly
2263                          >= - !ql * vly > - vy
2264                    };
2265           end;
2266           let ghost xc = Array.copy (x.data) in
2267           assert { forall j. x.offset <= j < x.offset + !i
2268                           -> (pelts x)[j] = xc.elts[j]
2269                    by
2270                    0 <= x.offset <= j /\ j < x.offset + !i <= xc.Array.length
2271                    so 0 <= j < xc.Array.length
2272                  } ;
2273           value_sub_frame (pelts x) xc.elts x.offset (x.offset + p2i !i);
2274           let c = wmpn_add_n_in_place xd y (sy - 1) in
2275           begin
2276           ensures { value x !i
2277                  = value (x at Adjust) !i }
2278           assert { forall j. x.offset <= j < x.offset + !i
2279                           -> (pelts x)[j] = xc.elts[j]
2280                    by
2281                    pelts (xd at Adjust) = pelts (x at Adjust)
2282                    so pelts x = pelts xd
2283                    so (pelts x)[j] = (pelts x at Adjust)[j]
2284                    so
2285                    ((pelts x at Adjust)[j] = xc.elts[j]
2286                    by
2287                    0 <= j /\ j < xc.Array.length
2288                    ) } ;
2289           value_sub_frame (pelts x) xc.elts x.offset (x.offset + p2i !i);
2290           end;
2291           label MidAdd in
2292           begin
2293             ensures { value xd (sy - 1) + power radix (sy - 1) * !x1
2294                       = value (xd at Adjust) (sy - 1)
2295                       + power radix (sy - 1) * (!x1 at Adjust)
2296                       + vy
2297                       - power radix sy }
2298             assert { 0 <= c <= 1
2299                      by
2300                      value xd (sy - 1) + c * power radix (sy - 1)
2301                      = value (xd at Adjust) (sy - 1)
2302                        + value y (sy - 1)
2303                      so
2304                         value (xd at Adjust) (sy - 1)
2305                         < power radix (sy - 1)
2306                      so value y (sy - 1) < power radix (sy - 1)
2307                      so value xd (sy - 1) >= 0
2308                      so c * power radix (sy - 1) < 2 * power radix (sy - 1)
2309                      so let p = power radix (sy - 1) in
2310                         (c < 2 by c * p < 2 * p so p > 0)
2311                    };
2312             let ghost c' = div (l2i !x1 + l2i dh + l2i c) radix in
2313             x1 := add_mod !x1 (add_mod dh c);
2314             assert { !x1 + c' * radix = !x1 at Adjust + dh + c
2315                      by
2316                         (!x1 = mod (!x1 at Adjust + dh + c) radix
2317                         by
2318                         !x1 = mod (!x1 at Adjust + (mod (dh + c) radix)) radix
2319                         so mod (div (dh + c) radix * radix + !x1 at Adjust
2320                                 + mod (dh + c) radix) radix
2321                            = mod (!x1 at Adjust + (mod (dh + c) radix)) radix
2322                         so !x1 = mod (div (dh + c) radix * radix + !x1 at Adjust
2323                                 + mod (dh + c) radix) radix
2324                                = mod (!x1 at Adjust + dh + c) radix
2325                         )
2326                      so (!x1 at Adjust) + dh + c
2327                         = div (!x1 at Adjust + dh + c) radix * radix
2328                           + mod (!x1 at Adjust + dh + c) radix
2329                         = c' * radix + !x1
2330                      };
2331             assert { 0 <= c' <= 1 };
2332             value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
2333             assert { value xd (sy - 1) + power radix (sy - 1) * !x1
2334                      = value (xd at Adjust) (sy - 1)
2335                        + power radix (sy - 1) * (!x1 at Adjust)
2336                        + vy
2337                        - power radix sy
2338                      by
2339                      value xd (sy - 1) + power radix (sy - 1) * c
2340                      = value (xd at Adjust) (sy - 1)
2341                        + value y (sy - 1)
2342                      so vy = value y (sy - 1)
2343                              + power radix (sy - 1) * dh
2344                      so value xd (sy - 1) + power radix (sy - 1) * c
2345                          + power radix (sy - 1) * (!x1 at Adjust)
2346                          + power radix (sy - 1) * dh
2347                         = value (xd at Adjust) (sy - 1)
2348                            + value y (sy - 1)
2349                            + power radix (sy - 1) * (!x1 at Adjust)
2350                            + power radix (sy - 1) * dh
2351                         = value (xd at Adjust) (sy - 1)
2352                            + power radix (sy - 1) * (!x1 at Adjust)
2353                            + vy
2354                      so value xd (sy - 1) + power radix (sy - 1) * c
2355                          + power radix (sy - 1) * (!x1 at Adjust)
2356                          + power radix (sy - 1) * dh
2357                         = value xd (sy - 1)
2358                           + power radix (sy - 1) * (c + dh + !x1 at Adjust)
2359                         = value xd (sy - 1)
2360                           + power radix (sy - 1) * (!x1 + radix * c')
2361                         = value xd (sy - 1)
2362                           + power radix (sy - 1) * !x1
2363                           + power radix sy * c'
2364                       so  value xd (sy - 1)
2365                             + power radix (sy - 1) * !x1
2366                             + power radix sy * c'
2367                           = value (xd at Adjust) (sy - 1)
2368                             + power radix (sy - 1) * (!x1 at Adjust)
2369                             + vy
2370                       so value (xd at Adjust) (sy - 1)
2371                             + power radix (sy - 1) * (!x1 at Adjust)
2372                          >= power radix sy - vy
2373                       so value xd (sy - 1) < power radix (sy - 1)
2374                       so !x1 <= radix - 1
2375                       so power radix (sy - 1) * !x1
2376                          <= power radix (sy - 1) * (radix - 1)
2377                       so value xd (sy - 1)
2378                             + power radix (sy - 1) * !x1
2379                          <= value xd (sy - 1)
2380                             + power radix (sy - 1) * (radix - 1)
2381                          < power radix (sy - 1)
2382                            + power radix (sy - 1) * (radix - 1)
2383                          = power radix sy
2384                       so c' <> 0
2385                       so c' = 1
2386                       };
2387           end;
2388           ql := !ql - 1;
2389           (* todo refl *)
2390           assert { value xd (sy - 1) + power radix (sy - 1) * !x1
2391                    =  value (xd at SubProd) sy
2392                    + power radix sy * (!x1 at StartLoop)
2393                    - !ql * vy
2394                    by
2395                    value xd (sy - 1) + power radix (sy - 1) * !x1
2396                    = value (xd at Adjust) (sy - 1)
2397                        + power radix (sy - 1) * (!x1 at Adjust)
2398                        + vy
2399                        - power radix sy
2400                    = value (xd at SubProd) sy
2401                      + power radix sy * (!x1 at StartLoop)
2402                      - (!ql at Adjust) * vy
2403                      + vy
2404                    = value (xd at SubProd) sy
2405                      + power radix sy * (!x1 at StartLoop)
2406                      - (!ql + 1) * vy
2407                      + vy
2408                    = value (xd at SubProd) sy
2409                      + power radix sy * (!x1 at StartLoop)
2410                      - !ql * vy };
2411           qp.contents <- C.incr !qp (-1);
2412           value_sub_update_no_change (pelts q) (!qp).offset
2413                             ((!qp).offset + 1)
2414                             ((!qp).offset + p2i sx - p2i sy - p2i !i)
2415                             !ql;
2416           C.set !qp !ql;
2417           value_sub_head (pelts q) (!qp).offset
2418             ((!qp).offset + p2i sx - p2i sy - p2i !i);
2419           value_sub_tail (pelts x) x.offset (x.offset + p2i sy + p2i !i - 1);
2420           value_sub_concat (pelts x) x.offset xd.offset (x.offset + s);
2421           (* todo refl *)
2422           assert { value (old x) sx =
2423                   (value !qp (sx - sy - !i)
2424                    + qh * power radix (sx - sy - !i))
2425                   * vy * power radix !i
2426                   + value x (sy + !i - 1)
2427                   + power radix (sy + !i - 1) * !x1
2428                  by
2429                     value !qp (sx - sy - !i)
2430                     = !ql + radix *
2431                           value_sub (pelts q) ((!qp).offset + 1)
2432                                     ((!qp).offset + sx - sy - !i)
2433                  so (value_sub (pelts q) ((!qp).offset + 1)
2434                               ((!qp).offset + sx - sy - !i)
2435                     = value (!qp at StartLoop)
2436                               (sx - sy - k)
2437                      by
2438                      (!qp at StartLoop).offset = (!qp).offset + 1
2439                      so ((!qp).offset + sx - sy - !i)
2440                          - ((!qp).offset + 1)
2441                         = sx - sy - k
2442                      )
2443                  so value !qp (sx - sy - !i)
2444                     = !ql + radix * value (!qp at StartLoop)
2445                                     (sx - sy - k)
2446                  so (value x s
2447                     = value x !i
2448                       + power radix !i
2449                         * value xd (sy - 1)
2450                      by
2451                         xd.offset = x.offset + !i
2452                      so x.offset + s = xd.offset + sy - 1
2453                      so pelts x = pelts xd
2454                      so x.offset + s - xd.offset = sy - 1
2455                      so value_sub (pelts x) xd.offset (x.offset + s)
2456                         = value xd (sy - 1)
2457                      so value x s
2458                         = value_sub (pelts x) x.offset xd.offset
2459                           + power radix !i * value_sub (pelts x) xd.offset (x.offset + s)
2460                         = value x !i
2461                           + power radix !i * value xd (sy - 1)
2462                         )
2463                  so (power radix s
2464                     = power radix !i * power radix (sy - 1)
2465                     by
2466                     let n = !i in
2467                     let m = sy - 1 in
2468                     let x = radix in
2469                     power x s = power x (n + m)
2470                     so (power x (n + m) = power x n * power x m
2471                        by 0 <= n
2472                        so 0 <= m
2473                        so forall x:int, n:int, m:int.
2474                                  0 <= n -> 0 <= m -> power x (n + m) = (power x n * power x m)))
2475                  so (value x s + power radix s * !x1
2476                     = value x !i
2477                       + power radix !i *
2478                         (value (xd at SubProd) sy
2479                         + power radix sy * (!x1 at StartLoop)
2480                         - !ql * vy)
2481                     by  value xd (sy - 1)
2482                          + power radix (sy - 1) * !x1
2483                        = value (xd at SubProd) sy
2484                          + power radix sy * (!x1 at StartLoop)
2485                          - !ql * vy
2486                     so value x s + power radix s * !x1
2487                     = value x !i
2488                       + power radix !i
2489                         * value xd (sy - 1)
2490                       + power radix (!i + sy - 1) * !x1
2491                     = value x !i
2492                       + power radix !i
2493                         * value xd (sy - 1)
2494                       + power radix !i
2495                         * power radix (sy - 1) * !x1
2496                     = value x !i
2497                       + power radix !i *
2498                         (value xd (sy - 1)
2499                          + power radix (sy - 1) * !x1)
2500                     = value x !i
2501                       + power radix !i *
2502                         (value (xd at SubProd) sy
2503                         + power radix sy * (!x1 at StartLoop)
2504                         - !ql * vy)
2505                     )
2506                  so (value (x at StartLoop) (sy + k - 1)
2507                     = value (x at SubProd) !i
2508                       + power radix !i
2509                         * value (xd at SubProd) sy
2510                     by
2511                      value (x at StartLoop) (sy + k - 1)
2512                      = value_sub (pelts x at SubProd) (x at SubProd).offset
2513                                  ((x at SubProd).offset + sy + k - 1)
2514                      = value_sub (pelts x at SubProd) (x at SubProd).offset xd.offset
2515                        + power radix (xd.offset - (x at SubProd).offset)
2516                          * value_sub (pelts x at SubProd) xd.offset
2517                            ((x at SubProd).offset + sy + k - 1)
2518                      so (x at SubProd).offset = x.offset
2519                      so xd.offset = x.offset + !i
2520                      so value_sub (pelts x at SubProd) (x at SubProd).offset xd.offset
2521                         = value (x at SubProd) !i
2522                      so power radix (xd.offset - x.offset) = power radix !i
2523                      so x.offset + sy + k - 1 - xd.offset = p2i sy
2524                      so value_sub (pelts x at SubProd) xd.offset
2525                                   (x.offset + sy + k - 1)
2526                         = value (xd at SubProd) sy
2527                     )
2528                  so (value x !i
2529                     = value (x at SubProd) !i
2530                     by
2531                     value x !i
2532                     = value (x at Adjust) !i
2533                     = value (x at SubProd) !i
2534                     )
2535                  so power radix !i * power radix sy = power radix (!i + sy)
2536                  so value x s + power radix s * !x1
2537                     - value (x at StartLoop) (sy + k - 1)
2538                     = value x !i
2539                       + power radix !i *
2540                         (value (xd at SubProd) sy
2541                         + power radix sy * (!x1 at StartLoop)
2542                         - !ql * vy)
2543                       - (value (x at SubProd) !i
2544                          + power radix !i
2545                          * value (xd at SubProd) sy)
2546                     = value x !i
2547                       + power radix !i *
2548                         (value (xd at SubProd) sy
2549                         + power radix sy * (!x1 at StartLoop)
2550                         - !ql * vy)
2551                       - (value x !i
2552                          + power radix !i
2553                          * value (xd at SubProd) sy)
2554                     =  power radix !i
2555                        * (power radix sy * (!x1 at StartLoop)
2556                           - !ql * vy)
2557                     = power radix !i * power radix sy * (!x1 at StartLoop)
2558                       - power radix !i * !ql * vy
2559                     = power radix (!i + sy) * (!x1 at StartLoop)
2560                       - power radix !i * !ql * vy
2561                     = power radix (sy + k - 1) * (!x1 at StartLoop)
2562                       - power radix !i * !ql * vy
2563                  so value x s + power radix s * !x1
2564                     = value (x at StartLoop) (sy + k - 1)
2565                       + power radix (sy + k - 1) * (!x1 at StartLoop)
2566                       - power radix !i * !ql * vy
2567                  so power radix (sx - sy - !i)
2568                     = radix * power radix (sx - sy - k)
2569                  so radix * power radix !i = power radix k
2570                  so  (value !qp (sx - sy - !i)
2571                        + qh * power radix (sx - sy - !i))
2572                        * vy * power radix !i
2573                        + value x (sy + !i - 1)
2574                        + power radix (sy + !i - 1) * !x1
2575                     = (!ql + radix * value (!qp at StartLoop)
2576                                                       (sx - sy - k)
2577                        + qh * power radix (sx - sy - !i))
2578                        * vy * power radix !i
2579                        + value x (sy + !i - 1)
2580                        + power radix (sy + !i - 1) * !x1
2581                     = (!ql + radix * value (!qp at StartLoop)
2582                                                       (sx - sy - k)
2583                        + qh * radix * power radix (sx - sy - k))
2584                        * vy * power radix !i
2585                        + value x (sy + !i - 1)
2586                        + power radix (sy + !i - 1) * !x1
2587                     = !ql * vy * power radix !i
2588                       + (value (!qp at StartLoop)
2589                                                  (sx - sy - k)
2590                                  + qh * power radix (sx - sy - k))
2591                         * vy * radix * power radix !i
2592                       + value x (sy + !i - 1)
2593                       + power radix (sy + !i - 1) * !x1
2594                     = !ql * vy * power radix !i
2595                       + (value (!qp at StartLoop)
2596                                                  (sx - sy - k)
2597                                  + qh * power radix (sx - sy - k))
2598                         * vy * power radix k
2599                       + value x (sy + !i - 1)
2600                       + power radix (sy + !i - 1) * !x1
2601                     = !ql * vy * power radix !i
2602                       + (value (!qp at StartLoop)
2603                                                  (sx - sy - k)
2604                                  + qh * power radix (sx - sy - k))
2605                         * vy * power radix k
2606                       + value x s
2607                       + power radix s * !x1
2608                     = !ql * vy * power radix !i
2609                       + (value (!qp at StartLoop)
2610                                                  (sx - sy - k)
2611                                  + qh * power radix (sx - sy - k))
2612                         * vy * power radix k
2613                       +  value (x at StartLoop) (sy + k - 1)
2614                       + power radix (sy + k - 1) * (!x1 at StartLoop)
2615                       - power radix !i * !ql * vy
2616                     = (value (!qp at StartLoop)
2617                                                  (sx - sy - k)
2618                                  + qh * power radix (sx - sy - k))
2619                         * vy * power radix k
2620                       +  value (x at StartLoop) (sy + k - 1)
2621                       + power radix (sy + k - 1) * (!x1 at StartLoop)
2622                     = value (old x) sx
2623                 };
2624         assert { value_sub (pelts x) (!xp.offset + mdn)
2625                  (!xp.offset + mdn + sy - 1)
2626                  + power radix (sy - 1) * !x1
2627                  < vy
2628                  by
2629                  (value xd (sy - 1) + power radix (sy - 1) * !x1 < vy
2630                    by
2631                    value xd (sy - 1) + power radix (sy - 1) * !x1
2632                    = value (xd at Adjust) (sy - 1)
2633                      + power radix (sy - 1) * (!x1 at Adjust)
2634                      + vy
2635                      - power radix sy
2636                    so value (xd at Adjust) (sy - 1)
2637                       < power radix (sy - 1)
2638                    so 1 + (!x1 at Adjust) <= radix
2639                    so value (xd at Adjust) (sy - 1)
2640                       + power radix (sy - 1) * (!x1 at Adjust)
2641                       + vy
2642                       - power radix sy
2643                     < power radix (sy - 1)
2644                       + power radix (sy - 1) * (!x1 at Adjust)
2645                       + vy
2646                       - power radix sy
2647                     = power radix (sy - 1) * (1 + !x1 at Adjust)
2648                       + vy
2649                       - power radix sy
2650                     <= power radix (sy - 1) * radix
2651                       + vy
2652                       - power radix sy
2653                     = vy
2654                  )
2655                  so pelts x = pelts xd
2656                  so xd.offset = !xp.offset + mdn
2657                  so value xd (sy - 1)
2658                     = value_sub (pelts x) (!xp.offset + mdn)
2659                                           (!xp.offset + mdn + sy - 1)
2660                };
2661         assert { dl + radix * dh
2662                  >= (pelts x)[(!xp).offset] + radix * !x1
2663                  by
2664                     vy = vly + power radix (sy - 2)
2665                              * (dl + radix * dh)
2666                  so value_sub (pelts x) (!xp.offset + mdn)
2667                               (!xp.offset + mdn + sy - 1)
2668                     + power radix (sy - 1) * !x1
2669                     < vy
2670                  so !xp.offset + mdn + sy - 1 = !xp.offset + 1
2671                  so power radix (sy - 1) = power radix (sy - 2) * radix
2672                  so - mdn = sy - 2
2673                  so vy
2674                     > value_sub (pelts x) (!xp.offset + mdn)
2675                                           (!xp.offset + mdn + sy - 1)
2676                       + power radix (sy - 1) * !x1
2677                     = value_sub (pelts x) (!xp.offset + mdn)
2678                                           (!xp.offset + 1)
2679                       + power radix (sy - 1) * !x1
2680                     = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2681                       + power radix (- mdn) * (pelts x)[(!xp).offset]
2682                       + power radix (sy - 1) * !x1
2683                     = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2684                       + power radix (sy - 2) * (pelts x)[(!xp).offset]
2685                       + power radix (sy - 1) * !x1
2686                     = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2687                       + power radix (sy - 2) * (pelts x)[(!xp).offset]
2688                       + power radix (sy - 2) * radix * !x1
2689                     = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2690                       + power radix (sy - 2)
2691                         * ((pelts x)[(!xp).offset] + radix * !x1)
2692                     >= power radix (sy - 2)
2693                         * ((pelts x)[(!xp).offset] + radix * !x1)
2694                  so vly < power radix (sy - 2)
2695                  so vy < power radix (sy - 2)
2696                          + power radix (sy - 2)
2697                            * (dl + radix * dh)
2698                          = power radix (sy - 2)
2699                            * (1 + dl + radix * dh)
2700                  so power radix (sy - 2)
2701                     * ((pelts x)[(!xp).offset] + radix * !x1)
2702                     < power radix (sy - 2) * (1 + dl + radix * dh)
2703                  so power radix (sy - 2)
2704                     * ((pelts x)[(!xp).offset] + radix * !x1
2705                       - (1 + dl + radix * dh))
2706                     < 0
2707                  so (pelts x)[(!xp).offset] + radix * !x1
2708                       - (1 + dl + radix * dh)
2709                     < 0
2710                };
2711         end
2712         else begin
2713           qp.contents <- C.incr !qp (-1);
2714           value_sub_update_no_change (pelts q) (!qp).offset
2715                             ((!qp).offset + 1)
2716                             ((!qp).offset + p2i sx - p2i sy - p2i !i)
2717                             !ql;
2718           C.set !qp !ql;
2719           value_sub_head (pelts q) (!qp).offset
2720             ((!qp).offset + p2i sx - p2i sy - p2i !i);
2721           assert { value !qp (sx - sy - !i) * vy
2722                     = !ql * vy + radix *
2723                           (value_sub (pelts q) ((!qp).offset + 1)
2724                                     ((!qp).offset + sx - sy - !i) * vy) }; (*nonlinear*)
2725           assert { value_sub (pelts q) ((!qp).offset + 1)
2726                                     ((!qp).offset + sx - sy - !i) * vy
2727                    =  (value !qp (sx - sy - !i) * vy at StartLoop) }; (*nonlinear*)
2728           value_tail x (sy + !i - 1);
2729           value_sub_concat (pelts x) x.offset xd.offset (x.offset + s);
2730           (* todo refl *)
2731           assert { cy2 = 0 };
2732           assert { value x !i = value (x at SubProd) !i };
2733           assert { value x s = value x !i + power radix !i * value xd (sy-1)
2734                      by xd.offset = x.offset + !i
2735                      so x.offset + s = xd.offset + sy - 1
2736                      so pelts x = pelts xd
2737                      so x.offset + s - xd.offset = sy - 1
2738                      so value_sub (pelts x) xd.offset (x.offset + s)
2739                         = value xd (sy - 1)
2740                      so value x s
2741                         = value_sub (pelts x) x.offset xd.offset
2742                           + power radix !i * value_sub (pelts x) xd.offset (x.offset + s)
2743                         = value x !i
2744                           + power radix !i * value xd (sy - 1)}; (*lifted from assertion*)
2745           assert { (value !qp (sx - sy - !i) + qh * power radix (sx - sy - !i))
2746                    * vy
2747                    = value !qp (sx - sy - !i) * vy
2748                      + qh * vy * power radix (sx - sy - !i)  }; (*nonlinear*)
2749           assert { ((value !qp (sx - sy - !i) + qh * power radix (sx - sy - !i))
2750                    * vy at StartLoop)
2751                    = (value !qp (sx - sy - !i) * vy
2752                      + qh * vy * power radix (sx - sy - !i) at StartLoop)  }; (*nonlinear*)
2753           assert { value x s = value x (sy + !i - 1) };
2754           assert { value (xd at SmallDiv) sy =
2755                     vlx + power radix (sy - 2) * xp0
2756                     + power radix (sy - 1) * xp1 };  (*nonlinear*)
2757           assert { value (x at SubProd) (sy + (!i at StartLoop) - 1)
2758                    = value (x at SubProd) !i + power radix !i * value (xd at SubProd) sy };
2759           assert { value (old x) sx =
2760                   (value !qp (sx - sy - !i)
2761                    + qh * power radix (sx - sy - !i))
2762                   * vy * power radix !i
2763                   + value x (sy + !i - 1)
2764                   + power radix (sy + !i - 1) * !x1 };
2765           value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
2766           value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 2);
2767           let ghost vly = value y (p2i sy - 2) in
2768           assert { vy = vly + power radix (sy - 2) * dl
2769                             + power radix (sy - 1) * dh
2770                    by (pelts y)[y.offset + sy - 1] = dh
2771                    so (pelts y)[y.offset + sy - 2] = dl
2772                    so
2773                       vy = value y (sy - 1)
2774                             + power radix (sy - 1) * dh
2775                          = vly + power radix (sy - 2) * dl
2776                             + power radix (sy - 1) * dh };
2777           assert { value_sub (pelts x) (!xp.offset + mdn)
2778                    (!xp.offset + mdn + sy - 1)
2779                    + power radix (sy - 1) * !x1
2780                    < vy
2781                    by
2782                      pelts x = pelts xd
2783                    so xd.offset = !xp.offset + mdn
2784                    so !xp.offset + mdn + sy - 1 = xd.offset + sy - 1
2785                    so
2786                      value xd (sy - 1)
2787                      = value_sub (pelts xd) xd.offset (xd.offset + sy - 1)
2788                      = value_sub (pelts x) (!xp.offset + mdn)
2789                        (!xp.offset + mdn + sy - 1)
2790                    so value xd (sy - 1)
2791                       + power radix (sy - 1) * !x1
2792                       - power radix sy * cy2
2793                       = value (xd at SubProd) sy
2794                       + power radix sy * (!x1 at StartLoop)
2795                       - !ql * vy
2796                    so cy2 = 0
2797                    so value xd (sy - 1)
2798                       + power radix (sy - 1) * !x1
2799                       = value (xd at SubProd) sy
2800                       + power radix sy * (!x1 at StartLoop)
2801                       - !ql * vy
2802                    so !ql * (dl + radix * dh)
2803                       + (rl + radix * rh)
2804                       = xp0
2805                         + radix * xp1
2806                         + radix * radix * (!x1 at StartLoop)
2807                    so vy = vly + power radix (sy - 2)
2808                                  * (dl + radix * dh)
2809                    so !ql * vy
2810                       = power radix (sy - 2) *
2811                           (xp0
2812                           + radix * xp1
2813                           + radix * radix * (!x1 at StartLoop))
2814                         - power radix (sy - 2) * (rl + radix * rh)
2815                         + !ql * vly
2816                    so value (xd at SubProd) sy
2817                        = vlx
2818                         + power radix (sy - 2) * (xp0 + radix * xp1)
2819                    so power radix sy
2820                       = power radix (sy - 2) * radix * radix
2821                    so (value (xd at SubProd) sy
2822                       + power radix sy * (!x1 at StartLoop)
2823                       - !ql * vy
2824                       < vy
2825                       by
2826                       (!ql * vly >= 0
2827                        by !ql >= 0 so vly >= 0)
2828                       so (power radix (sy - 2) * (rl + radix * rh)
2829                          <= power radix (sy - 2)
2830                             * (dl + radix * dh)
2831                             - power radix (sy - 2)
2832                          by
2833                          rl + radix * rh <= dl + radix * dh - 1
2834                          so power radix (sy - 2) >= 0
2835                          so power radix (sy - 2) * (rl + radix * rh)
2836                             <= power radix (sy - 2)
2837                                * (dl + radix * dh - 1)
2838                             = power radix (sy - 2)
2839                               * (dl + radix * dh)
2840                               - power radix (sy - 2)
2841                          )
2842                       so vlx < power radix (sy - 2)
2843                       so value (xd at SubProd) sy
2844                         + power radix sy * (!x1 at StartLoop)
2845                         - !ql * vy
2846                       = vlx
2847                         + power radix (sy - 2) * (xp0 + radix * xp1)
2848                         + power radix sy * (!x1 at StartLoop)
2849                         - !ql * vy
2850                       = vlx
2851                         + power radix (sy - 2) *
2852                                 (xp0 + radix * xp1
2853                                   + radix * radix * (!x1 at StartLoop))
2854                         - !ql * vy
2855                       = vlx
2856                         + power radix (sy - 2) *
2857                                 (xp0 + radix * xp1
2858                                   + radix * radix * (!x1 at StartLoop))
2859                         -  (power radix (sy - 2) *
2860                            (xp0
2861                             + radix * xp1
2862                             + radix * radix * (!x1 at StartLoop))
2863                           - power radix (sy - 2) * (rl + radix * rh)
2864                           + !ql * vly)
2865                       = vlx
2866                         + power radix (sy - 2) * (rl + radix * rh)
2867                         - !ql * vly
2868                       <= vlx
2869                         + power radix (sy - 2) * (rl + radix * rh)
2870                       <= vlx
2871                          +  power radix (sy - 2)
2872                             * (dl + radix * dh)
2873                          - power radix (sy - 2)
2874                       <  power radix (sy - 2)
2875                          +  power radix (sy - 2)
2876                             * (dl + radix * dh)
2877                          - power radix (sy - 2)
2878                       = power radix (sy - 2) * (dl + radix * dh)
2879                       = vy - vly <= vy
2880                       )
2881                     so value_sub (pelts x) (!xp.offset + mdn)
2882                        (!xp.offset + mdn + sy - 1)
2883                        + power radix (sy - 1) * !x1
2884                        = value xd (sy - 1)
2885                          + power radix (sy - 1) * !x1
2886                        = value (xd at SubProd) sy
2887                         + power radix sy * (!x1 at StartLoop)
2888                         - !ql * vy
2889                        < vy
2890                  };
2891           value_sub_tail (pelts x) (!xp.offset + p2i mdn) (!xp.offset);
2892           value_sub_upper_bound (pelts y) (y.offset) (y.offset + p2i sy - 2);
2893           value_sub_lower_bound (pelts x) (!xp.offset + p2i mdn) (!xp.offset);
2894           assert { dl + radix * dh
2895                    >= (pelts x)[(!xp).offset] + radix * !x1
2896                    by
2897                    vy = vly + power radix (sy - 2)
2898                               * (dl + radix * dh)
2899                    so value_sub (pelts x) (!xp.offset + mdn)
2900                    (!xp.offset + mdn + sy - 1)
2901                    + power radix (sy - 1) * !x1
2902                    < vy
2903                    so !xp.offset + mdn + sy - 1 = !xp.offset + 1
2904                    so power radix (sy - 1) = power radix (sy - 2) * radix
2905                    so - mdn = sy - 2
2906                    so vy
2907                       > value_sub (pelts x) (!xp.offset + mdn)
2908                                             (!xp.offset + mdn + sy - 1)
2909                         + power radix (sy - 1) * !x1
2910                       = value_sub (pelts x) (!xp.offset + mdn)
2911                                             (!xp.offset + 1)
2912                         + power radix (sy - 1) * !x1
2913                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2914                         + power radix (- mdn) * (pelts x)[(!xp).offset]
2915                         + power radix (sy - 1) * !x1
2916                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2917                         + power radix (sy - 2) * (pelts x)[(!xp).offset]
2918                         + power radix (sy - 1) * !x1
2919                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2920                         + power radix (sy - 2) * (pelts x)[(!xp).offset]
2921                         + power radix (sy - 2) * radix * !x1
2922                       = value_sub (pelts x) (!xp.offset + mdn) (!xp.offset)
2923                         + power radix (sy - 2)
2924                           * ((pelts x)[(!xp).offset] + radix * !x1)
2925                       >= power radix (sy - 2)
2926                           * ((pelts x)[(!xp).offset] + radix * !x1)
2927                    so vly < power radix (sy - 2)
2928                    so vy < power radix (sy - 2)
2929                            + power radix (sy - 2)
2930                               * (dl + radix * dh)
2931                            = power radix (sy - 2)
2932                              * (1 + dl + radix * dh)
2933                    so power radix (sy - 2)
2934                       * ((pelts x)[(!xp).offset] + radix * !x1)
2935                       < power radix (sy - 2) * (1 + dl + radix * dh)
2936                    so power radix (sy - 2)
2937                       * ((pelts x)[(!xp).offset] + radix * !x1
2938                         - (1 + dl + radix * dh))
2939                       < 0
2940                    so (pelts x)[(!xp).offset] + radix * !x1
2941                         - (1 + dl + radix * dh)
2942                       < 0
2943                  };
2944         end;
2945       end;
2946     done;
2947     label EndLoop in
2948     assert { !i = 0 };
2949     assert { !xp.offset = x.offset + sy - 2 };
2950     value_sub_update_no_change (pelts x) (!xp.offset + 1)
2951                                x.offset (!xp.offset) !x1;
2952     C.set_ofs !xp 1 !x1;
2953     assert { value x (sy - 1) =
2954              value (x at EndLoop) (sy - 1)
2955              by pelts x = Map.set (pelts x at EndLoop) (x.offset + sy - 1) !x1 };
2956     value_sub_tail (pelts x) x.offset (!xp.offset+1);
2957     (* todo refl *)
2958     assert { value (old x) sx =
2959               (value q (sx - sy)
2960               + power radix (sx - sy) * qh)
2961                 * value y sy
2962               + value x sy
2963              by
2964                value x sy
2965                 = value x (sy - 1)
2966                   + power radix (sy - 1) * !x1
2967              so vy = value y sy
2968              so value (old x) sx
2969                 = (value !qp (sx - sy - !i)
2970                   + qh * power radix (sx - sy - !i))
2971                   * vy * power radix !i
2972                   + value x (sy + !i - 1)
2973                   + power radix (sy + !i - 1) * !x1
2974                 = (value !qp (sx - sy)
2975                   + qh * power radix (sx - sy))
2976                   * vy * 1
2977                   + value x (sy - 1)
2978                   + power radix (sy - 1) * !x1
2979                 = (value !qp (sx - sy)
2980                   + qh * power radix (sx - sy))
2981                   * value y sy
2982                   + value x sy };
2983     qh
2984
2985   let wmpn_divrem_2 (q x y:t) (sx:int32) : limb
2986     requires { 2 <= sx }
2987     requires { valid x sx }
2988     requires { valid y 2 }
2989     requires { valid q (sx - 2) }
2990     requires { (pelts y)[y.offset + 1] >= div radix 2 }
2991     requires { writable q /\ writable x }
2992     ensures { value (old x) sx =
2993               (value q (sx - 2)
2994                + power radix (sx - 2) * result)
2995               * value y 2
2996               + value x 2 }
2997     ensures { value x 2 < value y 2 }
2998     ensures { 0 <= result <= 1 }
2999   =
3000     let xp = ref (C.incr x (sx - 2)) in
3001     let dh = C.get_ofs y 1 in
3002     let dl = C.get y in
3003     let rh = ref (C.get_ofs !xp 1) in
3004     let rl = ref (C.get !xp) in
3005     let qh = ref 0 in
3006     let lx = ref 0 in
3007     assert { value y 2 = dl + radix * dh };
3008     let i = ref (sx - 2) in
3009     let dinv = reciprocal_word_3by2 dh dl in
3010     ([@vc:sp] if (!rh >= dh && ([@vc:sp] !rh > dh || !rl >= dl))
3011     then
3012       label Adjust in
3013       begin
3014         ensures { value x sx
3015                   = (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3016                       + !qh * power radix (sx - 2 - !i))
3017                     * value y 2 * power radix !i
3018                     + value x !i
3019                     + power radix !i * (!rl + radix * !rh) }
3020         ensures { !rl + radix * !rh < dl + radix * dh }
3021         ensures { !qh = 1 }
3022         let (r0, b) = sub_with_borrow !rl dl 0 in
3023         let (r1, ghost b') = sub_with_borrow !rh dh b in
3024         assert { b' = 0 };
3025         assert { r0 + radix * r1 = !rl + radix * !rh - (dl + radix * dh) };
3026         value_sub_tail (pelts x) x.offset (x.offset + p2i sx - 1);
3027         value_sub_tail (pelts x) x.offset (x.offset + p2i sx - 2);
3028         rh := r1;
3029         rl := r0;
3030         qh := 1;
3031         assert { value x sx
3032                   = (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3033                       + !qh * power radix (sx - 2 - !i))
3034                     * value y 2 * power radix !i
3035                     + value x !i
3036                     + power radix !i * (!rl + radix * !rh)
3037                  by
3038                  value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2) = 0
3039                  so (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3040                       + !qh * power radix (sx - 2 - !i))
3041                     * value y 2 * power radix !i
3042                     + value x !i
3043                     + power radix !i * (!rl + radix * !rh)
3044                   = value y 2 * power radix !i
3045                     + value x !i
3046                     + power radix !i * (!rl + radix * !rh)
3047                   = value x !i
3048                     + power radix !i * (dl + radix * dh + !rl + radix * !rh)
3049                   = value x !i
3050                     + power radix !i * (!rl at Adjust + radix * !rh at Adjust)
3051                   = value x !i
3052                     + power radix !i * !rl at Adjust
3053                     + power radix (!i+1) * !rh at Adjust
3054                   = value x sx
3055                };
3056       end
3057     else
3058       begin
3059         ensures { value x sx
3060                   = (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3061                       + !qh * power radix (sx - 2 - !i))
3062                     * value y 2 * power radix !i
3063                     + value x !i
3064                     + power radix !i * (!rl + radix * !rh) }
3065         ensures { !rl + radix * !rh < dl + radix * dh }
3066         ensures { !qh = 0 }
3067         value_sub_tail (pelts x) x.offset (x.offset + p2i sx - 1);
3068         value_sub_tail (pelts x) x.offset (x.offset + p2i sx - 2);
3069       end);
3070     while (!i > 0) do
3071       variant { !i }
3072       invariant { 0 <= !i <= sx - 2 }
3073       invariant { !xp.offset = x.offset + !i }
3074       invariant { plength !xp = plength x }
3075       invariant { !xp.min = x.min }
3076       invariant { !xp.max = x.max }
3077       invariant { pelts !xp = pelts x }
3078       invariant { writable !xp }
3079       invariant { value x sx
3080                   = (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3081                       + !qh * power radix (sx - 2 - !i))
3082                     * value y 2 * power radix !i
3083                     + value x !i
3084                     + power radix !i * (!rl + radix * !rh) }
3085       invariant { !rl + radix * !rh < dl + radix * dh }
3086       label StartLoop in
3087       let ghost k = p2i !i in
3088       xp.contents <- C.incr !xp (-1);
3089       lx := C.get !xp;
3090       label Got in
3091       let (qu, r0, r1) = div3by2_inv !rh !rl !lx dh dl dinv in
3092       rh := r1;
3093       rl := r0;
3094       i := !i - 1;
3095       C.set_ofs q !i qu;
3096       assert { qu * (dl + radix * dh) + r0 + radix * r1
3097                = !lx + radix * (!rl at StartLoop)
3098                      + radix * radix * (!rh at StartLoop)
3099                by
3100                  radix * ((!rl at StartLoop) + radix * (!rh at StartLoop))
3101                  = radix * (!rl at StartLoop) + radix * radix * (!rh at StartLoop)
3102                so
3103                  qu * (dl + radix * dh) + r0 + radix * r1
3104                = !lx + radix * ((!rl at StartLoop) + radix * (!rh at StartLoop))
3105                = !lx + radix * (!rl at StartLoop)
3106                      + radix * radix * (!rh at StartLoop)
3107              };
3108       value_sub_head (pelts q) (q.offset + p2i !i) (q.offset + p2i sx - 2);
3109       value_sub_tail (pelts x) x.offset (x.offset + p2i !i);
3110       assert { value x sx
3111                   = (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3112                       + !qh * power radix (sx - 2 - !i))
3113                     * value y 2 * power radix !i
3114                     + value x !i
3115                     + power radix !i * (!rl + radix * !rh)
3116                by
3117                   value x k = value x !i + power radix !i * !lx
3118                so value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3119                  = qu + radix
3120                         * value_sub (pelts q) (q.offset + k) (q.offset + sx - 2)
3121                so power radix (sx - 2 - !i) = radix * power radix (sx - 2 - k)
3122                so
3123                  (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3124                         + !qh * power radix (sx - 2 - !i))
3125                  = qu + radix
3126                         * (value_sub (pelts q) (q.offset + k) (q.offset + sx - 2)
3127                            + !qh * power radix (sx - 2 - k))
3128                so power radix !i * radix  = power radix k
3129                so ((value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3130                       + !qh * power radix (sx - 2 - !i))
3131                     * value y 2 * power radix !i
3132                   = power radix !i * qu * (dl + radix * dh)
3133                     + (value_sub (pelts q) (q.offset + k)
3134                                                    (q.offset + sx - 2)
3135                                + !qh * power radix (sx - 2 - k))
3136                       * value y 2 * power radix k
3137                   by
3138                   (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3139                       + !qh * power radix (sx - 2 - !i))
3140                     * value y 2 * power radix !i
3141                   = (qu + radix
3142                          * (value_sub (pelts q) (q.offset + k)
3143                                                 (q.offset + sx - 2)
3144                              + !qh * power radix (sx - 2 - k)))
3145                     * value y 2 * power radix !i
3146                   = power radix !i * qu * (dl + radix * dh)
3147                     + radix * (value_sub (pelts q) (q.offset + k)
3148                                                    (q.offset + sx - 2)
3149                                + !qh * power radix (sx - 2 - k))
3150                       * value y 2 * power radix !i
3151                   = power radix !i * qu * (dl + radix * dh)
3152                     + (value_sub (pelts q) (q.offset + k)
3153                                                    (q.offset + sx - 2)
3154                                + !qh * power radix (sx - 2 - k))
3155                       * value y 2 * power radix k)
3156                so (value_sub (pelts q) (q.offset + !i) (q.offset + sx - 2)
3157                       + !qh * power radix (sx - 2 - !i))
3158                     * value y 2 * power radix !i
3159                     + value x !i
3160                     + power radix !i * (!rl + radix * !rh)
3161                   = power radix !i * qu * (dl + radix * dh)
3162                     + (value_sub (pelts q) (q.offset + k)
3163                                            (q.offset + sx - 2)
3164                                + !qh * power radix (sx - 2 - k))
3165                       * value y 2 * power radix k
3166                     + value x !i
3167                     + power radix !i * (!rl + radix * !rh)
3168                   = (value_sub (pelts q) (q.offset + k)
3169                                            (q.offset + sx - 2)
3170                                + !qh * power radix (sx - 2 - k))
3171                       * value y 2 * power radix k
3172                     + value x !i
3173                     + power radix !i * (qu * (dl + radix * dh)
3174                                         + !rl + radix * !rh)
3175                   = (value_sub (pelts q) (q.offset + k)
3176                                            (q.offset + sx - 2)
3177                                + !qh * power radix (sx - 2 - k))
3178                       * value y 2 * power radix k
3179                     + value x !i
3180                     + power radix !i
3181                       * (!lx + radix * (!rl at StartLoop)
3182                          + radix * radix * (!rh at StartLoop))
3183                   = (value_sub (pelts q) (q.offset + k)
3184                                            (q.offset + sx - 2)
3185                                + !qh * power radix (sx - 2 - k))
3186                       * value y 2 * power radix k
3187                     + value x !i
3188                     + power radix !i * !lx
3189                     + power radix !i * (radix * (!rl at StartLoop
3190                                                  + radix * !rh at StartLoop))
3191                   = (value_sub (pelts q) (q.offset + k)
3192                                            (q.offset + sx - 2)
3193                                + !qh * power radix (sx - 2 - k))
3194                       * value y 2 * power radix k
3195                     + value x k
3196                     + power radix !i * (radix * (!rl at StartLoop
3197                                                  + radix * !rh at StartLoop))
3198                   = (value_sub (pelts q) (q.offset + k)
3199                                            (q.offset + sx - 2)
3200                                + !qh * power radix (sx - 2 - k))
3201                       * value y 2 * power radix k
3202                     + value x k
3203                     + power radix k * (!rl at StartLoop
3204                                        + radix * !rh at StartLoop)
3205                   = value x sx
3206           };
3207     done;
3208     assert { !i = 0 };
3209     assert { value x sx
3210              = (value_sub (pelts q) q.offset (q.offset + sx - 2)
3211                       + !qh * power radix (sx - 2))
3212                     * value y 2
3213                     + !rl + radix * !rh
3214              by power radix !i = 1 };
3215     C.set_ofs x 1 !rh;
3216     C.set x !rl;
3217     assert { value x 2 = !rl + radix * !rh
3218              by (pelts x)[x.offset] = !rl
3219                 /\ (pelts x)[x.offset + 1] = !rh};
3220     !qh
3221
3222   (** `div_qr q r x y sx sy` divides `(x, sx)` by `(y, sy)`, writes the quotient
3223     in `(q, (sx-sy))` and the remainder in `(r, sy)`. Corresponds to
3224     `mpn_tdiv_qr`. *)
3225   let div_qr (q r x y nx ny:t) (sx sy:int32) : unit
3226     requires { 1 <= sy <= sx <= (Int32.max_int32 - 1) }
3227     requires { valid x sx }
3228     requires { valid y sy }
3229     requires { valid q (sx - sy + 1) }
3230     requires { valid r sy }
3231     requires { valid nx (sx + 1) }
3232     requires { valid ny sy }
3233     requires { writable nx /\ writable ny }
3234     requires { (pelts y)[y.offset + sy - 1] > 0 }
3235     requires { writable q /\ writable r }
3236     ensures { value x sx
3237               = value q (sx - sy + 1) * value y sy
3238                 + value r sy }
3239     ensures { value r sy < value y sy }
3240   =
3241     label Start in
3242     value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
3243     value_sub_lower_bound (pelts y) y.offset (y.offset + p2i sy - 1);
3244     assert { value y sy >= power radix (sy - 1) };
3245     if (sy = 1)
3246     then
3247       let lr = wmpn_divrem_1 q x sx (C.get y) in
3248       C.set r lr
3249     else
3250     if (sy = 2)
3251     then
3252       let clz = clz_ext (C.get_ofs y (sy - 1)) in
3253       let ghost p = power 2 (p2i clz) in
3254       if clz = 0
3255       then begin
3256         wmpn_copyi nx x sx;
3257         value_sub_shift_no_change (pelts x) x.offset (p2i sx) (p2i sx) 0;
3258         C.set_ofs nx sx 0;
3259         value_sub_frame_shift (pelts x) (pelts nx) x.offset nx.offset (p2i sx);
3260         label Div2_ns in
3261         let ghost _qh = wmpn_divrem_2 q nx y (sx + 1) in
3262         wmpn_copyi r nx sy;
3263         assert { value x sx
3264                  = value q (sx - sy + 1) * value y sy
3265                  + value r sy
3266                by value r sy = value nx sy
3267                so value (nx at Div2_ns) (sx + 1) < power radix sx
3268                so value (nx at Div2_ns) (sx + 1)
3269                   = value (nx at Div2_ns) sx
3270                so (_qh = 0
3271                    by
3272                    power radix sx
3273                    > value (nx at Div2_ns) (sx + 1)
3274                    = (value q (sx - 1) + power radix (sx - 1) * _qh)
3275                     * value y 2
3276                     + value nx 2
3277                    so value nx 2 >= 0
3278                    so value y 2 >= radix
3279                    so value q (sx - 1) >= 0
3280                    so _qh >= 0
3281                    so (value q (sx - 1)
3282                            + power radix (sx - 1) * _qh) >= 0
3283                    so (value q (sx - 1) + power radix (sx - 1) * _qh)
3284                       * value y 2
3285                       + value nx 2
3286                       >= (value q (sx - 1)
3287                            + power radix (sx - 1) * _qh)
3288                          * value y 2
3289                       >= (value q (sx - 1)
3290                            + power radix (sx - 1) * _qh)
3291                          * radix
3292                       >= power radix (sx - 1) * _qh * radix
3293                       = power radix sx * _qh
3294                    so power radix sx > power radix sx * _qh
3295                   )
3296                so value x sx = value (nx at Div2_ns) sx
3297                };
3298         ()
3299       end
3300       else begin
3301         let ghost _c = wmpn_lshift ny y sy (Limb.of_int32 clz) in
3302         begin
3303           ensures { normalized ny sy }
3304           ensures { value ny sy = power 2 clz * value y sy }
3305           let ghost dh = (pelts y)[y.offset + p2i sy - 1] in
3306           assert { value y sy
3307                    = value y (sy - 1) + power radix (sy - 1) * dh };
3308           value_sub_upper_bound (pelts y) y.offset (y.offset + p2i sy - 1);
3309           value_sub_tail (pelts ny) ny.offset (ny.offset + p2i sy - 1);
3310           value_sub_upper_bound (pelts ny) ny.offset (ny.offset + p2i sy - 1);
3311           let ghost ndh = (pelts ny)[ny.offset + p2i sy - 1] in
3312           assert { normalized ny sy
3313                    /\ value ny sy = power 2 clz * value y sy
3314                    by
3315                    value y sy < (dh + 1) * power radix (sy - 1)
3316                    so value ny sy + (power radix sy) * _c
3317                       = power 2 clz * value y sy
3318                       = power 2 clz
3319                         * (value y (sy - 1)
3320                            + dh * power radix (sy - 1))
3321                    so power 2 clz * dh <= radix - power 2 clz
3322                    so value ny sy + (power radix sy) * _c
3323                       = power 2 clz * value y (sy - 1)
3324                         + power 2 clz * dh * power radix (sy - 1)
3325                       < power 2 clz * power radix (sy - 1)
3326                         + power 2 clz * dh * power radix (sy - 1)
3327                       <= power 2 clz * power radix (sy - 1)
3328                         + (radix - power 2 clz) * power radix (sy - 1)
3329                        = radix * power radix (sy - 1)
3330                        = power radix sy
3331                    so _c = 0
3332                    so value ny sy
3333                       = power 2 clz * value y sy
3334                    so value y sy >= dh * power radix (sy - 1)
3335                    so value ny sy
3336                       >= power 2  clz * dh * power radix (sy - 1)
3337                    so value ny sy =
3338                       value ny (sy - 1) + power radix (sy - 1) * ndh
3339                       < power radix (sy - 1) + power radix (sy - 1) * ndh
3340                       = power radix (sy - 1) * (ndh + 1)
3341                    so power radix (sy - 1) * (ndh + 1)
3342                       > power radix (sy - 1) * (power 2 clz * dh)
3343                    so ndh + 1 > power 2 clz * dh
3344                    so ndh >= power 2 clz * dh
3345                    so 2 * power 2 clz * dh >= radix
3346                    so 2 * ndh >= radix
3347                    so ndh >= div radix 2
3348                  };
3349         end;
3350         let h = wmpn_lshift nx x sx (Limb.of_int32 clz) in
3351         C.set_ofs nx sx h;
3352         begin
3353           ensures { value nx (sx + 1)
3354                     = p * value x sx  }
3355           value_sub_tail (pelts nx) nx.offset (nx.offset + p2i sx);
3356           assert { value nx (sx + 1)
3357                     = p * value x sx
3358                    by
3359                    value nx sx + power radix sx * h
3360                     = p * value x sx
3361                    so value nx (sx + 1)
3362                       = value nx sx + power radix sx * h
3363                  }
3364         end;
3365         label Div2_s in
3366         (* TODO don't add 1 when not needed, cf "adjust" in GMP algo *)
3367         let ghost _qh = wmpn_divrem_2 q nx ny (sx + 1) in
3368         let ghost _l = wmpn_rshift r nx sy (Limb.of_int32 clz) in
3369         begin ensures { value nx 2 = p * value r 2 }
3370           assert { _l = 0
3371                  by
3372                    (mod (value nx sy) p = 0
3373                     by
3374                       value (nx at Div2_s) (sx + 1)
3375                    = (value q (sx - 1) + power radix (sx - 1) * _qh)
3376                       * value ny sy
3377                       + value nx sy
3378                    so value (nx at Div2_s) (sx + 1)
3379                     = p * value x sx
3380                    so value ny sy = p * value y sy
3381                    so value nx sy
3382                       = value (nx at Div2_s) (sx + 1)
3383                         - (value q (sx - 1)
3384                           + power radix (sx - 1) * _qh)
3385                           * value ny sy
3386                       = p * value x sx
3387                         - p * (value q (sx - 1)
3388                                + power radix (sx - 1) * _qh)
3389                             * value y sy
3390                       = p * (value x sx
3391                              - (value q (sx - 1)
3392                                + power radix (sx - 1) * _qh)
3393                                * value y sy)
3394                    so let n = (value x sx
3395                                 - (value q (sx - 1)
3396                                    + power radix (sx - 1) * _qh)
3397                                    * value y sy)
3398                       in
3399                       value nx sy = p * n
3400                       so value nx sy >= 0
3401                       so p > 0
3402                       so n >= 0
3403                       so mod (value nx sy) p
3404                          = mod (p * n) p
3405                          = mod ((p*n)+0) p
3406                          = mod 0 p
3407                          = 0
3408                  )
3409                  so _l + radix * value r sy
3410                     = power 2 (Limb.length - clz) * (value nx sy)
3411                  so let a = div (value nx sy) p in
3412                     value nx sy = p * a
3413                  so power 2 (Limb.length - clz) * p = radix
3414                  so power 2 (Limb.length - clz) * (value nx sy)
3415                     = power 2 (Limb.length - clz) * (p * a)
3416                     = (power 2 (Limb.length - clz) * p) * a
3417                     = radix * a
3418                  so mod (radix * value r sy + _l) radix
3419                     = mod _l radix
3420                  so mod (radix * value r sy + _l) radix
3421                     = mod (radix * a) radix = 0
3422                  so mod _l radix = 0
3423                  so 0 <= _l < radix
3424                };
3425           assert { value nx 2 = p * value r 2
3426                    by
3427                    radix * value r 2
3428                    = power 2 (Limb.length - clz) * value nx 2
3429                    so p * power 2 (Limb.length - clz)
3430                       = radix
3431                    so p * radix * value r 2
3432                       = p * power 2 (Limb.length - clz) * value nx 2
3433                       = radix * value nx 2
3434                    so p * value r 2 = value nx 2
3435                  }
3436         end;
3437         assert { value x sx
3438                  = value q (sx - sy + 1) * value y sy
3439                  + value r sy
3440                 by
3441                  value (nx at Div2_s) (sx + 1)
3442                  = (value q (sx - 1) + power radix (sx - 1) * _qh)
3443                    * value ny 2
3444                    + value nx 2
3445                 so value (nx at Div2_s) (sx + 1)
3446                     = p * value x sx
3447                 so value ny 2 = p * value y 2
3448                 so (_qh = 0
3449                     by
3450                     value x sx < power radix sx
3451                     so value y 2 >= radix
3452                     so value ny 2 >= p * radix
3453                     so value q (sx - 1) >= 0
3454                     so value nx 2 >= 0
3455                     so (value q (sx - 1) + power radix (sx - 1) * _qh)
3456                         >= 0
3457                     so (value q (sx - 1) + power radix (sx - 1) * _qh)
3458                          * value ny 2
3459                        + value nx 2
3460                        >= (value q (sx - 1)
3461                            + power radix (sx - 1) * _qh)
3462                           * value ny 2
3463                        >= (value q (sx - 1)
3464                            + power radix (sx - 1) * _qh)
3465                           * (p * radix)
3466                        >= power radix (sx - 1) * _qh * p * radix
3467                        = power radix sx * p * _qh
3468                     so power radix sx * p
3469                        > value (nx at Div2_s) (sx + 1)
3470                        >= power radix sx * p * _qh
3471                     )
3472                 so value nx 2 = p * value r 2
3473                 so p * value x sx
3474                    = value q (sx - 1) * p * value y 2
3475                      + p * value r 2
3476                    = p * (value q (sx - 1) * value y 2
3477                           + value r 2)
3478                };
3479         ()
3480       end
3481     else
3482      (* let qn = ref (Int32.(-) (Int32.(+) sx 1) sy) in
3483       if (Int32.(>=)  (Int32.(+) !qn !qn) sx)
3484       then*) begin
3485         let adjust =
3486           if (get_ofs x (sx - 1)) >= (get_ofs y (sy - 1))
3487           then 1
3488           else 0
3489         in
3490         let clz = clz_ext (C.get_ofs y (sy - 1)) in
3491         let ghost p = power 2 (p2i clz) in
3492         if clz = 0
3493         then begin
3494           wmpn_copyi nx x sx;
3495           value_sub_shift_no_change (pelts x) x.offset
3496                                     (p2i sx) (p2i sx) 0;
3497           C.set_ofs nx sx 0;
3498           value_sub_frame_shift (pelts x) (pelts nx) x.offset nx.offset (p2i sx);
3499           assert { value y sy * (power radix (sx - sy + adjust))
3500                    > value nx (sx + adjust)
3501                    by
3502                    let dh = (pelts y)[y.offset + sy - 1] in
3503                    value y sy >= dh * power radix (sy - 1)
3504                    so value nx (sx + adjust) = value nx sx = value x sx
3505                    so [@case_split]
3506                       ((adjust = 1
3507                            so value x sx < power radix sx
3508                            so value y sy * power radix (sx - sy + adjust)
3509                               >= dh * power radix (sy - 1)
3510                                     * power radix (sx - sy + adjust)
3511                               = dh * power radix ((sy - 1) + (sx - sy + adjust))
3512                               = dh * power radix sx
3513                            so dh >= div radix 2 > 1
3514                            so dh * power radix sx > power radix sx )
3515                       \/
3516                       (adjust = 0
3517                            so let ah = (pelts x)[x.offset + sx - 1] in
3518                               value x sx < (ah + 1) * power radix (sx - 1)
3519                               so ah + 1 <= dh
3520                               so value x sx < dh * power radix (sx - 1)
3521                               so value y sy * power radix (sx - sy + adjust)
3522                                  = value y sy * power radix (sx - sy)
3523                                  >= dh * power radix (sy - 1)
3524                                        * power radix (sx - sy)
3525                                  = dh * power radix (sy - 1 + sx - sy)
3526                                  = dh * power radix (sx - 1))) };
3527           label Div_ns in
3528           let ghost _qh = div_sb_qr q nx (sx + adjust) y sy in
3529           wmpn_copyi r nx sy;
3530           assert { value x sx
3531                    = value q (sx - sy + adjust) * value y sy
3532                    + value r sy
3533                  by value r sy = value nx sy
3534                  so value (nx at Div_ns) (sx + adjust) = value x sx < power radix sx
3535                  so value (nx at Div_ns) (sx + adjust)
3536                     = value (nx at Div_ns) sx
3537                  so (_qh = 0
3538                      by
3539                      value (nx at Div_ns) (sx + adjust)
3540                      = (value q (sx - sy + adjust)
3541                         + power radix (sx - sy + adjust) * _qh)
3542                       * value y sy
3543                       + value nx sy
3544                      so value nx sy >= 0
3545                      so value q (sx - sy + adjust) >= 0
3546                      so _qh >= 0
3547                      so (value q (sx - sy + adjust)
3548                              + power radix (sx - sy + adjust) * _qh) >= 0
3549                      so (value q (sx - sy + adjust)
3550                          + power radix (sx - sy + adjust) * _qh)
3551                         * value y sy
3552                         + value nx sy
3553                         >= (value q (sx - sy + adjust)
3554                              + power radix (sx - sy + adjust) * _qh)
3555                            * value y sy
3556                         >= power radix (sx - sy + adjust) * _qh * value y sy
3557                      so _qh <> 1)
3558                  so value x sx = value (nx at Div_ns) sx
3559                  };
3560            label Ret_ns in
3561            begin
3562              ensures { value q (sx - sy + 1)
3563                        = value (q at Ret_ns) (sx - sy + adjust) }
3564              if (adjust = 0)
3565              then begin
3566                value_sub_shift_no_change (pelts x) x.offset
3567                                          (p2i sx) (p2i sx) 0;
3568                set_ofs q (sx - sy) 0;
3569                value_sub_tail (pelts q) q.offset (q.offset + p2i sx - p2i sy);
3570                ()
3571              end
3572            end
3573         end
3574         else begin
3575           let ghost _c = wmpn_lshift ny y sy (Limb.of_int32 clz) in
3576           begin
3577             ensures { normalized ny sy }
3578             ensures { value ny sy
3579                       = power 2 clz * value y sy }
3580             let ghost dh = (pelts y)[y.offset + p2i sy - 1] in
3581             assert { value y sy
3582                      = value y (sy - 1) + power radix (sy - 1) * dh };
3583             value_sub_upper_bound (pelts y) y.offset (y.offset + p2i sy - 1);
3584             value_sub_tail (pelts ny) ny.offset (ny.offset + p2i sy - 1);
3585             value_sub_upper_bound (pelts ny) ny.offset (ny.offset + p2i sy - 1);
3586             let ghost ndh = (pelts ny)[ny.offset + p2i sy - 1] in
3587             assert { normalized ny sy
3588                      /\ value ny sy
3589                         = power 2 clz * value y sy
3590                      by
3591                      value y sy < (dh + 1) * power radix (sy - 1)
3592                      so value ny sy + (power radix sy) * _c
3593                         = power 2 clz * value y sy
3594                         = power 2 clz
3595                           * (value y (sy - 1)
3596                              + dh * power radix (sy - 1))
3597                      so power 2 clz * dh <= radix - power 2 clz
3598                      so (_c = 0
3599                         by
3600                           value ny sy + (power radix sy) * _c
3601                           = power 2 clz * value y (sy - 1)
3602                             + power 2 clz * dh * power radix (sy - 1)
3603                           < power 2 clz * power radix (sy - 1)
3604                             + power 2 clz * dh * power radix (sy - 1)
3605                           <= power 2 clz * power radix (sy - 1)
3606                             + (radix - power 2 clz) * power radix (sy - 1)
3607                            = radix * power radix (sy - 1)
3608                            = power radix sy
3609                         so value ny sy >= 0
3610                         so power radix sy * _c < power radix sy
3611                         so power radix sy > 0
3612                         so _c >= 0
3613                         )
3614                      so value ny sy
3615                         = power 2 clz * value y sy
3616                      so value y sy >= dh * power radix (sy - 1)
3617                      so value ny sy
3618                         >= power 2  clz * dh * power radix (sy - 1)
3619                      so value ny sy =
3620                         value ny (sy - 1) + power radix (sy - 1) * ndh
3621                         < power radix (sy - 1) + power radix (sy - 1) * ndh
3622                         = power radix (sy - 1) * (ndh + 1)
3623                      so power radix (sy - 1) * (ndh + 1)
3624                         > power radix (sy - 1) * (power 2 clz * dh)
3625                      so ndh + 1 > power 2 clz * dh
3626                      so ndh >= power 2 clz * dh
3627                      so 2 * power 2 clz * dh >= radix
3628                      so 2 * ndh >= radix
3629                      so ndh >= div radix 2
3630                    };
3631           end;
3632           let h = wmpn_lshift nx x sx (Limb.of_int32 clz) in
3633           label Shifted in
3634           C.set_ofs nx sx h;
3635           begin
3636             ensures { value nx (sx + adjust)
3637                       = p * value x sx  }
3638             if (adjust = 1)
3639             then begin
3640               value_sub_tail (pelts nx) nx.offset (nx.offset + p2i sx);
3641               assert { value nx (sx + 1)
3642                         = p * value x sx
3643                        by
3644                        value nx sx + power radix sx * h
3645                         = p * value x sx
3646                        so value nx (sx + 1)
3647                           = value nx sx + power radix sx * h
3648                      } end
3649             else begin
3650               assert { adjust = 0 };
3651               assert { h = 0
3652                        by
3653                        let dh = (pelts y)[y.offset + sy - 1] in
3654                        let ah = (pelts x)[x.offset + sx - 1] in
3655                        p * dh < radix
3656                        so 0 <= ah <= dh
3657                        so p * ah < radix
3658                        so (p * ah <= radix - p
3659                            by
3660                            let q = power 2 (Limb.length - clz) in
3661                            radix = p * q
3662                            so p * ah < p * q
3663                            so ah < q
3664                            so ah <= q - 1
3665                            so p * ah <= p * (q - 1) = radix - p
3666                            )
3667                        so p * (ah + 1) <= radix
3668                        so let s = power radix (sx - 1) in
3669                           value x sx < (ah + 1) * s
3670                        so p * value x sx < p * (ah + 1) * s
3671                        so (p * (ah + 1) * s
3672                            <= radix * s
3673                            by
3674                            [@case_split]
3675                            (p * (ah + 1) = radix
3676                            \/ (p * (ah + 1) < radix
3677                                so s > 0
3678                                so p * (ah + 1) * s
3679                                  < radix * s)))
3680                        so radix * power radix (sx - 1) = power radix sx
3681                        so value (nx at Shifted) sx + power radix sx * h
3682                           < power radix sx
3683                        so power radix sx * h < power radix sx * 1
3684                        so (h < 1 by power radix sx > 0)
3685                       }
3686             end
3687           end;
3688           label Div_s in
3689           assert { value ny sy * (power radix (sx - sy + adjust))
3690                    > value nx (sx + adjust)
3691                    by
3692                    let dh = (pelts y)[y.offset + sy - 1] in
3693                    value ny sy >= p * dh * power radix (sy - 1)
3694                    so value nx (sx + adjust) = p * value x sx
3695                    so p > 0
3696                    so [@case_split]
3697                       ((adjust = 1
3698                            so value x sx < power radix sx
3699                            so p * value x sx < p * power radix sx
3700                            so value ny sy * power radix (sx - sy + adjust)
3701                               >= p * dh * power radix (sy - 1)
3702                                     * power radix (sx - sy + adjust)
3703                               = p * dh * power radix ((sy - 1) + (sx - sy + adjust))
3704                               = p * dh * power radix sx
3705                            so dh >= 1
3706                            so p * dh * power radix sx >= p * power radix sx )
3707                       \/
3708                       (adjust = 0
3709                            so let ah = (pelts x)[x.offset + sx - 1] in
3710                               value x sx < (ah + 1) * power radix (sx - 1)
3711                               so ah + 1 <= dh
3712                               so value x sx < dh * power radix (sx - 1)
3713                               so p * value x sx < p * dh * power radix (sx - 1)
3714                               so value ny sy * power radix (sx - sy + adjust)
3715                                  = value ny sy * power radix (sx - sy)
3716                                  >= p * dh * power radix (sy - 1)
3717                                        * power radix (sx - sy)
3718                                  = p * dh * power radix (sy - 1 + sx - sy)
3719                                  = p * dh * power radix (sx - 1))) };
3720           let ghost _qh = div_sb_qr q nx (sx + adjust) ny sy in
3721           let ghost _l = wmpn_rshift r nx sy (Limb.of_int32 clz) in
3722           begin ensures { value nx sy = p * value r sy }
3723             assert { _l = 0
3724                    by
3725                      (mod (value nx sy) p = 0
3726                       by
3727                         value (nx at Div_s) (sx + adjust)
3728                      = (value q (sx - sy + adjust)
3729                          + power radix (sx - sy + adjust) * _qh)
3730                         * value ny sy
3731                         + value nx sy
3732                      so value (nx at Div_s) (sx + adjust)
3733                       = p * value x sx
3734                      so value ny sy = p * value y sy
3735                      so value nx sy
3736                         = value (nx at Div_s) (sx + adjust)
3737                           - (value q (sx - sy + adjust)
3738                             + power radix (sx - sy + adjust) * _qh)
3739                             * value ny sy
3740                         = p * value x sx
3741                           - p * (value q (sx - sy + adjust)
3742                                  + power radix (sx - sy + adjust) * _qh)
3743                               * value y sy
3744                         = p * (value x sx
3745                                - (value q (sx - sy + adjust)
3746                                  + power radix (sx - sy + adjust) * _qh)
3747                                  * value y sy)
3748                      so let n = (value x sx
3749                                   - (value q (sx - sy + adjust)
3750                                      + power radix (sx - sy + adjust) * _qh)
3751                                      * value y sy)
3752                         in
3753                         value nx sy = p * n
3754                         so value nx sy >= 0
3755                         so p > 0
3756                         so n >= 0
3757                         so mod (value nx sy) p
3758                            = mod (p * n) p
3759                            = mod ((p*n)+0) p
3760                            = mod 0 p
3761                            = 0
3762                    )
3763                    so _l + radix * value r sy
3764                       = power 2 (Limb.length - clz) * (value nx sy)
3765                    so let a = div (value nx sy) p in
3766                       value nx sy = p * a
3767                    so power 2 (Limb.length - clz) * p = radix
3768                    so power 2 (Limb.length - clz) * (value nx sy)
3769                       = power 2 (Limb.length - clz) * (p * a)
3770                       = (power 2 (Limb.length - clz) * p) * a
3771                       = radix * a
3772                    so mod (radix * value r sy + _l) radix
3773                       = mod _l radix
3774                    so mod (radix * value r sy + _l) radix
3775                       = mod (radix * a) radix = 0
3776                    so mod _l radix = 0
3777                    so 0 <= _l < radix
3778                  };
3779             assert { value nx sy = p * value r sy
3780                      by
3781                      radix * value r sy
3782                      = power 2 (Limb.length - clz) * value nx sy
3783                      so p * power 2 (Limb.length - clz)
3784                         = radix
3785                      so p * radix * value r sy
3786                         = p * power 2 (Limb.length - clz) * value nx sy
3787                         = radix * value nx sy
3788                      so p * value r sy = value nx sy
3789                    }
3790           end;
3791           assert { value x sx
3792                    = value q (sx - sy + adjust) * value y sy
3793                    + value r sy
3794                   by
3795                    value (nx at Div_s) (sx + adjust)
3796                    = (value q (sx - sy + adjust)
3797                        + power radix (sx - sy + adjust) * _qh)
3798                      * value ny sy
3799                      + value nx sy
3800                   so value (nx at Div_s) (sx + adjust)
3801                       = p * value x sx
3802                   so power radix (sx - sy + 1) * power radix (sy - 1)
3803                      = power radix sx
3804                   so value ny sy = p * value y sy
3805                   so (_qh = 0
3806                      by
3807                      value (nx at Div_s) (sx + adjust)
3808                      = (value q (sx - sy + adjust)
3809                         + power radix (sx - sy + adjust) * _qh)
3810                       * value ny sy
3811                       + value nx sy
3812                      so value nx sy >= 0
3813                      so value q (sx - sy + adjust) >= 0
3814                      so _qh >= 0
3815                      so (value q (sx - sy + adjust)
3816                              + power radix (sx - sy + adjust) * _qh) >= 0
3817                      so (value q (sx - sy + adjust)
3818                          + power radix (sx - sy + adjust) * _qh)
3819                         * value ny sy
3820                         + value nx sy
3821                         >= (value q (sx - sy + adjust)
3822                              + power radix (sx - sy + adjust) * _qh)
3823                            * value ny sy
3824                         >= power radix (sx - sy + adjust) * _qh * value ny sy
3825                      so _qh <> 1)
3826                   so value nx sy = p * value r sy
3827                   so p * value x sx
3828                      = value q (sx - sy + adjust) * p * value y sy
3829                        + p * value r sy
3830                      = p * (value q (sx - sy + adjust)
3831                        * value y sy
3832                             + value r sy)
3833                  };
3834           label Ret_s in
3835           begin
3836             ensures { value q (sx - sy + 1)
3837                       = value (q at Ret_s) (sx - sy + adjust) }
3838             if (adjust = 0)
3839             then begin
3840               value_sub_shift_no_change (pelts x) x.offset
3841                                         (p2i sx) (p2i sx) 0;
3842               set_ofs q (sx - sy) 0;
3843               value_sub_tail (pelts q) q.offset (q.offset + p2i sx - p2i sy);
3844               assert { value q (sx - sy + 1) = value (q at Ret_s) (sx - sy)
3845                        by value q (sx - sy + 1)
3846                           = value (q at Ret_s) (sx - sy) + power radix (sx - sy) * 0
3847                           = value (q at Ret_s) (sx - sy) }
3848             end
3849           end;
3850           ()
3851         end
3852         end
3853
3854   let wmpn_tdiv_qr (q r:t) (qxn:int32) (x:t) (sx:int32) (y:t) (sy:int32) : unit
3855     requires { 1 <= sy <= sx <= (Int32.max_int32 - 1) }
3856     requires { valid x sx }
3857     requires { valid y sy }
3858     requires { valid q (sx - sy + 1) }
3859     requires { valid r sy }
3860     requires { writable q /\ writable r }
3861     requires { qxn = 0 }
3862     requires { (pelts y)[y.offset + sy - 1] > 0 }
3863     ensures { value x sx
3864               = value q (sx - sy + 1) * value y sy
3865                 + value r sy }
3866     ensures { value r sy < value y sy }
3867   =
3868     let nx = malloc (UInt32.(+) (UInt32.of_int32 sx) 1) in
3869     c_assert (is_not_null nx);
3870     let ny = malloc (UInt32.of_int32 sy) in
3871     c_assert (is_not_null ny);
3872     div_qr q r x y nx ny sx sy;
3873     free nx;
3874     free ny
3875
3876   let div_qr_in_place (q x y nx ny:t) (sx sy:int32) : unit
3877     requires { 1 <= sy <= sx <= (Int32.max_int32 - 1) }
3878     requires { valid x sx }
3879     requires { valid y sy }
3880     requires { valid q (sx - sy + 1) }
3881     requires { valid nx (sx + 1) }
3882     requires { valid ny sy }
3883     requires { writable q /\ writable x }
3884     requires { writable nx /\ writable ny }
3885     requires { (pelts y)[y.offset + sy - 1] > 0 }
3886     ensures { value (old x) sx
3887               = value q (sx - sy + 1) * value y sy
3888                 + value x sy }
3889     ensures { value x sy < value y sy }
3890   =
3891     label Start in
3892     value_sub_tail (pelts y) y.offset (y.offset + p2i sy - 1);
3893     value_sub_lower_bound (pelts y) y.offset (y.offset + p2i sy - 1);
3894     assert { value y sy >= power radix (sy - 1) };
3895     let ghost ox = pure { x } in
3896     if (sy = 1)
3897     then
3898       let lr = wmpn_divrem_1 q x sx (C.get y) in
3899       C.set x lr
3900     else
3901     if (sy = 2)
3902     then
3903       let clz = clz_ext (C.get_ofs y (sy - 1)) in
3904       let ghost p = power 2 (p2i clz) in
3905       if clz = 0
3906       then begin
3907         wmpn_copyi nx x sx;
3908         value_sub_shift_no_change (pelts x) x.offset (p2i sx) (p2i sx) 0;
3909         C.set_ofs nx sx 0;
3910         value_sub_frame_shift (pelts x) (pelts nx) x.offset nx.offset (p2i sx);
3911         label Div2_ns in
3912         let ghost _qh = wmpn_divrem_2 q nx y (sx + 1) in
3913         wmpn_copyi x nx sy;
3914         assert { value ox sx
3915                  = value q (sx - sy + 1) * value y sy
3916                  + value x sy
3917                by value x sy = value nx sy
3918                so value (nx at Div2_ns) (sx + 1) < power radix sx
3919                so value (nx at Div2_ns) (sx + 1)
3920                   = value (nx at Div2_ns) sx
3921                so (_qh = 0
3922                    by
3923                    power radix sx
3924                    > value (nx at Div2_ns) (sx + 1)
3925                    = (value q (sx - 1) + power radix (sx - 1) * _qh)
3926                     * value y 2
3927                     + value nx 2
3928                    so value nx 2 >= 0
3929                    so value y 2 >= radix
3930                    so value q (sx - 1) >= 0
3931                    so _qh >= 0
3932                    so (value q (sx - 1)
3933                            + power radix (sx - 1) * _qh) >= 0
3934                    so (value q (sx - 1) + power radix (sx - 1) * _qh)
3935                       * value y 2
3936                       + value nx 2
3937                       >= (value q (sx - 1)
3938                            + power radix (sx - 1) * _qh)
3939                          * value y 2
3940                       >= (value q (sx - 1)
3941                            + power radix (sx - 1) * _qh)
3942                          * radix
3943                       >= power radix (sx - 1) * _qh * radix
3944                       = power radix sx * _qh
3945                    so power radix sx > power radix sx * _qh
3946                   )
3947                so value ox sx = value (nx at Div2_ns) sx
3948                };
3949         ()
3950       end
3951       else begin
3952         let ghost _c = wmpn_lshift ny y sy (Limb.of_int32 clz) in
3953         begin
3954           ensures { normalized ny sy }
3955           ensures { value ny sy = power 2 clz * value y sy }
3956           let ghost dh = (pelts y)[y.offset + p2i sy - 1] in
3957           assert { value y sy
3958                    = value y (sy - 1) + power radix (sy - 1) * dh };
3959           value_sub_upper_bound (pelts y) y.offset (y.offset + p2i sy - 1);
3960           value_sub_tail (pelts ny) ny.offset (ny.offset + p2i sy - 1);
3961           value_sub_upper_bound (pelts ny) ny.offset (ny.offset + p2i sy - 1);
3962           let ghost ndh = (pelts ny)[ny.offset + p2i sy - 1] in
3963           assert { normalized ny sy
3964                    /\ value ny sy = power 2 clz * value y sy
3965                    by
3966                    value y sy < (dh + 1) * power radix (sy - 1)
3967                    so value ny sy + (power radix sy) * _c
3968                       = power 2 clz * value y sy
3969                       = power 2 clz
3970                         * (value y (sy - 1)
3971                            + dh * power radix (sy - 1))
3972                    so power 2 clz * dh <= radix - power 2 clz
3973                    so value ny sy + (power radix sy) * _c
3974                       = power 2 clz * value y (sy - 1)
3975                         + power 2 clz * dh * power radix (sy - 1)
3976                       < power 2 clz * power radix (sy - 1)
3977                         + power 2 clz * dh * power radix (sy - 1)
3978                       <= power 2 clz * power radix (sy - 1)
3979                         + (radix - power 2 clz) * power radix (sy - 1)
3980                        = radix * power radix (sy - 1)
3981                        = power radix sy
3982                    so _c = 0
3983                    so value ny sy
3984                       = power 2 clz * value y sy
3985                    so value y sy >= dh * power radix (sy - 1)
3986                    so value ny sy
3987                       >= power 2  clz * dh * power radix (sy - 1)
3988                    so value ny sy =
3989                       value ny (sy - 1) + power radix (sy - 1) * ndh
3990                       < power radix (sy - 1) + power radix (sy - 1) * ndh
3991                       = power radix (sy - 1) * (ndh + 1)
3992                    so power radix (sy - 1) * (ndh + 1)
3993                       > power radix (sy - 1) * (power 2 clz * dh)
3994                    so ndh + 1 > power 2 clz * dh
3995                    so ndh >= power 2 clz * dh
3996                    so 2 * power 2 clz * dh >= radix
3997                    so 2 * ndh >= radix
3998                    so ndh >= div radix 2
3999                  };
4000         end;
4001         let h = wmpn_lshift nx x sx (Limb.of_int32 clz) in
4002         C.set_ofs nx sx h;
4003         begin
4004           ensures { value nx (sx + 1)
4005                     = p * value ox sx  }
4006           value_sub_tail (pelts nx) nx.offset (nx.offset + p2i sx);
4007           assert { value nx (sx + 1)
4008                     = p * value ox sx
4009                    by
4010                    value nx sx + power radix sx * h
4011                     = p * value ox sx
4012                    so value nx (sx + 1)
4013                       = value nx sx + power radix sx * h
4014                  }
4015         end;
4016         label Div2_s in
4017         (* TODO don't add 1 when not needed, cf "adjust" in GMP algo *)
4018         let ghost _qh = wmpn_divrem_2 q nx ny (sx + 1) in
4019         let ghost _l = wmpn_rshift x nx sy (Limb.of_int32 clz) in
4020         begin ensures { value nx 2 = p * value x 2 }
4021           assert { _l = 0
4022                  by
4023                    (mod (value nx sy) p = 0
4024                     by
4025                       value (nx at Div2_s) (sx + 1)
4026                    = (value q (sx - 1) + power radix (sx - 1) * _qh)
4027                       * value ny sy
4028                       + value nx sy
4029                    so value (nx at Div2_s) (sx + 1)
4030                     = p * value ox sx
4031                    so value ny sy = p * value y sy
4032                    so value nx sy
4033                       = value (nx at Div2_s) (sx + 1)
4034                         - (value q (sx - 1)
4035                           + power radix (sx - 1) * _qh)
4036                           * value ny sy
4037                       = p * value ox sx
4038                         - p * (value q (sx - 1)
4039                                + power radix (sx - 1) * _qh)
4040                             * value y sy
4041                       = p * (value ox sx
4042                              - (value q (sx - 1)
4043                                + power radix (sx - 1) * _qh)
4044                                * value y sy)
4045                    so let n = (value ox sx
4046                                 - (value q (sx - 1)
4047                                    + power radix (sx - 1) * _qh)
4048                                    * value y sy)
4049                       in
4050                       value nx sy = p * n
4051                       so value nx sy >= 0
4052                       so p > 0
4053                       so n >= 0
4054                       so mod (value nx sy) p
4055                          = mod (p * n) p
4056                          = mod ((p*n)+0) p
4057                          = mod 0 p
4058                          = 0
4059                  )
4060                  so _l + radix * value x sy
4061                     = power 2 (Limb.length - clz) * (value nx sy)
4062                  so let a = div (value nx sy) p in
4063                     value nx sy = p * a
4064                  so power 2 (Limb.length - clz) * p = radix
4065                  so power 2 (Limb.length - clz) * (value nx sy)
4066                     = power 2 (Limb.length - clz) * (p * a)
4067                     = (power 2 (Limb.length - clz) * p) * a
4068                     = radix * a
4069                  so mod (radix * value x sy + _l) radix
4070                     = mod _l radix
4071                  so mod (radix * value x sy + _l) radix
4072                     = mod (radix * a) radix = 0
4073                  so mod _l radix = 0
4074                  so 0 <= _l < radix
4075                };
4076           assert { value nx 2 = p * value x 2
4077                    by
4078                    radix * value x 2
4079                    = power 2 (Limb.length - clz) * value nx 2
4080                    so p * power 2 (Limb.length - clz)
4081                       = radix
4082                    so p * radix * value x 2
4083                       = p * power 2 (Limb.length - clz) * value nx 2
4084                       = radix * value nx 2
4085                    so p * value x 2 = value nx 2
4086                  }
4087         end;
4088         assert { value ox sx
4089                  = value q (sx - sy + 1) * value y sy
4090                  + value x sy
4091                 by
4092                  value (nx at Div2_s) (sx + 1)
4093                  = (value q (sx - 1) + power radix (sx - 1) * _qh)
4094                    * value ny 2
4095                    + value nx 2
4096                 so value (nx at Div2_s) (sx + 1)
4097                     = p * value ox sx
4098                 so value ny 2 = p * value y 2
4099                 so (_qh = 0
4100                     by
4101                     value ox sx < power radix sx
4102                     so value y 2 >= radix
4103                     so value ny 2 >= p * radix
4104                     so value q (sx - 1) >= 0
4105                     so value nx 2 >= 0
4106                     so (value q (sx - 1) + power radix (sx - 1) * _qh)
4107                         >= 0
4108                     so (value q (sx - 1) + power radix (sx - 1) * _qh)
4109                          * value ny 2
4110                        + value nx 2
4111                        >= (value q (sx - 1)
4112                            + power radix (sx - 1) * _qh)
4113                           * value ny 2
4114                        >= (value q (sx - 1)
4115                            + power radix (sx - 1) * _qh)
4116                           * (p * radix)
4117                        >= power radix (sx - 1) * _qh * p * radix
4118                        = power radix sx * p * _qh
4119                     so power radix sx * p
4120                        > value (nx at Div2_s) (sx + 1)
4121                        >= power radix sx * p * _qh
4122                     )
4123                 so value nx 2 = p * value x 2
4124                 so p * value ox sx
4125                    = value q (sx - 1) * p * value y 2
4126                      + p * value x 2
4127                    = p * (value q (sx - 1) * value y 2
4128                           + value x 2)
4129                };
4130         ()
4131       end
4132     else
4133      (* let qn = ref (Int32.(-) (Int32.(+) sx 1) sy) in
4134       if (Int32.(>=)  (Int32.(+) !qn !qn) sx)
4135       then*) begin
4136         let adjust =
4137           if (get_ofs x (sx - 1)) >= (get_ofs y (sy - 1))
4138           then 1
4139           else 0
4140         in
4141         let clz = clz_ext (C.get_ofs y (sy - 1)) in
4142         let ghost p = power 2 (p2i clz) in
4143         if clz = 0
4144         then begin
4145           wmpn_copyi nx x sx;
4146           value_sub_shift_no_change (pelts x) x.offset
4147                                     (p2i sx) (p2i sx) 0;
4148           C.set_ofs nx sx 0;
4149           value_sub_frame_shift (pelts x) (pelts nx) x.offset nx.offset (p2i sx);
4150           assert { value y sy * (power radix (sx - sy + adjust))
4151                    > value nx (sx + adjust)
4152                    by
4153                    let dh = (pelts y)[y.offset + sy - 1] in
4154                    value y sy >= dh * power radix (sy - 1)
4155                    so value nx (sx + adjust) = value nx sx = value ox sx
4156                    so [@case_split]
4157                       ((adjust = 1
4158                            so value ox sx < power radix sx
4159                            so value y sy * power radix (sx - sy + adjust)
4160                               >= dh * power radix (sy - 1)
4161                                     * power radix (sx - sy + adjust)
4162                               = dh * power radix ((sy - 1) + (sx - sy + adjust))
4163                               = dh * power radix sx
4164                            so dh >= div radix 2 > 1
4165                            so dh * power radix sx > power radix sx )
4166                       \/
4167                       (adjust = 0
4168                            so let ah = (pelts x)[x.offset + sx - 1] in
4169                               value ox sx < (ah + 1) * power radix (sx - 1)
4170                               so ah + 1 <= dh
4171                               so value ox sx < dh * power radix (sx - 1)
4172                               so value y sy * power radix (sx - sy + adjust)
4173                                  = value y sy * power radix (sx - sy)
4174                                  >= dh * power radix (sy - 1)
4175                                        * power radix (sx - sy)
4176                                  = dh * power radix (sy - 1 + sx - sy)
4177                                  = dh * power radix (sx - 1))) };
4178           label Div_ns in
4179           let ghost _qh = div_sb_qr q nx (sx + adjust) y sy in
4180           wmpn_copyi x nx sy;
4181           assert { value ox sx
4182                    = value q (sx - sy + adjust) * value y sy
4183                    + value x sy
4184                  by value x sy = value nx sy
4185                  so value (nx at Div_ns) (sx + adjust) = value ox sx < power radix sx
4186                  so value (nx at Div_ns) (sx + adjust)
4187                     = value (nx at Div_ns) sx
4188                  so (_qh = 0
4189                      by
4190                      value (nx at Div_ns) (sx + adjust)
4191                      = (value q (sx - sy + adjust)
4192                         + power radix (sx - sy + adjust) * _qh)
4193                       * value y sy
4194                       + value nx sy
4195                      so value nx sy >= 0
4196                      so value q (sx - sy + adjust) >= 0
4197                      so _qh >= 0
4198                      so (value q (sx - sy + adjust)
4199                              + power radix (sx - sy + adjust) * _qh) >= 0
4200                      so (value q (sx - sy + adjust)
4201                          + power radix (sx - sy + adjust) * _qh)
4202                         * value y sy
4203                         + value nx sy
4204                         >= (value q (sx - sy + adjust)
4205                              + power radix (sx - sy + adjust) * _qh)
4206                            * value y sy
4207                         >= power radix (sx - sy + adjust) * _qh * value y sy
4208                      so _qh <> 1)
4209                  so value ox sx = value (nx at Div_ns) sx
4210                  };
4211            label Ret_ns in
4212            begin
4213              ensures { value q (sx - sy + 1)
4214                        = value (q at Ret_ns) (sx - sy + adjust) }
4215              if (adjust = 0)
4216              then begin
4217                value_sub_shift_no_change (pelts x) x.offset
4218                                          (p2i sx) (p2i sx) 0;
4219                set_ofs q (sx - sy) 0;
4220                value_sub_tail (pelts q) q.offset (q.offset + p2i sx - p2i sy);
4221                ()
4222              end
4223            end
4224         end
4225         else begin
4226           let ghost _c = wmpn_lshift ny y sy (Limb.of_int32 clz) in
4227           begin
4228             ensures { normalized ny sy }
4229             ensures { value ny sy
4230                       = power 2 clz * value y sy }
4231             let ghost dh = (pelts y)[y.offset + p2i sy - 1] in
4232             assert { value y sy
4233                      = value y (sy - 1) + power radix (sy - 1) * dh };
4234             value_sub_upper_bound (pelts y) y.offset (y.offset + p2i sy - 1);
4235             value_sub_tail (pelts ny) ny.offset (ny.offset + p2i sy - 1);
4236             value_sub_upper_bound (pelts ny) ny.offset (ny.offset + p2i sy - 1);
4237             let ghost ndh = (pelts ny)[ny.offset + p2i sy - 1] in
4238             assert { normalized ny sy
4239                      /\ value ny sy
4240                         = power 2 clz * value y sy
4241                      by
4242                      value y sy < (dh + 1) * power radix (sy - 1)
4243                      so value ny sy + (power radix sy) * _c
4244                         = power 2 clz * value y sy
4245                         = power 2 clz
4246                           * (value y (sy - 1)
4247                              + dh * power radix (sy - 1))
4248                      so power 2 clz * dh <= radix - power 2 clz
4249                      so (_c = 0
4250                         by
4251                           value ny sy + (power radix sy) * _c
4252                           = power 2 clz * value y (sy - 1)
4253                             + power 2 clz * dh * power radix (sy - 1)
4254                           < power 2 clz * power radix (sy - 1)
4255                             + power 2 clz * dh * power radix (sy - 1)
4256                           <= power 2 clz * power radix (sy - 1)
4257                             + (radix - power 2 clz) * power radix (sy - 1)
4258                            = radix * power radix (sy - 1)
4259                            = power radix sy
4260                         so value ny sy >= 0
4261                         so power radix sy * _c < power radix sy
4262                         so power radix sy > 0
4263                         so _c >= 0
4264                         )
4265                      so value ny sy
4266                         = power 2 clz * value y sy
4267                      so value y sy >= dh * power radix (sy - 1)
4268                      so value ny sy
4269                         >= power 2  clz * dh * power radix (sy - 1)
4270                      so value ny sy =
4271                         value ny (sy - 1) + power radix (sy - 1) * ndh
4272                         < power radix (sy - 1) + power radix (sy - 1) * ndh
4273                         = power radix (sy - 1) * (ndh + 1)
4274                      so power radix (sy - 1) * (ndh + 1)
4275                         > power radix (sy - 1) * (power 2 clz * dh)
4276                      so ndh + 1 > power 2 clz * dh
4277                      so ndh >= power 2 clz * dh
4278                      so 2 * power 2 clz * dh >= radix
4279                      so 2 * ndh >= radix
4280                      so ndh >= div radix 2
4281                    };
4282           end;
4283           let h = wmpn_lshift nx x sx (Limb.of_int32 clz) in
4284           label Shifted in
4285           C.set_ofs nx sx h;
4286           begin
4287             ensures { value nx (sx + adjust)
4288                       = p * value ox sx  }
4289             if (adjust = 1)
4290             then begin
4291               value_sub_tail (pelts nx) nx.offset (nx.offset + p2i sx);
4292               assert { value nx (sx + 1)
4293                         = p * value ox sx
4294                        by
4295                        value nx sx + power radix sx * h
4296                         = p * value ox sx
4297                        so value nx (sx + 1)
4298                           = value nx sx + power radix sx * h
4299                      } end
4300             else begin
4301               assert { adjust = 0 };
4302               assert { h = 0
4303                        by
4304                        let dh = (pelts y)[y.offset + sy - 1] in
4305                        let ah = (pelts x)[x.offset + sx - 1] in
4306                        p * dh < radix
4307                        so 0 <= ah <= dh
4308                        so p * ah < radix
4309                        so (p * ah <= radix - p
4310                            by
4311                            let q = power 2 (Limb.length - clz) in
4312                            radix = p * q
4313                            so p * ah < p * q
4314                            so ah < q
4315                            so ah <= q - 1
4316                            so p * ah <= p * (q - 1) = radix - p
4317                            )
4318                        so p * (ah + 1) <= radix
4319                        so let s = power radix (sx - 1) in
4320                           value ox sx < (ah + 1) * s
4321                        so p * value ox sx < p * (ah + 1) * s
4322                        so (p * (ah + 1) * s
4323                            <= radix * s
4324                            by
4325                            [@case_split]
4326                            (p * (ah + 1) = radix
4327                            \/ (p * (ah + 1) < radix
4328                                so s > 0
4329                                so p * (ah + 1) * s
4330                                  < radix * s)))
4331                        so radix * power radix (sx - 1) = power radix sx
4332                        so value (nx at Shifted) sx + power radix sx * h
4333                           < power radix sx
4334                        so power radix sx * h < power radix sx * 1
4335                        so (h < 1 by power radix sx > 0)
4336                       }
4337             end
4338           end;
4339           label Div_s in
4340           assert { value ny sy * (power radix (sx - sy + adjust))
4341                    > value nx (sx + adjust)
4342                    by
4343                    let dh = (pelts y)[y.offset + sy - 1] in
4344                    value ny sy >= p * dh * power radix (sy - 1)
4345                    so value nx (sx + adjust) = p * value ox sx
4346                    so p > 0
4347                    so [@case_split]
4348                       ((adjust = 1
4349                            so value ox sx < power radix sx
4350                            so p * value ox sx < p * power radix sx
4351                            so value ny sy * power radix (sx - sy + adjust)
4352                               >= p * dh * power radix (sy - 1)
4353                                     * power radix (sx - sy + adjust)
4354                               = p * dh * power radix ((sy - 1) + (sx - sy + adjust))
4355                               = p * dh * power radix sx
4356                            so dh >= 1
4357                            so p * dh * power radix sx >= p * power radix sx )
4358                       \/
4359                       (adjust = 0
4360                            so let ah = (pelts x)[x.offset + sx - 1] in
4361                               value ox sx < (ah + 1) * power radix (sx - 1)
4362                               so ah + 1 <= dh
4363                               so value ox sx < dh * power radix (sx - 1)
4364                               so p * value ox sx < p * dh * power radix (sx - 1)
4365                               so value ny sy * power radix (sx - sy + adjust)
4366                                  = value ny sy * power radix (sx - sy)
4367                                  >= p * dh * power radix (sy - 1)
4368                                        * power radix (sx - sy)
4369                                  = p * dh * power radix (sy - 1 + sx - sy)
4370                                  = p * dh * power radix (sx - 1))) };
4371           let ghost _qh = div_sb_qr q nx (sx + adjust) ny sy in
4372           let ghost _l = wmpn_rshift x nx sy (Limb.of_int32 clz) in
4373           begin ensures { value nx sy = p * value x sy }
4374             assert { _l = 0
4375                    by
4376                      (mod (value nx sy) p = 0
4377                       by
4378                         value (nx at Div_s) (sx + adjust)
4379                      = (value q (sx - sy + adjust)
4380                          + power radix (sx - sy + adjust) * _qh)
4381                         * value ny sy
4382                         + value nx sy
4383                      so value (nx at Div_s) (sx + adjust)
4384                       = p * value ox sx
4385                      so value ny sy = p * value y sy
4386                      so value nx sy
4387                         = value (nx at Div_s) (sx + adjust)
4388                           - (value q (sx - sy + adjust)
4389                             + power radix (sx - sy + adjust) * _qh)
4390                             * value ny sy
4391                         = p * value ox sx
4392                           - p * (value q (sx - sy + adjust)
4393                                  + power radix (sx - sy + adjust) * _qh)
4394                               * value y sy
4395                         = p * (value ox sx
4396                                - (value q (sx - sy + adjust)
4397                                  + power radix (sx - sy + adjust) * _qh)
4398                                  * value y sy)
4399                      so let n = (value ox sx
4400                                   - (value q (sx - sy + adjust)
4401                                      + power radix (sx - sy + adjust) * _qh)
4402                                      * value y sy)
4403                         in
4404                         value nx sy = p * n
4405                         so value nx sy >= 0
4406                         so p > 0
4407                         so n >= 0
4408                         so mod (value nx sy) p
4409                            = mod (p * n) p
4410                            = mod ((p*n)+0) p
4411                            = mod 0 p
4412                            = 0
4413                    )
4414                    so _l + radix * value x sy
4415                       = power 2 (Limb.length - clz) * (value nx sy)
4416                    so let a = div (value nx sy) p in
4417                       value nx sy = p * a
4418                    so power 2 (Limb.length - clz) * p = radix
4419                    so power 2 (Limb.length - clz) * (value nx sy)
4420                       = power 2 (Limb.length - clz) * (p * a)
4421                       = (power 2 (Limb.length - clz) * p) * a
4422                       = radix * a
4423                    so mod (radix * value x sy + _l) radix
4424                       = mod _l radix
4425                    so mod (radix * value x sy + _l) radix
4426                       = mod (radix * a) radix = 0
4427                    so mod _l radix = 0
4428                    so 0 <= _l < radix
4429                  };
4430             assert { value nx sy = p * value x sy
4431                      by
4432                      radix * value x sy
4433                      = power 2 (Limb.length - clz) * value nx sy
4434                      so p * power 2 (Limb.length - clz)
4435                         = radix
4436                      so p * radix * value x sy
4437                         = p * power 2 (Limb.length - clz) * value nx sy
4438                         = radix * value nx sy
4439                      so p * value x sy = value nx sy
4440                    }
4441           end;
4442           assert { value ox sx
4443                    = value q (sx - sy + adjust) * value y sy
4444                    + value x sy
4445                   by
4446                    value (nx at Div_s) (sx + adjust)
4447                    = (value q (sx - sy + adjust)
4448                        + power radix (sx - sy + adjust) * _qh)
4449                      * value ny sy
4450                      + value nx sy
4451                   so value (nx at Div_s) (sx + adjust)
4452                       = p * value ox sx
4453                   so power radix (sx - sy + 1) * power radix (sy - 1)
4454                      = power radix sx
4455                   so value ny sy = p * value y sy
4456                   so (_qh = 0
4457                      by
4458                      value (nx at Div_s) (sx + adjust)
4459                      = (value q (sx - sy + adjust)
4460                         + power radix (sx - sy + adjust) * _qh)
4461                       * value ny sy
4462                       + value nx sy
4463                      so value nx sy >= 0
4464                      so value q (sx - sy + adjust) >= 0
4465                      so _qh >= 0
4466                      so (value q (sx - sy + adjust)
4467                              + power radix (sx - sy + adjust) * _qh) >= 0
4468                      so (value q (sx - sy + adjust)
4469                          + power radix (sx - sy + adjust) * _qh)
4470                         * value ny sy
4471                         + value nx sy
4472                         >= (value q (sx - sy + adjust)
4473                              + power radix (sx - sy + adjust) * _qh)
4474                            * value ny sy
4475                         >= power radix (sx - sy + adjust) * _qh * value ny sy
4476                      so _qh <> 1)
4477                   so value nx sy = p * value x sy
4478                   so p * value ox sx
4479                      = value q (sx - sy + adjust) * p * value y sy
4480                        + p * value x sy
4481                      = p * (value q (sx - sy + adjust)
4482                        * value y sy
4483                             + value x sy)
4484                  };
4485           label Ret_s in
4486           begin
4487             ensures { value q (sx - sy + 1)
4488                       = value (q at Ret_s) (sx - sy + adjust) }
4489             if (adjust = 0)
4490             then begin
4491               value_sub_shift_no_change (pelts x) x.offset
4492                                         (p2i sx) (p2i sx) 0;
4493               set_ofs q (sx - sy) 0;
4494               value_sub_tail (pelts q) q.offset (q.offset + p2i sx - p2i sy);
4495               assert { value q (sx - sy + 1) = value (q at Ret_s) (sx - sy)
4496                        by value q (sx - sy + 1)
4497                           = value (q at Ret_s) (sx - sy) + power radix (sx - sy) * 0
4498                           = value (q at Ret_s) (sx - sy) }
4499             end
4500           end;
4501           ()
4502         end
4503         end
4504
4505   let wmpn_tdiv_qr_in_place (q:t) (qxn:int32) (x:t) (sx:int32) (y:t) (sy:int32)
4506       : unit
4507     requires { 1 <= sy <= sx <= (Int32.max_int32 - 1) }
4508     requires { valid x sx }
4509     requires { valid y sy }
4510     requires { valid q (sx - sy + 1) }
4511     requires { writable x /\ writable q }
4512     requires { qxn = 0 }
4513     requires { (pelts y)[y.offset + sy - 1] > 0 }
4514     ensures { value (old x) sx
4515               = value q (sx - sy + 1) * value y sy
4516                 + value x sy }
4517     ensures { value x sy < value y sy }
4518   =
4519     let nx = malloc (UInt32.(+) (UInt32.of_int32 sx) 1) in
4520     c_assert (is_not_null nx);
4521     let ny = malloc (UInt32.of_int32 sy) in
4522     c_assert (is_not_null ny);
4523     div_qr_in_place q x y nx ny sx sy;
4524     free nx;
4525     free ny
4526
4527
4528 end