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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
/*
* Copyright (c) 1995 David I. Bell
* Permission is granted to use, distribute, or modify this source,
* provided that this copyright notice remains intact.
*
* Attempt to factor numbers using elliptic functions.
* y^2 = x^3 + a*x + b (mod N).
*
* Many points (x,y) (mod N) are found that solve the above equation,
* starting from a trivial solution and 'multiplying' that point together
* to generate high powers of the point, looking for such a point whose
* order contains a common factor with N. The order of the group of points
* varies almost randomly within a certain interval for each choice of a
* and b, and thus each choice provides an independent opportunity to
* factor N. To generate a trivial solution, a is chosen and then b is
* selected so that (1,1) is a solution. The multiplication is done using
* the basic fact that the equation is a cubic, and so if a line hits the
* curve in two rational points, then the third intersection point must
* also be rational. Thus by drawing lines between known rational points
* the number of rational solutions can be made very large. When modular
* arithmetic is used, solving for the third point requires the taking of a
* modular inverse (instead of division), and if this fails, then the GCD
* of the failing value and N provides a factor of N. This description is
* only an approximation, read "A Course in Number Theory and Cryptography"
* by Neal Koblitz for a good explanation.
*
* factor(iN, ia, B, force)
* iN is the number to be factored.
* ia is the initial value of a in the equation, and each successive
* value of a is an independent attempt at factoring (default 1).
* B is the limit of the primes that make up the high power that the
* point is raised to for each factoring attempt (default 100).
* force is a flag to attempt to factor numbers even if they are
* thought to already be prime (default FALSE).
*
* Making B larger makes the power the point being raised to contain more
* prime factors, thus increasing the chance that the order of the point
* will be made up of those factors. The higher B is then, the greater
* the chance that any individual attempt will find a factor. However,
* a higher B also slows down the number of independent functions being
* examined. The order of the point for any particular function might
* contain a large prime and so won't succeed even for a really large B,
* whereas the next function might have an order which is quickly found.
* So you want to trade off the depth of a particular search with the
* number of searches made. For example, for factoring 30 digits, I make
* B be about 1000 (probably still too small).
*
* If you have lots of machines available, then you can run parallel
* factoring attempts for the same number by giving different starting
* values of ia for each machine (e.g. 1000, 2000, 3000).
*
* The output as the function is running is (occasionally) the value of a
* when a new function is started, the prime that is being included in the
* high power being calculated, and the current point which is the result
* of the powers so far.
*
* If a factor is found, it is returned and is also saved in the global
* variable f. The number being factored is also saved in the global
* variable N.
*/
obj point {x, y};
global N; /* number to factor */
global a; /* first coefficient */
global b; /* second coefficient */
global f; /* found factor */
define factor(iN, ia, B, force)
{
local C, x, p;
if (!force && ptest(iN, 50))
return 1;
if (isnull(B))
B = 100;
if (isnull(ia))
ia = 1;
obj point x;
a = ia;
b = -ia;
N = iN;
C = isqrt(N);
C = 2 * C + 2 * isqrt(C) + 1;
f = 0;
while (f == 0) {
print "A =", a;
x.x = 1;
x.y = 1;
print 2, x;
x = x ^ (2 ^ (highbit(C) + 1));
for (p = 3; ((p < B) && (f == 0)); p += 2) {
if (!ptest(p, 1))
continue;
print p, x;
x = x ^ (p ^ ((highbit(C) // highbit(p)) + 1));
}
a++;
b--;
}
return f;
}
define point_print(p)
{
print "(" : p.x : "," : p.y : ")" :;
}
define point_mul(p1, p2)
{
local r, m;
if (p2 == 1)
return p1;
if (p1 == p2)
return point_square(&p1);
obj point r;
m = (minv(p2.x - p1.x, N) * (p2.y - p1.y)) % N;
if (m == 0) {
if (f == 0)
f = gcd(p2.x - p1.x, N);
r.x = 1;
r.y = 1;
return r;
}
r.x = (m^2 - p1.x - p2.x) % N;
r.y = ((m * (p1.x - r.x)) - p1.y) % N;
return r;
}
define point_square(p)
{
local r, m;
obj point r;
m = ((3 * p.x^2 + a) * minv(p.y << 1, N)) % N;
if (m == 0) {
if (f == 0)
f = gcd(p.y << 1, N);
r.x = 1;
r.y = 1;
return r;
}
r.x = (m^2 - p.x - p.x) % N;
r.y = ((m * (p.x - r.x)) - p.y) % N;
return r;
}
define point_pow(p, pow)
{
local bit, r, t;
r = 1;
if (isodd(pow))
r = p;
t = p;
for (bit = 2; ((bit <= pow) && (f == 0)); bit <<= 1) {
t = point_square(&t);
if (bit & pow)
r = point_mul(&t, &r);
}
return r;
}
if (config("lib_debug") >= 0) {
print "factor(N, I, B, force) defined";
}
|
