mirror of
https://review.haiku-os.org/buildtools
synced 2025-02-07 14:34:51 +01:00
f133d007f9
Old version was from 2012-05-06, 6.1.2 is from 2016-12-16 A lot of support for newer processors and speedups since then See gmp/NEWS for details
448 lines
8.2 KiB
C
448 lines
8.2 KiB
C
/* 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);
|
|
}
|