1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
|
/*
* Copyright 1996 Thierry Bousch
* Licensed under the Gnu Public License, Version 2
*
* $Id: pollard.c,v 1.8 1996/07/18 17:52:27 bousch Exp $
*
* Factorization using Pollard's rho method. This method is called
* "Monte-Carlo factorization" by Knuth and appears as Algorithm B in
* "Seminumerical Algorithms". The running time is roughly p^(1/2)
* where p is the smallest prime factor.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "saml.h"
#include "factorint.h"
int find_factor_pollard (mref_t p, mref_t N, int max_loops)
{
mref_t fast, slow, one;
char buffer[12];
int loops, seed=2, l, k, found=0;
double log_N;
if (max_loops < -1) {
log_N = log10(atof(mref_string(N)));
if (log_N < 10.0)
max_loops = -1; /* don't bother with MPQS */
else if (log_N > 70.0)
max_loops = 1000000;
else
max_loops = pow(10.0, (50.0+log_N)/20.0);
}
fast = mref_new();
slow = mref_new();
one = mref_new();
mref_build(one, ST_INTEGER, "1");
new_seed:
sprintf(buffer, "%d", seed);
mref_build(slow, ST_INTEGER, buffer);
mref_copy(fast, slow);
k = l = 1; loops = 0;
while(1) {
if (max_loops >= 0 && --max_loops < 0) {
/* Too many loops, give up Pollard's method */
break;
}
if (0==(++loops & 1023) && !quiet) {
/* Display a cute rotating cursor */
buffer[0] = "-/|\\"[(loops>>10)&3];
buffer[1] = '\b';
fwrite(buffer, 1, 2, stderr);
}
/* Now apply z->z^2+1 to fast */
mref_mul(fast, fast, fast);
mref_add(fast, fast, one);
mref_mod(fast, fast, N);
/* Factor found? */
mref_sub(p, fast, slow);
if (!mref_notzero(p)) {
/*
* End of cycle: the algorithm failed. For instance
* when seed=2, N=10592225491. So we try again
* with another seed.
*/
++seed;
goto new_seed;
}
mref_gcd(p, p, N);
if (mref_differ(p, one)) {
/* We've found a non-trivial factor */
if (is_pseudo_prime(p)) {
found = 1;
break;
}
}
if (--k == 0) {
k = l = 2*l;
mref_copy(slow, fast);
}
}
mref_free(fast);
mref_free(slow);
mref_free(one);
return found;
}
|
