[BACK]Return to factorize.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib / gmp / demos

Diff for /OpenXM_contrib/gmp/demos/Attic/factorize.c between version 1.1 and 1.1.1.3

version 1.1, 2000/01/10 15:35:22 version 1.1.1.3, 2003/08/25 16:06:03
Line 1 
Line 1 
 /* Factoring with Pollard's rho method.  /* Factoring with Pollard's rho method.
   
    Copyright (C) 1995 Free Software Foundation, Inc.  Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002 Free Software Foundation,
   Inc.
   
 This program is free software; you can redistribute it and/or modify it  This program is free software; you can redistribute it and/or modify it under
 under the terms of the GNU General Public License as published by the  the terms of the GNU General Public License as published by the Free Software
 Free Software Foundation; either version 2, or (at your option) any  Foundation; either version 2, or (at your option) any later version.
 later version.  
   
 This program is distributed in the hope that it will be useful, but  This program is distributed in the hope that it will be useful, but WITHOUT ANY
 WITHOUT ANY WARRANTY; without even the implied warranty of  WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU  PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 General Public License for more details.  
   
 You should have received a copy of the GNU General Public License along  You should have received a copy of the GNU General Public License along with
 with this program; see the file COPYING.  If not, write to the Free Software  this program; see the file COPYING.  If not, write to the Free Software
 Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
   
   #include <stdlib.h>
 #include <stdio.h>  #include <stdio.h>
   #include <string.h>
   
 #include "gmp.h"  #include "gmp.h"
   
 int flag_mersenne = 0;  int flag_verbose = 0;
   
 static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};  static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
   
 factor_using_division (t, limit)  void
      mpz_t t;  factor_using_division (mpz_t t, unsigned int limit)
      unsigned int limit;  
 {  {
   mpz_t q, r;    mpz_t q, r;
   unsigned long int f;    unsigned long int f;
   int i, ai;    int ai;
   unsigned *addv = add;    unsigned *addv = add;
     unsigned int failures;
   
     if (flag_verbose)
       {
         printf ("[trial division (%u)] ", limit);
         fflush (stdout);
       }
   
   mpz_init (q);    mpz_init (q);
   mpz_init (r);    mpz_init (r);
   
   if (mpz_probab_prime_p (t, 50))    f = mpz_scan1 (t, 0);
     goto ready;    mpz_div_2exp (t, t, f);
     while (f)
   for (;;)  
     {      {
       mpz_tdiv_qr_ui (q, r, t, 2);  
       if (mpz_cmp_ui (r, 0) != 0)  
         break;  
       mpz_set (t, q);  
       printf ("2 ");        printf ("2 ");
       fflush (stdout);        fflush (stdout);
       if (mpz_probab_prime_p (t, 50))        --f;
         goto ready;  
     }      }
   
   for (;;)    for (;;)
Line 58  factor_using_division (t, limit)
Line 60  factor_using_division (t, limit)
       mpz_set (t, q);        mpz_set (t, q);
       printf ("3 ");        printf ("3 ");
       fflush (stdout);        fflush (stdout);
       if (mpz_probab_prime_p (t, 50))  
         goto ready;  
     }      }
   
   for (;;)    for (;;)
Line 70  factor_using_division (t, limit)
Line 70  factor_using_division (t, limit)
       mpz_set (t, q);        mpz_set (t, q);
       printf ("5 ");        printf ("5 ");
       fflush (stdout);        fflush (stdout);
       if (mpz_probab_prime_p (t, 50))  
         goto ready;  
     }      }
   
     failures = 0;
   f = 7;    f = 7;
   ai = 0;    ai = 0;
   for (;;)    while (mpz_cmp_ui (t, 1) != 0)
     {      {
       mpz_tdiv_qr_ui (q, r, t, f);        mpz_tdiv_qr_ui (q, r, t, f);
       if (mpz_cmp_ui (r, 0) != 0)        if (mpz_cmp_ui (r, 0) != 0)
         {          {
           f += addv[ai];            f += addv[ai];
           if (f > limit)            if (mpz_cmp_ui (q, f) < 0)
             goto ret;              break;
           ai = (ai + 1) & 7;            ai = (ai + 1) & 7;
             failures++;
             if (failures > limit)
               break;
         }          }
       else        else
         {          {
           mpz_set (t, q);            mpz_swap (t, q);
           printf ("%lu ", f);            printf ("%lu ", f);
           fflush (stdout);            fflush (stdout);
           if (mpz_probab_prime_p (t, 50))            failures = 0;
             goto ready;  
         }          }
     }      }
   
  ready:  
   mpz_out_str (stdout, 10, t);  
   fflush (stdout);  
   mpz_set_ui (t, 1);  
   fputc (' ', stdout);  
  ret:  
   mpz_clear (q);    mpz_clear (q);
   mpz_clear (r);    mpz_clear (r);
 }  }
   
 void  void
 factor_using_pollard_rho (m, a_int, x0, p)  factor_using_division_2kp (mpz_t t, unsigned int limit, unsigned long p)
      mpz_t m;  
      long a_int;  
      long x0;  
      unsigned long p;  
 {  {
   mpz_t x, y, q;    mpz_t r;
     mpz_t f;
     unsigned int k;
   
     mpz_init (r);
     mpz_init_set_ui (f, 2 * p);
     mpz_add_ui (f, f, 1);
     for (k = 1; k < limit; k++)
       {
         mpz_tdiv_r (r, t, f);
         while (mpz_cmp_ui (r, 0) == 0)
           {
             mpz_tdiv_q (t, t, f);
             mpz_tdiv_r (r, t, f);
             mpz_out_str (stdout, 10, f);
             fflush (stdout);
             fputc (' ', stdout);
           }
         mpz_add_ui (f, f, 2 * p);
       }
   
     mpz_clear (f);
     mpz_clear (r);
   }
   
   void
   factor_using_pollard_rho (mpz_t n, int a_int, unsigned long p)
   {
     mpz_t x, x1, y, P;
   mpz_t a;    mpz_t a;
   mpz_t d;    mpz_t g;
   mpz_t tmp;    mpz_t t1, t2;
   mpz_t n;    int k, l, c, i;
   int i = 1;  
   int j = 1;  
   
   mpz_init_set (n, m);    if (flag_verbose)
       {
         printf ("[pollard-rho (%d)] ", a_int);
         fflush (stdout);
       }
   
   mpz_init (d);    mpz_init (g);
   mpz_init_set_ui (q, 1);    mpz_init (t1);
   mpz_init (tmp);    mpz_init (t2);
   
   mpz_init_set_si (a, a_int);    mpz_init_set_si (a, a_int);
   mpz_init_set_si (x, x0);    mpz_init_set_si (y, 2);
   mpz_init_set_si (y, x0);    mpz_init_set_si (x, 2);
     mpz_init_set_si (x1, 2);
     k = 1;
     l = 1;
     mpz_init_set_ui (P, 1);
     c = 0;
   
   while (mpz_cmp_ui (n, 1) != 0)    while (mpz_cmp_ui (n, 1) != 0)
     {      {
       if (flag_mersenne)  S2:
         if (p != 0)
         {          {
           mpz_powm_ui (x, x, p, n);  mpz_add (x, x, a);            mpz_powm_ui (x, x, p, n); mpz_add (x, x, a);
           mpz_powm_ui (y, y, p, n);  mpz_add (y, y, a);  
           mpz_powm_ui (y, y, p, n);  mpz_add (y, y, a);  
         }          }
       else        else
         {          {
           mpz_mul (x, x, x);  mpz_add (x, x, a); mpz_mod (x, x, n);            mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
           mpz_mul (y, y, y);  mpz_add (y, y, a); mpz_mod (y, y, n);  
           mpz_mul (y, y, y);  mpz_add (y, y, a); mpz_mod (y, y, n);  
         }          }
         mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
         c++;
         if (c == 20)
           {
             c = 0;
             mpz_gcd (g, P, n);
             if (mpz_cmp_ui (g, 1) != 0)
               goto S4;
             mpz_set (y, x);
           }
   S3:
         k--;
         if (k != 0)
           goto S2;
   
       if (mpz_cmp (x, y) > 0)        mpz_gcd (g, P, n);
         mpz_sub (tmp, x, y);        if (mpz_cmp_ui (g, 1) != 0)
       else          goto S4;
         mpz_sub (tmp, y, x);  
       mpz_mul (q, q, tmp);  
       mpz_mod (q, q, n);  
   
       if (++i % j == 0)        mpz_set (x1, x);
         k = l;
         l = 2 * l;
         for (i = 0; i < k; i++)
         {          {
           j += 1;            if (p != 0)
           mpz_gcd (d, q, n);  
           if (mpz_cmp_ui (d, 1) != 0)  
             {              {
               if (!mpz_probab_prime_p (d, 50))                mpz_powm_ui (x, x, p, n); mpz_add (x, x, a);
                 factor_using_pollard_rho (d, (random () & 31) - 16,  
                                           (random () & 31), p);  
               else  
                 {  
                   mpz_out_str (stdout, 10, d);  
                   fflush (stdout);  
                   fputc (' ', stdout);  
                 }  
               mpz_div (n, n, d);  
               if (mpz_probab_prime_p (n, 50))  
                 {  
                   mpz_out_str (stdout, 10, n);  
                   fflush (stdout);  
                   fputc (' ', stdout);  
                   break;  
                 }  
             }              }
             else
               {
                 mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
               }
         }          }
         mpz_set (y, x);
         c = 0;
         goto S2;
   S4:
         do
           {
             if (p != 0)
               {
                 mpz_powm_ui (y, y, p, n); mpz_add (y, y, a);
               }
             else
               {
                 mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
               }
             mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
           }
         while (mpz_cmp_ui (g, 1) == 0);
   
         if (!mpz_probab_prime_p (g, 3))
           {
             do
               {
                 mp_limb_t a_limb;
                 mpn_random (&a_limb, (mp_size_t) 1);
                 a_int = (int) a_limb;
               }
             while (a_int == -2 || a_int == 0);
   
             if (flag_verbose)
               {
                 printf ("[composite factor--restarting pollard-rho] ");
                 fflush (stdout);
               }
             factor_using_pollard_rho (g, a_int, p);
             break;
           }
         else
           {
             mpz_out_str (stdout, 10, g);
             fflush (stdout);
             fputc (' ', stdout);
           }
         mpz_div (n, n, g);
         mpz_mod (x, x, n);
         mpz_mod (x1, x1, n);
         mpz_mod (y, y, n);
         if (mpz_probab_prime_p (n, 3))
           {
             mpz_out_str (stdout, 10, n);
             fflush (stdout);
             fputc (' ', stdout);
             break;
           }
     }      }
   
   mpz_clear (n);    mpz_clear (g);
   mpz_clear (d);    mpz_clear (P);
   mpz_clear (q);    mpz_clear (t2);
   mpz_clear (tmp);    mpz_clear (t1);
   mpz_clear (a);    mpz_clear (a);
     mpz_clear (x1);
   mpz_clear (x);    mpz_clear (x);
   mpz_clear (y);    mpz_clear (y);
 }  }
   
 factor (t, a, x0, p)  void
      mpz_t t;  factor (mpz_t t, unsigned long p)
      long a;  
      long x0;  
      unsigned long p;  
 {  {
   factor_using_division (t, 1000000);    unsigned int division_limit;
   factor_using_pollard_rho (t, a, x0, p);  
     if (mpz_sgn (t) == 0)
       return;
   
     /* Set the trial division limit according the size of t.  */
     division_limit = mpz_sizeinbase (t, 2);
     if (division_limit > 1000)
       division_limit = 1000 * 1000;
     else
       division_limit = division_limit * division_limit;
   
     if (p != 0)
       factor_using_division_2kp (t, division_limit / 10, p);
     else
       factor_using_division (t, division_limit);
   
     if (mpz_cmp_ui (t, 1) != 0)
       {
         if (flag_verbose)
           {
             printf ("[is number prime?] ");
             fflush (stdout);
           }
         if (mpz_probab_prime_p (t, 3))
           mpz_out_str (stdout, 10, t);
         else
           factor_using_pollard_rho (t, 1, p);
       }
 }  }
   
 main (argc, argv)  main (int argc, char *argv[])
      int argc;  
      char *argv[];  
 {  {
   mpz_t t;    mpz_t t;
   long x0, a;  
   unsigned long p;    unsigned long p;
   int i;    int i;
   
   for (i = 1; i < argc; i++)    if (argc > 1 && !strcmp (argv[1], "-v"))
     {      {
       if (!strncmp (argv[i], "-Mp", 3))        flag_verbose = 1;
         argv++;
         argc--;
       }
   
     mpz_init (t);
     if (argc > 1)
       {
         p = 0;
         for (i = 1; i < argc; i++)
         {          {
           p = atoi (argv[i] + 3);            if (!strncmp (argv[i], "-Mp", 3))
           mpz_init_set_ui (t, 1);              {
           mpz_mul_2exp (t, t, p);                p = atoi (argv[i] + 3);
           mpz_sub_ui (t, t, 1);                mpz_set_ui (t, 1);
           flag_mersenne = 1;                mpz_mul_2exp (t, t, p);
                 mpz_sub_ui (t, t, 1);
               }
             else if (!strncmp (argv[i], "-2kp", 4))
               {
                 p = atoi (argv[i] + 4);
                 continue;
               }
             else
               {
                 mpz_set_str (t, argv[i], 0);
               }
   
             if (mpz_cmp_ui (t, 0) == 0)
               puts ("-");
             else
               {
                 factor (t, p);
                 puts ("");
               }
         }          }
       else      }
     else
       {
         for (;;)
         {          {
           p = 0;            mpz_inp_str (t, stdin, 0);
           mpz_init_set_str (t, argv[i], 0);            if (feof (stdin))
               break;
             mpz_out_str (stdout, 10, t); printf (" = ");
             factor (t, 0);
             puts ("");
         }          }
   
       a = -1;  
       x0 = 3;  
   
       factor (t, a, x0, p);  
       puts ("");  
     }      }
   
     exit (0);
 }  }
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.1.1.3
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>