← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/gmp/demos/factorize.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter Van Eynde
  • Date: 2006-05-17 02:46:26 UTC
  • Revision ID: james.westby@ubuntu.com-20060517024626-lljr08ftv9g9vefl
Tags: upstream-0.9h-20060510
Import upstream version 0.9h-20060510

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/gmp/demos/factorize.c
 
1
/* Factoring with Pollard's rho method.
 
2
 
 
3
Copyright 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003 Free Software
 
4
Foundation, Inc.
 
5
 
 
6
This program is free software; you can redistribute it and/or modify it under
 
7
the terms of the GNU General Public License as published by the Free Software
 
8
Foundation; either version 2, or (at your option) any later version.
 
9
 
 
10
This program is distributed in the hope that it will be useful, but WITHOUT ANY
 
11
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
 
12
PARTICULAR PURPOSE.  See the GNU General Public License for more details.
 
13
 
 
14
You should have received a copy of the GNU General Public License along with
 
15
this program; see the file COPYING.  If not, write to the Free Software
 
16
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  */
 
17
 
 
18
#include <stdlib.h>
 
19
#include <stdio.h>
 
20
#include <string.h>
 
21
 
 
22
#include "gmp.h"
 
23
 
 
24
int flag_verbose = 0;
 
25
 
 
26
static unsigned add[] = {4, 2, 4, 2, 4, 6, 2, 6};
 
27
 
 
28
void
 
29
factor_using_division (mpz_t t, unsigned int limit)
 
30
{
 
31
  mpz_t q, r;
 
32
  unsigned long int f;
 
33
  int ai;
 
34
  unsigned *addv = add;
 
35
  unsigned int failures;
 
36
 
 
37
  if (flag_verbose)
 
38
    {
 
39
      printf ("[trial division (%u)] ", limit);
 
40
      fflush (stdout);
 
41
    }
 
42
 
 
43
  mpz_init (q);
 
44
  mpz_init (r);
 
45
 
 
46
  f = mpz_scan1 (t, 0);
 
47
  mpz_div_2exp (t, t, f);
 
48
  while (f)
 
49
    {
 
50
      printf ("2 ");
 
51
      fflush (stdout);
 
52
      --f;
 
53
    }
 
54
 
 
55
  for (;;)
 
56
    {
 
57
      mpz_tdiv_qr_ui (q, r, t, 3);
 
58
      if (mpz_cmp_ui (r, 0) != 0)
 
59
        break;
 
60
      mpz_set (t, q);
 
61
      printf ("3 ");
 
62
      fflush (stdout);
 
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
    }
 
74
 
 
75
  failures = 0;
 
76
  f = 7;
 
77
  ai = 0;
 
78
  while (mpz_cmp_ui (t, 1) != 0)
 
79
    {
 
80
      mpz_tdiv_qr_ui (q, r, t, f);
 
81
      if (mpz_cmp_ui (r, 0) != 0)
 
82
        {
 
83
          f += addv[ai];
 
84
          if (mpz_cmp_ui (q, f) < 0)
 
85
            break;
 
86
          ai = (ai + 1) & 7;
 
87
          failures++;
 
88
          if (failures > limit)
 
89
            break;
 
90
        }
 
91
      else
 
92
        {
 
93
          mpz_swap (t, q);
 
94
          printf ("%lu ", f);
 
95
          fflush (stdout);
 
96
          failures = 0;
 
97
        }
 
98
    }
 
99
 
 
100
  mpz_clear (q);
 
101
  mpz_clear (r);
 
102
}
 
103
 
 
104
void
 
105
factor_using_division_2kp (mpz_t t, unsigned int limit, unsigned long p)
 
106
{
 
107
  mpz_t r;
 
108
  mpz_t f;
 
109
  unsigned int k;
 
110
 
 
111
  mpz_init (r);
 
112
  mpz_init_set_ui (f, 2 * p);
 
113
  mpz_add_ui (f, f, 1);
 
114
  for (k = 1; k < limit; k++)
 
115
    {
 
116
      mpz_tdiv_r (r, t, f);
 
117
      while (mpz_cmp_ui (r, 0) == 0)
 
118
        {
 
119
          mpz_tdiv_q (t, t, f);
 
120
          mpz_tdiv_r (r, t, f);
 
121
          mpz_out_str (stdout, 10, f);
 
122
          fflush (stdout);
 
123
          fputc (' ', stdout);
 
124
        }
 
125
      mpz_add_ui (f, f, 2 * p);
 
126
    }
 
127
 
 
128
  mpz_clear (f);
 
129
  mpz_clear (r);
 
130
}
 
131
 
 
132
void
 
133
factor_using_pollard_rho (mpz_t n, int a_int, unsigned long p)
 
134
{
 
135
  mpz_t x, x1, y, P;
 
136
  mpz_t a;
 
137
  mpz_t g;
 
138
  mpz_t t1, t2;
 
139
  int k, l, c, i;
 
140
 
 
141
  if (flag_verbose)
 
142
    {
 
143
      printf ("[pollard-rho (%d)] ", a_int);
 
144
      fflush (stdout);
 
145
    }
 
146
 
 
147
  mpz_init (g);
 
148
  mpz_init (t1);
 
149
  mpz_init (t2);
 
150
 
 
151
  mpz_init_set_si (a, a_int);
 
152
  mpz_init_set_si (y, 2);
 
153
  mpz_init_set_si (x, 2);
 
154
  mpz_init_set_si (x1, 2);
 
155
  k = 1;
 
156
  l = 1;
 
157
  mpz_init_set_ui (P, 1);
 
158
  c = 0;
 
159
 
 
160
  while (mpz_cmp_ui (n, 1) != 0)
 
161
    {
 
162
S2:
 
163
      if (p != 0)
 
164
        {
 
165
          mpz_powm_ui (x, x, p, n); mpz_add (x, x, a);
 
166
        }
 
167
      else
 
168
        {
 
169
          mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 
170
        }
 
171
      mpz_sub (t1, x1, x); mpz_mul (t2, P, t1); mpz_mod (P, t2, n);
 
172
      c++;
 
173
      if (c == 20)
 
174
        {
 
175
          c = 0;
 
176
          mpz_gcd (g, P, n);
 
177
          if (mpz_cmp_ui (g, 1) != 0)
 
178
            goto S4;
 
179
          mpz_set (y, x);
 
180
        }
 
181
S3:
 
182
      k--;
 
183
      if (k > 0)
 
184
        goto S2;
 
185
 
 
186
      mpz_gcd (g, P, n);
 
187
      if (mpz_cmp_ui (g, 1) != 0)
 
188
        goto S4;
 
189
 
 
190
      mpz_set (x1, x);
 
191
      k = l;
 
192
      l = 2 * l;
 
193
      for (i = 0; i < k; i++)
 
194
        {
 
195
          if (p != 0)
 
196
            {
 
197
              mpz_powm_ui (x, x, p, n); mpz_add (x, x, a);
 
198
            }
 
199
          else
 
200
            {
 
201
              mpz_mul (x, x, x); mpz_add (x, x, a); mpz_mod (x, x, n);
 
202
            }
 
203
        }
 
204
      mpz_set (y, x);
 
205
      c = 0;
 
206
      goto S2;
 
207
S4:
 
208
      do
 
209
        {
 
210
          if (p != 0)
 
211
            {
 
212
              mpz_powm_ui (y, y, p, n); mpz_add (y, y, a); 
 
213
            }
 
214
          else
 
215
            {
 
216
              mpz_mul (y, y, y); mpz_add (y, y, a); mpz_mod (y, y, n);
 
217
            }
 
218
          mpz_sub (t1, x1, y); mpz_gcd (g, t1, n);
 
219
        }
 
220
      while (mpz_cmp_ui (g, 1) == 0);
 
221
 
 
222
      if (!mpz_probab_prime_p (g, 3))
 
223
        {
 
224
          do
 
225
            {
 
226
              mp_limb_t a_limb;
 
227
              mpn_random (&a_limb, (mp_size_t) 1);
 
228
              a_int = (int) a_limb;
 
229
            }
 
230
          while (a_int == -2 || a_int == 0);
 
231
 
 
232
          if (flag_verbose)
 
233
            {
 
234
              printf ("[composite factor--restarting pollard-rho] ");
 
235
              fflush (stdout);
 
236
            }
 
237
          factor_using_pollard_rho (g, a_int, p);
 
238
          break;
 
239
        }
 
240
      else
 
241
        {
 
242
          mpz_out_str (stdout, 10, g);
 
243
          fflush (stdout);
 
244
          fputc (' ', stdout);
 
245
        }
 
246
      mpz_div (n, n, g);
 
247
      mpz_mod (x, x, n);
 
248
      mpz_mod (x1, x1, n);
 
249
      mpz_mod (y, y, n);
 
250
      if (mpz_probab_prime_p (n, 3))
 
251
        {
 
252
          mpz_out_str (stdout, 10, n);
 
253
          fflush (stdout);
 
254
          fputc (' ', stdout);
 
255
          break;
 
256
        }
 
257
    }
 
258
 
 
259
  mpz_clear (g);
 
260
  mpz_clear (P);
 
261
  mpz_clear (t2);
 
262
  mpz_clear (t1);
 
263
  mpz_clear (a);
 
264
  mpz_clear (x1);
 
265
  mpz_clear (x);
 
266
  mpz_clear (y);
 
267
}
 
268
 
 
269
void
 
270
factor (mpz_t t, unsigned long p)
 
271
{
 
272
  unsigned int division_limit;
 
273
 
 
274
  if (mpz_sgn (t) == 0)
 
275
    return;
 
276
 
 
277
  /* Set the trial division limit according the size of t.  */
 
278
  division_limit = mpz_sizeinbase (t, 2);
 
279
  if (division_limit > 1000)
 
280
    division_limit = 1000 * 1000;
 
281
  else
 
282
    division_limit = division_limit * division_limit;
 
283
 
 
284
  if (p != 0)
 
285
    factor_using_division_2kp (t, division_limit / 10, p);
 
286
  else
 
287
    factor_using_division (t, division_limit);
 
288
 
 
289
  if (mpz_cmp_ui (t, 1) != 0)
 
290
    {
 
291
      if (flag_verbose)
 
292
        {
 
293
          printf ("[is number prime?] ");
 
294
          fflush (stdout);
 
295
        }
 
296
      if (mpz_probab_prime_p (t, 3))
 
297
        mpz_out_str (stdout, 10, t);
 
298
      else
 
299
        factor_using_pollard_rho (t, 1, p);
 
300
    }
 
301
}
 
302
 
 
303
main (int argc, char *argv[])
 
304
{
 
305
  mpz_t t;
 
306
  unsigned long p;
 
307
  int i;
 
308
 
 
309
  if (argc > 1 && !strcmp (argv[1], "-v"))
 
310
    {
 
311
      flag_verbose = 1;
 
312
      argv++;
 
313
      argc--;
 
314
    }
 
315
 
 
316
  mpz_init (t);
 
317
  if (argc > 1)
 
318
    {
 
319
      p = 0;
 
320
      for (i = 1; i < argc; i++)
 
321
        {
 
322
          if (!strncmp (argv[i], "-Mp", 3))
 
323
            {
 
324
              p = atoi (argv[i] + 3);
 
325
              mpz_set_ui (t, 1);
 
326
              mpz_mul_2exp (t, t, p);
 
327
              mpz_sub_ui (t, t, 1);
 
328
            }
 
329
          else if (!strncmp (argv[i], "-2kp", 4))
 
330
            {
 
331
              p = atoi (argv[i] + 4);
 
332
              continue;
 
333
            }
 
334
          else
 
335
            {
 
336
              mpz_set_str (t, argv[i], 0);
 
337
            }
 
338
 
 
339
          if (mpz_cmp_ui (t, 0) == 0)
 
340
            puts ("-");
 
341
          else
 
342
            {
 
343
              factor (t, p);
 
344
              puts ("");
 
345
            }
 
346
        }
 
347
    }
 
348
  else
 
349
    {
 
350
      for (;;)
 
351
        {
 
352
          mpz_inp_str (t, stdin, 0);
 
353
          if (feof (stdin))
 
354
            break;
 
355
          mpz_out_str (stdout, 10, t); printf (" = ");
 
356
          factor (t, 0);
 
357
          puts ("");
 
358
        }
 
359
    }
 
360
 
 
361
  exit (0);
 
362
}