1 /* mpn_divrem -- Divide natural numbers, producing both remainder and
2 quotient.
3
4 Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
5
6 This file is part of the GNU MP Library.
7
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 2.1 of the License, or (at your
11 option) any later version.
12
13 The GNU MP Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
17
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21 MA 02111-1307, USA. */
22
23 #include <config.h>
24 #include "gmp-impl.h"
25
26 /* Divide num (NP/NSIZE) by den (DP/DSIZE) and write
27 the NSIZE-DSIZE least significant quotient limbs at QP
28 and the DSIZE long remainder at NP. If QEXTRA_LIMBS is
29 non-zero, generate that many fraction bits and append them after the
30 other quotient limbs.
31 Return the most significant limb of the quotient, this is always 0 or 1.
32
33 Preconditions:
34 0. NSIZE >= DSIZE.
35 1. The most significant bit of the divisor must be set.
36 2. QP must either not overlap with the input operands at all, or
37 QP + DSIZE >= NP must hold true. (This means that it's
38 possible to put the quotient in the high part of NUM, right after the
39 remainder in NUM.
40 3. NSIZE >= DSIZE, even if QEXTRA_LIMBS is non-zero. */
41
42 mp_limb_t
43 #if __STDC__
44 mpn_divrem (mp_ptr qp, mp_size_t qextra_limbs,
45 mp_ptr np, mp_size_t nsize,
46 mp_srcptr dp, mp_size_t dsize)
47 #else
48 mpn_divrem (qp, qextra_limbs, np, nsize, dp, dsize)
49 mp_ptr qp;
50 mp_size_t qextra_limbs;
51 mp_ptr np;
52 mp_size_t nsize;
53 mp_srcptr dp;
54 mp_size_t dsize;
55 #endif
56 {
57 mp_limb_t most_significant_q_limb = 0;
58
59 switch (dsize)
60 {
61 case 0:
62 /* We are asked to divide by zero, so go ahead and do it! (To make
63 the compiler not remove this statement, return the value.) */
64 return 1 / dsize;
65
66 case 1:
67 {
68 mp_size_t i;
69 mp_limb_t n1;
70 mp_limb_t d;
71
72 d = dp[0];
73 n1 = np[nsize - 1];
74
75 if (n1 >= d)
76 {
77 n1 -= d;
78 most_significant_q_limb = 1;
79 }
80
81 qp += qextra_limbs;
82 for (i = nsize - 2; i >= 0; i--)
83 udiv_qrnnd (qp[i], n1, n1, np[i], d);
84 qp -= qextra_limbs;
85
86 for (i = qextra_limbs - 1; i >= 0; i--)
87 udiv_qrnnd (qp[i], n1, n1, 0, d);
88
89 np[0] = n1;
90 }
91 break;
92
93 case 2:
94 {
95 mp_size_t i;
96 mp_limb_t n1, n0, n2;
97 mp_limb_t d1, d0;
98
99 np += nsize - 2;
100 d1 = dp[1];
101 d0 = dp[0];
102 n1 = np[1];
103 n0 = np[0];
104
105 if (n1 >= d1 && (n1 > d1 || n0 >= d0))
106 {
107 sub_ddmmss (n1, n0, n1, n0, d1, d0);
108 most_significant_q_limb = 1;
109 }
110
111 for (i = qextra_limbs + nsize - 2 - 1; i >= 0; i--)
112 {
113 mp_limb_t q;
114 mp_limb_t r;
115
116 if (i >= qextra_limbs)
117 np--;
118 else
119 np[0] = 0;
120
121 if (n1 == d1)
122 {
123 /* Q should be either 111..111 or 111..110. Need special
124 treatment of this rare case as normal division would
125 give overflow. */
126 q = ~(mp_limb_t) 0;
127
128 r = n0 + d1;
129 if (r < d1) /* Carry in the addition? */
130 {
131 add_ssaaaa (n1, n0, r - d0, np[0], 0, d0);
132 qp[i] = q;
133 continue;
134 }
135 n1 = d0 - (d0 != 0);
136 n0 = -d0;
137 }
138 else
139 {
140 udiv_qrnnd (q, r, n1, n0, d1);
141 umul_ppmm (n1, n0, d0, q);
142 }
143
144 n2 = np[0];
145 q_test:
146 if (n1 > r || (n1 == r && n0 > n2))
147 {
148 /* The estimated Q was too large. */
149 q--;
150
151 sub_ddmmss (n1, n0, n1, n0, 0, d0);
152 r += d1;
153 if (r >= d1) /* If not carry, test Q again. */
154 goto q_test;
155 }
156
157 qp[i] = q;
158 sub_ddmmss (n1, n0, r, n2, n1, n0);
159 }
160 np[1] = n1;
161 np[0] = n0;
162 }
163 break;
164
165 default:
166 {
167 mp_size_t i;
168 mp_limb_t dX, d1, n0;
169
170 np += nsize - dsize;
171 dX = dp[dsize - 1];
172 d1 = dp[dsize - 2];
173 n0 = np[dsize - 1];
174
175 if (n0 >= dX)
176 {
177 if (n0 > dX || mpn_cmp (np, dp, dsize - 1) >= 0)
178 {
179 mpn_sub_n (np, np, dp, dsize);
180 n0 = np[dsize - 1];
181 most_significant_q_limb = 1;
182 }
183 }
184
185 for (i = qextra_limbs + nsize - dsize - 1; i >= 0; i--)
186 {
187 mp_limb_t q;
188 mp_limb_t n1, n2;
189 mp_limb_t cy_limb;
190
191 if (i >= qextra_limbs)
192 {
193 np--;
194 n2 = np[dsize];
195 }
196 else
197 {
198 n2 = np[dsize - 1];
199 MPN_COPY_DECR (np + 1, np, dsize);
200 np[0] = 0;
201 }
202
203 if (n0 == dX)
204 /* This might over-estimate q, but it's probably not worth
205 the extra code here to find out. */
206 q = ~(mp_limb_t) 0;
207 else
208 {
209 mp_limb_t r;
210
211 udiv_qrnnd (q, r, n0, np[dsize - 1], dX);
212 umul_ppmm (n1, n0, d1, q);
213
214 while (n1 > r || (n1 == r && n0 > np[dsize - 2]))
215 {
216 q--;
217 r += dX;
218 if (r < dX) /* I.e. "carry in previous addition?" */
219 break;
220 n1 -= n0 < d1;
221 n0 -= d1;
222 }
223 }
224
225 /* Possible optimization: We already have (q * n0) and (1 * n1)
226 after the calculation of q. Taking advantage of that, we
227 could make this loop make two iterations less. */
228
229 cy_limb = mpn_submul_1 (np, dp, dsize, q);
230
231 if (n2 != cy_limb)
232 {
233 mpn_add_n (np, np, dp, dsize);
234 q--;
235 }
236
237 qp[i] = q;
238 n0 = np[dsize - 1];
239 }
240 }
241 }
242
243 return most_significant_q_limb;
244 }