From 11da511c784eca003deb90c23570f0873954e0de Mon Sep 17 00:00:00 2001 From: Duncan Wilkie Date: Sat, 18 Nov 2023 06:11:09 -0600 Subject: Initial commit. --- gmp-6.3.0/demos/factorize.c | 447 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 447 insertions(+) create mode 100644 gmp-6.3.0/demos/factorize.c (limited to 'gmp-6.3.0/demos/factorize.c') 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 +#include +#include +#include + +#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); +} -- cgit v1.2.3