1 /* mpn_mu_div_qr, mpn_preinv_mu_div_qr.
2
3 Compute Q = floor(N / D) and R = N-QD. N is nn limbs and D is dn limbs and
4 must be normalized, and Q must be nn-dn limbs. The requirement that Q is
5 nn-dn limbs (and not nn-dn+1 limbs) was put in place in order to allow us to
6 let N be unmodified during the operation.
7
8 Contributed to the GNU project by Torbjorn Granlund.
9
10 THE FUNCTIONS IN THIS FILE ARE INTERNAL WITH MUTABLE INTERFACES. IT IS ONLY
11 SAFE TO REACH THEM THROUGH DOCUMENTED INTERFACES. IN FACT, IT IS ALMOST
12 GUARANTEED THAT THEY WILL CHANGE OR DISAPPEAR IN A FUTURE GMP RELEASE.
13
14 Copyright 2005, 2006, 2007, 2009, 2010 Free Software Foundation, Inc.
15
16 This file is part of the GNU MP Library.
17
18 The GNU MP Library is free software; you can redistribute it and/or modify
19 it under the terms of the GNU Lesser General Public License as published by
20 the Free Software Foundation; either version 3 of the License, or (at your
21 option) any later version.
22
23 The GNU MP Library is distributed in the hope that it will be useful, but
24 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
26 License for more details.
27
28 You should have received a copy of the GNU Lesser General Public License
29 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
30
31
32 /*
33 The idea of the algorithm used herein is to compute a smaller inverted value
34 than used in the standard Barrett algorithm, and thus save time in the
35 Newton iterations, and pay just a small price when using the inverted value
36 for developing quotient bits. This algorithm was presented at ICMS 2006.
37 */
38
39 /* CAUTION: This code and the code in mu_divappr_q.c should be edited in sync.
40
41 Things to work on:
42
43 * This isn't optimal when the quotient isn't needed, as it might take a lot
44 of space. The computation is always needed, though, so there is no time to
45 save with special code.
46
47 * The itch/scratch scheme isn't perhaps such a good idea as it once seemed,
48 demonstrated by the fact that the mpn_invertappr function's scratch needs
49 mean that we need to keep a large allocation long after it is needed.
50 Things are worse as mpn_mul_fft does not accept any scratch parameter,
51 which means we'll have a large memory hole while in mpn_mul_fft. In
52 general, a peak scratch need in the beginning of a function isn't
53 well-handled by the itch/scratch scheme.
54 */
55
56 #ifdef STAT
57 #undef STAT
58 #define STAT(x) x
59 #else
60 #define STAT(x)
61 #endif
62
63 #include <stdlib.h> /* for NULL */
64 #include "gmp.h"
65 #include "gmp-impl.h"
66
67
68 /* FIXME: The MU_DIV_QR_SKEW_THRESHOLD was not analysed properly. It gives a
69 speedup according to old measurements, but does the decision mechanism
70 really make sense? It seem like the quotient between dn and qn might be
71 what we really should be checking. */
72 #ifndef MU_DIV_QR_SKEW_THRESHOLD
73 #define MU_DIV_QR_SKEW_THRESHOLD 100
74 #endif
75
76 #ifdef CHECK /* FIXME: Enable in minithres */
77 #undef MU_DIV_QR_SKEW_THRESHOLD
78 #define MU_DIV_QR_SKEW_THRESHOLD 1
79 #endif
80
81
82 static mp_limb_t mpn_mu_div_qr2 (mp_ptr, mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t, mp_ptr);
83
84
85 mp_limb_t
mpn_mu_div_qr(mp_ptr qp,mp_ptr rp,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_ptr scratch)86 mpn_mu_div_qr (mp_ptr qp,
87 mp_ptr rp,
88 mp_srcptr np,
89 mp_size_t nn,
90 mp_srcptr dp,
91 mp_size_t dn,
92 mp_ptr scratch)
93 {
94 mp_size_t qn;
95 mp_limb_t cy, qh;
96
97 qn = nn - dn;
98 if (qn + MU_DIV_QR_SKEW_THRESHOLD < dn)
99 {
100 /* |______________|_ign_first__| dividend nn
101 |_______|_ign_first__| divisor dn
102
103 |______| quotient (prel) qn
104
105 |___________________| quotient * ignored-divisor-part dn-1
106 */
107
108 /* Compute a preliminary quotient and a partial remainder by dividing the
109 most significant limbs of each operand. */
110 qh = mpn_mu_div_qr2 (qp, rp + nn - (2 * qn + 1),
111 np + nn - (2 * qn + 1), 2 * qn + 1,
112 dp + dn - (qn + 1), qn + 1,
113 scratch);
114
115 /* Multiply the quotient by the divisor limbs ignored above. */
116 if (dn - (qn + 1) > qn)
117 mpn_mul (scratch, dp, dn - (qn + 1), qp, qn); /* prod is dn-1 limbs */
118 else
119 mpn_mul (scratch, qp, qn, dp, dn - (qn + 1)); /* prod is dn-1 limbs */
120
121 if (qh)
122 cy = mpn_add_n (scratch + qn, scratch + qn, dp, dn - (qn + 1));
123 else
124 cy = 0;
125 scratch[dn - 1] = cy;
126
127 cy = mpn_sub_n (rp, np, scratch, nn - (2 * qn + 1));
128 cy = mpn_sub_nc (rp + nn - (2 * qn + 1),
129 rp + nn - (2 * qn + 1),
130 scratch + nn - (2 * qn + 1),
131 qn + 1, cy);
132 if (cy)
133 {
134 qh -= mpn_sub_1 (qp, qp, qn, 1);
135 mpn_add_n (rp, rp, dp, dn);
136 }
137 }
138 else
139 {
140 qh = mpn_mu_div_qr2 (qp, rp, np, nn, dp, dn, scratch);
141 }
142
143 return qh;
144 }
145
146 static mp_limb_t
mpn_mu_div_qr2(mp_ptr qp,mp_ptr rp,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_ptr scratch)147 mpn_mu_div_qr2 (mp_ptr qp,
148 mp_ptr rp,
149 mp_srcptr np,
150 mp_size_t nn,
151 mp_srcptr dp,
152 mp_size_t dn,
153 mp_ptr scratch)
154 {
155 mp_size_t qn, in;
156 mp_limb_t cy, qh;
157 mp_ptr ip, tp;
158
159 ASSERT (dn > 1);
160
161 qn = nn - dn;
162
163 /* Compute the inverse size. */
164 in = mpn_mu_div_qr_choose_in (qn, dn, 0);
165 ASSERT (in <= dn);
166
167 #if 1
168 /* This alternative inverse computation method gets slightly more accurate
169 results. FIXMEs: (1) Temp allocation needs not analysed (2) itch function
170 not adapted (3) mpn_invertappr scratch needs not met. */
171 ip = scratch;
172 tp = scratch + in + 1;
173
174 /* compute an approximate inverse on (in+1) limbs */
175 if (dn == in)
176 {
177 MPN_COPY (tp + 1, dp, in);
178 tp[0] = 1;
179 mpn_invertappr (ip, tp, in + 1, NULL);
180 MPN_COPY_INCR (ip, ip + 1, in);
181 }
182 else
183 {
184 cy = mpn_add_1 (tp, dp + dn - (in + 1), in + 1, 1);
185 if (UNLIKELY (cy != 0))
186 MPN_ZERO (ip, in);
187 else
188 {
189 mpn_invertappr (ip, tp, in + 1, NULL);
190 MPN_COPY_INCR (ip, ip + 1, in);
191 }
192 }
193 #else
194 /* This older inverse computation method gets slightly worse results than the
195 one above. */
196 ip = scratch;
197 tp = scratch + in;
198
199 /* Compute inverse of D to in+1 limbs, then round to 'in' limbs. Ideally the
200 inversion function should do this automatically. */
201 if (dn == in)
202 {
203 tp[in + 1] = 0;
204 MPN_COPY (tp + in + 2, dp, in);
205 mpn_invertappr (tp, tp + in + 1, in + 1, NULL);
206 }
207 else
208 {
209 mpn_invertappr (tp, dp + dn - (in + 1), in + 1, NULL);
210 }
211 cy = mpn_sub_1 (tp, tp, in + 1, GMP_NUMB_HIGHBIT);
212 if (UNLIKELY (cy != 0))
213 MPN_ZERO (tp + 1, in);
214 MPN_COPY (ip, tp + 1, in);
215 #endif
216
217 qh = mpn_preinv_mu_div_qr (qp, rp, np, nn, dp, dn, ip, in, scratch + in);
218
219 return qh;
220 }
221
222 mp_limb_t
mpn_preinv_mu_div_qr(mp_ptr qp,mp_ptr rp,mp_srcptr np,mp_size_t nn,mp_srcptr dp,mp_size_t dn,mp_srcptr ip,mp_size_t in,mp_ptr scratch)223 mpn_preinv_mu_div_qr (mp_ptr qp,
224 mp_ptr rp,
225 mp_srcptr np,
226 mp_size_t nn,
227 mp_srcptr dp,
228 mp_size_t dn,
229 mp_srcptr ip,
230 mp_size_t in,
231 mp_ptr scratch)
232 {
233 mp_size_t qn;
234 mp_limb_t cy, cx, qh;
235 mp_limb_t r;
236 mp_size_t tn, wn;
237
238 #define tp scratch
239 #define scratch_out (scratch + tn)
240
241 qn = nn - dn;
242
243 np += qn;
244 qp += qn;
245
246 qh = mpn_cmp (np, dp, dn) >= 0;
247 if (qh != 0)
248 mpn_sub_n (rp, np, dp, dn);
249 else
250 MPN_COPY (rp, np, dn);
251
252 if (qn == 0)
253 return qh; /* Degenerate use. Should we allow this? */
254
255 while (qn > 0)
256 {
257 if (qn < in)
258 {
259 ip += in - qn;
260 in = qn;
261 }
262 np -= in;
263 qp -= in;
264
265 /* Compute the next block of quotient limbs by multiplying the inverse I
266 by the upper part of the partial remainder R. */
267 mpn_mul_n (tp, rp + dn - in, ip, in); /* mulhi */
268 cy = mpn_add_n (qp, tp + in, rp + dn - in, in); /* I's msb implicit */
269 ASSERT_ALWAYS (cy == 0);
270
271 qn -= in;
272
273 /* Compute the product of the quotient block and the divisor D, to be
274 subtracted from the partial remainder combined with new limbs from the
275 dividend N. We only really need the low dn+1 limbs. */
276
277 if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
278 mpn_mul (tp, dp, dn, qp, in); /* dn+in limbs, high 'in' cancels */
279 else
280 {
281 tn = mpn_mulmod_bnm1_next_size (dn + 1);
282 mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, in, scratch_out);
283 wn = dn + in - tn; /* number of wrapped limbs */
284 if (wn > 0)
285 {
286 cy = mpn_sub_n (tp, tp, rp + dn - wn, wn);
287 cy = mpn_sub_1 (tp + wn, tp + wn, tn - wn, cy);
288 cx = mpn_cmp (rp + dn - in, tp + dn, tn - dn) < 0;
289 ASSERT_ALWAYS (cx >= cy);
290 mpn_incr_u (tp, cx - cy);
291 }
292 }
293
294 r = rp[dn - in] - tp[dn];
295
296 /* Subtract the product from the partial remainder combined with new
297 limbs from the dividend N, generating a new partial remainder R. */
298 if (dn != in)
299 {
300 cy = mpn_sub_n (tp, np, tp, in); /* get next 'in' limbs from N */
301 cy = mpn_sub_nc (tp + in, rp, tp + in, dn - in, cy);
302 MPN_COPY (rp, tp, dn); /* FIXME: try to avoid this */
303 }
304 else
305 {
306 cy = mpn_sub_n (rp, np, tp, in); /* get next 'in' limbs from N */
307 }
308
309 STAT (int i; int err = 0;
310 static int errarr[5]; static int err_rec; static int tot);
311
312 /* Check the remainder R and adjust the quotient as needed. */
313 r -= cy;
314 while (r != 0)
315 {
316 /* We loop 0 times with about 69% probability, 1 time with about 31%
317 probability, 2 times with about 0.6% probability, if inverse is
318 computed as recommended. */
319 mpn_incr_u (qp, 1);
320 cy = mpn_sub_n (rp, rp, dp, dn);
321 r -= cy;
322 STAT (err++);
323 }
324 if (mpn_cmp (rp, dp, dn) >= 0)
325 {
326 /* This is executed with about 76% probability. */
327 mpn_incr_u (qp, 1);
328 cy = mpn_sub_n (rp, rp, dp, dn);
329 STAT (err++);
330 }
331
332 STAT (
333 tot++;
334 errarr[err]++;
335 if (err > err_rec)
336 err_rec = err;
337 if (tot % 0x10000 == 0)
338 {
339 for (i = 0; i <= err_rec; i++)
340 printf (" %d(%.1f%%)", errarr[i], 100.0*errarr[i]/tot);
341 printf ("\n");
342 }
343 );
344 }
345
346 return qh;
347 }
348
349 /* In case k=0 (automatic choice), we distinguish 3 cases:
350 (a) dn < qn: in = ceil(qn / ceil(qn/dn))
351 (b) dn/3 < qn <= dn: in = ceil(qn / 2)
352 (c) qn < dn/3: in = qn
353 In all cases we have in <= dn.
354 */
355 mp_size_t
mpn_mu_div_qr_choose_in(mp_size_t qn,mp_size_t dn,int k)356 mpn_mu_div_qr_choose_in (mp_size_t qn, mp_size_t dn, int k)
357 {
358 mp_size_t in;
359
360 if (k == 0)
361 {
362 mp_size_t b;
363 if (qn > dn)
364 {
365 /* Compute an inverse size that is a nice partition of the quotient. */
366 b = (qn - 1) / dn + 1; /* ceil(qn/dn), number of blocks */
367 in = (qn - 1) / b + 1; /* ceil(qn/b) = ceil(qn / ceil(qn/dn)) */
368 }
369 else if (3 * qn > dn)
370 {
371 in = (qn - 1) / 2 + 1; /* b = 2 */
372 }
373 else
374 {
375 in = (qn - 1) / 1 + 1; /* b = 1 */
376 }
377 }
378 else
379 {
380 mp_size_t xn;
381 xn = MIN (dn, qn);
382 in = (xn - 1) / k + 1;
383 }
384
385 return in;
386 }
387
388 mp_size_t
mpn_mu_div_qr_itch(mp_size_t nn,mp_size_t dn,int mua_k)389 mpn_mu_div_qr_itch (mp_size_t nn, mp_size_t dn, int mua_k)
390 {
391 mp_size_t itch_local = mpn_mulmod_bnm1_next_size (dn + 1);
392 mp_size_t in = mpn_mu_div_qr_choose_in (nn - dn, dn, mua_k);
393 mp_size_t itch_out = mpn_mulmod_bnm1_itch (itch_local, dn, in);
394
395 return in + itch_local + itch_out;
396 }
397
398 mp_size_t
mpn_preinv_mu_div_qr_itch(mp_size_t nn,mp_size_t dn,mp_size_t in)399 mpn_preinv_mu_div_qr_itch (mp_size_t nn, mp_size_t dn, mp_size_t in)
400 {
401 mp_size_t itch_local = mpn_mulmod_bnm1_next_size (dn + 1);
402 mp_size_t itch_out = mpn_mulmod_bnm1_itch (itch_local, dn, in);
403
404 return itch_local + itch_out;
405 }
406