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

[Jim Meyering]

Jim Meyering wrote:

> 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: Torbjorn Granlund <address@hidden>
> Date: Tue, 4 Sep 2012 16:22:47 +0200
> Subject: [PATCH] factor: don't ever declare composites to be prime

Torbjörn, I've just noticed that I misspelled your name above.

Here's the NEWS/tests addition.
Following is an adjusted commit that spells your name properly.

>From e561ff991b74dc19f6728aa1e6e61d1927055ac1 Mon Sep 17 00:00:00 2001
From: Jim Meyering <address@hidden>
Date: Tue, 4 Sep 2012 18:26:25 +0200
Subject: [PATCH] factor: doc and tests
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* NEWS (Bug fixes): Mention it.
* tests/misc/factor.pl: Add five of Torbjörn's tests.
---
 NEWS                 | 3 +++
 tests/misc/factor.pl | 5 +++++
 2 files changed, 8 insertions(+)

diff --git a/NEWS b/NEWS
index f3874fd..ffa7939 100644
--- a/NEWS
+++ b/NEWS
@@ -15,6 +15,9 @@ GNU coreutils NEWS                                    -*- 
outline -*-
   it detects this precise type of cycle, diagnoses it as such and
   eventually exits nonzero.

+  factor (when using gmp) would mistakenly declare some composite numbers
+  to be prime, e.g., 465658903, 2242724851, 6635692801.
+
   rm -i -d now prompts the user then removes an empty directory, rather
   than ignoring the -d option and failing with an 'Is a directory' error.
   [bug introduced in coreutils-8.19, with the addition of --dir (-d)]
diff --git a/tests/misc/factor.pl b/tests/misc/factor.pl
index 47f9343..38a5037 100755
--- a/tests/misc/factor.pl
+++ b/tests/misc/factor.pl
@@ -67,6 +67,11 @@ my @Tests =
       {OUT => "4: 2 2\n"},
       {ERR => "$prog: 'a' is not a valid positive integer\n"},
       {EXIT => 1}],
+     ['bug-2012-a', '465658903', {OUT => '15259 30517'}],
+     ['bug-2012-b', '2242724851', {OUT => '33487 66973'}],
+     ['bug-2012-c', '6635692801', {OUT => '57601 115201'}],
+     ['bug-2012-d', '17709149503', {OUT => '94099 188197'}],
+     ['bug-2012-e', '17754345703', {OUT => '94219 188437'}],
     );

 # Prepend the command line argument and append a newline to end
--
1.7.12.176.g3fc0e4c


>From 4c21a96443ee26eb0d4da31526ce4cf180ac7a4e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Torbj=C3=B6rn=20Granlund?= <address@hidden>
Date: Tue, 4 Sep 2012 18:38:29 +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,
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 test, 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.
---
 THANKS.in    | 1 -
 src/factor.c | 9 ++++++---
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/THANKS.in b/THANKS.in
index 1580151..2c3f83c 100644
--- a/THANKS.in
+++ b/THANKS.in
@@ -608,7 +608,6 @@ Tony Leneis                         address@hidden
 Tony Robinson                       address@hidden
 Toomas Soome                        address@hidden
 Toralf Förster                      address@hidden
-Torbjorn Granlund                   address@hidden
 Torbjorn Lindgren                   address@hidden
 Torsten Landschoff                  address@hidden
 Travis Gummels                      address@hidden
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