"Fossies" - the Fresh Open Source Software Archive

Member "NZMATH-1.2.0/nzmath/factor/ecm.py" (19 Nov 2012, 9237 Bytes) of package /linux/misc/old/NZMATH-1.2.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Python source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file. For more information about "ecm.py" see the Fossies "Dox" file reference documentation.

```    1 """
2 Elliptic Curve Method Factorization
3
4 It's using Montgomery's parametrization.
5 """
6
7 from __future__ import division
8 import logging
9 import random
10 import nzmath.arith1 as _arith1
11 import nzmath.gcd as _gcd
12 import nzmath.prime as _prime
13
14 _log = logging.getLogger('nzmath.factor.ecm')
15
16
17 # curve types
18 S = SUYAMA = 1
19 B = BERNSTEIN = 2
20 A1 = ASUNCION1 = 5
21 A2 = ASUNCION2 = 6
22 A3 = ASUNCION3 = 8
23 A4 = ASUNCION4 = 9
24 A5 = ASUNCION5 = 10
25 N3 = NAKAMURA3 = 13
26
27 class Curve (object):
28     """
29     Elliptic curves for factorization.
30     """
31
32     _CURVE_TYPES = (S, B, A1, A2, A3, A4, A5, N3)
33
34     def __init__(self, c):
35         """
36         Initialize a Curve object with Montgomery's parameter c.
37         """
38         self.c = c
39         self.c2 = (c + 2)//4
40
41     @classmethod
42     def get_random_curve_with_point(cls, curve_type, n, bounds):
43         """
44         Return the curve with parameter C and a point Q on the curve,
45         according to the curve_type, factorization target n and the
46         bounds for stages.
47
48         curve_type should be one of the module constants corresponding
49         to parameters:
50           S: Suyama's parameter selection strategy
51           B: Bernstein's [2:1], [16,18,4,2]
52           A1: Asuncion's [2:1], [4,14,1,1]
53           A2: Asuncion's [2:1], [16,174,4,41]
54           A3: Asuncion's [3:1], [9,48,1,2]
55           A4: Asuncion's [3:1], [9,39,1,1]
56           A5: Asuncion's [4:1], [16,84,1,1]
57           N3: Nakamura's [2:1], [28,22,7,3]
58           This is a class method.
59         """
60         bound = bounds.first
61         if curve_type not in cls._CURVE_TYPES:
62             raise ValueError("Input curve_type is wrong.")
63         if curve_type == SUYAMA:
64             t = n
65             while _gcd.gcd(t, n) != 1:
66                 sigma = random.randrange(6, bound + 1)
67                 u, v = (sigma**2 - 5) % n, (4*sigma) % n
68                 t = 4*(u**3)*v
69             d = _arith1.inverse(t, n)
70             curve = cls((((u - v)**3 * (3*u + v)) * d - 2) % n)
71             start_point = Point(pow(u, 3, n), pow(v, 3, n))
72         elif curve_type == BERNSTEIN:
73             d = random.randrange(1, bound + 1)
74             start_point = Point(2, 1)
75             curve = cls((4*d + 2) % n)
76         elif curve_type == ASUNCION1:
77             d = random.randrange(1, bound + 1)
78             start_point = Point(2, 1)
79             curve = cls((d + 1) % n)
80         elif curve_type == ASUNCION2:
81             d = random.randrange(1, bound + 1)
82             start_point = Point(2, 1)
83             curve = cls((4*d + 41) % n)
84         elif curve_type == ASUNCION3:
85             d = random.randrange(1, bound + 1)
86             start_point = Point(3, 1)
87             curve = cls((d + 2) % n)
88         elif curve_type == ASUNCION4:
89             d = random.randrange(1, bound + 1)
90             start_point = Point(3, 1)
91             curve = cls((d + 1) % n)
92         elif curve_type == ASUNCION5:
93             d = random.randrange(1, bound + 1)
94             start_point = Point(4, 1)
95             curve = cls((d + 1) % n)
96         elif curve_type == NAKAMURA3:
97             d = random.randrange(1, bound + 1)
98             start_point = Point(2, 1)
99             curve = cls((7*d + 3) % n)
100         return curve, start_point
101
102
103 class Point (tuple):
104     """
105     x and z projective coordinates of elliptic curves.
106     """
107     def __new__(cls, x, z):
108         """
109         Create a new Point object.
110         """
111         return super(Point, cls).__new__(Point, (x, z))
112
113     def _get_x(self):
114         "getter for x"
115         return self[0]
116
117     def _get_z(self):
118         "getter for z"
119         return self[1]
120
121     x = property(_get_x, None, None, "x")
122     z = property(_get_z, None, None, "z")
123
124 POINT_AT_INFINITY = Point(0, 0)
125
126
127 class Bounds (object):
128     """
129     Bounds for ECM.
130
131     public attributes:
132       first: bound for first stage
133       second: bound for second stage
134     """
135     def __init__(self, n):
136         """
137         Create a Bounds object from target number n.
138         """
139         if _arith1.log(n, 10) < 20:
140             self.first = 1000
141         elif _arith1.log(n, 10) < 25:
142             self.first = 10000
143         else:
144             self.first = 100000
145         self.second = self.first * 50
146
147     def increment(self, scale=10):
148         """
149         Multiply both bounds by optional argument scale (default to 10).
150         """
151         self.first *= scale
152         self.second *= scale
153
154     def __getitem__(self, key):
155         """
156         Return first by key 1, and second key 2.
157         Otherwise KeyError is raised.
158         """
159         if key == 1:
160             return self.first
161         if key == 2:
162             return self.second
163         raise KeyError("Key must be 1 or 2.")
164
165     def __hash__(self):
166         """
167         Bounds objects are mutable and thus unhashable.
168         """
169         raise TypeError("Bounds objects are unhashable")
170
171     def __str__(self):
172         """
173         (first, second)
174         """
175         return "(%d, %d)" % (self.first, self.second)
176
177
178 def ecm(n, curve_type=A1, incs=3, trials=20, **options):
179     """
180     Find a factor of n with Elliptic Curve Method.
181     An unsuccessful factorization returns 1.
182
183     There are a few optional arguments.
184
185     By 'curve_type', the function choose a family of curves.
186     Please use a module constant to specify the curve_type.
187
188     The second optional argument 'incs' specifies the number of times
189     for changing bounds.  The function repeats factorization trials
190     several times changing curves with a fixed bounds.
191
192     The last optional argument 'trials' can control how quickly move
193     on to the next higher bounds.
194     """
195     # verbosity
196     verbose = options.get('verbose', False)
197     if verbose:
198         _log.setLevel(logging.DEBUG)
199         _log.debug("verbose")
200     else:
201         _log.setLevel(logging.NOTSET)
202
203     # trivial checks
204     if _prime.primeq(n):
205         _log.info("%d is prime!" % n)
206         return n
207     if _gcd.gcd(n, 6) != 1:
208         _log.info("%d is not coprime to 6!" % n)
209         if n % 2 == 0:
210             return 2
211         if n % 3 == 0:
212             return 3
213
214     # main loop
215     bounds = Bounds(n)
216     for inc in range(incs):
217         _log.info("bounds increment %d times" % inc)
218         _log.debug("Bounds B1, B2 = %s" % (bounds,))
219         for trial in range(trials):
220             _log.info("Trial #: %d" % (trial + 1))
221             curve, point = Curve.get_random_curve_with_point(curve_type, n, bounds)
222             _log.debug("Curve param: %d" % curve.c)
223             g = stage1(n, bounds, curve, point)
224             if 1 < g < n:
225                 return g
226             g = stage2(n, bounds, curve, point)
227             if 1 < g < n:
228                 return g
229         bounds.increment()
230     _log.info("ECM 2 step test failed!")
231
232     if not verbose:
233         _log.setLevel(logging.DEBUG)
234     return 1
235
236 def stage1(n, bounds, C, Q):
237     """
238     ECM stage 1 for factoring n.
239     The upper bound for primes to be tested is bounds.first.
240     It uses curve C and starting point Q.
241     """
242     for p in _prime.generator_eratosthenes(bounds.first):
243         q = p
244         while q < bounds.first:
245             Q = mul(Q, p, C, n)
246             q *= p
247     g = _gcd.gcd(Q.z, n)
248     _log.debug("Stage 1: %d" % g)
249     return g
250
251 def stage2(n, bounds, C, Q):
252     """
253     ECM stage 2 for factoring n.
254     The upper bounds for primes to be tested are stored in bounds.
255     It uses curve C and starting point Q.
256     """
257     d1 = bounds.first + random.randrange(1, 16)
258     d2 = 0
259     while _gcd.gcd(d1, d2) != 1:
260         d2 = random.randrange(2, d1//5 + 1) # We want to keep d2 small
261     for i in range(1, bounds.second//d1):
262         if _gcd.gcd(i, d2) != 1:
263             continue
264         for j in range(1, d1//2):
265             if _gcd.gcd(j, d1) == 1:
266                 Q = mul(Q, i*d1 + j*d2, C, n)
267                 if i*d1 - j*d2 > bounds.first:
268                     Q = mul(Q, i*d1 - j*d2, C, n)
269         g = _gcd.gcd(Q.z, n)
270         if 1 < g < n:
271             _log.debug("Stage 2: %d" % g)
272             return g
273     _log.debug("Stage 2: %d" % g)
274     return 1
275
276 def mul(Q, x, C, n):
277     """
278     Return x*Q on the curve C mod n.
279
280     m*Q and (m+1)*Q are being tracked in the main loop.
281     """
282     if x == 0:
283         return POINT_AT_INFINITY
284     if x == 1:
285         return Q
286     if x == 2:
287         return double(Q, C, n)
288     minor = Q
289     major = double(Q, C, n)
290     binary = _arith1.expand(x, 2)
291     lastbit, binary = binary[0], binary[1:]
292     while binary:
293         if binary.pop() == 1:
294             minor = add(major, minor, Q, n)
295             major = double(major, C, n)
296         else:
297             major = add(minor, major, Q, n)
298             minor = double(minor, C, n)
299     if lastbit:
300         return add(minor, major, Q, n)
301     return double(minor, C, n)
302
303 # Montgomery's Arithmetic
304
305 def double(Q, C, n):
306     """
307     Return the doubled Q on the curve C mod n.
308     """
309     u = (Q.x + Q.z)**2
310     v = (Q.x - Q.z)**2
311     t = C.c2 * (u - v) + v
312     return Point(u * v % n, (u - v) * t % n)
313
314 def add(P, Q, M, n):
315     """
316     Return the sum of P and Q mod n with a help of auxiliary point M.
317     M must be P - Q. The curve is not explicitly necessary.
318     """
319     u = (P.x + P.z) * (Q.x - Q.z)
320     v = (P.x - P.z) * (Q.x + Q.z)
321     w = (u + v)**2
322     t = (u - v)**2
323     return Point(M.z * w % n, M.x * t % n)
```