[Gmp-commit] /var/hg/gmp: Make mpz_remove use mpn_remove.
mercurial at gmplib.org
mercurial at gmplib.org
Thu Oct 25 21:47:05 CEST 2012
details: /var/hg/gmp/rev/ce92190c7189
changeset: 15084:ce92190c7189
user: Torbjorn Granlund <tege at gmplib.org>
date: Thu Oct 25 21:47:02 2012 +0200
description:
Make mpz_remove use mpn_remove.
diffstat:
ChangeLog | 6 ++
mpn/generic/remove.c | 30 ++++++++---
mpz/remove.c | 133 ++++++++++++++++++++++++++------------------------
3 files changed, 97 insertions(+), 72 deletions(-)
diffs (233 lines):
diff -r d3a9efa9d582 -r ce92190c7189 ChangeLog
--- a/ChangeLog Wed Oct 17 23:03:05 2012 +0200
+++ b/ChangeLog Thu Oct 25 21:47:02 2012 +0200
@@ -1,3 +1,9 @@
+2012-10-25 Torbjorn Granlund <tege at gmplib.org>
+
+ * mpn/generic/remove.c (mpn_bdiv_qr_wrap): New functions.
+ (mpn_remove): Call mpn_bdiv_qr_wrap.
+ * mpz/remove.c: Enable suppressed mpn_remove call.
+
2012-10-17 Torbjorn Granlund <tege at gmplib.org>
* mpz/powm_ui.c (mpz_powm_ui): Deflect to mpz_powm for large exponent.
diff -r d3a9efa9d582 -r ce92190c7189 mpn/generic/remove.c
--- a/mpn/generic/remove.c Wed Oct 17 23:03:05 2012 +0200
+++ b/mpn/generic/remove.c Thu Oct 25 21:47:02 2012 +0200
@@ -47,6 +47,25 @@
* If we allow ourselves to clobber U, we could save the other of qp and qp2.
*/
+/* FIXME: We need to wrap mpn_bdiv_qr due to the itch interface. This need
+ indicates a flaw in the current itch mechanism: Which operands not greater
+ than un,un will incur the worst itch? We need a parallel foo_maxitch set
+ of functions. */
+void
+mpn_bdiv_qr_wrap (mp_ptr qp, mp_ptr rp,
+ mp_srcptr np, mp_size_t nn,
+ mp_srcptr dp, mp_size_t dn)
+{
+ mp_ptr scratch_out;
+ TMP_DECL;
+
+ TMP_MARK;
+ scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (nn, dn));
+ mpn_bdiv_qr (qp, rp, np, nn, dp, dn, scratch_out);
+
+ TMP_FREE;
+}
+
mp_bitcnt_t
mpn_remove (mp_ptr wp, mp_size_t *wn,
mp_ptr up, mp_size_t un, mp_ptr vp, mp_size_t vn,
@@ -55,7 +74,7 @@
mp_ptr pwpsp[LOG];
mp_size_t pwpsn[LOG];
mp_size_t npowers;
- mp_ptr tp, qp, np, pp, qp2, scratch_out;
+ mp_ptr tp, qp, np, pp, qp2;
mp_size_t pn, nn, qn, i;
mp_bitcnt_t pwr;
TMP_DECL;
@@ -74,11 +93,6 @@
pp = vp;
pn = vn;
- /* FIXME: This allocation need indicate a flaw in the current itch mechanism:
- Which operands not greater than un,un will incur the worst itch? We need
- a parallel foo_maxitch set of functions. */
- scratch_out = TMP_ALLOC_LIMBS (mpn_bdiv_qr_itch (un + 1, (un + 1) >> 1));
-
MPN_COPY (qp, up, un);
qn = un;
@@ -86,7 +100,7 @@
while (qn >= pn)
{
qp[qn] = 0;
- mpn_bdiv_qr (qp2, tp, qp, qn + 1, pp, pn, scratch_out);
+ mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pp, pn);
if (!mpn_zero_p (tp, pn))
break; /* could not divide by V^npowers */
@@ -125,7 +139,7 @@
continue; /* V^i would bring us past cap */
qp[qn] = 0;
- mpn_bdiv_qr (qp2, tp, qp, qn + 1, pp, pn, scratch_out);
+ mpn_bdiv_qr_wrap (qp2, tp, qp, qn + 1, pp, pn);
if (!mpn_zero_p (tp, pn))
continue; /* could not divide by V^i */
diff -r d3a9efa9d582 -r ce92190c7189 mpz/remove.c
--- a/mpz/remove.c Wed Oct 17 23:03:05 2012 +0200
+++ b/mpz/remove.c Thu Oct 25 21:47:02 2012 +0200
@@ -35,79 +35,84 @@
fp0 = fp[0];
if (UNLIKELY ((afn <= (fp0 == 1)) /* mpz_cmpabs_ui (f, 1) <= 0 */
- | (sn == 0))) {
- /* f = 0 or f = +- 1 or src = 0 */
- if (afn == 0)
- DIVIDE_BY_ZERO;
- mpz_set (dest, src);
- return 0;
- }
-
-#if 0
- if ((fp0 & 1) == 1) { /* f is odd */
- mp_ptr dp;
- mp_size_t dn;
-
- dn = ABS (sn);
- dp = MPZ_REALLOC (dest, dn);
-
- pwr = mpn_remove (dp, &dn, PTR(src), dn, PTR(f), afn, ~(mp_bitcnt_t) 0);
-
- SIZ (dest) = ((pwr & (fn < 0)) ^ (sn < 0)) ? -dn : dn;
- } else
-#endif
- if (afn == (fp0 == 2)) { /* mpz_cmpabs_ui (f, 2) == 0 */
- pwr = mpz_scan1 (src, 0);
- mpz_div_2exp (dest, src, pwr);
- if (pwr & (SIZ (f) < 0)) /*((pwr % 2 == 1) && (SIZ (f) < 0))*/
- mpz_neg (dest, dest);
- } else { /* f != +-2 */
-
- mpz_t fpow[GMP_LIMB_BITS]; /* Really MP_SIZE_T_BITS */
- mpz_t x, rem;
- int p;
-
- /* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0). It is an
- upper bound of the result we're seeking. We could also shift down the
- operands so that they become odd, to make intermediate values smaller. */
-
- mpz_init (rem);
- mpz_init (x);
-
- pwr = 0;
- mpz_init (fpow[0]);
- mpz_set (fpow[0], f);
- mpz_set (dest, src);
-
- /* Divide by f, f^2, ..., f^(2^k) until we get a remainder for f^(2^k). */
- for (p = 0;; p++) {
- mpz_tdiv_qr (x, rem, dest, fpow[p]);
- if (SIZ (rem) != 0)
- break;
- mpz_init (fpow[p + 1]);
- mpz_mul (fpow[p + 1], fpow[p], fpow[p]);
- mpz_set (dest, x);
+ | (sn == 0)))
+ {
+ /* f = 0 or f = +- 1 or src = 0 */
+ if (afn == 0)
+ DIVIDE_BY_ZERO;
+ mpz_set (dest, src);
+ return 0;
}
- pwr = (1L << p) - 1;
+ if ((fp0 & 1) == 1)
+ { /* f is odd */
+ mp_ptr dp;
+ mp_size_t dn;
- mpz_clear (fpow[p]);
+ dn = ABS (sn);
+ dp = MPZ_REALLOC (dest, dn);
- /* Divide by f^(2^(k-1)), f^(2^(k-2)), ..., f for all divisors that give a
- zero remainder. */
- while (--p >= 0) {
- mpz_tdiv_qr (x, rem, dest, fpow[p]);
- if (SIZ (rem) == 0)
+ pwr = mpn_remove (dp, &dn, PTR(src), dn, PTR(f), afn, ~(mp_bitcnt_t) 0);
+
+ SIZ (dest) = ((pwr & (fn < 0)) ^ (sn < 0)) ? -dn : dn;
+ }
+ else if (afn == (fp0 == 2))
+ { /* mpz_cmpabs_ui (f, 2) == 0 */
+ pwr = mpz_scan1 (src, 0);
+ mpz_div_2exp (dest, src, pwr);
+ if (pwr & (SIZ (f) < 0)) /*((pwr % 2 == 1) && (SIZ (f) < 0))*/
+ mpz_neg (dest, dest);
+ }
+ else
+ { /* f != +-2 */
+ mpz_t fpow[GMP_LIMB_BITS]; /* Really MP_SIZE_T_BITS */
+ mpz_t x, rem;
+ int p;
+
+ /* We could perhaps compute mpz_scan1(src,0)/mpz_scan1(f,0). It is an
+ upper bound of the result we're seeking. We could also shift down the
+ operands so that they become odd, to make intermediate values
+ smaller. */
+
+ mpz_init (rem);
+ mpz_init (x);
+
+ pwr = 0;
+ mpz_init (fpow[0]);
+ mpz_set (fpow[0], f);
+ mpz_set (dest, src);
+
+ /* Divide by f, f^2 ... f^(2^k) until we get a remainder for f^(2^k). */
+ for (p = 0;; p++)
{
- pwr += 1L << p;
+ mpz_tdiv_qr (x, rem, dest, fpow[p]);
+ if (SIZ (rem) != 0)
+ break;
+ mpz_init (fpow[p + 1]);
+ mpz_mul (fpow[p + 1], fpow[p], fpow[p]);
mpz_set (dest, x);
}
+
+ pwr = (1L << p) - 1;
+
mpz_clear (fpow[p]);
+
+ /* Divide by f^(2^(k-1)), f^(2^(k-2)), ..., f for all divisors that give
+ a zero remainder. */
+ while (--p >= 0)
+ {
+ mpz_tdiv_qr (x, rem, dest, fpow[p]);
+ if (SIZ (rem) == 0)
+ {
+ pwr += 1L << p;
+ mpz_set (dest, x);
+ }
+ mpz_clear (fpow[p]);
+ }
+
+ mpz_clear (x);
+ mpz_clear (rem);
}
- mpz_clear (x);
- mpz_clear (rem);
- }
-
return pwr;
}
More information about the gmp-commit
mailing list