Annotation of OpenXM_contrib/gmp/mpn/generic/divmod_1.c, Revision 1.1
1.1 ! maekawa 1: /* mpn_divmod_1(quot_ptr, dividend_ptr, dividend_size, divisor_limb) --
! 2: Divide (DIVIDEND_PTR,,DIVIDEND_SIZE) by DIVISOR_LIMB.
! 3: Write DIVIDEND_SIZE limbs of quotient at QUOT_PTR.
! 4: Return the single-limb remainder.
! 5: There are no constraints on the value of the divisor.
! 6:
! 7: QUOT_PTR and DIVIDEND_PTR might point to the same limb.
! 8:
! 9: Copyright (C) 1991, 1993, 1994, 1996 Free Software Foundation, Inc.
! 10:
! 11: This file is part of the GNU MP Library.
! 12:
! 13: The GNU MP Library is free software; you can redistribute it and/or modify
! 14: it under the terms of the GNU Library General Public License as published by
! 15: the Free Software Foundation; either version 2 of the License, or (at your
! 16: option) any later version.
! 17:
! 18: The GNU MP Library is distributed in the hope that it will be useful, but
! 19: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 20: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
! 21: License for more details.
! 22:
! 23: You should have received a copy of the GNU Library General Public License
! 24: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 25: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 26: MA 02111-1307, USA. */
! 27:
! 28: #include "gmp.h"
! 29: #include "gmp-impl.h"
! 30: #include "longlong.h"
! 31:
! 32: #ifndef UMUL_TIME
! 33: #define UMUL_TIME 1
! 34: #endif
! 35:
! 36: #ifndef UDIV_TIME
! 37: #define UDIV_TIME UMUL_TIME
! 38: #endif
! 39:
! 40: /* FIXME: We should be using invert_limb (or invert_normalized_limb)
! 41: here (not udiv_qrnnd). */
! 42:
! 43: mp_limb_t
! 44: #if __STDC__
! 45: mpn_divmod_1 (mp_ptr quot_ptr,
! 46: mp_srcptr dividend_ptr, mp_size_t dividend_size,
! 47: mp_limb_t divisor_limb)
! 48: #else
! 49: mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
! 50: mp_ptr quot_ptr;
! 51: mp_srcptr dividend_ptr;
! 52: mp_size_t dividend_size;
! 53: mp_limb_t divisor_limb;
! 54: #endif
! 55: {
! 56: mp_size_t i;
! 57: mp_limb_t n1, n0, r;
! 58: int dummy;
! 59:
! 60: /* ??? Should this be handled at all? Rely on callers? */
! 61: if (dividend_size == 0)
! 62: return 0;
! 63:
! 64: /* If multiplication is much faster than division, and the
! 65: dividend is large, pre-invert the divisor, and use
! 66: only multiplications in the inner loop. */
! 67:
! 68: /* This test should be read:
! 69: Does it ever help to use udiv_qrnnd_preinv?
! 70: && Does what we save compensate for the inversion overhead? */
! 71: if (UDIV_TIME > (2 * UMUL_TIME + 6)
! 72: && (UDIV_TIME - (2 * UMUL_TIME + 6)) * dividend_size > UDIV_TIME)
! 73: {
! 74: int normalization_steps;
! 75:
! 76: count_leading_zeros (normalization_steps, divisor_limb);
! 77: if (normalization_steps != 0)
! 78: {
! 79: mp_limb_t divisor_limb_inverted;
! 80:
! 81: divisor_limb <<= normalization_steps;
! 82:
! 83: /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
! 84: result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
! 85: most significant bit (with weight 2**N) implicit. */
! 86:
! 87: /* Special case for DIVISOR_LIMB == 100...000. */
! 88: if (divisor_limb << 1 == 0)
! 89: divisor_limb_inverted = ~(mp_limb_t) 0;
! 90: else
! 91: udiv_qrnnd (divisor_limb_inverted, dummy,
! 92: -divisor_limb, 0, divisor_limb);
! 93:
! 94: n1 = dividend_ptr[dividend_size - 1];
! 95: r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
! 96:
! 97: /* Possible optimization:
! 98: if (r == 0
! 99: && divisor_limb > ((n1 << normalization_steps)
! 100: | (dividend_ptr[dividend_size - 2] >> ...)))
! 101: ...one division less... */
! 102:
! 103: for (i = dividend_size - 2; i >= 0; i--)
! 104: {
! 105: n0 = dividend_ptr[i];
! 106: udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
! 107: ((n1 << normalization_steps)
! 108: | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
! 109: divisor_limb, divisor_limb_inverted);
! 110: n1 = n0;
! 111: }
! 112: udiv_qrnnd_preinv (quot_ptr[0], r, r,
! 113: n1 << normalization_steps,
! 114: divisor_limb, divisor_limb_inverted);
! 115: return r >> normalization_steps;
! 116: }
! 117: else
! 118: {
! 119: mp_limb_t divisor_limb_inverted;
! 120:
! 121: /* Compute (2**2N - 2**N * DIVISOR_LIMB) / DIVISOR_LIMB. The
! 122: result is a (N+1)-bit approximation to 1/DIVISOR_LIMB, with the
! 123: most significant bit (with weight 2**N) implicit. */
! 124:
! 125: /* Special case for DIVISOR_LIMB == 100...000. */
! 126: if (divisor_limb << 1 == 0)
! 127: divisor_limb_inverted = ~(mp_limb_t) 0;
! 128: else
! 129: udiv_qrnnd (divisor_limb_inverted, dummy,
! 130: -divisor_limb, 0, divisor_limb);
! 131:
! 132: i = dividend_size - 1;
! 133: r = dividend_ptr[i];
! 134:
! 135: if (r >= divisor_limb)
! 136: r = 0;
! 137: else
! 138: {
! 139: quot_ptr[i] = 0;
! 140: i--;
! 141: }
! 142:
! 143: for (; i >= 0; i--)
! 144: {
! 145: n0 = dividend_ptr[i];
! 146: udiv_qrnnd_preinv (quot_ptr[i], r, r,
! 147: n0, divisor_limb, divisor_limb_inverted);
! 148: }
! 149: return r;
! 150: }
! 151: }
! 152: else
! 153: {
! 154: if (UDIV_NEEDS_NORMALIZATION)
! 155: {
! 156: int normalization_steps;
! 157:
! 158: count_leading_zeros (normalization_steps, divisor_limb);
! 159: if (normalization_steps != 0)
! 160: {
! 161: divisor_limb <<= normalization_steps;
! 162:
! 163: n1 = dividend_ptr[dividend_size - 1];
! 164: r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
! 165:
! 166: /* Possible optimization:
! 167: if (r == 0
! 168: && divisor_limb > ((n1 << normalization_steps)
! 169: | (dividend_ptr[dividend_size - 2] >> ...)))
! 170: ...one division less... */
! 171:
! 172: for (i = dividend_size - 2; i >= 0; i--)
! 173: {
! 174: n0 = dividend_ptr[i];
! 175: udiv_qrnnd (quot_ptr[i + 1], r, r,
! 176: ((n1 << normalization_steps)
! 177: | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
! 178: divisor_limb);
! 179: n1 = n0;
! 180: }
! 181: udiv_qrnnd (quot_ptr[0], r, r,
! 182: n1 << normalization_steps,
! 183: divisor_limb);
! 184: return r >> normalization_steps;
! 185: }
! 186: }
! 187: /* No normalization needed, either because udiv_qrnnd doesn't require
! 188: it, or because DIVISOR_LIMB is already normalized. */
! 189:
! 190: i = dividend_size - 1;
! 191: r = dividend_ptr[i];
! 192:
! 193: if (r >= divisor_limb)
! 194: r = 0;
! 195: else
! 196: {
! 197: quot_ptr[i] = 0;
! 198: i--;
! 199: }
! 200:
! 201: for (; i >= 0; i--)
! 202: {
! 203: n0 = dividend_ptr[i];
! 204: udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);
! 205: }
! 206: return r;
! 207: }
! 208: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to divmod_1.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpn / generic