void Int::Div(Int *a, Int *mod) {
  if (a->IsGreater(this)) {
    if (mod) mod->Set(this);
    CLEAR();
    return;
  }

  if (a->IsZero()) {
    printf("Divide by 0!\n");
    return;
  }

  if (IsEqual(a)) {
    if (mod) mod->CLEAR();
    Set(&_ONE);
    return;
  }

  Int rem(this);
  Int d(a);
  Int dq;
  CLEAR();

  // Size
  uint32_t dSize = d.GetSize64();
  uint32_t tSize = rem.GetSize64();
  uint32_t qSize = tSize - dSize + 1;

  // D1 normalize the divisor (d!=0)
  uint32_t shift = (uint32_t)LZC(d.bits64[dSize - 1]);
  d.ShiftL(shift);
  rem.ShiftL(shift);

  uint64_t _dh = d.bits64[dSize - 1];
  uint64_t _dl = (dSize > 1) ? d.bits64[dSize - 2] : 0;
  int sb = tSize - 1;

  // D2 Initialize j
  for (int j = 0; j < (int)qSize; j++) {
    // D3 Estimate qhat
    uint64_t qhat = 0;
    uint64_t qrem = 0;
    int skipCorrection = false;
    uint64_t nh = rem.bits64[sb - j + 1];
    uint64_t nm = rem.bits64[sb - j];

    if (nh == _dh) {
      qhat = ~0;
      qrem = nh + nm;
      skipCorrection = qrem < nh;
    } else {
      __uint128_t dividend = ((__uint128_t)nh << 64) | nm;
      qhat = dividend / _dh;
      qrem = dividend % _dh;
    }

    if (qhat == 0)
      continue;

    if (!skipCorrection) {
      // Correct qhat
      uint64_t nl = rem.bits64[sb - j - 1];

      uint64_t estProH;
      uint64_t estProL = _umul128(_dl, qhat, &estProH);
      if (isStrictGreater128(estProH, estProL, qrem, nl)) {
        qhat--;
        qrem += _dh;
        if (qrem >= _dh) {
          estProL = _umul128(_dl, qhat, &estProH);
          if (isStrictGreater128(estProH, estProL, qrem, nl))
            qhat--;
        }
      }
    }
 
    dq.Mult(&d, qhat);
    rem.ShiftL64BitAndSub(&dq, qSize - j - 1);
    if (rem.IsNegative()) {
      // Overflow
      rem.Add(&d);
      qhat--;
    }

    bits64[qSize - j - 1] = qhat;
  }

  if (mod) {
    rem.ShiftR(shift);
    mod->Set(&rem);
  }
}