diff options
author | tege <tege@gmplib.org> | 2001-11-14 20:22:51 +0100 |
---|---|---|
committer | tege <tege@gmplib.org> | 2001-11-14 20:22:51 +0100 |
commit | 9c0d6df1a9572211f3b37dfeb77be4d3b42d1bff (patch) | |
tree | bb2a0189a099afa5b39b8e70243236aa3b0a0468 /mpz | |
parent | b6ba500c63df36ec936085dad2b8edf8a3ff88ca (diff) | |
download | gmp-9c0d6df1a9572211f3b37dfeb77be4d3b42d1bff.tar.gz |
Adjust for negative b, efter exponentiation done. Add
(still disabled) code for handling negative exponents. Misc cleanups.
Diffstat (limited to 'mpz')
-rw-r--r-- | mpz/powm.c | 94 |
1 files changed, 57 insertions, 37 deletions
diff --git a/mpz/powm.c b/mpz/powm.c index 2e8f058ef..68de435bc 100644 --- a/mpz/powm.c +++ b/mpz/powm.c @@ -144,6 +144,9 @@ phi (mp_limb_t t) #define POWM_THRESHOLD ((8 * KARATSUBA_SQR_THRESHOLD) / 3) #endif +#undef HANDLE_NEGATIVE_EXPONENT +#undef REDUCE_EXPONENT + void #ifndef BERKELEY_MP mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m) @@ -151,15 +154,18 @@ mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m) pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) #endif /* BERKELEY_MP */ { - mp_limb_t invm, c; - mp_size_t bn, mn, xn; - unsigned long int enb; mp_ptr xp, tp, qp, gp, this_gp; mp_srcptr bp, ep, mp; + mp_size_t bn, es, en, mn, xn; + mp_limb_t invm, c; + unsigned long int enb; mp_size_t i, K, j, l, k; int m_zero_cnt, e_zero_cnt; int sh; int use_redc; +#if HANDLE_NEGATIVE_EXPONENT + mpz_t new_b; +#endif #if REDUCE_EXPONENT mpz_t new_e; #endif @@ -170,18 +176,33 @@ pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) if (mn == 0) DIVIDE_BY_ZERO; - if (SIZ (e) <= 0) + es = SIZ (e); + if (es <= 0) { - /* Exponent is zero, result is 1 mod m, i.e., 1 or 0 depending on if m - equals 1. */ - SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1; - PTR(r)[0] = 1; - return; + if (es == 0) + { + /* Exponent is zero, result is 1 mod m, i.e., 1 or 0 depending on if + m equals 1. */ + SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1; + PTR(r)[0] = 1; + return; + } +#if HANDLE_NEGATIVE_EXPONENT + MPZ_TMP_INIT (new_b, ABSIZ (b) + 1); + + if (! mpz_invert (new_b, b, m)) + DIVIDE_BY_ZERO; + b = new_b; + es = -es; +#else + DIVIDE_BY_ZERO; +#endif } + en = es; #if REDUCE_EXPONENT /* Reduce exponent by dividing it by phi(m) when m small. */ - if (mn == 1 && mp[0] < 0xffffffffL && SIZ (e) * BITS_PER_MP_LIMB > 150) + if (mn == 1 && mp[0] < 0x7fffffffL && en * BITS_PER_MP_LIMB > 150) { MPZ_TMP_INIT (new_e, 2); mpz_mod_ui (new_e, e, phi (mp[0])); @@ -200,6 +221,8 @@ pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) } else { + /* Normalize m (i.e. make its most significant bit set) as required by + division functions below. */ count_leading_zeros (m_zero_cnt, mp[mn - 1]); if (m_zero_cnt != 0) { @@ -212,8 +235,8 @@ pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) /* Determine optimal value of k, the number of exponent bits we look at at a time. */ - count_leading_zeros (e_zero_cnt, PTR(e)[SIZ(e) - 1]); - enb = SIZ (e) * BITS_PER_MP_LIMB - e_zero_cnt; /* number of bits of exponent */ + count_leading_zeros (e_zero_cnt, PTR(e)[en - 1]); + enb = en * BITS_PER_MP_LIMB - e_zero_cnt; /* number of bits of exponent */ k = 1; K = 2; while (2 * enb > K * (2 + k * (3 + k))) @@ -225,56 +248,45 @@ pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) tp = TMP_ALLOC_LIMBS (2 * mn + 1); qp = TMP_ALLOC_LIMBS (mn + 1); - gp = __GMP_ALLOCATE_FUNC_TYPE (K / 2 * mn, mp_limb_t); + gp = __GMP_ALLOCATE_FUNC_LIMBS (K / 2 * mn); /* Compute x*R^n where R=2^BITS_PER_MP_LIMB. */ bn = ABSIZ (b); bp = PTR(b); - /* Handle |b| >= m by computing b mod m. While it is not strictly necessary + /* Handle |b| >= m by computing b mod m. FIXME: It is not strictly necessary for speed or correctness to do this when b and m have the same number of - limbs, it simplifies the "else" clause's handling of b < 0. */ + limbs, perhaps remove mpn_cmp call. */ if (bn > mn || (bn == mn && mpn_cmp (bp, mp, mn) >= 0)) { - /* Reduce possibly huge base while moving it to gp[0]. */ + /* Reduce possibly huge base while moving it to gp[0]. Use a function + call to reduce, since we don't want the quotient allocation to + live until function return. */ if (use_redc) { reduce (tp + mn, bp, bn, mp, mn); /* b mod m */ - if (SIZ (b) < 0) - mpn_sub_n (tp + mn, mp, tp + mn, mn); MPN_ZERO (tp, mn); mpn_tdiv_qr (qp, gp, 0L, tp, 2 * mn, mp, mn); /* unnormnalized! */ } else { reduce (gp, bp, bn, mp, mn); - if (SIZ (b) < 0) - mpn_sub_n (gp, mp, gp, mn); } } else { - /* |b| < m. Possibly, |b| << m. */ + /* |b| < m. We pad out operands to become mn limbs, which simplifies + the rest of the function, but slows things down when the |b| << m. */ if (use_redc) { MPN_ZERO (tp, mn); - if (SIZ (b) < 0) - mpn_sub (tp + mn, mp, mn, bp, bn); - else - { - MPN_COPY (tp + mn, bp, bn); - MPN_ZERO (tp + mn + bn, mn - bn); - } + MPN_COPY (tp + mn, bp, bn); + MPN_ZERO (tp + mn + bn, mn - bn); mpn_tdiv_qr (qp, gp, 0L, tp, 2 * mn, mp, mn); } else { - if (SIZ (b) < 0) - mpn_sub (gp, mp, mn, bp, bn); - else - { - MPN_COPY (gp, bp, bn); - MPN_ZERO (gp + bn, mn - bn); - } + MPN_COPY (gp, bp, bn); + MPN_ZERO (gp + bn, mn - bn); } } @@ -299,7 +311,7 @@ pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) /* Start the real stuff. */ ep = PTR (e); - i = SIZ (e) - 1; /* current index */ + i = en - 1; /* current index */ c = ep[i]; /* current limb */ sh = BITS_PER_MP_LIMB - e_zero_cnt; /* significant bits in ep[i] */ sh -= k; /* index of lower bit of ep[i] to take into account */ @@ -430,10 +442,18 @@ pow (mpz_srcptr b, mpz_srcptr e, mpz_srcptr m, mpz_ptr r) } xn = mn; MPN_NORMALIZE (xp, xn); + + if ((ep[0] & 1) && SIZ(b) < 0 && xn != 0) + { + mp = PTR(m); /* want original, unnormalized m */ + mpn_sub (xp, mp, mn, xp, xn); + xn = mn; + MPN_NORMALIZE (xp, xn); + } MPZ_REALLOC (r, xn); SIZ (r) = xn; MPN_COPY (PTR(r), xp, xn); - __GMP_FREE_FUNC_TYPE (gp, K / 2 * mn, mp_limb_t); + __GMP_FREE_FUNC_LIMBS (gp, K / 2 * mn); TMP_FREE (marker); } |