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 }