1from sympy.ntheory import sieve, isprime 2from sympy.core.power import integer_log 3from sympy.core.compatibility import as_int 4import random 5 6rgen = random.Random() 7 8#----------------------------------------------------------------------------# 9# # 10# Lenstra's Elliptic Curve Factorization # 11# # 12#----------------------------------------------------------------------------# 13 14 15class Point: 16 """Montgomery form of Points in an elliptic curve. 17 In this form, the addition and doubling of points 18 does not need any y-coordinate information thus 19 decreasing the number of operations. 20 Using Montgomery form we try to perform point addition 21 and doubling in least amount of multiplications. 22 23 The elliptic curve used here is of the form 24 (E : b*y**2*z = x**3 + a*x**2*z + x*z**2). 25 The a_24 parameter is equal to (a + 2)/4. 26 27 References 28 ========== 29 30 .. [1] http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf 31 """ 32 33 def __init__(self, x_cord, z_cord, a_24, mod): 34 """ 35 Initial parameters for the Point class. 36 37 Parameters 38 ========== 39 40 x_cord : X coordinate of the Point 41 z_cord : Z coordinate of the Point 42 a_24 : Parameter of the elliptic curve in Montgomery form 43 mod : modulus 44 """ 45 self.x_cord = x_cord 46 self.z_cord = z_cord 47 self.a_24 = a_24 48 self.mod = mod 49 50 def __eq__(self, other): 51 """Two points are equal if X/Z of both points are equal 52 """ 53 from sympy import mod_inverse 54 if self.a_24 != other.a_24 or self.mod != other.mod: 55 return False 56 return self.x_cord * mod_inverse(self.z_cord, self.mod) % self.mod ==\ 57 other.x_cord * mod_inverse(other.z_cord, self.mod) % self.mod 58 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 algorihtm 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. 69 70 71 Parameters 72 ========== 73 74 Q : point on the curve in Montgomery form 75 diff : self - Q 76 77 Examples 78 ======== 79 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) 95 96 def double(self): 97 """ 98 Doubles a point in an elliptic curve in Montgomery form. 99 This algorithm requires 5 multiplications. 100 101 Examples 102 ======== 103 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, v = self.x_cord + self.z_cord, self.x_cord - self.z_cord 113 u, v = u*u, v*v 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) 118 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. 125 126 Parameters 127 ========== 128 129 k : The positive integer multiplier 130 131 Examples 132 ======== 133 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 152 153 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. 160 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. 167 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). 175 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. 181 182 Parameters 183 ========== 184 185 n : Number to be Factored 186 B1 : Stage 1 Bound 187 B2 : Stage 2 Bound 188 max_curve : Maximum number of curves generated 189 190 References 191 ========== 192 193 .. [1] Carl Pomerance and Richard Crandall "Prime Numbers: 194 A Computational Perspective" (2nd Ed.), page 344 195 """ 196 from sympy import gcd, mod_inverse, sqrt 197 n = as_int(n) 198 if B1 % 2 != 0 or B2 % 2 != 0: 199 raise ValueError("The Bounds should be an even integer") 200 sieve.extend(B2) 201 202 if isprime(n): 203 return n 204 205 curve = 0 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 while(curve <= max_curve): 213 curve += 1 214 215 #Suyama's Paramatrization 216 sigma = rgen.randint(6, n - 1) 217 u = (sigma*sigma - 5) % n 218 v = (4*sigma) % n 219 diff = v - u 220 u_3 = pow(u, 3, n) 221 222 try: 223 C = (pow(diff, 3, n)*(3*u + v)*mod_inverse(4*u_3*v, n) - 2) % n 224 except ValueError: 225 #If the mod_inverse(4*u_3*v, n) doesn't exist 226 return gcd(4*u_3*v, n) 227 228 a24 = (C + 2)*mod_inverse(4, n) % n 229 Q = Point(u_3 , pow(v, 3, n), a24, n) 230 Q = Q.mont_ladder(k) 231 g = gcd(Q.z_cord, n) 232 233 #Stage 1 factor 234 if g != 1 and g != n: 235 return g 236 #Stage 1 failure. Q.z = 0, Try another curve 237 elif g == n: 238 continue 239 240 #Stage 2 - Improved Standard Continuation 241 S[1] = Q.double() 242 S[2] = S[1].double() 243 beta[1] = (S[1].x_cord*S[1].z_cord) % n 244 beta[2] = (S[2].x_cord*S[2].z_cord) % n 245 246 for d in range(3, D + 1): 247 S[d] = S[d - 1].add(S[1], S[d - 2]) 248 beta[d] = (S[d].x_cord*S[d].z_cord) % n 249 250 g = 1 251 B = B1 - 1 252 T = Q.mont_ladder(B - 2*D) 253 R = Q.mont_ladder(B) 254 255 for r in range(B, B2, 2*D): 256 alpha = (R.x_cord*R.z_cord) % n 257 for q in sieve.primerange(r + 2, r + 2*D + 1): 258 delta = (q - r) // 2 259 f = (R.x_cord - S[d].x_cord)*(R.z_cord + S[d].z_cord) -\ 260 alpha + beta[delta] 261 g = (g*f) % n 262 #Swap 263 T, R = R, R.add(S[D], T) 264 g = gcd(n, g) 265 266 #Stage 2 Factor found 267 if g != 1 and g != n: 268 return g 269 270 #ECM failed, Increase the bounds 271 raise ValueError("Increase the bounds") 272 273 274def ecm(n, B1=10000, B2=100000, max_curve=200, seed=1234): 275 """Performs factorization using Lenstra's Elliptic curve method. 276 277 This function repeatedly calls `ecm_one_factor` to compute the factors 278 of n. First all the small factors are taken out using trial division. 279 Then `ecm_one_factor` is used to compute one factor at a time. 280 281 Parameters 282 ========== 283 284 n : Number to be Factored 285 B1 : Stage 1 Bound 286 B2 : Stage 2 Bound 287 max_curve : Maximum number of curves generated 288 seed : Initialize pseudorandom generator 289 290 Examples 291 ======== 292 293 >>> from sympy.ntheory import ecm 294 >>> ecm(25645121643901801) 295 {5394769, 4753701529} 296 >>> ecm(9804659461513846513) 297 {4641991, 2112166839943} 298 """ 299 _factors = set() 300 for prime in sieve.primerange(1, 100000): 301 if n % prime == 0: 302 _factors.add(prime) 303 while(n % prime == 0): 304 n //= prime 305 rgen.seed(seed) 306 while(n > 1): 307 try: 308 factor = _ecm_one_factor(n, B1, B2, max_curve) 309 except ValueError: 310 raise ValueError("Increase the bounds") 311 _factors.add(factor) 312 n //= factor 313 314 factors = set() 315 for factor in _factors: 316 if isprime(factor): 317 factors.add(factor) 318 continue 319 factors |= ecm(factor) 320 return factors 321