aboutsummaryrefslogtreecommitdiff
path: root/gmp-6.3.0/demos/factorize.c
diff options
context:
space:
mode:
Diffstat (limited to 'gmp-6.3.0/demos/factorize.c')
-rw-r--r--gmp-6.3.0/demos/factorize.c447
1 files changed, 447 insertions, 0 deletions
diff --git a/gmp-6.3.0/demos/factorize.c b/gmp-6.3.0/demos/factorize.c
new file mode 100644
index 0000000..91e6455
--- /dev/null
+++ b/gmp-6.3.0/demos/factorize.c
@@ -0,0 +1,447 @@
+/* Factoring with Pollard's rho method.
+
+Copyright 1995, 1997-2003, 2005, 2009, 2012, 2015 Free Software
+Foundation, Inc.
+
+This program is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License as published by the Free Software
+Foundation; either version 3 of the License, or (at your option) any later
+version.
+
+This program is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
+PARTICULAR PURPOSE. See the GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along with
+this program. If not, see https://www.gnu.org/licenses/. */
+
+
+#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;
+
+/* 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_init (struct factors *factors)
+{
+ 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] ");
+ }
+
+ mpz_init (q);
+
+ p = mpz_scan1 (t, 0);
+ mpz_fdiv_q_2exp (t, t, p);
+ while (p)
+ {
+ factor_insert_ui (factors, 2);
+ --p;
+ }
+
+ p = 3;
+ for (i = 1; i <= PRIMES_PTAB_ENTRIES;)
+ {
+ if (! mpz_divisible_ui_p (t, p))
+ {
+ p += primes_diff[i++];
+ if (mpz_cmp_ui (t, p * p) < 0)
+ break;
+ }
+ else
+ {
+ mpz_tdiv_q_ui (t, t, p);
+ factor_insert_ui (factors, p);
+ }
+ }
+
+ 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. */
+ is_prime = (r == MR_REPS - 1);
+ }
+
+ if (is_prime)
+ goto ret1;
+
+ mpz_add_ui (a, a, primes_diff[r]); /* Establish new base. */
+
+ if (!mp_millerrabin (n, nm1, a, tmp, q, k))
+ {
+ is_prime = 0;
+ goto ret1;
+ }
+ }
+
+ fprintf (stderr, "Lucas prime test failure. This should not happen\n");
+ abort ();
+
+ ret1:
+ if (flag_prove_primality)
+ factor_clear (&factors);
+ ret2:
+ mpz_clears (q, a, nm1, tmp, NULL);
+
+ return is_prime;
+}
+
+void
+factor_using_pollard_rho (mpz_t n, unsigned long a, struct factors *factors)
+{
+ mpz_t x, z, y, P;
+ mpz_t t, t2;
+ unsigned long long k, l, i;
+
+ if (flag_verbose > 0)
+ {
+ printf ("[pollard-rho (%lu)] ", a);
+ }
+
+ mpz_inits (t, t2, NULL);
+ mpz_init_set_si (y, 2);
+ mpz_init_set_si (x, 2);
+ mpz_init_set_si (z, 2);
+ mpz_init_set_ui (P, 1);
+ k = 1;
+ l = 1;
+
+ while (mpz_cmp_ui (n, 1) != 0)
+ {
+ for (;;)
+ {
+ do
+ {
+ mpz_mul (t, x, x);
+ mpz_mod (x, t, n);
+ mpz_add_ui (x, x, a);
+
+ mpz_sub (t, z, x);
+ mpz_mul (t2, P, t);
+ mpz_mod (P, t2, n);
+
+ if (k % 32 == 1)
+ {
+ mpz_gcd (t, P, n);
+ if (mpz_cmp_ui (t, 1) != 0)
+ goto factor_found;
+ mpz_set (y, x);
+ }
+ }
+ while (--k != 0);
+
+ mpz_set (z, x);
+ k = l;
+ l = 2 * l;
+ for (i = 0; i < k; i++)
+ {
+ mpz_mul (t, x, x);
+ mpz_mod (x, t, n);
+ mpz_add_ui (x, x, a);
+ }
+ mpz_set (y, x);
+ }
+
+ factor_found:
+ do
+ {
+ mpz_mul (t, y, y);
+ mpz_mod (y, t, n);
+ mpz_add_ui (y, y, a);
+
+ mpz_sub (t, z, y);
+ mpz_gcd (t, t, n);
+ }
+ while (mpz_cmp_ui (t, 1) == 0);
+
+ mpz_divexact (n, n, t); /* divide by t, before t is overwritten */
+
+ if (!mp_prime_p (t))
+ {
+ if (flag_verbose > 0)
+ {
+ printf ("[composite factor--restarting pollard-rho] ");
+ }
+ factor_using_pollard_rho (t, a + 1, factors);
+ }
+ else
+ {
+ factor_insert (factors, t);
+ }
+
+ if (mp_prime_p (n))
+ {
+ factor_insert (factors, n);
+ break;
+ }
+
+ mpz_mod (x, x, n);
+ mpz_mod (z, z, n);
+ mpz_mod (y, y, n);
+ }
+
+ mpz_clears (P, t2, t, z, x, y, NULL);
+}
+
+void
+factor (mpz_t t, struct factors *factors)
+{
+ factor_init (factors);
+
+ if (mpz_sgn (t) != 0)
+ {
+ factor_using_division (t, factors);
+
+ if (mpz_cmp_ui (t, 1) != 0)
+ {
+ if (flag_verbose > 0)
+ {
+ printf ("[is number prime?] ");
+ }
+ if (mp_prime_p (t))
+ factor_insert (factors, t);
+ else
+ factor_using_pollard_rho (t, 1, factors);
+ }
+ }
+}
+
+int
+main (int argc, char *argv[])
+{
+ mpz_t t;
+ int i, j, k;
+ struct factors factors;
+
+ while (argc > 1)
+ {
+ if (!strcmp (argv[1], "-v"))
+ flag_verbose = 1;
+ else if (!strcmp (argv[1], "-w"))
+ flag_prove_primality = 0;
+ else
+ break;
+
+ argv++;
+ argc--;
+ }
+
+ mpz_init (t);
+ if (argc > 1)
+ {
+ for (i = 1; i < argc; i++)
+ {
+ mpz_set_str (t, argv[i], 0);
+
+ gmp_printf ("%Zd:", t);
+ factor (t, &factors);
+
+ for (j = 0; j < factors.nfactors; j++)
+ for (k = 0; k < factors.e[j]; k++)
+ gmp_printf (" %Zd", factors.p[j]);
+
+ puts ("");
+ factor_clear (&factors);
+ }
+ }
+ else
+ {
+ for (;;)
+ {
+ mpz_inp_str (t, stdin, 0);
+ if (feof (stdin))
+ break;
+
+ gmp_printf ("%Zd:", t);
+ factor (t, &factors);
+
+ for (j = 0; j < factors.nfactors; j++)
+ for (k = 0; k < factors.e[j]; k++)
+ gmp_printf (" %Zd", factors.p[j]);
+
+ puts ("");
+ factor_clear (&factors);
+ }
+ }
+
+ exit (0);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-09-19 17:03:21 +0000