third_party/gmp/demos/factorize.c - RealtimeRoboticsGroup/test - Gitiles

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