bug#12350: Composites identified as primes in factor.c (when HAVE_GMP)

[Jim Meyering]

Torbjorn Granlund wrote:
> The very old factoring code cut from an now obsolete version GMP does
> not pass proper arguments to the mpz_probab_prime_p function.  It ask
> for 3 Miller-Rabin tests only, which is not sufficient.

Hi Torbjorn

Thank you for the patch and explanation.
I've converted that into the commit below in your name.
Please proofread it and let me know if you'd like to change anything.
I tweaked the patch to change MR_REPS from a #define to an enum
and to add the comment just preceding.

I'll add NEWS and tests separately.


>From ea6dd126e6452504f9fa1d6708d25473e2c27e67 Mon Sep 17 00:00:00 2001
From: Torbjorn Granlund <address@hidden>
Date: Tue, 4 Sep 2012 16:22:47 +0200
Subject: [PATCH] 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
by 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 for N=pq, p,q prime and (p-1)/(q-1) = s,
where s is a small integer.  (There are other problem forms too,
incvolving 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 test, one would expect about 9008992/4^64 = 140766 to be
invalidly recognised 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.
---
 src/factor.c | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/factor.c b/src/factor.c
index 1d55805..e63e0e0 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);
--
1.7.12.176.g3fc0e4c