mpz/stronglucas.c (mpz_oddjacobi_ui): New helper function.

1 files changed, 39 insertions, 18 deletions

diff --git a/mpz/stronglucas.c b/mpz/stronglucas.c
index 350dd2a48..e9a71370e 100644
--- a/mpz/stronglucas.c
+++ b/mpz/stronglucas.c
@@ -55,6 +55,26 @@ limb_apprsqrt (mp_limb_t x)
   return ((CNST_LIMB(1) << s) + (x >> s)) >> 1;
 }
 
+static int
+mpz_oddjacobi_ui (mpz_t b, mp_limb_t a)
+{
+  mp_limb_t  b_rem;
+  int        result_bit1;
+
+  ASSERT (a & 1);
+  ASSERT (a > 1);
+  ASSERT (SIZ (b) > 0);
+  ASSERT ((*PTR (b) & 1) == 1);
+
+  result_bit1 = 0;
+  JACOBI_MOD_OR_MODEXACT_1_ODD (result_bit1, b_rem, PTR (b), SIZ (b), a);
+  if (UNLIKELY (b_rem == 0))
+    return 0;
+  else
+    return mpn_jacobi_base (b_rem, a, result_bit1);
+}
+
+
 /* Performs strong Lucas' test on x, with parameters suggested */
 /* for the BPSW test. Qk and V are passed to recycle variables. */
 /* Requires GCD (x,6) = 1.*/
@@ -64,6 +84,7 @@ mpz_stronglucas (mpz_srcptr x, mpz_ptr V, mpz_ptr Qk)
   mp_bitcnt_t b0;
   mpz_t n;
   mp_limb_t D; /* The absolute value is stored. */
+  mp_limb_t g;
   long Q;
   mpz_t T1, T2;
 
@@ -74,25 +95,25 @@ mpz_stronglucas (mpz_srcptr x, mpz_ptr V, mpz_ptr Qk)
   /* ASSERT (mpz_gcd_ui (NULL, n, 6) == 1);	*/
 #if GMP_NUMB_BITS % 16 == 0
   /* (2^12 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1)	*/
-  D = mpn_mod_34lsub1 (PTR (n), SIZ (n));
+  g = mpn_mod_34lsub1 (PTR (n), SIZ (n));
   /* (2^12 - 1) = 3^2 * 5 * 7 * 13		*/
-  ASSERT (D % 3 != 0 && D % 5 != 0 && D % 7 != 0);
-  if ((D % 5 & 2) != 0)
+  ASSERT (g % 3 != 0 && g % 5 != 0 && g % 7 != 0);
+  if ((g % 5 & 2) != 0)
     /* (5/n) = -1, iff n = 2 or 3 (mod 5)	*/
     /* D = 5; Q = -1 */
     return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
-  else if (! POW2_P (D % 7))
+  else if (! POW2_P (g % 7))
     /* (-7/n) = -1, iff n = 3,5 or 6 (mod 7)	*/
     D = 7; /* Q = 2 */
     /* (9/n) = -1, never: 9 = 3^2	*/
-  else if (mpz_kronecker_ui (n, 11) == -1)
+  else if (mpz_oddjacobi_ui (n, 11) == -1)
     /* (-11/n) = (n/11)	*/
     D = 11; /* Q = 3 */
-  else if ((((D % 13 - (D % 13 >> 3)) & 7) > 4) ||
-	   (((D % 13 - (D % 13 >> 3)) & 7) == 2))
+  else if ((((g % 13 - (g % 13 >> 3)) & 7) > 4) ||
+	   (((g % 13 - (g % 13 >> 3)) & 7) == 2))
     /* (13/n) = -1, iff n = 2,5,6,7,8 or 11 (mod 13)	*/
     D = 13; /* Q = -3 */
-  else if (D % 3 == 2)
+  else if (g % 3 == 2)
     /* (-15/n) = (n/15) = (n/5)*(n/3)	*/
     /* Here, (n/5) = 1, and		*/
     /* (n/3) = -1, iff n = 2 (mod 3)	*/
@@ -100,20 +121,20 @@ mpz_stronglucas (mpz_srcptr x, mpz_ptr V, mpz_ptr Qk)
 #if GMP_NUMB_BITS % 32 == 0
   /* (2^24 - 1) | (2^{GMP_NUMB_BITS*3/4} - 1)	*/
   /* (2^24 - 1) = (2^12 - 1) * 17 * 241		*/
-  else if (! POW2_P (D % 17) && ! POW2_P (17 - D % 17))
+  else if (! POW2_P (g % 17) && ! POW2_P (17 - g % 17))
     D = 17; /* Q = -4 */
 #endif
 #else
-  if (mpz_kronecker_ui (n, 5) == -1)
+  if (mpz_oddjacobi_ui (n, 5) == -1)
     return mpn_strongfibo (PTR (n), SIZ (n), PTR (V));
 #endif
   else
   {
-    mp_limb_t tl;
     mp_limb_t maxD;
-    int jac_bit1;
+    int jac;
 
-    if (UNLIKELY (mpz_perfect_square_p (n)))
+    /* n is odd, to possibly be a square, n % 8 = 1 is needed. */
+    if (((*PTR (n) & 6) == 0) && UNLIKELY (mpz_perfect_square_p (n)))
       return 0; /* A square is composite. */
 
     /* Check Ds up to square root (in case, n is prime)
@@ -137,12 +158,12 @@ mpz_stronglucas (mpz_srcptr x, mpz_ptr V, mpz_ptr Qk)
 	if (UNLIKELY (D >= maxD))
 	  return 1;
 	D += 2;
-	jac_bit1 = 0;
-	JACOBI_MOD_OR_MODEXACT_1_ODD (jac_bit1, tl, PTR (n), SIZ (n), D);
-	if (UNLIKELY (tl == 0))
-	  return 0;
+	jac = mpz_oddjacobi_ui (n, D);
       }
-    while (mpn_jacobi_base (tl, D, jac_bit1) == 1);
+    while (jac == 1);
+
+    if (UNLIKELY (jac == 0))
+      return 0;
   }
 
   /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */