mpz/stronglucas.c: New file

4 files changed, 188 insertions, 2 deletions

diff --git a/Makefile.am b/Makefile.am
index 563b23cee..2709ec02d 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -210,7 +210,8 @@ MPZ_OBJECTS = mpz/abs$U.lo mpz/add$U.lo mpz/add_ui$U.lo			\
   mpz/set_q$U.lo mpz/set_si$U.lo mpz/set_str$U.lo mpz/set_ui$U.lo	\
   mpz/setbit$U.lo							\
   mpz/size$U.lo mpz/sizeinbase$U.lo mpz/sqrt$U.lo			\
-  mpz/sqrtrem$U.lo mpz/sub$U.lo mpz/sub_ui$U.lo mpz/swap$U.lo		\
+  mpz/sqrtrem$U.lo mpz/stronglucas$U.lo mpz/sub$U.lo			\
+  mpz/sub_ui$U.lo mpz/swap$U.lo						\
   mpz/tdiv_ui$U.lo mpz/tdiv_q$U.lo mpz/tdiv_q_2exp$U.lo			\
   mpz/tdiv_q_ui$U.lo mpz/tdiv_qr$U.lo mpz/tdiv_qr_ui$U.lo		\
   mpz/tdiv_r$U.lo mpz/tdiv_r_2exp$U.lo mpz/tdiv_r_ui$U.lo		\
diff --git a/gmp-impl.h b/gmp-impl.h
index 55263ebf9..5c72c5968 100644
--- a/gmp-impl.h
+++ b/gmp-impl.h
@@ -1136,6 +1136,9 @@ __GMP_DECLSPEC mp_size_t mpn_fib2_ui (mp_ptr, mp_ptr, unsigned long);
 #define mpn_fib2m __MPN(fib2m)
 __GMP_DECLSPEC int mpn_fib2m (mp_ptr, mp_ptr, mp_srcptr, mp_size_t, mp_srcptr, mp_size_t);
 
+#define mpn_strongfibo __MPN(strongfibo)
+__GMP_DECLSPEC int mpn_strongfibo (mp_srcptr, mp_size_t, mp_ptr);
+
 /* Remap names of internal mpn functions.  */
 #define __clz_tab               __MPN(clz_tab)
 #define mpn_udiv_w_sdiv		__MPN(udiv_w_sdiv)
@@ -1704,6 +1707,9 @@ __GMP_DECLSPEC mp_size_t mpz_prodlimbs (mpz_ptr, mp_ptr, mp_size_t);
 #define mpz_oddfac_1  __gmpz_oddfac_1
 __GMP_DECLSPEC void mpz_oddfac_1 (mpz_ptr, mp_limb_t, unsigned);
 
+#define mpz_stronglucas  __gmpz_stronglucas
+__GMP_DECLSPEC int mpz_stronglucas (mpz_srcptr, mpz_ptr, mpz_ptr);
+
 #define mpz_lucas_mod  __gmpz_lucas_mod
 __GMP_DECLSPEC int mpz_lucas_mod (mpz_ptr, mpz_ptr, long, mp_bitcnt_t, mpz_srcptr, mpz_ptr, mpz_ptr);
 
diff --git a/mpz/Makefile.am b/mpz/Makefile.am
index c18622270..a5e1f5720 100644
--- a/mpz/Makefile.am
+++ b/mpz/Makefile.am
@@ -61,7 +61,8 @@ libmpz_la_SOURCES = aors.h aors_ui.h fits_s.h mul_i.h \
   powm_sec.c powm_ui.c pprime_p.c prodlimbs.c primorial_ui.c random.c random2.c \
   realloc.c realloc2.c remove.c roinit_n.c root.c rootrem.c rrandomb.c \
   scan0.c scan1.c set.c set_d.c set_f.c set_q.c set_si.c set_str.c \
-  set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c sub.c sub_ui.c \
+  set_ui.c setbit.c size.c sizeinbase.c sqrt.c sqrtrem.c stronglucas.c \
+  sub.c sub_ui.c \
   swap.c tdiv_ui.c tdiv_q.c tdiv_q_2exp.c tdiv_q_ui.c tdiv_qr.c \
   tdiv_qr_ui.c tdiv_r.c tdiv_r_2exp.c tdiv_r_ui.c tstbit.c ui_pow_ui.c \
   ui_sub.c urandomb.c urandomm.c xor.c
diff --git a/mpz/stronglucas.c b/mpz/stronglucas.c
new file mode 100644
index 000000000..1672d722f
--- /dev/null
+++ b/mpz/stronglucas.c
@@ -0,0 +1,178 @@
+/* mpz_stronglucas(n,reps) -- An implementation of the strong Lucas
+   primality test, using parameters as suggested by the BPSW test.
+
+   THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY.  THEY'RE ALMOST
+   CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN
+   FUTURE GNU MP RELEASES.
+
+Copyright 2018 Free Software Foundation, Inc.
+
+Contributed by Marco Bodrato.
+
+This file is part of the GNU MP Library.
+
+The GNU MP Library is free software; you can redistribute it and/or modify
+it under the terms of either:
+
+  * the GNU Lesser General Public License as published by the Free
+    Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+or
+
+  * the GNU General Public License as published by the Free Software
+    Foundation; either version 2 of the License, or (at your option) any
+    later version.
+
+or both in parallel, as here.
+
+The GNU MP Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+for more details.
+
+You should have received copies of the GNU General Public License and the
+GNU Lesser General Public License along with the GNU MP Library.  If not,
+see https://www.gnu.org/licenses/.  */
+
+#include "gmp-impl.h"
+#include "longlong.h"
+
+/* Returns an approximation of the sqare root of x.
+ * It gives:
+ *   limb_apprsqrt (x) ^ 2 <= x < (limb_apprsqrt (x)+1) ^ 2
+ * or
+ *   x <= limb_apprsqrt (x) ^ 2 <= x * 9/8
+ */
+static mp_limb_t
+limb_apprsqrt (mp_limb_t x)
+{
+  int s;
+
+  ASSERT (x > 2);
+  count_leading_zeros (s, x);
+  s = (GMP_LIMB_BITS - s) >> 1;
+  return ((CNST_LIMB(1) << s) + (x >> s)) >> 1;
+}
+
+/* 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.*/
+int
+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. */
+  long Q;
+  mpz_t T1, T2;
+
+  /* Test on the absolute value. */
+  mpz_roinit_n (n, PTR (x), ABSIZ (x));
+
+  ASSERT (mpz_odd_p (n));
+  /* 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));
+  /* (2^12 - 1) = 3^2 * 5 * 7 * 13		*/
+  ASSERT (D % 3 != 0 && D % 5 != 0 && D % 7 != 0);
+  if ((D % 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))
+    /* (-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)
+    /* (-11/n) = (n/11)	*/
+    D = 11; /* Q = 3 */
+  else if ((((D % 13 - (D % 13 >> 3)) & 7) > 4) ||
+	   (((D % 13 - (D % 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)
+    /* (-15/n) = (n/15) = (n/5)*(n/3)	*/
+    /* Here, (n/5) = 1, and		*/
+    /* (n/3) = -1, iff n = 2 (mod 3)	*/
+    D = 15; /* Q = 4 */
+#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))
+    D = 17; /* Q = -4 */
+#endif
+#else
+  if (mpz_kronecker_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;
+
+    if (mpz_perfect_square_p (n))
+      return 0; /* A square is composite. */
+
+    /* Check Ds up to square root (in case, n is prime)
+       or avoid overflows */
+    if (SIZ (n) == 1)
+      maxD = limb_apprsqrt (* PTR (n));
+    else if (BITS_PER_ULONG >= GMP_NUMB_BITS && SIZ (n) == 2)
+      mpn_sqrtrem (&maxD, (mp_ptr) NULL, PTR (n), 2);
+    else
+      maxD = GMP_NUMB_MAX;
+    maxD = MIN (maxD, ULONG_MAX);
+
+    D = GMP_NUMB_BITS % 16 == 0 ? (GMP_NUMB_BITS % 32 == 0 ? 17 : 15) : 5;
+
+    /* Search a D such that (D/n) = -1 in the sequence 5,-7,9,-11,.. */
+    /* For those Ds we have (D/n) = (n/|D|) */
+    /* FIXME: Should we loop only on prime Ds?	*/
+    /* The only interesting composite D is 15.	*/
+    do
+      {
+	if (D >= maxD)
+	  return 1;
+	D += 2;
+	jac_bit1 = 0;
+	JACOBI_MOD_OR_MODEXACT_1_ODD (jac_bit1, tl, PTR (n), SIZ (n), D);
+	if (tl == 0)
+	  return 0;
+      }
+    while (mpn_jacobi_base (tl, D, jac_bit1) == 1);
+  }
+
+  /* D= P^2 - 4Q; P = 1; Q = (1-D)/4 */
+  Q = (D & 2) ? (D >> 2) + 1 : -(long) (D >> 2);
+  /* ASSERT (mpz_si_kronecker ((D & 2) ? NEG_CAST (long, D) : D, n) == -1); */
+
+  /* n-(D/n) = n+1 = d*2^{b0}, with d = (n>>b0) | 1 */
+  b0 = mpz_scan0 (n, 0);
+
+  mpz_init (T1);
+  mpz_init (T2);
+
+  /* If Ud != 0 && Vd != 0 */
+  if (mpz_lucas_mod (V, Qk, Q, b0, n, T1, T2) == 0)
+    if (--b0 != 0)
+      do
+	{
+	  /* V_{2k} <- V_k ^ 2 - 2Q^k */
+	  mpz_mul (T2, V, V);
+	  mpz_submul_ui (T2, Qk, 2);
+	  mpz_tdiv_r (V, T2, n);
+	  if (SIZ (V) == 0 || --b0 == 0)
+	    break;
+	  /* Q^{2k} = (Q^k)^2 */
+	  mpz_mul (T2, Qk, Qk);
+	  mpz_tdiv_r (Qk, T2, n);
+	} while (1);
+
+  mpz_clear (T1);
+  mpz_clear (T2);
+
+  return (b0 != 0);
+}