Coverage for /usr/lib/python3/dist-packages/sympy/ntheory/ecm.py: 12%
118 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from sympy.ntheory import sieve, isprime
2from sympy.core.numbers import mod_inverse
3from sympy.core.power import integer_log
4from sympy.utilities.misc import as_int
5import random
7rgen = random.Random()
9#----------------------------------------------------------------------------#
10# #
11# Lenstra's Elliptic Curve Factorization #
12# #
13#----------------------------------------------------------------------------#
16class Point:
17 """Montgomery form of Points in an elliptic curve.
18 In this form, the addition and doubling of points
19 does not need any y-coordinate information thus
20 decreasing the number of operations.
21 Using Montgomery form we try to perform point addition
22 and doubling in least amount of multiplications.
24 The elliptic curve used here is of the form
25 (E : b*y**2*z = x**3 + a*x**2*z + x*z**2).
26 The a_24 parameter is equal to (a + 2)/4.
28 References
29 ==========
31 .. [1] https://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf
32 """
34 def __init__(self, x_cord, z_cord, a_24, mod):
35 """
36 Initial parameters for the Point class.
38 Parameters
39 ==========
41 x_cord : X coordinate of the Point
42 z_cord : Z coordinate of the Point
43 a_24 : Parameter of the elliptic curve in Montgomery form
44 mod : modulus
45 """
46 self.x_cord = x_cord
47 self.z_cord = z_cord
48 self.a_24 = a_24
49 self.mod = mod
51 def __eq__(self, other):
52 """Two points are equal if X/Z of both points are equal
53 """
54 if self.a_24 != other.a_24 or self.mod != other.mod:
55 return False
56 return self.x_cord * other.z_cord % self.mod ==\
57 other.x_cord * self.z_cord % self.mod
59 def add(self, Q, diff):
60 """
61 Add two points self and Q where diff = self - Q. Moreover the assumption
62 is self.x_cord*Q.x_cord*(self.x_cord - Q.x_cord) != 0. This algorithm
63 requires 6 multiplications. Here the difference between the points
64 is already known and using this algorithm speeds up the addition
65 by reducing the number of multiplication required. Also in the
66 mont_ladder algorithm is constructed in a way so that the difference
67 between intermediate points is always equal to the initial point.
68 So, we always know what the difference between the point is.
71 Parameters
72 ==========
74 Q : point on the curve in Montgomery form
75 diff : self - Q
77 Examples
78 ========
80 >>> from sympy.ntheory.ecm import Point
81 >>> p1 = Point(11, 16, 7, 29)
82 >>> p2 = Point(13, 10, 7, 29)
83 >>> p3 = p2.add(p1, p1)
84 >>> p3.x_cord
85 23
86 >>> p3.z_cord
87 17
88 """
89 u = (self.x_cord - self.z_cord)*(Q.x_cord + Q.z_cord)
90 v = (self.x_cord + self.z_cord)*(Q.x_cord - Q.z_cord)
91 add, subt = u + v, u - v
92 x_cord = diff.z_cord * add * add % self.mod
93 z_cord = diff.x_cord * subt * subt % self.mod
94 return Point(x_cord, z_cord, self.a_24, self.mod)
96 def double(self):
97 """
98 Doubles a point in an elliptic curve in Montgomery form.
99 This algorithm requires 5 multiplications.
101 Examples
102 ========
104 >>> from sympy.ntheory.ecm import Point
105 >>> p1 = Point(11, 16, 7, 29)
106 >>> p2 = p1.double()
107 >>> p2.x_cord
108 13
109 >>> p2.z_cord
110 10
111 """
112 u = pow(self.x_cord + self.z_cord, 2, self.mod)
113 v = pow(self.x_cord - self.z_cord, 2, self.mod)
114 diff = u - v
115 x_cord = u*v % self.mod
116 z_cord = diff*(v + self.a_24*diff) % self.mod
117 return Point(x_cord, z_cord, self.a_24, self.mod)
119 def mont_ladder(self, k):
120 """
121 Scalar multiplication of a point in Montgomery form
122 using Montgomery Ladder Algorithm.
123 A total of 11 multiplications are required in each step of this
124 algorithm.
126 Parameters
127 ==========
129 k : The positive integer multiplier
131 Examples
132 ========
134 >>> from sympy.ntheory.ecm import Point
135 >>> p1 = Point(11, 16, 7, 29)
136 >>> p3 = p1.mont_ladder(3)
137 >>> p3.x_cord
138 23
139 >>> p3.z_cord
140 17
141 """
142 Q = self
143 R = self.double()
144 for i in bin(k)[3:]:
145 if i == '1':
146 Q = R.add(Q, self)
147 R = R.double()
148 else:
149 R = Q.add(R, self)
150 Q = Q.double()
151 return Q
154def _ecm_one_factor(n, B1=10000, B2=100000, max_curve=200):
155 """Returns one factor of n using
156 Lenstra's 2 Stage Elliptic curve Factorization
157 with Suyama's Parameterization. Here Montgomery
158 arithmetic is used for fast computation of addition
159 and doubling of points in elliptic curve.
161 This ECM method considers elliptic curves in Montgomery
162 form (E : b*y**2*z = x**3 + a*x**2*z + x*z**2) and involves
163 elliptic curve operations (mod N), where the elements in
164 Z are reduced (mod N). Since N is not a prime, E over FF(N)
165 is not really an elliptic curve but we can still do point additions
166 and doubling as if FF(N) was a field.
168 Stage 1 : The basic algorithm involves taking a random point (P) on an
169 elliptic curve in FF(N). The compute k*P using Montgomery ladder algorithm.
170 Let q be an unknown factor of N. Then the order of the curve E, |E(FF(q))|,
171 might be a smooth number that divides k. Then we have k = l * |E(FF(q))|
172 for some l. For any point belonging to the curve E, |E(FF(q))|*P = O,
173 hence k*P = l*|E(FF(q))|*P. Thus kP.z_cord = 0 (mod q), and the unknownn
174 factor of N (q) can be recovered by taking gcd(kP.z_cord, N).
176 Stage 2 : This is a continuation of Stage 1 if k*P != O. The idea utilize
177 the fact that even if kP != 0, the value of k might miss just one large
178 prime divisor of |E(FF(q))|. In this case we only need to compute the
179 scalar multiplication by p to get p*k*P = O. Here a second bound B2
180 restrict the size of possible values of p.
182 Parameters
183 ==========
185 n : Number to be Factored
186 B1 : Stage 1 Bound
187 B2 : Stage 2 Bound
188 max_curve : Maximum number of curves generated
190 References
191 ==========
193 .. [1] Carl Pomerance and Richard Crandall "Prime Numbers:
194 A Computational Perspective" (2nd Ed.), page 344
195 """
196 n = as_int(n)
197 if B1 % 2 != 0 or B2 % 2 != 0:
198 raise ValueError("The Bounds should be an even integer")
199 sieve.extend(B2)
201 if isprime(n):
202 return n
204 from sympy.functions.elementary.miscellaneous import sqrt
205 from sympy.polys.polytools import gcd
206 D = int(sqrt(B2))
207 beta = [0]*(D + 1)
208 S = [0]*(D + 1)
209 k = 1
210 for p in sieve.primerange(1, B1 + 1):
211 k *= pow(p, integer_log(B1, p)[0])
212 for _ in range(max_curve):
213 #Suyama's Parametrization
214 sigma = rgen.randint(6, n - 1)
215 u = (sigma*sigma - 5) % n
216 v = (4*sigma) % n
217 u_3 = pow(u, 3, n)
219 try:
220 # We use the elliptic curve y**2 = x**3 + a*x**2 + x
221 # where a = pow(v - u, 3, n)*(3*u + v)*mod_inverse(4*u_3*v, n) - 2
222 # However, we do not declare a because it is more convenient
223 # to use a24 = (a + 2)*mod_inverse(4, n) in the calculation.
224 a24 = pow(v - u, 3, n)*(3*u + v)*mod_inverse(16*u_3*v, n) % n
225 except ValueError:
226 #If the mod_inverse(16*u_3*v, n) doesn't exist (i.e., g != 1)
227 g = gcd(16*u_3*v, n)
228 #If g = n, try another curve
229 if g == n:
230 continue
231 return g
233 Q = Point(u_3, pow(v, 3, n), a24, n)
234 Q = Q.mont_ladder(k)
235 g = gcd(Q.z_cord, n)
237 #Stage 1 factor
238 if g != 1 and g != n:
239 return g
240 #Stage 1 failure. Q.z = 0, Try another curve
241 elif g == n:
242 continue
244 #Stage 2 - Improved Standard Continuation
245 S[1] = Q.double()
246 S[2] = S[1].double()
247 beta[1] = (S[1].x_cord*S[1].z_cord) % n
248 beta[2] = (S[2].x_cord*S[2].z_cord) % n
250 for d in range(3, D + 1):
251 S[d] = S[d - 1].add(S[1], S[d - 2])
252 beta[d] = (S[d].x_cord*S[d].z_cord) % n
254 g = 1
255 B = B1 - 1
256 T = Q.mont_ladder(B - 2*D)
257 R = Q.mont_ladder(B)
259 for r in range(B, B2, 2*D):
260 alpha = (R.x_cord*R.z_cord) % n
261 for q in sieve.primerange(r + 2, r + 2*D + 1):
262 delta = (q - r) // 2
263 # We want to calculate
264 # f = R.x_cord * S[delta].z_cord - S[delta].x_cord * R.z_cord
265 f = (R.x_cord - S[delta].x_cord)*\
266 (R.z_cord + S[delta].z_cord) - alpha + beta[delta]
267 g = (g*f) % n
268 #Swap
269 T, R = R, R.add(S[D], T)
270 g = gcd(n, g)
272 #Stage 2 Factor found
273 if g != 1 and g != n:
274 return g
276 #ECM failed, Increase the bounds
277 raise ValueError("Increase the bounds")
280def ecm(n, B1=10000, B2=100000, max_curve=200, seed=1234):
281 """Performs factorization using Lenstra's Elliptic curve method.
283 This function repeatedly calls `ecm_one_factor` to compute the factors
284 of n. First all the small factors are taken out using trial division.
285 Then `ecm_one_factor` is used to compute one factor at a time.
287 Parameters
288 ==========
290 n : Number to be Factored
291 B1 : Stage 1 Bound
292 B2 : Stage 2 Bound
293 max_curve : Maximum number of curves generated
294 seed : Initialize pseudorandom generator
296 Examples
297 ========
299 >>> from sympy.ntheory import ecm
300 >>> ecm(25645121643901801)
301 {5394769, 4753701529}
302 >>> ecm(9804659461513846513)
303 {4641991, 2112166839943}
304 """
305 _factors = set()
306 for prime in sieve.primerange(1, 100000):
307 if n % prime == 0:
308 _factors.add(prime)
309 while(n % prime == 0):
310 n //= prime
311 rgen.seed(seed)
312 while(n > 1):
313 try:
314 factor = _ecm_one_factor(n, B1, B2, max_curve)
315 except ValueError:
316 raise ValueError("Increase the bounds")
317 _factors.add(factor)
318 n //= factor
320 factors = set()
321 for factor in _factors:
322 if isprime(factor):
323 factors.add(factor)
324 continue
325 factors |= ecm(factor)
326 return factors