diff options
author | Kevin Ryde <user42@zip.com.au> | 2004-03-03 21:51:11 +0100 |
---|---|---|
committer | Kevin Ryde <user42@zip.com.au> | 2004-03-03 21:51:11 +0100 |
commit | 8beba9fd010f794b9758e7130fe6b168c4ffe006 (patch) | |
tree | 54e2c147b8c3f65550923df32cfa7dfa2d24240e | |
parent | 9f5babcb6705b4949ca6b4c2aadc1a46f2bbb471 (diff) | |
download | gmp-8beba9fd010f794b9758e7130fe6b168c4ffe006.tar.gz |
* mpf/set_q.c: Use mpn_tdiv_qr rather than mpn_divrem, so no shifting.
Don't truncate the divisor, it can make the result inaccurate.
-rw-r--r-- | mpf/set_q.c | 147 |
1 files changed, 68 insertions, 79 deletions
diff --git a/mpf/set_q.c b/mpf/set_q.c index d13b42a63..075f3c3c0 100644 --- a/mpf/set_q.c +++ b/mpf/set_q.c @@ -23,17 +23,50 @@ MA 02111-1307, USA. */ #include "gmp-impl.h" #include "longlong.h" + +/* As usual the aim is to produce PREC(r) limbs, with the high non-zero. + The basic mpn_tdiv_qr produces a quotient of nsize-dsize+1 limbs, with + either the high or second highest limb non-zero. We arrange for + nsize-dsize+1 to equal prec+1, hence giving either prec or prec+1 result + limbs at PTR(r). + + nsize-dsize+1 == prec+1 is achieved by adjusting num(q), either dropping + low limbs if it's too big, or padding with low zeros if it's too small. + The full given den(q) is always used. + + We cannot truncate den(q), because even when it's much bigger than prec + the last limbs can still influence the final quotient. Often they don't, + but we leave optimization of that to a prospective quotient-only mpn + division. + + Not done: + + If den(q) is a power of 2 then we may end up with low zero limbs on the + result. But nothing is done about this, since it should be unlikely on + random data, and can be left to an application to call mpf_div_2exp if it + might occur with any frequency. + + Enhancements: + + The high quotient limb is non-zero when high{np,dsize} >= {dp,dsize}. We + could make that comparison and use qsize==prec instead of qsize==prec+1, + to save one limb in the division. + + Future: + + If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of + padding n with zeros in temporary space. + + If/when a quotient-only division exists it can be used here immediately. + remp is only to satisfy mpn_tdiv_qr, the remainder is not used. */ + void mpf_set_q (mpf_t r, mpq_srcptr q) { - mp_ptr np, dp; - mp_ptr rp; - mp_size_t nsize, dsize; - mp_size_t qsize, rsize; - mp_size_t sign_quotient; - mp_limb_t qlimb; - mp_ptr qp; - mp_size_t prec; + mp_srcptr np, dp; + mp_size_t prec, nsize, dsize, qsize, prospective_qsize, tsize, zeros; + mp_size_t sign_quotient, high_zero; + mp_ptr qp, tp, remp; mp_exp_t exp; TMP_DECL (marker); @@ -42,17 +75,16 @@ mpf_set_q (mpf_t r, mpq_srcptr q) nsize = SIZ (&q->_mp_num); dsize = SIZ (&q->_mp_den); - if (nsize == 0) + if (UNLIKELY (nsize == 0)) { SIZ (r) = 0; EXP (r) = 0; return; } - prec = PREC (r) + 1; - TMP_MARK (marker); + prec = PREC (r); qp = PTR (r); sign_quotient = nsize; @@ -60,83 +92,40 @@ mpf_set_q (mpf_t r, mpq_srcptr q) np = PTR (&q->_mp_num); dp = PTR (&q->_mp_den); - exp = nsize - dsize; + prospective_qsize = nsize - dsize + 1; /* q from using given n,d sizes */ + exp = prospective_qsize; /* ie. number of integer limbs */ + qsize = prec + 1; /* desired q */ - if (nsize > prec) - { - np += nsize - prec; - nsize = prec; - } - if (dsize > prec) - { - dp += dsize - prec; - dsize = prec; - } + /* how many zero limbs to pad n with to get desired qsize */ + zeros = qsize - prospective_qsize; + + /* space for copy of n, if need zeros */ + tsize = (zeros > 0 ? nsize + zeros : 1); - rsize = MAX (nsize, dsize); - rp = (mp_ptr) TMP_ALLOC ((rsize + 1) * BYTES_PER_MP_LIMB); + TMP_ALLOC_LIMBS_2 (remp, dsize, tp, tsize); - /* Normalize the denominator, i.e. make its most significant bit set by - shifting it NORMALIZATION_STEPS bits to the left. Also shift the - numerator the same number of steps (to keep the quotient the same!). */ - if ((dp[dsize - 1] & GMP_NUMB_HIGHBIT) == 0) + if (zeros > 0) { - mp_ptr tp; - mp_limb_t nlimb; - unsigned normalization_steps; - - count_leading_zeros (normalization_steps, dp[dsize - 1]); - normalization_steps -= GMP_NAIL_BITS; - - /* Shift up the denominator setting the most significant bit of - the most significant limb. Use temporary storage not to clobber - the original contents of the denominator. */ - tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB); - mpn_lshift (tp, dp, dsize, normalization_steps); - dp = tp; - - if (rsize != nsize) - { - MPN_ZERO (rp, rsize - nsize); - nlimb = mpn_lshift (rp + (rsize - nsize), - np, nsize, normalization_steps); - } - else - { - nlimb = mpn_lshift (rp, np, rsize, normalization_steps); - } - if (nlimb != 0) - { - rp[rsize] = nlimb; - exp++; - /* Don't just increase rsize, chop off rp at the low end instead. */ - if (rsize == prec) - rp++; - else - rsize++; - } + /* pad n with zeros into temporary space */ + MPN_ZERO (tp, zeros); + MPN_COPY (tp+zeros, np, nsize); + np = tp; + nsize = tsize; } else { - if (rsize != nsize) - { - MPN_ZERO (rp, rsize - nsize); - MPN_COPY (rp + (rsize - nsize), np, nsize); - } - else - { - MPN_COPY (rp, np, rsize); - } + /* shorten n to get desired qsize */ + nsize += zeros; + np -= zeros; } - qlimb = mpn_divrem (qp, prec - 1 - (rsize - dsize), rp, rsize, dp, dsize); - qsize = prec - 1; - if (qlimb) - { - qp[qsize] = qlimb; - qsize++; - exp++; - } + ASSERT (nsize-dsize+1 == qsize); + mpn_tdiv_qr (qp, remp, (mp_size_t) 0, np, nsize, dp, dsize); + + /* strip possible zero high limb */ + high_zero = (qp[qsize-1] == 0); + qsize -= high_zero; + exp -= high_zero; EXP (r) = exp; SIZ (r) = sign_quotient >= 0 ? qsize : -qsize; |