[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.1.1 and 1.1.1.2

version 1.1.1.1, 2000/01/10 15:35:22 version 1.1.1.2, 2000/09/09 14:13:19
Line 1 
Line 1 
 /* Factoring with Pollard's rho method.  /* Factoring with Pollard's rho method.
   
    Copyright (C) 1995 Free Software Foundation, Inc.     Copyright (C) 1995, 1997, 1998, 1999, 2000 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 the terms of the GNU General Public License as published by the     under the terms of the GNU General Public License as published by the
 Free Software Foundation; either version 2, or (at your option) any     Free Software 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 WARRANTY; without even the implied warranty of     WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU     MERCHANTABILITY or FITNESS FOR A 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 this program; see the file COPYING.  If not, write to the Free Software     with this program; see the file COPYING.  If not, write to the Free
 Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */     Software 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)  #if defined (__hpux) || defined (__alpha)  || defined (__svr4__) || defined (__SVR4)
      mpz_t t;  /* HPUX lacks random().  DEC OSF/1 1.2 random() returns a double.  */
      unsigned int limit;  long mrand48 ();
   static long
   random ()
 {  {
     return mrand48 ();
   }
   #else
   /* Glibc stdlib.h has "int32_t random();" which, on i386 at least, conflicts
      with a redeclaration as "long". */
   #ifndef __GLIBC__
   long random ();
   #endif
   #endif
   
   void
   factor_using_division (mpz_t t, 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 78  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 88  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
               a_int = random ();
             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);  
     /* 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;
   
     if (argc > 1 && !strcmp (argv[1], "-v"))
       {
         flag_verbose = 1;
         argv++;
         argc--;
       }
   
     p = 0;
   for (i = 1; i < argc; i++)    for (i = 1; i < argc; i++)
     {      {
       if (!strncmp (argv[i], "-Mp", 3))        if (!strncmp (argv[i], "-Mp", 3))
Line 216  main (argc, argv)
Line 333  main (argc, argv)
           mpz_init_set_ui (t, 1);            mpz_init_set_ui (t, 1);
           mpz_mul_2exp (t, t, p);            mpz_mul_2exp (t, t, p);
           mpz_sub_ui (t, t, 1);            mpz_sub_ui (t, t, 1);
           flag_mersenne = 1;  
         }          }
         else if (!strncmp (argv[i], "-2kp", 4))
           {
             p = atoi (argv[i] + 4);
             continue;
           }
       else        else
         {          {
           p = 0;  
           mpz_init_set_str (t, argv[i], 0);            mpz_init_set_str (t, argv[i], 0);
         }          }
   
       a = -1;        if (mpz_cmp_ui (t, 0) == 0)
       x0 = 3;          puts ("-");
         else
       factor (t, a, x0, p);          {
       puts ("");            factor (t, p);
             puts ("");
           }
     }      }
     exit (0);
   }
   
   void
   dmp (mpz_t x)
   {
     mpz_out_str (stdout, 10, x);
     puts ("");
 }  }
Legend:
Removed from v.1.1.1.1  
changed lines
  Added in v.1.1.1.2
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>