1 /* mpn_divexact(qp,np,nn,dp,dn,tp) -- Divide N = {np,nn} by D = {dp,dn} storing
2 the result in Q = {qp,nn-dn+1} expecting no remainder. Overlap allowed
3 between Q and N; all other overlap disallowed.
4
5 Contributed to the GNU project by Torbjorn Granlund.
6
7 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
8 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
9 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
10
11 Copyright 2006, 2007, 2009 Free Software Foundation, Inc.
12
13 This file is part of the GNU MP Library.
14
15 The GNU MP Library is free software; you can redistribute it and/or modify
16 it under the terms of the GNU Lesser General Public License as published by
17 the Free Software Foundation; either version 3 of the License, or (at your
18 option) any later version.
19
20 The GNU MP Library is distributed in the hope that it will be useful, but
21 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
22 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 License for more details.
24
25 You should have received a copy of the GNU Lesser General Public License
26 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
27
28
29 #include "gmp.h"
30 #include "gmp-impl.h"
31 #include "longlong.h"
32
33 #if 1
34 void
mpn_divexact(mp_ptr qp,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn)35 mpn_divexact (mp_ptr qp,
36 mp_srcptr np, mp_size_t nn,
37 mp_srcptr dp, mp_size_t dn)
38 {
39 unsigned shift;
40 mp_size_t qn;
41 mp_ptr tp, wp;
42 TMP_DECL;
43
44 ASSERT (dn > 0);
45 ASSERT (nn >= dn);
46 ASSERT (dp[dn-1] > 0);
47
48 while (dp[0] == 0)
49 {
50 ASSERT (np[0] == 0);
51 dp++;
52 np++;
53 dn--;
54 nn--;
55 }
56
57 if (dn == 1)
58 {
59 MPN_DIVREM_OR_DIVEXACT_1 (qp, np, nn, dp[0]);
60 return;
61 }
62
63 TMP_MARK;
64
65 qn = nn + 1 - dn;
66 count_trailing_zeros (shift, dp[0]);
67
68 if (shift > 0)
69 {
70 mp_size_t ss = (dn > qn) ? qn + 1 : dn;
71
72 tp = TMP_ALLOC_LIMBS (ss);
73 mpn_rshift (tp, dp, ss, shift);
74 dp = tp;
75
76 /* Since we have excluded dn == 1, we have nn > qn, and we need
77 to shift one limb beyond qn. */
78 wp = TMP_ALLOC_LIMBS (qn + 1);
79 mpn_rshift (wp, np, qn + 1, shift);
80 }
81 else
82 {
83 wp = TMP_ALLOC_LIMBS (qn);
84 MPN_COPY (wp, np, qn);
85 }
86
87 if (dn > qn)
88 dn = qn;
89
90 tp = TMP_ALLOC_LIMBS (mpn_bdiv_q_itch (qn, dn));
91 mpn_bdiv_q (qp, wp, qn, dp, dn, tp);
92 TMP_FREE;
93 }
94
95 #else
96
97 /* We use the Jebelean's bidirectional exact division algorithm. This is
98 somewhat naively implemented, with equal quotient parts done by 2-adic
99 division and truncating division. Since 2-adic division is faster, it
100 should be used for a larger chunk.
101
102 This code is horrendously ugly, in all sorts of ways.
103
104 * It was hacked without much care or thought, but with a testing program.
105 * It handles scratch space frivolously, and furthermore the itch function
106 is broken.
107 * Doesn't provide any measures to deal with mu_divappr_q's +3 error. We
108 have yet to provoke an error due to this, though.
109 * Algorithm selection leaves a lot to be desired. In particular, the choice
110 between DC and MU isn't a point, but we treat it like one.
111 * It makes the msb part 1 or 2 limbs larger than the lsb part, in spite of
112 that the latter is faster. We should at least reverse this, but perhaps
113 we should make the lsb part considerably larger. (How do we tune this?)
114 */
115
116 mp_size_t
mpn_divexact_itch(mp_size_t nn,mp_size_t dn)117 mpn_divexact_itch (mp_size_t nn, mp_size_t dn)
118 {
119 return nn + dn; /* FIXME this is not right */
120 }
121
122 void
mpn_divexact(mp_ptr qp,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_ptr scratch)123 mpn_divexact (mp_ptr qp,
124 mp_srcptr np, mp_size_t nn,
125 mp_srcptr dp, mp_size_t dn,
126 mp_ptr scratch)
127 {
128 mp_size_t qn;
129 mp_size_t nn0, qn0;
130 mp_size_t nn1, qn1;
131 mp_ptr tp;
132 mp_limb_t qml;
133 mp_limb_t qh;
134 int cnt;
135 mp_ptr xdp;
136 mp_limb_t di;
137 mp_limb_t cy;
138 gmp_pi1_t dinv;
139 TMP_DECL;
140
141 TMP_MARK;
142
143 qn = nn - dn + 1;
144
145 /* For small divisors, and small quotients, don't use Jebelean's algorithm. */
146 if (dn < DIVEXACT_JEB_THRESHOLD || qn < DIVEXACT_JEB_THRESHOLD)
147 {
148 tp = scratch;
149 MPN_COPY (tp, np, qn);
150 binvert_limb (di, dp[0]); di = -di;
151 dn = MIN (dn, qn);
152 mpn_sbpi1_bdiv_q (qp, tp, qn, dp, dn, di);
153 TMP_FREE;
154 return;
155 }
156
157 qn0 = ((nn - dn) >> 1) + 1; /* low quotient size */
158
159 /* If quotient is much larger than the divisor, the bidirectional algorithm
160 does not work as currently implemented. Fall back to plain bdiv. */
161 if (qn0 > dn)
162 {
163 if (BELOW_THRESHOLD (dn, DC_BDIV_Q_THRESHOLD))
164 {
165 tp = scratch;
166 MPN_COPY (tp, np, qn);
167 binvert_limb (di, dp[0]); di = -di;
168 dn = MIN (dn, qn);
169 mpn_sbpi1_bdiv_q (qp, tp, qn, dp, dn, di);
170 }
171 else if (BELOW_THRESHOLD (dn, MU_BDIV_Q_THRESHOLD))
172 {
173 tp = scratch;
174 MPN_COPY (tp, np, qn);
175 binvert_limb (di, dp[0]); di = -di;
176 mpn_dcpi1_bdiv_q (qp, tp, qn, dp, dn, di);
177 }
178 else
179 {
180 mpn_mu_bdiv_q (qp, np, qn, dp, dn, scratch);
181 }
182 TMP_FREE;
183 return;
184 }
185
186 nn0 = qn0 + qn0;
187
188 nn1 = nn0 - 1 + ((nn-dn) & 1);
189 qn1 = qn0;
190 if (LIKELY (qn0 != dn))
191 {
192 nn1 = nn1 + 1;
193 qn1 = qn1 + 1;
194 if (UNLIKELY (dp[dn - 1] == 1 && qn1 != dn))
195 {
196 /* If the leading divisor limb == 1, i.e. has just one bit, we have
197 to include an extra limb in order to get the needed overlap. */
198 /* FIXME: Now with the mu_divappr_q function, we should really need
199 more overlap. That indicates one of two things: (1) The test code
200 is not good. (2) We actually overlap too much by default. */
201 nn1 = nn1 + 1;
202 qn1 = qn1 + 1;
203 }
204 }
205
206 tp = TMP_ALLOC_LIMBS (nn1 + 1);
207
208 count_leading_zeros (cnt, dp[dn - 1]);
209
210 /* Normalize divisor, store into tmp area. */
211 if (cnt != 0)
212 {
213 xdp = TMP_ALLOC_LIMBS (qn1);
214 mpn_lshift (xdp, dp + dn - qn1, qn1, cnt);
215 }
216 else
217 {
218 xdp = (mp_ptr) dp + dn - qn1;
219 }
220
221 /* Shift dividend according to the divisor normalization. */
222 /* FIXME: We compute too much here for XX_divappr_q, but these functions'
223 interfaces want a pointer to the imaginative least significant limb, not
224 to the least significant *used* limb. Of course, we could leave nn1-qn1
225 rubbish limbs in the low part, to save some time. */
226 if (cnt != 0)
227 {
228 cy = mpn_lshift (tp, np + nn - nn1, nn1, cnt);
229 if (cy != 0)
230 {
231 tp[nn1] = cy;
232 nn1++;
233 }
234 }
235 else
236 {
237 /* FIXME: This copy is not needed for mpn_mu_divappr_q, except when the
238 mpn_sub_n right before is executed. */
239 MPN_COPY (tp, np + nn - nn1, nn1);
240 }
241
242 invert_pi1 (dinv, xdp[qn1 - 1], xdp[qn1 - 2]);
243 if (BELOW_THRESHOLD (qn1, DC_DIVAPPR_Q_THRESHOLD))
244 {
245 qp[qn0 - 1 + nn1 - qn1] = mpn_sbpi1_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, dinv.inv32);
246 }
247 else if (BELOW_THRESHOLD (qn1, MU_DIVAPPR_Q_THRESHOLD))
248 {
249 qp[qn0 - 1 + nn1 - qn1] = mpn_dcpi1_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, &dinv);
250 }
251 else
252 {
253 /* FIXME: mpn_mu_divappr_q doesn't handle qh != 0. Work around it with a
254 conditional subtraction here. */
255 qh = mpn_cmp (tp + nn1 - qn1, xdp, qn1) >= 0;
256 if (qh)
257 mpn_sub_n (tp + nn1 - qn1, tp + nn1 - qn1, xdp, qn1);
258 mpn_mu_divappr_q (qp + qn0 - 1, tp, nn1, xdp, qn1, scratch);
259 qp[qn0 - 1 + nn1 - qn1] = qh;
260 }
261 qml = qp[qn0 - 1];
262
263 binvert_limb (di, dp[0]); di = -di;
264
265 if (BELOW_THRESHOLD (qn0, DC_BDIV_Q_THRESHOLD))
266 {
267 MPN_COPY (tp, np, qn0);
268 mpn_sbpi1_bdiv_q (qp, tp, qn0, dp, qn0, di);
269 }
270 else if (BELOW_THRESHOLD (qn0, MU_BDIV_Q_THRESHOLD))
271 {
272 MPN_COPY (tp, np, qn0);
273 mpn_dcpi1_bdiv_q (qp, tp, qn0, dp, qn0, di);
274 }
275 else
276 {
277 mpn_mu_bdiv_q (qp, np, qn0, dp, qn0, scratch);
278 }
279
280 if (qml < qp[qn0 - 1])
281 mpn_decr_u (qp + qn0, 1);
282
283 TMP_FREE;
284 }
285 #endif
286