factor: don't ever declare composites to be prime

The multiple-precision factoring code (with HAVE_GMP) was copied from
a now-obsolete version of GMP that did not pass proper arguments to
the mpz_probab_prime_p function.  It makes that code perform no more
than 3 Miller-Rabin tests only, which is not sufficient.

A Miller-Rabin test will detect composites with at least a probability
of 3/4.  For a uniform random composite, the probability will actually
be much higher.

Or put another way, of the N-3 possible Miller-Rabin tests for checking
the composite N, there is no number N for which more than (N-3)/4 of the
tests will fail to detect the number as a composite.  For most numbers N
the number of "false witnesses" will be much, much lower.

Problem numbers are of the form N=pq, p,q prime and (p-1)/(q-1) = s,
where s is a small integer.  (There are other problem forms too,
involving 3 or more prime factors.)  When s = 2, we get the 3/4 factor.

It is easy to find numbers of that form that cause coreutils' factor to
fail:

  465658903
  2242724851
  6635692801
  17709149503
  17754345703
  20889169003
  42743470771
  54890944111
  72047131003
  85862644003
  98275842811
  114654168091
  117225546301
  ...

There are 9008992 composites of the form with s=2 below 2^64.  With 3
Miller-Rabin tests, one would expect about 9008992/64 = 140766 to be
invalidly recognized as primes in that range.

* src/factor.c (MR_REPS): Define to 25.
(factor_using_pollard_rho): Use MR_REPS, not 3.
(print_factors_multi): Likewise.
* THANKS.in: Remove my name, now that it will be automatically
included in the generated THANKS file.

2 files changed, 6 insertions, 4 deletions

diff --git a/THANKS.in b/THANKS.in
index 158015187..2c3f83cdc 100644
--- a/THANKS.in
+++ b/THANKS.in
@@ -608,7 +608,6 @@ Tony Leneis                         tony@plaza.ds.adp.com
 Tony Robinson                       ajr@eng.cam.ac.uk
 Toomas Soome                        Toomas.Soome@Elion.ee
 Toralf Förster                      toralf.foerster@gmx.de
-Torbjorn Granlund                   tege@nada.kth.se
 Torbjorn Lindgren                   tl@funcom.no
 Torsten Landschoff                  torsten@pclab.ifg.uni-kiel.de
 Travis Gummels                      tgummels@redhat.com
diff --git a/src/factor.c b/src/factor.c
index 1d5580507..e63e0e01d 100644
--- a/src/factor.c
+++ b/src/factor.c
@@ -153,6 +153,9 @@ factor_using_division (mpz_t t, unsigned int limit)
   mpz_clear (r);
 }
 
+/* The number of Miller-Rabin tests we require.  */
+enum { MR_REPS = 25 };
+
 static void
 factor_using_pollard_rho (mpz_t n, int a_int)
 {
@@ -222,7 +225,7 @@ S4:
 
       mpz_div (n, n, g);	/* divide by g, before g is overwritten */
 
-      if (!mpz_probab_prime_p (g, 3))
+      if (!mpz_probab_prime_p (g, MR_REPS))
         {
           do
             {
@@ -242,7 +245,7 @@ S4:
       mpz_mod (x, x, n);
       mpz_mod (x1, x1, n);
       mpz_mod (y, y, n);
-      if (mpz_probab_prime_p (n, 3))
+      if (mpz_probab_prime_p (n, MR_REPS))
         {
           emit_factor (n);
           break;
@@ -411,7 +414,7 @@ print_factors_multi (mpz_t t)
       if (mpz_cmp_ui (t, 1) != 0)
         {
           debug ("[is number prime?] ");
-          if (mpz_probab_prime_p (t, 3))
+          if (mpz_probab_prime_p (t, MR_REPS))
             emit_factor (t);
           else
             factor_using_pollard_rho (t, 1);