For testing

CMakeLists.txt

@@ -10,4 +10,4 @@ add_compile_options(-Wall -Wextra -Wpedantic -O2 -funroll-loops)
 
 add_executable(karatsuba src/karatsuba.cpp)
 add_executable(pseudorandom src/pseudorandom.cpp)
-add_executable(multivar src/multivar.cpp)
+add_executable(types src/types.cpp)

src/karatsuba.cpp

@@ -124,6 +124,7 @@ template <typename R> vector<R> poly_mult_basic(vector<R> &a, vector<R> &b) {
 }
 
 #define THRESHOLD 16
+// #define THRESHOLD 1
 
 /**
  * A step of the Karatsuba function.

src/multivar.cpp

@@ -1,147 +0,0 @@
-#include <list>
-#include <vector>
-#include <version>
-
-#include "types.hpp"
-
-/**
- * Constraint: degrees should have the same size across the board
- */
-template<typename R>
-struct monomial {
-  R coefficient;
-  vector<size_t> degrees;
-};
-
-template <typename R>
-std::ostream& operator<<(std::ostream& os, const monomial<R>& r) {
-  os << r.coefficient;
-  for (size_t i = 0; i < r.degrees.size(); i++) {
-    os << " " << "X_" << i << "^" << r.degrees[i];
-  }
-
-  return os;
-}
-
-template<typename R>
-int lexi_compare(const monomial<R> &lhs, const monomial<R> &rhs) {
-  for (int i = 0; i < lhs.degrees.size(); i++) {
-    if (lhs.degrees[i] > rhs.degrees[i]) {
-      return 1;
-    }
-    if (lhs.degrees[i] < rhs.degrees[i]) {
-      return -1;
-    }
-  }
-  // All degree terms are the same
-  return 0;
-}
-
-/**
- * Ordered according to lexicographic order
- * @tparam R ring
- */
-template<typename R>
-struct multipoly {
-  list<monomial<R>> monomials;
-};
-
-template <typename R>
-std::ostream& operator<<(std::ostream& os, const multipoly<R>& poly) {
-  if (poly.monomials.empty()) {
-    os << "0";
-  } else {
-    auto ite = poly.monomials.begin();
-    os << *ite;
-    ++ite;
-    for (; ite != poly.monomials.end(); ++ite) {
-      os << " + " << *ite;
-    }
-  }
-
-  return os;
-}
-
-template <typename R>
-multipoly<R> &add_monomial(multipoly<R> &poly, const monomial<R> &monomial) {
-  auto ite = poly.monomials.begin();
-  for (; ite != poly.monomials.end(); ++ite) {
-    int cmp = lexi_compare(monomial, *ite);
-    if (cmp < 0) {
-      poly.monomials.insert(ite, monomial);
-      return poly;
-    }
-    if (cmp == 0) {
-      ite->coefficient = ite->coefficient + monomial.coefficient;
-      if (ite->coefficient == 0) {
-        poly.monomials.erase(ite);
-      }
-      return poly;
-    }
-  }
-
-  // Fail to insert before the very back
-  poly.monomials.insert(ite, monomial);
-  return poly;
-}
-
-template <typename R>
-multipoly<R> operator+(const multipoly<R> &lhs, const multipoly<R> &rhs) {
-  auto result = lhs;
-  for (auto &rhs_mono : rhs.monomials) {
-    add_monomial(result, rhs_mono);
-  }
-  return result;
-}
-
-template <typename R>
-multipoly<R> operator-(const multipoly<R> &lhs, const multipoly<R> &rhs) {
-  auto result = lhs;
-  for (auto &rhs_mono : rhs.monomials) {
-    auto negate_rhs = rhs_mono;
-    negate_rhs.coefficient = -negate_rhs.coefficient;
-    add_monomial(result, negate_rhs);
-  }
-  return result;
-}
-
-template <typename R>
-multipoly<R> operator*(const multipoly<R> &lhs, const multipoly<R> &rhs) {
-  auto result = lhs;
-  // TODO Implement multiplication
-  return result;
-}
-
-template <typename R>
-bool is_zero(multipoly<R> &poly) {
-  return poly.monomials.empty();
-}
-
-int main() {
-  auto m1 = monomial{2, vector<size_t>{3, 4}};
-  auto m2 = monomial{2, vector<size_t>{4, 2}};
-  auto m3 = monomial{3, vector<size_t>{3, 2}};
-  cout << "m1, m2, m3: " << m1 << ", " << m2 << ", " << m3 << endl;
-  cout << "compare m1, m2: " << lexi_compare(m1, m2) << endl;
-  cout << "compare m1, m3: " << lexi_compare(m1, m3) << endl;
-  cout << "compare m1, m1: " << lexi_compare(m1, m1) << endl;
-
-  cout << endl;
-  auto poly = multipoly<int>();
-  cout << "poly - step 0: " << poly << endl;
-  add_monomial(poly, m1);
-  cout << "poly - step 1: " << poly << endl;
-  add_monomial(poly, m2);
-  cout << "poly - step 2: " << poly << endl;
-  add_monomial(poly, m3);
-  cout << "poly - step 3: " << poly << endl;
-  add_monomial(poly, m1);
-  cout << "poly - step 4: " << poly << endl;
-
-  cout << "poly + poly: " << poly + poly << endl;
-  auto sub = poly - poly;
-  cout << "poly - poly: " << sub << endl;
-  cout << "poly - poly is zero: " << is_zero(sub) << endl;
-
-  return 0;
-}

src/types.cpp

@@ -0,0 +1,76 @@
+/*
+ * Until next time, let's boilerplate later
+ */
+#include <cmath>
+#include <iostream>
+using namespace std;
+
+template<typename R>
+struct rational {
+  R num;
+  R denom;
+};
+
+template<typename R>
+rational<R> operator+(rational<R> &l, rational<R> &r) {
+  return rational(l.num * r.denom + r.num * l.denom, l.denom * r.denom);
+}
+
+template<typename R>
+rational<R> operator-(rational<R> &l, rational<R> &r) {
+  return rational(l.num * r.denom - r.num * l.denom, l.denom * r.denom);
+}
+
+template<typename R>
+rational<R> operator*(rational<R> &l, rational<R> &r) {
+  return rational(l.num * r.num, l.denom * r.denom);
+}
+
+template<typename R>
+rational<R> operator/(rational<R> &l, rational<R> &r) {
+  return rational(l.num * r.denom, r.denom * l.num);
+}
+
+// TODO Take care of normalization
+template<typename R>
+bool operator==(rational<R> &l, rational<R> &r) {
+  return l.num == r.num && l.denom == r.denom;
+}
+
+template<typename R>
+bool operator!=(rational<R> &l, rational<R> &r) {
+  return l.num != r.num || l.denom != r.denom;
+}
+
+/*
+template<typename R>
+ostream operator<< (ostream &os, rational<R> &r) {
+
+}
+*/
+
+template<unsigned int>
+struct prime_field {
+  size_t number;
+};
+
+template<unsigned int p>
+prime_field<p> operator+(const prime_field<p> &l, const prime_field<p> &r) {
+  return prime_field(l.number + r.number % p);
+}
+
+template<unsigned int p>
+prime_field<p> operator-(const prime_field<p> &l, const prime_field<p> &r) {
+  return prime_field(l.number - r.number % p);
+}
+
+template<unsigned int p>
+prime_field<p> operator*(const prime_field<p> &l, const prime_field<p> &r) {
+  return prime_field(l.number * r.number % p);
+}
+
+// TODO /, ==, !=
+
+int main() {
+  cout << "Hello World!" << endl;
+}

src/types.hpp

@@ -1,119 +0,0 @@
-#include <cmath>
-#include <iostream>
-using namespace std;
-
-/**
- * A generic function for computing the GCD
- * @tparam R : this template variable should be a type for which all the necessary
- * operations for Euclidean domains have been implemented
- * @param a
- * @param b
- * @return
- */
-template <typename R>
-R gcd(R a, R b) {
-  if (b == 0)
-    return a;
-  R r = a % b;
-  while (r != 0) {
-    a = b;
-    b = r;
-    r = a % b;
-  }
-  return b;
-}
-
-/**
- * Rational on R which is always normalized.
- * @tparam R : a type for Euclidean domain
- */
-template<typename R>
-struct ratio {
-  R num;
-  R denom;
-};
-
-template<typename R>
-ratio<R> &rational_of(R n) {
-  return ratio<R>(n, 1);
-}
-
-template<typename R>
-ratio<R> &normalize(ratio<R> &n) {
-  const R gcd = gcd(n.num, n.denom);
-  n.num = n.num / gcd;
-  n.denom = n.denom / gcd;
-  return n;
-}
-
-template<typename R>
-ratio<R> &operator+(const ratio<R> &l, const ratio<R> &r) {
-  return normalize(rational(l.num * r.denom + r.num * l.denom, l.denom * r.denom));
-}
-
-template<typename R>
-ratio<R> &operator-(const ratio<R> &l, const ratio<R> &r) {
-  return normalize(rational(l.num * r.denom - r.num * l.denom, l.denom * r.denom));
-}
-
-template<typename R>
-ratio<R> &operator*(const ratio<R> &l, const ratio<R> &r) {
-  return normalize(rational(l.num * r.num, l.denom * r.denom));
-}
-
-template<typename R>
-ratio<R> &operator/(const ratio<R> &l, const ratio<R> &r) {
-  return normalize(rational(l.num * r.denom, r.denom * l.num));
-}
-
-template<typename R>
-bool operator==(const ratio<R> &l, const ratio<R> &r) {
-  return l.num == r.num && l.denom == r.denom;
-}
-
-template<typename R>
-bool operator!=(const ratio<R> &l, const ratio<R> &r) {
-  return l.num != r.num || l.denom != r.denom;
-}
-
-template <typename R>
-std::ostream& operator<<(std::ostream& os, const ratio<R>& r) {
-  os << r.a << "/" << r.b;
-
-  return os;
-}
-
-/**
- * Prime field where number is kept to be in range [0, p).
- */
-template<unsigned int>
-struct prime_field {
-  size_t number;
-};
-
-template<unsigned int p>
-prime_field<p> &operator+(const prime_field<p> &l, const prime_field<p> &r) {
-  return prime_field((l.number + r.number) % p);
-}
-
-template<unsigned int p>
-prime_field<p> &operator-(const prime_field<p> &l, const prime_field<p> &r) {
-  return prime_field((l.number - r.number) % p);
-}
-
-template<unsigned int p>
-prime_field<p> &operator*(const prime_field<p> &l, const prime_field<p> &r) {
-  return prime_field((l.number * r.number) % p);
-}
-
-// TODO prime_field /
-
-template<unsigned int p>
-bool operator==(const prime_field<p> &l, const prime_field<p> &r) {
-  return l.number == r.number;
-}
-
-template<unsigned int p>
-bool operator!=(const prime_field<p> &l, const prime_field<p> &r) {
-  return l.number != r.number;
-}