Annotation of OpenXM_contrib/gmp/mpz/dmincl.c, Revision 1.1
1.1 ! maekawa 1: /* dmincl.c -- include file for tdiv_qr.c, tdiv_r.c.
! 2:
! 3: Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
! 4:
! 5: This file is part of the GNU MP Library.
! 6:
! 7: The GNU MP Library is free software; you can redistribute it and/or modify
! 8: it under the terms of the GNU Library General Public License as published by
! 9: the Free Software Foundation; either version 2 of the License, or (at your
! 10: option) any later version.
! 11:
! 12: The GNU MP Library is distributed in the hope that it will be useful, but
! 13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 15: License for more details.
! 16:
! 17: You should have received a copy of the GNU Library General Public License
! 18: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 20: MA 02111-1307, USA. */
! 21:
! 22: /* If den == quot, den needs temporary storage.
! 23: If den == rem, den needs temporary storage.
! 24: If num == quot, num needs temporary storage.
! 25: If den has temporary storage, it can be normalized while being copied,
! 26: i.e no extra storage should be allocated. */
! 27:
! 28: /* This is the function body of mdiv, mpz_divmod, and mpz_mod.
! 29:
! 30: If COMPUTE_QUOTIENT is defined, the quotient is put in the MP_INT
! 31: object quot, otherwise that variable is not referenced at all.
! 32:
! 33: The remainder is always computed, and the result is put in the MP_INT
! 34: object rem. */
! 35:
! 36: {
! 37: mp_ptr np, dp;
! 38: mp_ptr qp, rp;
! 39: mp_size_t nsize = num->_mp_size;
! 40: mp_size_t dsize = den->_mp_size;
! 41: mp_size_t qsize, rsize;
! 42: mp_size_t sign_remainder = nsize;
! 43: #ifdef COMPUTE_QUOTIENT
! 44: mp_size_t sign_quotient = nsize ^ dsize;
! 45: #endif
! 46: unsigned normalization_steps;
! 47: mp_limb_t q_limb;
! 48: TMP_DECL (marker);
! 49:
! 50: nsize = ABS (nsize);
! 51: dsize = ABS (dsize);
! 52:
! 53: /* Ensure space is enough for quotient and remainder. */
! 54:
! 55: /* We need space for an extra limb in the remainder, because it's
! 56: up-shifted (normalized) below. */
! 57: rsize = nsize + 1;
! 58: if (rem->_mp_alloc < rsize)
! 59: _mpz_realloc (rem, rsize);
! 60:
! 61: qsize = rsize - dsize; /* qsize cannot be bigger than this. */
! 62: if (qsize <= 0)
! 63: {
! 64: if (num != rem)
! 65: {
! 66: rem->_mp_size = num->_mp_size;
! 67: MPN_COPY (rem->_mp_d, num->_mp_d, nsize);
! 68: }
! 69: #ifdef COMPUTE_QUOTIENT
! 70: /* This needs to follow the assignment to rem, in case the
! 71: numerator and quotient are the same. */
! 72: quot->_mp_size = 0;
! 73: #endif
! 74: return;
! 75: }
! 76:
! 77: #ifdef COMPUTE_QUOTIENT
! 78: if (quot->_mp_alloc < qsize)
! 79: _mpz_realloc (quot, qsize);
! 80: #endif
! 81:
! 82: /* Read pointers here, when reallocation is finished. */
! 83: np = num->_mp_d;
! 84: dp = den->_mp_d;
! 85: rp = rem->_mp_d;
! 86:
! 87: /* Optimize division by a single-limb divisor. */
! 88: if (dsize == 1)
! 89: {
! 90: mp_limb_t rlimb;
! 91: #ifdef COMPUTE_QUOTIENT
! 92: qp = quot->_mp_d;
! 93: rlimb = mpn_divmod_1 (qp, np, nsize, dp[0]);
! 94: qsize -= qp[qsize - 1] == 0;
! 95: quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize;
! 96: #else
! 97: rlimb = mpn_mod_1 (np, nsize, dp[0]);
! 98: #endif
! 99: rp[0] = rlimb;
! 100: rsize = rlimb != 0;
! 101: rem->_mp_size = sign_remainder >= 0 ? rsize : -rsize;
! 102: return;
! 103: }
! 104:
! 105: TMP_MARK (marker);
! 106:
! 107: #ifdef COMPUTE_QUOTIENT
! 108: qp = quot->_mp_d;
! 109:
! 110: /* Make sure QP and NP point to different objects. Otherwise the
! 111: numerator would be gradually overwritten by the quotient limbs. */
! 112: if (qp == np)
! 113: {
! 114: /* Copy NP object to temporary space. */
! 115: np = (mp_ptr) TMP_ALLOC (nsize * BYTES_PER_MP_LIMB);
! 116: MPN_COPY (np, qp, nsize);
! 117: }
! 118:
! 119: #else
! 120: /* Put quotient at top of remainder. */
! 121: qp = rp + dsize;
! 122: #endif
! 123:
! 124: count_leading_zeros (normalization_steps, dp[dsize - 1]);
! 125:
! 126: /* Normalize the denominator, i.e. make its most significant bit set by
! 127: shifting it NORMALIZATION_STEPS bits to the left. Also shift the
! 128: numerator the same number of steps (to keep the quotient the same!). */
! 129: if (normalization_steps != 0)
! 130: {
! 131: mp_ptr tp;
! 132: mp_limb_t nlimb;
! 133:
! 134: /* Shift up the denominator setting the most significant bit of
! 135: the most significant word. Use temporary storage not to clobber
! 136: the original contents of the denominator. */
! 137: tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
! 138: mpn_lshift (tp, dp, dsize, normalization_steps);
! 139: dp = tp;
! 140:
! 141: /* Shift up the numerator, possibly introducing a new most
! 142: significant word. Move the shifted numerator in the remainder
! 143: meanwhile. */
! 144: nlimb = mpn_lshift (rp, np, nsize, normalization_steps);
! 145: if (nlimb != 0)
! 146: {
! 147: rp[nsize] = nlimb;
! 148: rsize = nsize + 1;
! 149: }
! 150: else
! 151: rsize = nsize;
! 152: }
! 153: else
! 154: {
! 155: /* The denominator is already normalized, as required. Copy it to
! 156: temporary space if it overlaps with the quotient or remainder. */
! 157: #ifdef COMPUTE_QUOTIENT
! 158: if (dp == rp || dp == qp)
! 159: #else
! 160: if (dp == rp)
! 161: #endif
! 162: {
! 163: mp_ptr tp;
! 164:
! 165: tp = (mp_ptr) TMP_ALLOC (dsize * BYTES_PER_MP_LIMB);
! 166: MPN_COPY (tp, dp, dsize);
! 167: dp = tp;
! 168: }
! 169:
! 170: /* Move the numerator to the remainder. */
! 171: if (rp != np)
! 172: MPN_COPY (rp, np, nsize);
! 173:
! 174: rsize = nsize;
! 175: }
! 176:
! 177: q_limb = mpn_divmod (qp, rp, rsize, dp, dsize);
! 178:
! 179: #ifdef COMPUTE_QUOTIENT
! 180: qsize = rsize - dsize;
! 181: if (q_limb)
! 182: {
! 183: qp[qsize] = q_limb;
! 184: qsize += 1;
! 185: }
! 186:
! 187: quot->_mp_size = sign_quotient >= 0 ? qsize : -qsize;
! 188: #endif
! 189:
! 190: rsize = dsize;
! 191: MPN_NORMALIZE (rp, rsize);
! 192:
! 193: if (normalization_steps != 0 && rsize != 0)
! 194: {
! 195: mpn_rshift (rp, rp, rsize, normalization_steps);
! 196: rsize -= rp[rsize - 1] == 0;
! 197: }
! 198:
! 199: rem->_mp_size = sign_remainder >= 0 ? rsize : -rsize;
! 200: TMP_FREE (marker);
! 201: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to dmincl.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz