ecm

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

ecm.cpp (2283B)


      1 #include "bigint.h"
      2 #include "zmodn.h"
      3 
      4 #include <cstdint>
      5 #include <iostream>
      6 #include <variant>
      7 
      8 constexpr BigInt N(NUMBER);
      9 
     10 class Point {
     11 public:
     12 	bool is_infinity;
     13 	Zmod<N> x;
     14 	Zmod<N> y;
     15 
     16 	Point(Zmod<N> a, Zmod<N> b) : is_infinity{false}, x{a}, y{b} {}
     17 	Point(bool inf) : is_infinity{inf}, x{0}, y{0} {}
     18 };
     19 
     20 std::variant<Point, BigInt<>>
     21 sum_or_factor(BigInt<> a, Point p, Point q) {
     22 	if (p.is_infinity) return q;
     23 	if (q.is_infinity) return p;
     24 	if (p.x == q.x && p.y + q.y == Zmod<N>(0)) return Point(true);
     25 
     26 	Zmod<N> l(0);
     27 	if (p.x != q.x)
     28 		if (auto inv_x = (p.x-q.x).inverse(); inv_x.has_value())
     29 			l = (p.y-q.y) * inv_x.value();
     30 		else
     31 			return std::get<0>(extended_gcd(N, (p.x-q.x).toint()));
     32 	else
     33 		if (auto inv_y = (p.y+q.y).inverse(); inv_y.has_value())
     34 			l = (Zmod<N>(3) * p.x * p.x + a) * inv_y.value();
     35 		else
     36 			return std::get<0>(extended_gcd(N, (p.y+q.y).toint()));
     37 
     38 	auto x = l*l - p.x - q.x;
     39 	auto y = l*(p.x-x) - p.y;
     40 
     41 	return Point(x, y);
     42 }
     43 
     44 std::variant<Point, BigInt<>>
     45 product_or_factor(BigInt<> a, std::variant<Point, BigInt<>> pi, BigInt<> m) {
     46 	if (std::holds_alternative<BigInt<>>(pi))
     47 		return std::get<BigInt<>>(pi);
     48 
     49 	auto p = std::get<Point>(pi);
     50 
     51 	if (m == 0)
     52 		return Point(true); // Anything multiplied by 0 is 0
     53 
     54 	// Divide-and-conquer multiplication (power) algorithm
     55 	if (m % 2 == 0) {
     56 		return product_or_factor(a, sum_or_factor(a, p, p), m/2);
     57 	} else {
     58 		auto pp = product_or_factor(a, p, m-1);
     59 		if (std::holds_alternative<BigInt<>>(pp))
     60 			return std::get<BigInt<>>(pp);
     61 		return sum_or_factor(a, p, std::get<Point>(pp));
     62 	}
     63 }
     64 
     65 BigInt<> find_factor_ecm() {
     66 	// Find a factor of the integer N.
     67 	// If N is prime, this method goes into an infinite loop.
     68 
     69 	while (true) {
     70 		BigInt a = BigInt<>::random(N);
     71 		Point p(BigInt<>::random(N), BigInt<>::random(N));
     72 		for (BigInt m = 2;
     73 		    (m < 256 || m*m*m*m*m*m*m*m < N) && !p.is_infinity; m += 1) {
     74 			auto x = product_or_factor(a, p, m);
     75 			if (std::holds_alternative<BigInt<>>(x))
     76 				return std::get<BigInt<>>(x);
     77 			else
     78 				p = std::get<Point>(x);
     79 		}
     80 	}
     81 }
     82 
     83 int main() {
     84 	// N is a compile-time constant
     85 	if (BigInt f = find_factor_ecm(); f > 1 && f < N)
     86 		std::cout << N << " = " << f << " * " << N/f << std::endl;
     87 	else
     88 		std::cout << N << " is prime" << std::endl;
     89 
     90 	return 0;
     91 }