ecm

Elliptic Curve factorization - code and slides
git clone https://git.tronto.net/ecm
Download | Log | Files | Refs | README

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)