/* 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);
}