diff options
author | Kevin Ryde <user42@zip.com.au> | 2004-03-13 23:58:07 +0100 |
---|---|---|
committer | Kevin Ryde <user42@zip.com.au> | 2004-03-13 23:58:07 +0100 |
commit | 726e41dd8e9b50f018ce66f21eaeb8513e0b8d58 (patch) | |
tree | 2accd97efabadaa4b5fe5696cd4d76cdf24e4253 /mpf | |
parent | fe360ada1150d595b3adc0b38ba9685833e89f2a (diff) | |
download | gmp-726e41dd8e9b50f018ce66f21eaeb8513e0b8d58.tar.gz |
* mpf/div.c: Use mpn_tdiv_qr. Use just one TMP_ALLOC. Use full
divisor, since truncating can lose accuracy.
Diffstat (limited to 'mpf')
-rw-r--r-- | mpf/div.c | 144 |
1 files changed, 76 insertions, 68 deletions
@@ -1,6 +1,7 @@ /* mpf_div -- Divide two floats. -Copyright 1993, 1994, 1996, 2000, 2001, 2002 Free Software Foundation, Inc. +Copyright 1993, 1994, 1996, 2000, 2001, 2002, 2004 Free Software Foundation, +Inc. This file is part of the GNU MP Library. @@ -19,21 +20,49 @@ along with the GNU MP Library; see the file COPYING.LIB. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ +#include <stdio.h> /* for NULL */ #include "gmp.h" #include "gmp-impl.h" #include "longlong.h" + +/* Not done: + + No attempt is made to identify an overlap u==v. The result will be + correct (1.0), but a full actual division is done whereas of course + x/x==1 needs no work. Such a call is not a sensible thing to make, and + it's left to an application to notice and optimize if it might arise + somehow through pointer aliasing or whatever. + + Enhancements: + + The high quotient limb is non-zero when high{up,vsize} >= {vp,vsize}. We + could make that comparison and use qsize==prec instead of qsize==prec+1, + to save one limb in the division. + + If r==u but the size is enough bigger than prec that there won't be an + overlap between quotient and dividend in mpn_tdiv_qr, then we can avoid + copying up,usize. This would only arise from a prec reduced with + mpf_set_prec_raw and will be pretty unusual, but might be worthwhile if + it could be worked into the copy_u decision cleanly. + + Future: + + If/when mpn_tdiv_qr supports its qxn parameter we can use that instead of + padding u 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_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) { mp_srcptr up, vp; - mp_ptr rp, tp, rtp; - mp_size_t usize, vsize; - mp_size_t rsize, tsize; - mp_size_t sign_quotient; - mp_size_t prec; - mp_limb_t q_limb; + mp_ptr rp, remp, tp, new_vp; + mp_size_t usize, vsize, rsize, prospective_rsize, tsize, zeros, copy_v_size; + mp_size_t sign_quotient, prec, high_zero, chop; mp_exp_t rexp; + int copy_u; TMP_DECL (marker); usize = u->_mp_size; @@ -54,86 +83,65 @@ mpf_div (mpf_ptr r, mpf_srcptr u, mpf_srcptr v) } TMP_MARK (marker); - rexp = u->_mp_exp - v->_mp_exp; + rexp = u->_mp_exp - v->_mp_exp + 1; rp = r->_mp_d; up = u->_mp_d; vp = v->_mp_d; - if (vsize > prec) - { - vp += vsize - prec; - vsize = prec; - } + prospective_rsize = usize - vsize + 1; /* quot from using given u,v sizes */ + rsize = prec + 1; /* desired quot */ - tsize = vsize + prec; - tp = (mp_ptr) TMP_ALLOC ((tsize + 1) * BYTES_PER_MP_LIMB); + zeros = rsize - prospective_rsize; /* padding u to give rsize */ + copy_u = (zeros > 0 || rp == up); /* copy u if overlap or padding */ - if (usize > tsize) + chop = MAX (-zeros, 0); /* negative zeros means shorten u */ + up += chop; + usize -= chop; + zeros += chop; /* now zeros >= 0 */ + + tsize = usize + zeros; /* size for possible copy of u */ + + if (WANT_TMP_DEBUG) { - up += usize - tsize; - usize = tsize; - rtp = tp; + /* separate blocks, for malloc debugging */ + remp = TMP_ALLOC_LIMBS (vsize); + tp = (copy_u ? TMP_ALLOC_LIMBS (tsize) : NULL); + new_vp = (rp == vp ? TMP_ALLOC_LIMBS (vsize) : NULL); } else { - MPN_ZERO (tp, tsize - usize); - rtp = tp + (tsize - usize); + /* one block with conditionalized size, for efficiency */ + copy_v_size = (rp == vp ? vsize : 0); + remp = TMP_ALLOC_LIMBS (vsize + copy_v_size + (copy_u ? tsize : 0)); + new_vp = remp + vsize; + tp = new_vp + copy_v_size; } - - /* Normalize the divisor and the dividend. */ - if ((vp[vsize-1] & GMP_NUMB_HIGHBIT) == 0) - { - mp_ptr tmp; - mp_limb_t nlimb; - unsigned normalization_steps; - - count_leading_zeros (normalization_steps, vp[vsize - 1]); - normalization_steps -= GMP_NAIL_BITS; - - /* Shift up the divisor setting the most significant bit of - the most significant limb. Use temporary storage not to clobber - the original contents of the divisor. */ - tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); - mpn_lshift (tmp, vp, vsize, normalization_steps); - vp = tmp; - - /* Shift up the dividend, possibly introducing a new most - significant word. Move the shifted dividend in the remainder - at the same time. */ - nlimb = mpn_lshift (rtp, up, usize, normalization_steps); - if (nlimb != 0) - { - rtp[usize] = nlimb; - tsize++; - rexp++; - } - } - else + /* copy and possibly extend u if necessary */ + if (copy_u) { - /* The divisor is already normalized, as required. - Copy it to temporary space if it overlaps with the quotient. */ - if (vp - rp <= tsize - vsize) - { - mp_ptr tmp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB); - MPN_COPY (tmp, vp, vsize); - vp = (mp_srcptr) tmp; - } - - /* Move the dividend to the remainder. */ - MPN_COPY (rtp, up, usize); + MPN_ZERO (tp, zeros); + MPN_COPY (tp+zeros, up, usize); + up = tp; + usize = tsize; } - q_limb = mpn_divmod (rp, tp, tsize, vp, vsize); - rsize = tsize - vsize; - if (q_limb) + /* ensure divisor doesn't overlap quotient */ + if (rp == vp) { - rp[rsize] = q_limb; - rsize++; - rexp++; + MPN_COPY (new_vp, vp, vsize); + vp = new_vp; } + ASSERT (usize-vsize+1 == rsize); + mpn_tdiv_qr (rp, remp, (mp_size_t) 0, up, usize, vp, vsize); + + /* strip possible zero high limb */ + high_zero = (rp[rsize-1] == 0); + rsize -= high_zero; + rexp -= high_zero; + r->_mp_size = sign_quotient >= 0 ? rsize : -rsize; r->_mp_exp = rexp; TMP_FREE (marker); |