Initial commit

7 files changed, 552 insertions, 0 deletions

diff --git a/code/cpp/a.out b/code/cpp/a.out
new file mode 100644
index 0000000..b52f011
--- /dev/null
+++ b/code/cpp/a.out
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 <cstdint>
+#include <iostream>
+#include <random>
+#include <string_view>
+
+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<uint64_t N = 50, uint64_t E = 9>
+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<uint64_t>(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<int> 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<N, E>& 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<BigInt, BigInt>
+        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<BigInt, BigInt>
+        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<uint64_t>(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 <cstdint>
+#include <iostream>
+#include <variant>
+
+constexpr BigInt N(NUMBER);
+
+class Point {
+public:
+        bool is_infinity;
+        Zmod<N> x;
+        Zmod<N> y;
+
+        Point(Zmod<N> a, Zmod<N> b) : is_infinity{false}, x{a}, y{b} {}
+        Point(bool inf) : is_infinity{inf}, x{0}, y{0} {}
+};
+
+std::variant<Point, BigInt<>>
+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<N>(0)) return Point(true);
+
+        Zmod<N> 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<N>(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<Point, BigInt<>>
+product_or_factor(BigInt<> a, std::variant<Point, BigInt<>> pi, BigInt<> m) {
+        if (std::holds_alternative<BigInt<>>(pi))
+                return std::get<BigInt<>>(pi);
+
+        auto p = std::get<Point>(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<BigInt<>>(pp))
+                        return std::get<BigInt<>>(pp);
+                return sum_or_factor(a, p, std::get<Point>(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<BigInt<>>(x))
+                                return std::get<BigInt<>>(x);
+                        else
+                                p = std::get<Point>(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 <cstdint>
+#include <iostream>
+
+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 <cstdint>
+#include <iostream>
+#include <optional>
+#include <tuple>
+#include <type_traits>
+
+template<typename T>
+concept Integer = requires(T a, T b, int i, std::ostream& os) {
+        {T(i)};
+
+        {a + b} -> std::same_as<T>;
+        {a - b} -> std::same_as<T>;
+        {a * b} -> std::same_as<T>;
+        {a / b} -> std::same_as<T>;
+        {a % b} -> std::same_as<T>;
+
+        {a == b} -> std::same_as<bool>;
+        {a != b} -> std::same_as<bool>;
+
+        {os << a} -> std::same_as<std::ostream&>;
+};
+
+template<Integer T>
+std::tuple<T, T, T> 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<Integer auto N>
+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<Zmod> inverse() const {
+                auto [g, a, _] = extended_gcd(value, N);
+                return g == 1 ? Zmod(a) : std::optional<Zmod>{};
+        }
+
+        std::optional<Zmod> operator/(const Zmod& d) const {
+                auto i = d.inverse();
+                return i ? (*this) * i.value() : i;
+        }
+
+        std::optional<Zmod> operator/=(const Zmod& d) {
+                auto q = *this / d;
+                return q ? (*this = q.value()) : q;
+        }
+
+        friend std::ostream& operator<<(std::ostream& os, const Zmod<N>& 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)