module Mul

  use int.Int
  use mach.int.Int32
  use import mach.int.UInt64GMP as Limb
  use int.Power
  use ref.Ref
  use mach.c.C
  use array.Array
  use map.Map
  use types.Types
  use types.Int32Eq
  use types.UInt64Eq
  use lemmas.Lemmas
  use util.Util
  use add.Add

  let wmpn_mul_1 (r x:t) (sz:int32) (y:limb) : limb
    requires { valid x sz }
    requires { valid r sz }
    requires { writable r }
    ensures { value r sz + (power radix sz) * result
                = value x sz * y }
    ensures { forall j. (j < offset r \/ offset r + sz <= j)
              -> (pelts r)[j] = old (pelts r)[j] }
    writes { r.data.elts }
  =
    let ref cl = 0 in
    let ref ul = 0 in
    let ref n = sz in
    let ref up = C.incr x 0 in
    let ref rp = C.incr r 0 in
    let ghost ref i : int32 = 0 in
    while n <> 0 do
      invariant { 0 <= n <= sz }
      invariant { i = sz - n }
      invariant { value r i + (power radix i) * cl =
                  value x i * y }
      invariant { rp.offset = r.offset + i }
      invariant { plength rp = plength r }
      invariant { rp.min = r.min }
      invariant { rp.max = r.max }
      invariant { writable rp }
      invariant { pelts rp = pelts r }
      invariant { up.offset = x.offset + i }
      invariant { plength up = plength x }
      invariant { up.min = x.min }
      invariant { up.max = x.max }
      invariant { pelts up = pelts x }
      invariant { forall j. (j < offset r \/ offset r + sz <= j)
                -> (pelts r)[j] = old (pelts r)[j] }
      variant { n }
      label StartLoop in
      ul <- C.get up;
      up <- C.incr up 1;
      let l, h = mul_double ul y in
      assert { h < radix - 1
               by ul * y <= (radix - 1) * (radix - 1)
               so radix * h <= ul * y };
      let lpl = add_mod l cl in
      begin ensures { lpl + radix * cl = ul * y + (cl at StartLoop) }
        cl <- (if lpl < cl then 1 else 0) + h;
      end;
      value_sub_update_no_change (pelts r) (r.offset + int32'int i)
                                 r.offset (r.offset + int32'int i) lpl;
      C.set rp lpl;
      assert { value r i = value r i at StartLoop };
      assert { (pelts r)[offset r + i] = lpl };
      value_tail r i;
      value_tail x i;
      assert { value x (i+1) = value x i + power radix i * ul };
      assert { value x (i+1) * y = value x i * y + power radix i * (ul * y) };
      assert { value r (i+1) + power radix (i+1) * cl = value x (i+1) * y };
      rp <- C.incr rp 1;
      n <- n-1;
      i <- i+1;
    done;
    cl

  let wmpn_addmul_1 (r x:t) (sz: int32) (y:limb) : limb
    requires { valid x sz }
    requires { valid r sz }
    requires { writable r }
    ensures { value r sz + (power radix sz) * result
            = value (old r) sz
              + value x sz * y }
    writes { r.data.elts }
    ensures { forall j. (j < r.offset \/ r.offset + sz <= j) ->
              (pelts r)[j] = (pelts (old r))[j] }
  =
    let ref ul = 0 in
    let ref rl = 0 in
    let ref cl = 0 in
    let ref n = sz in
    let ghost ref i = 0 in
    let ref rp = C.incr r 0 in
    let ref up = C.incr x 0 in
    while n <> 0 do
      invariant { 0 <= n <= sz }
      invariant { i = sz - n }
      invariant { value r i + (power radix i) * cl
                  = value (old r) i
                    + value x i * y }
      invariant { (rp).offset = r.offset + i }
      invariant { rp.min = r.min }
      invariant { rp.max = r.max }
      invariant { writable rp }
      invariant { pelts rp = pelts r }
      invariant { up.offset = x.offset + i }
      invariant { plength up = plength x }
      invariant { up.min = x.min }
      invariant { up.max = x.max }
      invariant { pelts up = pelts x }
      invariant { forall j. r.offset + i <= j < r.offset + sz ->
                 (pelts (old r)) [j] = (pelts r)[j]  }
      invariant { forall j. j < r.offset \/ r.offset + sz <= j ->
                 (pelts r)[j] = (pelts (old r))[j] }
      variant { sz - i }
      label StartLoop in
      ul <- C.get up;
      up <- C.incr up 1;
      let l, h = mul_double ul y in
      let ref lpl = add_mod l cl in
      assert { h < radix - 1
               by ul * y <= (radix - 1) * (radix - 1)
               so radix * h <= ul * y };
      begin ensures { lpl + radix * cl = ul * y + (cl at StartLoop) }
        cl <- (if lpl < cl then 1 else 0) + h;
      end;
      assert { cl = radix - 1 -> lpl = 0
               by ul * y <= (radix - 1) * (radix - 1)
               so cl at StartLoop <= radix - 1
               so lpl + radix * cl
                  = ul * y + (cl at StartLoop) = radix * (radix - 1) };
      rl <- C.get rp;
      assert { rl = (pelts (old r))[r.offset + i] };
      lpl <- add_mod rl lpl;
      begin ensures { lpl + radix * cl = ul * y + (cl at StartLoop) + rl }
        cl <- (if lpl < rl then 1 else 0) + cl;
      end;
      value_sub_update_no_change (pelts r) (r.offset + int32'int i)
                                 r.offset (r.offset + int32'int i) lpl;
      C.set rp lpl;
      assert { value r i = value r i at StartLoop };
      value_tail r i;
      value_tail x i;
      assert { value (old r) (i+1)
               = value (old r) i + power radix i * rl };
      assert { value x (i+1) = value x i + power radix i * ul };
      assert { value r (i+1) = value r i + power radix i * lpl };
      assert { value x (i+1) * y = value x i * y + power radix i * (ul * y) };
      assert { value r (i+1) + power radix (i+1) * cl =
               value (old r) (i+1) + value x (i+1) * y };
      rp <- C.incr rp 1;
      n <- n-1;
      i <- i+1;
    done;
    cl

 let wmpn_addmul_n (r x y:t) (sz:int32) : limb
    requires { sz > 0 }
    requires { valid x sz }
    requires { valid y sz }
    requires { valid r (sz + sz) }
    requires { writable r }
    writes { r.data.elts }
    ensures { value r (sz + sz)
                + power radix (sz + sz) * result
              = value (old r) (sz + sz)
                + value x sz * value y sz }
  = [@vc:do_not_keep_trace] (* traceability info breaks the proof *)
    let ref rp = C.incr r 0 in
    let ref vp = C.incr y 0 in
    let ref lr = 0 in
    let ref c = 0 in
    let ref vn = sz in
    let ghost ref i = 0 in
    while vn <> 0 do
      invariant { 0 <= i <= sz }
      invariant { i = sz - vn }
      invariant { value r (i + sz)
                    + (power radix (i + sz)) * c
                  = value (old r) (i + sz)
                    + value x sz * value y i }
      invariant { rp.offset = r.offset + i }
      invariant { rp.min = r.min }
      invariant { rp.max = r.max }
      invariant { writable rp }
      invariant { pelts rp = pelts r }
      invariant { plength rp = plength r }
      invariant { vp.offset = y.offset + i }
      invariant { plength vp = plength y }
      invariant { vp.min = y.min }
      invariant { vp.max = y.max }
      invariant { pelts vp = pelts y }
      invariant { 0 <= c <= 1 }
      invariant { forall j. rp.offset + sz <= j ->
                 (pelts (old r)) [j] = (pelts r)[j] }
      variant { sz - i }
      label StartLoop in
      value_concat r i (i+sz);
      assert { value rp sz
               = value_sub (pelts r) (r.offset + i) (r.offset + (i + sz)) };
      let ghost ly = pure { (pelts y)[y.offset + i] } in
      let c' =  wmpn_addmul_1 rp x sz (C.get vp) in
      assert { forall j. rp.offset + sz <= j ->
                 (pelts (old r)) [j] = (pelts r)[j]
                 by (pelts r)[j]
                 = (pelts rp)[j]
                 = (pelts rp)[j] at StartLoop
                 = (pelts (old r))[j]};
      assert { value rp sz + power radix sz * c'
              = value (rp at StartLoop) sz + value x sz * ly };
      assert { MapEq.map_eq_sub (pelts r) (pelts r at StartLoop)
                r.offset rp.offset
                by rp.offset = r.offset + i
                so forall j. r.offset <= j < rp.offset
                   ->
                   (j < rp.offset
                    so (pelts rp)[j] = (pelts rp at StartLoop)[j]
                         = (pelts r at StartLoop)[j]) };
      lr <- get_ofs rp sz;
      assert { lr = (pelts (old r))[(old r).offset + (i + sz)] };
      let (res, carry) = add_with_carry c' lr c in
      label BeforeCarry in
      value_sub_update_no_change (pelts r) (rp.offset + p2i sz)
                        r.offset  (r.offset + p2i i) res;
      set_ofs rp sz res;
      assert { value rp sz = value (rp at BeforeCarry) sz };
      c <- carry;
      assert { value r i = value (r at BeforeCarry) i
               = value (r at StartLoop) i};
      value_tail r (i+sz);
      value_tail y i;
      value_tail (pure { old r }) (i+sz);
      assert { value (old r) ((i+sz)+1)
               = value (old r) (i+sz) + power radix (i+sz) * lr };
      assert { (pelts r)[r.offset + (i + sz)] = res };
      value_concat r i (i+sz);
      assert { value_sub (pelts r) (r.offset + i) (r.offset+(i+sz))
               = value rp sz };
      assert { value r (i + sz + 1) =
               value r i + power radix i * (value rp sz) + power radix (i + sz) * res };
      assert { value x sz * value y (i+1)
               = value x sz * value y i + power radix i * (value x sz * ly) };
             (* nonlinear *)
      assert { value r (i + sz + 1)
                    + (power radix (i + sz + 1)) * c
                  = value (old r) (i + sz + 1)
                    + value x sz * value y (i + 1) };
      i <- i + 1;
      rp <- C.incr rp 1;
      vp <- C.incr vp 1;
      vn <- vn - 1;
      assert { forall j. rp.offset + sz <= j ->
                 (pelts (old r)) [j] = (pelts r)[j] };
    done;
    c

  let wmpn_mul_1_in_place (x:t) (sz:int32) (y:limb) : limb
    requires { valid x sz }
    requires { writable x }
    ensures  { value x sz + (power radix sz) * result
               = old value x sz * y }
    ensures  { forall j. (j < offset x \/ offset x + sz <= j)
              -> (pelts x)[j] = old (pelts x)[j] }
    writes   { x.data.elts }
  =
    let ghost ox = pure { x } in
    let ref cl = 0 in
    let ref ul = 0 in
    let ref n = sz in
    let ref up = C.incr x 0 in
    let ghost ref i : int32 = 0 in
    while n <> 0 do
      invariant { 0 <= n <= sz }
      invariant { i = sz - n }
      invariant { value x i + (power radix i) * cl =
                  value ox i * y }
      invariant { up.offset = x.offset + i }
      invariant { plength up = plength x }
      invariant { up.min = x.min }
      invariant { up.max = x.max }
      invariant { pelts up = pelts x }
      invariant { writable up }
      invariant { forall j. (j < offset x \/ offset x + i <= j)
                -> (pelts x)[j] = (pelts ox)[j] }
      variant { n }
      label StartLoop in
      ul <- C.get up;
      let l, h = mul_double ul y in
      assert { h < radix - 1
               by ul * y <= (radix - 1) * (radix - 1)
               so radix * h <= ul * y };
      let lpl = add_mod l cl in
      begin ensures { lpl + radix * cl = ul * y + (cl at StartLoop) }
        cl <- (if lpl < cl then 1 else 0) + h;
      end;
      value_sub_update_no_change (pelts x) (x.offset + int32'int i)
                                 x.offset (x.offset + int32'int i) lpl;
      C.set up lpl;
      assert { value x i = value x i at StartLoop };
      assert { (pelts x)[offset x + i] = lpl };
      value_tail x i;
      value_tail ox i;
      assert { value ox (i+1) = value ox i + power radix i * ul };
      assert { value ox (i+1) * y = value ox i * y + power radix i * (ul * y) };
      assert { value x (i+1) + power radix (i+1) * cl = value ox (i+1) * y };
      up <- C.incr up 1;
      n <- n-1;
      i <- i+1;
    done;
    cl

  let wmpn_submul_1 (r x:t) (sz:int32) (y:limb):limb
    requires { valid x sz }
    requires { valid r sz }
    requires { writable r }
    ensures { value r sz - (power radix sz) * result
            = value (old r) sz
              - value x sz * y }
    writes { r.data.elts }
    ensures { forall j. j < r.offset \/ r.offset + sz <= j ->
              (pelts r)[j] = (pelts (old r))[j] }
  =
    let ref ul = 0 in
    let ref rl = 0 in
    let ref cl = 0 in
    let ref n = sz in
    let ghost ref i = 0 in
    let ref rp = C.incr r 0 in
    let ref up = C.incr x 0 in
    while n <> 0 do
      invariant { 0 <= n <= sz }
      invariant { i = sz - n }
      invariant { value r i - (power radix i) * cl
                  = value (old r) i
                    - value x i * y }
      invariant { (rp).offset = r.offset + i }
      invariant { rp.min = r.min }
      invariant { rp.max = r.max }
      invariant { writable rp }
      invariant { pelts rp = pelts r }
      invariant { up.offset = x.offset + i }
      invariant { plength up = plength x }
      invariant { up.min = x.min }
      invariant { up.max = x.max }
      invariant { pelts up = pelts x }
      invariant { forall j. r.offset + i <= j < r.offset + sz ->
                 (pelts (old r)) [j] = (pelts r)[j]  }
      invariant { forall j. j < r.offset \/ r.offset + sz <= j ->
                 (pelts r)[j] = (pelts (old r))[j] }
      variant { sz - i }
      label StartLoop in
      ul <- C.get up;
      up <- C.incr up 1;
      let l, h = mul_double ul y in
      let ref lpl = add_mod l cl in
      assert { h < radix - 1
               by ul * y <= (radix - 1) * (radix - 1)
               so radix * h <= ul * y };
      begin ensures { lpl + radix * cl = ul * y + (cl at StartLoop) }
        cl <- (if lpl < cl then 1 else 0) + h;
      end;
      assert { cl = radix - 1 -> lpl = 0
               by ul * y <= (radix - 1) * (radix - 1)
               so cl at StartLoop <= radix - 1
               so lpl + radix * cl
                  = ul * y + (cl at StartLoop) = radix * (radix - 1) };
      rl <- C.get rp;
      assert { rl = (pelts (old r))[r.offset + i] };
      lpl <- sub_mod rl lpl;
      begin ensures { lpl - radix * cl = rl - ul * y - (cl at StartLoop) }
        cl <- (if lpl > rl then 1 else 0) + cl;
      end;
      value_sub_update_no_change (pelts r) (r.offset + int32'int i)
                                 r.offset (r.offset + int32'int i) lpl;
      C.set rp lpl;
      assert { value r i = value r i at StartLoop };
      value_tail r i;
      value_tail x i;
      assert { value (old r) (i+1)
               = value (old r) i + power radix i * rl };
      assert { value x (i+1) = value x i + power radix i * ul };
      assert { value r (i+1) = value r i + power radix i * lpl };
      assert { value x (i+1) * y = value x i * y + power radix i * (ul * y) };
      assert { value r (i+1) - power radix (i+1) * cl =
               value (old r) (i+1) - value x (i+1) * y };
      rp <- C.incr rp 1;
      n <- n-1;
      i <- i+1;
    done;
    cl

 let wmpn_addmul_2 (r x:t) (sz:int32) (y:t) : limb
    requires { sz > 0 }
    requires { valid x sz }
    requires { valid y 2 }
    requires { valid r (sz + 2) }
    requires { writable r }
    writes   { r.data.elts }
    ensures  { value r (sz + 1)
                 + power radix (sz + 1) * result
               = value (old r) sz
                 + value x sz * value y 2 }
    ensures  { forall j. (j < r.offset \/ r.offset + sz + 1 <= j) ->
               (pelts r)[j] = (pelts (old r))[j] }
  =
    let ghost or = pure { r } in
    let y0 = C.get y in
    let rn = wmpn_addmul_1 r x sz y0 in
    value_sub_update_no_change (pelts r) (offset r + int32'int sz) (offset r)
                               (offset r + int32'int sz) rn;
    C.set_ofs r sz rn;
    value_tail r sz;
    assert { value r (sz + 1) = value or sz + value x sz * y0 };
    value_concat r 1 (sz + 1);
    let r1 = C.incr r 1 in
    assert { value r (sz + 1) = value r 1 + radix * value r1 sz };
    let y1 = C.get_ofs y 1 in
    assert { value y 2 = y0 + radix * y1 };
    label A2 in
    let c = wmpn_addmul_1 r1 x sz y1 in
    value_concat r 1 (sz + 1);
    value_sub_frame (pelts r) (pure { pelts r at A2 })
                    (offset r) (offset r + 1);
    assert { value r1 sz + power radix sz * c
             = value r1 sz at A2 + value x sz * y1 };
    assert { value r (sz + 1) + power radix (sz + 1) * c
             = value or sz + value x sz * value y 2
             by value r (sz + 1) + power radix (sz + 1) * c
                = value r 1 + radix * value r1 sz + power radix (sz + 1) * c
                = value r 1 + radix * (value r1 sz + power radix sz * c)
                = value r 1 + radix * (value r1 sz at A2 + value x sz * y1)
                = value r 1 at A2 + radix * value r1 sz at A2
                  + radix * (value x sz * y1)
                = value r (sz + 1) at A2 + radix * (value x sz * y1)
                = value or sz + value x sz * y0 + radix * value x sz * y1
                = value or sz + value x sz * (y0 + radix * y1)
                = value or sz + value x sz * value y 2 };
    assert { forall j. (j < r.offset \/ r.offset + sz + 1 <= j) ->
             (pelts r)[j] = (pelts (old r))[j]
             by (pelts r)[j] = (pelts r)[j] at A2 = (pelts (old r))[j] };
    c

end

module Mul_basecase

  use int.Int
  use mach.int.Int32
  use import mach.int.UInt64GMP as Limb
  use int.Power
  use ref.Ref
  use mach.c.C
  use array.Array
  use map.Map
  use types.Types
  use types.Int32Eq
  use types.UInt64Eq
  use lemmas.Lemmas
  use util.Util
  use add.Add
  use Mul

  let wmpn_mul_basecase (r x:t) (sx:int32) (y:t) (sy:int32) : unit
    requires { 0 < sy <= sx }
    requires { valid x sx }
    requires { valid y sy }
    requires { valid r (sy + sx) }
    requires { writable r }
    writes   { r.data.elts }
    ensures  { value r (sy + sx) = value x sx * value y sy }
    ensures { forall j. (j < offset r \/ offset r + (sy + sx) <= j)
              -> (pelts r)[j] = old (pelts r)[j] }
    (*ensures  { result = (pelts r)[r.offset + sx + sy - 1] }*)
  =
    let c = wmpn_mul_1 r x sx (C.get y) in
    value_sub_update_no_change (pelts r) (r.offset + p2i sx)
                               r.offset (r.offset + p2i sx - 1) c;
    set_ofs r sx c;
    value_sub_tail (pelts r) r.offset (r.offset + p2i sx);
    assert { value r (sx + 1) = value x sx * value y 1
             by (pelts y)[offset y] = value y 1
             so value r sx + power radix sx * c = value x sx * value y 1 };
    let ref rp = C.incr r 1 in
    let ref vp = C.incr y 1 in
    let ghost ref i = 1 in
    let ref vn = sy - 1 in
    while vn >= 2 do
      invariant { 1 <= i <= sy }
      invariant { i = sy - vn }
      invariant { value r (i + sx) = value x sx * value y i }
      invariant { rp.offset = r.offset + i }
      invariant { plength rp = plength r }
      invariant { rp.min = r.min }
      invariant { rp.max = r.max }
      invariant { pelts rp = pelts r }
      invariant { writable rp }
      invariant { vp.offset = y.offset + i }
      invariant { plength vp = plength y }
      invariant { vp.min = y.min }
      invariant { vp.max = y.max }
      invariant { pelts vp = pelts y }
      invariant { forall j. (j < offset r \/ offset r + (sy + sx) <= j)
                -> (pelts r)[j] = old (pelts r)[j] }
      variant { vn }
      label StartLoop in
      value_concat r i (i + sx);
      assert { value rp sx =
               value_sub (pelts r) (r.offset + i) (r.offset + (i + sx)) };
      let res =  wmpn_addmul_2 rp x sx vp in
      assert { value rp (sx+1) + power radix (sx+1) * res
              = value (rp at StartLoop) sx + value x sx * value vp 2 };
      assert { MapEq.map_eq_sub (pelts r) (pelts r at StartLoop)
                r.offset rp.offset
                by rp.offset = r.offset + i
                so forall j. r.offset <= j < rp.offset
                   ->
                   (j < rp.offset
                    so (pelts rp)[j] = (pelts rp at StartLoop)[j]
                         = (pelts r at StartLoop)[j]) };
      label BeforeCarry in
      value_sub_update_no_change (pelts r) (rp.offset + p2i sx + 1)
                        r.offset  (r.offset + p2i i) res;
      set_ofs rp (sx+1) res;
      assert { value rp (sx+1) = value (rp at BeforeCarry) (sx+1) };
      assert { value r i = value (r at BeforeCarry) i
             = value (r at StartLoop) i };
      value_tail r (i + sx + 1);
      value_tail rp (sx + 1);
      value_concat y i (i+2);
      assert { value vp 2
               = value_sub (pelts y) (y.offset + i) (y.offset + (i+2)) };
      value_concat r i (i + sx + 2);
(*      assert { value_sub (pelts r) (r.offset + i) (r.offset + (i + sx + 1))
               = value rp (sx + 1) };
      assert { (pelts r)[r.offset + (i + sx + 1)] = res };*)
      assert { value rp (sx + 2) = value (rp at StartLoop) sx
               + value x sx * value vp 2 };
      assert { value x sx * value y (i+2)
               = value x sx * value y i
                 + power radix i * (value x sx * value vp 2) };
      assert { value r (i + sx + 2) = value x sx * value y (i+2)
               by value r (i + sx + 2)
                  = value r i
                    + power radix i * (value_sub (pelts r) (offset r + i)
                                                 (offset r + i + sx + 2))
                  = value r i + power radix i * value rp (sx + 2)
                  = value r i
                    + power radix i
                      * (value (rp at StartLoop) sx + value x sx * value vp 2)
                  = value (r at StartLoop) i
                    + power radix i
                      * (value (rp at StartLoop) sx + value x sx * value vp 2)
                  = value (r at StartLoop) i
                    + power radix i * value (rp at StartLoop) sx
                    + power radix i * (value x sx * value vp 2)
                  = value (r at StartLoop) (sx + i)
                    + power radix i * (value x sx * value vp 2)
                  = value x sx * value y i
                    + power radix i * value x sx * value vp 2
                  = value x sx * (value y i + power radix i * value vp 2)
                  = value x sx * value y (i+2) };
      i <- i + 2;
      rp <- C.incr rp 2;
      vp <- C.incr vp 2;
      vn <- vn - 2;
    done;
    while vn >= 1 do
      invariant { 1 <= i <= sy }
      invariant { i = sy - vn }
      invariant { value r (i + sx) = value x sx * value y i }
      invariant { rp.offset = r.offset + i }
      invariant { plength rp = plength r }
      invariant { rp.min = r.min }
      invariant { rp.max = r.max }
      invariant { pelts rp = pelts r }
      invariant { writable rp }
      invariant { vp.offset = y.offset + i }
      invariant { plength vp = plength y }
      invariant { vp.min = y.min }
      invariant { vp.max = y.max }
      invariant { pelts vp = pelts y }
      invariant { forall j. (j < offset r \/ offset r + (sy + sx) <= j)
                -> (pelts r)[j] = old (pelts r)[j] }
      variant { vn }
      label StartLoop in
      value_concat r i (i + sx);
      assert { value rp sx =
               value_sub (pelts r) (r.offset + i) (r.offset + (i + sx)) };
      let ghost ly = pure { (pelts y)[y.offset + i] } in
      let res =  wmpn_addmul_1 rp x sx (C.get vp) in
      assert { value rp sx + power radix sx * res
              = value (rp at StartLoop) sx + value x sx * ly };
      assert { MapEq.map_eq_sub (pelts r) (pelts r at StartLoop)
                r.offset rp.offset
                by rp.offset = r.offset + i
                so forall j. r.offset <= j < rp.offset
                   ->
                   (j < rp.offset
                    so (pelts rp)[j] = (pelts rp at StartLoop)[j]
                         = (pelts r at StartLoop)[j]) };
      label BeforeCarry in
      value_sub_update_no_change (pelts r) (rp.offset + p2i sx)
                        r.offset  (r.offset + p2i i) res;
      set_ofs rp sx res;
      assert { value rp sx = value (rp at BeforeCarry) sx };
      assert { value r i = value (r at BeforeCarry) i
             = value (r at StartLoop) i };
      value_tail r (i + sx);
      value_tail y i;
      value_concat r i (i+sx);
      assert { value_sub (pelts r) (r.offset + i) (r.offset+(i+sx))
               = value rp sx };
      assert { (pelts r)[r.offset + (i+sx)] = res };
      assert { value x sx * value y (i+1)
               = value x sx * value y i
                 + power radix i * (value x sx * ly) };
             (*nonlinear*)
      assert { value r (i + sx + 1) = value x sx * value y (i+1) };
      i <- i + 1;
      rp <- C.incr rp 1;
      vp <- C.incr vp 1;
      vn <- vn - 1;
    done

end