From 30b7851cdf97503a572b7e9a1ab3938f67df81e8 Mon Sep 17 00:00:00 2001
From: Paul Eggert <eggert@cs.ucla.edu>
Date: Sat, 14 Jun 2025 23:10:25 -0700
Subject: [PATCH] factor: always use Baillie-PSW

* src/factor.c (USE_BAILLIE_PSW): New constant.
(prime_p, prime2_p): Use it, i.e., always use Baillie-PSW.
Do so by using mp_prime_p.  Do not tell GCC these functions
are pure; the pure mark was present only to pacify GCC
and is no longer needed now that thes functions call mp_prime_p.
---
 src/factor.c | 30 ++++++++++++++++++++++++++----
 1 file changed, 26 insertions(+), 4 deletions(-)

diff --git a/src/factor.c b/src/factor.c
index 0542c0e6f1..82b283f6a2 100644
--- a/src/factor.c
+++ b/src/factor.c
@@ -1261,9 +1261,18 @@ mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
   return false;
 }
 
-/* Lucas' prime test.  The number of iterations vary greatly, up to a few dozen
-   have been observed.  The average seem to be about 2.  */
-static bool ATTRIBUTE_PURE
+static bool mp_prime_p (mpz_t);
+
+/* Use the Baillie-PSW primality test for all integers.
+   Code that is currently never executed because it is protected by
+   the equivalent of "if (!USE_BAILLIE_PSW)" dates back to when
+   factor.c used an inferior primality test.
+   FIXME: Either remove the !USE_BAILLIE_PSW code, or fix it to be
+   more efficient than the USE_BAILLIE_PSW code, by special casing
+   Baillie-PSW for single- and/or double-word args.  */
+enum { USE_BAILLIE_PSW = true };
+
+static bool
 prime_p (mp_limb_t n)
 {
   mp_limb_t a_prim, one, ni;
@@ -1276,6 +1285,12 @@ prime_p (mp_limb_t n)
   if (n < (mp_limb_t) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME)
     return true;
 
+  if (USE_BAILLIE_PSW)
+    {
+      mpz_t mn = MPZ_ROINIT_N (&n, 1);
+      return mp_prime_p (mn);
+    }
+
   /* Precomputation for Miller-Rabin.  */
   int k = stdc_trailing_zeros (n - 1);
   mp_limb_t q = (n - 1) >> k;
@@ -1344,7 +1359,7 @@ prime_p (mp_limb_t n)
     }
 }
 
-static bool ATTRIBUTE_PURE
+static bool
 prime2_p (mp_limb_t n1, mp_limb_t n0)
 {
   mp_limb_t q[2], nm1[2];
@@ -1358,6 +1373,13 @@ prime2_p (mp_limb_t n1, mp_limb_t n0)
   if (n1 == 0)
     return prime_p (n0);
 
+  if (USE_BAILLIE_PSW)
+    {
+      uuint un = make_uuint (n1, n0);
+      mpz_t mn = MPZ_ROINIT_N (un.uu, 2);
+      return mp_prime_p (mn);
+    }
+
   nm1[1] = n1 - (n0 == 0);
   nm1[0] = n0 - 1;
   if (nm1[0] == 0)
-- 
2.47.3