1 /* mpn_div -- Divide natural numbers, producing both remainder and
4 Copyright (C) 1991 Free Software Foundation, Inc.
6 This file is part of the GNU MP Library.
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2, or (at your option)
13 The GNU MP Library is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
18 You should have received a copy of the GNU General Public License
19 along with the GNU MP Library; see the file COPYING. If not, write to
20 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
26 /* Divide num (NUM_PTR/NUM_SIZE) by den (DEN_PTR/DEN_SIZE) and write
27 the quotient at QUOT_PTR and the remainder at NUM_PTR.
29 Return 0 or 1, depending on if the quotient size is (NSIZE - DSIZE)
30 or (NSIZE - DSIZE + 1).
33 1. The most significant bit of d must be set.
34 2. QUOT_PTR != DEN_PTR and QUOT_PTR != NUM_PTR, i.e. the quotient storage
35 area must be distinct from either input operands.
37 The exact sizes of the quotient and remainder must be determined
38 by the caller, in spite of the return value. The return value just
39 informs the caller about if the highest digit is written or not, and
40 it may very well be 0. */
42 /* THIS WILL BE IMPROVED SOON. MORE COMMENTS AND FASTER CODE. */
46 mpn_div (mp_ptr quot_ptr,
47 mp_ptr num_ptr, mp_size num_size,
48 mp_srcptr den_ptr, mp_size den_size)
50 mpn_div (quot_ptr, num_ptr, num_size, den_ptr, den_size)
58 mp_size q_is_long = 0;
63 /* We are asked to divide by zero, so go ahead and do it!
64 (To make the compiler not remove this statement, assign NUM_SIZE
66 num_size = 1 / den_size;
89 udiv_qrnnd (quot_ptr[i], n1, n1, n0, d);
102 num_ptr += num_size - 2;
117 for (i = num_size - den_size - 1; i >= 0; i--)
125 /* Q should be either 111..111 or 111..110. Need special
126 treatment of this rare case as normal division would
131 if (r < d0) /* Carry in the addition? */
135 add_ssaaaa (n0, n1, r - d1, n2, 0, d1);
144 udiv_qrnnd (q, r, n0, n1, d0);
145 umul_ppmm (n0, n1, d1, q);
150 if (n0 > r || (n0 == r && n1 > n2))
152 /* The estimated Q was too large. */
155 sub_ddmmss (n0, n1, n0, n1, 0, d1);
157 if (r >= d0) /* If not carry, test q again. */
162 sub_ddmmss (n0, n1, r, n2, n0, n1);
172 mp_limb d0 = den_ptr[den_size - 1];
173 mp_limb d1 = den_ptr[den_size - 2];
174 mp_limb n0 = num_ptr[num_size - 1];
175 int ugly_hack_flag = 0;
180 /* There's a problem with this case, which shows up later in the
181 code. q becomes 1 (and sometimes 0) the first time when
182 we've been here, and w_cy == 0 after the main do-loops below.
183 But c = num_ptr[j] reads rubbish outside the num_ptr vector!
184 Maybe I can solve this cleanly when I fix the early-end
185 optimization here in the default case. For now, I change the
186 add_back entering condition, to kludge. Leaving the stray
189 HACK: Added ugly_hack_flag to make it work. */
199 for (i = num_size - den_size - 1; i >= 0; i--)
209 /* This might over-estimate q, but it's probably not worth
210 the extra code here to find out. */
216 udiv_qrnnd (q, r, n0, num_ptr[-1], d0);
217 umul_ppmm (n1, n0, d1, q);
219 while (n1 > r || (n1 == r && n0 > num_ptr[-2]))
223 if (r < d0) /* I.e. "carry in previous addition?" */
236 umul_ppmm (n1, n0, d, q);
238 w_cy = (n0 < w_cy) + n1;
263 umul_ppmm (n1, n0, d, q);
265 w_cy = (n0 < w_cy) + n1;
299 abort (); /* We should always have a carry out! */
304 n0 = num_ptr[j] + d + 1;