[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