[Gmp-commit] /var/hg/gmp-proj/mini-gmp: Implemented the other mpz_[cft]div_[q...
mercurial at gmplib.org
mercurial at gmplib.org
Fri Jan 13 22:02:40 CET 2012
details: /var/hg/gmp-proj/mini-gmp/rev/b327ec321322
changeset: 77:b327ec321322
user: Niels M?ller <nisse at lysator.liu.se>
date: Fri Jan 13 22:02:38 2012 +0100
description:
Implemented the other mpz_[cft]div_[qr]_2exp functions.
diffstat:
mini-gmp.c | 222 ++++++++++++++++++++++++++++++++++++++++++++++------
mini-gmp.h | 5 +
tests/t-div_2exp.c | 20 ++--
3 files changed, 208 insertions(+), 39 deletions(-)
diffs (truncated from 304 to 300 lines):
diff -r ce396c502706 -r b327ec321322 mini-gmp.c
--- a/mini-gmp.c Fri Jan 13 19:21:44 2012 +0100
+++ b/mini-gmp.c Fri Jan 13 22:02:38 2012 +0100
@@ -37,8 +37,6 @@
mpz_divisible_p
mpz_divisible_ui_p
mpz_export
- mpz_fdiv_q_2exp
- mpz_fdiv_r_2exp
mpz_import.
*/
@@ -2167,35 +2165,201 @@
mpz_div_qr (NULL, r, n, d, DIV_TRUNC);
}
+static int
+mpz_divisible_2exp_p (const mpz_t x, mp_bitcnt_t bit_index)
+{
+ mp_size_t xn, limb_index;
+ mp_ptr xp;
+ mp_limb_t bit;
+ mp_size_t i;
+
+ xn = GMP_ABS (x->_mp_size);
+ if (xn == 0)
+ return 1;
+
+ limb_index = bit_index / GMP_LIMB_BITS;
+ if (limb_index >= xn)
+ return 0;
+
+ xp = x->_mp_d;
+ for (i = limb_index; i-- > 0; )
+ if (xp[i] != 0)
+ return 0;
+
+ bit = (mp_limb_t) 1 << (bit_index % GMP_LIMB_BITS);
+ return (xp[limb_index] & (bit - 1)) == 0;
+}
+
+static void
+mpz_div_q_2exp (mpz_t q, const mpz_t u, mp_bitcnt_t bit_index,
+ enum mpz_div_round_mode mode)
+{
+ mp_size_t un, qn;
+ mp_size_t limb_cnt;
+ mp_ptr qp;
+ mp_limb_t adjust;
+
+ un = u->_mp_size;
+ if (un == 0)
+ {
+ q->_mp_size = 0;
+ return;
+ }
+ limb_cnt = bit_index / GMP_LIMB_BITS;
+ qn = GMP_ABS (un) - limb_cnt;
+
+ if ( (mode == DIV_FLOOR && un < 0) || (mode == DIV_CEIL && un > 0))
+ adjust = !mpz_divisible_2exp_p (u, bit_index);
+ else
+ adjust = 0;
+
+ if (qn <= 0)
+ qn = 0;
+
+ else
+ {
+ qp = MPZ_REALLOC (q, qn);
+
+ bit_index %= GMP_LIMB_BITS;
+ if (bit_index != 0)
+ {
+ mpn_rshift (qp, u->_mp_d + limb_cnt, qn, bit_index);
+ qn -= qp[qn - 1] == 0;
+ }
+ else
+ {
+ mpn_copyi (qp, u->_mp_d + limb_cnt, qn);
+ }
+ }
+
+ q->_mp_size = qn;
+
+ mpz_add_ui (q, q, adjust);
+ if (un < 0)
+ mpz_neg (q, q);
+}
+
+static void
+mpz_div_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t bit_index,
+ enum mpz_div_round_mode mode)
+{
+ mp_size_t us, un, rn;
+ mp_ptr rp;
+ mp_limb_t mask;
+
+ us = u->_mp_size;
+ if (us == 0 || bit_index == 0)
+ {
+ r->_mp_size = 0;
+ return;
+ }
+ rn = (bit_index + GMP_LIMB_BITS - 1) / GMP_LIMB_BITS;
+ assert (rn > 0);
+
+ rp = MPZ_REALLOC (r, rn);
+ un = GMP_ABS (us);
+
+ mask = GMP_LIMB_MAX >> (rn * GMP_LIMB_BITS - bit_index);
+
+ if (rn > un)
+ {
+ /* Quotient (with truncation) is zero, and remainder is
+ non-zero */
+ if ( (mode == DIV_FLOOR && us < 0) || (mode == DIV_CEIL && us > 0))
+ {
+ /* Have to negate and sign extend. */
+ mp_size_t i;
+ mp_limb_t cy;
+
+ for (cy = 1, i = 0; i < un; i++)
+ {
+ mp_limb_t s = ~u->_mp_d[i] + cy;
+ cy = s < cy;
+ rp[i] = s;
+ }
+ assert (cy == 0);
+ for (; i < rn - 1; i++)
+ rp[i] = GMP_LIMB_MAX;
+
+ rp[rn-1] = mask;
+ us = -us;
+ }
+ else
+ {
+ /* Just copy */
+ if (r != u)
+ mpn_copyi (rp, u->_mp_d, un);
+
+ rn = un;
+ }
+ }
+ else
+ {
+ if (r != u)
+ mpn_copyi (rp, u->_mp_d, rn - 1);
+
+ rp[rn-1] = u->_mp_d[rn-1] & mask;
+
+ if ( (mode == DIV_FLOOR && us < 0) || (mode == DIV_CEIL && us > 0))
+ {
+ /* If r != 0, compute 2^{bit_count} - r. */
+ mp_size_t i;
+
+ for (i = 0; i < rn && rp[i] == 0; i++)
+ ;
+ if (i < rn)
+ {
+ /* r > 0, need to flip sign. */
+ rp[i] = ~rp[i] + 1;
+ for (i++; i < rn; i++)
+ rp[i] = ~rp[i];
+
+ rp[rn-1] &= mask;
+
+ /* us is not used for anything else, so we can modify it
+ here to indicate flipped sign. */
+ us = -us;
+ }
+ }
+ }
+ rn = mpn_normalized_size (rp, rn);
+ r->_mp_size = us < 0 ? -rn : rn;
+}
+
+void
+mpz_cdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
+{
+ mpz_div_q_2exp (r, u, cnt, DIV_CEIL);
+}
+
+void
+mpz_fdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
+{
+ mpz_div_q_2exp (r, u, cnt, DIV_FLOOR);
+}
+
void
mpz_tdiv_q_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
{
- mp_size_t un, rn;
- mp_size_t limb_cnt;
- mp_ptr rp;
-
- un = u->_mp_size;
- limb_cnt = cnt / GMP_LIMB_BITS;
- rn = GMP_ABS (un) - limb_cnt;
- if (rn <= 0)
- rn = 0;
- else
- {
- rp = MPZ_REALLOC (r, rn);
-
- cnt %= GMP_LIMB_BITS;
- if (cnt != 0)
- {
- mpn_rshift (rp, u->_mp_d + limb_cnt, rn, cnt);
- rn -= rp[rn - 1] == 0;
- }
- else
- {
- mpn_copyi (rp, u->_mp_d + limb_cnt, rn);
- }
- }
-
- r->_mp_size = un >= 0 ? rn : -rn;
+ mpz_div_q_2exp (r, u, cnt, DIV_TRUNC);
+}
+
+void
+mpz_cdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
+{
+ mpz_div_r_2exp (r, u, cnt, DIV_CEIL);
+}
+
+void
+mpz_fdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
+{
+ mpz_div_r_2exp (r, u, cnt, DIV_FLOOR);
+}
+
+void
+mpz_tdiv_r_2exp (mpz_t r, const mpz_t u, mp_bitcnt_t cnt)
+{
+ mpz_div_r_2exp (r, u, cnt, DIV_TRUNC);
}
void
@@ -3025,7 +3189,7 @@
assert (limb_index < dn);
gmp_assert_nocarry (mpn_sub_1 (dp + limb_index, dp + limb_index,
- dn - limb_index, bit));
+ dn - limb_index, bit));
dn -= (dp[dn-1] == 0);
d->_mp_size = (d->_mp_size < 0) ? - dn : dn;
}
diff -r ce396c502706 -r b327ec321322 mini-gmp.h
--- a/mini-gmp.h Fri Jan 13 19:21:44 2012 +0100
+++ b/mini-gmp.h Fri Jan 13 22:02:38 2012 +0100
@@ -126,7 +126,12 @@
void mpz_fdiv_r (mpz_t, const mpz_t, const mpz_t);
void mpz_tdiv_r (mpz_t, const mpz_t, const mpz_t);
+void mpz_cdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
+void mpz_fdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
void mpz_tdiv_q_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
+void mpz_cdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
+void mpz_fdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
+void mpz_tdiv_r_2exp (mpz_t, const mpz_t, mp_bitcnt_t);
void mpz_divexact (mpz_t q, const mpz_t n, const mpz_t d);
diff -r ce396c502706 -r b327ec321322 tests/t-div_2exp.c
--- a/tests/t-div_2exp.c Fri Jan 13 19:21:44 2012 +0100
+++ b/tests/t-div_2exp.c Fri Jan 13 22:02:38 2012 +0100
@@ -33,25 +33,25 @@
for (i = 0; i < COUNT; i++)
{
unsigned j;
- for (j = 0; j < 1; j++)
+ for (j = 0; j < 6; j++)
{
static const enum hex_random_op ops[6] =
{
- /* OP_CDIV_Q_2, OP_CDIV_R_2,
- OP_FDIV_Q_2, OP_FDIV_R_2, */
- OP_TDIV_Q_2, /* OP_TDIV_R_2 */
+ OP_CDIV_Q_2, OP_CDIV_R_2,
+ OP_FDIV_Q_2, OP_FDIV_R_2,
+ OP_TDIV_Q_2, OP_TDIV_R_2
};
static const char *name[6] =
{
- /* "cdiv_q", "cdiv_r",
- "fdiv_q", "fdiv_r", */
- "tdiv_q", /* "tdiv_r" */
+ "cdiv_q", "cdiv_r",
+ "fdiv_q", "fdiv_r",
+ "tdiv_q", "tdiv_r"
};
static div_func * const div [6] =
{
- /* mpz_cdiv_q_2exp, mpz_cdiv_r_2exp,
- mpz_fdiv_q_2exp, mpz_fdiv_r_2exp, */
- mpz_tdiv_q_2exp, /* mpz_tdiv_r_2exp */
+ mpz_cdiv_q_2exp, mpz_cdiv_r_2exp,
+ mpz_fdiv_q_2exp, mpz_fdiv_r_2exp,
More information about the gmp-commit
mailing list