_mpd_ptrswap(&big, &small);
swap++;
}
+ /* align the coefficients */
if (!mpd_iszerocoeff(big)) {
- /* Test for adjexp(small) + big->digits < adjexp(big), if big-digits > prec
- * Test for adjexp(small) + prec + 1 < adjexp(big), if big-digits <= prec
- * If true, the magnitudes of the numbers are so far apart that one can as
- * well add or subtract 1*10**big->exp. */
exp = big->exp - 1;
exp += (big->digits > ctx->prec) ? 0 : big->digits-ctx->prec-1;
if (mpd_adjexp(small) < exp) {
+ /*
+ * Avoid huge shifts by substituting a value for small that is
+ * guaranteed to produce the same results.
+ *
+ * adjexp(small) < exp if and only if:
+ *
+ * bdigits <= prec AND
+ * bdigits+shift >= prec+2+sdigits AND
+ * exp = bexp+bdigits-prec-2
+ *
+ * 1234567000000000 -> bdigits + shift
+ * ----------XX1234 -> sdigits
+ * ----------X1 -> tiny-digits
+ * |- prec -|
+ *
+ * OR
+ *
+ * bdigits > prec AND
+ * shift > sdigits AND
+ * exp = bexp-1
+ *
+ * 1234567892100000 -> bdigits + shift
+ * ----------XX1234 -> sdigits
+ * ----------X1 -> tiny-digits
+ * |- prec -|
+ *
+ * If tiny is zero, adding or subtracting is a no-op.
+ * Otherwise, adding tiny generates a non-zero digit either
+ * below the rounding digit or the least significant digit
+ * of big. When subtracting, tiny is in the same position as
+ * the carry that would be generated by subtracting sdigits.
+ */
mpd_copy_flags(&tiny, small);
tiny.exp = exp;
tiny.digits = 1;
tiny.data[0] = mpd_iszerocoeff(small) ? 0 : 1;
small = &tiny;
}
- /* this cannot wrap: the difference is positive and <= maxprec+1 */
+ /* This cannot wrap: the difference is positive and <= maxprec */
shift = big->exp - small->exp;
if (!mpd_qshiftl(&big_aligned, big, shift, status)) {
mpd_seterror(result, MPD_Malloc_error, status);
/* Internal function. */
static void
_mpd_qdivmod(mpd_t *q, mpd_t *r, const mpd_t *a, const mpd_t *b,
- const mpd_context_t *ctx, uint32_t *status)
+ const mpd_context_t *ctx, uint32_t *status)
{
MPD_NEW_STATIC(aligned,0,0,0,0);
mpd_ssize_t qsize, rsize;