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