static int do_div_mod()

in mysql_strings/decimal.cc [2197:2455]


static int do_div_mod(const decimal_t *from1, const decimal_t *from2,
                      decimal_t *to, decimal_t *mod, int scale_incr) {
  /*
    frac* - number of digits in fractional part of the number
    prec* - precision of the number
    intg* - number of digits in the integer part
    buf* - buffer having the actual number
    All variables ending with 0 - like frac0, intg0 etc are
    for the final result. Similarly frac1, intg1 etc are for
    the first number and frac2, intg2 etc are for the second number
   */
  int frac1 = ROUND_UP(from1->frac) * DIG_PER_DEC1, prec1 = from1->intg + frac1,
      frac2 = ROUND_UP(from2->frac) * DIG_PER_DEC1, prec2 = from2->intg + frac2,
      error = 0, i, intg0, frac0, len1, len2,
      dintg, /* Holds the estimate of number of integer digits in final result
              */
      div_mod = (!mod) /*true if this is division */;
  dec1 *buf0, *buf1 = from1->buf, *buf2 = from2->buf, *start1, *stop1, *start2,
              *stop2, *stop0, norm2, carry, dcarry, *tmp1;
  dec2 norm_factor, x, guess, y;

  if (mod) to = mod;

  sanity(to);

  /*
    removing all the leading zeroes in the second number. Leading zeroes are
    added later to the result.
   */
  i = ((prec2 - 1) % DIG_PER_DEC1) + 1;
  while (prec2 > 0 && *buf2 == 0) {
    prec2 -= i;
    i = DIG_PER_DEC1;
    buf2++;
  }
  if (prec2 <= 0) /* short-circuit everything: from2 == 0 */
    return E_DEC_DIV_ZERO;

  /*
    Remove the remanining zeroes . For ex: for 0.000000000001
    the above while loop removes 9 zeroes and the result will have 0.0001
    these remaining zeroes are removed here
   */
  prec2 -= count_leading_zeroes((prec2 - 1) % DIG_PER_DEC1, *buf2);
  assert(prec2 > 0);

  /*
   Do the same for the first number. Remove the leading zeroes.
   Check if the number is actually 0. Then remove the remaining zeroes.
   */

  i = ((prec1 - 1) % DIG_PER_DEC1) + 1;
  while (prec1 > 0 && *buf1 == 0) {
    prec1 -= i;
    i = DIG_PER_DEC1;
    buf1++;
  }
  if (prec1 <= 0) { /* short-circuit everything: from1 == 0 */
    decimal_make_zero(to);
    return E_DEC_OK;
  }
  prec1 -= count_leading_zeroes((prec1 - 1) % DIG_PER_DEC1, *buf1);
  assert(prec1 > 0);

  /* let's fix scale_incr, taking into account frac1,frac2 increase */
  if ((scale_incr -= frac1 - from1->frac + frac2 - from2->frac) < 0)
    scale_incr = 0;

  /* Calculate the integer digits in final result */
  dintg = (prec1 - frac1) - (prec2 - frac2) + (*buf1 >= *buf2);
  if (dintg < 0) {
    dintg /= DIG_PER_DEC1;
    intg0 = 0;
  } else
    intg0 = ROUND_UP(dintg);
  if (mod) {
    /* we're calculating N1 % N2.
       The result will have
         frac=max(frac1, frac2), as for subtraction
         intg=intg2
    */
    to->sign = from1->sign;
    to->frac = std::max(from1->frac, from2->frac);
    frac0 = 0;
  } else {
    /*
      we're calculating N1/N2. N1 is in the buf1, has prec1 digits
      N2 is in the buf2, has prec2 digits. Scales are frac1 and
      frac2 accordingly.
      Thus, the result will have
         frac = ROUND_UP(frac1+frac2+scale_incr)
      and
         intg = (prec1-frac1) - (prec2-frac2) + 1
         prec = intg+frac
    */
    frac0 = ROUND_UP(frac1 + frac2 + scale_incr);
    FIX_INTG_FRAC_ERROR(to->len, intg0, frac0, error);
    to->sign = from1->sign != from2->sign;
    to->intg = intg0 * DIG_PER_DEC1;
    to->frac = frac0 * DIG_PER_DEC1;
  }
  buf0 = to->buf;
  stop0 = buf0 + intg0 + frac0;
  if (likely(div_mod))
    while (dintg++ < 0 && buf0 < &to->buf[to->len]) {
      *buf0++ = 0;
    }

  len1 = (i = ROUND_UP(prec1)) + ROUND_UP(2 * frac2 + scale_incr + 1) + 1;
  len1 = std::max(len1, 3);
  if (!(tmp1 = (dec1 *)my_alloca(len1 * sizeof(dec1)))) return E_DEC_OOM;
  memcpy(tmp1, buf1, i * sizeof(dec1));
  memset(tmp1 + i, 0, (len1 - i) * sizeof(dec1));

  start1 = tmp1;
  stop1 = start1 + len1;
  start2 = buf2;
  stop2 = buf2 + ROUND_UP(prec2) - 1;

  /* removing end zeroes */
  while (*stop2 == 0 && stop2 >= start2) stop2--;
  len2 = (int)(stop2++ - start2);

  /*
    calculating norm2 (normalized *start2) - we need *start2 to be large
    (at least > DIG_BASE/2), but unlike Knuth's Alg. D we don't want to
    normalize input numbers (as we don't make a copy of the divisor).
    Thus we normalize first dec1 of buf2 only, and we'll normalize *start1
    on the fly for the purpose of guesstimation only.
    It's also faster, as we're saving on normalization of buf2
  */
  norm_factor = DIG_BASE / (*start2 + 1);
  norm2 = (dec1)(norm_factor * start2[0]);
  if (likely(len2 > 0)) norm2 += (dec1)(norm_factor * start2[1] / DIG_BASE);

  if (*start1 < *start2)
    dcarry = *start1++;
  else
    dcarry = 0;

  /* main loop */
  for (; buf0 < stop0; buf0++) {
    /* short-circuit, if possible */
    if (unlikely(dcarry == 0 && *start1 < *start2))
      guess = 0;
    else {
      /* D3: make a guess */
      x = start1[0] + ((dec2)dcarry) * DIG_BASE;
      y = start1[1];
      guess = (norm_factor * x + norm_factor * y / DIG_BASE) / norm2;
      if (unlikely(guess >= DIG_BASE)) guess = DIG_BASE - 1;
      if (likely(len2 > 0)) {
        /* hmm, this is a suspicious trick - I removed normalization here */
        if (start2[1] * guess > (x - guess * start2[0]) * DIG_BASE + y) guess--;
        if (unlikely(start2[1] * guess >
                     (x - guess * start2[0]) * DIG_BASE + y))
          guess--;
        assert(start2[1] * guess <= (x - guess * start2[0]) * DIG_BASE + y);
      }

      /* D4: multiply and subtract */
      buf2 = stop2;
      buf1 = start1 + len2;
      assert(buf1 < stop1);
      for (carry = 0; buf2 > start2; buf1--) {
        dec1 hi, lo;
        x = guess * (*--buf2);
        hi = (dec1)(x / DIG_BASE);
        lo = (dec1)(x - ((dec2)hi) * DIG_BASE);
        SUB2(*buf1, *buf1, lo, carry);
        carry += hi;
      }
      carry = dcarry < carry;

      /* D5: check the remainder */
      if (unlikely(carry)) {
        /* D6: correct the guess */
        guess--;
        buf2 = stop2;
        buf1 = start1 + len2;
        for (carry = 0; buf2 > start2; buf1--) {
          ADD(*buf1, *buf1, *--buf2, carry);
        }
      }
    }
    if (likely(div_mod)) {
      assert(buf0 < to->buf + to->len);
      *buf0 = (dec1)guess;
    }
    dcarry = *start1;
    start1++;
  }
  if (mod) {
    /*
      now the result is in tmp1, it has
      intg=prec1-frac1  if there were no leading zeroes.
                        If leading zeroes were present, they have been removed
                        earlier. We need to now add them back to the result.
      frac=max(frac1, frac2)=to->frac
     */
    if (dcarry) *--start1 = dcarry;
    buf0 = to->buf;
    /* Calculate the final result's integer digits */
    dintg = (prec1 - frac1) - ((start1 - tmp1) * DIG_PER_DEC1);
    if (dintg < 0) {
      /* If leading zeroes in the fractional part were earlier stripped */
      intg0 = dintg / DIG_PER_DEC1;
    } else
      intg0 = ROUND_UP(dintg);
    frac0 = ROUND_UP(to->frac);
    error = E_DEC_OK;
    if (unlikely(frac0 == 0 && intg0 == 0)) {
      decimal_make_zero(to);
      goto done;
    }
    if (intg0 <= 0) {
      /* Add back the leading zeroes that were earlier stripped */
      if (unlikely(-intg0 >= to->len)) {
        decimal_make_zero(to);
        error = E_DEC_TRUNCATED;
        goto done;
      }
      stop1 = start1 + frac0 + intg0;
      frac0 += intg0;
      to->intg = 0;
      while (intg0++ < 0) *buf0++ = 0;
    } else {
      if (unlikely(intg0 > to->len)) {
        frac0 = 0;
        intg0 = to->len;
        error = E_DEC_OVERFLOW;
        goto done;
      }
      assert(intg0 <= ROUND_UP(from2->intg));
      stop1 = start1 + frac0 + intg0;
      to->intg = std::min(intg0 * DIG_PER_DEC1, from2->intg);
    }
    if (unlikely(intg0 + frac0 > to->len)) {
      stop1 -= frac0 + intg0 - to->len;
      frac0 = to->len - intg0;
      to->frac = frac0 * DIG_PER_DEC1;
      error = E_DEC_TRUNCATED;
    }
    assert(buf0 + (stop1 - start1) <= to->buf + to->len);
    while (start1 < stop1) *buf0++ = *start1++;
  }
done:
  if (decimal_is_zero(to)) {
    // Return "0." rather than "0.000000"
    decimal_make_zero(to);
  } else {
    tmp1 = remove_leading_zeroes(to, &to->intg);
    if (to->buf != tmp1)
      memmove(to->buf, tmp1,
              (ROUND_UP(to->intg) + ROUND_UP(to->frac)) * sizeof(dec1));
  }
  assert(to->intg + to->frac > 0);
  return error;
}