module Sub

  use mach.c.C
  use lemmas.Lemmas
  use int.Int
  use mach.int.Int32
  use import mach.int.UInt64GMP as Limb
  use types.Types
  use array.Array
  use map.Map
  use ptralias.Alias
  use int.Power
  use map.MapEq

  (** `wmpn_sub_n r x y sz` subtracts `(y, sz)` from `(x, sz)` and
      writes the result to `(r, sz)`. Returns borrow, either 0 or
      1. Corresponds to `mpn_sub_n`. *)
  let wmpn_sub_n (r x y:ptr limb) (sz:int32) : limb
    requires { 0 <= sz }
    requires { valid r sz /\ valid x sz /\ valid y sz }
    requires { offset r = offset x \/ offset r + sz <= offset x
               \/ offset x + sz <= offset r }
    requires { offset r = offset y \/ offset r + sz <= offset y
               \/ offset y + sz <= offset r }
    requires { r.data = x.data = y.data }
    requires { writable r }
    alias    { r.data with x.data, r.data with y.data }
    ensures  { value r sz - power radix sz * result
               = old (value x sz - value y sz) }
    ensures  { 0 <= result <= 1 }
    ensures  { forall j. (j < offset r \/ offset r + sz <= j)
                         -> (pelts r)[j] = old (pelts r)[j] }
    ensures  { offset r = offset x \/ map_eq_sub (pelts x) (old pelts x)
                                                 (offset x) (offset x + sz) }
    ensures  { offset r = offset y \/ map_eq_sub (pelts y) (old pelts y)
                                                 (offset y) (offset y + sz) }
    writes { r.data.elts }
  =
    let ref lx = 0 in
    let ref ly = 0 in
    let ref b = 0 in
    let ref i = 0 in
    let ghost ox = pure { x } in
    let ghost oy = pure { y } in
    while i < sz do
      variant { sz - i }
      invariant { 0 <= i <= sz }
      invariant { value r i - power radix i * b = value ox i - value oy i }
      invariant { forall j. (j < offset r \/ offset r + i <= j)
                            -> (pelts r)[j] = old (pelts r)[j] }
      invariant { pelts x = pelts r = pelts y }
      invariant { offset r = offset x \/ map_eq_sub (pelts x) (old pelts x)
                                                   (offset x) (offset x + sz) }
      invariant { offset r = offset y \/ map_eq_sub (pelts y) (old pelts y)
                                                    (offset y) (offset y + sz) }
      invariant { 0 <= b <= 1 }
      label StartLoop in
      lx <- get_ofs x i;
      ly <- get_ofs y i;
      let ghost olx = get_ofs ox i in
      let ghost oly = get_ofs oy i in
      assert { lx = olx /\ ly = oly };
      let res, borrow = sub_with_borrow lx ly b in
      value_sub_update_no_change (pelts r) (offset r + int32'int i)
                                 (offset r) (offset r + int32'int i) res;
      set_ofs r i res;
      assert { value r i = value r i at StartLoop };
      assert { value r i - power radix i * b = value ox i - value oy i };
      b <- borrow;
      value_tail r i;
      value_tail ox i;
      value_tail oy i;
      assert { value r (i+1) - power radix (i+1) * b
               = value ox (i+1) - value oy (i+1) };
      i <- i+1;
    done;
    b

  (** `wmpn_sub r x sx y sy` subtracts `(y,sy)` from `(x, sx)` and writes the
      result in `(r, sx)`. `sx` must be greater than or equal to
      `sy`. Returns borrow, either 0 or 1. Corresponds to `mpn_sub`. *)
  let wmpn_sub (r x:ptr limb) (sx:int32) (y:ptr limb) (sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid r sx /\ valid x sx /\ valid y sy }
    requires { offset r = offset x \/ offset r + sx <= offset x
               \/ offset x + sx <= offset r }
    requires { offset r = offset y \/ offset r + sx <= offset y
               \/ offset y + sy <= offset r }
    requires { r.data = x.data = y.data }
    requires { writable r }
    alias    { r.data with x.data }
    alias    { r.data with y.data }
    alias    { x.data with y.data }
    ensures  { value r sx - power radix sx * result
               = old (value x sx - value y sy) }
    ensures  { 0 <= result <= 1 }
    ensures  { forall j. (j < offset r \/ offset r + sx <= j)
                         -> (pelts r)[j] = old (pelts r)[j] }
    ensures  { offset r = offset x \/ map_eq_sub (pelts x) (old pelts x)
                                                 (offset x) (offset x + sx) }
    ensures  { offset r = offset y \/ map_eq_sub (pelts y) (old pelts y)
                                                 (offset y) (offset y + sy) }
    writes { r.data.elts }
  =
    let ref lx = 0 in
    let ox = pure { x } in
    let oy = pure { y } in
    let ref b = wmpn_sub_n r x y sy in
    let ref i = sy in
    assert { offset r <> offset x ->
             forall j. offset x <= j < offset x + sx ->
                       (pelts x)[j] = (pelts ox)[j]
             by offset r + sy <= j \/ j < offset r };
    if (b <> 0)
    then begin
      while i < sx do
        variant { sx - i }
        invariant { sy <= i <= sx }
        invariant { value r i - power radix i * b
                    = value ox i - value oy sy }
        invariant { forall j. (j < offset r \/ offset r + i <= j)
                              -> (pelts r)[j] = old (pelts r)[j] }
        invariant { pelts x = pelts r }
        invariant { pelts y = pelts r }
        invariant { offset r = offset x \/ map_eq_sub (pelts x) (old pelts x)
                                                   (offset x) (offset x + sx) }
        invariant { offset r = offset y \/ map_eq_sub (pelts y) (old pelts y)
                                                   (offset y) (offset y + sy) }
        invariant { 0 <= b <= 1 }
        invariant { i = sx \/ b = 1 }
        label StartLoop in
        assert { b = 1 };
        lx <- get_ofs x i;
        let ghost olx = get_ofs ox i in
        assert { lx = olx };
        let res = sub_mod lx 1 in
        value_sub_update_no_change (pelts r) (r.offset + int32'int i)
                                   r.offset (r.offset + int32'int i) res;
        set_ofs r i res;
        assert { value r i = value r i at StartLoop };
        value_tail r i;
        value_tail ox i;
        i <- i+1;
        if lx <> 0
        then begin
          b <- 0;
          assert { res = lx - 1 };
          assert { value r i = value ox i - value oy sy };
          break
        end
        else begin
          assert { res = radix - 1 };
          assert { value r i - power radix i * b = value ox i - value oy sy };
        end
      done
    end;
    while i < sx do
        variant { sx - i }
        invariant { sy <= i <= sx }
        invariant { i = sx \/ b = 0 }
        invariant { value r i - power radix i * b
                    = value ox i - value oy sy }
        invariant { forall j. (j < offset r \/ offset r + i <= j)
                              -> (pelts r)[j] = old (pelts r)[j] }
        invariant { pelts x = pelts r }
        invariant { pelts y = pelts r }
        invariant { offset r = offset x \/ map_eq_sub (pelts x) (old pelts x)
                                                   (offset x) (offset x + sx) }
        invariant { offset r = offset y \/ map_eq_sub (pelts y) (old pelts y)
                                                   (offset y) (offset y + sy) }
        invariant { 0 <= b <= 1 }
        label StartLoop2 in
        assert { b = 0 by i < sx };
        lx <- get_ofs x i;
        let ghost olx = get_ofs ox i in
        assert { olx = lx };
        value_sub_update_no_change (pelts r) (r.offset + int32'int i)
                                   r.offset (r.offset + int32'int i) lx;
        set_ofs r i lx;
        assert { value r i = value r i at StartLoop2 };
        value_tail r i;
        value_tail ox i;
        assert { value r (i+1) - power radix (i+1) * b
                 = value ox (i+1) - value oy sy };
        i <- i+1
     done;
     b

  let sub_n (*[@extraction:inline]*) (r x y:ptr limb) (sz:int32) : limb
    requires { 0 <= sz }
    requires { valid r sz /\ valid x sz /\ valid y sz }
    requires { writable r }
    ensures  { value r sz - power radix sz * result
               = old (value x sz - value y sz) }
    ensures  { 0 <= result <= 1 }
    ensures { forall j. (j < offset r \/ offset r + sz <= j)
              -> (pelts r)[j] = old (pelts r)[j] }
    ensures { forall j. (pelts x)[j] = old (pelts x)[j] }
    ensures { forall j. (pelts y)[j] = old (pelts y)[j] }
    ensures { value x sz = old value x sz }
    ensures { value y sz = old value y sz }
    ensures  { min x = old min x /\ max x = old max x
               /\ plength x = old plength x }
    ensures  { min y = old min y /\ max y = old max y
               /\ plength y = old plength y }
    ensures  { min r = old min r /\ max r = old max r
               /\ plength r = old plength r }
  =
    let ghost ox = pure { x } in
    let ghost oy = pure { y } in
    let nr, nx, ny, m = open_sep r x sz y sz in
    label Sub in
    let res = wmpn_sub_n nr nx ny sz in
    let ghost onx = pure { nx } in
    let ghost ony = pure { ny } in
    close_sep r x sz y sz nr nx ny m;
    assert { forall j. 0 <= j < sz
                       -> (pelts x)[offset x + j] = (pelts ox)[offset x + j]
             by (pelts x)[offset x + j]
                = (pelts onx)[offset onx + j]
                = (pelts nx at Sub)[offset onx + j]
                = (pelts ox)[offset x + j] };
    assert { forall j. 0 <= j < sz
                       -> (pelts y)[offset y + j] = (pelts oy)[offset y + j]
             by (pelts y)[offset y + j]
                = (pelts ony)[offset ony + j]
                = (pelts ny at Sub)[offset ony + j]
                = (pelts oy)[offset y + j] };
    res

  let sub_n_rx (*[@extraction:inline]*) (x y:ptr limb) (sz:int32) : limb
    requires { 0 <= sz }
    requires { valid x sz /\ valid y sz }
    requires { writable x }
    ensures  { value x sz - power radix sz * result
               = old (value x sz - value y sz) }
    ensures  { 0 <= result <= 1 }
    ensures  { forall j. (j < offset x \/ offset x + sz <= j)
                      -> (pelts x)[j] = old (pelts x)[j] }
    ensures  { forall j. (pelts y)[j] = old (pelts y)[j] }
    ensures  { value y sz = old value y sz }
    ensures  { min x = old min x /\ max x = old max x
               /\ plength x = old plength x }
    ensures  { min y = old min y /\ max y = old max y
               /\ plength y = old plength y }
  = let ghost oy = pure { y } in
    let nr, nx, ny, m = open_rx x sz y sz in
    label Sub in
    let res = wmpn_sub_n nr nx ny sz in
    let ghost ony = pure { ny } in
    close_rx x sz y sz nr nx ny m;
    assert { forall j. 0 <= j < sz
                       -> (pelts y)[offset y + j] = (pelts oy)[offset y + j]
             by (pelts y)[offset y + j]
                = (pelts ony)[offset ony + j]
                = (pelts ny at Sub)[offset ony + j]
                = (pelts oy)[offset y + j] };
    value_sub_frame_shift (pelts y) (pelts oy) (offset y) (offset y)
                          (int32'int sz);
    res

  let sub_n_ry (*[@extraction:inline]*) (x y:ptr limb) (sz:int32) : limb
    requires { 0 <= sz }
    requires { valid x sz /\ valid y sz }
    requires { writable y }
    ensures  { value y sz - power radix sz * result
               = old (value x sz - value y sz) }
    ensures  { 0 <= result <= 1 }
    ensures  { forall j. (j < offset y \/ offset y + sz <= j)
               -> (pelts y)[j] = old (pelts y)[j] }
    ensures  { forall j. (pelts x)[j] = (old pelts x)[j] }
    ensures  { value x sz = old value x sz }
    ensures  { min x = old min x /\ max x = old max x
               /\ plength x = old plength x }
    ensures  { min y = old min y /\ max y = old max y
               /\ plength y = old plength y }
  = let ghost ox = pure { x } in
    let nr, ny, nx, m = open_rx y sz x sz in
    label Sub in
    let res = wmpn_sub_n nr nx ny sz in
    let ghost onx = pure { nx } in
    close_rx y sz x sz nr ny nx m;
    assert { forall j. 0 <= j < sz
                       -> (pelts x)[offset x + j] = (pelts ox)[offset x + j]
             by (pelts x)[offset x + j]
                = (pelts onx)[offset onx + j]
                = (pelts nx at Sub)[offset onx + j]
                = (pelts ox)[offset x + j] };
    value_sub_frame_shift (pelts x) (pelts ox)
                          (offset x) (offset x) (int32'int sz);
    res

  let sub (*[@extraction:inline]*) (r x:ptr limb) (sx:int32)
                                      (y:ptr limb) (sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid r sx /\ valid x sx /\ valid y sy }
    requires { writable r }
    ensures  { value r sx - power radix sx * result
               = old (value x sx - value y sy) }
    ensures  { 0 <= result <= 1 }
    ensures { forall j. (j < offset r \/ offset r + sx <= j)
              -> (pelts r)[j] = old (pelts r)[j] }
    ensures { forall j. (pelts x)[j] = old (pelts x)[j] }
    ensures { forall j. (pelts y)[j] = old (pelts y)[j] }
    ensures { value x sx = old value x sx }
    ensures { value y sy = old value y sy }
    ensures  { min x = old min x /\ max x = old max x
               /\ plength x = old plength x }
    ensures  { min y = old min y /\ max y = old max y
               /\ plength y = old plength y }
    ensures  { min r = old min r /\ max r = old max r
               /\ plength r = old plength r }
  =
    let ghost ox = pure { x } in
    let ghost oy = pure { y } in
    let nr, nx, ny, m = open_sep r x sx y sy in
    label Sub in
    let res = wmpn_sub nr nx sx ny sy in
    let ghost onx = pure { nx } in
    let ghost ony = pure { ny } in
    close_sep r x sx y sy nr nx ny m;
    assert { forall j. 0 <= j < sx
                       -> (pelts x)[offset x + j] = (pelts ox)[offset x + j]
             by (pelts x)[offset x + j]
                = (pelts onx)[offset onx + j]
                = (pelts nx at Sub)[offset onx + j]
                = (pelts ox)[offset x + j] };
    assert { forall j. 0 <= j < sy
                       -> (pelts y)[offset y + j] = (pelts oy)[offset y + j]
             by (pelts y)[offset y + j]
                = (pelts ony)[offset ony + j]
                = (pelts ny at Sub)[offset ony + j]
                = (pelts oy)[offset y + j] };
    res

  let sub_rx (*[@extraction:inline]*) (x:ptr limb) (sx:int32)
                                  (y:ptr limb) (sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid x sx /\ valid y sy }
    requires { writable x }
    ensures  { value x sx - power radix sx * result
               = old (value x sx - value y sy) }
    ensures  { 0 <= result <= 1 }
    ensures  { forall j. (j < offset x \/ offset x + sx <= j)
               -> (pelts x)[j] = old (pelts x)[j] }
    ensures  { forall j. (pelts y)[j] = (old pelts y)[j] }
    ensures  { value y sy = old value y sy }
    ensures  { min x = old min x /\ max x = old max x
               /\ plength x = old plength x }
    ensures  { min y = old min y /\ max y = old max y
               /\ plength y = old plength y }
  =
    let ghost oy = pure { y } in
    let nr, nx, ny, m = open_rx x sx y sy in
    label Sub in
    let res = wmpn_sub nr nx sx ny sy in
    let ghost ony = pure { ny } in
    close_rx x sx y sy nr nx ny m;
    assert { forall j. 0 <= j < sy
                       -> (pelts y)[offset y + j] = (pelts oy)[offset y + j]
             by (pelts y)[offset y + j]
                = (pelts ony)[offset ony + j]
                = (pelts ny at Sub)[offset ony + j]
                = (pelts oy)[offset y + j] };
    value_sub_frame_shift (pelts y) (pelts oy)
                          (offset y) (offset oy) (int32'int sy);
    res

  let sub_ry (*[@extraction:inline]*) (x:ptr limb) (sx:int32)
                                  (y:ptr limb) (sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid x sx /\ valid y sx }
    requires { writable y }
    ensures  { value y sx - power radix sx * result
               = old (value x sx - value y sy) }
    ensures  { 0 <= result <= 1 }
    ensures  { forall j. (j < offset y \/ offset y + sx <= j)
               -> (pelts y)[j] = old (pelts y)[j] }
    ensures  { forall j. (pelts x)[j] = (old pelts x)[j] }
    ensures  { value x sx = old value x sx }
    ensures  { min x = old min x /\ max x = old max x
               /\ plength x = old plength x }
    ensures  { min y = old min y /\ max y = old max y
               /\ plength y = old plength y }
  =
    let ghost ox = pure { x } in
    let ghost oy = pure { y } in
    let nr, ny, nx, m = open_rx y sx x sx in
    label Sub in
    value_sub_frame_shift (pelts ny) (pelts oy)
                          (offset ny) (offset y) (int32'int sy);
    assert { value ny sy = old (value y sy) };
    assert { value nx sx = old (value x sx) };
    let res = wmpn_sub nr nx sx ny sy in
    let ghost onx = pure { nx } in
    close_rx y sx x sx nr ny nx m;
    assert { forall j. 0 <= j < sx
                       -> (pelts x)[offset x + j] = (pelts ox)[offset x + j]
             by (pelts x)[offset x + j]
                = (pelts onx)[offset onx + j]
                = (pelts nx at Sub)[offset onx + j]
                = (pelts ox)[offset x + j] };
    value_sub_frame_shift (pelts x) (pelts ox)
                          (offset x) (offset x) (int32'int sx);
    res


end

(* TODO remove when sub_rx can be used instead *)
module SubOld

  use mach.c.C
  use lemmas.Lemmas
  use int.Int
  use mach.int.Int32
  use import mach.int.UInt64GMP as Limb
  use types.Types
  use array.Array
  use map.Map
  use int.Power
  use map.MapEq
  use ptralias.Alias
  use ref.Ref

  let wmpn_sub_n (r x y:t) (sz:int32) : limb
    requires { valid x sz }
    requires { valid y sz }
    requires { valid r sz }
    requires { writable r }
    ensures { 0 <= result <= 1 }
    ensures { value r sz - power radix sz * result
              = value x sz - value y sz }
    ensures { forall j. (j < offset r \/ offset r + sz <= j)
              -> (pelts r)[j] = old (pelts r)[j] }
    writes { r.data.elts }
  =
    let lx = ref 0 in
    let ly = ref 0 in
    let b = ref 0 in
    let i = ref 0 in
    while !i < sz do
      variant { sz - !i }
      invariant { 0 <= !i <= sz }
      invariant { value r !i - (power radix !i) * !b
                  = value x !i - value y !i }
      invariant { forall j. (j < offset r \/ offset r + sz <= j)
                -> (pelts r)[j] = old (pelts r)[j] }
      invariant { 0 <= !b <= 1 }
      label StartLoop in
      lx := get_ofs x !i;
      ly := get_ofs y !i;
      let res, borrow = sub_with_borrow !lx !ly !b in
      set_ofs r !i res;
      assert { value r !i - (power radix !i) * !b =
      value x !i - value y !i };
      b := borrow;
      value_tail r !i;
      value_tail x !i;
      value_tail y !i;
      assert { value r (!i+1) - (power radix (!i+1)) * !b
                  = value x (!i+1) - value y (!i+1) };
      i := !i + 1;
      done;
      !b

  let wmpn_sub (r x:t) (sx:int32) (y:t) (sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid x sx }
    requires { valid y sy }
    requires { valid r sx }
    requires { writable r }
    ensures { value r sx  - power radix sx * result
              = value x sx - value y sy }
    ensures { forall j. (j < offset r \/ offset r + sx <= j)
              -> (pelts r)[j] = old (pelts r)[j] }
    ensures { 0 <= result <= 1 }
    writes { r.data.elts }
  =
    let ref lx = 0 in
    let ref b = wmpn_sub_n r x y sy in
    let ref i = sy in
    if b <> 0
    then begin
      while i < sx do
        variant { sx - i }
        invariant { sy <= i <= sx }
        invariant { value r i - power radix i * b =
                    value x i - value y sy }
        invariant { forall j. (j < offset r \/ offset r + sx <= j)
                  -> (pelts r)[j] = old (pelts r)[j] }
        invariant { 0 <= b <= 1 }
        invariant { i = sx \/ b = 1 }
        assert { b = 1 };
        lx <- get_ofs x i;
        let res = sub_mod lx 1 in
        value_sub_update_no_change (pelts x) (x.offset + int32'int i)
                                             (x.offset + int32'int i + 1)
                                             (x.offset + int32'int sx) res;
        set_ofs r i res;
        assert { value r i - (power radix i) * b = value x i - value y sy };
        value_tail r i;
        value_tail x i;
        i <- i + 1;
        if (lx <> 0)
        then begin
          b <- 0;
          assert { res = lx - 1 };
          assert { value r i = value x i - value y sy};
          break;
        end
        else begin
          assert { res = radix - 1 };
          assert { value r i - power radix i = value x i - value y sy }
        end
      done
    end;
    while i < sx do
      variant { sx - i }
      invariant { sy <= i <= sx }
      invariant { i = sx \/ b = 0 }
      invariant { value r i - power radix i * b =
                value x i - value y sy }
      invariant { forall j. (j < offset r \/ offset r + sx <= j)
                -> (pelts r)[j] = old (pelts r)[j] }
      assert { b = 0 by i < sx };
      lx <- get_ofs x i;
      set_ofs r i lx;
      value_tail r i;
      value_tail x i;
      assert { value r i = value x i - value y sy };
      assert { value r (i+1) - power radix (i+1) * b
               = value x (i+1) - value y sy };
      i <- i + 1;
    done;
    b

  let wmpn_sub_n_in_place (x y:t) (sz:int32) : limb
    requires { 0 <= sz }
    requires { valid x sz }
    requires { valid y sz }
    requires { writable x }
    ensures { value x sz  - power radix sz * result
              = value (old x) sz - value y sz }
    ensures { 0 <= result <= 1 }
    writes { x.data.elts }
    ensures { forall j. j < x.offset \/ x.offset + sz <= j ->
              (pelts x)[j] = (pelts (old x))[j] }
  =
    let ghost ox = { x } in
    let ref lx = 0 in
    let ref ly = 0 in
    let ref b = 0 in
    let ref i = 0 in
    while i < sz do
      variant { sz - i }
      invariant { 0 <= i <= sz }
      invariant { value x i - power radix i * b =
                  value ox i - value y i }
      invariant { 0 <= b <= 1 }
      invariant { forall j. i <= j < sz ->
                  (pelts x)[x.offset + j] = (pelts ox)[x.offset + j] }
      invariant { forall j. j < x.offset \/ x.offset + sz <= j ->
                  (pelts x)[j] = (pelts (old x))[j] }
      label StartLoop in
      lx <- get_ofs x i;
      assert { lx = (pelts ox)[ox.offset + i] };
      ly <- get_ofs y i;
      let res, borrow = sub_with_borrow lx ly b in
      value_sub_update_no_change (pelts x) (offset x + int32'int i)
                                 (offset x) (offset x + int32'int i) res;
      set_ofs x i res;
      assert { forall j. i < j < sz ->
                 (pelts x)[x.offset + j]
                 = (pelts ox)[x.offset + j]
                 by (pelts x)[x.offset + j]
                 = (pelts (x at StartLoop))[x.offset + j]
                 = (pelts ox)[x.offset + j]};
      assert { value x i - power radix i * b = value ox i - value y i };
      b <- borrow;
      value_tail ox i;
      value_tail x i;
      value_tail y i;
      assert { value x (i+1) - power radix (i+1) * b =
              value ox (i+1) - value y (i+1) };
      i <- i + 1;
    done;
    b

  let wmpn_sub_in_place (x:t) (sx:int32) (y:t) (sy:int32) : limb
    requires { 0 <= sy <= sx }
    requires { valid x sx }
    requires { valid y sy }
    requires { writable x }
    ensures { value x sx  - power radix sx * result
              = value (old x) sx - value y sy }
    ensures { 0 <= result <= 1 }
    writes { x.data.elts }
    ensures { forall j. j < x.offset \/ x.offset + sx <= j ->
              (pelts x)[j] = (pelts (old x))[j] }
  =
    let ghost ox = { x } in
    let ref lx = 0 in
    let ref b = wmpn_sub_n_in_place x y sy in
    let ref i = sy in
    if not (b = 0)
    then begin
      while i < sx do
        variant { sx - i }
        invariant { sy <= i <= sx }
        invariant { value x i - power radix i * b =
                    value ox i - value y sy }
        invariant { 0 <= b <= 1 }
        invariant { i = sx \/ b = 1 }
        invariant { forall j. i <= j < sx ->
                    (pelts x)[x.offset + j] = (pelts ox) [x.offset + j] }
        invariant { forall j. j < x.offset \/ x.offset + sx <= j ->
                    (pelts x)[j] = (pelts (old x))[j] }
        assert { b = 1 };
        lx <- get_ofs x i;
        assert { lx = (pelts ox)[ox.offset + i] };
        let res = sub_mod lx 1 in
        value_sub_update_no_change (pelts x) (x.offset + int32'int i)
                                             (x.offset + int32'int i + 1)
                                             (x.offset + int32'int sx) res;
        set_ofs x i res;
        assert { value x i - power radix i * b = value ox i - value y sy };
        assert { forall j. i < j < sx ->
                 (pelts x)[x.offset + j] = (pelts ox) [x.offset + j] };
        value_tail ox i;
        value_tail x i;
        i <- i + 1;
        if not (lx =  0)
        then begin
          b <- 0;
          assert { res = lx - 1 };
          assert { value x i = value ox i - value y sy };
          break;
        end
        else begin
          assert { res = radix - 1 };
          assert { value x i - power radix i = value ox i - value y sy }
        end
      done
    end;
    assert { forall j. x.offset + i <= j < x.offset + sx
             -> (pelts x)[j] = (pelts ox)[j]
             by i <= j - x.offset < sx
             so (pelts x)[x.offset + (j - x.offset)]
                = (pelts ox)[x.offset + (j - x.offset)] };
    value_sub_frame (pelts x) (pelts ox) (x.offset + int32'int i)
                                         (x.offset + int32'int sx);
    value_concat x i sx;
    value_concat ox i sx;
    b

end