[Gmp-commit] /var/hg/gmp-proj/mini-gmp: 3 new changesets
mercurial at gmplib.org
mercurial at gmplib.org
Wed Jan 4 21:26:36 CET 2012
details: /var/hg/gmp-proj/mini-gmp/rev/4bde2db4ad2f
changeset: 43:4bde2db4ad2f
user: Niels M?ller <nisse at lysator.liu.se>
date: Wed Jan 04 19:53:34 2012 +0100
description:
Implemented mpz_invert and mpz_powm.
details: /var/hg/gmp-proj/mini-gmp/rev/f776169353f3
changeset: 44:f776169353f3
user: Niels M?ller <nisse at lysator.liu.se>
date: Wed Jan 04 19:56:16 2012 +0100
description:
Tests for powm, and logops (forgotten in earlier commit).
details: /var/hg/gmp-proj/mini-gmp/rev/b753004ece2c
changeset: 45:b753004ece2c
user: Niels M?ller <nisse at lysator.liu.se>
date: Wed Jan 04 21:25:37 2012 +0100
description:
Improved preinv division, and use it in mpz_powm.
struct gmp_div_qr_1_inverse: generalized and renamed to gmp_div_inverse.
mpn_div_qr_2_invert, mpn_div_qr_invert: New functions.
mpn_div_qr_2_preinv: New function.
mpn_div_qr_preinv: New function. For dn > 2, requires that caller
normalizes d. Still shifts n as needed.
mpn_div_qr: Rewritten to use mpn_div_qr_preinv.
mpz_powm: Use mpn_div_qr_preinv.
diffstat:
.hgignore | 1 +
mini-gmp.c | 310 +++++++++++++++++++++++++++++++++++++++++++--------
mini-gmp.h | 3 +
tests/Makefile | 2 +-
tests/hex-random.c | 118 +++++++++++++------
tests/hex-random.h | 11 +-
tests/mini-random.c | 18 +-
tests/mini-random.h | 4 +-
tests/t-div.c | 2 +-
tests/t-logops.c | 75 ++++++++++++
tests/t-powm.c | 54 +++++++++
11 files changed, 492 insertions(+), 106 deletions(-)
diffs (truncated from 856 to 300 lines):
diff -r 0a5e174f762f -r b753004ece2c .hgignore
--- a/.hgignore Tue Jan 03 22:31:57 2012 +0100
+++ b/.hgignore Wed Jan 04 21:25:37 2012 +0100
@@ -10,5 +10,6 @@
^tests/t-div$
^tests/t-double$
^tests/t-gcd$
+^tests/t-powm$
^tests/t-logops$
^tests/t-str$
diff -r 0a5e174f762f -r b753004ece2c mini-gmp.c
--- a/mini-gmp.c Tue Jan 03 22:31:57 2012 +0100
+++ b/mini-gmp.c Wed Jan 04 21:25:37 2012 +0100
@@ -674,31 +674,83 @@
return m;
}
-struct gmp_div_qr_1_inverse
+struct gmp_div_inverse
{
/* Normalization shift count. */
unsigned shift;
- /* Normalized divisor */
- mp_limb_t d;
- /* Inverse */
+ /* Normalized divisor (d0 unused for mpn_div_qr_1) */
+ mp_limb_t d1, d0;
+ /* Inverse, for 2/1 or 3/2. */
mp_limb_t di;
};
static void
-mpn_div_qr_1_invert (struct gmp_div_qr_1_inverse *inv, mp_limb_t d)
+mpn_div_qr_1_invert (struct gmp_div_inverse *inv, mp_limb_t d)
{
unsigned shift;
+
+ assert (d > 0);
gmp_clz (shift, d);
inv->shift = shift;
- inv->d = d << shift;
- inv->di = mpn_invert_limb (inv->d);
+ inv->d1 = d << shift;
+ inv->di = mpn_invert_limb (inv->d1);
+}
+
+static void
+mpn_div_qr_2_invert (struct gmp_div_inverse *inv,
+ mp_limb_t d1, mp_limb_t d0)
+{
+ unsigned shift;
+
+ assert (d1 > 0);
+ gmp_clz (shift, d1);
+ inv->shift = shift;
+ if (shift > 0)
+ {
+ d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
+ d0 <<= shift;
+ }
+ inv->d1 = d1;
+ inv->d0 = d0;
+ inv->di = mpn_invert_3by2 (d1, d0);
+}
+
+static void
+mpn_div_qr_invert (struct gmp_div_inverse *inv,
+ mp_srcptr dp, mp_size_t dn)
+{
+ assert (dn > 0);
+
+ if (dn == 1)
+ mpn_div_qr_1_invert (inv, dp[0]);
+ else if (dn == 2)
+ mpn_div_qr_2_invert (inv, dp[1], dp[0]);
+ else
+ {
+ unsigned shift;
+ mp_limb_t d1, d0;
+
+ d1 = dp[dn-1];
+ d0 = dp[dn-2];
+ assert (d1 > 0);
+ gmp_clz (shift, d1);
+ inv->shift = shift;
+ if (shift > 0)
+ {
+ d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
+ d0 = (d0 << shift) | (dp[dn-3] >> (GMP_LIMB_BITS - shift));
+ }
+ inv->d1 = d1;
+ inv->d0 = d0;
+ inv->di = mpn_invert_3by2 (d1, d0);
+ }
}
/* Not matching current public gmp interface, rather corresponding to
the sbpi1_div_* functions. */
static mp_limb_t
mpn_div_qr_1_preinv (mp_ptr qp, mp_srcptr np, mp_size_t nn,
- const struct gmp_div_qr_1_inverse *inv)
+ const struct gmp_div_inverse *inv)
{
mp_limb_t d, di;
mp_limb_t r;
@@ -713,8 +765,7 @@
else
r = 0;
-
- d = inv->d;
+ d = inv->d1;
di = inv->di;
while (nn-- > 0)
{
@@ -735,40 +786,42 @@
{
assert (d > 0);
+ /* Special case for powers of two. */
if (d > 1 && (d & (d-1)) == 0)
{
- /* Power of two */
unsigned shift;
mp_limb_t r = np[0] & (d-1);
gmp_ctz (shift, d);
- mpn_rshift (qp, np, nn, shift);
+ if (qp)
+ mpn_rshift (qp, np, nn, shift);
return r;
}
else
{
- struct gmp_div_qr_1_inverse inv;
+ struct gmp_div_inverse inv;
mpn_div_qr_1_invert (&inv, d);
return mpn_div_qr_1_preinv (qp, np, nn, &inv);
}
}
static void
-mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
- mp_limb_t d1, mp_limb_t d0)
+mpn_div_qr_2_preinv (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
+ const struct gmp_div_inverse *inv)
{
unsigned shift;
mp_size_t i;
- mp_limb_t di, r1, r0;
+ mp_limb_t d1, d0, di, r1, r0;
mp_ptr tp;
assert (nn >= 2);
-
- gmp_clz (shift, d1);
+ shift = inv->shift;
+ d1 = inv->d1;
+ d0 = inv->d0;
+ di = inv->di;
+
if (shift > 0)
{
- d1 = (d1 << shift) | (d0 >> (GMP_LIMB_BITS - shift));
- d0 <<= shift;
tp = gmp_xalloc_limbs (nn);
r1 = mpn_lshift (tp, np, nn, shift);
np = tp;
@@ -776,9 +829,6 @@
else
r1 = 0;
- di = mpn_invert_3by2 (d1, d0);
-
- /* np += nn - 1; */
r0 = np[nn - 1];
for (i = nn - 2; i >= 0; i--)
@@ -804,6 +854,19 @@
rp[0] = r0;
}
+#if 0
+static void
+mpn_div_qr_2 (mp_ptr qp, mp_ptr rp, mp_srcptr np, mp_size_t nn,
+ mp_limb_t d1, mp_limb_t d0)
+{
+ struct gmp_div_inverse inv;
+ assert (nn >= 2);
+
+ mpn_div_qr_2_invert (&inv, d1, d0);
+ mpn_div_qr_2_preinv (qp, rp, np, nn, &inv);
+}
+#endif
+
static void
mpn_div_qr_pi1 (mp_ptr qp,
mp_ptr np, mp_size_t nn, mp_limb_t n1,
@@ -866,48 +929,61 @@
}
static void
-mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
+mpn_div_qr_preinv (mp_ptr qp, mp_ptr np, mp_size_t nn,
+ mp_srcptr dp, mp_size_t dn,
+ const struct gmp_div_inverse *inv)
{
assert (dn > 0);
assert (nn >= dn);
if (dn == 1)
- np[0] = mpn_div_qr_1 (qp, np, nn, dp[0]);
+ np[0] = mpn_div_qr_1_preinv (qp, np, nn, inv);
else if (dn == 2)
- mpn_div_qr_2 (qp, np, np, nn, dp[1], dp[0]);
+ mpn_div_qr_2_preinv (qp, np, np, nn, inv);
else
{
- mp_ptr tp = NULL;
- mp_limb_t d1, d0, di;
mp_limb_t nh;
unsigned shift;
- d1 = dp[dn - 1];
-
- gmp_clz (shift, d1);
+ assert (dp[dn-1] & GMP_LIMB_HIGHBIT);
+
+ shift = inv->shift;
if (shift > 0)
- {
- tp = gmp_xalloc_limbs (dn);
- gmp_assert_nocarry (mpn_lshift (tp, dp, dn, shift));
- nh = mpn_lshift (np, np, nn, shift);
- dp = tp;
- d1 = dp[dn - 1];
- }
+ nh = mpn_lshift (np, np, nn, shift);
else
nh = 0;
- d0 = dp[dn-2];
- di = mpn_invert_3by2 (d1, d0);
- mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, di);
+ assert (inv->d1 == dp[dn-1]);
+ assert (inv->d0 == dp[dn-2]);
+
+ mpn_div_qr_pi1 (qp, np, nn, nh, dp, dn, inv->di);
if (shift > 0)
- {
- free (tp);
- gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
- }
+ gmp_assert_nocarry (mpn_rshift (np, np, dn, shift));
}
}
+static void
+mpn_div_qr (mp_ptr qp, mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn)
+{
+ struct gmp_div_inverse inv;
+ mp_ptr tp = NULL;
+
+ assert (dn > 0);
+ assert (nn >= dn);
+
+ mpn_div_qr_invert (&inv, dp, dn);
+ if (dn > 2 && inv.shift > 0)
+ {
+ tp = gmp_xalloc_limbs (dn);
+ gmp_assert_nocarry (mpn_lshift (tp, dp, dn, inv.shift));
+ dp = tp;
+ }
+ mpn_div_qr_preinv (qp, np, nn, dp, dn, &inv);
+ if (tp)
+ free (tp);
+}
+
static unsigned
mpn_base_power_of_two_p (unsigned b)
{
@@ -992,7 +1068,7 @@
the end. */
static size_t
mpn_limb_get_str (unsigned char *sp, mp_limb_t w,
- const struct gmp_div_qr_1_inverse *binv)
+ const struct gmp_div_inverse *binv)
{
mp_size_t i;
for (i = 0; w > 0; i++)
@@ -1002,7 +1078,7 @@
h = w >> (GMP_LIMB_BITS - binv->shift);
l = w << binv->shift;
- gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d, binv->di);
+ gmp_udiv_qrnnd_preinv (w, r, h, l, binv->d1, binv->di);
assert ( (r << (GMP_LIMB_BITS - binv->shift)) == 0);
r >>= binv->shift;
More information about the gmp-commit
mailing list