diff --git a/CMakeLists.txt b/CMakeLists.txt
index 39a0fe4..21ed6d7 100644
--- a/CMakeLists.txt
+++ b/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(types src/types.cpp)
+add_executable(multivar src/multivar.cpp)
diff --git a/src/multivar.cpp b/src/multivar.cpp
new file mode 100644
index 0000000..46f6ab8
--- /dev/null
+++ b/src/multivar.cpp
@@ -0,0 +1,147 @@
+#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;
+}
\ No newline at end of file
diff --git a/src/types.cpp b/src/types.cpp
deleted file mode 100644
index 6261624..0000000
--- a/src/types.cpp
+++ /dev/null
@@ -1,76 +0,0 @@
-/*
- * 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;
-}
diff --git a/src/types.hpp b/src/types.hpp
new file mode 100644
index 0000000..18be62d
--- /dev/null
+++ b/src/types.hpp
@@ -0,0 +1,119 @@
+#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;
+}
\ No newline at end of file