mpn/generic/mu_div_qr.c

5    nn-dn limbs (and not nn-dn+1 limbs) was put in place in order to allow us to
31 or both in parallel, as here.
33 The GNU MP Library is distributed in the hope that it will be useful, but
45    than used in the standard Barrett algorithm, and thus save time in the
50 /* CAUTION: This code and the code in mu_divappr_q.c should be edited in sync.
62     which means we'll have a large memory hole while in mpn_mul_fft.  In
63     general, a peak scratch need in the beginning of a function isn't
86 #ifdef CHECK				/* FIXME: Enable in minithres */
166   mp_size_t qn, in;
175   in = mpn_mu_div_qr_choose_in (qn, dn, 0);
176   ASSERT (in <= dn);
183   tp = scratch + in + 1;
185   /* compute an approximate inverse on (in+1) limbs */
186   if (dn == in)
188       MPN_COPY (tp + 1, dp, in);
190       mpn_invertappr (ip, tp, in + 1, tp + in + 1);
191       MPN_COPY_INCR (ip, ip + 1, in);
195       cy = mpn_add_1 (tp, dp + dn - (in + 1), in + 1, 1);
197 	MPN_ZERO (ip, in);
200 	  mpn_invertappr (ip, tp, in + 1, tp + in + 1);
201 	  MPN_COPY_INCR (ip, ip + 1, in);
208   tp = scratch + in;
210   /* Compute inverse of D to in+1 limbs, then round to 'in' limbs.  Ideally the
212   if (dn == in)
214       tp[in + 1] = 0;
215       MPN_COPY (tp + in + 2, dp, in);
216       mpn_invertappr (tp, tp + in + 1, in + 1, NULL);
220       mpn_invertappr (tp, dp + dn - (in + 1), in + 1, NULL);
222   cy = mpn_sub_1 (tp, tp, in + 1, GMP_NUMB_HIGHBIT);
224     MPN_ZERO (tp + 1, in);
225   MPN_COPY (ip, tp + 1, in);
228   qh = mpn_preinv_mu_div_qr (qp, rp, np, nn, dp, dn, ip, in, scratch + in);
241 		      mp_size_t in,
268       if (qn < in)
270 	  ip += in - qn;
271 	  in = qn;
273       np -= in;
274       qp -= in;
278       mpn_mul_n (tp, rp + dn - in, ip, in);		/* mulhi  */
279       cy = mpn_add_n (qp, tp + in, rp + dn - in, in);	/* I's msb implicit */
282       qn -= in;
288       if (BELOW_THRESHOLD (in, MUL_TO_MULMOD_BNM1_FOR_2NXN_THRESHOLD))
289 	mpn_mul (tp, dp, dn, qp, in);		/* dn+in limbs, high 'in' cancels */
293 	  mpn_mulmod_bnm1 (tp, tn, dp, dn, qp, in, scratch_out);
294 	  wn = dn + in - tn;			/* number of wrapped limbs */
299 	      cx = mpn_cmp (rp + dn - in, tp + dn, tn - dn) < 0;
305       r = rp[dn - in] - tp[dn];
309       if (dn != in)
311 	  cy = mpn_sub_n (tp, np, tp, in);	/* get next 'in' limbs from N */
312 	  cy = mpn_sub_nc (tp + in, rp, tp + in, dn - in, cy);
317 	  cy = mpn_sub_n (rp, np, tp, in);	/* get next 'in' limbs from N */
361    (a) dn < qn:         in = ceil(qn / ceil(qn/dn))
362    (b) dn/3 < qn <= dn: in = ceil(qn / 2)
363    (c) qn < dn/3:       in = qn
364    In all cases we have in <= dn.
369   mp_size_t in;
378 	  in = (qn - 1) / b + 1;	/* ceil(qn/b) = ceil(qn / ceil(qn/dn)) */
382 	  in = (qn - 1) / 2 + 1;	/* b = 2 */
386 	  in = (qn - 1) / 1 + 1;	/* b = 1 */
393       in = (xn - 1) / k + 1;
396   return in;
402   mp_size_t in = mpn_mu_div_qr_choose_in (nn - dn, dn, mua_k);
403   mp_size_t itch_preinv = mpn_preinv_mu_div_qr_itch (nn, dn, in);
404   mp_size_t itch_invapp = mpn_invertappr_itch (in + 1) + in + 2; /* 3in + 4 */
407   return in + MAX (itch_invapp, itch_preinv);
411 mpn_preinv_mu_div_qr_itch (mp_size_t nn, mp_size_t dn, mp_size_t in)
414   mp_size_t itch_out = mpn_mulmod_bnm1_itch (itch_local, dn, in);