Subject: Re: bug#78476: GNU 'factor' problems with 128-bit word

From: Paul Eggert <[email-masked]>
Date: Wed, 9 Jul 2025 22:44:32 -0700

On 2025-05-20 00:04, Torbjörn Granlund wrote:
> I am too busy to take care of this.  Sorry.

Understandable. Thanks for taking the time to look at my email.

> Preferably, this problem should be fully understood.

Yes, that'd be better but it appears that nobody has the time to do so. 
However, we can address this in a different way, by removing the buggy 
code instead of understanding it. I did so by re-reading this bug report 
along with your January proposal[1], looking at factor.c again, and 
committing changes mostly along the lines you suggested, as follows:

* Merge the January proposed factor.c into GNU coreutils.

* Port to mini-gmp.c, which lacks mpn_tdiv_qr.

* Don't use uintmax_t or unsigned __int128 as a word; just use the same 
mp_limb_t type that GMP uses.

* Remove the USE_SQUFOF code (which was unused).

* Don't try to prove primality; trust Baillie-PSW as being good enough. 
Use mpz_probab_prime_p for that, even for numbers that fit into two 
limbs or less.

* Port to unlikely platforms lacking 32- or 64-bit words.

* Minor speedups, e.g., convert strings to uuint faster, and decrement 
MR_REPS from 25 to 24 as that's a sweet spot in GMP.

* Statically check for unsupported (and unlikely) build configurations.

* Add some comments.

Two FIXMEs remain in factor.c. To address the first one, would GMP 
accept longlong.h changes to make it standalone even with 
-Wunused-macros? Currently longlong.h is not standalone in that sense, 
even if LONGLONG_STANDALONE is defined. The changes should be minor.

You can grab all the coreutils updates from Savannah by cloning [2] as 
usual. For convenience in reviewing I am attaching the factors-related 
changes between the old and new Coreutils as a single patch. Comments 
are of course welcome.

[1]: https://lists.gnu.org/r/coreutils/2025-01/msg00000.html
[2]: https://git.savannah.gnu.org/git/coreutils.git

diff --git a/NEWS b/NEWS
index 87ba2659c..eda971bb6 100644
--- a/NEWS
+++ b/NEWS
@@ -2,6 +2,10 @@ GNU coreutils NEWS                                    -*- outline -*-
 
 * Noteworthy changes in release ?.? (????-??-??) [?]
 
+** Changes in behavior
+
+  'factor' is now much faster at identifying large prime numbers.
+
 ** Bug fixes
 
   cksum was not compilable by Apple LLVM 10.0.0 x86-64, which
diff --git a/doc/coreutils.texi b/doc/coreutils.texi
index 89a0b5972..7fc2db571 100644
--- a/doc/coreutils.texi
+++ b/doc/coreutils.texi
@@ -19158,7 +19158,7 @@ $ factor --exponents 3000
 
 If the number to be factored is small (less than @math{2^{127}} on
 typical machines), @command{factor} uses a faster algorithm.
-For example, on a circa-2017 Intel Xeon Silver 4116, factoring the
+For example, on a circa-2021 Intel Xeon W-1350, factoring the
 product of the eighth and ninth Mersenne primes (approximately
 @math{2^{92}}) takes about 4 ms of CPU time:
 
@@ -19169,22 +19169,32 @@ $ n=$(echo "$M8 * $M9" | bc)
 $ bash -c "time factor $n"
 4951760154835678088235319297: 2147483647 2305843009213693951
 
-real	0m0.004s
+real	0m0.006s
 user	0m0.004s
-sys	0m0.000s
+sys	0m0.002s
 @end example
 
 For larger numbers, @command{factor} uses a slower algorithm.  On the
 same platform, factoring the eighth Fermat number @math{2^{256} + 1}
-takes about 14 seconds, and the slower algorithm would have taken
-about 750 ms to factor @math{2^{127} - 3} instead of the 50 ms needed by
-the faster algorithm.
+takes about 6400 ms.  However, large primes are identified quickly:
+it takes just 420 ms to factor the Mersenne prime @math{2^{11213} - 1}
+into itself and 1.
 
-Factoring large numbers is, in general, hard.  The Pollard-Brent rho
+Factoring large numbers is, in general, hard.  The Pollard--Brent rho
 algorithm used by @command{factor} is particularly effective for
-numbers with relatively small factors.  If you wish to factor large
-numbers which do not have small factors (for example, numbers which
-are the product of two large primes), other methods are far better.
+numbers with relatively small factors.  Other methods are far better
+for factoring large composite numbers that lack relatively small
+factors, such as numbers that are the product of two large primes.
+
+When testing whether a number is prime, @command{factor} uses the
+Baillie--PSW primality heuristic for speed.  Although Baillie--PSW has
+not been proven to reject all composite numbers, the probability of a
+Baillie--PSW pseudoprime is astronomically small and no such
+pseudoprime has been discovered despite decades of searching by
+mathematicians.  Any such pseudoprime is greater than @math{2^{64}}.
+For more, see: Baillie R, Fiori A, Wagstaff Jr SS.@: Strengthening the
+Baillie--PSW primality test.@: @i{Math Comp}.@: 2021;90:1931--1955.@:
+@uref{https://doi.org/10.1090/mcom/3616}.
 
 @exitstatus
 
diff --git a/src/factor.c b/src/factor.c
index 948bf0527..44c0deb6a 100644
--- a/src/factor.c
+++ b/src/factor.c
@@ -19,9 +19,10 @@
    Arbitrary-precision code adapted by James Youngman from Torbjörn
    Granlund's factorize.c, from GNU MP version 4.2.2.
    In 2012, the core was rewritten by Torbjörn Granlund and Niels Möller.
+   In 2025, Torbjörn Granlund and Paul Eggert sped up primality checking.
    Contains code from GNU MP.  */
 
-/* Efficiently factor numbers that fit in one or two words (word = wide_uint),
+/* Efficiently factor numbers that fit in one or two words (word = mp_limb_t),
    or, with GMP, numbers of any size.
 
   Code organization:
@@ -41,38 +42,29 @@
         division since the primes table store inverses modulo the word base.
         (The GMP variant of this code doesn't make use of the precomputed
         inverses, but instead relies on GMP for fast divisibility testing.)
-    (2) Check the nature of any non-factored part using Miller-Rabin for
-        detecting composites, and Lucas for detecting primes.
-    (3) Factor any remaining composite part using the Pollard-Brent rho
-        algorithm or if USE_SQUFOF is defined to 1, try that first.
-        Status of found factors are checked again using Miller-Rabin and Lucas.
+    (2) Check the nature of any non-factored part using Baillie-PSW.
+    (3) Factor any remaining composite part using Pollard-Brent rho.
+        Status of found factors are checked again using Baillie-PSW.
 
-    We prefer using Hensel norm in the divisions, not the more familiar
+    Prefer using Hensel norm in the divisions, not the more familiar
     Euclidean norm, since the former leads to much faster code.  In the
-    Pollard-Brent rho code and the prime testing code, we use Montgomery's
-    trick of multiplying all n-residues by the word base, allowing cheap Hensel
-    reductions mod n.
+    Pollard-Brent rho code, use Montgomery's trick of multiplying
+    all n-residues by the word base, allowing cheap Hensel reductions mod n.
 
-    The GMP code uses an algorithm that can be considerably slower;
-    for example, on a circa-2017 Intel Xeon Silver 4116, factoring
-    2^{127}-3 takes about 50 ms with the two-word algorithm but would
-    take about 750 ms with the GMP code.
+    The GMP code uses an algorithm that can be considerably slower.
 
   Improvements:
 
-    * Use modular inverses also for exact division in the Lucas code, and
-      elsewhere.  A problem is to locate the inverses not from an index, but
-      from a prime.  We might instead compute the inverse on-the-fly.
+    * Use modular inverses also for exact division.  A problem is to
+      locate the inverses not from an index, but from a prime.
+      We might instead compute the inverse on-the-fly.
 
     * Tune trial division table size (not forgetting that this is a standalone
       program where the table will be read from secondary storage for
       each invocation).
 
-    * Implement less naive powm, using k-ary exponentiation for k = 3 or
-      perhaps k = 4.
-
-    * Try to speed trial division code for single wide_uint numbers, i.e., the
-      code using DIVBLOCK.  It currently runs at 2 cycles per prime (Intel SBR,
+    * Try to speed trial division code for single word numbers, i.e., the
+      code using divblock.  It currently runs at 2 cycles per prime (Intel SBR,
       IBR), 3 cycles per prime (AMD Stars) and 5 cycles per prime (AMD BD) when
       using gcc 4.6 and 4.7.  Some software pipelining should help; 1, 2, and 4
       respectively cycles ought to be possible.
@@ -80,26 +72,9 @@
     * The redcify function could be vastly improved by using (plain Euclidean)
       pre-inversion (such as GMP's invert_limb) and udiv_qrnnd_preinv (from
       GMP's gmp-impl.h).  The redcify2 function could be vastly improved using
-      similar methods.  These functions currently dominate run time when using
-      the -w option.
+      similar methods.
 */
 
-/* Whether to recursively factor to prove primality,
-   or run faster probabilistic tests.  */
-#ifndef PROVE_PRIMALITY
-# define PROVE_PRIMALITY 1
-#endif
-
-/* Faster for certain ranges but less general.  */
-#ifndef USE_SQUFOF
-# define USE_SQUFOF 0
-#endif
-
-/* Output SQUFOF statistics.  */
-#ifndef STAT_SQUFOF
-# define STAT_SQUFOF 0
-#endif
-
 
 #include <config.h>
 #include <getopt.h>
@@ -109,7 +84,6 @@
 
 #include "system.h"
 #include "assure.h"
-#include "c-ctype.h"
 #include "full-write.h"
 #include "quote.h"
 #include "readtokens.h"
@@ -124,56 +98,63 @@
   proper_name_lite ("Niels Moller", "Niels M\303\266ller")
 
 /* Token delimiters when reading from a file.  */
-#define DELIM "\n\t "
-
-/* __int128 is experimental; to use it, compile with -DUSE_INT128.  */
-#ifndef USE_INT128
-# define USE_INT128 false
-#endif
-
-/* Typedefs and macros related to an unsigned type that is no narrower
-   than 32 bits and no narrower than unsigned int.  For efficiency,
-   use the widest hardware-supported type.  */
-#if USE_INT128
-typedef unsigned __int128 wide_uint;
-typedef __int128 wide_int;
-# define W_TYPE_SIZE 128
+static char const DELIM[] = "\n\t ";
+
+/* GMP uses the unsigned integer type mp_limb_t as its word in
+   multiprecision arithmetic.  This code uses the same word for single
+   and double precision integer arithmetic.  Although previous
+   versions of this code used uintmax_t for single and double
+   precision, that introduced opportunities for bugs and was not worth
+   the hassle, as mp_limb_t and uintmax_t are invariably the same on
+   64-bit platforms, and 32-bit platforms are less important now.
+
+   Although GMP can be built with GMP_NUMB_BITS < GMP_LIMB_BITS,
+   so that some high-order bits of a word are not used, do not
+   do this in single and double precision integer arithmetic.
+   Instead, always use the full word.  */
+
+/* The word size in bits.
+   In the rest of this file's comments, B = 2^W_TYPE_SIZE is the base,
+   and the notation (a1,a0) stands for B*a1 + a0.  */
+#ifdef GMP_LIMB_BITS
+# define W_TYPE_SIZE GMP_LIMB_BITS
 #else
-typedef uintmax_t wide_uint;
-typedef intmax_t wide_int;
-# define W_TYPE_SIZE UINTMAX_WIDTH
+/* An older GMP, or mini-gmp; guess the usual value.  */
+# define W_TYPE_SIZE ULONG_WIDTH
 #endif
-#define WIDE_UINT_MAX ((wide_uint) -1)
+
+/* The maximum value of a word.  */
+#define MP_LIMB_MAX ((mp_limb_t) -1)
+
+/* Check W_TYPE_SIZE's value, as it might be a guess.  */
+static_assert (MP_LIMB_MAX >> (W_TYPE_SIZE - 1) == 1);
+
+/* Check that the builder didn't specify something perverse like
+   "-DMINI_GMP_LIMB_TYPE=short -DW_TYPE_SIZE=USHRT_WIDTH".
+   This could result in undefined behavior due to signed integer
+   overflow if a word promotes to int.  */
+static_assert (INT_MAX < MP_LIMB_MAX);
 
 #ifndef USE_LONGLONG_H
-/* With the way we use longlong.h, it's only safe to use
-   when UWtype = UHWtype, as there were various cases
-   (as can be seen in the history for longlong.h) where
-   for example, _LP64 was required to enable W_TYPE_SIZE==64 code,
-   to avoid compile time or run time issues.  */
-# define USE_LONGLONG_H (W_TYPE_SIZE == ULONG_WIDTH)
+# if (defined INT32_MAX && defined UINT32_MAX \
+      && defined INT64_MAX && defined UINT64_MAX)
+#  define USE_LONGLONG_H true
+# else
+#  define USE_LONGLONG_H false
+# endif
 #endif
 
 #if USE_LONGLONG_H
 
 /* Make definitions for longlong.h to make it do what it can do for us */
 
-# define UWtype  wide_uint
-# define UHWtype unsigned long int
+typedef mp_limb_t UWtype;
+typedef unsigned int UHWtype;
 # undef UDWtype
-# if HAVE_ATTRIBUTE_MODE
-typedef unsigned int UQItype    __attribute__ ((mode (QI)));
-typedef          int SItype     __attribute__ ((mode (SI)));
-typedef unsigned int USItype    __attribute__ ((mode (SI)));
-typedef          int DItype     __attribute__ ((mode (DI)));
-typedef unsigned int UDItype    __attribute__ ((mode (DI)));
-# else
-typedef unsigned char UQItype;
-typedef          long SItype;
-typedef unsigned long int USItype;
-typedef long long int DItype;
-typedef unsigned long long int UDItype;
-# endif
+typedef int32_t SItype;
+typedef uint32_t USItype;
+typedef int64_t DItype;
+typedef uint64_t UDItype;
 # define LONGLONG_STANDALONE     /* Don't require GMP's longlong.h mdep files */
 
 /* FIXME make longlong.h really standalone, so that ASSERT, __GMP_DECLSPEC
@@ -184,31 +165,31 @@ typedef unsigned long long int UDItype;
 #  if defined __GNUC__ && defined __GNUC_MINOR__
 #   define __GMP_GNUC_PREREQ(a, b) ((a) < __GNUC__ + ((b) <= __GNUC_MINOR__))
 #  else
-#   define __GMP_GNUC_PREREQ(a, b) 0
+#   define __GMP_GNUC_PREREQ(a, b) false
 #  endif
 # endif
 
 /* longlong.h uses these macros only in certain system compiler combinations.
    Ensure usage to pacify -Wunused-macros.  */
-# if (defined ASSERT || defined UHWtype \
-      || defined __GMP_DECLSPEC || defined __GMP_GNUC_PREREQ)
+# if defined ASSERT || defined __GMP_DECLSPEC || defined __GMP_GNUC_PREREQ
 # endif
 
 # if _ARCH_PPC
 #  define HAVE_HOST_CPU_FAMILY_powerpc 1
 # endif
 # include "longlong.h"
-
-#else /* not USE_LONGLONG_H */
-
-# define __ll_B ((wide_uint) 1 << (W_TYPE_SIZE / 2))
-# define __ll_lowpart(t)  ((wide_uint) (t) & (__ll_B - 1))
-# define __ll_highpart(t) ((wide_uint) (t) >> (W_TYPE_SIZE / 2))
-
 #endif
 
-/* 2*3*5*7*11...*101 is 128 bits, and has 26 prime factors */
-#define MAX_NFACTS 26
+/* 2*3*5*7*11...*101 fits in 128 bits, and has 26 prime factors.
+   This code can be extended in the future as needed; show as an example
+   2*3*5*7*11...*193 which fits in 257 bits, and has 44 prime factors.  */
+#if 2 * W_TYPE_SIZE <= 128
+enum { MAX_NFACTS = 26 };
+#elif 2 * W_TYPE_SIZE <= 257
+enum { MAX_NFACTS = 44 };
+#else
+# error "configuration has a wide word; please update MAX_NFACTS definition"
+#endif
 
 enum
 {
@@ -227,79 +208,129 @@ static struct option const long_options[] =
 /* If true, use p^e output format.  */
 static bool print_exponents;
 
-/* This represents an unsigned integer twice as wide as wide_uint.  */
-typedef struct { wide_uint uu[2]; } uuint;
+/* This represents an unsigned integer twice as wide as a word.  */
+typedef struct { mp_limb_t uu[2]; } uuint;
 
 /* Accessors and constructors for the type.  Programs should not
    access the type's internals directly, in case some future version
    replaces the type with unsigned __int256 or whatever.  */
-static wide_uint lo (uuint u) { return u.uu[0]; }
-static wide_uint hi (uuint u) { return u.uu[1]; }
-static void hiset (uuint *u, wide_uint hi) { u->uu[1] = hi; }
+static mp_limb_t lo (uuint u) { return u.uu[0]; }
+static mp_limb_t hi (uuint u) { return u.uu[1]; }
+static void hiset (uuint *u, mp_limb_t hi) { u->uu[1] = hi; }
+static bool hi_is_set (uuint const *pu) { return pu->uu[1] != 0; }
 static void
-uuset (wide_uint *phi, wide_uint *plo, uuint uu)
+uuset (mp_limb_t *phi, mp_limb_t *plo, uuint uu)
 {
   *phi = hi (uu);
   *plo = lo (uu);
 }
 static uuint
-make_uuint (wide_uint hi, wide_uint lo)
+make_uuint (mp_limb_t hi, mp_limb_t lo)
 {
   return (uuint) {{lo, hi}};
 }
 
-/* BIG_POWER_OF_10 is a positive power of 10 that does not exceed WIDE_UINT_MAX.
+/* BIG_POWER_OF_10 is a positive power of 10 that fits in a word.
    The larger it is, the more efficient the code will likely be.
    LOG_BIG_POWER_OF_10 = log (BIG_POWER_OF_10).  */
-#if W_TYPE_SIZE < 64
-# error "platform does not support 64-bit integers"
+#if W_TYPE_SIZE < 4
+# error "Configuration error; platform word is impossibly narrow"
+#elif W_TYPE_SIZE < 30
+/* An unusually narrow word.  */
+static mp_limb_t const BIG_POWER_OF_10 = 10;
+enum { LOG_BIG_POWER_OF_10 = 1 };
+#elif W_TYPE_SIZE < 64
+/* Almost surely a 32-bit word.  */
+static mp_limb_t const BIG_POWER_OF_10 = 1000000000;
+enum { LOG_BIG_POWER_OF_10 = 9 };
 #elif W_TYPE_SIZE < 128
-/* A mainstream platforms, with at-least-64-bit uintmax_t.  */
-static wide_uint const BIG_POWER_OF_10 = 10000000000000000000llu;
+/* Almost surely a 64-bit word.  */
+static mp_limb_t const BIG_POWER_OF_10 = 10000000000000000000llu;
 enum { LOG_BIG_POWER_OF_10 = 19 };
 #else
-/* If built with -DUSE_INT128, a platform that supports __int128; otherwise,
-   a so-far-only-theoretical platform with at-least-128-bit uintmax_t.
-   This is for performance; the 64-bit mainstream code will still work.  */
-static wide_uint const BIG_POWER_OF_10 =
-  (wide_uint) 10000000000000000000llu * 10000000000000000000llu;
+/* An unusually wide word.  */
+static mp_limb_t const BIG_POWER_OF_10 =
+  (mp_limb_t) 10000000000000000000llu * 10000000000000000000llu;
 enum { LOG_BIG_POWER_OF_10 = 38 };
 #endif
 
+/* Check that struct factors can use unsigned char to record a uuint's
+   prime factor's multiplicity, which is at most 2 * W_TYPE_SIZE - 1.  */
+static_assert (2 * W_TYPE_SIZE - 1 <= UCHAR_MAX);
+
+/* Likewise for recording the number of prime factors.  */
+static_assert (MAX_NFACTS <= UCHAR_MAX);
+
+/* Prime factors of a uuint.  At most one is a uuint.  */
 struct factors
 {
-  uuint plarge; /* Can have a single large factor */
-  wide_uint     p[MAX_NFACTS];
+  /* If PLARGE.uu[1] is nonzero, PLARGE is a double-word prime factor
+     with multiplicity 1; otherwise, PLARGE is not a factor and
+     PLARGE.uu[0] might not be initialized.  */
+  uuint plarge;
+
+  /* Distinct single-word prime factors, their multiplicities,
+     and their number.  */
+  mp_limb_t     p[MAX_NFACTS];
   unsigned char e[MAX_NFACTS];
   unsigned char nfactors;
 };
 
+/* An mpz's prime factor and multiplicity.  */
+struct mp_factor
+{
+  mpz_t p;
+  mp_bitcnt_t e;
+};
+
+/* Prime factors of an mpz_t.  */
 struct mp_factors
 {
-  mpz_t             *p;
-  unsigned long int *e;
+  /* A vector of distinct prime factors, a count of the factors,
+     and the number of allocated slots in the vector.  */
+  struct mp_factor *f;
   idx_t nfactors;
   idx_t nalloc;
 };
 
-static void factor (wide_uint, wide_uint, struct factors *);
+static void factor (struct factors *, mp_limb_t, mp_limb_t);
+static void factor_up (struct factors *, mp_limb_t, mp_limb_t, idx_t);
 
+/* Set (w1,w0) = u * v.  */
+#ifndef umul_ppmm
+/* Speed things up if there is an unsigned type uuroom_t that is wide
+   enough to hold two words.  */
+# if W_TYPE_SIZE <= UINTMAX_WIDTH / 2
+#  define uuroom_t uintmax_t
+# elif W_TYPE_SIZE <= 64 && defined __SIZEOF_INT128__
+#  define uuroom_t unsigned __int128
+# endif
+# ifdef uuroom_t
+#  define umul_ppmm(w1, w0, u, v)		\
+    do {					\
+      uuroom_t _u = u, _w = _u * (v);		\
+      (w1) = _w >> W_TYPE_SIZE;			\
+      (w0) = _w;				\
+    } while (false)
+# endif
+#endif
 #ifndef umul_ppmm
+static mp_limb_t const __ll_B = (mp_limb_t) 1 << (W_TYPE_SIZE / 2);
+static mp_limb_t __ll_lowpart (mp_limb_t t) { return t & (__ll_B - 1); }
+static mp_limb_t __ll_highpart (mp_limb_t t) { return t >> (W_TYPE_SIZE / 2); }
 # define umul_ppmm(w1, w0, u, v)                                        \
   do {                                                                  \
-    wide_uint __x0, __x1, __x2, __x3;                                   \
-    unsigned long int __ul, __vl, __uh, __vh;                           \
-    wide_uint __u = (u), __v = (v);                                     \
+    mp_limb_t __u = u, __v = v,						\
                                                                         \
-    __ul = __ll_lowpart (__u);                                          \
-    __uh = __ll_highpart (__u);                                         \
-    __vl = __ll_lowpart (__v);                                          \
-    __vh = __ll_highpart (__v);                                         \
+      __ul = __ll_lowpart (__u),					\
+      __uh = __ll_highpart (__u),					\
+      __vl = __ll_lowpart (__v),					\
+      __vh = __ll_highpart (__v),					\
                                                                         \
-    __x0 = (wide_uint) __ul * __vl;                                     \
-    __x1 = (wide_uint) __ul * __vh;                                     \
-    __x2 = (wide_uint) __uh * __vl;                                     \
-    __x3 = (wide_uint) __uh * __vh;                                     \
+      __x0 = __ul * __vl,						\
+      __x1 = __ul * __vh,						\
+      __x2 = __uh * __vl,						\
+      __x3 = __uh * __vh;						\
                                                                         \
     __x1 += __ll_highpart (__x0);/* This can't give carry.  */		\
     if (ckd_add (&__x1, __x1, __x2))	/* Did this give a carry?  */	\
@@ -307,9 +338,10 @@ static void factor (wide_uint, wide_uint, struct factors *);
                                                                         \
     (w1) = __x3 + __ll_highpart (__x1);                                 \
     (w0) = (__x1 << W_TYPE_SIZE / 2) + __ll_lowpart (__x0);             \
-  } while (0)
+  } while (false)
 #endif
 
+/* Set (q,r) to the quotient and remainder of dividing (n1,n0) by d.  */
 #if !defined udiv_qrnnd || defined UDIV_NEEDS_NORMALIZATION
 /* Define our own, not needing normalization.  This function is
    currently not performance critical, so keep it simple.  Similar to
@@ -317,12 +349,9 @@ static void factor (wide_uint, wide_uint, struct factors *);
 # undef udiv_qrnnd
 # define udiv_qrnnd(q, r, n1, n0, d)                                    \
   do {                                                                  \
-    wide_uint __d1, __d0, __q, __r1, __r0;                              \
+    mp_limb_t __d1 = d, __d0 = 0, __q = 0, __r1 = n1, __r0 = n0;	\
                                                                         \
-    __d1 = (d); __d0 = 0;                                               \
-    __r1 = (n1); __r0 = (n0);                                           \
     affirm (__r1 < __d1);                                               \
-    __q = 0;                                                            \
     for (int __i = W_TYPE_SIZE; __i > 0; __i--)                         \
       {                                                                 \
         rsh2 (__d1, __d0, __d1, __d0, 1);                               \
@@ -335,16 +364,13 @@ static void factor (wide_uint, wide_uint, struct factors *);
       }                                                                 \
     (r) = __r0;                                                         \
     (q) = __q;                                                          \
-  } while (0)
+  } while (false)
 #endif
 
+/* Set (sh,sl) = (ah,al) + (bh,bl).  Overflow wraps around.  */
 #if !defined add_ssaaaa
 # define add_ssaaaa(sh, sl, ah, al, bh, bl)                             \
-  do {                                                                  \
-    wide_uint _add_x;							\
-    (sh) = (ah) + (bh) + ckd_add (&_add_x, al, bl);			\
-    (sl) = _add_x;                                                      \
-  } while (0)
+   ((sh) = (ah) + (bh) + ckd_add (&(sl), al, bl))
 #endif
 
 /* Set (rh,rl) = (ah,al) >> cnt, where 0 < cnt < W_TYPE_SIZE.  */
@@ -352,62 +378,77 @@ static void factor (wide_uint, wide_uint, struct factors *);
   do {                                                                  \
     (rl) = ((ah) << (W_TYPE_SIZE - (cnt))) | ((al) >> (cnt));           \
     (rh) = (ah) >> (cnt);                                               \
-  } while (0)
+  } while (false)
 
-/* Set (rh,rl) = (ah,al) << cnt, where 0 < cnt < W_TYPE_SIZE.  */
+/* Set (rh,rl) = (ah,al) << cnt, where 0 < cnt < W_TYPE_SIZE.
+   Overflow wraps around.  */
 #define lsh2(rh, rl, ah, al, cnt)                                       \
   do {                                                                  \
     (rh) = ((ah) << cnt) | ((al) >> (W_TYPE_SIZE - (cnt)));             \
     (rl) = (al) << (cnt);                                               \
-  } while (0)
+  } while (false)
 
-#define ge2(ah, al, bh, bl)                                             \
-  ((ah) > (bh) || ((ah) == (bh) && (al) >= (bl)))
+/* (ah,hl) < (bh,bl)?  */
+static bool
+lt2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl)
+{
+  mp_limb_t dh, dl;
+  bool vh = ckd_sub (&dh, ah, bh);
+  mp_limb_t vl = ckd_sub (&dl, al, bl);
+  return vh | ckd_sub (&dh, dh, vl);
+}
 
-#define gt2(ah, al, bh, bl)                                             \
-  ((ah) > (bh) || ((ah) == (bh) && (al) > (bl)))
+/* (ah,hl) >= (bh,bl)?  */
+static bool
+ge2 (mp_limb_t ah, mp_limb_t al, mp_limb_t bh, mp_limb_t bl)
+{
+  return !lt2 (ah, al, bh, bl);
+}
 
+/* Set (rh,rl) = (ah,al) - (bh,bl).  Overflow wraps around.  */
 #ifndef sub_ddmmss
 # define sub_ddmmss(rh, rl, ah, al, bh, bl)                             \
    ((rh) = (ah) - (bh) - ckd_sub (&(rl), al, bl))
 #endif
 
-/* Requires that a < n and b <= n */
+/* Set r = (a - b) mod n, where a < n & b <= n.  */
 #define submod(r,a,b,n)                                                 \
   do {                                                                  \
-    wide_uint _s, _t = -ckd_sub (&_s, a, b);				\
+    mp_limb_t _s, _t = -ckd_sub (&_s, a, b);				\
     (r) = ((n) & _t) + _s;						\
-  } while (0)
+  } while (false)
 
+/* Set r = (a + b) mod n, where a < n & b <= n.  */
 #define addmod(r,a,b,n)                                                 \
-  submod ((r), (a), ((n) - (b)), (n))
+  submod (r, a, (n) - (b), n)
 
 /* Modular two-word addition and subtraction.  For performance reasons, the
    most significant bit of n1 must be clear.  The destination variables must be
    distinct from the mod operand.  */
 #define addmod2(r1, r0, a1, a0, b1, b0, n1, n0)                         \
   do {                                                                  \
-    add_ssaaaa ((r1), (r0), (a1), (a0), (b1), (b0));                    \
-    if (ge2 ((r1), (r0), (n1), (n0)))                                   \
-      sub_ddmmss ((r1), (r0), (r1), (r0), (n1), (n0));                  \
-  } while (0)
+    add_ssaaaa (r1, r0, a1, a0, b1, b0);				\
+    if (ge2 (r1, r0, n1, n0))						\
+      sub_ddmmss (r1, r0, r1, r0, n1, n0);				\
+  } while (false)
 #define submod2(r1, r0, a1, a0, b1, b0, n1, n0)                         \
   do {                                                                  \
-    sub_ddmmss ((r1), (r0), (a1), (a0), (b1), (b0));                    \
-    if ((wide_int) (r1) < 0)                                            \
-      add_ssaaaa ((r1), (r0), (r1), (r0), (n1), (n0));                  \
-  } while (0)
-
-#define HIGHBIT_TO_MASK(x)                                              \
-  (((wide_int)-1 >> 1) < 0                                              \
-   ? (wide_uint)((wide_int)(x) >> (W_TYPE_SIZE - 1))                    \
-   : ((x) & ((wide_uint) 1 << (W_TYPE_SIZE - 1))                        \
-      ? WIDE_UINT_MAX : (wide_uint) 0))
-
-/* Return r = a mod d, where a = <a1,a0>, d = <d1,d0>.
-   Requires that d1 != 0.  */
+    bool _v1 = ckd_sub (&(r1), a1, b1);					\
+    mp_limb_t _v0 = ckd_sub (&(r0), a0, b0);				\
+    if (_v1 | ckd_sub (&(r1), r1, _v0))					\
+      add_ssaaaa (r1, r0, r1, r0, n1, n0);				\
+  } while (false)
+
+/* Return 0 if x < B/2, MP_LIMB_MAX otherwise.  */
+static mp_limb_t
+highbit_to_mask (mp_limb_t x)
+{
+  return - (x >> (W_TYPE_SIZE - 1));
+}
+
+/* Return (a1,a0) mod (d1,d0), where d1 != 0.  */
 ATTRIBUTE_PURE static uuint
-mod2 (wide_uint a1, wide_uint a0, wide_uint d1, wide_uint d0)
+mod2 (mp_limb_t a1, mp_limb_t a0, mp_limb_t d1, mp_limb_t d0)
 {
   affirm (d1 != 0);
 
@@ -431,16 +472,12 @@ mod2 (wide_uint a1, wide_uint a0, wide_uint d1, wide_uint d0)
   return make_uuint (a1, a0);
 }
 
+/* Return the greatest common divisor of a and b,
+   where b is odd.  */
 ATTRIBUTE_CONST
-static wide_uint
-gcd_odd (wide_uint a, wide_uint b)
+static mp_limb_t
+gcd_odd (mp_limb_t a, mp_limb_t b)
 {
-  if ((b & 1) == 0)
-    {
-      wide_uint t = b;
-      b = a;
-      a = t;
-    }
   if (a == 0)
     return b;
 
@@ -449,18 +486,17 @@ gcd_odd (wide_uint a, wide_uint b)
 
   for (;;)
     {
-      wide_uint t;
-      wide_uint bgta;
+      mp_limb_t t;
+      mp_limb_t bgta;
 
       assume (a);
-      a >>= stdc_trailing_zeros (a);
-      a >>= 1;
+      mp_limb_t ao = a >> stdc_trailing_zeros (a);
 
-      t = a - b;
+      t = (ao >> 1) - b;
       if (t == 0)
-        return (a << 1) + 1;
+        return ao;
 
-      bgta = HIGHBIT_TO_MASK (t);
+      bgta = highbit_to_mask (t);
 
       /* b <-- min (a, b) */
       b += (bgta & t);
@@ -470,94 +506,100 @@ gcd_odd (wide_uint a, wide_uint b)
     }
 }
 
+/* Return the greatest common divisor of (a1,a0) and (b1,b0),
+   where (b1,b0) is odd.  */
 ATTRIBUTE_PURE static uuint
-gcd2_odd (wide_uint a1, wide_uint a0, wide_uint b1, wide_uint b0)
+gcd2_odd (mp_limb_t a1, mp_limb_t a0, mp_limb_t b1, mp_limb_t b0)
 {
   affirm (b0 & 1);
 
-  if ((a0 | a1) == 0)
-    return make_uuint (b1, b0);
   if (!a0)
-    a0 = a1, a1 = 0;
-  assume (a0);
-  int ctz = stdc_trailing_zeros (a0);
-  if (ctz)
-    rsh2 (a1, a0, a1, a0, ctz);
+    {
+      a0 = a1, a1 = 0;
+      if (!a0)
+        return make_uuint (b1, b0);
+    }
+  goto make_A_odd;
 
   for (;;)
     {
       if ((b1 | a1) == 0)
         return make_uuint (0, gcd_odd (b0, a0));
 
-      if (gt2 (a1, a0, b1, b0))
+      if (lt2 (b1, b0, a1, a0))
         {
           sub_ddmmss (a1, a0, a1, a0, b1, b0);
           if (!a0)
             a0 = a1, a1 = 0;
+        make_A_odd:
           assume (a0);
-          ctz = stdc_trailing_zeros (a0);
+          int ctz = stdc_trailing_zeros (a0);
           if (ctz)
             rsh2 (a1, a0, a1, a0, ctz);
         }
-      else if (gt2 (b1, b0, a1, a0))
+      else
         {
           sub_ddmmss (b1, b0, b1, b0, a1, a0);
           if (!b0)
-            b0 = b1, b1 = 0;
-          assume (b0);
-          ctz = stdc_trailing_zeros (b0);
+            {
+              b0 = b1, b1 = 0;
+              if (!b0)
+                break;
+            }
+          int ctz = stdc_trailing_zeros (b0);
           if (ctz)
             rsh2 (b1, b0, b1, b0, ctz);
         }
-      else
-        break;
     }
 
   return make_uuint (a1, a0);
 }
 
+/* Store into FACTORS a prime factor PRIME with multiplicity M.  */
 static void
 factor_insert_multiplicity (struct factors *factors,
-                            wide_uint prime, int m)
+                            mp_limb_t prime, int m)
 {
   int nfactors = factors->nfactors;
-  wide_uint *p = factors->p;
+  mp_limb_t *p = factors->p;
   unsigned char *e = factors->e;
 
   /* Locate position for insert new or increment e.  */
   int i;
-  for (i = nfactors - 1; i >= 0; i--)
+  for (i = nfactors; 0 < i; i--)
     {
-      if (p[i] <= prime)
+      if (p[i - 1] < prime)
         break;
-    }
-
-  if (i < 0 || p[i] != prime)
-    {
-      for (int j = nfactors - 1; j > i; j--)
+      if (p[i - 1] == prime)
         {
-          p[j + 1] = p[j];
-          e[j + 1] = e[j];
+          e[i - 1] += m;
+          return;
         }
-      p[i + 1] = prime;
-      e[i + 1] = m;
-      factors->nfactors = nfactors + 1;
-    }
-  else
-    {
-      e[i] += m;
     }
+
+  factors->nfactors = nfactors + 1;
+  memmove (&p[i + 1], &p[i], (nfactors - i) * sizeof *p);
+  memmove (&e[i + 1], &e[i], (nfactors - i) * sizeof *e);
+  e[i] = m;
+  p[i] = prime;
 }
 
-#define factor_insert(f, p) factor_insert_multiplicity (f, p, 1)
+/* Store into FACTORS a prime factor PRIME.  */
+static void
+factor_insert (struct factors *factors, mp_limb_t prime)
+{
+  factor_insert_multiplicity (factors, prime, 1);
+}
 
+/* Store into FACTORS a prime factor (P1,P0).  If P1 != 0,
+   FACTORS must not already contain a large prime factor.  */
 static void
 factor_insert_large (struct factors *factors,
-                     wide_uint p1, wide_uint p0)
+                     mp_limb_t p1, mp_limb_t p0)
 {
   if (p1 > 0)
     {
-      affirm (hi (factors->plarge) == 0);
+      affirm (!hi_is_set (&factors->plarge));
       factors->plarge = make_uuint (p1, p0);
     }
   else
@@ -586,109 +628,118 @@ mpz_va_init (void (*mpz_single_init)(mpz_t), ...)
 }
 #endif
 
-static void mp_factor (mpz_t, struct mp_factors *);
+#ifndef mpn_tdiv_qr
 
 static void
-mp_factor_init (struct mp_factors *factors)
+copy_mpn_from_mpz (mp_limb_t *p, mp_size_t n, mpz_t z)
 {
-  factors->p = nullptr;
-  factors->e = nullptr;
-  factors->nfactors = 0;
-  factors->nalloc = 0;
+  mp_size_t zsize = mpz_size (z);
+  mpn_copyi (p, mpz_limbs_read (z), zsize);
+  mpn_zero (p + zsize, n - zsize);
+}
+
+static void
+mpn_tdiv_qr (mp_limb_t *qp, mp_limb_t *rp, MAYBE_UNUSED mp_size_t qxn,
+             mp_limb_t const *np, mp_size_t nn,
+             mp_limb_t const *dp, mp_size_t dn)
+{
+  mpz_t q, r, n, d;
+  mpz_inits (q, r, nullptr);
+  mpz_tdiv_qr (q, r, mpz_roinit_n (n, np, nn), mpz_roinit_n (d, dp, dn));
+  copy_mpn_from_mpz (qp, nn - dn + 1, q);
+  copy_mpn_from_mpz (rp, dn, r);
+  mpz_clears (q, r, nullptr);
+}
+#endif
+
+static struct mp_factors mp_factor (mpz_t);
+
+/* Return an empty set of factors.  */
+static struct mp_factors
+mp_no_factors (void)
+{
+  return (struct mp_factors) {0,};
 }
 
+/* Free storage allocated for FACTORS, making it uninitialized.  */
 static void
 mp_factor_clear (struct mp_factors *factors)
 {
+  struct mp_factor *f = factors->f;
   for (idx_t i = 0; i < factors->nfactors; i++)
-    mpz_clear (factors->p[i]);
-
-  free (factors->p);
-  free (factors->e);
+    mpz_clear (f[i].p);
+  free (f);
 }
 
+/* Store into FACTORS a prime factor PRIME with multiplicity M.  */
 static void
-mp_factor_insert (struct mp_factors *factors, mpz_t prime)
+mp_factor_insert (struct mp_factors *factors, mpz_t prime, mp_bitcnt_t m)
 {
   idx_t nfactors = factors->nfactors;
-  mpz_t *p = factors->p;
-  unsigned long int *e = factors->e;
-  ptrdiff_t i;
+  struct mp_factor *f = factors->f;
+  idx_t i;
 
   /* Locate position for insert new or increment e.  */
-  for (i = nfactors - 1; i >= 0; i--)
+  for (i = nfactors; 0 < i; i--)
     {
-      if (mpz_cmp (p[i], prime) <= 0)
+      int sgn = mpz_cmp (f[i - 1].p, prime);
+      if (sgn < 0)
         break;
-    }
-
-  if (i < 0 || mpz_cmp (p[i], prime) != 0)
-    {
-      if (factors->nfactors == factors->nalloc)
+      if (sgn == 0)
         {
-          p = xpalloc (p, &factors->nalloc, 1, -1, sizeof *p);
-          e = xireallocarray (e, factors->nalloc, sizeof *e);
+          f[i - 1].e += m;
+          return;
         }
+    }
 
-      mpz_init (p[nfactors]);
-      for (ptrdiff_t j = nfactors - 1; j > i; j--)
-        {
-          mpz_set (p[j + 1], p[j]);
-          e[j + 1] = e[j];
-        }
-      mpz_set (p[i + 1], prime);
-      e[i + 1] = 1;
+  if (nfactors == factors->nalloc)
+    factors->f = f = xpalloc (f, &factors->nalloc, 1, -1, sizeof *f);
 
-      factors->p = p;
-      factors->e = e;
-      factors->nfactors = nfactors + 1;
-    }
-  else
-    {
-      e[i] += 1;
-    }
+  factors->nfactors = nfactors + 1;
+  memmove (&f[i + 1], &f[i], (nfactors - i) * sizeof *f);
+  f[i].e = m;
+  mpz_init_set (f[i].p, prime);
 }
 
+/* Store into FACTORS a prime factor PRIME with multiplicity M.  */
 static void
-mp_factor_insert_ui (struct mp_factors *factors, unsigned long int prime)
+mp_factor_insert1 (struct mp_factors *factors, mp_limb_t prime, mp_bitcnt_t m)
 {
-  mpz_t pz;
-
-  mpz_init_set_ui (pz, prime);
-  mp_factor_insert (factors, pz);
-  mpz_clear (pz);
+  mpz_t pz = MPZ_ROINIT_N (&prime, 1);
+  mp_factor_insert (factors, pz, m);
 }

-#define P(a,b,c,d) a,
-static const unsigned char primes_diff[] = {
+/* primes_ptab[i] is prime i, where the zeroth prime is 3.
+   However, the last 7 entries are sentinels.  */
+#define P(a,b,c) a,
+static int_least16_t const primes_ptab[] = {
 #include "primes.h"
 0,0,0,0,0,0,0                           /* 7 sentinels for 8-way loop */
 };
 #undef P
 
-#define PRIMES_PTAB_ENTRIES (ARRAY_CARDINALITY (primes_diff) - 8 + 1)
-
-#define P(a,b,c,d) b,
-static const unsigned char primes_diff8[] = {
-#include "primes.h"
-0,0,0,0,0,0,0                           /* 7 sentinels for 8-way loop */
-};
-#undef P
+enum { PRIMES_PTAB_ENTRIES = ARRAY_CARDINALITY (primes_ptab) - 8 + 1 };
 
 struct primes_dtab
 {
-  wide_uint binv, lim;
+  mp_limb_t binv, lim;
 };
 
-#define P(a,b,c,d) {c,d},
+/* primes_dtab[i].binv is the multiplicative inverse of prime i
+   modulo B, where the zeroth prime is 3.  That is,
+   ((primes_dtab[i].binv) * (prime i)) mod B = 1.
+   primes_dtab[i].lim is the maximum value that won't overflow in
+   mp_limb_t arithmetic when multiplied by the ith prime.
+   However, the last 7 entries of primes_dtab are sentinels.  */
+#define P(a,b,c) {b,c},
 static const struct primes_dtab primes_dtab[] = {
 #include "primes.h"
 {1,0},{1,0},{1,0},{1,0},{1,0},{1,0},{1,0} /* 7 sentinels for 8-way loop */
 };
 #undef P
 
-/* Verify that wide_uint is not wider than
+/* Verify that a word is not wider than
    the integers used to generate primes.h.  */
 static_assert (W_TYPE_SIZE <= WIDE_UINT_BITS);
 
@@ -696,24 +747,30 @@ static_assert (W_TYPE_SIZE <= WIDE_UINT_BITS);
    This flag is used only in the GMP code.  */
 static bool dev_debug = false;
 
-/* Prove primality or run probabilistic tests.  */
-static bool flag_prove_primality = PROVE_PRIMALITY;
+/* Number of Miller-Rabin tests to run.  For more, see:
 
-/* Number of Miller-Rabin tests to run when not proving primality.  */
-#define MR_REPS 25
+     Ishmukhametov ST, Mubarakov BG, Rubtsova RG.
+     On the Number of Witnesses in the Miller-Rabin Primality Test.
+     Symmetry. 2020;12(6):890. https://doi.org/10.3390/sym12060890
 
-static void
-factor_insert_refind (struct factors *factors, wide_uint p, int i, int off)
-{
-  for (int j = 0; j < off; j++)
-    p += primes_diff[i + j];
-  factor_insert (factors, p);
-}
+   Its Corollary 1 suggests that the probability of error on random inputs
+   is less than 16^-MR_REPS, an improvement on the 4^-MR_REPS commonly cited.
+   If MR_REPS is 24 this means the probability of error is less than 1.26e-29,
+   which is much less than the likelihood of hardware error and so can
+   be treated as essentially zero.
+
+   For adversarial inputs, no known false positives exist for Baillie-PSW,
+   which mpz_probab_prime_p always uses.  So default MR_REPS to 24,
+   the maximum value for which mpz_probab_prime_p does not do extra
+   Miller-Rabin tests.  */
+#ifndef MR_REPS
+# define MR_REPS 24
+#endif
 
 /* Trial division with odd primes uses the following trick.
 
-   Let p be an odd prime, and B = 2^{W_TYPE_SIZE}.  For simplicity,
-   consider the case t < B (this is the second loop below).
+   Let p be an odd prime.  For simplicity, consider the case t < B;
+   this is the second loop below.
 
    From our tables we get
 
@@ -743,9 +800,30 @@ factor_insert_refind (struct factors *factors, wide_uint p, int i, int off)
    order, and the non-multiples of p onto the range lim < q < B.
  */
 
+/* The kernel of factor_using_division.  This function is called in a
+   loop unrolled by hand.  */
+static inline mp_limb_t
+divblock (struct factors *factors, mp_limb_t t0, struct primes_dtab const *pd,
+          idx_t i, int ioff)
+{
+  for (;;)
+    {
+      mp_limb_t q = t0 * pd[ioff].binv;
+      if (LIKELY (pd[ioff].lim < q))
+        return t0;
+      t0 = q;
+      factor_insert (factors, primes_ptab[i + ioff]);
+    }
+}
+
+/* Insert into FACTORS the factors of (T1,T0) found via trial division.
+   The candidate factors are 2 and the primes in the primes table.
+   However, primes less than prime I have
+   already been considered, and need not be looked at again.
+   Return (T1,T0) divided by the factors found.  */
 static uuint
-factor_using_division (wide_uint t1, wide_uint t0,
-                       struct factors *factors)
+factor_using_division (struct factors *factors, mp_limb_t t1, mp_limb_t t0,
+                       idx_t i)
 {
   if (t0 % 2 == 0)
     {
@@ -768,14 +846,13 @@ factor_using_division (wide_uint t1, wide_uint t0,
       factor_insert_multiplicity (factors, 2, cnt);
     }
 
-  wide_uint p = 3;
-  idx_t i;
-  for (i = 0; t1 > 0 && i < PRIMES_PTAB_ENTRIES; i++)
+  for (; t1 > 0 && i < PRIMES_PTAB_ENTRIES; i++)
     {
+      mp_limb_t p = primes_ptab[i];
       for (;;)
         {
-          wide_uint q1, q0, hi;
-          MAYBE_UNUSED wide_uint lo;
+          mp_limb_t q1, q0, hi;
+          MAYBE_UNUSED mp_limb_t lo;
 
           q0 = t0 * primes_dtab[i].binv;
           umul_ppmm (hi, lo, q0, p);
@@ -788,35 +865,21 @@ factor_using_division (wide_uint t1, wide_uint t0,
           t1 = q1; t0 = q0;
           factor_insert (factors, p);
         }
-      p += primes_diff[i + 1];
     }
 
-#define DIVBLOCK(I)                                                     \
-  do {                                                                  \
-    for (;;)                                                            \
-      {                                                                 \
-        q = t0 * pd[I].binv;                                            \
-        if (LIKELY (q > pd[I].lim))                                     \
-          break;                                                        \
-        t0 = q;                                                         \
-        factor_insert_refind (factors, p, i + 1, I);                    \
-      }                                                                 \
-  } while (0)
-
   for (; i < PRIMES_PTAB_ENTRIES; i += 8)
     {
-      wide_uint q;
       const struct primes_dtab *pd = &primes_dtab[i];
-      DIVBLOCK (0);
-      DIVBLOCK (1);
-      DIVBLOCK (2);
-      DIVBLOCK (3);
-      DIVBLOCK (4);
-      DIVBLOCK (5);
-      DIVBLOCK (6);
-      DIVBLOCK (7);
-
-      p += primes_diff8[i];
+      t0 = divblock (factors, t0, pd, i, 0);
+      t0 = divblock (factors, t0, pd, i, 1);
+      t0 = divblock (factors, t0, pd, i, 2);
+      t0 = divblock (factors, t0, pd, i, 3);
+      t0 = divblock (factors, t0, pd, i, 4);
+      t0 = divblock (factors, t0, pd, i, 5);
+      t0 = divblock (factors, t0, pd, i, 6);
+      t0 = divblock (factors, t0, pd, i, 7);
+
+      int_least32_t p = primes_ptab[i + 8];
       if (p * p > t0)
         break;
     }
@@ -824,45 +887,100 @@ factor_using_division (wide_uint t1, wide_uint t0,
   return make_uuint (t1, t0);
 }
 
-static void
-mp_factor_using_division (mpz_t t, struct mp_factors *factors)
+/* Return the number of limbs in positive N.  */
+static mp_size_t
+mp_size (mpz_t n)
 {
-  mpz_t q;
-  mp_bitcnt_t p;
+  /* Tell the compiler that N is positive; this can speed up access to N.  */
+  assume (0 < mpz_sgn (n));
 
-  devmsg ("[trial division] ");
+  return mpz_size (n);
+}
 
-  mpz_init (q);
+/* Insert into MP_FACTORS the factors of N if N < B^2 / 2,
+   and return true.  Otherwise, return false.
+   Primes less than prime PRIME_IDX have
+   already been considered, and need not be looked at again.  */
+static bool
+mp_finish_up_in_single (struct mp_factors *mp_factors, mpz_t n,
+                        idx_t prime_idx)
+{
+  if (2 < mp_size (n))
+    return false;
+  mp_limb_t n1 = mpz_getlimbn (n, 1);
+  if (n1 >> (W_TYPE_SIZE - 1))
+    return false;
+  mp_limb_t n0 = mpz_getlimbn (n, 0);
+  mpz_set_ui (n, 1);
 
-  p = mpz_scan1 (t, 0);
-  mpz_fdiv_q_2exp (t, t, p);
-  while (p)
+  struct factors factors;
+  factor_up (&factors, n1, n0, prime_idx);
+
+  if (hi_is_set (&factors.plarge))
     {
-      mp_factor_insert_ui (factors, 2);
-      --p;
+      mpz_t p = MPZ_ROINIT_N (factors.plarge.uu, 2);
+      mp_factor_insert (mp_factors, p, 1);
     }
 
-  unsigned long int d = 3;
-  for (idx_t i = 1; i <= PRIMES_PTAB_ENTRIES;)
+  for (int i = factors.nfactors; 0 < i; i--)
     {
-      if (! mpz_divisible_ui_p (t, d))
-        {
-          d += primes_diff[i++];
-          if (mpz_cmp_ui (t, d * d) < 0)
-            break;
-        }
-      else
+      mpz_t p = MPZ_ROINIT_N (&factors.p[i - 1], 1);
+      mp_factor_insert (mp_factors, p, factors.e[i - 1]);
+    }
+
+  return true;
+}
+
+/* Insert into MP_FACTORS the factors of N if N < B^2 / 2, and
+   return true.  Otherwise, return false.  N must be odd.  */
+static bool
+mp_finish_in_single (struct mp_factors *mp_factors, mpz_t n)
+{
+  return mp_finish_up_in_single (mp_factors, n, 0);
+}
+
+/* Return some or all factors of T.
+   Divide T by the factors found.  */
+static struct mp_factors
+mp_factor_using_division (mpz_t t)
+{
+  devmsg ("[trial division] ");
+  struct mp_factors factors = mp_no_factors ();
+
+  mp_bitcnt_t m = mpz_scan1 (t, 0);
+  if (m)
+    {
+      mpz_fdiv_q_2exp (t, t, m);
+      mp_factor_insert1 (&factors, 2, m);
+      if (mp_finish_in_single (&factors, t))
+        return factors;
+    }
+
+  for (idx_t i = 0; i < PRIMES_PTAB_ENTRIES; i++)
+    {
+      mp_limb_t d = primes_ptab[i];
+      for (m = 0; mpz_divisible_ui_p (t, d); m++)
         {
-          mpz_tdiv_q_ui (t, t, d);
-          mp_factor_insert_ui (factors, d);
+          mpz_divexact_ui (t, t, d);
+          if (mp_finish_up_in_single (&factors, t, i))
+            {
+              mp_factor_insert1 (&factors, d, m + 1);
+              return factors;
+            }
         }
+      if (m)
+        mp_factor_insert1 (&factors, d, m);
+      static_assert (SQUARE_OF_FIRST_OMITTED_PRIME
+                     <= MIN (MP_LIMB_MAX, ULONG_MAX));
+      if (mpz_cmp_ui (t, d * d) < 0)
+        break;
     }
 
-  mpz_clear (q);
+  return factors;
 }
 
 /* Entry i contains (2i+1)^(-1) mod 2^8.  */
-static const unsigned char  binvert_table[128] =
+static const unsigned char binvert_table[128] =
 {
   0x01, 0xAB, 0xCD, 0xB7, 0x39, 0xA3, 0xC5, 0xEF,
   0xF1, 0x1B, 0x3D, 0xA7, 0x29, 0x13, 0x35, 0xDF,
@@ -882,39 +1000,29 @@ static const unsigned char  binvert_table[128] =
   0x11, 0x3B, 0x5D, 0xC7, 0x49, 0x33, 0x55, 0xFF
 };
 
-/* Compute n^(-1) mod B, using a Newton iteration.  */
-#define binv(inv,n)                                                     \
-  do {                                                                  \
-    wide_uint  __n = (n);                                               \
-    wide_uint  __inv;                                                   \
-                                                                        \
-    __inv = binvert_table[(__n / 2) & 0x7F]; /*  8 */                   \
-    if (W_TYPE_SIZE > 8)   __inv = 2 * __inv - __inv * __inv * __n;     \
-    if (W_TYPE_SIZE > 16)  __inv = 2 * __inv - __inv * __inv * __n;     \
-    if (W_TYPE_SIZE > 32)  __inv = 2 * __inv - __inv * __inv * __n;     \
-                                                                        \
-    if (W_TYPE_SIZE > 64)                                               \
-      {                                                                 \
-        int  __invbits = 64;                                            \
-        do {                                                            \
-          __inv = 2 * __inv - __inv * __inv * __n;                      \
-          __invbits *= 2;                                               \
-        } while (__invbits < W_TYPE_SIZE);                              \
-      }                                                                 \
-                                                                        \
-    (inv) = __inv;                                                      \
-  } while (0)
+/* Compute n^(-1) mod B.  n must be odd.  */
+static mp_limb_t
+binv_limb (mp_limb_t n)
+{
+  mp_limb_t inv = binvert_table[(n / 2) & 0x7F];
+  if ( 8 < W_TYPE_SIZE) inv = 2 * inv - inv * inv * n;
+  if (16 < W_TYPE_SIZE) inv = 2 * inv - inv * inv * n;
+  if (32 < W_TYPE_SIZE) inv = 2 * inv - inv * inv * n;
+  for (int invbits = 64; invbits < W_TYPE_SIZE; invbits *= 2)
+    inv = 2 * inv - inv * inv * n;
+  return inv;
+}
 
 /* q = u / d, assuming d|u.  */
 #define divexact_21(q1, q0, u1, u0, d)                                  \
   do {                                                                  \
-    wide_uint _di, _q0;                                                 \
-    binv (_di, (d));                                                    \
+    mp_limb_t _di, _q0;							\
+    _di = binv_limb (d);						\
     _q0 = (u0) * _di;                                                   \
     if ((u1) >= (d))                                                    \
       {                                                                 \
-        wide_uint _p1;                                                  \
-        MAYBE_UNUSED wide_int _p0;                                      \
+        mp_limb_t _p1;							\
+        MAYBE_UNUSED mp_limb_t _p0;					\
         umul_ppmm (_p1, _p0, _q0, d);                                   \
         (q1) = ((u1) - _p1) * _di;                                      \
         (q0) = _q0;                                                     \
@@ -924,19 +1032,19 @@ static const unsigned char  binvert_table[128] =
         (q0) = _q0;                                                     \
         (q1) = 0;                                                       \
       }                                                                 \
-  } while (0)
+  } while (false)
 
 /* x B (mod n).  */
 #define redcify(r_prim, r, n)                                           \
   do {                                                                  \
-    MAYBE_UNUSED wide_uint _redcify_q;					\
+    MAYBE_UNUSED mp_limb_t _redcify_q;					\
     udiv_qrnnd (_redcify_q, r_prim, r, 0, n);                           \
-  } while (0)
+  } while (false)
 
 /* x B^2 (mod n).  Requires x > 0, n1 < B/2.  */
 #define redcify2(r1, r0, x, n1, n0)                                     \
   do {                                                                  \
-    wide_uint _r1, _r0, _i;                                             \
+    mp_limb_t _r1, _r0, _i;						\
     if ((x) < (n1))                                                     \
       {                                                                 \
         _r1 = (x); _r0 = 0;                                             \
@@ -950,20 +1058,20 @@ static const unsigned char  binvert_table[128] =
     while (_i-- > 0)                                                    \
       {                                                                 \
         lsh2 (_r1, _r0, _r1, _r0, 1);                                   \
-        if (ge2 (_r1, _r0, (n1), (n0)))                                 \
-          sub_ddmmss (_r1, _r0, _r1, _r0, (n1), (n0));                  \
+        if (ge2 (_r1, _r0, n1, n0))					\
+          sub_ddmmss (_r1, _r0, _r1, _r0, n1, n0);			\
       }                                                                 \
     (r1) = _r1;                                                         \
     (r0) = _r0;                                                         \
-  } while (0)
+  } while (false)
 
 /* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B.
    Both a and b must be in redc form, the result will be in redc form too.  */
-static inline wide_uint
-mulredc (wide_uint a, wide_uint b, wide_uint m, wide_uint mi)
+static inline mp_limb_t
+mulredc (mp_limb_t a, mp_limb_t b, mp_limb_t m, mp_limb_t mi)
 {
-  wide_uint rh, rl, q, th, xh;
-  MAYBE_UNUSED wide_uint tl;
+  mp_limb_t rh, rl, q, th, xh;
+  MAYBE_UNUSED mp_limb_t tl;
 
   umul_ppmm (rh, rl, a, b);
   q = rl * mi;
@@ -977,17 +1085,17 @@ mulredc (wide_uint a, wide_uint b, wide_uint m, wide_uint mi)
 /* Modular two-word multiplication, r = a * b mod m, with mi = m^(-1) mod B.
    Both a and b must be in redc form, the result will be in redc form too.
    For performance reasons, the most significant bit of m must be clear.  */
-static wide_uint
-mulredc2 (wide_uint *r1p,
-          wide_uint a1, wide_uint a0, wide_uint b1, wide_uint b0,
-          wide_uint m1, wide_uint m0, wide_uint mi)
+static mp_limb_t
+mulredc2 (mp_limb_t *r1p,
+          mp_limb_t a1, mp_limb_t a0, mp_limb_t b1, mp_limb_t b0,
+          mp_limb_t m1, mp_limb_t m0, mp_limb_t mi)
 {
-  wide_uint r1, r0, q, p1, t1, t0, s1, s0;
-  MAYBE_UNUSED wide_uint p0;
+  mp_limb_t r1, r0, q, p1, t1, t0, s1, s0;
+  MAYBE_UNUSED mp_limb_t p0;
   mi = -mi;
   affirm ((m1 >> (W_TYPE_SIZE - 1)) == 0);
 
-  /* First compute a0 * <b1, b0> B^{-1}
+  /* First compute a0 * (b1,b0) B^{-1}
         +-----+
         |a0 b0|
      +--+--+--+
@@ -1010,7 +1118,7 @@ mulredc2 (wide_uint *r1p,
   add_ssaaaa (r1, r0, r1, r0, 0, t1);
   add_ssaaaa (r1, r0, r1, r0, s1, s0);
 
-  /* Next, (a1 * <b1, b0> + <r1, r0> B^{-1}
+  /* Next, (a1 * (b1,b0) + (r1,r0) B^{-1}
         +-----+
         |a1 b0|
         +--+--+
@@ -1044,428 +1152,61 @@ mulredc2 (wide_uint *r1p,
   return r0;
 }
 
-ATTRIBUTE_CONST
-static wide_uint
-powm (wide_uint b, wide_uint e, wide_uint n, wide_uint ni, wide_uint one)
-{
-  wide_uint y = one;
-
-  if (e & 1)
-    y = b;
-
-  while (e != 0)
-    {
-      b = mulredc (b, b, n, ni);
-      e >>= 1;
-
-      if (e & 1)
-        y = mulredc (y, b, n, ni);
-    }
-
-  return y;
-}
-
-ATTRIBUTE_PURE static uuint
-powm2 (const wide_uint *bp, const wide_uint *ep, const wide_uint *np,
-       wide_uint ni, const wide_uint *one)
-{
-  wide_uint r1, r0, b1, b0, n1, n0;
-  int i;
-  wide_uint e;
-
-  b0 = bp[0];
-  b1 = bp[1];
-  n0 = np[0];
-  n1 = np[1];
-
-  r0 = one[0];
-  r1 = one[1];
-
-  for (e = ep[0], i = W_TYPE_SIZE; i > 0; i--, e >>= 1)
-    {
-      if (e & 1)
-        {
-          wide_uint r1m1;
-          r0 = mulredc2 (&r1m1, r1, r0, b1, b0, n1, n0, ni);
-          r1 = r1m1;
-        }
-      wide_uint r1m;
-      b0 = mulredc2 (&r1m, b1, b0, b1, b0, n1, n0, ni);
-      b1 = r1m;
-    }
-  for (e = ep[1]; e > 0; e >>= 1)
-    {
-      if (e & 1)
-        {
-          wide_uint r1m1;
-          r0 = mulredc2 (&r1m1, r1, r0, b1, b0, n1, n0, ni);
-          r1 = r1m1;
-        }
-      wide_uint r1m;
-      b0 = mulredc2 (&r1m, b1, b0, b1, b0, n1, n0, ni);
-      b1 = r1m;
-    }
-  return make_uuint (r1, r0);
-}
-
-ATTRIBUTE_CONST
-static bool
-millerrabin (wide_uint n, wide_uint ni, wide_uint b, wide_uint q,
-             int k, wide_uint one)
-{
-  wide_uint y = powm (b, q, n, ni, one);
-
-  wide_uint nm1 = n - one;      /* -1, but in redc representation.  */
-
-  if (y == one || y == nm1)
-    return true;
-
-  for (int i = 1; i < k; i++)
-    {
-      y = mulredc (y, y, n, ni);
-
-      if (y == nm1)
-        return true;
-      if (y == one)
-        return false;
-    }
-  return false;
-}
-
-ATTRIBUTE_PURE static bool
-millerrabin2 (const wide_uint *np, wide_uint ni, const wide_uint *bp,
-              const wide_uint *qp, int k, const wide_uint *one)
-{
-  wide_uint y1, y0, nm1_1, nm1_0, r1m;
-
-  uuset (&y1, &y0, powm2 (bp, qp, np, ni, one));
-
-  if (y0 == one[0] && y1 == one[1])
-    return true;
-
-  sub_ddmmss (nm1_1, nm1_0, np[1], np[0], one[1], one[0]);
-
-  if (y0 == nm1_0 && y1 == nm1_1)
-    return true;
-
-  for (int i = 1; i < k; i++)
-    {
-      y0 = mulredc2 (&r1m, y1, y0, y1, y0, np[1], np[0], ni);
-      y1 = r1m;
-
-      if (y0 == nm1_0 && y1 == nm1_1)
-        return true;
-      if (y0 == one[0] && y1 == one[1])
-        return false;
-    }
-  return false;
-}
+static bool mp_prime_p (mpz_t);
 
+/* Is N prime?  N cannot be even or be a composite number less than
+   SQUARE_OF_FIRST_OMITTED_PRIME.  */
 static bool
-mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
-                mpz_srcptr q, mp_bitcnt_t k)
+prime_p (mp_limb_t n)
 {
-  mpz_powm (y, x, q, n);
-
-  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
-    return true;
-
-  for (mp_bitcnt_t i = 1; i < k; i++)
-    {
-      mpz_powm_ui (y, y, 2, n);
-      if (mpz_cmp (y, nm1) == 0)
-        return true;
-      if (mpz_cmp_ui (y, 1) == 0)
-        return false;
-    }
-  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
-prime_p (wide_uint n)
-{
-  bool is_prime;
-  wide_uint a_prim, one, ni;
-  struct factors factors;
-
   if (n <= 1)
     return false;
 
-  wide_uint cast_out_limit
-    = (wide_uint) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME;
-
-#ifndef EXHIBIT_INT128_BUG
-  /* FIXME: Do the small-prime performance improvement only if
-     wide_uint is exactly 64 bits wide.  We don't know why the code
-     misbehaves when wide_uint is wider; e.g., when compiled with
-     'gcc -DUSE_INT128 -DEXHIBIT_INT128_BUG', 'factor' mishandles
-     340282366920938463463374607431768211355.  */
-  if (W_TYPE_SIZE != 64)
-    cast_out_limit = 2;
-#endif
-
   /* We have already cast out small primes.  */
-  if (n < cast_out_limit)
+  if (n < SQUARE_OF_FIRST_OMITTED_PRIME)
     return true;
 
-  /* Precomputation for Miller-Rabin.  */
-  int k = stdc_trailing_zeros (n - 1);
-  wide_uint q = (n - 1) >> k;
-
-  wide_uint a = 2;
-  binv (ni, n);                 /* ni <- 1/n mod B */
-  redcify (one, 1, n);
-  addmod (a_prim, one, one, n); /* i.e., redcify a = 2 */
-
-  /* Perform a Miller-Rabin test, finds most composites quickly.  */
-  if (!millerrabin (n, ni, a_prim, q, k, one))
-    return false;
-
-  if (flag_prove_primality)
-    {
-      /* Factor n-1 for Lucas.  */
-      factor (0, n - 1, &factors);
-    }
-
-  /* Loop until Lucas proves our number prime, or Miller-Rabin proves our
-     number composite.  */
-  for (idx_t r = 0; r < PRIMES_PTAB_ENTRIES; r++)
-    {
-      if (flag_prove_primality)
-        {
-          is_prime = true;
-          for (int i = 0; i < factors.nfactors && is_prime; i++)
-            {
-              is_prime
-                = powm (a_prim, (n - 1) / factors.p[i], n, ni, one) != one;
-            }
-        }
-      else
-        {
-          /* After enough Miller-Rabin runs, be content.  */
-          is_prime = (r == MR_REPS - 1);
-        }
-
-      if (is_prime)
-        return true;
-
-      a += primes_diff[r];      /* Establish new base.  */
-
-      /* The following is equivalent to redcify (a_prim, a, n).  It runs faster
-         on most processors, since it avoids udiv_qrnnd.  If we go down the
-         udiv_qrnnd_preinv path, this code should be replaced.  */
-      {
-        wide_uint s1, s0;
-        umul_ppmm (s1, s0, one, a);
-        if (LIKELY (s1 == 0))
-          a_prim = s0 % n;
-        else
-          {
-            MAYBE_UNUSED wide_uint dummy;
-            udiv_qrnnd (dummy, a_prim, s1, s0, n);
-          }
-      }
-
-      if (!millerrabin (n, ni, a_prim, q, k, one))
-        return false;
-    }
-
-  affirm (!"Lucas prime test failure.  This should not happen");
+  mpz_t mn = MPZ_ROINIT_N (&n, 1);
+  return mp_prime_p (mn);
 }
 
-static bool ATTRIBUTE_PURE
-prime2_p (wide_uint n1, wide_uint n0)
+/* Is (n1,n0) prime?  (n1,n0) cannot be even or be a composite number
+   less than SQUARE_OF_FIRST_OMITTED_PRIME.  */
+static bool
+prime2_p (mp_limb_t n1, mp_limb_t n0)
 {
-  wide_uint q[2], nm1[2];
-  wide_uint a_prim[2];
-  wide_uint one[2];
-  wide_uint na[2];
-  wide_uint ni;
-  int k;
-  struct factors factors;
-
   if (n1 == 0)
     return prime_p (n0);
 
-  nm1[1] = n1 - (n0 == 0);
-  nm1[0] = n0 - 1;
-  if (nm1[0] == 0)
-    {
-      assume (nm1[1]);
-      k = stdc_trailing_zeros (nm1[1]);
-
-      q[0] = nm1[1] >> k;
-      q[1] = 0;
-      k += W_TYPE_SIZE;
-    }
-  else
-    {
-      k = stdc_trailing_zeros (nm1[0]);
-      rsh2 (q[1], q[0], nm1[1], nm1[0], k);
-    }
-
-  wide_uint a = 2;
-  binv (ni, n0);
-  redcify2 (one[1], one[0], 1, n1, n0);
-  addmod2 (a_prim[1], a_prim[0], one[1], one[0], one[1], one[0], n1, n0);
-
-  /* FIXME: Use scalars or pointers in arguments?  Some consistency needed.  */
-  na[0] = n0;
-  na[1] = n1;
-
-  if (!millerrabin2 (na, ni, a_prim, q, k, one))
-    return false;
-
-  if (flag_prove_primality)
-    {
-      /* Factor n-1 for Lucas.  */
-      factor (nm1[1], nm1[0], &factors);
-    }
-
-  /* Loop until Lucas proves our number prime, or Miller-Rabin proves our
-     number composite.  */
-  for (idx_t r = 0; r < PRIMES_PTAB_ENTRIES; r++)
-    {
-      bool is_prime;
-      wide_uint e[2];
-      uuint y;
-
-      if (flag_prove_primality)
-        {
-          is_prime = true;
-          if (hi (factors.plarge))
-            {
-              wide_uint pi;
-              binv (pi, lo (factors.plarge));
-              e[0] = pi * nm1[0];
-              e[1] = 0;
-              y = powm2 (a_prim, e, na, ni, one);
-              is_prime = (lo (y) != one[0] || hi (y) != one[1]);
-            }
-          for (int i = 0; i < factors.nfactors && is_prime; i++)
-            {
-              /* FIXME: We always have the factor 2.  Do we really need to
-                 handle it here?  We have done the same powering as part
-                 of millerrabin.  */
-              if (factors.p[i] == 2)
-                rsh2 (e[1], e[0], nm1[1], nm1[0], 1);
-              else
-                divexact_21 (e[1], e[0], nm1[1], nm1[0], factors.p[i]);
-              y = powm2 (a_prim, e, na, ni, one);
-              is_prime = (lo (y) != one[0] || hi (y) != one[1]);
-            }
-        }
-      else
-        {
-          /* After enough Miller-Rabin runs, be content.  */
-          is_prime = (r == MR_REPS - 1);
-        }
-
-      if (is_prime)
-        return true;
-
-      a += primes_diff[r];      /* Establish new base.  */
-      redcify2 (a_prim[1], a_prim[0], a, n1, n0);
-
-      if (!millerrabin2 (na, ni, a_prim, q, k, one))
-        return false;
-    }
-
-  affirm (!"Lucas prime test failure.  This should not happen");
+  uuint n = make_uuint (n1, n0);
+  mpz_t mn = MPZ_ROINIT_N (n.uu, 2);
+  return mp_prime_p (mn);
 }
 
+/* Is N prime?  N cannot be even or be a composite number less than
+   SQUARE_OF_FIRST_OMITTED_PRIME.  */
 static bool
 mp_prime_p (mpz_t n)
 {
-  bool is_prime;
-  mpz_t q, a, nm1, tmp;
-  struct mp_factors factors;
-
   if (mpz_cmp_ui (n, 1) <= 0)
     return false;
 
   /* We have already cast out small primes.  */
-  if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0)
+  if (mpz_cmp_ui (n, SQUARE_OF_FIRST_OMITTED_PRIME) < 0)
     return true;
 
-  mpz_inits (q, a, nm1, tmp, nullptr);
-
-  /* Precomputation for Miller-Rabin.  */
-  mpz_sub_ui (nm1, n, 1);
-
-  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
-  mp_bitcnt_t k = mpz_scan1 (nm1, 0);
-  mpz_tdiv_q_2exp (q, nm1, k);
-
-  mpz_set_ui (a, 2);
-
-  /* Perform a Miller-Rabin test, finds most composites quickly.  */
-  if (!mp_millerrabin (n, nm1, a, tmp, q, k))
-    {
-      is_prime = false;
-      goto ret2;
-    }
-
-  if (flag_prove_primality)
-    {
-      /* Factor n-1 for Lucas.  */
-      mpz_set (tmp, nm1);
-      mp_factor (tmp, &factors);
-    }
-
-  /* Loop until Lucas proves our number prime, or Miller-Rabin proves our
-     number composite.  */
-  for (idx_t r = 0; r < PRIMES_PTAB_ENTRIES; r++)
-    {
-      if (flag_prove_primality)
-        {
-          is_prime = true;
-          for (idx_t i = 0; i < factors.nfactors && is_prime; i++)
-            {
-              mpz_divexact (tmp, nm1, factors.p[i]);
-              mpz_powm (tmp, a, tmp, n);
-              is_prime = mpz_cmp_ui (tmp, 1) != 0;
-            }
-        }
-      else
-        {
-          /* After enough Miller-Rabin runs, be content.  */
-          is_prime = (r == MR_REPS - 1);
-        }
-
-      if (is_prime)
-        goto ret1;
-
-      mpz_add_ui (a, a, primes_diff[r]);        /* Establish new base.  */
-
-      if (!mp_millerrabin (n, nm1, a, tmp, q, k))
-        {
-          is_prime = false;
-          goto ret1;
-        }
-    }
-
-  affirm (!"Lucas prime test failure.  This should not happen");
-
- ret1:
-  if (flag_prove_primality)
-    mp_factor_clear (&factors);
- ret2:
-  mpz_clears (q, a, nm1, tmp, nullptr);
-
-  return is_prime;
+  return !!mpz_probab_prime_p (n, MR_REPS);
 }
 
+/* Insert into FACTORS the result of factoring N,
+   using Pollard-rho with starting value A.  N must be odd.  */
 static void
-factor_using_pollard_rho (wide_uint n, unsigned long int a,
-                          struct factors *factors)
+factor_using_pollard_rho (struct factors *factors, mp_limb_t n, mp_limb_t a)
 {
-  wide_uint x, z, y, P, t, ni, g;
+  mp_limb_t x, z, y, P, t, ni, g;
 
-  unsigned long int k = 1;
-  unsigned long int l = 1;
+  int_fast64_t k = 1;
+  int_fast64_t l = 1;
 
   redcify (P, 1, n);
   addmod (x, P, P, n);          /* i.e., redcify(2) */
@@ -1475,7 +1216,7 @@ factor_using_pollard_rho (wide_uint n, unsigned long int a,
     {
       affirm (a < n);
 
-      binv (ni, n);             /* FIXME: when could we use old 'ni' value?  */
+      ni = binv_limb (n);	/* FIXME: when could we use old 'ni' value?  */
 
       for (;;)
         {
@@ -1487,7 +1228,7 @@ factor_using_pollard_rho (wide_uint n, unsigned long int a,
               submod (t, z, x, n);
               P = mulredc (P, t, n, ni);
 
-              if (k % 32 == 1)
+              if ((k & 31) == 1)
                 {
                   if (gcd_odd (P, n) != 1)
                     goto factor_found;
@@ -1499,7 +1240,7 @@ factor_using_pollard_rho (wide_uint n, unsigned long int a,
           z = x;
           k = l;
           l = 2 * l;
-          for (unsigned long int i = 0; i < k; i++)
+          for (int_fast64_t i = 0; i < k; i++)
             {
               x = mulredc (x, x, n, ni);
               addmod (x, x, a, n);
@@ -1521,14 +1262,14 @@ factor_using_pollard_rho (wide_uint n, unsigned long int a,
       if (n == g)
         {
           /* Found n itself as factor.  Restart with different params.  */
-          factor_using_pollard_rho (n, a + 1, factors);
+          factor_using_pollard_rho (factors, n, a + 1);
           return;
         }
 
       n = n / g;
 
       if (!prime_p (g))
-        factor_using_pollard_rho (g, a + 1, factors);
+        factor_using_pollard_rho (factors, g, a + 1);
       else
         factor_insert (factors, g);
 
@@ -1545,13 +1286,13 @@ factor_using_pollard_rho (wide_uint n, unsigned long int a,
 }
 
 static void
-factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
-                           struct factors *factors)
+factor_using_pollard_rho2 (struct factors *factors,
+                           mp_limb_t n1, mp_limb_t n0, mp_limb_t a)
 {
-  wide_uint x1, x0, z1, z0, y1, y0, P1, P0, t1, t0, ni, g1, g0, r1m;
+  mp_limb_t x1, x0, z1, z0, y1, y0, P1, P0, t1, t0, g1, g0, r1m;
 
-  unsigned long int k = 1;
-  unsigned long int l = 1;
+  int_fast64_t k = 1;
+  int_fast64_t l = 1;
 
   redcify2 (P1, P0, 1, n1, n0);
   addmod2 (x1, x0, P1, P0, P1, P0, n1, n0); /* i.e., redcify(2) */
@@ -1560,7 +1301,7 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
 
   while (n1 != 0 || n0 != 1)
     {
-      binv (ni, n0);
+      mp_limb_t ni = binv_limb (n0);
 
       for (;;)
         {
@@ -1568,13 +1309,13 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
             {
               x0 = mulredc2 (&r1m, x1, x0, x1, x0, n1, n0, ni);
               x1 = r1m;
-              addmod2 (x1, x0, x1, x0, 0, (wide_uint) a, n1, n0);
+              addmod2 (x1, x0, x1, x0, 0, a, n1, n0);
 
               submod2 (t1, t0, z1, z0, x1, x0, n1, n0);
               P0 = mulredc2 (&r1m, P1, P0, t1, t0, n1, n0, ni);
               P1 = r1m;
 
-              if (k % 32 == 1)
+              if ((k & 31) == 1)
                 {
                   uuset (&g1, &g0, gcd2_odd (P1, P0, n1, n0));
                   if (g1 != 0 || g0 != 1)
@@ -1587,11 +1328,11 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
           z1 = x1; z0 = x0;
           k = l;
           l = 2 * l;
-          for (unsigned long int i = 0; i < k; i++)
+          for (int_fast64_t i = 0; i < k; i++)
             {
               x0 = mulredc2 (&r1m, x1, x0, x1, x0, n1, n0, ni);
               x1 = r1m;
-              addmod2 (x1, x0, x1, x0, 0, (wide_uint) a, n1, n0);
+              addmod2 (x1, x0, x1, x0, 0, a, n1, n0);
             }
           y1 = x1; y0 = x0;
         }
@@ -1601,7 +1342,7 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
         {
           y0 = mulredc2 (&r1m, y1, y0, y1, y0, n1, n0, ni);
           y1 = r1m;
-          addmod2 (y1, y0, y1, y0, 0, (wide_uint) a, n1, n0);
+          addmod2 (y1, y0, y1, y0, 0, a, n1, n0);
 
           submod2 (t1, t0, z1, z0, y1, y0, n1, n0);
           uuset (&g1, &g0, gcd2_odd (t1, t0, n1, n0));
@@ -1614,7 +1355,7 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
           divexact_21 (n1, n0, n1, n0, g0);     /* n = n / g */
 
           if (!prime_p (g0))
-            factor_using_pollard_rho (g0, a + 1, factors);
+            factor_using_pollard_rho (factors, g0, a + 1);
           else
             factor_insert (factors, g0);
         }
@@ -1622,24 +1363,22 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
         {
           /* The found factor is two words.  This is highly unlikely, thus hard
              to trigger.  Please be careful before you change this code!  */
-          wide_uint ginv;
 
           if (n1 == g1 && n0 == g0)
             {
               /* Found n itself as factor.  Restart with different params.  */
-              factor_using_pollard_rho2 (n1, n0, a + 1, factors);
+              factor_using_pollard_rho2 (factors, n1, n0, a + 1);
               return;
             }
 
           /* Compute n = n / g.  Since the result will fit one word,
              we can compute the quotient modulo B, ignoring the high
              divisor word.  */
-          binv (ginv, g0);
-          n0 = ginv * n0;
+          n0 = binv_limb (g0) * n0;
           n1 = 0;
 
           if (!prime2_p (g1, g0))
-            factor_using_pollard_rho2 (g1, g0, a + 1, factors);
+            factor_using_pollard_rho2 (factors, g1, g0, a + 1);
           else
             factor_insert_large (factors, g1, g0);
         }
@@ -1652,7 +1391,7 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
               break;
             }
 
-          factor_using_pollard_rho (n0, a, factors);
+          factor_using_pollard_rho (factors, n0, a);
           return;
         }
 
@@ -1668,551 +1407,169 @@ factor_using_pollard_rho2 (wide_uint n1, wide_uint n0, unsigned long int a,
     }
 }
 
+/* Set RP = (AP + BP) mod MP.  All values are nonnegative and take up
+   N>0 words, and AP and BP are both less than MP.  */
 static void
-mp_factor_using_pollard_rho (mpz_t n, unsigned long int a,
-                             struct mp_factors *factors)
+mp_modadd (mp_limb_t *rp, mp_limb_t const *ap, mp_limb_t const *bp,
+           mp_limb_t const *mp, mp_size_t n)
 {
-  mpz_t x, z, y, P;
-  mpz_t t, t2;
-
-  devmsg ("[pollard-rho (%lu)] ", a);
-
-  mpz_inits (t, t2, nullptr);
-  mpz_init_set_si (y, 2);
-  mpz_init_set_si (x, 2);
-  mpz_init_set_si (z, 2);
-  mpz_init_set_ui (P, 1);
-
-  unsigned long long int k = 1;
-  unsigned long long int l = 1;
-
-  while (mpz_cmp_ui (n, 1) != 0)
-    {
-      for (;;)
-        {
-          do
-            {
-              mpz_mul (t, x, x);
-              mpz_mod (x, t, n);
-              mpz_add_ui (x, x, a);
-
-              mpz_sub (t, z, x);
-              mpz_mul (t2, P, t);
-              mpz_mod (P, t2, n);
-
-              if (k % 32 == 1)
-                {
-                  mpz_gcd (t, P, n);
-                  if (mpz_cmp_ui (t, 1) != 0)
-                    goto factor_found;
-                  mpz_set (y, x);
-                }
-            }
-          while (--k != 0);
-
-          mpz_set (z, x);
-          k = l;
-          l = 2 * l;
-          for (unsigned long long int i = 0; i < k; i++)
-            {
-              mpz_mul (t, x, x);
-              mpz_mod (x, t, n);
-              mpz_add_ui (x, x, a);
-            }
-          mpz_set (y, x);
-        }
-
-    factor_found:
-      do
-        {
-          mpz_mul (t, y, y);
-          mpz_mod (y, t, n);
-          mpz_add_ui (y, y, a);
-
-          mpz_sub (t, z, y);
-          mpz_gcd (t, t, n);
-        }
-      while (mpz_cmp_ui (t, 1) == 0);
-
-      mpz_divexact (n, n, t);   /* divide by t, before t is overwritten */
-
-      if (!mp_prime_p (t))
-        {
-          devmsg ("[composite factor--restarting pollard-rho] ");
-          mp_factor_using_pollard_rho (t, a + 1, factors);
-        }
-      else
-        {
-          mp_factor_insert (factors, t);
-        }
-
-      if (mp_prime_p (n))
-        {
-          mp_factor_insert (factors, n);
-          break;
-        }
-
-      mpz_mod (x, x, n);
-      mpz_mod (z, z, n);
-      mpz_mod (y, y, n);
-    }
-
-  mpz_clears (P, t2, t, z, x, y, nullptr);
+  mp_limb_t cy = mpn_add_n (rp, ap, bp, n);
+  if (cy || mpn_cmp (rp, mp, n) >= 0)
+    mpn_sub_n (rp, rp, mp, n);
 }
 
-#if USE_SQUFOF
-/* FIXME: Maybe better to use an iteration converging to 1/sqrt(n)?  If
-   algorithm is replaced, consider also returning the remainder.  */
-ATTRIBUTE_CONST
-static wide_uint
-isqrt (wide_uint n)
+/* Set RP = (AP - BP) mod MP.  All values are nonnegative and take up
+   N>0 words, and AP and BP are both less than MP.  */
+static void
+mp_modsub (mp_limb_t *rp, mp_limb_t const *ap, mp_limb_t const *bp,
+           mp_limb_t const *mp, mp_size_t n)
 {
-  if (n == 0)
-    return 0;
-
-  int c = stdc_leading_zeros (n);
-
-  /* Make x > sqrt(n).  This will be invariant through the loop.  */
-  wide_uint x = (wide_uint) 1 << ((W_TYPE_SIZE + 1 - c) >> 1);
-
-  for (;;)
-    {
-      wide_uint y = (x + n / x) / 2;
-      if (y >= x)
-        return x;
-
-      x = y;
-    }
+  mp_limb_t cy = mpn_sub_n (rp, ap, bp, n);
+  if (cy)
+    mpn_add_n (rp, rp, mp, n);
 }
 
-ATTRIBUTE_CONST
-static wide_uint
-isqrt2 (wide_uint nh, wide_uint nl)
+/* Set RP = (AP - B0) mod MP.  All values are nonnegative, AP and MP
+   both take up N>0 words, and AP < MP.  */
+static void
+mp_modadd_1 (mp_limb_t *rp, mp_limb_t const *ap, mp_limb_t b0,
+             mp_limb_t const *mp, mp_size_t n)
 {
-  /* Ensures the remainder fits in a wide_uint.  */
-  affirm (nh < ((wide_uint) 1 << (W_TYPE_SIZE - 2)));
-
-  if (nh == 0)
-    return isqrt (nl);
-
-  int shift = stdc_leading_zeros (nh) & ~1;
-
-  /* Make x > sqrt (n).  */
-  wide_uint x = isqrt ((nh << shift) + (nl >> (W_TYPE_SIZE - shift))) + 1;
-  x <<= (W_TYPE_SIZE - shift) >> 1;
-
-  /* Do we need more than one iteration?  */
-  for (;;)
-    {
-      MAYBE_UNUSED wide_uint r;
-      wide_uint q, y;
-      udiv_qrnnd (q, r, nh, nl, x);
-      y = (x + q) / 2;
-
-      if (y >= x)
-        {
-          wide_uint hi, lo;
-          umul_ppmm (hi, lo, x + 1, x + 1);
-          affirm (gt2 (hi, lo, nh, nl));
-
-          umul_ppmm (hi, lo, x, x);
-          affirm (ge2 (nh, nl, hi, lo));
-          sub_ddmmss (hi, lo, nh, nl, hi, lo);
-          affirm (hi == 0);
-
-          return x;
-        }
-
-      x = y;
-    }
+  mp_limb_t cy = mpn_add_1 (rp, ap, n, b0);
+  if (cy || mpn_cmp (rp, mp, n) >= 0)
+    mpn_sub_n (rp, rp, mp, n);
 }
 
-/* MAGIC[N] has a bit i set iff i is a quadratic residue mod N.  */
-# define MAGIC64 0x0202021202030213ULL
-# define MAGIC63 0x0402483012450293ULL
-# define MAGIC65 0x218a019866014613ULL
-# define MAGIC11 0x23b
-
-/* Return the square root if the input is a square, otherwise 0.  */
-ATTRIBUTE_CONST
-static wide_uint
-is_square (wide_uint x)
+static void
+mp_mulredc (mp_limb_t *rp, mp_limb_t const *ap, mp_limb_t const *bp,
+            mp_limb_t const *mp, mp_size_t n, mp_limb_t m0inv, mp_limb_t *tp)
 {
-  /* Uses the tests suggested by Cohen.  Excludes 99% of the non-squares before
-     computing the square root.  */
-  if (((MAGIC64 >> (x & 63)) & 1)
-      && ((MAGIC63 >> (x % 63)) & 1)
-      /* Both 0 and 64 are squares mod (65).  */
-      && ((MAGIC65 >> ((x % 65) & 63)) & 1)
-      && ((MAGIC11 >> (x % 11) & 1)))
+  tp[n] = mpn_mul_1 (tp, ap, n, bp[0]);
+  tp[0] = mpn_addmul_1 (tp, mp, n, tp[0] * m0inv);
+
+  for (mp_size_t i = 1; i < n; i++)
     {
-      wide_uint r = isqrt (x);
-      if (r * r == x)
-        return r;
+      tp[n + i] = mpn_addmul_1 (tp + i, ap, n, bp[i]);
+      tp[i] = mpn_addmul_1 (tp + i, mp, n, tp[i] * m0inv);
     }
-  return 0;
+  mp_modadd (rp, tp, tp + n, mp, n);
 }
 
-/* invtab[i] = floor (0x10000 / (0x100 + i) */
-static short const invtab[0x81] =
-  {
-    0x200,
-    0x1fc, 0x1f8, 0x1f4, 0x1f0, 0x1ec, 0x1e9, 0x1e5, 0x1e1,
-    0x1de, 0x1da, 0x1d7, 0x1d4, 0x1d0, 0x1cd, 0x1ca, 0x1c7,
-    0x1c3, 0x1c0, 0x1bd, 0x1ba, 0x1b7, 0x1b4, 0x1b2, 0x1af,
-    0x1ac, 0x1a9, 0x1a6, 0x1a4, 0x1a1, 0x19e, 0x19c, 0x199,
-    0x197, 0x194, 0x192, 0x18f, 0x18d, 0x18a, 0x188, 0x186,
-    0x183, 0x181, 0x17f, 0x17d, 0x17a, 0x178, 0x176, 0x174,
-    0x172, 0x170, 0x16e, 0x16c, 0x16a, 0x168, 0x166, 0x164,
-    0x162, 0x160, 0x15e, 0x15c, 0x15a, 0x158, 0x157, 0x155,
-    0x153, 0x151, 0x150, 0x14e, 0x14c, 0x14a, 0x149, 0x147,
-    0x146, 0x144, 0x142, 0x141, 0x13f, 0x13e, 0x13c, 0x13b,
-    0x139, 0x138, 0x136, 0x135, 0x133, 0x132, 0x130, 0x12f,
-    0x12e, 0x12c, 0x12b, 0x129, 0x128, 0x127, 0x125, 0x124,
-    0x123, 0x121, 0x120, 0x11f, 0x11e, 0x11c, 0x11b, 0x11a,
-    0x119, 0x118, 0x116, 0x115, 0x114, 0x113, 0x112, 0x111,
-    0x10f, 0x10e, 0x10d, 0x10c, 0x10b, 0x10a, 0x109, 0x108,
-    0x107, 0x106, 0x105, 0x104, 0x103, 0x102, 0x101, 0x100,
-  };
-
-/* Compute q = [u/d], r = u mod d.  Avoids slow hardware division for the case
-   that q < 0x40; here it instead uses a table of (Euclidean) inverses.  */
-# define div_smallq(q, r, u, d)                                          \
-  do {                                                                  \
-    if ((u) / 0x40 < (d))                                               \
-      {                                                                 \
-        wide_uint _dinv, _mask, _q, _r;                                 \
-        int _cnt = stdc_leading_zeros (d);				\
-        _r = (u);                                                       \
-        if (UNLIKELY (_cnt > (W_TYPE_SIZE - 8)))                        \
-          {                                                             \
-            _dinv = invtab[((d) << (_cnt + 8 - W_TYPE_SIZE)) - 0x80];   \
-            _q = _dinv * _r >> (8 + W_TYPE_SIZE - _cnt);                \
-          }                                                             \
-        else                                                            \
-          {                                                             \
-            _dinv = invtab[((d) >> (W_TYPE_SIZE - 8 - _cnt)) - 0x7f];   \
-            _q = _dinv * (_r >> (W_TYPE_SIZE - 3 - _cnt)) >> 11;        \
-          }                                                             \
-        _r -= _q * (d);                                                 \
-                                                                        \
-        _mask = -(wide_uint) (_r >= (d));                               \
-        (r) = _r - (_mask & (d));                                       \
-        (q) = _q - _mask;                                               \
-        affirm ((q) * (d) + (r) == u);					\
-      }                                                                 \
-    else                                                                \
-      {                                                                 \
-        wide_uint _q = (u) / (d);                                       \
-        (r) = (u) - _q * (d);                                           \
-        (q) = _q;                                                       \
-      }                                                                 \
-  } while (0)
-
-/* Notes: Example N = 22117019.  After first phase we find Q1 = 6314, Q
-   = 3025, P = 1737, representing F_{18} = (-6314, 2 * 1737, 3025),
-   with 3025 = 55^2.
-
-   Constructing the square root, we get Q1 = 55, Q = 8653, P = 4652,
-   representing G_0 = (-55, 2 * 4652, 8653).
-
-   In the notation of the paper:
-
-   S_{-1} = 55, S_0 = 8653, R_0 = 4652
+/* Maximum value for N in mp_factor_using_pollard_rho,
+   to avoid integer overflow when it calls xinmalloc.  */
+#define MP_FACTOR_USING_POLLARD_RHO_N_MAX \
+  ((MIN (IDX_MAX, TYPE_MAXIMUM (mp_size_t)) - 3) / 10)
 
-   Put
-
-     t_0 = floor([q_0 + R_0] / S0) = 1
-     R_1 = t_0 * S_0 - R_0 = 4001
-     S_1 = S_{-1} +t_0 (R_0 - R_1) = 706
-*/
-
-/* Multipliers, in order of efficiency:
-   0.7268  3*5*7*11 = 1155 = 3 (mod 4)
-   0.7317  3*5*7    =  105 = 1
-   0.7820  3*5*11   =  165 = 1
-   0.7872  3*5      =   15 = 3
-   0.8101  3*7*11   =  231 = 3
-   0.8155  3*7      =   21 = 1
-   0.8284  5*7*11   =  385 = 1
-   0.8339  5*7      =   35 = 3
-   0.8716  3*11     =   33 = 1
-   0.8774  3        =    3 = 3
-   0.8913  5*11     =   55 = 3
-   0.8972  5        =    5 = 1
-   0.9233  7*11     =   77 = 1
-   0.9295  7        =    7 = 3
-   0.9934  11       =   11 = 3
-*/
-# define QUEUE_SIZE 50
-#endif
-
-#if STAT_SQUFOF
-# define Q_FREQ_SIZE 50
-/* Element 0 keeps the total */
-static int q_freq[Q_FREQ_SIZE + 1];
-#endif
-
-#if USE_SQUFOF
-/* Return true on success.  Expected to fail only for numbers
-   >= 2^{2*W_TYPE_SIZE - 2}, or close to that limit.  */
-static bool
-factor_using_squfof (wide_uint n1, wide_uint n0, struct factors *factors)
+/* Insert into FACTORS the result of factoring MP, of size N,
+   using Pollard-rho with starting value A.  MP must be odd. */
+static void
+mp_factor_using_pollard_rho (struct mp_factors *factors,
+                             mp_limb_t const *mp, mp_size_t n,
+                             mp_limb_t a)
 {
-  /* Uses algorithm and notation from
+  devmsg ("[pollard-rho (%lu)] ", (unsigned long int) a);
 
-     SQUARE FORM FACTORIZATION
-     JASON E. GOWER AND SAMUEL S. WAGSTAFF, JR.
+  static_assert (10 * MP_FACTOR_USING_POLLARD_RHO_N_MAX + 3
+                 <= MIN (IDX_MAX, TYPE_MAXIMUM (mp_size_t)));
+  mp_limb_t *scratch = xinmalloc (10 * n + 3, sizeof *scratch);
+  mp_limb_t *qp = scratch + 2 * n + 1, *pp = qp + n + 2,
+    *xp = pp + n, *yp = xp + n, *zp = yp + n,
+    *tp = zp + n, *sp = tp + n, *gp = sp + n;
+  mp_size_t gn;
 
-     https://homes.cerias.purdue.edu/~ssw/squfof.pdf
-   */
+  mpn_zero (scratch, n);
+  scratch[n] = 1;
+  mpn_tdiv_qr (qp, pp, 0, scratch, n + 1, mp, n);
 
-  static short const multipliers_1[] =
-    { /* = 1 (mod 4) */
-      105, 165, 21, 385, 33, 5, 77, 1, 0
-    };
-  static short const multipliers_3[] =
-    { /* = 3 (mod 4) */
-      1155, 15, 231, 35, 3, 55, 7, 11, 0
-    };
+  mp_modadd (xp, pp, pp, mp, n);
+  mpn_copyi (yp, xp, n);
+  mpn_copyi (zp, xp, n);
 
-  struct { wide_uint Q; wide_uint P; } queue[QUEUE_SIZE];
-
-  if (n1 >= ((wide_uint) 1 << (W_TYPE_SIZE - 2)))
-    return false;
+  mp_limb_t m0inv = binv_limb (-mp[0]);
 
-  wide_uint sqrt_n = isqrt2 (n1, n0);
-
-  if (n0 == sqrt_n * sqrt_n)
+  for (int_fast64_t k = 1; ; k *= 2)
     {
-      wide_uint p1, p0;
+      for (int_fast64_t i = k; 0 < i; i--)
+        {
+          mp_mulredc (tp, xp, xp, mp, n, m0inv, scratch);
+          mp_modadd_1 (xp, tp, a, mp, n);
 
-      umul_ppmm (p1, p0, sqrt_n, sqrt_n);
-      affirm (p0 == n0);
+          mp_modsub (tp, zp, xp, mp, n);
+          mp_mulredc (pp, pp, tp, mp, n, m0inv, scratch);
 
-      if (n1 == p1)
-        {
-          if (prime_p (sqrt_n))
-            factor_insert_multiplicity (factors, sqrt_n, 2);
-          else
+          if (i % 128 == 1)
             {
-              struct factors f;
-
-              f.nfactors = 0;
-              if (!factor_using_squfof (0, sqrt_n, &f))
+              if (mpn_zero_p (pp, n))
                 {
-                  /* Try pollard rho instead */
-                  factor_using_pollard_rho (sqrt_n, 1, &f);
+                  mp_factor_using_pollard_rho (factors, mp, n, a + 1);
+                  goto finish;
                 }
-              /* Duplicate the new factors */
-              for (unsigned int i = 0; i < f.nfactors; i++)
-                factor_insert_multiplicity (factors, f.p[i], 2 * f.e[i]);
+              mpn_copyi (tp, pp, n);
+              mpn_copyi (sp, mp, n);
+              gn = mpn_gcd (gp, tp, n, sp, n);
+              if (gn != 1 || gp[0] != 1)
+                goto factor_found;
+              mpn_copyi (yp, xp, n);
             }
-          return true;
         }
+
+      mpn_copyi (zp, xp, n);
+      for (int_fast64_t i = 2 * k; 0 < i; i--)
+        {
+          mp_mulredc (tp, xp, xp, mp, n, m0inv, scratch);
+          mp_modadd_1 (xp, tp, a, mp, n);
+        }
+      mpn_copyi (yp, xp, n);
     }
 
-  /* Select multipliers so we always get n * mu = 3 (mod 4) */
-  for (short const *m = (n0 % 4 == 1) ? multipliers_3 : multipliers_1;
-       *m; m++)
+ factor_found:
+  do
     {
-      wide_uint S, Dh, Dl, Q1, Q, P, L, L1, B;
-      unsigned int i;
-      unsigned int mu = *m;
-      int qpos = 0;
-
-      affirm (mu * n0 % 4 == 3);
+      mp_mulredc (tp, yp, yp, mp, n, m0inv, scratch);
+      mp_modadd_1 (yp, tp, a, mp, n);
+      mp_modsub (tp, zp, yp, mp, n);
+      mpn_copyi (sp, mp, n);
+      gn = mpn_gcd (gp, tp, n, sp, n);
+    }
+  while (gn == 1 && gp[0] == 1);
 
-      /* In the notation of the paper, with mu * n == 3 (mod 4), we
-         get \Delta = 4 mu * n, and the paper's \mu is 2 mu.  As far as
-         I understand it, the necessary bound is 4 \mu^3 < n, or 32
-         mu^3 < n.
+  mpz_t g = MPZ_ROINIT_N (gp, gn);
+  mpz_t m = MPZ_ROINIT_N ((mp_limb_t *) mp, n), q;
+  mpz_init (q);
+  mpz_divexact (q, m, g);
 
-         However, this seems insufficient: With n = 37243139 and mu =
-         105, we get a trivial factor, from the square 38809 = 197^2,
-         without any corresponding Q earlier in the iteration.
+  if (!mp_finish_in_single (factors, g))
+    {
+      if (mp_prime_p (g))
+        mp_factor_insert (factors, g, 1);
+      else
+        mp_factor_using_pollard_rho (factors, gp, gn, a + 1);
+    }
 
-         Requiring 64 mu^3 < n seems sufficient.  */
-      if (n1 == 0)
-        {
-          if ((wide_uint) mu * mu * mu >= n0 / 64)
-            continue;
-        }
+  if (!mp_finish_in_single (factors, q))
+    {
+      if (mp_prime_p (q))
+        mp_factor_insert (factors, q, 1);
       else
         {
-          if (n1 > ((wide_uint) 1 << (W_TYPE_SIZE - 2)) / mu)
-            continue;
+          devmsg ("[composite factor--restarting pollard-rho] ");
+          mp_factor_using_pollard_rho (factors, mpz_limbs_read (q),
+                                       mp_size (q), a + 1);
         }
-      umul_ppmm (Dh, Dl, n0, mu);
-      Dh += n1 * mu;
-
-      affirm (Dl % 4 != 1);
-      affirm (Dh < (wide_uint) 1 << (W_TYPE_SIZE - 2));
-
-      S = isqrt2 (Dh, Dl);
-
-      Q1 = 1;
-      P = S;
-
-      /* Square root remainder fits in one word, so ignore high part.  */
-      Q = Dl - P * P;
-      /* FIXME: When can this differ from floor (sqrt (2 * sqrt (D)))?  */
-      L = isqrt (2 * S);
-      B = 2 * L;
-      L1 = mu * 2 * L;
-
-      /* The form is (+/- Q1, 2P, -/+ Q), of discriminant 4 (P^2 + Q Q1) =
-         4 D.  */
-
-      for (i = 0; i <= B; i++)
-        {
-          wide_uint q, P1, t, rem;
-
-          div_smallq (q, rem, S + P, Q);
-          P1 = S - rem; /* P1 = q*Q - P */
-
-          affirm (q > 0 && Q > 0);
-
-# if STAT_SQUFOF
-          q_freq[0]++;
-          q_freq[MIN (q, Q_FREQ_SIZE)]++;
-# endif
-
-          if (Q <= L1)
-            {
-              wide_uint g = Q;
-
-              if ((Q & 1) == 0)
-                g /= 2;
-
-              g /= gcd_odd (g, mu);
-
-              if (g <= L)
-                {
-                  if (qpos >= QUEUE_SIZE)
-                    error (EXIT_FAILURE, 0, _("squfof queue overflow"));
-                  queue[qpos].Q = g;
-                  queue[qpos].P = P % g;
-                  qpos++;
-                }
-            }
+    }
 
-          /* I think the difference can be either sign, but mod
-             2^W_TYPE_SIZE arithmetic should be fine.  */
-          t = Q1 + q * (P - P1);
-          Q1 = Q;
-          Q = t;
-          P = P1;
+  mpz_clear (q);
 
-          if ((i & 1) == 0)
-            {
-              wide_uint r = is_square (Q);
-              if (r)
-                {
-                  for (int j = 0; j < qpos; j++)
-                    {
-                      if (queue[j].Q == r)
-                        {
-                          if (r == 1)
-                            /* Traversed entire cycle.  */
-                            goto next_multiplier;
-
-                          /* Need the absolute value for divisibility test.  */
-                          if (P >= queue[j].P)
-                            t = P - queue[j].P;
-                          else
-                            t = queue[j].P - P;
-                          if (t % r == 0)
-                            {
-                              /* Delete entries up to and including entry
-                                 j, which matched.  */
-                              memmove (queue, queue + j + 1,
-                                       (qpos - j - 1) * sizeof (queue[0]));
-                              qpos -= (j + 1);
-                            }
-                          goto next_i;
-                        }
-                    }
-
-                  /* We have found a square form, which should give a
-                     factor.  */
-                  Q1 = r;
-                  affirm (S >= P); /* What signs are possible?  */
-                  P += r * ((S - P) / r);
-
-                  /* Note: Paper says (N - P*P) / Q1, that seems incorrect
-                     for the case D = 2N.  */
-                  /* Compute Q = (D - P*P) / Q1, but we need double
-                     precision.  */
-                  wide_uint hi, lo;
-                  umul_ppmm (hi, lo, P, P);
-                  sub_ddmmss (hi, lo, Dh, Dl, hi, lo);
-                  udiv_qrnnd (Q, rem, hi, lo, Q1);
-                  affirm (rem == 0);
-
-                  for (;;)
-                    {
-                      /* Note: There appears to by a typo in the paper,
-                         Step 4a in the algorithm description says q <--
-                         floor([S+P]/\hat Q), but looking at the equations
-                         in Sec. 3.1, it should be q <-- floor([S+P] / Q).
-                         (In this code, \hat Q is Q1).  */
-                      div_smallq (q, rem, S + P, Q);
-                      P1 = S - rem;     /* P1 = q*Q - P */
-
-# if STAT_SQUFOF
-                      q_freq[0]++;
-                      q_freq[MIN (q, Q_FREQ_SIZE)]++;
-# endif
-                      if (P == P1)
-                        break;
-                      t = Q1 + q * (P - P1);
-                      Q1 = Q;
-                      Q = t;
-                      P = P1;
-                    }
-
-                  if ((Q & 1) == 0)
-                    Q /= 2;
-                  Q /= gcd_odd (Q, mu);
-
-                  affirm (Q > 1 && (n1 || Q < n0));
-
-                  if (prime_p (Q))
-                    factor_insert (factors, Q);
-                  else if (!factor_using_squfof (0, Q, factors))
-                    factor_using_pollard_rho (Q, 2, factors);
-
-                  divexact_21 (n1, n0, n1, n0, Q);
-
-                  if (prime2_p (n1, n0))
-                    factor_insert_large (factors, n1, n0);
-                  else
-                    {
-                      if (!factor_using_squfof (n1, n0, factors))
-                        {
-                          if (n1 == 0)
-                            factor_using_pollard_rho (n0, 1, factors);
-                          else
-                            factor_using_pollard_rho2 (n1, n0, 1, factors);
-                        }
-                    }
-
-                  return true;
-                }
-            }
-        next_i:;
-        }
-    next_multiplier:;
-    }
-  return false;
+ finish:
+  free (scratch);
 }
-#endif
 
-/* Compute the prime factors of the 128-bit number (T1,T0), and put the
-   results in FACTORS.  */
+/* Insert into FACTORS the prime factors of the two-word number (T1,T0).
+   Primes less than the prime with index PRIME_IDX
+   have already been considered, and need not be looked at again.  */
 static void
-factor (wide_uint t1, wide_uint t0, struct factors *factors)
+factor_up (struct factors *factors, mp_limb_t t1, mp_limb_t t0,
+           idx_t prime_idx)
 {
   factors->nfactors = 0;
   hiset (&factors->plarge, 0);
@@ -2220,7 +1577,7 @@ factor (wide_uint t1, wide_uint t0, struct factors *factors)
   if (t1 == 0 && t0 < 2)
     return;
 
-  uuset (&t1, &t0, factor_using_division (t1, t0, factors));
+  uuset (&t1, &t0, factor_using_division (factors, t1, t0, prime_idx));
 
   if (t1 == 0 && t0 < 2)
     return;
@@ -2229,102 +1586,81 @@ factor (wide_uint t1, wide_uint t0, struct factors *factors)
     factor_insert_large (factors, t1, t0);
   else
     {
-#if USE_SQUFOF
-      if (factor_using_squfof (t1, t0, factors))
-        return;
-#endif
-
       if (t1 == 0)
-        factor_using_pollard_rho (t0, 1, factors);
+        factor_using_pollard_rho (factors, t0, 1);
       else
-        factor_using_pollard_rho2 (t1, t0, 1, factors);
+        factor_using_pollard_rho2 (factors, t1, t0, 1);
     }
 }
 
-/* Use Pollard-rho to compute the prime factors of
-   arbitrary-precision T, and put the results in FACTORS.  */
+/* Compute the prime factors of the two-word number (T1,T0),
+   and put the results in FACTORS.  */
 static void
-mp_factor (mpz_t t, struct mp_factors *factors)
+factor (struct factors *factors, mp_limb_t t1, mp_limb_t t0)
 {
-  mp_factor_init (factors);
+  return factor_up (factors, t1, t0, 0);
+}
+
+/* Return the prime factors of T.  */
+static struct mp_factors
+mp_factor (mpz_t t)
+{
+  struct mp_factors factors = mp_no_factors ();
 
   if (mpz_sgn (t) != 0)
     {
-      mp_factor_using_division (t, factors);
+      factors = mp_factor_using_division (t);
 
       if (mpz_cmp_ui (t, 1) != 0)
         {
           devmsg ("[is number prime?] ");
           if (mp_prime_p (t))
-            mp_factor_insert (factors, t);
+            mp_factor_insert (&factors, t, 1);
           else
-            mp_factor_using_pollard_rho (t, 1, factors);
+            mp_factor_using_pollard_rho (&factors, mpz_limbs_read (t),
+                                         mp_size (t), 1);
         }
     }
+
+  return factors;
 }
 
+/* Convert to *U the value represented by S, and return LONGINT_OK.
+   However, on error simply return a value other than LONGINT_OK.  */
 static strtol_error
-strto2wide_uint (wide_uint *hip, wide_uint *lop, char const *s)
+strtouuint (uuint *u, char const *s)
 {
-  int lo_carry;
-  wide_uint hi = 0, lo = 0;
-
-  strtol_error err = LONGINT_INVALID;
-
-  /* Initial scan for invalid digits.  */
-  char const *p = s;
-  for (;;)
-    {
-      unsigned char c = *p++;
-      if (c == 0)
-        break;
-
-      if (UNLIKELY (!c_isdigit (c)))
-        {
-          err = LONGINT_INVALID;
-          break;
-        }
+  mp_limb_t hi = 0, lo = *s++ - '0';
 
-      err = LONGINT_OK;           /* we've seen at least one valid digit */
-    }
+  if (UNLIKELY (9 < lo))
+    return LONGINT_INVALID;
 
-  while (err == LONGINT_OK)
+  for (; LIKELY (0 <= *s - '0' && *s - '0' <= 9); s++)
     {
-      unsigned char c = *s++;
-      if (c == 0)
-        break;
-
-      c -= '0';
-
       if (UNLIKELY (ckd_mul (&hi, hi, 10)))
-        {
-          err = LONGINT_OVERFLOW;
-          break;
-        }
+        return LONGINT_OVERFLOW;
 
-      lo_carry = (lo >> (W_TYPE_SIZE - 3)) + (lo >> (W_TYPE_SIZE - 1));
+      int lo_carry = (lo >> (W_TYPE_SIZE - 3)) + (lo >> (W_TYPE_SIZE - 1));
       lo_carry += 10 * lo < 2 * lo;
 
       lo = 10 * lo;
-      lo_carry += ckd_add (&lo, lo, c);
+      lo_carry += ckd_add (&lo, lo, *s - '0');
       if (UNLIKELY (ckd_add (&hi, hi, lo_carry)))
-        {
-          err = LONGINT_OVERFLOW;
-          break;
-        }
+        return LONGINT_OVERFLOW;
     }
 
-  *hip = hi;
-  *lop = lo;
+  if (UNLIKELY (*s))
+    return LONGINT_INVALID;
 
-  return err;
+  *u = make_uuint (hi, lo);
+  return LONGINT_OK;
 }
 
 /* FACTOR_PIPE_BUF is chosen to give good performance,
    and also is the max guaranteed size that
    consumers can read atomically through pipes.
    Also it's big enough to cater for max line length
-   even with 128 bit wide_uint.  */
+   even with 128 bit word.  */
 #ifndef _POSIX_PIPE_BUF
 # define _POSIX_PIPE_BUF 512
 #endif
@@ -2412,7 +1748,7 @@ lbuf_putnl (void)
    everything from BUFEND to lbuf_buf's end.  Use the area just before
    BUFEND temporarily.  */
 static void
-lbuf_putint_append (wide_uint i, char *bufend)
+lbuf_putint_append (mp_limb_t i, char *bufend)
 {
   char *istr = bufend;
   do
@@ -2432,21 +1768,29 @@ lbuf_putint_append (wide_uint i, char *bufend)
 
 /* Append the string representation of I to lbuf_buf.  */
 static void
-lbuf_putint (wide_uint i)
+lbuf_putint (mp_limb_t i)
 {
   lbuf_putint_append (i, lbuf_buf + sizeof lbuf_buf);
 }
+static void
+lbuf_putbitcnt (mp_bitcnt_t i)
+{
+  char *bufend = lbuf_buf + sizeof lbuf_buf;
+  for (; MP_LIMB_MAX < i; i /= 10)
+    *--bufend = '0' + i % 10;
+  lbuf_putint_append (i, bufend);
+}
 
 /* Append the string representation of T to lbuf_buf.  */
 static void
 print_uuint (uuint t)
 {
-  wide_uint t1 = hi (t), t0 = lo (t);
+  mp_limb_t t1 = hi (t), t0 = lo (t);
   char *bufend = lbuf_buf + sizeof lbuf_buf;
 
   while (t1)
     {
-      wide_uint r = t1 % BIG_POWER_OF_10;
+      mp_limb_t r = t1 % BIG_POWER_OF_10;
       t1 /= BIG_POWER_OF_10;
       udiv_qrnnd (t0, r, r, t0, BIG_POWER_OF_10);
       for (int i = 0; i < LOG_BIG_POWER_OF_10; i++)
@@ -2459,7 +1803,7 @@ print_uuint (uuint t)
   lbuf_putint_append (t0, bufend);
 }
 
-/* Buffer an mpz to the internal LBUF, possibly writing if it is long.  */
+/* Buffer the mpz I to lbuf_buf, possibly writing if it is long.  */
 static void
 lbuf_putmpz (mpz_t const i)
 {
@@ -2486,16 +1830,16 @@ lbuf_putmpz (mpz_t const i)
     }
 }
 
-/* Single-precision factoring */
+/* Emit the factors of T, which is less than B^2 / 2.  */
 static void
-print_factors_single (wide_uint t1, wide_uint t0)
+print_factors_single (uuint t)
 {
   struct factors factors;
 
-  print_uuint (make_uuint (t1, t0));
+  print_uuint (t);
   lbuf_putc (':');
 
-  factor (t1, t0, &factors);
+  factor (&factors, hi (t), lo (t));
 
   for (int j = 0; j < factors.nfactors; j++)
     for (int k = 0; k < factors.e[j]; k++)
@@ -2510,7 +1854,7 @@ print_factors_single (wide_uint t1, wide_uint t0)
           }
       }
 
-  if (hi (factors.plarge))
+  if (hi_is_set (&factors.plarge))
     {
       lbuf_putc (' ');
       print_uuint (factors.plarge);
@@ -2533,25 +1877,23 @@ print_factors (char const *input)
     str++;
   str += *str == '+';
 
-  wide_uint t1, t0;
-
   /* Try converting the number to one or two words.  If it fails, use GMP or
      print an error message.  The 2nd condition checks that the most
      significant bit of the two-word number is clear, in a typesize neutral
      way.  */
-  strtol_error err = strto2wide_uint (&t1, &t0, str);
+  uuint u;
+  strtol_error err = strtouuint (&u, str);
 
   switch (err)
     {
     case LONGINT_OK:
-      if (((t1 << 1) >> 1) == t1)
+      if (hi (u) >> (W_TYPE_SIZE - 1) == 0)
         {
           devmsg ("[using single-precision arithmetic] ");
-          print_factors_single (t1, t0);
+          print_factors_single (u);
           return true;
         }
-      break;
-
+      FALLTHROUGH;
     case LONGINT_OVERFLOW:
       /* Try GMP.  */
       break;
@@ -2566,23 +1908,23 @@ print_factors (char const *input)
 
   devmsg ("[using arbitrary-precision arithmetic] ");
   mpz_t t;
-  struct mp_factors factors;
-
   mpz_init_set_str (t, str, 10);
+  if (MP_FACTOR_USING_POLLARD_RHO_N_MAX < mp_size (t))
+    xalloc_die ();
 
   lbuf_putmpz (t);
   lbuf_putc (':');
-  mp_factor (t, &factors);
+  struct mp_factors factors = mp_factor (t);
 
   for (idx_t j = 0; j < factors.nfactors; j++)
-    for (unsigned long int k = 0; k < factors.e[j]; k++)
+    for (mp_bitcnt_t k = 0; k < factors.f[j].e; k++)
       {
         lbuf_putc (' ');
-        lbuf_putmpz (factors.p[j]);
-        if (print_exponents && factors.e[j] > 1)
+        lbuf_putmpz (factors.f[j].p);
+        if (print_exponents && factors.f[j].e > 1)
           {
             lbuf_putc ('^');
-            lbuf_putint (factors.e[j]);
+            lbuf_putbitcnt (factors.f[j].e);
             break;
           }
       }
@@ -2593,6 +1935,7 @@ print_factors (char const *input)
   return true;
 }
 
+/* Output a usage diagnostic and exit with STATUS.  */
 void
 usage (int status)
 {
@@ -2619,6 +1962,7 @@ are specified on the command line, read them from standard input.\n\
   exit (status);
 }
 
+/* Read numbers from stdin, one per line, and output their factors.  */
 static bool
 do_stdin (void)
 {
@@ -2680,10 +2024,6 @@ main (int argc, char **argv)
 
   atexit (lbuf_flush);
 
-#if STAT_SQUFOF
-  memset (q_freq, 0, sizeof (q_freq));
-#endif
-
   bool ok;
   if (argc <= optind)
     ok = do_stdin ();
@@ -2695,20 +2035,5 @@ main (int argc, char **argv)
           ok = false;
     }
 
-#if STAT_SQUFOF
-  if (q_freq[0] > 0)
-    {
-      double acc_f;
-      printf ("q  freq.  cum. freq.(total: %d)\n", q_freq[0]);
-      for (int i = 1, acc_f = 0.0; i <= Q_FREQ_SIZE; i++)
-        {
-          double f = (double) q_freq[i] / q_freq[0];
-          acc_f += f;
-          printf ("%s%d %.2f%% %.2f%%\n", i == Q_FREQ_SIZE ? ">=" : "", i,
-                  100.0 * f, 100.0 * acc_f);
-        }
-    }
-#endif
-
   return ok ? EXIT_SUCCESS : EXIT_FAILURE;
 }
diff --git a/src/make-prime-list.c b/src/make-prime-list.c
index 6b53a624e..2aef09106 100644
--- a/src/make-prime-list.c
+++ b/src/make-prime-list.c
@@ -108,7 +108,7 @@ print_wide_uint (wide_uint n, int nesting, unsigned wide_uint_bits)
     }
   else if (nesting)
     {
-      printf ("(wide_uint) ");
+      printf ("(mp_limb_t) ");
       hex_digits_per_literal
         = ((wide_uint_bits - 1) % bits_per_literal) % 4 + 1;
     }
@@ -139,20 +139,21 @@ output_primes (const struct prime *primes, unsigned nprimes)
   puts ("/* Generated file -- DO NOT EDIT */\n");
   printf ("#define WIDE_UINT_BITS %u\n", wide_uint_bits);
 
-  for (i = 0, p = 2; i < nprimes; i++)
+  for (i = 0; i < nprimes; i++)
     {
-      unsigned int d8 = i + 8 < nprimes ? primes[i + 8].p - primes[i].p : 0xff;
-      if (255 < d8) /* this happens at 668221 */
+      /* Check that primes[i].p fits in int_least16_t on all platforms,
+         and that its square fits in int_least32_t on all platforms,
+         as factor.c relies on this.  */
+      if ((1 << (16 - 1)) <= primes[i].p)
         abort ();
-      printf ("P (%u, %u,\n   (", primes[i].p - p, d8);
+
+      printf ("P (%u,\n   (", primes[i].p);
       print_wide_uint (primes[i].pinv, 0, wide_uint_bits);
-      printf ("),\n   UINTMAX_MAX / %u)\n", primes[i].p);
-      p = primes[i].p;
+      printf ("),\n   (mp_limb_t) -1 / %u)\n", primes[i].p);
     }
 
-  printf ("\n#undef FIRST_OMITTED_PRIME\n");
-
   /* Find next prime */
+  p = primes[nprimes - 1].p;
   do
     {
       p += 2;
@@ -169,7 +170,7 @@ output_primes (const struct prime *primes, unsigned nprimes)
     }
   while (!is_prime);
 
-  printf ("#define FIRST_OMITTED_PRIME %u\n", p);
+  printf ("#define SQUARE_OF_FIRST_OMITTED_PRIME %u\n", p * p);
 }
 
 ATTRIBUTE_MALLOC
diff --git a/tests/factor/create-test.sh b/tests/factor/create-test.sh
index fc84b1484..7fad40cdc 100755
--- a/tests/factor/create-test.sh
+++ b/tests/factor/create-test.sh
@@ -27,6 +27,10 @@ t2=170141183460469229545748130981302223887
 # https://bugs.gnu.org/73474
 bug73474=22222222222222222202111121111
 
+# 2^400 - 593
+bigprime=25822498780869085896559191720030118743297057928292235128306593565406\
+47622016841194629645353280137831435903171972747492783
+
 # Each test is a triple: lo, hi, sha1 of result.
 # The test script, run.sh, runs seq lo hi|factor|sha1sum
 # and verifies that the actual and expected checksums are the same.
@@ -70,6 +74,7 @@ case $t in
   t35) set   ${q}958336   ${q}960335 2374919a89196e1fce93adfe779cb4664556d4b6 ;;
   t36) set   ${q}960336   ${q}962335 569e4363e8d9e8830a187d9ab27365eef08abde1 ;;
   t37) set    $bug73474    $bug73474 61d04aaf757acc5a37eb1d5581a98eea78ef50e8 ;;
+  t38) set    $bigprime    $bigprime 02f3c51a2896ff4524fd76de5f5854029879a179 ;;
   *)
     echo "$0: error: unknown test: '$test_name' -> '$t'" >&2
     exit 1
diff --git a/tests/factor/factor.pl b/tests/factor/factor.pl
index a4edb7acb..782deefe6 100755
--- a/tests/factor/factor.pl
+++ b/tests/factor/factor.pl
@@ -86,12 +86,15 @@ my @Tests =
       '115792089237316195423570985008687907853'
       . '269984665640564039457584007913129639936',
       {OUT => '2 'x255 . '2'}],
-     ['bug-gmp-plus_2_sup_128_plus_1',
+     ['bug-gmp-plus_2_sup_127_plus_1',
       '+170141183460469231731687303715884105729',
       {OUT => '3 56713727820156410577229101238628035243'}],
      ['bug-with-128-bit-uintmax_t',
       '340282366920938463463374607431768211355',
       {OUT => '5 31 2195370109167344925570158757624311041'}],
+     ['lt-5000-times-128-bit',
+      '1701071552237771378853409662551409288273703',
+      {OUT => '4999 340282366920938463463374607431768211297'}],
      ['h-1', '-h 3000', {OUT => '2^3 3 5^3'}],
      ['h-2', '3000 --exponents', {OUT => '2^3 3 5^3'}],
     );
diff --git a/tests/local.mk b/tests/local.mk
index 03114f759..dd07032da 100644
--- a/tests/local.mk
+++ b/tests/local.mk
@@ -765,7 +765,7 @@ factor_tests = \
   $(tf)/t20.sh $(tf)/t21.sh $(tf)/t22.sh $(tf)/t23.sh $(tf)/t24.sh \
   $(tf)/t25.sh $(tf)/t26.sh $(tf)/t27.sh $(tf)/t28.sh $(tf)/t29.sh \
   $(tf)/t30.sh $(tf)/t31.sh $(tf)/t32.sh $(tf)/t33.sh $(tf)/t34.sh \
-  $(tf)/t35.sh $(tf)/t36.sh $(tf)/t37.sh
+  $(tf)/t35.sh $(tf)/t36.sh $(tf)/t37.sh $(tf)/t38.sh
 
 $(factor_tests): $(tf)/run.sh $(tf)/create-test.sh
 	$(AM_V_GEN)$(MKDIR_P) $(tf)