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.
7 QUOT_PTR and DIVIDEND_PTR might point to the same limb.
9 Copyright (C) 1991, 1993 Free Software Foundation, Inc.
11 This file is part of the GNU MP Library.
13 The GNU MP Library is free software; you can redistribute it and/or modify
14 it under the terms of the GNU General Public License as published by
15 the Free Software Foundation; either version 2, or (at your option)
18 The GNU MP Library is distributed in the hope that it will be useful,
19 but WITHOUT ANY WARRANTY; without even the implied warranty of
20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 GNU General Public License for more details.
23 You should have received a copy of the GNU General Public License
24 along with the GNU MP Library; see the file COPYING. If not, write to
25 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
36 #define UDIV_TIME UMUL_TIME
39 #if UDIV_TIME > 2 * UMUL_TIME
40 #undef UDIV_NEEDS_NORMALIZATION
41 #define UDIV_NEEDS_NORMALIZATION 1
44 #define udiv_qrnnd_preinv(q, r, nh, nl, d, di) \
46 unsigned long int _q, _ql, _r; \
47 unsigned long int _xh, _xl; \
48 umul_ppmm (_q, _ql, (nh), (di)); \
49 _q += (nh); /* DI is 2**BITS_PER_MP_LIMB too small. */\
50 umul_ppmm (_xh, _xl, _q, (d)); \
51 sub_ddmmss (_xh, _r, (nh), (nl), _xh, _xl); \
54 sub_ddmmss (_xh, _r, _xh, _r, 0, (d)); \
58 sub_ddmmss (_xh, _r, _xh, _r, 0, (d)); \
73 mpn_divmod_1 (mp_ptr quot_ptr,
74 mp_srcptr dividend_ptr, mp_size dividend_size,
75 unsigned long int divisor_limb)
77 mpn_divmod_1 (quot_ptr, dividend_ptr, dividend_size, divisor_limb)
79 mp_srcptr dividend_ptr;
80 mp_size dividend_size;
81 unsigned long int divisor_limb;
87 /* Botch: Should this be handled at all? Rely on callers? */
88 if (dividend_size == 0)
91 if (UDIV_NEEDS_NORMALIZATION)
93 int normalization_steps;
95 count_leading_zeros (normalization_steps, divisor_limb);
96 if (normalization_steps != 0)
98 divisor_limb <<= normalization_steps;
100 n1 = dividend_ptr[dividend_size - 1];
101 r = n1 >> (BITS_PER_MP_LIMB - normalization_steps);
103 /* Possible optimization:
105 && divisor_limb > ((n1 << normalization_steps)
106 | (dividend_ptr[dividend_size - 2] >> ...)))
107 ...one division less...
108 [Don't forget to zero most sign. quotient limb!] */
110 /* If multiplication is much faster than division, and the
111 dividend is large, pre-invert the divisor, and use
112 only multiplications in the inner loop. */
113 if (UDIV_TIME > 2 * UMUL_TIME && dividend_size >= 4)
115 mp_limb divisor_limb_inverted;
117 /* Compute (2**64 - 2**32 * DIVISOR_LIMB) / DIVISOR_LIMB.
118 The result is an 33-bit approximation to 1/DIVISOR_LIMB,
119 with the most significant bit (weight 2**32) implicit. */
121 /* Special case for DIVISOR_LIMB == 100...000. */
122 if (divisor_limb << 1 == 0)
123 divisor_limb_inverted = ~0;
125 udiv_qrnnd (divisor_limb_inverted, dummy,
126 -divisor_limb, 0, divisor_limb);
128 for (i = dividend_size - 2; i >= 0; i--)
130 n0 = dividend_ptr[i];
131 udiv_qrnnd_preinv (quot_ptr[i + 1], r, r,
132 ((n1 << normalization_steps)
133 | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
134 divisor_limb, divisor_limb_inverted);
137 udiv_qrnnd_preinv (quot_ptr[0], r, r,
138 n1 << normalization_steps,
139 divisor_limb, divisor_limb_inverted);
140 return r >> normalization_steps;
144 for (i = dividend_size - 2; i >= 0; i--)
146 n0 = dividend_ptr[i];
147 udiv_qrnnd (quot_ptr[i + 1], r, r,
148 ((n1 << normalization_steps)
149 | (n0 >> (BITS_PER_MP_LIMB - normalization_steps))),
153 udiv_qrnnd (quot_ptr[0], r, r,
154 n1 << normalization_steps,
156 return r >> normalization_steps;
161 /* No normalization needed, either because udiv_qrnnd doesn't require
162 it, or because DIVISOR_LIMB is already normalized. */
164 i = dividend_size - 1;
167 if (r >= divisor_limb)
173 /* Callers expect the quotient to be DIVIDEND_SIZE limbs. Store
174 a leading zero to make that expectation come true. */
181 n0 = dividend_ptr[i];
182 udiv_qrnnd (quot_ptr[i], r, r, n0, divisor_limb);