From 7f1b8a358515b45d717c3c4b5baf2fbc0f568170 Mon Sep 17 00:00:00 2001 From: Sebastiano Tronto Date: Wed, 26 Feb 2025 17:12:06 +0100 Subject: Initial commit --- README.md | 63 +++++ code/cpp/a.out | Bin 0 -> 43664 bytes code/cpp/bigint.h | 266 ++++++++++++++++++++ code/cpp/ecm.cpp | 91 +++++++ code/cpp/naive.cpp | 23 ++ code/cpp/zmodn.h | 79 ++++++ code/python/ecm.py | 75 ++++++ code/python/naive.py | 18 ++ images/beer.jpg | Bin 0 -> 508889 bytes images/clock.png | Bin 0 -> 14805 bytes images/clock2.png | Bin 0 -> 30188 bytes images/demo.jpg | Bin 0 -> 24404 bytes images/ec1.png | Bin 0 -> 17744 bytes images/ec2.png | Bin 0 -> 18917 bytes images/ec3.png | Bin 0 -> 14559 bytes images/euclid.png | Bin 0 -> 1054492 bytes images/factorization.webp | Bin 0 -> 14496 bytes images/number-line.svg | 262 +++++++++++++++++++ images/numbers.jpg | Bin 0 -> 164839 bytes images/questions.png | Bin 0 -> 125274 bytes images/sum-1a.png | Bin 0 -> 18339 bytes images/sum-1b.png | Bin 0 -> 19456 bytes images/sum-1c.png | Bin 0 -> 20134 bytes images/sum-2a.png | Bin 0 -> 18120 bytes images/sum-2b.png | Bin 0 -> 24334 bytes images/sum-2c.png | Bin 0 -> 25172 bytes images/sum-3a.png | Bin 0 -> 17859 bytes images/sum-3b.png | Bin 0 -> 24262 bytes images/sum-3c.png | Bin 0 -> 25047 bytes images/sum-4a.png | Bin 0 -> 17644 bytes images/sum-4b.png | Bin 0 -> 21832 bytes images/sum-4c.png | Bin 0 -> 22574 bytes index.html | 622 ++++++++++++++++++++++++++++++++++++++++++++++ 33 files changed, 1499 insertions(+) create mode 100644 README.md create mode 100644 code/cpp/a.out create mode 100644 code/cpp/bigint.h create mode 100644 code/cpp/ecm.cpp create mode 100644 code/cpp/naive.cpp create mode 100644 code/cpp/zmodn.h create mode 100755 code/python/ecm.py create mode 100755 code/python/naive.py create mode 100644 images/beer.jpg create mode 100644 images/clock.png create mode 100644 images/clock2.png create mode 100644 images/demo.jpg create mode 100644 images/ec1.png create mode 100644 images/ec2.png create mode 100644 images/ec3.png create mode 100644 images/euclid.png create mode 100644 images/factorization.webp create mode 100644 images/number-line.svg create mode 100644 images/numbers.jpg create mode 100644 images/questions.png create mode 100644 images/sum-1a.png create mode 100644 images/sum-1b.png create mode 100644 images/sum-1c.png create mode 100644 images/sum-2a.png create mode 100644 images/sum-2b.png create mode 100644 images/sum-2c.png create mode 100644 images/sum-3a.png create mode 100644 images/sum-3b.png create mode 100644 images/sum-3c.png create mode 100644 images/sum-4a.png create mode 100644 images/sum-4b.png create mode 100644 images/sum-4c.png create mode 100644 index.html diff --git a/README.md b/README.md new file mode 100644 index 0000000..75f495a --- /dev/null +++ b/README.md @@ -0,0 +1,63 @@ +# Elliptic curve factorization method + +Slides and code for the a presentation about +[Lenstra's elliptic-curve factorization](https://en.wikipedia.org/wiki/Lenstra_elliptic-curve_factorization). + +See also [this blog post](https://sebastiano.tronto.net/blog/2025-02-27-ecm). + +## Abstract + +Elliptic curves are mathematical objects that have both a geometric and an +arithmetic side. They turn out to be useful for real-world applications +because they sit in a sweet spot: they are complicated enough to have +interesting and useful arithmetic properties, but simple enough to be +implemented in software in an efficient way. For example, they are used in +cryptographic schemes, such as the Elliptic-curve Diffie-Hellman scheme, +to obtain greater security with smaller keys. + +After introducing elliptic curves and modular arithmetic, we will take a +look at the elliptic curve factorization method (ECM), one of the most +efficient method to find the prime factors of an integer number. We +will see in practice how much faster this method is compared to a naive +algorithm, and we'll see that the implementation of this method is not +that hard at all. + +## Slides + +The slides are a single html file, `index.html`. They rely on a couple +of external JavaScript libraries. They are also hosted +[here](https://sebastiano.tronto.net/talks/ecm). + +## Code + +The folder `code/python` contains two files: + +* `ecm.py`: an implementation of the ECM algorithm. +* `naive.py`: an implementation of the simple O(√n) algorithm for finding + a factor of a number, for comparison. + +To use any of the two, pass the number to factor as a command-line argument, +for example: + +``` +$ ./ecm.py 255000007030000033 +255000007030000033 = 510000011 * 500000003 +``` + +Some benchmarks (note: the ECM is randomized, the time can vary a lot): + +``` +$ time ./ecm.py 255000007030000033 +255000007030000033 = 510000011 * 500000003 + 0m01.45s real 0m01.43s user 0m00.01s system +$ time ./naive.py 255000007030000033 +255000007030000033 = 500000003 * 510000011 + 0m26.62s real 0m26.50s user 0m00.01s system +``` + +### C++ code (experimental) + +The folder `code/cpp` contains an experimental implementation of the ECM +algorithm in C++. It works, but it is very slow: it requires suppport +for compile-time big integers, which I implement in an inefficient way. I +may optimize this code in the future. diff --git a/code/cpp/a.out b/code/cpp/a.out new file mode 100644 index 0000000..b52f011 Binary files /dev/null and b/code/cpp/a.out differ diff --git a/code/cpp/bigint.h b/code/cpp/bigint.h new file mode 100644 index 0000000..096ed9c --- /dev/null +++ b/code/cpp/bigint.h @@ -0,0 +1,266 @@ +#ifndef BIGUNSIGNED_H +#define BIGUNSIGNED_H + +#include +#include +#include +#include + +constexpr uint64_t abs64(int64_t); +constexpr uint64_t pow10(uint64_t); + +// Big integer class for numbers of at most N decimal digits. +// The number E is used to tune the size of each digit, mostly for +// testing purposes. + +template +requires (E < 10) +class BigInt { +public: + // The member variables sign and digits are declared public so that + // BigInt becomes a structural type and can be used in templates. + + static constexpr uint64_t M = pow10(E); + static constexpr uint64_t D = (N / E) + 1; + + bool sign; + uint64_t digits[D]; + + constexpr BigInt() : sign{true} { + std::fill(digits, digits+D, 0); + } + + constexpr BigInt(int64_t n) : sign{n >= 0} { + std::fill(digits, digits+D, 0); + digits[0] = abs64(n); + carryover(); + } + + constexpr BigInt(const std::string_view s) : sign{true} { + std::fill(digits, digits+D, 0); + if (s.size() == 0) + return; + for (int i = s.size()-1, j = 0; i >= 0; i--, j++) { + if (s[i] == '\'') + continue; + if (i == 0 && s[i] == '-') { + sign = false; + break; + } + digits[j/E] += (pow10(j % E)) + * static_cast(s[i] - '0'); + } + } + + constexpr auto operator<=>(const BigInt& other) const { + if (sign != other.sign) + return sign <=> other.sign; + + for (int i = D-1; i >= 0; i--) + if (digits[i] != other.digits[i]) + return sign ? + digits[i] <=> other.digits[i] : + other.digits[i] <=> digits[i]; + + return 0 <=> 0; + } + + constexpr bool operator==(const BigInt& other) const = default; + + constexpr BigInt abs() const { + BigInt ret = *this; + ret.sign = true; + return ret; + } + + constexpr BigInt operator-() const { + if (*this == 0) + return 0; + BigInt ret = *this; + ret.sign = !ret.sign; + return ret; + } + + constexpr BigInt operator+(const BigInt& z) const { + if (sign && z.sign) + return positive_sum(*this, z); + else if (sign && !z.sign) + return positive_diff(*this, -z); + else if (!sign && z.sign) + return positive_diff(z, -*this); + else + return -positive_sum(-*this, -z); + } + + constexpr BigInt operator-(const BigInt& z) const { + return *this + (-z); + } + + constexpr BigInt operator*(const BigInt& z) const { + BigInt ret; + ret.sign = !(sign ^ z.sign); + for (int i = 0; i < D; i++) + for (int j = 0; i+j < D; j++) + ret.digits[i+j] += digits[i] * z.digits[j]; + ret.carryover(); + return ret; + } + + constexpr BigInt operator/(const BigInt& z) const { + auto [q, r] = euclidean_division(*this, z); + return q; + } + + constexpr BigInt operator%(const BigInt& z) const { + auto [q, r] = euclidean_division(*this, z); + return r; + } + + constexpr BigInt operator+=(const BigInt& z) { return *this = *this + z; } + constexpr BigInt operator++() { return *this += 1; } + constexpr BigInt operator-=(const BigInt& z) { return *this = *this - z; } + constexpr BigInt operator--() { return *this -= 1; } + constexpr BigInt operator*=(const BigInt& z) { return *this = *this * z; } + constexpr BigInt operator/=(const BigInt& z) { return *this = *this / z; } + constexpr BigInt operator%=(const BigInt& z) { return *this = *this % z; } + + static BigInt random(BigInt r) { + std::random_device rd; + std::default_random_engine rng(rd()); + std::uniform_int_distribution distribution(0, M-1); + + BigInt ret; + for (uint64_t i = 0; i < D; i++) + ret.digits[i] = distribution(rng); + + return ret % r; + } + + friend std::ostream& operator<<(std::ostream& os, const BigInt& z) { + if (z == 0) { + os << "0"; + return os; + } + + if (!z.sign) + os << "-"; + + int j; + for (j = z.D-1; z.digits[j] == 0; j--) ; + os << z.digits[j]; // Top digit is not padded + + for (int i = j-1; i >= 0; i--) { + std::string num = std::to_string(z.digits[i]); + os << std::string(E - num.length(), '0') << num; + } + return os; + } + +private: + constexpr void carryover() { + for (int i = 1; i < D; i++) { + auto c = digits[i-1] / M; + digits[i-1] -= c * M; + digits[i] += c; + } + } + + constexpr BigInt half() const { + BigInt ret; + uint64_t carry = 0; + for (int i = D-1; i >= 0; i--) { + ret.digits[i] += (digits[i] + M * carry) / 2; + carry = digits[i] % 2; + } + return ret; + } + + static constexpr BigInt powM(uint64_t e) { + BigInt ret; + ret.digits[e] = 1; + return ret; + } + + // Sum of non-negative integers + static constexpr BigInt positive_sum(const BigInt& x, const BigInt& y) { + BigInt ret; + for (int i = 0; i < D; i++) + ret.digits[i] = x.digits[i] + y.digits[i]; + ret.carryover(); + return ret; + } + + // Difference of non-negative integers (result may be negative) + static constexpr BigInt positive_diff(const BigInt& x, const BigInt& y) { + if (y > x) + return -positive_diff(y, x); + + BigInt ret; + uint64_t carry = 0; + for (int i = 0; i < D; i++) { + uint64_t oldcarry = carry; + if (x.digits[i] < y.digits[i] + oldcarry) { + ret.digits[i] = M; + carry = 1; + } else { + carry = 0; + } + ret.digits[i] += x.digits[i]; + ret.digits[i] -= y.digits[i] + oldcarry; + } + ret.carryover(); + return ret; + } + + // Division with remainder, UB if y == 0 + static constexpr std::pair + euclidean_division(const BigInt& x, const BigInt& y) { + auto [q, r] = positive_div(x.abs(), y.abs()); + if (x.sign && y.sign) + return std::pair(q, r); + else if (x.sign && !y.sign) + return r == 0 ? std::pair(-q, 0) : std::pair(-q-1, y+r); + else if (!x.sign && y.sign) + return r == 0 ? std::pair(-q, r) : std::pair(-q-1, y-r); + else + return std::pair(q, -r); + } + + // Division with remainder of non-negative integers, UB if y == 0 + // This method is inefficient (O(log(x/y)) BigInt multiplications) + static constexpr std::pair + positive_div(const BigInt& x, const BigInt& y) { + BigInt q = 0; + BigInt r = x; + + if (y > x) + return std::pair(q, r); + + BigInt lb = 0; + BigInt ub = x; + while (true) { + BigInt q = (ub + lb).half(); + BigInt r = x - y*q; + + if (r < 0) + ub = q; + else if (r >= y) + lb = q+1; + else + return std::pair(q, r); + } + } +}; + +constexpr uint64_t abs64(int64_t x) { + return static_cast(x > 0 ? x : -x); +} + +constexpr uint64_t pow10(uint64_t e) { + if (e == 0) + return 1; + else + return 10 * pow10(e-1); +} + +#endif diff --git a/code/cpp/ecm.cpp b/code/cpp/ecm.cpp new file mode 100644 index 0000000..64dcb6b --- /dev/null +++ b/code/cpp/ecm.cpp @@ -0,0 +1,91 @@ +#include "bigint.h" +#include "zmodn.h" + +#include +#include +#include + +constexpr BigInt N(NUMBER); + +class Point { +public: + bool is_infinity; + Zmod x; + Zmod y; + + Point(Zmod a, Zmod b) : is_infinity{false}, x{a}, y{b} {} + Point(bool inf) : is_infinity{inf}, x{0}, y{0} {} +}; + +std::variant> +sum_or_factor(BigInt<> a, Point p, Point q) { + if (p.is_infinity) return q; + if (q.is_infinity) return p; + if (p.x == q.x && p.y + q.y == Zmod(0)) return Point(true); + + Zmod l(0); + if (p.x != q.x) + if (auto inv_x = (p.x-q.x).inverse(); inv_x.has_value()) + l = (p.y-q.y) * inv_x.value(); + else + return std::get<0>(extended_gcd(N, (p.x-q.x).toint())); + else + if (auto inv_y = (p.y+q.y).inverse(); inv_y.has_value()) + l = (Zmod(3) * p.x * p.x + a) * inv_y.value(); + else + return std::get<0>(extended_gcd(N, (p.y+q.y).toint())); + + auto x = l*l - p.x - q.x; + auto y = l*(p.x-x) - p.y; + + return Point(x, y); +} + +std::variant> +product_or_factor(BigInt<> a, std::variant> pi, BigInt<> m) { + if (std::holds_alternative>(pi)) + return std::get>(pi); + + auto p = std::get(pi); + + if (m == 0) + return Point(true); // Anything multiplied by 0 is 0 + + // Divide-and-conquer multiplication (power) algorithm + if (m % 2 == 0) { + return product_or_factor(a, sum_or_factor(a, p, p), m/2); + } else { + auto pp = product_or_factor(a, p, m-1); + if (std::holds_alternative>(pp)) + return std::get>(pp); + return sum_or_factor(a, p, std::get(pp)); + } +} + +BigInt<> find_factor_ecm() { + // Find a factor of the integer N. + // If N is prime, this method goes into an infinite loop. + + while (true) { + BigInt a = BigInt<>::random(N); + Point p(BigInt<>::random(N), BigInt<>::random(N)); + for (BigInt m = 2; + (m < 256 || m*m*m*m*m*m*m*m < N) && !p.is_infinity; m += 1) { + auto x = product_or_factor(a, p, m); + if (std::holds_alternative>(x)) + return std::get>(x); + else + p = std::get(x); + } + } +} + +int main() { + // N is a compile-time constant + if (BigInt f = find_factor_ecm(); f > 1 && f < N) + std::cout << N << " = " << f << " * " << N/f << std::endl; + else + std::cout << N << " is prime" << std::endl; + + return 0; +} diff --git a/code/cpp/naive.cpp b/code/cpp/naive.cpp new file mode 100644 index 0000000..ed741e6 --- /dev/null +++ b/code/cpp/naive.cpp @@ -0,0 +1,23 @@ +#include "bigint.h" + +#include +#include + +constexpr BigInt N(NUMBER); + +BigInt<> find_factor() { + for (BigInt i = 2; i*i < N; i += 1) + if (N % i == 0) + return i; + return -1; +} + +int main() { + // N is a compile-time constant + if (auto f = find_factor(); f > 1 && f < N) + std::cout << N << " = " << f << " * " << N/f << std::endl; + else + std::cout << N << " is prime" << std::endl; + + return 0; +} diff --git a/code/cpp/zmodn.h b/code/cpp/zmodn.h new file mode 100644 index 0000000..a24f293 --- /dev/null +++ b/code/cpp/zmodn.h @@ -0,0 +1,79 @@ +#ifndef ZMODN_H +#define ZMODN_H + +#include +#include +#include +#include +#include + +template +concept Integer = requires(T a, T b, int i, std::ostream& os) { + {T(i)}; + + {a + b} -> std::same_as; + {a - b} -> std::same_as; + {a * b} -> std::same_as; + {a / b} -> std::same_as; + {a % b} -> std::same_as; + + {a == b} -> std::same_as; + {a != b} -> std::same_as; + + {os << a} -> std::same_as; +}; + +template +std::tuple extended_gcd(T a, T b) { + if (b == 0) return {a, 1, 0}; + auto [g, x, y] = extended_gcd(b, a%b); + return {g, y, x - y*(a/b)}; +} + +template +requires(N > 1) +class Zmod { +public: + Zmod(decltype(N) z) : value{(z%N + N) % N} {} + decltype(N) toint() const { return value; } + + Zmod operator+(const Zmod& z) const { return value + z.value; } + Zmod operator-(const Zmod& z) const { return value - z.value; } + Zmod operator*(const Zmod& z) const { return value * z.value; } + + Zmod operator+=(const Zmod& z) { return (*this) = value + z.value; } + Zmod operator-=(const Zmod& z) { return (*this) = value - z.value; } + Zmod operator*=(const Zmod& z) { return (*this) = value * z.value; } + + Zmod operator^(decltype(N) z) const { + if (z == 0) return 1; + if (z % 2 == 0) return (((*this) * (*this)) ^ (z/2)); + return (*this) * ((*this) ^ (z-1)); + } + + bool operator==(const Zmod& z) const { return value == z.value; } + bool operator!=(const Zmod& z) const { return value != z.value; } + + std::optional inverse() const { + auto [g, a, _] = extended_gcd(value, N); + return g == 1 ? Zmod(a) : std::optional{}; + } + + std::optional operator/(const Zmod& d) const { + auto i = d.inverse(); + return i ? (*this) * i.value() : i; + } + + std::optional operator/=(const Zmod& d) { + auto q = *this / d; + return q ? (*this = q.value()) : q; + } + + friend std::ostream& operator<<(std::ostream& os, const Zmod& z) { + return os << "(" << z.value << " mod " << N << ")"; + } +private: + decltype(N) value; +}; + +#endif diff --git a/code/python/ecm.py b/code/python/ecm.py new file mode 100755 index 0000000..662facf --- /dev/null +++ b/code/python/ecm.py @@ -0,0 +1,75 @@ +#! /bin/env python + +from sys import argv +from random import randint +from math import sqrt +from dataclasses import dataclass + +@dataclass +class Point: + x: int = 0 + y: int = 0 + is_zero: bool = False + +@dataclass +class FactorFound(Exception): + factor: int = 0 + +# Returns gcd(a, b) and x, y such that ax + by = gcd(a,b) +def extended_gcd(a: int, b: int) -> (int, int, int): + if b == 0: + return a, 1, 0 + g, x, y = extended_gcd(b, a % b) + return g, y, x - y * (a // b) + +def inverse_modulo(a: int, n: int) -> int: + g, x, y = extended_gcd(a, n) + if g != 1: + raise FactorFound(g) + return x + +def ec_sum(p: Point, q: Point, A: int, n: int) -> Point: + if p.is_zero: + return q + if q.is_zero: + return p + if (p.x - q.x) % n == 0 and (p.y + q.y) % n == 0: + return Point(is_zero = True) + + if (p.x - q.x) % n != 0: + k = ((p.y - q.y) * inverse_modulo(p.x - q.x, n)) % n + else: + k = ((3 * p.x**2 + A) * inverse_modulo(p.y + q.y, n)) % n + + x = (k**2 - p.x - q.x) % n + y = (k * (p.x - x) - p.y) % n + + return Point(x = x, y = y) + +def ec_mul(M: int, p: Point, A: int, n: int) -> Point: + if M == 0: + return Point(is_zero = True) + if M % 2 == 0: + return ec_mul(M // 2, ec_sum(p, p, A, n), A, n) + return ec_sum(p, ec_mul(M - 1, p, A, n), A, n) + +# Elliptic curve factorization method +# If n is prime, this method goes into an infinite loop +def find_factor(n: int) -> int: + bound = max(int(sqrt(sqrt(sqrt(n)))), 256) + + while True: + A = randint(0,n) + p = Point(x = randint(0,n), y = randint(0,n)) + M = 2 + while M < bound and not p.is_zero: + try: + p = ec_mul(M, p, A, n) + except FactorFound as ff: + #print("Found with A =", A, "M =", M, "and P =", p) + return ff.factor + M += 1 + +N = int(argv[-1]) +f = find_factor(N) +print(N, "=", f, "*", N//f) diff --git a/code/python/naive.py b/code/python/naive.py new file mode 100755 index 0000000..7f6ec00 --- /dev/null +++ b/code/python/naive.py @@ -0,0 +1,18 @@ +#!/bin/env python + +from sys import argv +from math import floor, sqrt + +def naive_factor(n: int) -> int: + for i in range(2,floor(sqrt(n))+1): + if n%i == 0: + return i + else: + return -1 + +N = int(argv[-1]) +f = naive_factor(N) +if f == -1: + print(N, "is prime") +else: + print(N, "=", f, "*", N//f) diff --git a/images/beer.jpg b/images/beer.jpg new file mode 100644 index 0000000..5a67225 Binary files /dev/null and b/images/beer.jpg differ diff --git a/images/clock.png b/images/clock.png new file mode 100644 index 0000000..631ee42 Binary files /dev/null and b/images/clock.png differ diff --git a/images/clock2.png b/images/clock2.png new file mode 100644 index 0000000..bc2706b Binary files /dev/null and b/images/clock2.png differ diff --git a/images/demo.jpg b/images/demo.jpg new file mode 100644 index 0000000..eaf2cd9 Binary files /dev/null and b/images/demo.jpg differ diff --git a/images/ec1.png b/images/ec1.png new file mode 100644 index 0000000..1c31a1e Binary files /dev/null and b/images/ec1.png differ diff --git a/images/ec2.png b/images/ec2.png new file mode 100644 index 0000000..c44f3c0 Binary files /dev/null and b/images/ec2.png differ diff --git a/images/ec3.png b/images/ec3.png new file mode 100644 index 0000000..6db47e4 Binary files /dev/null and b/images/ec3.png differ diff --git a/images/euclid.png b/images/euclid.png new file mode 100644 index 0000000..ccfec81 Binary files /dev/null and b/images/euclid.png differ diff --git a/images/factorization.webp b/images/factorization.webp new file mode 100644 index 0000000..892a53d Binary files /dev/null and b/images/factorization.webp differ diff --git a/images/number-line.svg b/images/number-line.svg new file mode 100644 index 0000000..bc765a7 --- /dev/null +++ b/images/number-line.svg @@ -0,0 +1,262 @@ + + + + + + + + + + + + style="overflow:visible"> + + + + + + + + + + + + + + + + + + + + + + + + + + image/svg+xml + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/images/numbers.jpg b/images/numbers.jpg new file mode 100644 index 0000000..9b0f835 Binary files /dev/null and b/images/numbers.jpg differ diff --git a/images/questions.png b/images/questions.png new file mode 100644 index 0000000..eee3df2 Binary files /dev/null and b/images/questions.png differ diff --git a/images/sum-1a.png b/images/sum-1a.png new file mode 100644 index 0000000..61fb531 Binary files /dev/null and b/images/sum-1a.png differ diff --git a/images/sum-1b.png b/images/sum-1b.png new file mode 100644 index 0000000..f3a7972 Binary files /dev/null and b/images/sum-1b.png differ diff --git a/images/sum-1c.png b/images/sum-1c.png new file mode 100644 index 0000000..c534cce Binary files /dev/null and b/images/sum-1c.png differ diff --git a/images/sum-2a.png b/images/sum-2a.png new file mode 100644 index 0000000..d925e7a Binary files /dev/null and b/images/sum-2a.png differ diff --git a/images/sum-2b.png b/images/sum-2b.png new file mode 100644 index 0000000..dbd6f03 Binary files /dev/null and b/images/sum-2b.png differ diff --git a/images/sum-2c.png b/images/sum-2c.png new file mode 100644 index 0000000..4637bb8 Binary files /dev/null and b/images/sum-2c.png differ diff --git a/images/sum-3a.png b/images/sum-3a.png new file mode 100644 index 0000000..f85d7ff Binary files /dev/null and b/images/sum-3a.png differ diff --git a/images/sum-3b.png b/images/sum-3b.png new file mode 100644 index 0000000..537c5b2 Binary files /dev/null and b/images/sum-3b.png differ diff --git a/images/sum-3c.png b/images/sum-3c.png new file mode 100644 index 0000000..ca093dc Binary files /dev/null and b/images/sum-3c.png differ diff --git a/images/sum-4a.png b/images/sum-4a.png new file mode 100644 index 0000000..82eee13 Binary files /dev/null and b/images/sum-4a.png differ diff --git a/images/sum-4b.png b/images/sum-4b.png new file mode 100644 index 0000000..827c0c2 Binary files /dev/null and b/images/sum-4b.png differ diff --git a/images/sum-4c.png b/images/sum-4c.png new file mode 100644 index 0000000..bf323fd Binary files /dev/null and b/images/sum-4c.png differ diff --git a/index.html b/index.html new file mode 100644 index 0000000..29d9e58 --- /dev/null +++ b/index.html @@ -0,0 +1,622 @@ + + + + Elliptic Curves and the ECM algorithm + + + + + + + + + + + + + + + + +
+ +

+ +

Elliptic curves and the ECM algorithm

+

Sebastiano Tronto

+

ALTEN Scientific Software Evening

+ +
+ +
+

+

Part I: Numbers

+

+
+ +
+

The integers

+ +The number line + +
    +
  • Operations: sum \(+\), difference \(-\) and multiplication \(\times\)
    • +
  • Various properties: associativity, commutativity, etc...
    • +
  • What about division (without remainder)?
    +If \(\frac ab\) is an integer we say that +\(a\) divides \(b\) (in code: +a % b == 0)
    • +
+
+ +
+

The integers modulo N

+ +
+ +
    + +
  • Two numbers are the same if they give the same remainder +when divided by N
    • +
  • Think of int, but with % N +after every operation
    • +
  • Examples with \(N=12\): +\[9+5\equiv 14\equiv 2\pmod{12}\] +\[7-11\equiv-4\equiv 8\pmod{12}\] +\[3\times 4\equiv 12\equiv 0\pmod{12}\] +
    • +
+ +The number clock + +
+ +
+ +
+

The integers modulo N - Division

+ +

What about division?

+ +
    +
  • +Sometimes it works +\[ +\frac 37\equiv 9 \pmod{12} \qquad +\text{because} \qquad 9\times 7 \equiv 63 \equiv 3 \pmod{12} +\] +
    • + +
  • +Sometimes it does not +\[ +\frac 32\equiv \; ? \pmod{4} \qquad {\color{red}\text{Impossible!}} +\] +
    • +
+
+ +
+

Integers modulo N - Division

+ +
    +
  • +Sometimes it's... weird? +\[ +\frac 62 \equiv \; ? \pmod{8} \quad +\begin{array}{l} +\rightarrow {\color{red}3} \times 2 \equiv 6\pmod{8}\\ +\rightarrow {\color{red}7} \times 2 \equiv 14 \equiv 6 \pmod{8} +\end{array} +\] +
    • +
+ +
+

Integers modulo N - Division

+
+
    +
  • Can divide by \(a\) when \[\operatorname{GCD}(a, N)=1\]
    • +
  • With the + +extended GCD algorithm find \(x\) and \(y\) such that +\[ +ax+Ny=1 +\] +
  • +This means \(\frac{1}{a}\equiv x\pmod{N}\) +
    • +
+

+def extended_gcd(a, b):
+    if b == 0:
+        return a, 1, 0
+    g, x, y = extended_gcd(b, a % b)
+    return g, y, x - y*(a // b)
+
+
+

+Division always works if \(N\) is a prime number!

+
+ +
+

Modular arithmetic - recap

+
+
    +
  • Integers modulo \(N\) are (almost) like numbers
    • +
  • Normal operations like \(+\), \(-\) and \(\times\) work
    • +
  • Division sometimes works, sometimes not
    • +
  • Division always works if \(N\) is prime
    • +
+The number line +
+
+ +
+

+

Part II: Elliptic Curves

+

+
+ +
+

Elliptic curves

+ +
+ +

+An elliptic curve is a curve with equation +\[ y^2 = x^3+Ax+B \] +Where \(A\) and \(B\) are numbers with \[4A^3+27B^2\neq 0\] +

+ +
+
\(y^2=x^3-x+1\)
\((A=-1, B=1\))
+An elliptic curve +
+ +
+
+ +
+

Elliptic curves

+ +
+ +
+
\(y^2=x^3+13x-34\)
\((A=13, B=-34\))
+An elliptic curve +
+ +
+
\(y^2=x^3-x\)
\((A=-1, B=0\))
+An elliptic curve +
+ +
+
+ +
+

Elliptic curves

+ +
+
    +
  • There is a "sum" operation for points of a curve +(NOT the sum of coordinates)
    • +
  • To make things work out nicely, we pretend the curve has +a point at infinity that acts as \(0\): +\[P+0=0\qquad 0+P=0\qquad P-P=0\]
    • +
+ +Elliptic curve sum + +
+
+ +
+

Elliptic curves - sum operation - example 1

+Elliptic curve sum +
+
+

Elliptic curves - sum operation - example 1

+Elliptic curve sum +
+
+

Elliptic curves - sum operation - example 1

+Elliptic curve sum +
+ +
+

Elliptic curves - sum operation - example 2

+Elliptic curve sum +
+
+

Elliptic curves - sum operation - example 2

+Elliptic curve sum +
+
+

Elliptic curves - sum operation - example 2

+Elliptic curve sum +
+ +
+

Elliptic curves - sum operation - example 3

+Elliptic curve sum +
+
+

Elliptic curves - sum operation - example 3

+Elliptic curve sum +
+
+

Elliptic curves - sum operation - example 3

+Elliptic curve sum +
+ +
+

Elliptic curves - sum operation - code

+ +
+

+# Computes p+q on the elliptic curve y^2 = x^3 + Ax + B
+def ec_sum(p: Point, q: Point, A: double) -> Point:
+    if p.is_zero:
+        return q
+    if q.is_zero:
+        return p
+    if p.x == q.x and p.y == -q.y:
+        return Point(is_zero = True)
+
+    if p.x != q.x:
+        k = (p.y - q.y) / (p.x - q.x)
+    else:
+        k = (3 * p.x**2 + A) / (p.y + q.y)
+
+    new_x = k**2 - p.x - q.x
+    new_y = k * (p.x - new_x) - p.y
+    return Point(x = new_x, y = new_y)
+
+ +

+@dataclass
+class Point:
+    x: int = 0
+    y: int = 0
+    is_zero: bool = False
+
+ +
+
+ +
+

Elliptic curves - recap

+
+
    +
  • Curves of equation \(y^2=x^3+Ax+B\)
    • +
  • Can "sum" points of the same curve
    • +
  • Nice properties: associativity, commutativity...
    • +
  • The sum operation can be easily implemented
    • +
+The number line +
+
+ +
+

+

Part III: The Elliptic Curve Factorization Method

+

+
+ +
+

Integer factorization

+

+Every positive integer can be written as the product of prime numbers +

+
+
    +
  • Example: \(69420 = 2\times 2\times 3\times 5\times 13\times 89\)
    • +
  • Computationally hard
    • +
  • Important for cryptography!
    • +
+ +
+
+ +
+

Integer factorization - high-level procedure

+ +
+

+# Returns the list of prime factors of n
+def factorize(n: int) -> list:
+    if n == 1:
+        return []
+
+    if is_prime(n):
+        return [n]
+
+    f = find_factor(n)
+
+    return factorize(n) + factorize(n//f)
+
+ +
    +
  • is_prime(n) can be implemented +efficiently +(Miller-Rabin, +AKS +or ECPP) +
    • +
  • Naive implementation of find_factor(n): +
    
+def find_factor(n: int) -> int:
+    for i in range(2,floor(sqrt(n))+1):
+        if n % i == 0:
+            return i
+
    +
    • +
+ +
+
+ +
+

Find Factor - Elliptic Curve Method

+
+

To find a factor of \(n\):

+

    +
  1. Take a random Elliptic Curve \(E\) +and a random point \(P\) of \(E\)
    2. +
  2. Take a suitable number \(m\)
    3. +
  3. Try to compute \(m\cdot P = P+P+\cdots+P\quad\) (\(m\) times) +with coordinates modulo \(n\)
    4. +
  4. If you attempt an impossible division by some number \(d\), +return \(\operatorname{GCD}(n,d)\)
    5. +
  5. Go back to 1.
    6. +
+
+
+ +
+

Find Factor - Elliptic Curve Method - Example

+
    +
  • Take \(n = 91\)
    • +
  • Take \(E: y^2 = x^3 + 51x -371\) and \(P = (11, 39)\) and \(M = 2\)
    • +
  • +Try to compute \(M\cdot P=P+P\pmod n\): +\[ k = \frac{3x_p^2+51}{2y_p} \pmod n\] +Is \(2y_p=78\) invertible modulo \(91\)? +\[ \operatorname{GCD}(78, 91) = 13 \neq 1 \quad \implies \quad \text{NO!} \] +
    • +
  • Found factor: \(13\)
    • +
+ +
+

+

Demo time!

+git.tronto.net/ecm +

+
+ +
+

Elliptic Curve Method - Questions

+
+

+Q: Aren't we just computing the \(\operatorname{GCD}\) with random numbers? +

+

+A: Yes, but Elliptic Curve operations produce "good candidates" +for these random numbers. +

+
+
+ +
+

Elliptic Curve Method - Questions

+
+

Q: Can we do the same without elliptic curves?

+

A: Yes, with + +Pollard's \(p-1\) Algorithm, but ECM is faster.

+
+
+ +
+

Elliptic Curve Method - Questions

+
+

+Q: Are there objects that are more complicated than Elliptic Curves +and can make the method even faster? +

+

A: Yes, there are higher-dimensional + +Abelian Varieties and other + +Algebraic Groups, but they are much harder (if not impossible) +to implement efficiently.

+
+
+ +
+

+

More questions?

+

+
+ +
+

+

Drinks!

+

+
+ + + + + -- cgit v1.3