examples/multiprecision/toom.mlw

   1 module Toom
   2
   3   use array.Array
   4   use map.Map
   5   use mach.c.C
   6   use ref.Ref
   7   use mach.int.Int32
   8   use import mach.int.UInt64GMP as Limb
   9   use int.Int
  10   use int.Power
  11   use valuation.Valuation
  12   use int.ComputerDivision
  13   use types.Types
  14   use types.Int32Eq
  15   use types.UInt64Eq
  16   use lemmas.Lemmas
  17   use compare.Compare
  18   use util.Util
  19   use util.UtilOld
  20   use add_1.Add_1
  21   use add.AddOld
  22   use sub_1.Sub_1
  23   use sub.SubOld
  24   use mul.Mul
  25   use mul.Mul_basecase
  26   use logical.Logical
  27
  28 let constant toom22_threshold : int32 = 29
  29
  30 let lemma no_borrow (x y r b m:int)
  31   requires { 0 <= y <= x }
  32   requires { 0 <= r < m }
  33   requires { r - m * b = x - y }
  34   requires { 0 <= b }
  35   ensures  { b = 0 }
  36 = assert { b <= 0 by m * b < m * 1 }
  37
  38 let lemma no_borrow_ptr (x y r: ptr limb) (nx ny:int) (b:limb)
  39   requires { 0 < ny <= nx }
  40   requires { value y ny <= value x nx }
  41   requires { value r nx - (power radix nx) * b = value x nx - value y ny }
  42   requires { 0 <= b }
  43   ensures  { b = 0 }
  44 = no_borrow (value x nx) (value y ny) (value r nx) (l2i b) (power radix nx)
  45
  46 let rec wmpn_toom22_mul (r x:ptr limb) (sx:int32) (y:ptr limb) (sy:int32)
  47                         (scratch: ptr limb) (ghost k: int) : unit
  48   requires { valid x sx }
  49   requires { valid y sy }
  50   requires { valid r (sx + sy) }
  51   requires { toom22_threshold < sy }
  52   requires { 0 < k <= 64 }
  53   requires { sx <= toom22_threshold * power 2 k }
  54   requires { valid scratch (2 * (sx + k)) } (*TODO faire en fonction de sy *)
  55   requires { writable r /\ writable scratch }
  56   requires { 8 * sx < max_int32 }
  57   requires { 2 < sy <= sx < sy + sy - 1 }
  58   requires { 4 * sx < 5 * sy }
  59   ensures  { min r = old min r }
  60   ensures  { max r = old max r }
  61   ensures  { plength r = old plength r }
  62   ensures  { min scratch = old min scratch }
  63   ensures  { max scratch = old max scratch }
  64   ensures  { plength scratch = old plength scratch }
  65   ensures  { value r (sx + sy) = value x sx * value y sy }
  66   ensures  { forall j. min r <= j < offset r \/ offset r + sx + sy <= j < max r
  67                        -> (pelts r)[j] = old (pelts r)[j] }
  68   ensures  { forall j. min scratch <= j < offset scratch
  69                        -> (pelts scratch)[j] = old (pelts scratch)[j] }
  70   variant { k + k }
  71 =
  72   let s = sx / 2 in (* TODO sx >> 1 *)
  73   let n = sx - s in
  74   let t = sy - n in
  75   assert { 0 < s };
  76   assert { n-1 <= s <= n };
  77   assert { 0 < t <= s };
  78   let x0 = x in
  79   let x1 = C.incr x n in
  80   let y0 = y in
  81   let y1 = C.incr y n in
  82   let ghost a0 = value x0 (int32'int n) in
  83   let ghost a1 = value x1 (int32'int s) in
  84   let ghost b0 = value y0 (int32'int n) in
  85   let ghost b1 = value y1 (int32'int t) in
  86   let ghost m = power radix (int32'int n) in
  87   value_concat x n sx;
  88   value_concat y n sy;
  89   assert { value x sx = a0 + m * a1 };
  90   assert { value y sy = b0 + m * b1 };
  91   let r' = decr_split r 0 in
  92   let ro = C.incr_split r (sx + sy) in
  93   let scratch' = decr_split scratch 0 in
  94   assert { min r = offset r /\ max r = offset r + sx + sy };
  95   let s_out = C.incr_split scratch (n + n) in
  96   let vinf = C.incr_split r (n + n) in
  97   label ASM1 in
  98   let xsm1 = r in
  99   let ysm1 = C.incr_split r n in
 100   let vm1_neg = ref false in
 101   begin ensures { (!vm1_neg /\ value xsm1 n = a1 - a0)
 102                   \/ (not !vm1_neg /\ value xsm1 n = a0 - a1) }
 103         ensures { min scratch = old min scratch }
 104         ensures { max scratch = old max scratch }
 105         ensures { plength scratch = old plength scratch }
 106     if (s = n)
 107     then
 108       if begin ensures { result <->  value x0 n < value x1 n }
 109            (wmpn_cmp x0 x1 n) < 0
 110          end
 111       then begin
 112         let ghost b = wmpn_sub_n xsm1 x1 x0 n in
 113         no_borrow_ptr x1 x0 xsm1 (p2i n) (p2i n) b;
 114         vm1_neg := true end
 115       else let ghost b = wmpn_sub_n xsm1 x0 x1 n in
 116            no_borrow_ptr x0 x1 xsm1 (p2i n) (p2i n) b
 117     else
 118       (* n-s=1*)
 119       if ((get_ofs x0 s) = 0) &&
 120          ((wmpn_cmp x0 x1 s) < 0)
 121       then begin
 122         assert { value x0 s < value x1 s };
 123         value_tail x0 s;
 124         assert { value x0 n = value x0 s };
 125         let ghost b = wmpn_sub_n xsm1 x1 x0 s in
 126         no_borrow_ptr x1 x0 xsm1 (p2i s) (p2i s) b;
 127         value_sub_shift_no_change (pelts xsm1) (xsm1.offset) (p2i s) (p2i s) 0;
 128         set_ofs xsm1 s 0;
 129         vm1_neg := true;
 130         value_tail xsm1 s
 131         end
 132       else begin
 133         let b = wmpn_sub_n xsm1 x0 x1 s in
 134         label Borrow in
 135         value_sub_shift_no_change (pelts xsm1) (xsm1.offset) (p2i s) (p2i s) b;
 136         let lx = get_ofs x0 s in
 137         set_ofs xsm1 s (lx - b);
 138         assert { value xsm1 s = value (xsm1 at Borrow) s };
 139         value_tail x0 s;
 140         value_tail xsm1 s;
 141         end
 142   end;
 143   label BSM1 in
 144   begin ensures { (!vm1_neg = (!vm1_neg at BSM1) /\ value ysm1 n = b0 - b1)
 145                   \/ (!vm1_neg = not (!vm1_neg at BSM1) /\ value ysm1 n = b1 - b0) }
 146         ensures { value xsm1 n = (value xsm1 n at BSM1) }
 147         ensures { min scratch = old min scratch }
 148         ensures { max scratch = old max scratch }
 149         ensures { plength scratch = old plength scratch }
 150     if (t = n)
 151     then
 152       if begin ensures { result <-> value y0 n < value y1 n }
 153            (wmpn_cmp y0 y1 n) < 0
 154          end
 155       then begin
 156         let ghost b = wmpn_sub_n ysm1 y1 y0 n in
 157         no_borrow_ptr y1 y0 ysm1 (p2i n) (p2i n) b;
 158         vm1_neg := not !vm1_neg end
 159       else
 160         let ghost b = wmpn_sub_n ysm1 y0 y1 n in
 161         no_borrow_ptr y0 y1 ysm1 (p2i n) (p2i n) b;
 162     else
 163       let y0t = C.incr y0 t in
 164       let c0 = ((wmpn_zero_p y0t (n - t)) = 1) in
 165       let c1 = ((wmpn_cmp y0 y1 t) < 0) in
 166       if c0 && c1
 167       then begin
 168         assert { value y0 t < value y1 t };
 169         value_concat y0 t n;
 170         assert { value y0 n = value y0 t };
 171         let ghost b = wmpn_sub_n ysm1 y1 y0 t in
 172         no_borrow_ptr y1 y0 ysm1 (p2i t) (p2i t) b;
 173         assert { forall j. (j < offset r \/ offset r + sx + sy <= j)
 174                              -> (pelts r)[j] = (pelts r)[j] at BSM1 };
 175         label Zero in
 176         let ghost ysm1z = { ysm1 } in
 177         let ysm1t = C.incr ysm1 t in
 178         wmpn_zero ysm1t (n - t);
 179         value_sub_frame_shift (pelts ysm1) (pelts ysm1z)
 180                               (offset ysm1) (offset ysm1z) (p2i t);
 181         assert { value ysm1 t = value ysm1 t at Zero };
 182         assert { value xsm1 n = value xsm1 n at Zero };
 183         value_concat ysm1 t n;
 184         assert { value ysm1 n = value ysm1 t };
 185         vm1_neg := not !vm1_neg end
 186       else begin
 187         value_concat y0 t n;
 188         assert { value y0 n = value y0 t + power radix t * value y0t (n-t) };
 189         assert { value y1 t <= value y0 n
 190                  by (not c0
 191                      so 1 <= value y0t (n - t)
 192                      so power radix t * 1 <= power radix t * value y0t (n-t)
 193                      so power radix t <= value y0 n
 194                      so value y1 t < power radix t)
 195                   \/ (not c1
 196                       so value y1 t <= value y0 t
 197                       so value y0 t <= value y0 n) };
 198         let ghost b = wmpn_sub ysm1 y0 n y1 t in
 199         no_borrow_ptr y0 y1 ysm1 (p2i n) (p2i t) b;
 200         end
 201   end;
 202   let ghost asm1_abs = value xsm1 (int32'int n) in
 203   let ghost bsm1_abs = value ysm1 (int32'int n) in
 204   label RecM1 in
 205   wmpn_toom22_mul_n_rec scratch xsm1 ysm1 s_out n (k-1);
 206   assert { value scratch (n+n) = asm1_abs * bsm1_abs };
 207   join r ysm1;
 208   assert { min scratch = offset scratch };
 209   assert { max scratch = max scratch at ASM1 };
 210   assert { plength scratch = old plength scratch };
 211   assert { min r = offset r };
 212   assert { max r = max r at ASM1 };
 213   assert { plength r = old plength r };
 214   assert { forall j. min scratch <= j < offset scratch -> (pelts scratch)[j]
 215             = old (pelts scratch)[j] };
 216   let v0 = r in
 217   label Rec in
 218   begin ensures { value scratch (n+n) = asm1_abs * bsm1_abs }
 219         ensures { value v0 (n+n) = a0 * b0 }
 220         ensures { value vinf (s+t) = a1 * b1 }
 221         ensures { min scratch = old min scratch }
 222         ensures { max scratch = old max scratch }
 223         ensures { plength scratch = old plength scratch }
 224         ensures { min s_out = old min s_out }
 225         ensures { max s_out = old max s_out }
 226         ensures { min vinf = old min vinf }
 227         ensures { max vinf = old max vinf }
 228         ensures { plength vinf = old plength vinf }
 229         ensures { min r = old min r }
 230         ensures { max r = old max r }
 231         ensures { plength r = old plength r }
 232     let lemma valid_monotonous (s n:int32)
 233       requires { valid s_out (2*(n+(k-1))) }
 234       requires { 0 <= s <= n }
 235       ensures  { valid s_out (2*(s+(k-1))) }
 236     = assert { 0 <= 2*(s+(k-1)) <= 2*(n+(k-1)) };
 237       assert { forall p: ptr limb, n1 n2. 0 <= n1 <= n2 -> valid p n2 -> valid p n1 } in
 238     valid_monotonous s n;
 239     valid_monotonous t n;
 240     (if s > t
 241      then wmpn_toom22_mul_rec vinf x1 s y1 t s_out (k-1)
 242      else wmpn_toom22_mul_n_rec vinf x1 y1 s_out s (k-1));
 243     wmpn_toom22_mul_n_rec v0 x0 y0 s_out n (k-1);
 244   end;
 245   label Adds in
 246   value_concat v0 n (n + n);
 247   value_concat vinf n (s + t);
 248   let v0n = incr_split v0 n in
 249   let vinfn = incr_split vinf n in
 250   let ghost lv0 = value v0 (int32'int n) in
 251   let ghost hv0 = value v0n (int32'int n) in
 252   let ghost lvinf = value vinf (int32'int n) in
 253   let ghost hvinf = value vinfn (int32'int s + int32'int t- int32'int n) in
 254   assert { lv0 + m * hv0 = a0 * b0 };
 255   assert { lvinf + m * hvinf = a1 * b1 };
 256   let cy = ref (wmpn_add_n_in_place vinf v0n n) in
 257   assert { value vinf n = lvinf + hv0 - m * !cy };
 258   let c = wmpn_add_n v0n vinf v0 n in
 259   let cy2 = c + !cy in
 260   assert { value v0n n = lvinf + hv0 + lv0 - m * cy2
 261            by value v0n n = lv0 + value vinf n - m * c
 262                           = lvinf + hv0 + lv0 - m * !cy - m * c
 263                           = lvinf + hv0 + lv0 - m * cy2 };
 264   label Add3 in
 265   let c' = wmpn_add_in_place vinf n vinfn ((s+t) - n) in
 266   cy := !cy + c';
 267   assert { value vinf n = hvinf + lvinf + hv0 - m * !cy
 268            by m * (!cy at Add3) + m * c' = m * !cy
 269            so value vinf n = value vinf n at Add3 + hvinf - m * c'
 270               = lvinf + hv0 - m * (!cy at Add3) + hvinf - m * c'
 271               = hvinf + lvinf + hv0 - m * !cy };
 272   assert { !cy <= 2 };
 273   label JoinMid in
 274   let ghost v0nj = { v0n } in
 275   let ghost vinfj = { vinf } in
 276   join v0n vinf;
 277   value_sub_frame (pelts v0n) (pelts v0nj) (offset v0n) (offset v0n + p2i n);
 278   assert { value v0n n = value v0nj n };
 279   value_sub_frame (pelts v0n) (pelts vinfj) (offset vinf) (offset vinf + p2i n);
 280   assert { value_sub (pelts v0n) (offset v0n + n) (offset v0n + n + n)
 281            = value vinfj n };
 282   value_concat v0n n (n + n);
 283   assert { value v0n (n+n) = a1 * b1 + a0 * b0 + hv0 + m * lvinf
 284              - m * cy2 - m * m * !cy
 285            by value v0n (n+n)
 286            = value v0n n at JoinMid + m * value vinf n at JoinMid
 287            = lvinf + hv0 + lv0 - m * cy2
 288              + m * (hvinf + lvinf + hv0 - m * !cy)
 289            = lvinf + hv0 + lv0 - m * cy2 + m * hvinf + m * lvinf
 290              + m * hv0 - m * m * !cy
 291            = a1 * b1 + a0 * b0 + hv0 + m * lvinf
 292              - m * cy2 - m * m * !cy };
 293   label AddSub in
 294   begin ensures { !cy <= 3 (*2?*)
 295                    /\ value v0n (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 296                                       + hv0 + m * lvinf - m * cy2 - m * m * !cy
 297                   \/ !cy = radix - 1
 298                       /\ value v0n (n+n) = a1 * b1 + a0 * b0  - (a0 - a1)*(b0 - b1)
 299                          + hv0 + m * lvinf - m * cy2 + m * m
 300                       /\ !cy at AddSub = 0 }
 301     assert { !vm1_neg /\ value scratch (n+n) = - (a0-a1)*(b0-b1)
 302              \/ not !vm1_neg /\ value scratch (n+n) = (a0-a1)*(b0-b1)
 303              by value scratch (n+n) = asm1_abs * bsm1_abs
 304              so [@case_split]
 305                  (!vm1_neg at BSM1 /\ !vm1_neg
 306                \/ !vm1_neg at BSM1 /\ not !vm1_neg
 307                \/ not !vm1_neg at BSM1 /\ !vm1_neg
 308                \/ not !vm1_neg at BSM1 /\ not !vm1_neg) };
 309     assert { power radix (n+n) = m * m };
 310     if !vm1_neg
 311     then
 312       let c'' = wmpn_add_n_in_place v0n scratch (n+n) in
 313       assert { value v0n (n+n)
 314                = value v0n (n+n) at AddSub + value scratch (n+n)
 315                  - power radix (n+n) * c'' };
 316       cy := !cy + c'';
 317       assert { value v0n (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 318                  + hv0 + m * lvinf - m * cy2 - m * m * !cy
 319                by - m * m * c'' - m * m * !cy at AddSub = - m * m * !cy
 320                so value scratch (n+n) = -(a0-a1)*(b0-b1)
 321                so value v0n (n+n)
 322                   = value v0n (n+n) at AddSub + value scratch (n+n)
 323                     - power radix (n+n) * c''
 324                   = value v0n (n+n) at AddSub - (a0 - a1)*(b0-b1) - m * m * c''
 325                   = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 326                     + hv0 + m * lvinf - m * cy2 - m * m * !cy }
 327     else
 328       let b = wmpn_sub_n_in_place v0n scratch (n+n) in
 329       assert { value v0n (n+n)
 330                = value v0n (n+n) at AddSub - value scratch (n+n)
 331                  + power radix (n+n) * b };
 332       cy := sub_mod !cy b;
 333       assert { !cy <= 2 /\ !cy = !cy at AddSub - b
 334                \/ !cy = radix - 1 /\ !cy at AddSub = 0 /\ b = 1
 335                by [@case_split]
 336                   ((!cy at AddSub = 0 /\ b = 1
 337                    so !cy = EuclideanDivision.mod (-1) radix = radix - 1)
 338                    \/
 339                    (1 <= !cy at AddSub \/ b = 0
 340                    so 0 <= !cy at AddSub - b < radix
 341                    so !cy = !cy at AddSub - b)) };
 342       assert { !cy <= 2 ->
 343                (value v0n (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 344                                   + hv0 + m * lvinf - m * cy2 - m * m * !cy
 345                by !cy = !cy at AddSub - b
 346                so m * m * b - m * m * !cy at AddSub = - m * m * !cy
 347                so value scratch (n+n) = (a0-a1)*(b0-b1)
 348                so value v0n (n+n)
 349                   = value v0n (n+n) at AddSub - value scratch (n+n)
 350                     + power radix (n+n) * b
 351                   = value v0n (n+n) at AddSub - (a0 - a1)*(b0-b1) + m * m * b
 352                   = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 353                     + hv0 + m * lvinf - m * cy2 - m * m * !cy) };
 354       assert { !cy = radix - 1 ->
 355                (value v0n (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 356                                   - m * cy2 + hv0 + m * lvinf + m * m
 357                 by b = 1 so power radix (n+n) * b = m * m
 358                 so m * m * !cy at AddSub = 0
 359                 so value scratch (n+n) = (a0-a1)*(b0-b1)
 360                 so value v0n (n+n)
 361                    = value v0n (n+n) at AddSub - value scratch (n+n)
 362                      + power radix (n+n) * b
 363                    = value v0n (n+n) at AddSub - (a0 - a1)*(b0-b1) + m*m
 364                    = a1 * b1 + a0 * b0 + hv0 + m * lvinf
 365                      - m * cy2 - m * m * (!cy at AddSub) - (a0 - a1)*(b0-b1) + m*m
 366                    = a1 * b1 + a0 * b0 + hv0 - (a0 - a1)*(b0 - b1)
 367                      + m * lvinf - m * cy2 + m * m) }
 368   end;
 369   label Join in
 370   let ghost rj = { r } in
 371   let ghost v0nj = { v0n } in
 372   let ghost vinfnj = { vinfn } in
 373   join r v0n;
 374   value_sub_frame (pelts r) (pelts rj) (offset r) (offset r + p2i n);
 375   assert { value r n = value rj n = lv0 };
 376   value_concat r n (3 * n);
 377   value_sub_frame (pelts r) (pelts v0nj) (offset r + p2i n)
 378                                          (offset r + 3 * p2i n);
 379   assert { value r (3*n) = value r n + m * value (v0n at Join) (n+n)
 380            by offset v0nj = offset r + n
 381            so offset v0nj + (n + n) = offset r + 3 * n
 382            so value_sub (pelts r) (offset r + n) (offset r + 3*n)
 383               = value v0nj (n+n) };
 384   label JoinH in
 385   let ghost rh = { r } in
 386   join r vinfn;
 387   value_sub_frame (pelts r) (pelts rh) (offset r) (offset r + 3 * p2i n);
 388   assert { value r (3*n) = value r (3*n) at JoinH };
 389   value_sub_frame (pelts r) (pelts rh) (offset r) (offset r + p2i n);
 390   assert { value r n = value r n at JoinH };
 391   value_concat r (3*n) ((3 * n) + ((s + t) - n));
 392   assert { forall i. offset r + 3 * n <= i < offset r + 3 * n + s + t - n ->
 393            min vinfnj <= i < max vinfnj
 394            by max vinfnj >= offset r + sx + sy
 395               = offset r + n + s + n + t
 396               = offset r + 3 * n + s + t - n };
 397   value_sub_frame (pelts r) (pelts vinfnj) (offset r + 3 * p2i n)
 398                             (offset r + 3 * p2i n + p2i s + p2i t - p2i n);
 399   assert { value r (sx + sy)
 400            = value r (3*n + s + t - n)
 401            = value r (3*n) + m*m*m* value (vinfn at Join) (s+t-n)
 402            by offset vinfnj = offset r + 3*n
 403            so offset vinfnj + s + t - n = offset r + 3*n + s + t - n
 404            so m * m * m = power radix n * power radix n * power radix n
 405               = power radix (n+n+n) = power radix (3 * n)
 406            so value_sub (pelts r) (offset r + 3*n) (offset r + 3*n + s + t - n)
 407               = value vinfnj (s+t-n) };
 408   join scratch s_out;
 409   assert { a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1) = a0 * b1 + a1 * b0 };
 410   assert { value x sx * value y sy
 411            = (a0 +  m * a1)*(b0 + m * b1)
 412            = a0 * b0 + m * (a0 * b1 + a1 * b0) + m * m * (a1 * b1) };
 413   assert { !cy <= 3 /\ value r (sx + sy) = value x sx * value y sy
 414                                             - m * m * cy2 - m * m * m * !cy
 415            \/ !cy = radix - 1 /\ value r (sx + sy) = value x sx * value y sy
 416                                   - m * m * cy2 + m * m * m
 417             by value r n = lv0
 418                so value vinfnj (s+t-n) = hvinf
 419                so value r (sx + sy)
 420                = value r (3 * n) + m * m * m * value vinfnj (s+t-n)
 421                = value r n + m * value v0nj (n+n)
 422                  + m*m*m * value vinfnj (s+t-n)
 423                = lv0 + m * value v0nj (n+n) + m * m * m * hvinf
 424             so (lv0 + m * (a1*b1 + a0*b0 - (a0 - a1)*(b0 - b1) + hv0 + m*lvinf)
 425                 + m * m * m * hvinf = value x sx * value y sy
 426                by lv0 + m * (a0 * b1 + a1 * b0 + hv0 + m * lvinf)
 427                       + m * m * m * hvinf
 428                   = lv0 + m * hv0 + m * (a0 * b1 + a1 * b0)
 429                      + m * m * lvinf + m * m * m * hvinf
 430                   = lv0 + m * hv0 + m * (a0 * b1 + a1 * b0)
 431                     + m * m * (lvinf + m * hvinf)
 432                   = a0 * b0 + m * (a0 * b1 + a1 * b0)
 433                     + m * m * a1 * b1
 434                   = value x sx * value y sy)
 435             so [@case_split] (* TODO *)
 436                ((!cy <= 3
 437                 so value v0nj (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 438                                      + hv0 + m * lvinf - m * cy2 - m * m * !cy
 439                    = a0 * b1 + a1 * b0 + hv0 + m * lvinf - m * cy2 - m * m * !cy
 440                 so value r (sx + sy)
 441                    = lv0 + m * value v0nj (n+n) + m * m * m * hvinf
 442                    = lv0 + m * (a0 * b1 + a1 * b0 + hv0 + m * lvinf
 443                                - m * cy2 - m * m * !cy)
 444                                + m * m * m * hvinf
 445                    = lv0 + m * (a0 * b1 + a1 * b0 + hv0 + m * lvinf)
 446                       + m * m * m * hvinf
 447                       - m * m * cy2 - m * m * m * !cy
 448                    = value x sx * value y sy
 449                      - m * m * cy2 - m * m * m * !cy)
 450                  \/
 451                  (!cy = radix - 1
 452                   so value v0nj (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 453                                      + hv0 + m * lvinf - m * cy2 + m * m
 454                    = a0 * b1 + a1 * b0 + hv0 + m * lvinf - m * cy2 + m * m
 455                   so value r (sx + sy)
 456                    = lv0 + m * value v0nj (n+n) + m * m * m * hvinf
 457                    = lv0 + m * (a0 * b1 + a1 * b0 + hv0 + m * lvinf
 458                                - m * cy2 + m * m)
 459                                + m * m * m * hvinf
 460                    = lv0 + m * (a0 * b1 + a1 * b0 + hv0 + m * lvinf)
 461                       + m * m * m * hvinf
 462                       - m * m * cy2 + m * m * m
 463                    = value x sx * value y sy
 464                      - m * m * cy2 + m * m * m)) };
 465   let vinf0 = C.incr r (n+n) in
 466   value_sub_upper_bound (pelts x1) (offset x1) (offset x1 + int32'int s);
 467   value_sub_upper_bound (pelts y1) (offset y1) (offset y1 + int32'int t);
 468   assert { s + t - n > 0
 469            by s >= n-1
 470            so 4 * sx < 5 * sy
 471            so 4 * n + 4 * s < 5 * n + 5 * t
 472            so 5 * t > 4 * s - n >= 4 * n - 4 - n = 3 * n - 4
 473            so 5 * t > 3 * n - 4
 474            so n > 3
 475            so t > 1 };
 476   let ghost m' = power radix (p2i s + p2i t - p2i n) in
 477   begin ensures { value r (sx + sy) + m * m * cy2  < m * m * m * m' }
 478         ensures { value x sx * value y sy < m * m * m * m' }
 479   assert { power radix (s+t) = m * m' };
 480   value_sub_upper_bound (pelts r) (offset r) (offset r + int32'int n);
 481   value_sub_upper_bound (pelts v0nj) (offset v0nj)
 482                         (offset v0nj + int32'int n + int32'int n);
 483   value_sub_upper_bound (pelts x) (offset x) (offset x + int32'int sx);
 484   value_sub_upper_bound (pelts y) (offset y) (offset y + int32'int sy);
 485   assert { !cy = radix - 1 -> hvinf <= m' - 2
 486            by value v0nj (n+n) < power radix (n+n) = m * m
 487            so value v0nj (n+n) = a1 * b1 + a0 * b0 - (a0 - a1)*(b0 - b1)
 488                                      + hv0 + m * lvinf - m * cy2 + m * m
 489                    = a0 * b1 + a1 * b0 + hv0 + m * lvinf - m * cy2 + m * m
 490                               so (0 <= a0 * b1 by 0 <= a0 /\ 0 <= b1)
 491            so (0 <= a1 * b0 by 0 <= a1 /\ 0 <= b0)
 492            so hv0 >= 0
 493            so value v0nj (n+n) >= m * lvinf - m * cy2 + m * m
 494            so m * lvinf - m * cy2 + m * m < m * m
 495            so m * (lvinf - cy2) = m * lvinf - m * cy2 < 0
 496            so lvinf - cy2 < 0
 497            so (cy2 <= 1 by !cy at AddSub = 0 so !cy at Add3 = 0)
 498            so (lvinf = 0 by 0 <= lvinf)
 499            so a1 * b1 = m * hvinf
 500            so if a1*b1 <= 0
 501            then (hvinf <= m' - 2
 502                  by m*hvinf=0 so 0 <> m so hvinf = 0
 503                  so 1 <= s+t-n
 504                  so radix = power radix 1 <= power radix (s+t-n)
 505                  so m' > 2
 506                  so hvinf <= m' - 2)
 507            else hvinf <= m'-2 by
 508                 (forall s. 0 <= s -> power radix s = power 2 (64*s)
 509                 by radix = power 2 64)
 510            so (m = power 2 (64 * n)
 511                by m = power radix n)
 512            so m >= 1
 513            so (hvinf >= 1 by m*hvinf > 0)
 514            so valuation m 2 = 64*n
 515            so 64*n <= valuation (a1*b1) 2 = valuation a1 2 + valuation b1 2
 516            so if 64*n < valuation (a1*b1) 2
 517            then hvinf <= m'-2
 518                  by valuation (m*hvinf) 2 > 64*n
 519                  so valuation m 2 = 64*n
 520                  so if valuation hvinf 2 = 0
 521                     then false by valuation (m*hvinf) 2 = valuation m 2
 522                     else hvinf <= m'-2
 523                     by even hvinf
 524                     so (odd (m'-1)
 525                        by m' = power radix (s+t-n)
 526                              = power 2 (64*(s+t-n))
 527                              = 2 * power 2 (64*(s+t-n) - 1)
 528                           so even m')
 529                     so hvinf <> m'-1
 530                     so (hvinf < m'
 531                         by hvinf = value vinfnj (s+t-n)
 532                                  < power radix ((offset vinfnj +(s+t-n))
 533                                                  - offset vinfnj)
 534                                  = power radix (s+t-n))
 535                     so hvinf < m'-1
 536            else hvinf <= m'-2
 537            by power radix s = power 2 (64*s)
 538            so power radix t = power 2 (64*t)
 539            so m' = power radix (s+t-n) = power 2 (64*(s+t-n))
 540            so let k = valuation a1 2 in
 541               let l = valuation b1 2 in
 542               let a1' = div a1 (power 2 k) in
 543               let b1' = div b1 (power 2 l) in
 544               a1 = (power 2 k) * a1'
 545            so b1 = (power 2 l) * b1'
 546            so 64*n = k + l
 547            so 1 <= a1 /\ 1 <= b1
 548            so 0 <= k /\ 0 <= l
 549            so 1 <= a1' /\ 1 <= b1'
 550            so (k <= 64*s by power 2 k <= a1 < power 2 (64*s)
 551                          so power 2 k < power 2 (64*s)
 552                          so not 64*s < k)
 553            so (l <= 64*t by power 2 l <= b1 < power 2 (64*t)
 554                          so power 2 l < power 2 (64*t)
 555                          so not 64*t < l)
 556            so (forall a b c m:int. 0 <= a -> 0 <= b -> 0 < c -> 0 < m
 557               -> a * c < m -> b * c >= m -> a < b)
 558            so a1' < power 2 (64*s - k)
 559               (*(by a1' * (power 2 k) = a1 < power radix s = power 2 (64*s)
 560               so a1' * power 2 k < power 2 (64*s)
 561               so power 2 (64*s-k) * power 2 k = power 2 (64*s))*)
 562            so a1' <= power 2 (64*s - k) - 1
 563            so b1' < power 2 (64*t - l)
 564            so b1' <= power 2 (64*t - l) - 1
 565            so a1' * b1'
 566               <= (power 2 (64*s - k) - 1) * b1'
 567               <= (power 2 (64*s - k) - 1) * (power 2 (64*t - l) - 1)
 568               = (power 2 (64*s-k))*(power 2 (64*t -l))
 569                 - power 2 (64*s-k)
 570                 - power 2 (64*t-l)
 571                 + 1
 572               <= (power 2 (64*s-k))*(power 2 (64*t -l)) - 2
 573               = power 2 (64*(s+t) - (k+l)) - 2
 574               = power 2 (64*(s+t) - 64*n) - 2
 575               = power 2 (64*(s+t-n)) - 2
 576               = power radix (s+t-n) - 2
 577               = m' - 2
 578            so a1 * b1 = (power 2 k) * a1' * (power 2 l) * b1'
 579               = power 2 k * power 2 l * a1' * b1'
 580               = (power 2 (k+l)) * a1' * b1'
 581               = power 2 (64*n) * a1' * b1'
 582               = power radix n * a1' * b1'
 583               = m * (a1' * b1')
 584            so a1 * b1 = m * hvinf = m * (a1' * b1')
 585            so hvinf = a1' * b1' <= m' - 2 };
 586   assert { value x sx * value y sy < m * m * m * m'
 587            by (m * m * m * m' = power radix (n+n+s+t)
 588               by m * m * m * m'
 589                  = power radix n * power radix n * power radix n
 590                    * power radix (s+t-n)
 591                  = power radix n * power radix n * power radix (s+t)
 592                  = power radix n * power radix (n+s+t)
 593                  = power radix (n+n+s+t))
 594            so value x sx < power radix sx = power radix (n+s)
 595            so value y sy < power radix sy = power radix (n+t)
 596            so (forall a b c d. 0 <= a < b -> 0 <= c < d -> a * c < b * d
 597                by [@case_split]
 598                   (((a=0\/c=0) so a*c=0 so b*d>0 so a*c < b*d)
 599                   \/
 600                   ((0<a /\ 0<c) so a * c < a * d = d*a < d * b = b * d)))
 601            so value x sx * value y sy < power radix (n+s) * power radix (n+t)
 602               = power radix (n+n+s+t) };
 603   assert { value r (sx + sy) + m * m * cy2  < m * m * m * m'
 604            by [@case_split]
 605            ((!cy <= 3
 606            so 0 <= m * m * m * !cy
 607            so value r (sx + sy) + m * m * cy2
 608            = value x sx * value y sy - m * m * cy2 - m * m * m * !cy
 609                + m * m * cy2
 610            = value x sx * value y sy - m * m * m * !cy
 611            <= value x sx * value y sy)
 612            \/
 613            (!cy = radix - 1
 614            so hvinf <= m' - 2
 615            so value r (sx + sy) = value r n + m * (value v0nj (n+n))
 616                                    + m*m*m * value vinfnj (s+t-n)
 617               = value r n + m * (value v0nj (n+n)) + m*m*m * hvinf
 618            so value r n < power radix n
 619            so value r n < m
 620            so value v0nj (n+n) < power radix (n+n) = m * m
 621            so value v0nj (n+n) <= m * m - 1
 622            so m * value v0nj (n+n) <= m * (m * m - 1)
 623            so value r n + m * (value v0nj (n+n))
 624               < m + m * (m * m -1) = m * m * m
 625            so m * m * m * hvinf <= m * m * m * (m'-2)
 626            so value r (sx + sy) < m * m * m + m * m * m * (m'-2)
 627               = m * m * m * m' - m * m * m
 628            so m * m * cy2 < m * m * m)) };
 629   end;
 630   value_concat r (n + n) (sx + sy);
 631   assert { value_sub (pelts r) (offset r + n + n) (offset r + sx + sy)
 632            = value vinf0 (s+t) };
 633   assert { value r (sx + sy) = value r (n+n) + m * m * value vinf0 (s+t) };
 634   assert { value vinf0 (s+t) + cy2 < m * m'
 635            by value r (n+n) >= 0
 636            so m * m * m * m' > value r (sx + sy) + m * m * cy2
 637               = value r (n+n) + m * m * value vinf0 (s+t) + m * m * cy2
 638               >= m * m * value vinf0 (s+t) + m * m * cy2
 639               = m * m * (value vinf0 (s+t) + cy2)
 640               so (m * m) * (value vinf0 (s+t) + cy2) < (m * m) * (m * m') };
 641   let ghost ri = { r } in
 642   label IncrM in
 643   wmpn_incr vinf0 (s + t) cy2;
 644   value_concat r (n + n) (sx + sy);
 645   assert { value_sub (pelts r) (offset r + n + n) (offset r + sx + sy)
 646            = value vinf0 (s+t) };
 647   assert { value r (sx + sy) = value r (n+n) + m * m * value vinf0 (s+t) };
 648   assert { forall j. offset r <= j < offset r + (n+n)
 649            -> (pelts r)[j] = (pelts r)[j] at IncrM };
 650   value_sub_frame (pelts r) (pelts ri) (offset r)
 651                   (offset r + (int32'int n + int32'int n));
 652   assert { value r (sx + sy) = value r (sx+sy) at IncrM + m * m * cy2
 653            by value vinf0 (s+t) = value vinf0 (s+t) at IncrM + cy2
 654            so value r (n+n) = value r (n+n) at IncrM
 655            so value r (sx + sy) = value r (n+n) + m * m * value vinf0 (s+t)
 656               = (value r (n+n) at IncrM)
 657                 + m * m * (value vinf0 (s+t) at IncrM + cy2)
 658               = (value r (n+n) + m * m * value vinf0 (s+t) at IncrM)
 659                 + m * m * cy2
 660               = value r (sx+sy) at IncrM + m * m * cy2 };
 661   assert { forall j. min r <= j < offset r \/ offset r + sx + sy <= j < max r
 662              -> (pelts r)[j] = (pelts r)[j] at IncrM
 663                 by (pelts r)[j] = (pelts vinf0)[j]
 664                    = (pelts vinf0)[j] at IncrM };
 665   assert { !cy <= 3 /\ value r (sx + sy)
 666                        = value x sx * value y sy - m * m * m * !cy
 667            \/ value r (sx + sy) = value x sx * value y sy + m * m * m };
 668   let rh = { r } in
 669   let vinfn = C.incr r (3 * n) in
 670   label IncrH in
 671   assert { valid vinfn (s+t-n) };
 672   value_concat r (3 * n) (sx + sy);
 673   assert { value_sub (pelts r) (offset r + 3*n) (offset r + sx + sy)
 674            = value vinfn (s+t-n)
 675            by pelts r = pelts vinfn
 676            so offset r + 3*n = offset vinfn
 677            so offset r + sx + sy = offset vinfn + s + t - n };
 678   assert { value r (sx + sy) = value r (3*n)
 679                                + power radix (3*n) * value vinfn (s+t-n)};
 680   assert { power radix (3*n) = m * m * m
 681            by power radix (3*n) = power radix (n+n+n)
 682               = power radix (n+n) * power radix n
 683               = power radix n * power radix n * power radix n };
 684   if ([@likely] !cy <= 3)
 685   then begin
 686     assert { value r (sx+sy) = value x sx * value y sy
 687                                - power radix (3*n) * !cy
 688               by value r (sx+sy) = value x sx * value y sy - m * m * m * !cy };
 689     assert { value r (sx+sy) + power radix (3*n)* !cy < power radix (3*n) * m'
 690              by value r (sx+sy) + power radix (3*n) * !cy
 691                 = value x sx * value y sy
 692                 < m*m*m*m' = power radix (3*n) * m' };
 693     assert { value vinfn (s+t-n) + !cy < m'
 694              by power radix (3*n) * m'
 695                 > value r (sx+sy) + power radix (3*n) * !cy
 696                 = value r (3*n) + power radix (3*n) * value vinfn (s+t-n)
 697                                 + power radix (3*n) * !cy
 698                 = value r (3*n) + power radix (3*n) * (value vinfn (s+t-n) + !cy)
 699                 >= power radix (3*n) * (value vinfn (s+t-n) + !cy)
 700                 so power radix (3*n) * (value vinfn (s+t-n) + !cy)
 701                    < power radix (3*n) * m'};
 702     wmpn_incr vinfn ((s + t) - n) !cy;
 703     value_concat r (3 * n) (sx + sy);
 704     assert { value_sub (pelts r) (offset r + 3*n) (offset r + sx + sy)
 705              = value vinfn (s+t-n)
 706              by pelts r = pelts vinfn
 707              so offset r + 3*n = offset vinfn
 708              so offset r + sx + sy = offset vinfn + s + t - n };
 709     assert { value r (sx + sy) = value r (3*n)
 710                                 + power radix (3*n) * value vinfn (s+t-n)};
 711     assert { forall j. offset r <= j < offset r + (3*n)
 712              -> (pelts r)[j] = (pelts r)[j] at IncrH };
 713     value_sub_frame (pelts r) (pelts rh) (offset r)
 714                     (offset r + (3 * int32'int n));
 715     assert { value r (sx + sy) = value x sx * value y sy
 716              by value vinfn (s+t-n) = (value vinfn (s+t-n) at IncrH) + !cy
 717              so value r (3*n) = value r (3*n) at IncrH
 718              so value r (sx + sy)
 719                 = value r (3*n) + power radix (3*n) * value vinfn (s+t-n)
 720                 = value r (3*n) at IncrH
 721                   + power radix (3*n) * (value vinfn (s+t-n) at IncrH + !cy)
 722                 = (value r (3*n) + power radix (3*n) * value vinfn (s+t-n)
 723                          at IncrH)
 724                   + power radix (3*n) * !cy
 725                 = value r (sx+sy) at IncrH + power radix (3*n) * !cy
 726                 = value x sx * value y sy - power radix (3*n) * !cy
 727                                           + power radix (3*n) * !cy };
 728   end
 729   else begin
 730     assert { !cy = radix - 1 };
 731     assert { value r (sx+sy) = value x sx * value y sy
 732                                + power radix (3*n)
 733               by value r (sx+sy) = value x sx * value y sy + m * m * m };
 734     value_sub_upper_bound (pelts r) (offset r) (offset r + 3 * int32'int n);
 735     assert { 0 <= value vinfn (s+t-n) - 1
 736              by 0 <= value x sx so 0 <= value y sy
 737              so 0 <= value x sx * value y sy
 738              so value r (sx + sy) >= power radix (3*n)
 739              so value r (sx + sy)
 740                 = value r (3*n) + power radix (3*n) * value vinfn (s+t-n)
 741              so value r (3*n) < power radix (3*n)
 742              so power radix (3*n) * value vinfn (s+t-n) > 0
 743              so value vinfn (s+t-n) > 0 };
 744     wmpn_decr_1 vinfn ((s + t) - n);
 745     value_concat r (3 * n) (sx + sy);
 746     assert { value_sub (pelts r) (offset r + 3*n) (offset r + sx + sy)
 747              = value vinfn (s+t-n)
 748              by pelts r = pelts vinfn
 749              so offset r + 3*n = offset vinfn
 750              so offset r + sx + sy = offset vinfn + s + t - n };
 751     assert { value r (sx + sy) = value r (3*n)
 752                                 + power radix (3*n) * value vinfn (s+t-n)};
 753     assert { forall j. offset r <= j < offset r + (3*n)
 754              -> (pelts r)[j] = (pelts r)[j] at IncrH };
 755     value_sub_frame (pelts r) (pelts rh) (offset r)
 756                     (offset r + (3 * int32'int n));
 757     assert { value r (sx + sy) = value x sx * value y sy
 758              by value vinfn (s+t-n) = (value vinfn (s+t-n) at IncrH) - 1
 759              so value r (3*n) = value r (3*n) at IncrH
 760              so value r (sx + sy)
 761                 = value r (3*n) + power radix (3*n) * value vinfn (s+t-n)
 762                 = value r (3*n) at IncrH
 763                   + power radix (3*n) * (value vinfn (s+t-n) at IncrH - 1)
 764                 = (value r (3*n) + power radix (3*n) * value vinfn (s+t-n)
 765                          at IncrH)
 766                   - power radix (3*n)
 767                 = value r (sx+sy) at IncrH - power radix (3*n)
 768                 = value x sx * value y sy };
 769   end;
 770   label JoinRight in
 771   let rf = { r } in
 772   C.join r ro;
 773   label JoinLeft in
 774   C.join_r r' r;
 775   assert { forall j. offset r <= j < offset r + sx + sy ->
 776            (pelts r)[j] = (pelts rf)[j]
 777            by (pelts r)[j] = (pelts r)[j] at JoinLeft
 778               = (pelts rf)[j] };
 779   value_sub_frame (pelts r) (pelts rf) (offset r) (offset r + p2i sx + p2i sy);
 780   assert { value r (sx + sy) = value r (sx + sy) at JoinRight };
 781   C.join_r scratch' scratch
 782
 783
 784   (* Choose which multiplication algorithm is called recursively.
 785      Equivalent to the macros WMPN_TOOM22_MUL_REC and WMPN_TOOM22_MUL_N_REC
 786      respectively, in wmpn_toom22_mul.c *)
 787 with wmpn_toom22_mul_rec (r x: ptr limb) (sx:int32) (y:ptr limb) (sy: int32)
 788                          (scratch:ptr limb) (ghost k: int)
 789      : unit
 790   requires { valid x sx }
 791   requires { valid y sy }
 792   requires { valid r (sx + sy) }
 793   requires { writable r /\ writable scratch }
 794   requires { 0 < sy <= sx <= sy + sy }
 795   requires { 8 * sx < max_int32 }
 796   requires { 0 <= k <= 64 }
 797   requires { sx <= toom22_threshold * power 2 k }
 798   requires { valid scratch (2 * (sx + k)) }
 799   ensures  { value r (sx + sy) = value x sx * value y sy }
 800   ensures  { forall j. min r <= j < offset r \/ offset r + sx + sy <= j < max r
 801                        -> (pelts r)[j] = old (pelts r)[j] }
 802   ensures  { forall j. min scratch <= j < offset scratch
 803              -> (pelts scratch)[j] = old (pelts scratch)[j] }
 804   ensures  { min r = old min r }
 805   ensures  { max r = old max r }
 806   ensures  { plength r = old plength r }
 807   ensures  { min scratch = old min scratch }
 808   ensures  { max scratch = old max scratch }
 809   ensures  { plength scratch = old plength scratch }
 810   variant  { k + k + 1 }
 811 =
 812   if sy <= toom22_threshold
 813   then wmpn_mul_basecase r x sx y sy
 814   else
 815     if (4 * sx < 5 * sy) (* ? *)
 816     then wmpn_toom22_mul r x sx y sy scratch k
 817     else wmpn_toom32_mul r x sx y sy scratch k
 818
 819 with wmpn_toom22_mul_n_rec (r x y scratch: ptr limb) (sz:int32) (ghost k: int) : unit
 820   requires { valid x sz }
 821   requires { valid y sz }
 822   requires { valid r (sz + sz) }
 823   requires { writable r /\ writable scratch }
 824   requires { 0 < sz }
 825   requires { 8 * sz < max_int32 }
 826   requires { 0 <= k <= 64 }
 827   requires { sz <= toom22_threshold * power 2 k }
 828   requires { valid scratch (2 * (sz + k)) }
 829   ensures  { value r (sz + sz) = value x sz * value y sz }
 830   ensures  { forall j. min r <= j < offset r \/ offset r + sz + sz <= j < max r
 831                        -> (pelts r)[j] = old (pelts r)[j] }
 832   ensures  { forall j. min scratch <= j < offset scratch
 833              -> (pelts scratch)[j] = old (pelts scratch)[j] }
 834   ensures  { min r = old min r }
 835   ensures  { max r = old max r }
 836   ensures  { plength r = old plength r }
 837   ensures  { min scratch = old min scratch }
 838   ensures  { max scratch = old max scratch }
 839   ensures  { plength scratch = old plength scratch }
 840   variant  { k + k + 1 }
 841 =
 842   if sz <= toom22_threshold
 843   then wmpn_mul_basecase r x sz y sz
 844   else wmpn_toom22_mul r x sz y sz scratch k
 845
 846 with wmpn_toom32_mul (r x:ptr limb) (sx:int32) (y:ptr limb) (sy:int32)
 847                      (scratch: ptr limb) (ghost k: int) : unit
 848   requires { valid x sx }
 849   requires { valid y sy }
 850   requires { valid r (sx + sy) }
 851   requires { writable r /\ writable scratch }
 852   requires { toom22_threshold < sy }
 853   requires { 0 < k <= 64 }
 854   requires { sx <= toom22_threshold * power 2 k }
 855   requires { valid scratch (2 * (sx + k)) }
 856   requires { 8 * sx < max_int32 }
 857   requires { 4 < sy + 2 <= sx }
 858   requires { sx + 6 <= 3 * sy }
 859   ensures  { min r = old min r }
 860   ensures  { max r = old max r }
 861   ensures  { plength r = old plength r }
 862   ensures  { min scratch = old min scratch }
 863   ensures  { max scratch = old max scratch }
 864   ensures  { plength scratch = old plength scratch }
 865   ensures  { value r (sx + sy) = value x sx * value y sy }
 866   ensures  { forall j. min r <= j < offset r \/ offset r + sx + sy <= j < max r
 867                        -> (pelts r)[j] = old (pelts r)[j] }
 868   ensures  { forall j. min scratch <= j < offset scratch
 869                        -> (pelts scratch)[j] = old (pelts scratch)[j] }
 870   variant { k + k }
 871 =
 872   let n = 1 + (if 2 * sx >= 3 * sy
 873                then (sx - 1) / 3
 874                else (sy - 1) / 2) in
 875   let s = sx - 2 * n in
 876   let t = sy - n in
 877   assert { 0 < s <= n };
 878   assert { 0 < t <= n };
 879   assert { s + t >= n };
 880   let x0 = x in
 881   let x1 = C.incr x n in
 882   let x2 = C.incr x1 n in
 883   let y0 = y in
 884   let y1 = C.incr y n in
 885   let ghost a0 = value x0 (int32'int n) in
 886   let ghost a1 = value x1 (int32'int n) in
 887   let ghost a2 = value x2 (int32'int s) in
 888   let ghost b0 = value y0 (int32'int n) in
 889   let ghost b1 = value y1 (int32'int t) in
 890   let ghost m  = power radix (int32'int n) in
 891   value_concat x1 n (sx - n);
 892   value_concat x n sx;
 893   value_concat y n sy;
 894   assert { value y sy = b0 + m * b1 };
 895   assert { value x sx = a0 + m * a1 + m * m * a2 };
 896   let rol = decr_split r 0 in
 897   let ror = incr_split r (sx + sy) in
 898   let sol = decr_split scratch 0 in
 899   let sor = incr_split scratch ((n + n) + 1) in
 900   (* xp1 | yp1 | xm1 | ym1 *)
 901   let xp1 = r in                  (* x(1) = a0 + a1 + a2  *)
 902   let yp1 = incr_split r n in     (* y(1) = b0 + b1       *)
 903   let xm1 = incr_split yp1 n in   (* x(-1) = a0 - a1 + a2 *)
 904   let ym1 = incr_split xm1 n in   (* y(-1) = b0 - b1      *)
 905   let v1  = scratch in            (* x(1)*y(1)            *)
 906   let vm1 = r in                  (* x(-1)*y(-1)          *)
 907   let xp1_hi = ref 0 in           (* high limb of xp1     *)
 908   let yp1_hi = ref 0 in           (* high bit of yp1      *)
 909   let hi = ref 0 in               (* high bit of xm1      *)
 910   let vm1_neg = ref false in      (* negate vm1 ?         *)
 911   begin ensures { value xp1 n + m * !xp1_hi = a0 + a1 + a2 }
 912         ensures { (!vm1_neg /\ value xm1 n + m * !hi = a1 - (a0 + a2))
 913                   \/ (not !vm1_neg /\ value xm1 n + m * !hi = a0 - a1 + a2) }
 914         ensures { 0 <= !xp1_hi <= 2 }
 915         ensures { 0 <= !hi <= 1 }
 916     xp1_hi := wmpn_add xp1 x0 n x2 s;
 917     begin
 918       let cmp = wmpn_cmp xp1 x1 n in
 919       if (*begin ensures { result <-> a0 + a2 < a1 }*)
 920            !xp1_hi = 0
 921            && (cmp < 0)
 922          (*end*)
 923       then begin
 924         assert { value xp1 n < value x1 n };
 925         assert { value xp1 n = a0 + a2 };
 926         let ghost b = wmpn_sub_n xm1 x1 xp1 n in
 927         assert { b = 0 };
 928         hi := 0;
 929         vm1_neg := true end
 930       else begin
 931         let b = wmpn_sub_n xm1 xp1 x1 n in
 932         assert { !xp1_hi = 0 -> b = 0 };
 933         hi := !xp1_hi - b;
 934         assert { value xm1 n + m * !hi = a0 - a1 + a2
 935                  by value xm1 n + m * !hi
 936                   = value xm1 n - m * b + m * !xp1_hi
 937                   = value xp1 n - a1 + m * !xp1_hi
 938                   = a0 - a1 + a2 };
 939         end
 940     end;
 941     label Addx1 in
 942     assert { value xp1 n + m * !xp1_hi = a0 + a2 };
 943     let c = wmpn_add_n_in_place xp1 x1 n in
 944     xp1_hi := !xp1_hi + c;
 945     assert { value xp1 n + m * !xp1_hi = a0 + a2 + a1
 946              by value xp1 n + m * !xp1_hi
 947                 = value xp1 n + m * c + m * (!xp1_hi at Addx1)
 948                 = (value xp1 n at Addx1) + value x1 n + m * (!xp1_hi at Addx1)
 949                 = a0 + a2 + value x1 n = a0 + a2 + a1 };
 950   end;
 951   label B1 in
 952   begin ensures { value yp1 n + m * !yp1_hi = b0 + b1 }
 953         ensures { (!vm1_neg = (!vm1_neg at B1) /\ value ym1 n = b0 - b1)
 954                   \/ (!vm1_neg = not (!vm1_neg at B1) /\ value ym1 n = b1 - b0) }
 955         ensures { 0 <= !yp1_hi <= 1 }
 956     if (t = n)
 957     then begin
 958       yp1_hi := wmpn_add_n yp1 y0 y1 n;
 959       let cmp = wmpn_cmp y0 y1 n in
 960       if (cmp < 0)
 961       then begin
 962         let ghost b = wmpn_sub_n ym1 y1 y0 n in
 963         assert { b = 0 };
 964         vm1_neg := not !vm1_neg end
 965       else begin
 966         let ghost b = wmpn_sub_n ym1 y0 y1 n in
 967         assert { b = 0 }
 968         end
 969       end
 970     else begin
 971       yp1_hi := wmpn_add yp1 y0 n y1 t;
 972       let y0t = C.incr y0 t in
 973       let c0 = ((wmpn_zero_p y0t (n - t)) = 1) in
 974       let cmp = wmpn_cmp y0 y1 t in
 975       let c1 = (cmp < 0) in
 976       if c0 && c1
 977       then begin
 978         value_concat y0 t n;
 979         assert { value y0 t = value y0 n };
 980         assert { value y0 t < value y1 t };
 981         let ghost b = wmpn_sub_n ym1 y1 y0 t in
 982         assert { b = 0 };
 983         let ghost ym1z = { ym1 } in
 984         let ym1t = C.incr ym1 t in
 985         label Zero in
 986         wmpn_zero ym1t (n - t);
 987         assert { forall i. 0 <= i < t ->
 988                  (pelts ym1)[offset ym1 + i] = (pelts ym1z)[offset ym1z + i]
 989                  by offset ym1 + i < offset ym1t
 990                  so (pelts ym1)[offset ym1 + i]
 991                      = (pelts ym1t)[offset ym1 + i]
 992                      = (pelts ym1t at Zero)[offset ym1 + i]
 993                      = (pelts ym1 at Zero)[offset ym1 + i] };
 994         value_sub_frame_shift (pelts ym1) (pelts ym1z)
 995                               (offset ym1) (offset ym1z) (p2i t);
 996         assert { value ym1 t = value ym1 t at Zero };
 997         value_concat ym1 t n;
 998         assert { value ym1 n = value ym1 t
 999                  by value_sub (pelts ym1) (offset ym1 + t) (offset ym1 + n)
1000                     = value ym1t (n-t) = 0 };
1001         vm1_neg := not !vm1_neg end
1002       else begin
1003         value_concat y0 t n;
1004         assert { value y0 n = value y0 t + power radix t * value y0t (n-t) };
1005         assert { value y1 t <= value y0 n
1006                  by (not c0
1007                      so 1 <= value y0t (n - t)
1008                      so power radix t * 1 <= power radix t * value y0t (n-t)
1009                      so power radix t <= value y0 n
1010                      so value y1 t < power radix t)
1011                   \/ (not c1
1012                       so value y1 t <= value y0 t
1013                       so value y0 t <= value y0 n) };
1014         let ghost b = wmpn_sub ym1 y0 n y1 t in
1015         assert { b = 0 }
1016         end;
1017       end
1018   end;
1019   let ghost am1_abs = value xm1 (int32'int n) + m * (uint64'int !hi) in
1020   let ghost bm1_abs = value ym1 (int32'int n) in
1021   label RecP1 in
1022   wmpn_toom22_mul_n_rec v1 xp1 yp1 sor n (k-1);
1023   let cy = ref 0 in
1024   begin ensures { value scratch (2 * n) + power radix (n + n) * !cy
1025                   = (a0 + a1 + a2) * (b0 + b1) }
1026     assert { value scratch (n+n) = value xp1 n * value yp1 n };
1027     begin ensures { value scratch (n + n) + power radix (n + n) * !cy
1028                     = (a0 + a1 + a2) * (value yp1 n) }
1029           ensures { 0 <= !cy <= 3 }
1030           (* actually 2, but this is enough to prove there is no overflow *)
1031     if (!xp1_hi = 1)
1032     then begin
1033       let sa = { scratch } in
1034       let sn = C.incr scratch n in
1035       label Adjust1 in
1036       value_concat scratch n (n+n);
1037       assert { pelts scratch = pelts sn };
1038       let c = wmpn_add_n_in_place sn yp1 n in
1039       assert { pelts scratch = pelts sn };
1040       value_concat scratch n (n+n);
1041       value_sub_frame_shift (pelts scratch) (pelts sa)
1042                             (offset scratch) (offset sa) (p2i n);
1043       assert { value scratch n = value (scratch at Adjust1) n };
1044       assert { m * m = power radix (n+n) };
1045       assert { value scratch (n+n) + m * m * c = (a0 + a1 + a2) * value yp1 n
1046                by value sn n = value_sub (pelts scratch) (offset scratch + n)
1047                                          (offset scratch + (n+n))
1048                so value scratch (n+n) = value scratch n + m * value sn n
1049                so value (sn at Adjust1) n
1050                    = value_sub (pelts sa) (offset scratch + n)
1051                                           (offset scratch + (n+n))
1052                so value (scratch at Adjust1) (n+n)
1053                    = value scratch n + m * value (sn at Adjust1) n
1054                so value sn n + m * c = value (sn at Adjust1) n + value yp1 n
1055                so value scratch (n+n) + m * m * c
1056                   = value scratch n + m * (value sn n + m * c)
1057                   = value scratch n + m * value (sn at Adjust1) n + m * value yp1 n
1058                   = value (scratch at Adjust1) (n+n) + m * value yp1 n
1059                   = value xp1 n * value yp1 n + m * value yp1 n
1060                   = (value xp1 n + m * !xp1_hi) * value yp1 n
1061                   = (a0 + a1 + a2) * value yp1 n };
1062       cy := c end
1063     else begin
1064       if (!xp1_hi = 2)
1065       then begin
1066       let sa = { scratch } in
1067       let sn = C.incr scratch n in
1068       label Adjust2 in
1069       value_concat scratch n (n+n);
1070       assert { pelts scratch = pelts sn };
1071       let c = wmpn_addmul_1 sn yp1 n 2 in
1072       assert { pelts scratch = pelts sn };
1073       value_concat scratch n (n+n);
1074       value_sub_frame_shift (pelts scratch) (pelts sa)
1075                             (offset scratch) (offset sa) (p2i n);
1076       assert { value scratch n = value (scratch at Adjust2) n };
1077       assert { m * m = power radix (n+n) };
1078       assert { value scratch (n+n) + m * m * c = (a0 + a1 + a2) * value yp1 n
1079                by value sn n = value_sub (pelts scratch) (offset scratch + n)
1080                                          (offset scratch + (n+n))
1081                so value scratch (n+n) = value scratch n + m * value sn n
1082                so value (sn at Adjust2) n
1083                    = value_sub (pelts sa) (offset scratch + n)
1084                                           (offset scratch + (n+n))
1085                so value (scratch at Adjust2) (n+n)
1086                    = value scratch n + m * value (sn at Adjust2) n
1087                so value sn n + m * c = value (sn at Adjust2) n + !xp1_hi * value yp1 n
1088                so value scratch (n+n) + m * m * c
1089                   = value scratch n + m * (value sn n + m * c)
1090                   = value scratch n + m * value (sn at Adjust2) n
1091                                     + m * !xp1_hi * value yp1 n
1092                   = value (scratch at Adjust2) (n+n) + m * !xp1_hi * value yp1 n
1093                   = value xp1 n * value yp1 n + m * !xp1_hi * value yp1 n
1094                   = (value xp1 n + m * !xp1_hi) * value yp1 n
1095                   = (a0 + a1 + a2) * value yp1 n };
1096       assert { c <= 3
1097                by value sn n + m * c = value (sn at Adjust2) n + (value yp1 n) * 2
1098                so 0 <= value sn n
1099                so value (sn at Adjust2) n < m
1100                so value yp1 n < m
1101                so m * c < m * 3
1102                so c <= 3 };
1103       cy := c end
1104     else assert { !xp1_hi = 0 }  end
1105     end;
1106     begin ensures { value scratch (n + n) + power radix (n + n) * !cy
1107                     = (a0 + a1 + a2) * (b0 + b1) }
1108     if not (!yp1_hi = 0)
1109     then begin
1110       let sa = { scratch } in
1111       let sn = C.incr scratch n in
1112       label Adjust3 in
1113       value_concat scratch n (n+n);
1114       assert { pelts scratch = pelts sn };
1115       let c = wmpn_add_n_in_place sn xp1 n in
1116       value_concat scratch n (n+n);
1117       assert { pelts scratch = pelts sn };
1118       value_sub_frame_shift (pelts scratch) (pelts sa)
1119                             (offset scratch) (offset sa) (p2i n);
1120       assert { value scratch n = value (scratch at Adjust3) n };
1121       cy := !xp1_hi * !yp1_hi + c + !cy;
1122       assert { m * m = power radix (n+n) };
1123       assert { value scratch (n + n) + power radix (n + n) * !cy
1124                     = (a0 + a1 + a2) * (b0 + b1)
1125                by value sn n = value_sub (pelts scratch) (offset scratch + n)
1126                                          (offset scratch + (n+n))
1127                so value scratch (n+n) = value scratch n + m * value sn n
1128                so value (sn at Adjust3) n
1129                    = value_sub (pelts sa) (offset scratch + n)
1130                                           (offset scratch + (n+n))
1131                so value (scratch at Adjust3) (n+n)
1132                    = value scratch n + m * value (sn at Adjust3) n
1133                so !yp1_hi = 1
1134                so value sn n + m * c
1135                   = value (sn at Adjust3) n + !yp1_hi * value xp1 n
1136                so !yp1_hi * value xp1 n + m * !xp1_hi * !yp1_hi
1137                   = (a0 + a1 + a2) * !yp1_hi
1138                so value scratch (n+n) + m * m * !cy
1139                   = value scratch (n+n) + m * m * c
1140                        + m * m * !xp1_hi * !yp1_hi + m * m * (!cy at Adjust3)
1141                   = value scratch n + m * (value sn n + m * c)
1142                        + m * m * !xp1_hi * !yp1_hi + m * m * (!cy at Adjust3)
1143                   = value scratch n
1144                        + m * (value (sn at Adjust3) n + !yp1_hi * value xp1 n)
1145                        + m * m * !xp1_hi * !yp1_hi + m * m * (!cy at Adjust3)
1146                   = value scratch n + m * value (sn at Adjust3) n
1147                        + m * !yp1_hi * value xp1 n
1148                        + m * m * !xp1_hi * !yp1_hi + m * m * (!cy at Adjust3)
1149                   = value (scratch at Adjust3) (n+n)
1150                        + m * !yp1_hi * value xp1 n
1151                        + m * m * !xp1_hi * !yp1_hi + m * m * (!cy at Adjust3)
1152                   = (a0 + a1 + a2) * value yp1 n
1153                        + m * !yp1_hi * value xp1 n
1154                        + m * m * !xp1_hi * !yp1_hi
1155                   = (a0 + a1 + a2) * value yp1 n + m * (a0 + a1 + a2) * !yp1_hi
1156                   = (a0 + a1 + a2) * (value yp1 n + m * !yp1_hi)
1157                   = (a0 + a1 + a2) * (b0 + b1) } end
1158       else assert { b0 + b1 = value yp1 n }
1159     end;
1160   end;
1161   label RecM1 in
1162   join vm1 yp1;
1163   wmpn_toom22_mul_n_rec vm1 xm1 ym1 sor n (k-1);
1164   begin ensures { value vm1 (2*n) + m * m * !hi = am1_abs * bm1_abs }
1165         ensures { min r = old min r }
1166         ensures { max r = old max r }
1167         ensures { plength r = old plength r }
1168         ensures { 0 <= !hi <= 1 }
1169   if (not (!hi = 0))
1170   then begin
1171     assert { !hi = 1 };
1172     value_concat vm1 n (2*n);
1173     label HiSplit in
1174     let vm1n = incr_split vm1 n in
1175     let vm1l = { vm1 } in
1176     assert { value vm1 (2*n) at HiSplit = value vm1l n + m * value vm1n n
1177              by value vm1n n = value_sub (pelts vm1 at HiSplit)
1178                                (offset vm1 + n) (offset vm1 + 2*n) };
1179     label HiAdd in
1180     let c = wmpn_add_n_in_place vm1n ym1 n in
1181     label HiJoin in
1182     let vm1nj = { vm1n } in
1183     join vm1 vm1n;
1184     value_concat vm1 n (2*n);
1185     value_sub_frame (pelts vm1) (pelts vm1l) (offset vm1) (offset vm1 + p2i n);
1186     value_sub_frame (pelts vm1) (pelts vm1nj) (offset vm1 + p2i n)
1187                                               (offset vm1 + 2 * p2i n);
1188     assert { value vm1 n = value vm1l n };
1189     assert { value vm1 (2*n) = value vm1l n + m * value vm1nj n
1190              by value vm1nj n = value_sub (pelts vm1) (offset vm1 + p2i n)
1191                                                   (offset vm1 + 2 * p2i n) };
1192     assert { value vm1 (2*n) + m * m * c = am1_abs * bm1_abs
1193              by value vm1 (2 * n) + m * m * c
1194                 = value vm1l n + m * (value (vm1n at HiJoin) n + m * c)
1195                 = value vm1l n + m * (value (vm1n at HiAdd) n + value ym1 n)
1196                 = value vm1 (2*n) at HiSplit + m * value ym1 n
1197                 = value xm1 n * value ym1 n + m * value ym1 n
1198                 = (value xm1 n + m * !hi) * value ym1 n
1199                 = am1_abs * bm1_abs };
1200     hi := c;
1201     end
1202   end;
1203   let ghost vx0 = a0 * b0 in
1204   let ghost vx1 = a1 * b0 + a0 * b1 in
1205   let ghost vx2 = a2 * b0 + a1 * b1 in
1206   let ghost vx3 = a2 * b1 in
1207   assert { 0 <= vx0 /\ 0 <= vx1 /\ 0 <= vx2 /\ 0 <= vx3
1208            by 0 <= a0 /\ 0 <= a1 /\ 0 <= a2 /\ 0 <= b0 /\ 0 <= b1 };
1209   assert { (a0 + m * a1 + m * m * a2) * (b0 + m * b1)
1210            = vx0 + m * vx1 + m * m * vx2 + m * m * m * vx3 };
1211   assert { (a0 + a1 + a2) * (b0 + b1) = vx0 + vx1 + vx2 + vx3 };
1212   assert { vx0 + vx2 < 3 * m * m
1213             by (0 <= a0 < m /\ 0 <= a1 < m /\ 0 <= a2 < m
1214                 /\ 0 <= b0 < m /\ 0 <= b1 < m
1215                 by a2 < power radix s <= m
1216                 so b1 < power radix t <= m)
1217             so (vx0 < m * m
1218                 by vx0 = a0 * b0 <= a0 * m = m * a0 < m * m)
1219             so (vx2 < 2 * m * m
1220                 by vx2 = a2 * b0 + a1 * b1
1221                 so a2 * b0 <= a2 * m = m * a2 < m * m
1222                 so a1 * b1 <= a1 * m = m * a1 < m * m) };
1223   begin ensures { value scratch (2*n + 1) = vx0 + vx2 }
1224         ensures { if !vm1_neg
1225                   then am1_abs*bm1_abs = (a0 - a1 + a2) * (b1 - b0)
1226                   else am1_abs*bm1_abs = (a0 - a1 + a2) * (b0 - b1) }
1227     begin ensures { value scratch (2*n + 1) = 2 * (vx0 + vx2) }
1228           ensures { if !vm1_neg
1229                     then am1_abs*bm1_abs = (a0 - a1 + a2) * (b1 - b0)
1230                     else am1_abs*bm1_abs = (a0 - a1 + a2) * (b0 - b1) }
1231     if !vm1_neg
1232     then begin
1233       assert { am1_abs * bm1_abs = (a0 - a1 + a2) * (b1 - b0)
1234                by if !vm1_neg at B1
1235                then am1_abs * bm1_abs = (a0 - a1 + a2) * (b1 - b0)
1236                     by am1_abs = a1 - (a0 + a2) /\ bm1_abs = b0 - b1
1237                else am1_abs * bm1_abs = (a0 - a1 + a2) * (b1 - b0)
1238                     by am1_abs = a0 - a1 + a2 /\ bm1_abs = b1 - b0 };
1239       assert { (a0 - a1 + a2) * (b1 - b0) = vx1 -vx0 + vx3 - vx2 };
1240       label Sub in
1241       let b = wmpn_sub_n_in_place scratch vm1 (2*n) in
1242       let r, ghost b' = sub_with_borrow !cy !hi b in
1243       assert { (b' = 0 /\ value scratch (2*n) + m * m * r = 2 * (vx0 + vx2))
1244                by r - radix * b' = !cy - !hi - b
1245                so m * m * r - m * m * radix * b'
1246                   = m * m * !cy - m * m * !hi - m * m * b
1247                so (m * m * b
1248                     = value scratch (2*n) - value (scratch at Sub) (2*n)
1249                       + value vm1 (2*n)
1250                    by value scratch (2*n) - (m * m) * b
1251                       = value (scratch at Sub) (2*n) - value vm1 (2*n))
1252                so m * m * r - m * m * radix * b' + value scratch (2 * n)
1253                   = m * m * !cy - m * m * !hi + value (scratch at Sub) (2*n)
1254                     - value vm1 (2*n)
1255                   = (value (scratch at Sub) (2*n) + m * m * !cy)
1256                     - (value vm1 (2*n) + m * m * !hi)
1257                   = (vx0 + vx1 + vx2 + vx3) - (vx1 - vx0 + vx3 - vx2)
1258                   = 2 * (vx0 + vx2)
1259                   >= 0
1260                so m * m * r - m * m * radix * b' + value scratch (2 * n) >= 0
1261                so value scratch (2 * n) < m * m
1262                so (m * m) * r <= (m * m) * (radix - 1)
1263                so value scratch (2 * n) + m * m * r < m * m * radix
1264                so m * m * radix * (1 - b') > 0
1265                so b' < 1
1266                so b' = 0
1267                };
1268       label Set in
1269       value_sub_shift_no_change (pelts scratch) (offset scratch)
1270                                 (2 * p2i n)
1271                                 (2 * p2i n) !cy;
1272       set_ofs scratch (2*n) r;
1273       assert { value scratch (2 * n) = value (scratch at Set) (2 * n) };
1274       value_tail scratch (2*n);
1275       assert { value scratch (2*n + 1) = 2 * (vx0 + vx2)
1276                by value scratch (2*n + 1) = value scratch (2*n) + (m*m) * r };
1277       end
1278     else begin
1279       assert { am1_abs * bm1_abs = (a0 - a1 + a2) * (b0 - b1)
1280                by if !vm1_neg at B1
1281                then am1_abs * bm1_abs = (a0 - a1 + a2) * (b0 - b1)
1282                     by am1_abs = a1 - (a0 + a2) /\ bm1_abs = b1 - b0
1283                else am1_abs * bm1_abs = (a0 - a1 + a2) * (b0 - b1)
1284                     by am1_abs = a0 - a1 + a2 /\ bm1_abs = b0 - b1 };
1285       assert { (a0 - a1 + a2) * (b0 - b1) = vx0 - vx1 + vx2 - vx3 };
1286       label Add in
1287       let c = wmpn_add_n_in_place scratch vm1 (2*n) in
1288       let r, ghost c' = add_with_carry !cy !hi c in
1289       assert { (c' = 0 /\ value scratch (2*n) + (m*m) * r = 2 * (vx0 + vx2))
1290                 by r + radix * c' = !cy + !hi + c
1291                 so m * m * r + m * m * radix * c'
1292                    = m * m * !cy + m * m * !hi + m * m * c
1293                 so (m * m * c
1294                     = value (scratch at Add) (2*n) - value scratch (2*n)
1295                       + value vm1 (2*n)
1296                     by value scratch (2*n) + (m*m)*c
1297                        = value (scratch at Add) (2*n) + value vm1 (2*n))
1298                 so m * m * r + m * m * radix * c' + value scratch (2*n)
1299                    = m * m * !cy + m * m * !hi + value (scratch at Add) (2*n)
1300                      + value vm1 (2*n)
1301                    = (value (scratch at Add) (2*n) + m * m * !cy)
1302                      + (value vm1 (2*n) + m * m * !hi)
1303                    = (vx0 + vx1 + vx2 + vx3) + (am1_abs * bm1_abs)
1304                    = (vx0 + vx1 + vx2 + vx3) + (vx0 - vx1 + vx2 - vx3)
1305                    = 2 * (vx0 + vx2)
1306                 so 2 * (vx0 + vx2) < 6 * m * m < m * m * radix
1307                 so (m * m * radix * c' < m * m * radix
1308                     by 0 <= r so 0 <= m * m * r
1309                     so 0 <= value scratch (2*n)
1310                     so m * m * radix * c'
1311                        <= m * m * r + m * m * radix * c' + value scratch (2*n)
1312                        < m * m * radix)
1313                 so c' < 1
1314                 so c' = 0
1315       };
1316       label Set in
1317       value_sub_shift_no_change (pelts scratch) (offset scratch)
1318                                 (2 * p2i n)
1319                                 (2 * p2i n) !cy;
1320       set_ofs scratch (2*n) r;
1321       assert { value scratch (2 * n) = value (scratch at Set) (2 * n) };
1322       value_tail scratch (2*n);
1323       assert { value scratch (2*n + 1) = 2 * (vx0 + vx2)
1324                by value scratch (2*n + 1) = value scratch (2*n) + (m*m) * r };
1325     end
1326     end;
1327     label Shift in
1328     let s = (2 * n) + 1 in
1329     let ghost low = wmpn_rshift scratch scratch s 1 in
1330     assert { low = 0 /\ value scratch s = vx0 + vx2
1331              by (low + radix * value scratch s)
1332                 = value (scratch at Shift) s * power 2 (Limb.length - 1)
1333                 = 2 * (vx0 + vx2) * power 2 (Limb.length - 1)
1334                 = (vx0 + vx2) * power 2 Limb.length
1335                 = (vx0 + vx2) * radix
1336              so divides radix ((vx0 + vx2) * radix)
1337              so divides radix (low + radix * value scratch s)
1338              so divides radix low }
1339   end;
1340   let ghost vy = vx1 + vx3 + (vx0 + vx2) * m in
1341   assert { vy = (vx0 + vx2) * m + vx0 + vx2 - (vx0 - vx1 + vx2 - vx3) };
1342   join xm1 ym1;
1343   (* (    r    |    xm1  ) *)
1344   let ghost ss = (2 * n) + 1 in
1345   assert { value scratch ss = vx0 + vx2 };
1346   let vy0 = scratch in
1347   let ghost l02 = value scratch (int32'int n) in
1348   let vy1 = xm1 in
1349   let vy2 = incr_split scratch n in
1350   let ghost h02 = value vy2 (int32'int n) in
1351   let t02 = get_ofs vy2 n in
1352   begin ensures { value vy0 n + m * value vy1 n + m * m * value vy2 (n+1) = vy }
1353         ensures { value vy2 (n+1) < (power radix n) * 6 }
1354   label Vy in
1355   let os = { vy0 } in
1356   value_tail vy2 n;
1357   value_concat scratch n ss;
1358   assert { l02 + m * h02 + m * m * t02 = vx0 + vx2
1359            by vx0 + vx2 = value scratch ss
1360               = value scratch n + m * value_sub (pelts scratch)
1361                                       (offset scratch + n) (offset scratch + ss)
1362               = value scratch n + m * value vy2 (n+1)
1363               = l02 + m * (h02 + m * t02) };
1364   assert { t02 < 3
1365            by l02 + m * h02 + m * m * t02 < m * m * 3
1366            so 0 <= l02 /\ 0 <= h02 /\ 0 <= m
1367            so 0 <= l02 + m * h02
1368            so m * m * t02 < m * m * 3 };
1369   let c = wmpn_add_n vy1 vy0 vy2 n in
1370   assert { value vy2 (n+1) < (power radix n) * 4
1371            by value vy2 (n+1) = value vy2 n + (power radix n) * t02
1372               < power radix n + (power radix n) * t02
1373               <= (power radix n) * 4 };
1374   assert { value vy2 (n+1) + (c + t02) < (power radix n) * 5
1375            by c + t02 <= power radix n
1376            so value vy2 (n+1) + (c+t02) <= (power radix n) * 5 };
1377   wmpn_incr vy2 (n+1) (c + t02);
1378   assert { value vy2 (n+1) < (power radix n) * 5 };
1379   value_sub_frame (pelts vy0) (pelts os) (offset scratch) (offset scratch + int32'int n);
1380   assert { value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)
1381            = l02 + m * (l02 + h02) + m * m * (h02 + t02 + m * t02)
1382            = (vx0 + vx2) * m + vx0 + vx2
1383            by value vy0 n = l02
1384            so value vy1 n + m * c = l02 + h02
1385            so value vy2 (n+1) = value vy2 (n+1) at Vy + c + t02
1386               = h02 + t02 + m * t02 + c
1387            so value vy0 n + m * value vy1 n + m * m * value vy2 (n+1) + m * c
1388               = l02 + m * (l02 + h02) + m * m * (h02 + t02 + m * t02) + m * c };
1389   let vm1n = incr vm1 n in
1390   value_concat vm1 n (2*n);
1391   assert { value vm1 n + m * value vm1n n + m * m * !hi = am1_abs * bm1_abs
1392            by value vm1 n + m * value vm1n n = value vm1 (2*n) };
1393   begin ensures { value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)
1394                   = old (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1))
1395                     - (vx0 - vx1 + vx2 - vx3) }
1396         ensures { value vy2 (n+1) < (power radix n) * 6 }
1397   label Vm1 in
1398   let ovy2 = { vy2 } in
1399   if !vm1_neg
1400   then begin
1401     assert { am1_abs*bm1_abs = (a0 - a1 + a2) * (b1 - b0)
1402              = - (vx0 - vx1 + vx2 - vx3) };
1403     let c1 = wmpn_add_n_in_place scratch vm1 n in
1404     assert { value scratch n = value scratch n at Vm1 + value vm1 n - m * c1 };
1405     let c2 = wmpn_add_n_in_place vy1 vm1n n in
1406     assert { value vy1 n = value vy1 n at Vm1 + value vm1n n - m * c2};
1407     hi := !hi + c2;
1408     let c3 = wmpn_add_1_in_place vy1 n c1 in
1409     assert { value vy1 n
1410              = value vy1 n at Vm1 + value vm1n n + c1 - m * (c2 + c3) };
1411     hi := !hi + c3;
1412     wmpn_incr vy2 (n+1) !hi;
1413     assert { value vy2 (n+1) = value ovy2 (n+1) + c2 + c3 + !hi at Vm1 };
1414     assert { value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)
1415              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1416                - (vx0 - vx1 + vx2 - vx3)
1417              by value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)
1418              = value vy0 n at Vm1 + value vm1 n - m * c1
1419                + m * (value vy1 n at Vm1 + value vm1n n + c1 -  m * (c2 + c3))
1420                + m * m * (value ovy2 (n+1) + c2 + c3 + !hi at Vm1)
1421              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1422                + (value vm1 n + m * value vm1n n + m * m * !hi at Vm1)
1423              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1424                + am1_abs * bm1_abs
1425              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1426                 - (vx0 - vx1 + vx2 - vx3) };
1427     end
1428   else begin
1429     assert { am1_abs*bm1_abs = (a0 - a1 + a2) * (b0 - b1)
1430              = vx0 - vx1 + vx2 - vx3 };
1431     let b1 = wmpn_sub_n_in_place scratch vm1 n in
1432     assert { value scratch n = value scratch n at Vm1 - value vm1 n + m * b1 };
1433     let b2 = wmpn_sub_n_in_place vy1 vm1n n in
1434     assert { value vy1 n = value vy1 n at Vm1 - value vm1n n + m * b2 };
1435     hi := !hi + b2;
1436     let b3 = wmpn_sub_1_in_place vy1 n b1 in
1437     assert { value vy1 n
1438              = value vy1 n at Vm1 - value vm1n n - b1 + m * (b2 + b3) };
1439     hi := !hi + b3;
1440     assert { value vy0 n + m * value vy1 n + m * m * (value vy2 (n+1) - !hi)
1441              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1442                - (vx0 - vx1 + vx2 - vx3)
1443              by value vy0 n + m * value vy1 n + m * m * (value vy2 (n+1) - !hi)
1444              = value vy0 n + m * value vy1 n
1445                + m * m * (value ovy2 (n+1) - b2 - b3 - !hi at Vm1)
1446              = value vy0 n at Vm1 - value vm1 n + m * b1
1447                + m * (value vy1 n at Vm1 - value vm1n n - b1 +  m * (b2 + b3))
1448                + m * m * (value ovy2 (n+1) - b2 - b3 - !hi at Vm1)
1449              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1450                - (value vm1 n + m * value vm1n n + m * m * !hi at Vm1)
1451              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1452                - am1_abs * bm1_abs
1453              = (value vy0 n + m * value vy1 n + m * m * value vy2 (n+1)) at Vm1
1454                - (vx0 - vx1 + vx2 - vx3) };
1455     assert { value vy2 (n+1) - !hi >= 0
1456              by value vy0 n + m * value vy1 n + m * m * (value vy2 (n+1) - !hi)
1457                 = vy
1458              so (vy >= 0 by 0 <= vx1 /\ 0 <= vx1 /\ 0 <= vx2 /\ 0 <= vx3)
1459              so value vy0 n < m
1460              so value vy1 n <= m - 1
1461              so m * value vy1 n <= m * (m-1)
1462              so value vy0 n + m * value vy1 n <= value vy0 n + m * (m-1)
1463                 < m + m * (m-1) = m * m
1464              so 0 <= value vy0 n + m * value vy1 n + m * m * (value vy2 (n+1) - !hi)
1465                 < m * m * (1 + value vy2 (n+1) - !hi)
1466              so 0 < (m * m) * (1 + value vy2 (n+1) - !hi)
1467              so 0 < (1 + value vy2 (n+1) - !hi) };
1468     wmpn_decr vy2 (n+1) !hi;
1469     assert { value vy2 (n+1) = value ovy2 (n+1) - b2 - b3 - !hi at Vm1 };
1470     end;
1471   end
1472   end;
1473   label Split3 in
1474   let ghost ovy1 = { vy1 } in
1475   wmpn_toom22_mul_n_rec r x y sor n (k-1);
1476   let r3n = incr_split xm1 n in
1477   value_sub_frame (pelts vy1) (pelts ovy1) (offset vy1)
1478                               (offset vy1 + int32'int n);
1479   assert { value vy1 n = value ovy1 n };
1480   assert { value r (n+n) = vx0
1481            by value x n = a0 so value y n = b0 };
1482   begin
1483     ensures { value r3n (s+t) = vx3 }
1484     ensures { min r3n = old min r3n /\ max r3n = old max r3n }
1485     ensures { plength r3n = old plength r3n }
1486     if (s > t)
1487     then let ghost _ = wmpn_mul r3n x2 s y1 t (k-1) in ()
1488     else let ghost _ = wmpn_mul r3n y1 t x2 s (k-1) in ()
1489   end;
1490   assert { (a0 + m * a1 + m * m * a2) * (b0 + m * b1)
1491            = m * vy + vx0 + m * m * m * vx3 - m * m * vx0 - m * vx3 };
1492   assert { 0 <= vx0 < m * m /\ 0 <= vx3 < m * m
1493            by a0 * b0 <= a0 * m = m * a0 < m * m
1494            so a2 < power radix s <= m
1495            so b1 < power radix t <= m
1496            so a2 * b1 <= a2 * m = m * a2 < m * m };
1497   let ghost or = { r } in
1498   let r1n = incr_split r n in
1499   value_sub_frame (pelts r) (pelts or) (offset r) (offset r + int32'int n);
1500   value_sub_frame (pelts r1n) (pelts or) (offset r1n) (offset r1n + int32'int n);
1501   let ghost lvx0 = value r (int32'int n) in
1502   let ghost hvx0 = value r1n (int32'int n) in
1503   value_concat or n (n + n);
1504   assert { vx0 = lvx0 + m * hvx0 };
1505   let ghost or3n = { r3n } in
1506   let r4n = incr_split r3n n in
1507   let r2n = xm1 in
1508   value_sub_frame (pelts r3n) (pelts or3n) (offset r3n)
1509                                            (offset r3n + int32'int n);
1510   value_sub_frame (pelts r4n) (pelts or3n) (offset r4n)
1511                   (offset r4n + int32'int s + int32'int t - int32'int n);
1512   let ghost lvx3 = value r3n (int32'int n) in
1513   let ghost hvx3 = value r4n (int32'int s + int32'int t- int32'int n) in
1514   value_concat or3n n (s + t);
1515   assert { vx3 = lvx3 + m * hvx3 };
1516   let ghost vvy0 = value vy0 (int32'int n) in
1517   let ghost vvy1 = value vy1 (int32'int n) in
1518   let ghost vvy2 = value vy2 (int32'int n+1) in
1519   assert { vy = vvy0 + m * vvy1 + m * m * vvy2 };
1520   assert { m * vy + vx0 + m * m * m * vx3 - m * m * vx0 - m * vx3
1521            = lvx0 + m * (vvy0 + hvx0 - lvx3) + m * m * (vvy1 - lvx0 - hvx3)
1522              + m * m * m * (vvy2 - (hvx0 - lvx3)) + m * m * m * m * hvx3 };
1523   label R1 in
1524   let bo1 = wmpn_sub_n_in_place r1n r3n n in
1525   let bo = ref bo1 in
1526   assert { value r1n n - m * bo1 = hvx0 - lvx3 };
1527   let ly2 = get_ofs vy2 n in
1528   value_tail vy2 n;
1529   assert { 0 <= ly2 < 6
1530            by 0 <= value vy2 n
1531            so value vy2 n + (power radix n) * ly2 = value vy2 (n+1)
1532            so (power radix n) * ly2 <= value vy2 (n+1) < (power radix n) * 6 };
1533   let h = ref (Limb.to_int64 (ly2 + bo1)) in
1534   label R2 in
1535   let bo2 = wmpn_sub_n_in_place r2n r n in
1536   let bo2' = wmpn_sub_1_in_place r2n n !bo in
1537   bo := bo2 + bo2';
1538   assert { value r2n n - m * !bo = vvy1 - lvx0 - (!bo at R2) };
1539   assert { value r1n n + m * value r2n n - m * m * !bo
1540            = hvx0 - lvx3 + m * (vvy1 - lvx0) };
1541   label R3 in
1542   let bo3 = wmpn_sub_n r3n vy2 r1n n in
1543   let bo3' = wmpn_sub_1_in_place r3n n !bo in
1544   bo := bo3 + bo3';
1545   assert { value r3n n - m * !bo = value vy2 n - value r1n n - (!bo at R3) };
1546   assert { value r1n n + m * value r2n n + m * m * value r3n n - m * m *m * !bo
1547            = hvx0 - lvx3 + m * (vvy1 - lvx0) + m * m * value vy2 n
1548              - m * m * (hvx0 - lvx3 + m * bo1) };
1549   h := Int64.(-) !h (Limb.to_int64 !bo);
1550   label Join3 in
1551   let ghost or1n = { r1n } in
1552   let ghost or2n = { r2n } in
1553   let ghost or3n = { r3n } in
1554   join r2n r3n;
1555   value_concat r2n n (n + n);
1556   join r1n r2n;
1557   value_sub_frame (pelts r1n) (pelts or1n)
1558                   (offset r1n) (offset r1n + int32'int n);
1559   value_sub_frame (pelts r1n) (pelts or2n)
1560                   (offset r1n + int32'int n) (offset r1n + 2 * int32'int n);
1561   value_sub_frame (pelts r1n) (pelts or3n)
1562                   (offset r1n + 2 * int32'int n) (offset r1n + 3 * int32'int n);
1563   value_concat r1n n (3 * n);
1564   value_sub_concat (pelts r1n) (offset r1n + int32'int n)
1565                    (offset r1n + int32'int (2 * n))
1566                    (offset r1n + int32'int (3 * n));
1567   assert { value r1n (3*n)
1568            = value r1n n + m * value r2n n + m * m * value r3n n at Join3
1569            by offset r2n at Join3 = offset r1n + n
1570            so offset r3n at Join3 = offset r1n + 2 * n
1571            so value r1n n = value r1n n at Join3
1572            so value_sub (pelts r1n) (offset r1n + n) (offset r1n + 2*n)
1573               = value r2n n at Join3
1574            so value_sub (pelts r1n) (offset r1n + 2*n) (offset r1n + 3*n)
1575               = value_sub (pelts or3n) (offset r1n + 2*n) (offset r1n + 3*n)
1576               = value r3n n at Join3
1577            so value r1n (3*n)
1578               = value r1n n
1579                 + m * value_sub (pelts r1n) (offset r1n + n) (offset r1n + 3*n)
1580               = value r1n n
1581                 + m *
1582                 (value_sub (pelts r1n) (offset r1n + n) (offset r1n + 2*n)
1583                  + m * value_sub (pelts r1n) (offset r1n + 2*n)
1584                                              (offset r1n + 3*n))
1585               = value r1n n
1586                 + m * value_sub (pelts r1n) (offset r1n + n) (offset r1n + 2*n)
1587                 + m * m * value_sub (pelts r1n) (offset r1n + 2*n)
1588                                                 (offset r1n + 3*n) };
1589   assert { value r1n (3*n) + m * m * m * !h
1590            = hvx0 - lvx3 + m * (vvy1 - lvx0) + m * m * vvy2
1591              - m * m * (hvx0 - lvx3)
1592            by value r1n (3*n) - m * m * m * !bo
1593               = hvx0 - lvx3 + m * (vvy1 - lvx0) + m * m * value vy2 n
1594                - m * m * (hvx0 - lvx3 + m * bo1)
1595            so !h = ly2 + bo1 - !bo
1596            so value r1n (3*n) + m * m * m * !h
1597               = value r1n (3*n) - m * m * m * !bo
1598                 + m * m * m * ly2 + m * m * m * bo1
1599            so m * m * value vy2 n + m * m * m * ly2 = m * m * vvy2 };
1600   label Addy0 in
1601   let c = wmpn_add_in_place r1n (3 * n) scratch n in
1602   h := Int64.(+) !h (Limb.to_int64 c);
1603   assert { power radix (3*n) = m * m * m };
1604   assert { value r1n (3*n) + m * m * m * !h
1605            = vy + hvx0 - lvx3 - m * lvx0 - m * m * (hvx0 - lvx3)
1606            = vy + hvx0 - lvx3 - m * vx0 + m * m * lvx3
1607            by value r1n (3*n) + m * m * m * !h
1608               = value r1n (3*n) + (m * m * m * !h at Addy0) + m * m * m * c
1609               = (value r1n (3*n) + m * m * m * !h at Addy0)
1610                  + vvy0
1611               = vvy0 + hvx0 - lvx3 + m * (vvy1 - lvx0) + m * m * vvy2
1612                 - m * m * (hvx0 - lvx3) };
1613   label Join1 in
1614   let ghost or = { r } in
1615   let ghost or1n = { r1n } in
1616   join r r1n;
1617   value_sub_frame (pelts r) (pelts or) (offset r) (offset r + int32'int n);
1618   value_sub_frame (pelts r) (pelts or1n)
1619                   (offset r + int32'int n)
1620                   (offset r + 4 * int32'int n);
1621   value_concat r n (4*n);
1622   assert { value r (4*n)
1623            = lvx0 + m * value r1n (3*n) at Join1
1624            by offset or1n = offset r + n
1625            so value r n = lvx0
1626            so value r1n (3*n) at Join1
1627               = value_sub (pelts r) (offset r + n) (offset r + 4*n) };
1628   assert { m * m * m * m = power radix (4*n)
1629            by m * m = power radix (2*n) };
1630   assert { value r (4*n) + m * m * m * m * !h
1631            = vx0 + m * vy - m * lvx3 - m * m * vx0 + m * m * m * lvx3 };
1632   let rs = s + t - n in
1633   begin ensures { value r (4*n) + m * m * m * m * value r4n rs
1634                   = (value x sx) * (value y sy) }
1635   if [@extraction:likely] (s + t > n)
1636   then begin
1637     let r2n = incr r (2 * n) in
1638     label Sub2 in
1639     value_concat r (2 * n) (4 * n);
1640     assert { value r (4*n) = value r (2*n) + m*m * value r2n (2*n)
1641              by m * m = power radix (2*n)
1642              so value_sub (pelts r) (offset r + 2*n) (offset r + 4*n)
1643                 = value r2n (2*n) };
1644     assert { value r4n rs = hvx3 };
1645     let ghost or = { r } in
1646     let b = wmpn_sub_in_place r2n (2 * n) r4n rs in
1647     value_sub_frame (pelts r) (pelts or) (offset r) (offset r + 2 * int32'int n);
1648     assert { value r (2*n) = value or (2*n) };
1649     value_concat r (2 * n) (4 * n);
1650     assert { value r (4*n) = value r (2*n) + m*m * value r2n (2*n)
1651              by m * m = power radix (2*n)
1652              so pelts r2n = pelts r
1653              so offset r2n = offset r + 2*n
1654              so offset r2n + 2*n = offset r + 4*n
1655              so value_sub (pelts r) (offset r + 2*n) (offset r + 4*n)
1656                 = value r2n (2*n) };
1657     h := Int64.(-) !h (Limb.to_int64 b);
1658     assert { value r2n (2*n) - m * m * b = (value r2n (2*n) at Sub2) - hvx3
1659              by m * m = power radix (2*n) };
1660     assert { value r (4*n) - m * m * m * m * b
1661              = value r (4*n) at Sub2 - m * m * hvx3 };
1662     assert { value r (4*n) + m * m * m * m * !h + m * m * m * m * value r4n rs
1663              = value x sx * value y sy
1664              by value r (4*n) + m * m * m * m * !h
1665              = value r (4*n) + m * m * m * m * (!h at Sub2) - m * m * m * m * b
1666              = value r (4*n) at Sub2 + m * m * m * m * !h at Sub2 - m * m * hvx3
1667              = vx0 + m * vy - m * lvx3 - m * m * vx0 + m * m * m * lvx3
1668                - m * m * hvx3
1669              = vx0 + m * vy - m * vx3 - m * m * vx0 + m * m * m * lvx3
1670              so value r (4*n) + m * m * m * m * !h
1671                 + m * m * m * m * value r4n rs
1672              = vx0 + m * vy - m * vx3 - m * m * vx0 + m * m * m * lvx3
1673                 + m * m * m * m * value r4n rs
1674              = vx0 + m * vy - m * vx3 - m * m * vx0 + m * m * m * vx3
1675              = (a0 + m * a1 + m * m * a2) * (b0 + m * b1)
1676              = value x sx * value y sy };
1677     (if Int64.(<) !h 0
1678     then begin
1679       assert { 0 <= value r4n rs + !h
1680                by 0 <= value x sx so 0 <= value y sy
1681                so 0 <= value x sx * value y sy
1682                so value r (4*n) < power radix (4*n) = m * m * m * m
1683                so 0 <= value r (4*n) + m * m * m * m * (value r4n rs + !h)
1684                     < (m * m * m * m) * (1 + value r4n rs + !h)
1685                so 0 < (m * m * m * m) * (1 + value r4n rs + !h)
1686                so 0 < (m * m * m * m)
1687                so 0 < 1 + value r4n rs + !h };
1688       wmpn_decr r4n rs (Limb.of_int64 (Int64.(-_) !h))
1689       end
1690     else begin
1691       assert { value r4n rs + !h < power radix rs
1692                by value x sx < power radix sx
1693                so value y sy < power radix sy
1694                so value x sx * value y sy
1695                   <= value x sx * power radix sy
1696                   < power radix sx * power radix sy
1697                   = power radix (sx + sy)
1698                   = power radix (3*n+s+t)
1699                so 0 <= value r (4*n)
1700                so power radix (4*n) * (value r4n rs + !h)
1701                   = m * m * m * m * (value r4n rs + !h)
1702                   <= value r (4*n) + m * m * m * m * (value r4n rs + !h)
1703                   < power radix (3*n + s + t) = power radix (4*n + rs)
1704                   = power radix (4*n) * power radix rs
1705                so power radix (4*n) * (value r4n rs + !h)
1706                   < power radix (4*n) * power radix rs
1707                so 0 < power radix (4*n) };
1708       wmpn_incr r4n rs (Limb.of_int64 !h)
1709       end)
1710   end
1711   else begin
1712     assert { !h = 0
1713              by rs = 0
1714              so sx + sy = 4 * n
1715              so hvx3 = 0
1716              so vx3 = lvx3
1717              so value r (4*n) + power radix (4*n) * !h = value x sx * value y sy
1718              so 0 <= value x sx < power radix sx
1719              so 0 <= value y sy < power radix sy
1720              so 0 <= value x sx * value y sy
1721              so value x sx * value y sy
1722                 <= value x sx * power radix sy
1723                 < power radix sx * power radix sy
1724                 = power radix (sx + sy)
1725                 = power radix (4*n)
1726              so 0 <= value r (4*n) < power radix (4*n)
1727              so 0 <= value x sx * value y sy
1728                 < power radix (4*n) + power radix (4*n) * !h
1729                 = power radix (4*n) * (1 + !h)
1730              so 0 < power radix (4*n) * (1 + !h)
1731              so 0 < 1 + !h
1732              so power radix (4*n) * 1 > value x sx * value y sy
1733                 >= power radix (4*n) * !h
1734              so power radix (4*n) * !h < power radix (4*n) * 1
1735              so !h < 1 };
1736   end
1737   end;
1738   label Join4 in
1739   let ghost or = { r } in
1740   let ghost or4n = { r4n } in
1741   join r r4n;
1742   value_sub_frame (pelts r) (pelts or) (offset r) (offset r + 4 * int32'int n);
1743   value_sub_frame (pelts r) (pelts or4n) (offset r + 4 * int32'int n)
1744                   (offset r + 3 * int32'int n + int32'int s + int32'int t);
1745   value_concat r (4*n) (3*n + s + t);
1746   assert { value r (3*n + s + t) = value x sx * value y sy
1747            by offset or4n = offset r + 4*n
1748            so offset r4n + rs = offset r + 3*n + s + t
1749            so value_sub (pelts r) (offset r + 4*n) (offset r + 3*n + s + t)
1750               = value_sub (pelts or4n) (offset r + 4*n) (offset r + 3*n + s + t)
1751               = value or4n rs
1752            so value r (3*n + s + t)
1753               = value or (4*n) + m * m * m * m * value or4n rs };
1754   join scratch vy2;
1755   join scratch sor;
1756   join_r sol scratch;
1757   label JoinR in
1758   let ghost or = { r } in
1759   join r ror;
1760   value_sub_frame (pelts r) (pelts or) (offset r)
1761                   (offset r + int32'int sx + int32'int sy);
1762   label JoinL in
1763   join_r rol r;
1764   value_sub_frame (pelts r) (pelts or) (offset r)
1765                   (offset r + int32'int sx + int32'int sy);
1766   assert { value r (sx + sy) = value x sx * value y sy
1767            by value r (sx + sy) = value r (sx + sy) at JoinL
1768               = value r (sx + sy) at JoinR }
1769
1770   with wmpn_mul (r x:ptr limb) (sx:int32) (y:ptr limb) (sy:int32) (ghost k: int) : limb
1771     requires { valid x sx }
1772     requires { valid y sy }
1773     requires { valid r (sx + sy) }
1774     requires { writable r }
1775     requires { 0 < sy <= sx }
1776     requires { 8 * sx < max_int32 }
1777     requires { sx <= toom22_threshold * power 2 k }
1778     requires { 0 <= k <= 64 }
1779     ensures  { min r = old min r }
1780     ensures  { max r = old max r }
1781     ensures  { plength r = old plength r }
1782     ensures  { value r (sx + sy) = value x sx * value y sy }
1783     ensures  { result = (pelts r)[offset r + sx + sy - 1] }
1784     ensures  { forall j. min r <= j < offset r \/ offset r + sx + sy <= j < max r
1785                          -> (pelts r)[j] = old (pelts r)[j] }
1786     variant  { k+k+1 }
1787   =
1788     if sy <= toom22_threshold
1789     then
1790     (* TODO block product if sx large, for memory locality according to GMP *)
1791      wmpn_mul_basecase r x sx y sy
1792     else begin
1793       (* this would be faster with salloc *)
1794       (*let ghost k = 64 in*) (* is always enough *)
1795       let scratch = salloc (UInt32.of_int32 ((5 * sy) + 128)) in
1796       (* c_assert (is_not_null scratch); *)
1797       let rol = decr_split r 0 in
1798       let ror = incr_split r (sx + sy) in
1799       if ((2 * sx) >= (5 * sy))
1800       then begin
1801         let un = ref sx in
1802         let su = (3 * sy) / 2 in
1803         assert { 0 < su };
1804         let ghost sr = su + sy in
1805         let ws = salloc (UInt32.of_int32 (4 * sy)) in
1806         (* c_assert (is_not_null ws); *)
1807         wmpn_toom32_mul r x su y sy scratch k;
1808         un := !un - su;
1809         let up = ref (C.incr x su) in
1810         let rp = ref (C.incr r su) in
1811         let ghost ou = ref su in
1812         let ghost or = ref sr in
1813         while (!un >= (2 * sy)) do (* 5/2?*)
1814           invariant { min_int32 <= 2 * !un <= max_int32 }
1815           invariant { !ou + !un = sx }
1816           invariant { !or = !ou + sy }
1817           invariant { su <= !ou < sx }
1818           invariant { !un < sx }
1819           invariant { 2 * sy - su <= !un }
1820           invariant { value r !or = value x !ou * value y sy }
1821           invariant { offset !rp = offset r + !ou }
1822           invariant { offset !up = offset x + !ou }
1823           invariant { min !up = min x /\ max !up = max x }
1824           invariant { plength !up = plength x }
1825           invariant { min !rp = min r /\ max !rp = max r }
1826           invariant { plength !rp = plength r }
1827           invariant { writable !rp }
1828           invariant { min ws = 0 /\ max ws = plength ws = 4 * sy }
1829           invariant { min scratch = 0 /\ max scratch = plength scratch }
1830           invariant { plength scratch = 5 * sy + 128 }
1831           invariant { pelts !rp = pelts r }
1832           invariant { pelts !up = pelts x }
1833           variant { p2i !un }
1834           (*wmpn_copyi ws !rp sy;
1835           let rr = rp.contents in
1836           wmpn_toom32_mul rr !up y scratch su sy k;
1837           let cy = wmpn_add_in_place rr ws sy sy in
1838           let rpn = C.incr rr sy in
1839           wmpn_incr rpn cy su;
1840           un := !un - su;
1841           up.contents <- C.incr !up su;
1842           ou := !ou + su;
1843           or := !or + su;
1844           rp.contents <- C.incr !rp su;*)
1845           label StartLoop in
1846           let ghost o_r = { r } in
1847           let ghost rrp = !rp in
1848           let ghost o_rp = { rrp } in (* TODO why not { !rp } ? *)
1849           value_concat r !ou !or;
1850           assert { value r !or = value r !ou + power radix !ou * value !rp sy };
1851           wmpn_toom32_mul ws !up su y sy scratch k;
1852           let cy = wmpn_add_n_in_place !rp ws sy in
1853           value_sub_frame (pelts r) (pelts o_r) (offset r) (offset r + p2i !ou);
1854           assert { value r !ou = value o_r !ou };
1855           assert { value !rp sy + (power radix sy) * cy
1856                    = value o_rp sy + value ws sy };
1857           let rpn = C.incr !rp sy in
1858           let wsy = C.incr ws sy in
1859           let orp = { rpn } in
1860           label Copy in
1861           wmpn_copyi rpn wsy su;
1862           value_sub_frame_shift (pelts rpn) (pelts wsy)
1863                                 (offset rpn) (offset wsy) (int32'int su);
1864           value_sub_frame (pelts r) (pelts orp) (offset r) (offset r + p2i !ou);
1865           value_sub_frame (pelts rpn) (pelts orp)
1866                           (offset !rp) (offset !rp + p2i sy);
1867           assert { value rpn su = value wsy su };
1868           assert { value !rp sy = value (!rp at Copy) sy };
1869           assert { value r !ou = value o_r !ou };
1870           value_concat ws sy sr;
1871           assert { value ws sr
1872                    = value ws (sy + su)
1873                    = value ws sy + power radix sy * value wsy su };
1874           value_concat r !ou (!ou + sr);
1875           assert { value r (!or + su) = value r (!ou + sr)
1876                    = value r !ou + power radix !ou * value !rp sr };
1877           value_concat !rp sy sr;
1878           assert { value !rp sr = value !rp sy + power radix sy * value rpn su };
1879           value_concat x !ou (!ou + su);
1880           assert { value x (!ou + su)
1881                    = value x !ou + power radix !ou * value !up su };
1882           assert { value r (!ou + sr) + (power radix !or) * cy
1883                    = value x (!ou + su) * value y sy
1884                    by value r (!ou + sr) + (power radix !or) * cy
1885                       = value r !ou + power radix !ou * value !rp sr
1886                         + power radix !ou * (power radix sy) * cy
1887                       = value r !ou + power radix !ou * value !rp sy
1888                         + power radix !ou * power radix sy * value wsy su
1889                         + power radix !ou * power radix sy * cy
1890                       = value r !ou
1891                         + power radix !ou * (value !rp sy + power radix sy * cy)
1892                         + power radix !ou * power radix sy * value wsy su
1893                       = value r !ou
1894                         + power radix !ou * (value o_rp sy + value ws sy)
1895                         + power radix !ou * power radix sy * value wsy su
1896                       = value r !ou + power radix !ou * value o_rp sy
1897                         + power radix !ou
1898                           * (value ws sy + power radix sy * value wsy su)
1899                       = value o_r !or
1900                         + power radix !ou * value ws sr
1901                       = value x !ou * value y sy
1902                         + power radix !ou * value !up su * value y sy
1903                       = (value x !ou + power radix !ou * value !up su)
1904                         * value y sy
1905                       = value x (!ou + su) * value y sy };
1906           value_concat r !or (!ou + sr);
1907           assert { value r (!ou + sr)
1908                       = value r !or + power radix !or * value rpn su
1909                    by value rpn su = value_sub (pelts r) (offset r + !or) (offset r + !ou + sr) };
1910           assert { value rpn su + cy < power radix su
1911                    by value x (!ou + su) < power radix (!ou + su)
1912                    so value y sy < power radix sy
1913                    so value x (!ou + su) * value y sy
1914                       < power radix (!ou + su) * power radix sy
1915                       = power radix (!ou + su + sy)
1916                       = power radix (!or + su)
1917                    so value r (!ou + sr) + power radix !or * cy
1918                       = value r !or + power radix !or * (value rpn su)
1919                         + power radix !or * cy
1920                       = value r !or + power radix !or * (value rpn su + cy)
1921                    so value r !or + power radix !or * (value rpn su + cy)
1922                       < power radix (!or + su)
1923                    so 0 <= value r !or
1924                    so power radix !or * (value rpn su + cy)
1925                       < power radix (!or + su)
1926                       = power radix !or * power radix su };
1927           label Incr in
1928           let orp = { rpn } in
1929           wmpn_incr rpn su cy;
1930           value_sub_frame (pelts r) (pelts orp) (offset r) (offset r + p2i !or);
1931           assert { value r !or = value r !or at Incr };
1932           value_concat r !or (!ou + sr);
1933           assert { value r (!ou + sr) = value x (!ou + su) * value y sy
1934                    by value r (!ou + sr)
1935                    = value r !or + power radix !or * value rpn su
1936                    = value r !or at Incr + power radix !or * value rpn su
1937                    = value r !or at Incr
1938                      + power radix !or * (value rpn su at Incr + cy)
1939                    = value r !or at Incr
1940                      + power radix !or * value rpn su at Incr
1941                      + power radix !or * cy
1942                    = value r (!ou + sr) at Incr + power radix !or * cy };
1943           un := !un - su;
1944           up.contents <- C.incr !up su;
1945           ou := !ou + su;
1946           or := !or + su;
1947           rp.contents <- C.incr !rp su;
1948         done;
1949  (*       wmpn_copyi ws !rp sy;*)
1950         (* sy <= !un <= 2.5 * sy *)
1951         value_concat r !ou !or;
1952         assert { value r !or = value r !ou + power radix !ou * value !rp sy };
1953         let ghost o_r = { r } in
1954         let ghost rrp = !rp in
1955         let ghost o_rp = { rrp } in (* TODO why not { !rp } ? *)
1956         begin
1957           ensures { value ws (!un + sy) = value !up !un * value y sy }
1958           ensures { min ws = old min ws /\ max ws = old max ws
1959                     /\ plength ws = old plength ws }
1960           ensures { min scratch = old min scratch
1961                     /\ max scratch = old max scratch
1962                     /\ plength scratch = old plength scratch }
1963           if sy <= !un
1964           then begin
1965             if ((4 * !un) < (5 * sy))
1966             then wmpn_toom22_mul ws !up !un y sy scratch k
1967             else wmpn_toom32_mul ws !up !un y sy scratch k
1968           end
1969           else let _ = wmpn_mul ws y sy !up !un (k-1) in ()
1970         end;
1971         let cy = wmpn_add_n_in_place !rp ws sy in
1972         value_sub_frame (pelts r) (pelts o_r) (offset r) (offset r + p2i !ou);
1973         assert { value r !ou = value o_r !ou };
1974         assert { value !rp sy + (power radix sy) * cy
1975                  = value o_rp sy + value ws sy };
1976         assert { valid !rp sy };
1977         let rpn = C.incr !rp sy in
1978         let wsy = C.incr ws sy in
1979         let orp = { rpn } in
1980         label Copy in
1981         wmpn_copyi rpn wsy !un;
1982         value_sub_frame_shift (pelts rpn) (pelts wsy)
1983                               (offset rpn) (offset wsy) (int32'int !un);
1984         value_sub_frame (pelts r) (pelts orp) (offset r) (offset r + p2i !ou);
1985         value_sub_frame (pelts rpn) (pelts orp)
1986                         (offset !rp) (offset !rp + p2i sy);
1987         assert { value rpn !un = value wsy !un };
1988         assert { value !rp sy = value (!rp at Copy) sy };
1989         assert { value r !ou = value o_r !ou };
1990         let ghost sr = sy + !un in
1991         value_concat ws sy sr;
1992         assert { value ws sr
1993                  = value ws sy + power radix sy * value wsy !un };
1994         value_concat r !ou (!ou + sr);
1995         assert { value r (!or + !un) = value r (!ou + sr)
1996                  = value r !ou + power radix !ou * value !rp sr };
1997         value_concat !rp sy sr;
1998         assert { value !rp sr = value !rp sy + power radix sy * value rpn !un };
1999         value_concat x !ou (!ou + !un);
2000         assert { value x (!ou + !un)
2001                  = value x !ou + power radix !ou * value !up !un };
2002         assert { value r (!ou + sr) + (power radix !or) * cy
2003                  = value x (!ou + !un) * value y sy
2004                  by value r (!ou + sr) + (power radix !or) * cy
2005                     = value r !ou + power radix !ou * value !rp sr
2006                       + power radix !ou * (power radix sy) * cy
2007                     = value r !ou + power radix !ou * value !rp sy
2008                       + power radix !ou * power radix sy * value wsy !un
2009                       + power radix !ou * power radix sy * cy
2010                     = value r !ou
2011                       + power radix !ou * (value !rp sy + power radix sy * cy)
2012                       + power radix !ou * power radix sy * value wsy !un
2013                     = value r !ou
2014                       + power radix !ou * (value o_rp sy + value ws sy)
2015                       + power radix !ou * power radix sy * value wsy !un
2016                     = value r !ou + power radix !ou * value o_rp sy
2017                       + power radix !ou
2018                         * (value ws sy + power radix sy * value wsy !un)
2019                     = value o_r !or
2020                       + power radix !ou * value ws sr
2021                     = value x !ou * value y sy
2022                       + power radix !ou * value !up !un * value y sy
2023                     = (value x !ou + power radix !ou * value !up !un)
2024                       * value y sy
2025                     = value x (!ou + !un) * value y sy };
2026         value_concat r !or (!ou + sr);
2027         assert { value r (!ou + sr)
2028                     = value r !or + power radix !or * value rpn !un };
2029         assert { value rpn !un + cy < power radix !un
2030                  by value x (!ou + !un) < power radix (!ou + !un)
2031                  so value y sy < power radix sy
2032                  so value x (!ou + !un) * value y sy
2033                     < power radix (!ou + !un) * power radix sy
2034                     = power radix (!ou + !un + sy)
2035                     = power radix (!or + !un)
2036                  so value r (!ou + sr) + power radix !or * cy
2037                     = value r !or + power radix !or * (value rpn !un)
2038                       + power radix !or * cy
2039                     = value r !or + power radix !or * (value rpn !un + cy)
2040                  so value r !or + power radix !or * (value rpn !un + cy)
2041                     < power radix (!or + !un)
2042                  so 0 <= value r !or
2043                  so power radix !or * (value rpn !un + cy)
2044                     < power radix (!or + !un)
2045                     = power radix !or * power radix !un };
2046         let orp = { rpn } in
2047         label Incr in
2048         wmpn_incr rpn !un cy;
2049         value_sub_frame (pelts r) (pelts orp) (offset r) (offset r + p2i !or);
2050         assert { value r !or = value r !or at Incr };
2051         value_concat r !or (!ou + sr);
2052         assert { value r (!ou + sr) = value x (!ou + !un) * value y sy
2053                  by value r (!ou + sr)
2054                  = value r !or + power radix !or * value rpn !un
2055                  = value r !or at Incr + power radix !or * value rpn !un
2056                  = value r !or at Incr
2057                    + power radix !or * (value rpn !un at Incr + cy)
2058                  = value r !or at Incr
2059                    + power radix !or * value rpn !un at Incr
2060                    + power radix !or * cy
2061                  = value r (!ou + sr) at Incr + power radix !or * cy };
2062         assert { value r (sx + sy) = value x sx * value y sy
2063                  by !ou + sr = sx + sy /\ !ou + !un = sx };
2064         sfree ws;
2065       end
2066       else begin
2067         if ((4 * sx) < (5 * sy))
2068         then wmpn_toom22_mul r x sx y sy scratch k
2069         else wmpn_toom32_mul r x sx y sy scratch k
2070       end;
2071       sfree scratch;
2072       label JoinR in
2073       let ghost or = { r } in
2074       join r ror;
2075       value_sub_frame (pelts r) (pelts or) (offset r)
2076                       (offset r + int32'int sx + int32'int sy);
2077       label JoinL in
2078       join_r rol r;
2079       value_sub_frame (pelts r) (pelts or) (offset r)
2080                       (offset r + int32'int sx + int32'int sy);
2081       assert { value r (sx + sy) = value x sx * value y sy
2082                by value r (sx + sy) = value r (sx + sy) at JoinL
2083                 = value r (sx + sy) at JoinR }
2084     end;
2085   C.get_ofs r (sx + sy - 1)
2086         (* sy <= !un <= 2.5 * sy *)
2087        (* if sy <= !un
2088         then begin
2089            if ((4 * !un) < (5 * sy))
2090            then wmpn_toom22_mul ws !up y scratch !un sy k
2091            else wmpn_toom32_mul ws !up y scratch !un sy k
2092         end
2093         else wmpn_mul ws y !up sy !un (k-1);
2094         let cy = wmpn_add_in_place !rp ws sy sy in
2095         let rpn = C.incr !rp sy in
2096         wmpn_copyi rpn (C.incr ws sy) !un;
2097         wmpn_incr rpn cy !un;
2098         sfree ws;
2099       end
2100       else begin
2101         if ((4 * sx) < (5 * sy))
2102         then wmpn_toom22_mul r x y scratch sx sy k
2103         else wmpn_toom32_mul r x y scratch sx sy k
2104       end;
2105       sfree scratch
2106     end*)
2107
2108   (** `wmpn_mul_n r x y sz` multiplies `(x, sz)` and `(y, sz)` and
2109   writes the result to `(r, 2*sz)`. `r` must not overlap with either
2110   `x` or `y`. Corresponds to `mpn_mul_n`.  *)
2111   let wmpn_mul_n (r x y:t) (sz:int32) (ghost k: int) : unit
2112     requires { sz > 0 }
2113     requires { valid x sz }
2114     requires { valid y sz }
2115     requires { valid r (sz + sz) }
2116     requires { writable r }
2117     requires { 8 * sz < max_int32 }
2118     requires { sz <= toom22_threshold * power 2 k }
2119     requires { 0 <= k <= 64 }
2120     ensures { value r (sz + sz) = value x sz * value y sz }
2121     ensures  { min r = old min r }
2122     ensures  { max r = old max r }
2123     ensures  { plength r = old plength r }
2124     ensures  { forall j. min r <= j < offset r \/ offset r + sz + sz <= j < max r
2125                          -> (pelts r)[j] = old (pelts r)[j] }
2126   =
2127     if sz <= toom22_threshold
2128     then wmpn_mul_basecase r x sz y sz
2129     else
2130       let ws = salloc (UInt32.of_int32 (2 * (sz + 64))) in
2131       wmpn_toom22_mul r x sz y sz ws 64
2132
2133   meta remove_prop axiom no_borrow
2134   meta remove_prop axiom no_borrow_ptr
2135   meta remove_prop lemma power_ge_1
2136   meta remove_prop axiom valuation'def
2137   meta remove_prop axiom valuation'spec
2138   meta remove_prop lemma valuation_mul
2139   meta remove_prop axiom valuation_times_pow
2140   meta remove_prop axiom valuation_lower_bound
2141   meta remove_prop axiom valuation_monotonous
2142   meta remove_prop lemma valuation_nondiv
2143   meta remove_prop lemma valuation_zero_prod
2144   meta remove_prop axiom valuation_times_nondiv
2145   meta remove_prop lemma valuation_split
2146   meta remove_prop lemma valuation_prod
2147   meta remove_prop lemma valuation_mod
2148   meta remove_prop lemma valuation_decomp
2149   meta remove_prop lemma valuation_pow
2150
2151 end