contrib/libgmp/demos/factorize.c

   1 /* Factoring with Pollard's rho method.
   2
   3    Copyright (C) 1995 Free Software Foundation, Inc.
   4
   5 This program is free software; you can redistribute it and/or modify it
   6 under the terms of the GNU General Public License as published by the
   7 Free Software Foundation; either version 2, or (at your option) any
   8 later version.
   9
  10 This program is distributed in the hope that it will be useful, but
  11 WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License along
  16 with this program; see the file COPYING.  If not, write to the Free Software
  17 Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
  18
  19 #include <stdio.h>
  20 #include "gmp.h"
  21
  22 int flag_mersenne = 0;
  23
  24 static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
  25
  26 factor_using_division (t, limit)
  27      mpz_t t;
  28      unsigned int limit;
  29 {
  30   mpz_t q, r;
  31   unsigned long int f;
  32   int i, ai;
  33   unsigned *addv = add;
  34
  35   mpz_init (q);
  36   mpz_init (r);
  37
  38   if (mpz_probab_prime_p (t, 50))
  39     goto ready;
  40
  41   for (;;)
  42     {
  43       mpz_tdiv_qr_ui (q, r, t, 2);
  44       if (mpz_cmp_ui (r, 0) != 0)
  45         break;
  46       mpz_set (t, q);
  47       printf ("2 ");
  48       fflush (stdout);
  49       if (mpz_probab_prime_p (t, 50))
  50         goto ready;
  51     }
  52
  53   for (;;)
  54     {
  55       mpz_tdiv_qr_ui (q, r, t, 3);
  56       if (mpz_cmp_ui (r, 0) != 0)
  57         break;
  58       mpz_set (t, q);
  59       printf ("3 ");
  60       fflush (stdout);
  61       if (mpz_probab_prime_p (t, 50))
  62         goto ready;
  63     }
  64
  65   for (;;)
  66     {
  67       mpz_tdiv_qr_ui (q, r, t, 5);
  68       if (mpz_cmp_ui (r, 0) != 0)
  69         break;
  70       mpz_set (t, q);
  71       printf ("5 ");
  72       fflush (stdout);
  73       if (mpz_probab_prime_p (t, 50))
  74         goto ready;
  75     }
  76
  77   f = 7;
  78   ai = 0;
  79   for (;;)
  80     {
  81       mpz_tdiv_qr_ui (q, r, t, f);
  82       if (mpz_cmp_ui (r, 0) != 0)
  83         {
  84           f += addv[ai];
  85           if (f > limit)
  86             goto ret;
  87           ai = (ai + 1) & 7;
  88         }
  89       else
  90         {
  91           mpz_set (t, q);
  92           printf ("%lu ", f);
  93           fflush (stdout);
  94           if (mpz_probab_prime_p (t, 50))
  95             goto ready;
  96         }
  97     }
  98
  99  ready:
 100   mpz_out_str (stdout, 10, t);
 101   fflush (stdout);
 102   mpz_set_ui (t, 1);
 103   fputc (' ', stdout);
 104  ret:
 105   mpz_clear (q);
 106   mpz_clear (r);
 107 }
 108
 109 void
 110 factor_using_pollard_rho (m, a_int, x0, p)
 111      mpz_t m;
 112      long a_int;
 113      long x0;
 114      unsigned long p;
 115 {
 116   mpz_t x, y, q;
 117   mpz_t a;
 118   mpz_t d;
 119   mpz_t tmp;
 120   mpz_t n;
 121   int i = 1;
 122   int j = 1;
 123
 124   mpz_init_set (n, m);
 125
 126   mpz_init (d);
 127   mpz_init_set_ui (q, 1);
 128   mpz_init (tmp);
 129
 130   mpz_init_set_si (a, a_int);
 131   mpz_init_set_si (x, x0);
 132   mpz_init_set_si (y, x0);
 133
 134   while (mpz_cmp_ui (n, 1) != 0)
 135     {
 136       if (flag_mersenne)
 137         {
 138           mpz_powm_ui (x, x, p, n);  mpz_add (x, x, a);
 139           mpz_powm_ui (y, y, p, n);  mpz_add (y, y, a);
 140           mpz_powm_ui (y, y, p, n);  mpz_add (y, y, a);
 141         }
 142       else
 143         {
 144           mpz_mul (x, x, x);  mpz_add (x, x, a); mpz_mod (x, x, n);
 145           mpz_mul (y, y, y);  mpz_add (y, y, a); mpz_mod (y, y, n);
 146           mpz_mul (y, y, y);  mpz_add (y, y, a); mpz_mod (y, y, n);
 147         }
 148
 149       if (mpz_cmp (x, y) > 0)
 150         mpz_sub (tmp, x, y);
 151       else
 152         mpz_sub (tmp, y, x);
 153       mpz_mul (q, q, tmp);
 154       mpz_mod (q, q, n);
 155
 156       if (++i % j == 0)
 157         {
 158           j += 1;
 159           mpz_gcd (d, q, n);
 160           if (mpz_cmp_ui (d, 1) != 0)
 161             {
 162               if (!mpz_probab_prime_p (d, 50))
 163                 factor_using_pollard_rho (d, (random () & 31) - 16,
 164                                           (random () & 31), p);
 165               else
 166                 {
 167                   mpz_out_str (stdout, 10, d);
 168                   fflush (stdout);
 169                   fputc (' ', stdout);
 170                 }
 171               mpz_div (n, n, d);
 172               if (mpz_probab_prime_p (n, 50))
 173                 {
 174                   mpz_out_str (stdout, 10, n);
 175                   fflush (stdout);
 176                   fputc (' ', stdout);
 177                   break;
 178                 }
 179             }
 180         }
 181     }
 182
 183   mpz_clear (n);
 184   mpz_clear (d);
 185   mpz_clear (q);
 186   mpz_clear (tmp);
 187   mpz_clear (a);
 188   mpz_clear (x);
 189   mpz_clear (y);
 190 }
 191
 192 factor (t, a, x0, p)
 193      mpz_t t;
 194      long a;
 195      long x0;
 196      unsigned long p;
 197 {
 198   factor_using_division (t, 1000000);
 199   factor_using_pollard_rho (t, a, x0, p);
 200 }
 201
 202 main (argc, argv)
 203      int argc;
 204      char *argv[];
 205 {
 206   mpz_t t;
 207   long x0, a;
 208   unsigned long p;
 209   int i;
 210
 211   for (i = 1; i < argc; i++)
 212     {
 213       if (!strncmp (argv[i], "-Mp", 3))
 214         {
 215           p = atoi (argv[i] + 3);
 216           mpz_init_set_ui (t, 1);
 217           mpz_mul_2exp (t, t, p);
 218           mpz_sub_ui (t, t, 1);
 219           flag_mersenne = 1;
 220         }
 221       else
 222         {
 223           p = 0;
 224           mpz_init_set_str (t, argv[i], 0);
 225         }
 226
 227       a = -1;
 228       x0 = 3;
 229
 230       factor (t, a, x0, p);
 231       puts ("");
 232     }
 233 }