[Gmp-commit] /var/hg/gmp: Rewrite no more current form. Implement Lucas prim...

mercurial at gmplib.org mercurial at gmplib.org
Mon Sep 10 13:42:47 CEST 2012


details:   /var/hg/gmp/rev/069dbc6582cc
changeset: 15079:069dbc6582cc
user:      Torbjorn Granlund <tege at gmplib.org>
date:      Mon Sep 10 13:42:44 2012 +0200
description:
Rewrite no more current form.  Implement Lucas prime proving, and make its use the default.

diffstat:

 ChangeLog         |    6 +
 demos/factorize.c |  503 ++++++++++++++++++++++++++++--------------------
 demos/primes.h    |  552 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 854 insertions(+), 207 deletions(-)

diffs (truncated from 1209 to 300 lines):

diff -r 7677276bdf92 -r 069dbc6582cc ChangeLog
--- a/ChangeLog	Fri Aug 24 18:43:40 2012 +0200
+++ b/ChangeLog	Mon Sep 10 13:42:44 2012 +0200
@@ -1,3 +1,9 @@
+2012-09-10  Torbjorn Granlund  <tege at gmplib.org>
+
+	* demos/factorize.c: Rewrite no more current form.  Implement Lucas
+	prime proving, and make its use the default.
+	* demos/primes.h: New file.
+
 2012-08-24  Torbjorn Granlund  <tege at gmplib.org>
 
 	* demos/factorize.c: Overhaul.
diff -r 7677276bdf92 -r 069dbc6582cc demos/factorize.c
--- a/demos/factorize.c	Fri Aug 24 18:43:40 2012 +0200
+++ b/demos/factorize.c	Mon Sep 10 13:42:44 2012 +0200
@@ -1,6 +1,6 @@
 /* Factoring with Pollard's rho method.
 
-Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2009
+Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2005, 2009, 2012
 Free Software Foundation, Inc.
 
 This program is free software; you can redistribute it and/or modify it under
@@ -19,138 +19,268 @@
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
+#include <inttypes.h>
 
 #include "gmp.h"
 
+static unsigned char primes_diff[] = {
+#define P(a,b,c) a,
+#include "primes.h"
+#undef P
+};
+#define PRIMES_PTAB_ENTRIES (sizeof(primes_diff) / sizeof(primes_diff[0]))
+
 int flag_verbose = 0;
 
-static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
+/* Prove primality or run probabilistic tests.  */
+int flag_prove_primality = 1;
+
+/* Number of Miller-Rabin tests to run when not proving primality. */
+#define MR_REPS 25
+
+struct factors
+{
+  mpz_t         *p;
+  unsigned long *e;
+  long nfactors;
+};
+
+void factor (mpz_t, struct factors *);
 
 void
-factor_using_division (mpz_t t, unsigned int limit)
+factor_init (struct factors *factors)
 {
-  mpz_t q, r;
-  unsigned long int f;
-  int ai;
-  unsigned *addv = add;
-  unsigned int failures;
+  factors->p = malloc (1);
+  factors->e = malloc (1);
+  factors->nfactors = 0;
+}
+
+void
+factor_clear (struct factors *factors)
+{
+  int i;
+
+  for (i = 0; i < factors->nfactors; i++)
+    mpz_clear (factors->p[i]);
+
+  free (factors->p);
+  free (factors->e);
+}
+
+void
+factor_insert (struct factors *factors, mpz_t prime)
+{
+  long    nfactors  = factors->nfactors;
+  mpz_t         *p  = factors->p;
+  unsigned long *e  = factors->e;
+  long i, j;
+
+  /* Locate position for insert new or increment e.  */
+  for (i = nfactors - 1; i >= 0; i--)
+    {
+      if (mpz_cmp (p[i], prime) <= 0)
+	break;
+    }
+
+  if (i < 0 || mpz_cmp (p[i], prime) != 0)
+    {
+      p = realloc (p, (nfactors + 1) * sizeof p[0]);
+      e = realloc (e, (nfactors + 1) * sizeof e[0]);
+
+      mpz_init (p[nfactors]);
+      for (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;
+
+      factors->p = p;
+      factors->e = e;
+      factors->nfactors = nfactors + 1;
+    }
+  else
+    {
+      e[i] += 1;
+    }
+}
+
+void
+factor_insert_ui (struct factors *factors, unsigned long prime)
+{
+  mpz_t pz;
+
+  mpz_init_set_ui (pz, prime);
+  factor_insert (factors, pz);
+  mpz_clear (pz);
+}
+
+
+void
+factor_using_division (mpz_t t, struct factors *factors)
+{
+  mpz_t q;
+  unsigned long int p;
+  int i;
 
   if (flag_verbose > 0)
     {
-      printf ("[trial division (%u)] ", limit);
-      fflush (stdout);
+      printf ("[trial division] ");
     }
 
   mpz_init (q);
-  mpz_init (r);
 
-  f = mpz_scan1 (t, 0);
-  mpz_div_2exp (t, t, f);
-  while (f)
+  p = mpz_scan1 (t, 0);
+  mpz_div_2exp (t, t, p);
+  while (p)
     {
-      printf ("2 ");
-      fflush (stdout);
-      --f;
+      factor_insert_ui (factors, 2);
+      --p;
     }
 
-  for (;;)
+  p = 3;
+  for (i = 1; i <= PRIMES_PTAB_ENTRIES;)
     {
-      mpz_tdiv_qr_ui (q, r, t, 3);
-      if (mpz_cmp_ui (r, 0) != 0)
-	break;
-      mpz_set (t, q);
-      printf ("3 ");
-      fflush (stdout);
-    }
-
-  for (;;)
-    {
-      mpz_tdiv_qr_ui (q, r, t, 5);
-      if (mpz_cmp_ui (r, 0) != 0)
-	break;
-      mpz_set (t, q);
-      printf ("5 ");
-      fflush (stdout);
-    }
-
-  failures = 0;
-  f = 7;
-  ai = 0;
-  while (mpz_cmp_ui (t, 1) != 0)
-    {
-      mpz_tdiv_qr_ui (q, r, t, f);
-      if (mpz_cmp_ui (r, 0) != 0)
+      if (! mpz_divisible_ui_p (t, p))
 	{
-	  f += addv[ai];
-	  if (mpz_cmp_ui (q, f) < 0)
-	    break;
-	  ai = (ai + 1) & 7;
-	  failures++;
-	  if (failures > limit)
+	  p += primes_diff[i++];
+	  if (mpz_cmp_ui (t, p * p) < 0)
 	    break;
 	}
       else
 	{
-	  mpz_swap (t, q);
-	  printf ("%lu ", f);
-	  fflush (stdout);
-	  failures = 0;
+	  mpz_tdiv_q_ui (t, t, p);
+	  factor_insert_ui (factors, p);
 	}
     }
 
-  mpz_clears (q, r, NULL);
+  mpz_clear (q);
+}
+
+static int
+mp_millerrabin (mpz_srcptr n, mpz_srcptr nm1, mpz_ptr x, mpz_ptr y,
+		mpz_srcptr q, unsigned long int k)
+{
+  unsigned long int i;
+
+  mpz_powm (y, x, q, n);
+
+  if (mpz_cmp_ui (y, 1) == 0 || mpz_cmp (y, nm1) == 0)
+    return 1;
+
+  for (i = 1; i < k; i++)
+    {
+      mpz_powm_ui (y, y, 2, n);
+      if (mpz_cmp (y, nm1) == 0)
+	return 1;
+      if (mpz_cmp_ui (y, 1) == 0)
+	return 0;
+    }
+  return 0;
+}
+
+int
+mp_prime_p (mpz_t n)
+{
+  int k, r, is_prime;
+  mpz_t q, a, nm1, tmp;
+  struct factors factors;
+
+  if (mpz_cmp_ui (n, 1) <= 0)
+    return 0;
+
+  /* We have already casted out small primes. */
+  if (mpz_cmp_ui (n, (long) FIRST_OMITTED_PRIME * FIRST_OMITTED_PRIME) < 0)
+    return 1;
+
+  mpz_inits (q, a, nm1, tmp, NULL);
+
+  /* Precomputation for Miller-Rabin.  */
+  mpz_sub_ui (nm1, n, 1);
+
+  /* Find q and k, where q is odd and n = 1 + 2**k * q.  */
+  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 = 0;
+      goto ret2;
+    }
+
+  if (flag_prove_primality)
+    {
+      /* Factor n-1 for Lucas.  */
+      mpz_set (tmp, nm1);
+      factor (tmp, &factors);
+    }
+
+  /* Loop until Lucas proves our number prime, or Miller-Rabin proves our
+     number composite.  */
+  for (r = 0; r < PRIMES_PTAB_ENTRIES; r++)
+    {
+      int i;
+
+      if (flag_prove_primality)
+	{
+	  is_prime = 1;
+	  for (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. */


More information about the gmp-commit mailing list