Commit 5d0d2e2b authored by Stefan Krah's avatar Stefan Krah

Explain the strategy to avoid huge alignment shifts in _mpd_qadd() in detail.

parent ed4b21ff
......@@ -3103,14 +3103,43 @@ _mpd_qaddsub(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b,
_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;
......@@ -3118,7 +3147,7 @@ _mpd_qaddsub(mpd_t *result, const mpd_t *a, const mpd_t *b, uint8_t sign_b,
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);
......@@ -3521,7 +3550,7 @@ mpd_qdiv(mpd_t *q, const mpd_t *a, const mpd_t *b,
/* 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;
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment