From 7f1b8a358515b45d717c3c4b5baf2fbc0f568170 Mon Sep 17 00:00:00 2001 From: Sebastiano Tronto Date: Wed, 26 Feb 2025 17:12:06 +0100 Subject: Initial commit --- 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 ++++++++++++++++ 5 files changed, 459 insertions(+) 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 (limited to 'code/cpp') 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 -- cgit v1.3