module Zdiv

use int.Int
use int.ComputerDivision
use int.Power
use map.Map
use mach.int.Int32GMP
use ref.Ref
use mach.c.C
use lemmas.Lemmas
use util.Util
use util.UtilOld as OU
use ptralias.Alias
use compare.Compare
use import mach.int.UInt64GMP as Limb
use types.Types
use types.Int32Eq
use types.UInt64Eq
use int.Abs
use div.Div
use mpz.Z
use mpz.Zutil

let wmpz_tdiv_qr (quot rem num den: mpz_ptr) : unit
  requires { mpz.alloc[num] > 0 /\ mpz.alloc[den] > 0 }
  requires { mpz.alloc[quot] > 0 /\ mpz.alloc[rem] > 0 }
  requires { mpz.readers[num] = 0 /\ mpz.readers[den] = 0
             /\ mpz.readers[quot] = 0 /\ mpz.readers[rem] = 0 }
  requires { quot <> rem }
  requires { value_of den mpz <> 0 }
  requires { mpz.abs_size[num] <= max_int32 - 1 }
  ensures  { value_of quot mpz * (old value_of den mpz) + value_of rem mpz
             = old value_of num mpz }
  ensures  { 0 <= mpz.abs_value_of rem < old mpz.abs_value_of den }
  ensures  { forall x. x <> quot -> x <> rem -> mpz_unchanged x mpz (old mpz) }
  ensures  { mpz.readers[num] = 0 /\ mpz.readers[den] = 0
             /\ mpz.readers[quot] = 0 /\ mpz.readers[rem] = 0 }
=
  let ns = size_of num in
  let ds = size_of den in
  let nl = abs ns in
  let ref dl = abs ds in
  let ref ql = nl - dl + 1 in
  assert { dl <> 0 };
  if ql <= 0
  then begin
    assert { old mpz.abs_value_of num < old mpz.abs_value_of den };
    if not mpz_eq num rem
    then begin
      let ghost ompz = pure { mpz } in
      let np = get_read_ptr num in
      label Realloc in
      let rp = wmpz_realloc rem dl in
      wmpn_copyd_sep rp np nl;
      assert { value rp nl = value np nl at Realloc
               = old mpz.abs_value_of num };
      set_size rem ns rp;
      unchanged num mpz (pure { mpz at Realloc });
      release_reader num np;
      unchanged num mpz ompz;
      release_writer rem rp;
      assert { value_of rem mpz = old value_of num mpz };
      ghost (if not mpz_eq num quot then unchanged quot mpz ompz);
    end;
    let ghost mpz0 = pure { mpz } in
    set_size_0 quot;
    unchanged rem mpz mpz0;
    assert { value_of quot mpz = 0 };
    assert { value_of rem mpz = old value_of num mpz };
    return
  end;
  let ref d' = den in
  let ref n' = num in
  let ghost mpz' = pure { mpz } in
  let ghost ref cmemd = any mpz_mem in
  let ghost ref cmemn = any mpz_mem in
  let ghost ref clear_d = false in
  let ghost ref clear_n = false in
  assert { mpz.readers[d'] = 0 /\ mpz.readers[n'] = 0
           /\ mpz.readers[quot] = 0 /\ mpz.readers[rem] = 0 };
  begin
    (*ensures { d' <> rem }
    ensures { d' <> quot }
    ensures { n' <> rem }
    ensures { n' <> quot }
    ensures { forall x. mpz_unchanged x mpz mpz' }
    ensures { mpz.readers[d'] = 0 /\ mpz.readers[n'] = 0
              /\ mpz.readers[quot] = 0 /\ mpz.readers[rem] = 0 }
    ensures { mpz.alloc[d'] > 0 }
    ensures { mpz.abs_size[d'] = old mpz.abs_size[den] }
    ensures { mpz.sgn[d'] = old mpz.sgn[den] }
    ensures { mpz.abs_value_of[d'] = old mpz.abs_value_of[den] }
    ensures { mpz.alloc[n'] > 0 }
    ensures { mpz.abs_size[n'] = old mpz.abs_size[num] }
    ensures { mpz.sgn[n'] = old mpz.sgn[num] }
    ensures { mpz.abs_value_of[n'] = old mpz.abs_value_of[num] }*)
  if mpz_eq den rem || mpz_eq den quot
  then begin
    let dp = get_read_ptr den in
    let ghost mpzd = pure { mpz } in
    let nd, memd = wmpz_ptr_decl () in
    d' <- nd;
    cmemd <- memd;
    clear_d <- true;
    let tdp = salloc (UInt32.of_int32 dl) in
    set_alloc d' dl;
    wmpn_copyd_sep tdp dp dl;
    assert { value tdp dl = mpz'.abs_value_of den };
    set_ptr d' tdp;
    set_size d' ds tdp;
    release_writer d' tdp;
    assert { mpz.readers[d'] = 0 };
    unchanged den mpz mpzd;
    release_reader den dp;
    unchanged den mpz mpz';
    assert { forall x. x <> d' ->  mpz_unchanged x mpz mpz' };
    unchanged num mpz mpz';
    unchanged rem mpz mpz';
    unchanged quot mpz mpz';
  end;
  let ghost mpz'' = pure { mpz } in
  if mpz_eq num rem || mpz_eq num quot
  then begin
    let np = get_read_ptr num in
    let ghost mpzn = pure { mpz } in
    let nn, memn = wmpz_ptr_decl () in
    cmemn <- memn;
    n' <- nn;
    clear_n <- true;
    let tnp = salloc (UInt32.of_int32 nl) in
    set_alloc n' nl;
    wmpn_copyd_sep tnp np nl;
    assert { value tnp nl = mpz'.abs_value_of num };
    set_ptr n' tnp;
    set_size n' ns tnp;
    release_writer n' tnp;
    unchanged num mpz mpzn;
    release_reader num np;
    unchanged num mpz mpz'';
    assert { forall x. x <> n' -> mpz_unchanged x mpz mpz'' };
    unchanged den mpz mpz'';
    unchanged rem mpz mpz'';
    unchanged quot mpz mpz'';
    unchanged d' mpz mpz'';
  end
  end;
  assert { d' <> quot /\ d' <> rem };
  assert { n' <> quot /\ n' <> rem };
  assert { forall x. x <> n' -> x <> d' -> mpz_unchanged x mpz mpz' };
  let ghost mpz' = pure { mpz } in
  let qp = wmpz_realloc quot ql in
  unchanged rem mpz mpz';
  let ghost mpzq = pure { mpz } in
  let rp = wmpz_realloc rem dl in
  unchanged n' mpz mpz';
  unchanged d' mpz mpz';
  let np = get_read_ptr n' in
  let dp = get_read_ptr d' in
  let ghost mpzn = pure { mpz } in
  let ghost mpzd = pure { mpz } in
  value_tail dp (dl - 1);
  assert { (pelts dp)[offset dp + dl - 1] > 0
           by value dp dl = mpz.abs_value_of[d']
              >= power radix (dl - 1)
           so value dp dl
              = value dp (dl - 1)
                + power radix (dl - 1) * (pelts dp)[offset dp + dl - 1]
              < power radix (dl - 1)
                + power radix (dl - 1) * (pelts dp)[offset dp + dl - 1]
              = power radix (dl - 1) * (1 + (pelts dp)[offset dp + dl - 1])
           so 1 + (pelts dp)[offset dp + dl - 1] > 1 };
  wmpn_tdiv_qr qp rp 0 np nl dp dl;
  label Norm in
  assert { value qp ql * value dp dl + value rp dl = value np nl };
  assert { nl - dl - 1 >= 0 -> value qp ql >= power radix (nl - dl - 1)
           by value rp dl < value dp dl
           so value qp ql * value dp dl + value rp dl
              < value qp ql * value dp dl + value dp dl
              = (value qp ql + 1) * value dp dl
              < (value qp ql + 1) * power radix dl
           so value np nl >= power radix (nl - 1)
           so value dp dl < power radix dl
           so power radix (nl - 1) < power radix dl * (value qp ql + 1)
           so power radix (nl - 1) = power radix dl * power radix (nl - dl - 1)
           so power radix (nl - dl - 1) < value qp ql + 1 };
  assert { value qp ql = value qp (nl - dl + 1) };
  let qh = C.get_ofs qp (ql - 1) in
  value_tail qp (ql - 1);
  ql <- ql - (if qh = 0 then 1 else 0);
  assert { value qp ql = value qp (nl - dl + 1) };
  assert { ql >= 1 -> value qp ql >= power radix (ql - 1)
           by if qh = 0
              then ql = nl - dl
                   so value qp ql >= power radix (ql - 1)
              else value qp ql = value qp (ql - 1) + power radix (ql - 1) * qh
                   >= power radix (ql - 1) * qh
                   >= power radix (ql - 1) };
  normalize rp dl;
  let qs = if (bxor ns ds >= 0) then ql else - ql in
  let rs = if ns >= 0 then dl else -dl in
  assert { sgn_value qp qs * sgn_value dp ds + sgn_value rp rs = sgn_value np ns
           by if ns >= 0
              then if ds >= 0
                   then qs = ql /\ rs = dl
                        so sgn_value qp qs = value qp ql
                        so sgn_value rp rs = value rp dl
                   else qs = - ql /\ rs = dl
                        so sgn_value qp qs = - value qp ql
                        so sgn_value rp rs = value rp dl
              else if ds >= 0
                   then qs = -ql /\ rs = -dl
                        so sgn_value qp qs = - value qp ql
                        so sgn_value rp rs = - value rp dl
                   else qs = ql /\ rs = -dl
                        so sgn_value qp qs = value qp ql
                        so sgn_value rp rs = - value rp dl };
  assert { sgn_value np ns = old value_of num mpz };
  assert { sgn_value dp ds = old value_of den mpz };
  unchanged quot mpz mpzq;
  set_size quot qs qp;
  let ghost mpzr = pure { mpz } in
  assert { sgn_value qp qs = value_of quot mpz };
  set_size rem rs rp;
  unchanged quot mpz mpzr;
  assert { sgn_value rp rs = value_of rem mpz };
  assert { value_of quot mpz * (old value_of den mpz) + value_of rem mpz
           = old value_of num mpz };
  release_writer rem rp;
  release_writer quot qp;
  let ghost mpzrq = pure { mpz } in
  assert { value_of quot mpz * (old value_of den mpz) + value_of rem mpz
           = old value_of num mpz };
  unchanged n' mpz mpzn;
  unchanged d' mpz mpzd;
  release_reader n' np;
  release_reader d' dp;
  unchanged d' mpz mpz';
  unchanged n' mpz mpz';
  ghost (if clear_d then wmpz_tmp_clear d' cmemd);
  ghost (if clear_n then wmpz_tmp_clear n' cmemn);
  unchanged quot mpz mpzrq;
  unchanged rem mpz mpzrq;
  return

let wmpz_tdiv_qr_ui (quot rem dividend: mpz_ptr) (divisor: limb) : limb
  requires { mpz.alloc[dividend] > 0 }
  requires { mpz.alloc[quot] > 0 /\ mpz.alloc[rem] > 0 }
  requires { mpz.readers[dividend] = 0 /\ mpz.readers[quot] = 0
             /\ mpz.readers[rem] = 0 }
  requires { quot <> rem }
  requires { divisor > 0 }
  ensures  { value_of quot mpz * divisor + value_of rem mpz
             = old value_of dividend mpz }
  ensures  { result = mpz.abs_value_of rem }
  ensures  { 0 <= mpz.abs_value_of rem < divisor }
  ensures  { forall x. x <> quot -> x <> rem -> mpz_unchanged x mpz (old mpz) }
  ensures  { mpz.readers[dividend] = 0 /\ mpz.readers[quot] = 0
             /\ mpz.readers[rem] = 0 }
=
  let ghost ompz = pure { mpz } in
  let ns = size_of dividend in
  let ghost avn = mpz.abs_value_of dividend in
  if ns = 0
  then begin
    set_size_0 quot;
    set_size_0 rem;
    assert { value_of quot mpz = 0 /\ value_of rem mpz = 0 };
    return 0;
  end;
  let nn = abs ns in
  let qp = wmpz_realloc quot nn in
  ghost (if not mpz_eq quot dividend
         then unchanged dividend mpz ompz);
  unchanged rem mpz ompz;
  let ref qn = nn in
  let ref rl = 0 in
  if mpz_eq quot dividend
  then begin
    let np = salloc (UInt32.of_int32 nn) in
    wmpn_copyd_sep np qp nn;
    rl <- wmpn_divrem_1 qp np nn divisor
  end else begin
    let np = get_read_ptr dividend in
    assert { plength np >= nn };
    rl <- wmpn_divrem_1 qp np nn divisor;
    release_reader dividend np;
    ghost (if not mpz_eq rem dividend
           then unchanged rem mpz ompz);
  end;
  assert { value qp qn * divisor + rl = avn };
  if rl = 0
  then set_size_0 rem
  else begin
    let rs = if ns >= 0 then 1 else -1 in
    let rp = get_write_ptr rem in
    assert { plength rp >= 1 };
    C.set rp rl;
    assert { value rp 1 = rl };
    set_size rem rs rp;
    release_writer rem rp;
  end;
  assert { nn - 2 >= 0 -> value qp qn >= power radix (nn - 2)
           by 0 <= rl < divisor
           so value qp qn * divisor + rl
              < value qp qn * divisor + divisor
              = (value qp qn + 1) * divisor
              < (value qp qn + 1) * radix
           so avn >= power radix (nn - 1)
           so power radix (nn - 1) < radix * (value qp qn + 1)
           so power radix (nn - 1) = radix * power radix (nn - 2)
           so power radix (nn - 2) < value qp qn + 1 };
  let qh = C.get_ofs qp (nn - 1) in
  value_tail qp (qn - 1);
  qn <- qn - (if qh = 0 then 1 else 0);
  assert { value qp qn = value qp nn };
  assert { qn >= 1 -> value qp qn >= power radix (qn - 1)
           by if qh = 0
           then qn = nn - 1
                so value qp qn >= power radix (nn - 2)
           else power radix (qn - 1) * 1 <= power radix (qn - 1) * qh
                so value qp qn = value qp (qn - 1) + power radix (qn - 1) * qh
                >= power radix (qn - 1) * qh
                >= power radix (qn - 1) };
  let qs = if ns >= 0 then qn else -qn in
  assert { sgn_value qp qs * divisor
           + value_of rem mpz = value_of dividend ompz
           by if ns >= 0
              then qs = qn so value_of dividend ompz = avn
                   so value_of rem mpz = rl
                   so sgn_value qp qs * divisor
                       + value_of rem mpz = value_of dividend ompz
              else qs = - qn /\ value_of dividend ompz = - avn
                   so value_of rem mpz = - rl
                   so sgn_value qp qs = - value qp qn
                   so sgn_value qp qs * divisor
                       + value_of rem mpz
                       = - (value qp qn * divisor + rl)
                       = - avn
                       = value_of dividend ompz };
  let ghost mpzq = pure { mpz } in
  set_size quot qs qp;
  release_writer quot qp;
  unchanged rem mpz mpzq;
  return rl

end