Added BigInt

4 files changed, 621 insertions, 7 deletions

diff --git a/README.md b/README.md
index 900a296..043b626 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,17 @@
 # ZmodN - A simple library for modular arithmetic
 
+## ZmodN
+
 Usage:
 
 1. `#include "zmodn.h"` in your project
 2. enjoy
 
-# Development
+## BigInt
+
+This repo contains also a class for big integer. Performance are not great,
+but you have a look at it for educational purposes.
+
+## Development
 
 Run `chmod +x test` and then `./test` to run tests.
diff --git a/bigint.h b/bigint.h
new file mode 100644
index 0000000..a8199b9
--- /dev/null
+++ b/bigint.h
@@ -0,0 +1,245 @@
+#ifndef BIGUNSIGNED_H
+#define BIGUNSIGNED_H
+
+#include <cstdint>
+#include <iostream>
+#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; }
+
+        friend std::ostream& operator<<(std::ostream& os, const BigInt<N, E>& z) {
+                bool fl = false;
+                if (!z.sign)
+                        os << "-";
+                for (int i = z.D-1; i >= 0; i--)
+                        if (fl = fl || z.digits[i] != 0; fl)
+                                os << z.digits[i];
+                if (z == 0)
+                        os << "0";
+                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/test b/test
index 4b426c2..46ed7c3 100755
--- a/test
+++ b/test
@@ -2,13 +2,16 @@
 
 cc=${CC:-g++}
 bin="$(mktemp)"
-${cc} -x c++ -std=c++20 -o "$bin" "$(realpath $0)"
+${cc} -x c++ -std=c++20 -o "$bin" -g -O0 "$(realpath $0)"
+echo "Running $bin"
 "$bin"
 
 exit 0
 #endif
 
 #include "zmodn.h"
+#include "bigint.h"
+
 #include <concepts>
 #include <functional>
 #include <iostream>
@@ -141,6 +144,350 @@ public:
                 assert_equal(p, 92);
         }
 },
+{
+        .name = "BigInt constructor zero",
+        .f = []() {
+                BigInt x;
+                BigInt y(0);
+
+                assert_equal(x, y);
+        }
+},
+{
+        .name = "BigInt constructor one digit",
+        .f = []() {
+                BigInt x(12345);
+                BigInt y("12345");
+
+                assert_equal(x, y);
+        }
+},
+{
+        .name = "BigInt constructor many small digits",
+        .f = []() {
+                BigInt<20, 2> x(123456789);
+                BigInt<20, 2> y("123456789");
+
+                assert_equal(x, y);
+        }
+},
+{
+        .name = "BigInt constructor negative many small digits",
+        .f = []() {
+                BigInt<20, 2> x(-123456789);
+                BigInt<20, 2> y("-123456789");
+
+                assert_equal(x, y);
+        }
+},
+{
+        .name = "BigInt operator==",
+        .f = []() {
+                BigInt<20, 2> x(123456789);
+                BigInt<20, 2> y("123456789");
+                BigInt<20, 2> z("12456789");
+
+                bool eq = (x == y);
+                bool diff = (x == z);
+
+                assert_equal(eq, true);
+                assert_equal(diff, false);
+        },
+},
+{
+        .name = "BigInt operator== negative",
+        .f = []() {
+                BigInt<20, 2> x("-123456789");
+                BigInt<20, 2> z("123456789");
+
+                bool diff = (x == z);
+
+                assert_equal(diff, false);
+        },
+},
+{
+        .name = "BigInt operator!= true",
+        .f = []() {
+                BigInt<20, 2> x(12345678);
+                BigInt<20, 2> y("123456789");
+                BigInt<20, 2> z("123456789");
+
+                bool diff = (x != y);
+                bool eq = (y != z);
+
+                assert_equal(diff, true);
+                assert_equal(eq, false);
+        },
+},
+{
+        .name = "BigInt operator< and operator>",
+        .f = []() {
+                BigInt<20, 2> x(7891);
+                BigInt<20, 2> y(7881);
+        
+                bool t = (y < x);
+                bool f = (x < y);
+
+                assert_equal(t, true);
+                assert_equal(f, false);
+        }
+},
+{
+        .name = "BigInt operator< both negative",
+        .f = []() {
+                BigInt<20, 2> x(-7891);
+                BigInt<20, 2> y(-7881);
+        
+                bool cmp = (x < y);
+
+                assert_equal(cmp, true);
+        }
+},
+{
+        .name = "BigInt operator< different sign",
+        .f = []() {
+                BigInt<20, 2> x(-7);
+                BigInt<20, 2> y(7);
+        
+                bool cmp = (x < y);
+
+                assert_equal(cmp, true);
+        }
+},
+{
+        .name = "BigInt abs",
+        .f = []() {
+                BigInt<20, 2> x(-1234567);
+                BigInt<20, 2> y(7654321);
+
+                assert_equal(x.abs(), BigInt<20, 2>(1234567));
+                assert_equal(y.abs(), y);
+        }
+},
+{
+        .name = "BigInt opposite",
+        .f = []() {
+                BigInt<20, 2> x(-1234567);
+                BigInt<20, 2> y(7654321);
+
+                assert_equal(-x, BigInt<20, 2>(1234567));
+                assert_equal(-y, BigInt<20, 2>(-7654321));
+        }
+},
+{
+        .name = "BigInt -0 == 0",
+        .f = []() {
+                BigInt z(0);
+
+                assert_equal(-z, z);
+        }
+},
+{
+        .name = "BigInt sum",
+        .f = []() {
+                BigInt<20, 2> x("987608548588589");
+                BigInt<20, 2> y("6793564545455289");
+                BigInt<20, 2> z("7781173094043878");
+
+                assert_equal(x+y, z);
+        }
+},
+{
+        .name = "BigInt sum both negative",
+        .f = []() {
+                BigInt<20, 2> x("-987608548588589");
+                BigInt<20, 2> y("-6793564545455289");
+                BigInt<20, 2> z("-7781173094043878");
+
+                assert_equal(x+y, z);
+        }
+},
+{
+        .name = "BigInt sum negative + positive, result positive",
+        .f = []() {
+                BigInt<20, 2> x("-987608548588589");
+                BigInt<20, 2> y("6793564545455289");
+                BigInt<20, 2> z("5805955996866700");
+
+                assert_equal(x+y, z);
+        }
+},
+{
+        .name = "BigInt sum positive + negative, result negative",
+        .f = []() {
+                BigInt<20, 2> x("987608548588589");
+                BigInt<20, 2> y("-6793564545455289");
+                BigInt<20, 2> z("-5805955996866700");
+
+                assert_equal(x+y, z);
+        }
+},
+{
+        .name = "BigInt difference",
+        .f = []() {
+                BigInt<20, 2> x("2342442323434134");
+                BigInt<20, 2> y("2524342523342342");
+                BigInt<20, 2> z("-181900199908208");
+
+                assert_equal(x-y, z);
+        }
+},
+{
+        .name = "BigInt product",
+        .f = []() {
+                BigInt<100, 3> x("134142345244134");
+                BigInt<100, 3> y("-56543047058245");
+                BigInt<100, 3> z("-7584816939642416135042584830");
+
+                assert_equal(x*y, z);
+        }
+},
+{
+        .name = "BigInt operator+=",
+        .f = []() {
+                BigInt<20, 2> x("987608548588589");
+                BigInt<20, 2> y("6793564545455289");
+                BigInt<20, 2> z("7781173094043878");
+
+                x += y;
+
+                assert_equal(x, z);
+        }
+},
+{
+        .name = "BigInt 14 / 3",
+        .f = []() {
+                BigInt x(14);
+                BigInt y(3);
+
+                assert_equal(x / y, 4);
+        }
+},
+{
+        .name = "BigInt 14 % 3",
+        .f = []() {
+                BigInt x(14);
+                BigInt y(3);
+
+                assert_equal(x % y, 2);
+        }
+},
+{
+        .name = "BigInt 14 / -3",
+        .f = []() {
+                BigInt x(14);
+                BigInt y(-3);
+
+                assert_equal(x / y, -5);
+        }
+},
+{
+        .name = "BigInt 14 % -3",
+        .f = []() {
+                BigInt x(14);
+                BigInt y(-3);
+
+                assert_equal(x % y, -1);
+        }
+},
+{
+        .name = "BigInt -14 / 3",
+        .f = []() {
+                BigInt x(-14);
+                BigInt y(3);
+
+                assert_equal(x / y, -5);
+        }
+},
+{
+        .name = "BigInt -14 % 3",
+        .f = []() {
+                BigInt x(-14);
+                BigInt y(3);
+
+                assert_equal(x % y, 1);
+        }
+},
+{
+        .name = "BigInt -14 / -3",
+        .f = []() {
+                BigInt x(-14);
+                BigInt y(-3);
+
+                assert_equal(x / y, 4);
+        }
+},
+{
+        .name = "BigInt -14 % -3",
+        .f = []() {
+                BigInt x(-14);
+                BigInt y(-3);
+
+                assert_equal(x % y, -2);
+        }
+},
+{
+        .name = "BigInt division large numbers, quotient = 0",
+        .f = []() {
+                BigInt<50, 3> x("4534435234134244242");
+                BigInt<50, 3> y("7832478748237487343");
+
+                assert_equal(x / y, 0);
+        }
+},
+{
+        .name = "BigInt division large numbers",
+        .f = []() {
+                BigInt<50, 3> x("12344534435234134244242");
+                BigInt<50, 3> y("7832478748237487343");
+                BigInt<50, 3> z(1576);
+
+                assert_equal(x / y, z);
+        }
+},
+{
+        .name = "BigInt modulo large numbers",
+        .f = []() {
+                BigInt<50, 3> x("12344534435234134244242");
+                BigInt<50, 3> y("7832478748237487343");
+                BigInt<50, 3> z("547928011854191674");
+
+                assert_equal(x % y, z);
+        }
+},
+{
+        .name = "Zmod with BigInt constructor",
+        .f = []() {
+                constexpr BigInt<50, 3> N("78923471");
+                constexpr BigInt<50, 3> x("145452451");
+                Zmod<N> xmodN(x);
+
+                assert_equal(xmodN.toint(), x % N);
+        }
+},
+{
+        .name = "Zmod with BigInt big inverse",
+        .f = []() {
+                constexpr BigInt<50, 3> N("7520824651249795349285");
+                constexpr BigInt<50, 3> x("234589234599896924596");
+                constexpr BigInt<50, 3> expected("5901181270843786267351");
+                Zmod<N> xmodN(x);
+
+                auto inv = xmodN.inverse();
+
+                assert_equal(inv.has_value(), true);
+                assert_equal(inv.value().toint(), expected);
+        }
+},
+/*
+{
+        .name = "This does not compile",
+        .f = []() {
+                constexpr double N = 1.2;
+                Zmod<N> x;
+        }
+}
+*/
 };
 
 int main() {
diff --git a/zmodn.h b/zmodn.h
index fdb1ee6..a24f293 100644
--- a/zmodn.h
+++ b/zmodn.h
@@ -7,16 +7,31 @@
 #include <tuple>
 #include <type_traits>
 
-template<typename INT>
-requires std::is_integral_v<INT>
-std::tuple<INT, INT, INT> extended_gcd(INT a, INT b) {
+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<auto N>
-requires(N > 1) && std::is_integral_v<decltype(N)>
+template<Integer auto N>
+requires(N > 1)
 class Zmod {
 public:
         Zmod(decltype(N) z) : value{(z%N + N) % N} {}