Quelle factorize.c Sprache: 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);
}

quality96%

¤ Dauer der Verarbeitung: 0.12 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.