[Gmp-commit] /var/hg/gmp: Rewrote broken handling of the case asize == 1 in m...
mercurial at gmplib.org
mercurial at gmplib.org
Sat May 21 23:49:13 CEST 2011
details: /var/hg/gmp/rev/fcdc50208fad
changeset: 14196:fcdc50208fad
user: Niels M?ller <nisse at lysator.liu.se>
date: Sat May 21 23:49:09 2011 +0200
description:
Rewrote broken handling of the case asize == 1 in mpz_jacobi.
diffstat:
ChangeLog | 10 ++++++
mpz/jacobi.c | 87 +++++++++++++++++-------------------------------------
tests/mpz/t-jac.c | 17 +++++++++-
3 files changed, 52 insertions(+), 62 deletions(-)
diffs (161 lines):
diff -r ddd56c4f2b69 -r fcdc50208fad ChangeLog
--- a/ChangeLog Fri May 20 20:23:55 2011 +0200
+++ b/ChangeLog Sat May 21 23:49:09 2011 +0200
@@ -1,3 +1,13 @@
+2011-05-21 Niels Möller <nisse at lysator.liu.se>
+
+ * mpz/jacobi.c (mpz_jacobi): The handling of asize == 1 was
+ broken. Rewrote it.
+
+ * tests/mpz/t-jac.c (mpz_nextprime_step): Sanity check that prime
+ candidate and step has no common factor.
+ (check_data): Added some test cases related to the asize == 1 case
+ in mpz_jacobi.
+
2011-05-20 Niels Möller <nisse at lysator.liu.se>
* gmp-impl.h: Jacobi-related prototypes.
diff -r ddd56c4f2b69 -r fcdc50208fad mpz/jacobi.c
--- a/mpz/jacobi.c Fri May 20 20:23:55 2011 +0200
+++ b/mpz/jacobi.c Sat May 21 23:49:09 2011 +0200
@@ -125,72 +125,39 @@
return mpn_jacobi_base (alow, blow, result_bit1);
}
- bits = mpn_jacobi_init (alow, blow, (result_bit1>>1) & 1);
-
if (asize == 1)
{
- /* We need least significant bits of the quotient */
- mp_limb_t b1 = mpn_mod_1 (bsrcp+1, bsize-1, alow);
- mp_limb_t q;
- udiv_qrnnd (q, blow, b1, blow, alow);
+ /* Logic copied from mpz_ui_kronecker */
+ if (alow == 1)
+ return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */
+
+ if (btwos > 0)
+ {
+ /* Only bit 1 of blow is used below. */
+ if (btwos == GMP_NUMB_BITS - 1)
+ blow = bsrcp[1] << 1;
+ else
+ blow >> btwos;
+ }
- bits = mpn_jacobi_update (bits, 0, q & 3);
+ else if ( (alow & 1) == 0)
+ {
+ unsigned atwos;
+ count_trailing_zeros (atwos, alow);
+ alow >>= atwos;
+ result_bit1 ^= JACOBI_TWOS_U_BIT1 (atwos, blow);
+ }
- if (bits >= 16)
- MP_LIMB_T_SWAP (alow, blow);
+ if (alow == 1)
+ return JACOBI_BIT1_TO_PN (result_bit1); /* (1/b)=1 */
- if (blow == 1)
- {
- /* FIXME: Do this in some better way? Currently,
- mpn_jacobi_finish requires that one number is reduced to
- zero. */
- bits = mpn_jacobi_update (bits, 1, alow & 3);
- return mpn_jacobi_finish (bits);
- }
- else
- return mpn_jacobi_base (alow, blow, bits << 1);
-
-#if 0
- /* FIXME: Is it better to avoid the use of mpn_jacobi_update for
- this special case? The below should work, if it hasn't
- suffered bit rot. */
- if (alow & 1)
- {
- result_bit1 ^= JACOBI_RECIP_UU_BIT1(alow, blow);
+ /* (a/b*2^n) = (b*2^n / a) * RECIP(a,b) */
+ result_bit1 ^= JACOBI_RECIP_UU_BIT1 (alow, blow);
+ JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, blow, bsrcp, bsize, alow);
+ return mpn_jacobi_base (blow, alow, result_bit1);
+ }
- if (alow == 1)
- return JACOBI_BIT1_TO_PN (result_bit1);
-
- if (bsize > 1)
- JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, blow, bsrcp, bsize, alow);
-
- return mpn_jacobi_base (blow, alow, result_bit1);
- }
- else
- {
- if (alow & 2)
- {
- /* We need least significant bits of the quotient */
- if (bsize > 1)
- {
- mp_limb_t b1 = mpn_mod_1 (bsrcp+1, bsize-1, alow);
- mp_limb_t q;
- udiv_qrnnd (q, blow, b1, blow, alow);
- /* Sign change is q (b-1)/2 + q (q-1) / 2
- = q (r-1)/2 + q (q+1)/2 */
- result_bit1 ^= ((q << 1) & blow) ^ (q << 1) ^q;
- }
- }
- else
- blow = mpn_mod_1 (bsrcp, bsize, alow);
-
- if (blow == 1)
- return JACOBI_BIT1_TO_PN (result_bit1);
-
- return mpn_jacobi_base (alow, blow, result_bit1);
- }
-#endif
- }
+ bits = mpn_jacobi_init (alow, blow, (result_bit1>>1) & 1);
/* Allocation strategy: When one operand is much larger than the
other, we currently have to allocate space for the entire
diff -r ddd56c4f2b69 -r fcdc50208fad tests/mpz/t-jac.c
--- a/tests/mpz/t-jac.c Fri May 20 20:23:55 2011 +0200
+++ b/tests/mpz/t-jac.c Sat May 21 23:49:09 2011 +0200
@@ -632,7 +632,15 @@
"451716845976689892447895811408978421929", -1 },
{ "67674091930576781943923596701346271058970643542491743605048620644676477275152701774960868941561652032482173612421015",
"4902678867794567120224500687210807069172039735", 0 },
- { "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683283672592", "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683290481773", 1 }
+ { "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683283672592", "2666617146103764067061017961903284334497474492754652499788571378062969111250584288683585223600172138551198546085281683290481773", 1 },
+
+ /* Exersizes the case asize == 1, btwos > 0 in mpz_jacobi. */
+ { "804609", "421248363205206617296534688032638102314410556521742428832362659824", 1 } ,
+ { "4190209", "2239744742177804210557442048984321017460028974602978995388383905961079286530650825925074203175536427000", 1 },
+
+ /* Exersizes the case asize == 1, btwos = 63 in mpz_jacobi
+ (relevant when GMP_LIMB_BITS == 64). */
+ { "17311973299000934401", "1675975991242824637446753124775689449936871337036614677577044717424700351103148799107651171694863695242089956242888229458836426332300124417011114380886016", 1 },
};
int i;
@@ -862,7 +870,7 @@
mp_size_t pn;
mp_bitcnt_t nbits;
unsigned incr;
- mpz_t step;
+ mpz_t step, gcd;
TMP_SDECL;
ASSERT_ALWAYS (mpz_sgn (step_in) > 0);
@@ -904,6 +912,11 @@
return;
}
+ mpz_init (gcd);
+ mpz_gcd (gcd, p, step);
+ ASSERT_ALWAYS (mpz_cmp_ui (gcd, 1) == 0);
+ mpz_clear (gcd);
+
pn = SIZ(p);
count_leading_zeros (cnt, PTR(p)[pn - 1]);
nbits = pn * GMP_NUMB_BITS - (cnt - GMP_NAIL_BITS);
More information about the gmp-commit
mailing list