ecm.py (1737B)
1 #! /bin/env python 2 3 from sys import argv 4 from random import randint 5 from math import sqrt 6 from dataclasses import dataclass 7 8 @dataclass 9 class Point: 10 x: int = 0 11 y: int = 0 12 is_zero: bool = False 13 14 @dataclass 15 class FactorFound(Exception): 16 factor: int = 0 17 18 # Returns gcd(a, b) and x, y such that ax + by = gcd(a,b) 19 def extended_gcd(a: int, b: int) -> (int, int, int): 20 if b == 0: 21 return a, 1, 0 22 g, x, y = extended_gcd(b, a % b) 23 return g, y, x - y * (a // b) 24 25 def inverse_modulo(a: int, n: int) -> int: 26 g, x, y = extended_gcd(a, n) 27 if g != 1: 28 raise FactorFound(g) 29 return x 30 31 def ec_sum(p: Point, q: Point, A: int, n: int) -> Point: 32 if p.is_zero: 33 return q 34 if q.is_zero: 35 return p 36 if (p.x - q.x) % n == 0 and (p.y + q.y) % n == 0: 37 return Point(is_zero = True) 38 39 if (p.x - q.x) % n != 0: 40 k = ((p.y - q.y) * inverse_modulo(p.x - q.x, n)) % n 41 else: 42 k = ((3 * p.x**2 + A) * inverse_modulo(p.y + q.y, n)) % n 43 44 x = (k**2 - p.x - q.x) % n 45 y = (k * (p.x - x) - p.y) % n 46 47 return Point(x = x, y = y) 48 49 def ec_mul(M: int, p: Point, A: int, n: int) -> Point: 50 if M == 0: 51 return Point(is_zero = True) 52 if M % 2 == 0: 53 return ec_mul(M // 2, ec_sum(p, p, A, n), A, n) 54 return ec_sum(p, ec_mul(M - 1, p, A, n), A, n) 55 56 # Elliptic curve factorization method 57 # If n is prime, this method goes into an infinite loop 58 def find_factor(n: int) -> int: 59 bound = max(int(sqrt(sqrt(sqrt(n)))), 256) 60 61 while True: 62 A = randint(0,n) 63 p = Point(x = randint(0,n), y = randint(0,n)) 64 M = 2 65 while M < bound and not p.is_zero: 66 try: 67 p = ec_mul(M, p, A, n) 68 except FactorFound as ff: 69 #print("Found with A =", A, "M =", M, "and P =", p) 70 return ff.factor 71 M += 1 72 73 N = int(argv[-1]) 74 f = find_factor(N) 75 print(N, "=", f, "*", N//f)