diff --git a/include/poly.h b/include/poly.h
index 65c5bb67..99b43cbf 100644
--- a/include/poly.h
+++ b/include/poly.h
@@ -68,6 +68,7 @@ typedef struct lp_interval_struct lp_interval_t;
 typedef struct lp_feasibility_set_struct lp_feasibility_set_t;
 typedef struct lp_polynomial_hash_set_struct lp_polynomial_hash_set_t;
 typedef struct lp_polynomial_vector_struct lp_polynomial_vector_t;
+typedef struct lp_polynomial_heap_struct lp_polynomial_heap_t;
 
 /** Enable a given tag for tracing */
 void lp_trace_enable(const char* tag);
diff --git a/include/polynomial_hash_set.h b/include/polynomial_hash_set.h
index dde568d7..597e0147 100644
--- a/include/polynomial_hash_set.h
+++ b/include/polynomial_hash_set.h
@@ -32,12 +32,18 @@ struct lp_polynomial_hash_set_struct {
   size_t data_size;
   /** Number of set elements */
   size_t size;
-  /** Treshold for resize */
+  /** Threshold for resize */
   size_t resize_threshold;
   /** Has the set been closed */
   int closed;
 };
 
+/** Allocate and construct a new hash set */
+lp_polynomial_hash_set_t* lp_polynomial_hash_set_new(void);
+
+/** Destructs and deletes the hash set */
+void lp_polynomial_hash_set_delete(lp_polynomial_hash_set_t* set);
+
 /** Construct a new set */
 void lp_polynomial_hash_set_construct(lp_polynomial_hash_set_t* set);
 
@@ -45,17 +51,36 @@ void lp_polynomial_hash_set_construct(lp_polynomial_hash_set_t* set);
 void lp_polynomial_hash_set_destruct(lp_polynomial_hash_set_t* set);
 
 /** Returns true if empty */
-int lp_polynomial_hash_set_is_empty(lp_polynomial_hash_set_t* set);
+int lp_polynomial_hash_set_is_empty(const lp_polynomial_hash_set_t* set);
+
+/** Returns the number of elements */
+size_t lp_polynomial_hash_set_size(const lp_polynomial_hash_set_t* set);
 
 /** Check whether p is in set. The set must not be closed). */
-int lp_polynomial_hash_set_contains(lp_polynomial_hash_set_t* set, const lp_polynomial_t* p);
+int lp_polynomial_hash_set_contains(const lp_polynomial_hash_set_t* set, const lp_polynomial_t* p);
 
-/** Add polynomial p to set. Returns true if p was added (not already in the set).  */
+/** Add polynomial p to the set. Returns true if p was added (not already in the set). */
 int lp_polynomial_hash_set_insert(lp_polynomial_hash_set_t* set, const lp_polynomial_t* p);
 
+/** Add polynomial p to the set. Returns true if p was added (not already in the set) and p becomes a 0 polynomial.
+ *  Returns false if p was found in the set and p remained unchanged. */
+int lp_polynomial_hash_set_insert_move(lp_polynomial_hash_set_t* set, lp_polynomial_t* p);
+
+/** Add all polynomials from the vector to the hash map. Returns the number of inserted polynomials */
+int lp_polynomial_hash_set_insert_vector(lp_polynomial_hash_set_t* set, const lp_polynomial_vector_t* v);
+
+/** Add polynomial p to set. Returns true if p was removed (was in the set). */
+int lp_polynomial_hash_set_remove(lp_polynomial_hash_set_t* set, const lp_polynomial_t* p);
+
+/** Removes all elements from set that are not in other */
+void lp_polynomial_hash_set_intersect(lp_polynomial_hash_set_t* set, const lp_polynomial_hash_set_t* other);
+
 /** Close the set: compact the data so that all elements get stored in data[0..size]. No addition after close! */
 void lp_polynomial_hash_set_close(lp_polynomial_hash_set_t* set);
 
+/** Returns one element at index. Hashset must be closed. */
+const lp_polynomial_t* lp_polynomial_hash_set_at(const lp_polynomial_hash_set_t* set, size_t n);
+
 /** Clear the set. */
 void lp_polynomial_hash_set_clear(lp_polynomial_hash_set_t* set);
 
diff --git a/include/polynomial_heap.h b/include/polynomial_heap.h
new file mode 100644
index 00000000..cd4cc197
--- /dev/null
+++ b/include/polynomial_heap.h
@@ -0,0 +1,88 @@
+/**
+ * Copyright 2015, SRI International.
+ *
+ * This file is part of LibPoly.
+ *
+ * LibPoly is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LibPoly is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with LibPoly.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "poly.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef int (*lp_polynomial_heap_compare_f)(const lp_polynomial_t* A, const lp_polynomial_t* B);
+
+struct lp_polynomial_heap_struct {
+    /** The data */
+    lp_polynomial_t** data;
+    /** Size of the data */
+    size_t data_size;
+    /** Number of set elements */
+    size_t size;
+    /** The compare function */
+    lp_polynomial_heap_compare_f cmp;
+};
+
+/** Allocates a new heap and constructs it */
+lp_polynomial_heap_t* lp_polynomial_heap_new(lp_polynomial_heap_compare_f cmp);
+
+/** Destructs a heap and frees the memory */
+void lp_polynomial_heap_delete(lp_polynomial_heap_t* heap);
+
+/** Construct a new heap */
+void lp_polynomial_heap_construct(lp_polynomial_heap_t* heap, lp_polynomial_heap_compare_f cmp);
+
+/** Destruct the heap */
+void lp_polynomial_heap_destruct(lp_polynomial_heap_t* heap);
+
+/** Returns true if empty */
+int lp_polynomial_heap_is_empty(const lp_polynomial_heap_t* heap);
+
+/** Returns the number of elements */
+size_t lp_polynomial_heap_size(const lp_polynomial_heap_t* heap);
+
+/** Add polynomial p to the heap.  */
+void lp_polynomial_heap_push(lp_polynomial_heap_t* heap, const lp_polynomial_t* p);
+
+/** Add polynomial p to the heap. Moves the content of p and p becomes a 0 polynomial. */
+void lp_polynomial_heap_push_move(lp_polynomial_heap_t* heap, lp_polynomial_t* p);
+
+/** Add all polynomials from the vector to the heap */
+void lp_polynomial_heap_push_vector(lp_polynomial_heap_t* heap, const lp_polynomial_vector_t* v);
+
+/** Removes and returns the top element of the heap, returns NULL if the heap is empty  */
+lp_polynomial_t* lp_polynomial_heap_pop(lp_polynomial_heap_t* heap);
+
+/** Removes an element from the heap. Returns number of removed polynomials. */
+int lp_polynomial_heap_remove(lp_polynomial_heap_t* heap, const lp_polynomial_t *p);
+
+/** Returns the top element without removing it. */
+const lp_polynomial_t* lp_polynomial_heap_peek(const lp_polynomial_heap_t* heap);
+
+/** Returns one element */
+const lp_polynomial_t* lp_polynomial_heap_at(const lp_polynomial_heap_t* heap, size_t n);
+
+/** Clear the heap. */
+void lp_polynomial_heap_clear(lp_polynomial_heap_t* heap);
+
+/** Prints the heap */
+void lp_polynomial_heap_print(lp_polynomial_heap_t* heap, FILE *);
+
+#ifdef __cplusplus
+} /* close extern "C" { */
+#endif
diff --git a/include/polynomial_vector.h b/include/polynomial_vector.h
index a53b2bc9..13f539b4 100644
--- a/include/polynomial_vector.h
+++ b/include/polynomial_vector.h
@@ -28,12 +28,22 @@ extern "C" {
 /** Allocate and construct a new vector */
 lp_polynomial_vector_t* lp_polynomial_vector_new(const lp_polynomial_context_t* ctx);
 
-/* Delete the vector */
+/** Allocate and construct a new vector, copies all elements from v */
+lp_polynomial_vector_t* lp_polynomial_vector_copy(const lp_polynomial_vector_t *v);
+
+/** Delete the vector */
 void lp_polynomial_vector_delete(lp_polynomial_vector_t* v);
 
+/** Swap two vectors */
+void lp_polynomial_vector_swap(lp_polynomial_vector_t *v1, lp_polynomial_vector_t *v2);
+
 /** Add to back (makes a copy, should be in the context of the vector) */
 void lp_polynomial_vector_push_back(lp_polynomial_vector_t* v, const lp_polynomial_t* p);
 
+/** Add to back by move. More efficient than copying, p becomes a (constructed) zero-polynomial
+ *  and must be in the context of the vector. */
+void lp_polynomial_vector_push_back_move(lp_polynomial_vector_t* v, lp_polynomial_t* p);
+
 /** Reset the vector to 0 elements */
 void lp_polynomial_vector_reset(lp_polynomial_vector_t* v);
 
@@ -43,6 +53,12 @@ size_t lp_polynomial_vector_size(const lp_polynomial_vector_t* v);
 /** Returns the polynomial at i (newly constructed each time) */
 lp_polynomial_t* lp_polynomial_vector_at(const lp_polynomial_vector_t* v, size_t i);
 
+/** Returns the context of the polynomial vector */
+const lp_polynomial_context_t* lp_polynomial_vector_get_context(const lp_polynomial_vector_t *v);
+
+/** Prints the polynomial vector to the given stream. */
+void lp_polynomial_vector_print(const lp_polynomial_vector_t* v, FILE* out);
+
 #ifdef __cplusplus
 } /* close extern "C" { */
 #endif
diff --git a/include/upolynomial.h b/include/upolynomial.h
index f7bc9247..51ad71e1 100644
--- a/include/upolynomial.h
+++ b/include/upolynomial.h
@@ -109,6 +109,12 @@ int lp_upolynomial_print(const lp_upolynomial_t* p, FILE* out);
  */
 char* lp_upolynomial_to_string(const lp_upolynomial_t* p);
 
+/**
+ * Converts the upolynomial to a polynomial in variable var.
+ * ctx must have the same ring than upolynomial.
+ */
+lp_polynomial_t* lp_upolynomial_to_polynomial(const lp_upolynomial_t* p, const lp_polynomial_context_t* ctx, lp_variable_t var);
+
 /**
  * Returns true if this is a zero polynomial
  */
diff --git a/include/variable_order.h b/include/variable_order.h
index 73a5a9c5..bc7ba71e 100644
--- a/include/variable_order.h
+++ b/include/variable_order.h
@@ -45,6 +45,22 @@ void lp_variable_order_detach(lp_variable_order_t* var_order);
   */
 int lp_variable_order_cmp(const lp_variable_order_t* var_order, lp_variable_t x, lp_variable_t y);
 
+/** returns the bigger of the two variables wrt. the order */
+static inline
+lp_variable_t lp_variable_order_max(const lp_variable_order_t* var_order, lp_variable_t x, lp_variable_t y) {
+  if (x == lp_variable_null) return y;
+  if (y == lp_variable_null) return x;
+  return lp_variable_order_cmp(var_order, x, y) < 0 ? x : y;
+}
+
+/** returns the smaller of the two variables wrt. the order */
+static inline
+lp_variable_t lp_variable_order_min(const lp_variable_order_t* var_order, lp_variable_t x, lp_variable_t y) {
+  if (x == lp_variable_null) return y;
+  if (y == lp_variable_null) return x;
+  return lp_variable_order_cmp(var_order, x, y) > 0 ? x : y;
+}
+
 /** Get the size of the order */
 size_t lp_variable_order_size(const lp_variable_order_t* var_order);
 
diff --git a/python/polypyPolynomial3.c b/python/polypyPolynomial3.c
index 8273a16d..026a80a8 100644
--- a/python/polypyPolynomial3.c
+++ b/python/polypyPolynomial3.c
@@ -182,7 +182,7 @@ PyMethodDef Polynomial_methods[] = {
     {"coefficients", (PyCFunction)Polynomial_coefficients, METH_NOARGS, "Returns a dictionary from degrees to coefficients"},
     {"reductum", (PyCFunction)Polynomial_reductum, METH_VARARGS, "Returns the reductum of the polynomial"},
     {"sgn", (PyCFunction)Polynomial_sgn, METH_VARARGS, "Returns the sign of the polynomials in the given model"},
-    {"sgn_check", (PyCFunction)Polynomial_sgn_check, METH_VARARGS, "Returns true if the sign of the polynomail respects the sign condition."},
+    {"sgn_check", (PyCFunction)Polynomial_sgn_check, METH_VARARGS, "Returns true if the sign of the polynomial respects the sign condition."},
     {"rem", (PyCFunction)Polynomial_rem, METH_VARARGS, "Returns the remainder of current and given polynomial"},
     {"prem", (PyCFunction)Polynomial_prem, METH_VARARGS, "Returns the pseudo remainder of current and given polynomial"},
     {"sprem", (PyCFunction)Polynomial_sprem, METH_VARARGS, "Returns the sparse pseudo remainder of current and given polynomial"},
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ac20a969..366cb1f6 100755
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -34,6 +34,7 @@ set(poly_SOURCES
   polynomial/polynomial_context.c
   polynomial/feasibility_set.c
   polynomial/polynomial_hash_set.c
+  polynomial/polynomial_heap.c
   polynomial/polynomial_vector.c
   poly.c
 )
diff --git a/src/polynomial/polynomial_hash_set.c b/src/polynomial/polynomial_hash_set.c
index 078d2420..92fa6692 100644
--- a/src/polynomial/polynomial_hash_set.c
+++ b/src/polynomial/polynomial_hash_set.c
@@ -19,6 +19,7 @@
 
 #include <polynomial_hash_set.h>
 #include <polynomial.h>
+#include <polynomial_vector.h>
 
 #include <stdlib.h>
 #include <string.h>
@@ -30,6 +31,17 @@
 /** Resize threshold: the size is doubled when nelems >= size * RESIZE_RATIO */
 #define LP_POLYNOMIAL_HASH_SET_RESIZE_RATIO 0.7
 
+lp_polynomial_hash_set_t* lp_polynomial_hash_set_new(void) {
+  lp_polynomial_hash_set_t *set = malloc(sizeof(lp_polynomial_hash_set_t));
+  lp_polynomial_hash_set_construct(set);
+  return set;
+}
+
+void lp_polynomial_hash_set_delete(lp_polynomial_hash_set_t* set) {
+  lp_polynomial_hash_set_destruct(set);
+  free(set);
+}
+
 void lp_polynomial_hash_set_construct(lp_polynomial_hash_set_t* set) {
   set->data = malloc(LP_POLYNOMIAL_HASH_SET_DEFAULT_SIZE*sizeof(lp_polynomial_t*));
   memset(set->data, 0, LP_POLYNOMIAL_HASH_SET_DEFAULT_SIZE*sizeof(lp_polynomial_t*));
@@ -52,12 +64,28 @@ void lp_polynomial_hash_set_destruct(lp_polynomial_hash_set_t* set) {
   set->data = NULL;
 }
 
-int lp_polynomial_hash_set_is_empty(lp_polynomial_hash_set_t* set) {
+int lp_polynomial_hash_set_is_empty(const lp_polynomial_hash_set_t* set) {
   return set->size == 0;
 }
 
+size_t lp_polynomial_hash_set_size(const lp_polynomial_hash_set_t* set) {
+    return set->size;
+}
+
 static
-void lp_polynomial_hash_set_insert_move(lp_polynomial_t** data, size_t mask, lp_polynomial_t* p) {
+int lp_polynomial_hash_set_insert_move_check(lp_polynomial_t** data, size_t mask, lp_polynomial_t* p) {
+  size_t i = lp_polynomial_hash(p) & mask;
+  while (data[i] != 0) {
+    if (lp_polynomial_eq(data[i], p)) { return 0; }
+    i ++;
+    i &= mask;
+  }
+  data[i] = p;
+  return 1;
+}
+
+static
+void lp_polynomial_hash_set_insert_move_no_check(lp_polynomial_t** data, size_t mask, lp_polynomial_t* p) {
   size_t i = lp_polynomial_hash(p) & mask;
   while (data[i] != 0) {
     i ++;
@@ -79,7 +107,20 @@ int lp_polynomial_hash_set_insert_copy(lp_polynomial_t** data, size_t mask, cons
 }
 
 static
-int lp_polynomial_hash_set_search(lp_polynomial_t** data, size_t mask, const lp_polynomial_t* p) {
+int lp_polynomial_hash_set_insert_swap(lp_polynomial_t** data, size_t mask, lp_polynomial_t* p) {
+  size_t i = lp_polynomial_hash(p) & mask;
+  while (data[i] != 0) {
+    if (lp_polynomial_eq(data[i], p)) { return 0; }
+    i ++;
+    i &= mask;
+  }
+  data[i] = lp_polynomial_new(lp_polynomial_get_context(p));
+  lp_polynomial_swap(data[i], p);
+  return 1;
+}
+
+static
+int lp_polynomial_hash_set_search(lp_polynomial_t* const* data, size_t mask, const lp_polynomial_t* p) {
   size_t i = lp_polynomial_hash(p) & mask;
   for (;;) {
     if (data[i] == 0) return 0;
@@ -90,6 +131,28 @@ int lp_polynomial_hash_set_search(lp_polynomial_t** data, size_t mask, const lp_
   return 0;
 }
 
+#define SWAP(A, B) ({lp_polynomial_t *tmp = A; A = B; B = tmp;})
+
+static
+int lp_polynomial_hash_set_search_and_remove(lp_polynomial_t** data, size_t mask, const lp_polynomial_t* p) {
+  size_t i = lp_polynomial_hash(p) & mask;
+  for (;;) {
+    if (data[i] == 0) return 0;
+    if (lp_polynomial_eq(data[i], p)) break;
+    i ++;
+    i &= mask;
+  }
+  lp_polynomial_delete(data[i]);
+  data[i] = 0;
+  for (;;) {
+    size_t j = (i + 1) & mask;
+    if (data[j] == 0 || (lp_polynomial_hash(data[j]) & mask) == j) break;
+    SWAP(data[i], data[j]);
+    i = j;
+  }
+  return 1;
+}
+
 /** Double the size of the set. */
 static
 void lp_polynomial_hash_set_extend(lp_polynomial_hash_set_t* set) {
@@ -105,7 +168,7 @@ void lp_polynomial_hash_set_extend(lp_polynomial_hash_set_t* set) {
   for (i = 0; i < old_data_size; ++ i) {
     lp_polynomial_t* p = set->data[i];
     if (p != 0) {
-      lp_polynomial_hash_set_insert_move(new_data, mask, p);
+      lp_polynomial_hash_set_insert_move_no_check(new_data, mask, p);
     }
   }
 
@@ -115,7 +178,7 @@ void lp_polynomial_hash_set_extend(lp_polynomial_hash_set_t* set) {
   set->resize_threshold = new_data_size*LP_POLYNOMIAL_HASH_SET_RESIZE_RATIO;
 }
 
-int lp_polynomial_hash_set_contains(lp_polynomial_hash_set_t* set, const lp_polynomial_t* p) {
+int lp_polynomial_hash_set_contains(const lp_polynomial_hash_set_t* set, const lp_polynomial_t* p) {
   assert(p);
   assert(!set->closed);
   return lp_polynomial_hash_set_search(set->data, set->data_size-1, p);
@@ -137,6 +200,80 @@ int lp_polynomial_hash_set_insert(lp_polynomial_hash_set_t* set, const lp_polyno
   return result;
 }
 
+int lp_polynomial_hash_set_insert_move(lp_polynomial_hash_set_t* set, lp_polynomial_t* p) {
+  assert(p);
+  assert(set->data_size > set->size);
+  assert(!set->closed);
+
+  int result = lp_polynomial_hash_set_insert_swap(set->data, set->data_size-1, p);
+  if (result) {
+    set->size ++;
+    if (set->size > set->resize_threshold) {
+      lp_polynomial_hash_set_extend(set);
+    }
+  }
+
+  return result;
+}
+
+int lp_polynomial_hash_set_insert_vector(lp_polynomial_hash_set_t* set, const lp_polynomial_vector_t* v) {
+  assert(v);
+  assert(set->data_size > set->size);
+  assert(!set->closed);
+
+  int inserted = 0;
+  size_t size = lp_polynomial_vector_size(v);
+  for (size_t i = 0; i < size; ++i) {
+    lp_polynomial_t *p = lp_polynomial_vector_at(v, i);
+    int result = lp_polynomial_hash_set_insert_move_check(set->data, set->data_size - 1, p);
+    if (result) {
+      inserted ++;
+      set->size ++;
+      if (set->size > set->resize_threshold) {
+        lp_polynomial_hash_set_extend(set);
+      }
+    } else {
+      lp_polynomial_delete(p);
+    }
+  }
+  return inserted;
+}
+
+int lp_polynomial_hash_set_remove(lp_polynomial_hash_set_t* set, const lp_polynomial_t* p) {
+  assert(p);
+  assert(!set->closed);
+
+  int result = lp_polynomial_hash_set_search_and_remove(set->data, set->data_size-1, p);
+  if (result) {
+    set->size --;
+  }
+
+  return result;
+}
+
+void lp_polynomial_hash_set_intersect(lp_polynomial_hash_set_t* set, const lp_polynomial_hash_set_t* other) {
+  assert(!set->closed);
+
+  size_t mask = set->data_size - 1;
+  for (size_t i = 0; i < set->data_size;) {
+    if (set->data[i] == NULL) {
+      ++i;
+      continue;
+    }
+    if (!lp_polynomial_hash_set_contains(other, set->data[i])) {
+      lp_polynomial_delete(set->data[i]);
+      set->data[i] = NULL;
+      size_t k = i;
+      for (;;) {
+        size_t j = (k + 1) & mask;
+        if (set->data[j] == 0 || (lp_polynomial_hash(set->data[j]) & mask) == j) break;
+        SWAP(set->data[k], set->data[j]);
+        k = j;
+      }
+    }
+  }
+}
+
 void lp_polynomial_hash_set_close(lp_polynomial_hash_set_t* set) {
 
   if (set->closed) {
@@ -160,6 +297,14 @@ void lp_polynomial_hash_set_close(lp_polynomial_hash_set_t* set) {
   assert(i == set->size && i < data_size);
 }
 
+const lp_polynomial_t* lp_polynomial_hash_set_at(const lp_polynomial_hash_set_t* set, size_t n) {
+  assert(set->closed);
+  if (n >= set->size) {
+    return NULL;
+  }
+  return set->data[n];
+}
+
 void lp_polynomial_hash_set_clear(lp_polynomial_hash_set_t* set) {
   lp_polynomial_hash_set_destruct(set);
   lp_polynomial_hash_set_construct(set);
diff --git a/src/polynomial/polynomial_heap.c b/src/polynomial/polynomial_heap.c
new file mode 100644
index 00000000..7bf6fd43
--- /dev/null
+++ b/src/polynomial/polynomial_heap.c
@@ -0,0 +1,197 @@
+/**
+ * Copyright 2015, SRI International.
+ *
+ * This file is part of LibPoly.
+ *
+ * LibPoly is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * LibPoly is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with LibPoly.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <polynomial_heap.h>
+#include <polynomial.h>
+#include <polynomial_vector.h>
+
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+
+/** Default initial size (must be a power of 2) */
+#define LP_POLYNOMIAL_HEAP_DEFAULT_SIZE 32
+
+lp_polynomial_heap_t* lp_polynomial_heap_new(lp_polynomial_heap_compare_f cmp) {
+  lp_polynomial_heap_t *heap = malloc(sizeof(lp_polynomial_heap_t));
+  lp_polynomial_heap_construct(heap, cmp);
+  return heap;
+}
+
+void lp_polynomial_heap_delete(lp_polynomial_heap_t* heap) {
+  lp_polynomial_heap_destruct(heap);
+  free(heap);
+}
+
+void lp_polynomial_heap_construct(lp_polynomial_heap_t *heap, lp_polynomial_heap_compare_f cmp) {
+  heap->data = malloc(LP_POLYNOMIAL_HEAP_DEFAULT_SIZE * sizeof(lp_polynomial_t*));
+  heap->data_size = LP_POLYNOMIAL_HEAP_DEFAULT_SIZE;
+  heap->size = 0;
+  heap->cmp = cmp;
+}
+
+void lp_polynomial_heap_destruct(lp_polynomial_heap_t* heap) {
+  // remove all the polynomials
+  for (size_t i = 0; i < heap->size; ++i) {
+    lp_polynomial_delete(heap->data[i]);
+  }
+  // free the space
+  free(heap->data);
+  heap->data = NULL;
+}
+
+int lp_polynomial_heap_is_empty(const lp_polynomial_heap_t *heap) {
+    return heap->size == 0;
+}
+
+size_t lp_polynomial_heap_size(const lp_polynomial_heap_t* heap) {
+    return heap->size;
+}
+
+#define SWAP(A,B) ({ lp_polynomial_t *tmp = B; B = A; A = tmp; })
+// macros assume a heap indices 1...size
+#define HEAP_SWAP(heap, i, j) SWAP(heap->data[(i)-1], heap->data[(j)-1])
+#define HEAP_CMP(heap, i, j) (heap->cmp(heap->data[(i)-1], heap->data[(j)-1]))
+
+static
+void lp_polynomial_heap_extend(lp_polynomial_heap_t *heap) {
+  // double the size
+  heap->data_size <<= 1;
+  heap->data = realloc(heap->data, heap->data_size * sizeof(lp_polynomial_t*));
+}
+
+static
+void lp_polynomial_heap_heapify_up(lp_polynomial_heap_t *heap) {
+  for (size_t pos = heap->size;
+       pos > 1 && HEAP_CMP(heap, pos / 2, pos) < 0;
+       pos /= 2) {
+    // if data[pos] is smaller or equal than data[parent] swap
+    HEAP_SWAP(heap, pos / 2, pos);
+  }
+}
+
+static
+void lp_polynomial_heap_heapify_down(lp_polynomial_heap_t *heap, size_t pos) {
+  // we're using a 1-based index to enable heap index calculation
+  ++ pos;
+  // while there is still at least a left child
+  while (2 * pos <= heap->size) {
+    size_t l = 2 * pos;
+    size_t r = 2 * pos + 1;
+
+    // find larger child
+    size_t o = (r <= heap->size && HEAP_CMP(heap, l, r) < 0) ? r : l;
+
+    // in case current element is larger than larger child, we're done
+    if (HEAP_CMP(heap, pos, o) >= 0) break;
+
+    // swap larger child with current pos
+    HEAP_SWAP(heap, pos, o);
+    pos = o;
+  }
+}
+
+/** inserts a polynomial, does not copy the polynomial */
+static
+void lp_polynomial_heap_insert(lp_polynomial_heap_t* heap, lp_polynomial_t* p) {
+  assert(p);
+  heap->size++;
+  if (heap->size > heap->data_size) {
+    lp_polynomial_heap_extend(heap);
+  }
+  heap->data[heap->size - 1] = p;
+  lp_polynomial_heap_heapify_up(heap);
+}
+
+void lp_polynomial_heap_push(lp_polynomial_heap_t* heap, const lp_polynomial_t* p) {
+  lp_polynomial_heap_insert(heap, lp_polynomial_new_copy(p));
+}
+
+void lp_polynomial_heap_push_move(lp_polynomial_heap_t* heap, lp_polynomial_t* p) {
+  lp_polynomial_t *tmp = lp_polynomial_new(lp_polynomial_get_context(p));
+  lp_polynomial_swap(tmp, p);
+  lp_polynomial_heap_insert(heap, tmp);
+}
+
+void lp_polynomial_heap_push_vector(lp_polynomial_heap_t* heap, const lp_polynomial_vector_t* v) {
+  size_t size = lp_polynomial_vector_size(v);
+  for (size_t i = 0; i < size; ++i) {
+    // lp_polynomial_vector_at gives a fresh copied polynomial, we reuse it
+    lp_polynomial_heap_insert(heap, lp_polynomial_vector_at(v, i));
+  }
+}
+
+lp_polynomial_t* lp_polynomial_heap_pop(lp_polynomial_heap_t* heap) {
+  // empty heap does not have a top
+  if (heap->size == 0) {
+    return NULL;
+  }
+  // get the top element
+  lp_polynomial_t *ret = heap->data[0];
+  // reduce size and moves last element to top
+  heap->data[0] = heap->data[--heap->size];
+  // restore heap condition
+  lp_polynomial_heap_heapify_down(heap, 0);
+
+  return ret;
+}
+
+int lp_polynomial_heap_remove(lp_polynomial_heap_t* heap, const lp_polynomial_t *p){
+  int result = 0;
+  for (size_t i = 0; i < heap->size; ++i) {
+    if (lp_polynomial_eq(p, heap->data[i])) {
+      heap->data[i] = heap->data[--heap->size];
+      lp_polynomial_heap_heapify_down(heap, i);
+      result++;
+    }
+  }
+  return result;
+}
+
+const lp_polynomial_t* lp_polynomial_heap_peek(const lp_polynomial_heap_t* heap) {
+  // empty heap does not have a top
+  if (heap->size == 0) {
+    return NULL;
+  }
+  // return the top element
+  return heap->data[0];
+}
+
+const lp_polynomial_t* lp_polynomial_heap_at(const lp_polynomial_heap_t* heap, size_t n) {
+  if (n >= heap->size) {
+    return NULL;
+  }
+  return heap->data[n];
+}
+
+void lp_polynomial_heap_clear(lp_polynomial_heap_t* heap) {
+  lp_polynomial_heap_compare_f cmp = heap->cmp;
+  lp_polynomial_heap_destruct(heap);
+  lp_polynomial_heap_construct(heap, cmp);
+}
+
+void lp_polynomial_heap_print(lp_polynomial_heap_t* heap, FILE *stream) {
+  fputc('[', stream);
+  for (size_t i = 0; i < heap->size; ++i) {
+    lp_polynomial_print(heap->data[i], stream);
+    if (i < heap->size - 1)
+      fputc(',', stream);
+  }
+  fputc(']', stream);
+}
\ No newline at end of file
diff --git a/src/polynomial/polynomial_vector.c b/src/polynomial/polynomial_vector.c
index 0de1ffb7..802bd184 100644
--- a/src/polynomial/polynomial_vector.c
+++ b/src/polynomial/polynomial_vector.c
@@ -19,6 +19,7 @@
 
 #include "polynomial_vector.h"
 #include "polynomial.h"
+#include "output.h"
 
 #include "polynomial/gcd.h"
 
@@ -32,6 +33,19 @@ lp_polynomial_vector_t* lp_polynomial_vector_new(const lp_polynomial_context_t*
   return result;
 }
 
+lp_polynomial_vector_t* lp_polynomial_vector_copy(const lp_polynomial_vector_t *v) {
+  lp_polynomial_vector_t* result = (lp_polynomial_vector_t*) malloc(sizeof(lp_polynomial_vector_t));
+  // copy all fields
+  *result = *v;
+  // deep copy of data
+  result->data = (coefficient_t*) malloc(v->capacity * sizeof(coefficient_t));
+  for (size_t i = 0; i < v->size; ++i) {
+    coefficient_construct_copy(result->ctx, result->data + i, v->data + i);
+  }
+  lp_polynomial_context_attach((lp_polynomial_context_t*)result->ctx);
+  return result;
+}
+
 void lp_polynomial_vector_delete(lp_polynomial_vector_t* v) {
   lp_polynomial_vector_destruct(v);
   free(v);
@@ -51,6 +65,14 @@ void lp_polynomial_vector_destruct(lp_polynomial_vector_t* v) {
   lp_polynomial_context_detach((lp_polynomial_context_t*)v->ctx);
 }
 
+const lp_polynomial_context_t* lp_polynomial_vector_get_context(const lp_polynomial_vector_t *v) {
+    return v->ctx;
+}
+
+void lp_polynomial_vector_swap(lp_polynomial_vector_t *v1, lp_polynomial_vector_t *v2) {
+    lp_polynomial_vector_t tmp = *v1; *v1 = *v2; *v2 = tmp;
+}
+
 static inline
 void lp_polynomial_vector_check_size_for_add(lp_polynomial_vector_t* v) {
   if (v->size == v->capacity) {
@@ -66,6 +88,15 @@ void lp_polynomial_vector_push_back(lp_polynomial_vector_t* v, const lp_polynomi
   v->size ++;
 }
 
+void lp_polynomial_vector_push_back_move(lp_polynomial_vector_t* v, lp_polynomial_t* p) {
+  assert(lp_polynomial_context_equal(v->ctx, p->ctx));
+  lp_polynomial_vector_check_size_for_add(v);
+  coefficient_t *v_p = v->data + v->size;
+  coefficient_construct_from_int(v->ctx, v_p, 0);
+  coefficient_swap(&p->data, v_p);
+  v->size ++;
+}
+
 void lp_polynomial_vector_push_back_coeff(lp_polynomial_vector_t* v, const coefficient_t* C) {
   lp_polynomial_vector_check_size_for_add(v);
   coefficient_construct_copy(v->ctx, v->data + v->size, C);
@@ -73,7 +104,7 @@ void lp_polynomial_vector_push_back_coeff(lp_polynomial_vector_t* v, const coeff
 }
 
 void lp_polynomial_vector_reset(lp_polynomial_vector_t* v) {
-  size_t i = 0;
+  size_t i;
   for (i = 0; i < v->size; ++ i) {
     coefficient_destruct(v->data + i);
   }
@@ -88,6 +119,17 @@ lp_polynomial_t* lp_polynomial_vector_at(const lp_polynomial_vector_t* v, size_t
   return lp_polynomial_new_from_coefficient(v->ctx, v->data + i);
 }
 
+void lp_polynomial_vector_print(const lp_polynomial_vector_t* v, FILE* out) {
+    fputc('[', out);
+    size_t i;
+    for (i = 0; i < v->size; ++ i) {
+        coefficient_print(v->ctx, v->data + i, out);
+        if (i != v->size - 1)
+            fputc(',', out);
+    }
+    fputc(']', out);
+}
+
 void lp_polynomial_vector_push_back_coeff_prime(lp_polynomial_vector_t* v, const coefficient_t* C) {
 
   size_t old_size = v->size;
diff --git a/src/upolynomial/upolynomial.c b/src/upolynomial/upolynomial.c
index 57a67818..cbf048ed 100644
--- a/src/upolynomial/upolynomial.c
+++ b/src/upolynomial/upolynomial.c
@@ -31,11 +31,11 @@
 
 #include <stdlib.h>
 #include <assert.h>
-#include <string.h>
 
-#include <upolynomial.h>
 #include "upolynomial/upolynomial.h"
 
+#include <polynomial.h>
+
 size_t lp_upolynomial_degree(const lp_upolynomial_t* p) {
   assert(p);
   assert(p->size > 0);
@@ -43,6 +43,7 @@ size_t lp_upolynomial_degree(const lp_upolynomial_t* p) {
 }
 
 // Construct the polynomial, but don't construct monomials
+static
 lp_upolynomial_t* lp_upolynomial_construct_empty(const lp_int_ring_t* K, size_t size) {
   size_t malloc_size = sizeof(lp_upolynomial_t) + size*sizeof(ulp_monomial_t);
   lp_upolynomial_t* new_p = (lp_upolynomial_t*) malloc(malloc_size);
@@ -1206,3 +1207,22 @@ lp_upolynomial_t* lp_upolynomial_subst_x_neg(const lp_upolynomial_t* f) {
 
   return neg;
 }
+
+lp_polynomial_t* lp_upolynomial_to_polynomial(const lp_upolynomial_t* p, const lp_polynomial_context_t* ctx, lp_variable_t var) {
+  assert(p);
+  assert(lp_int_ring_equal(p->K, ctx->K));
+
+  lp_polynomial_t *poly = lp_polynomial_new(ctx);
+
+  lp_monomial_t mono;
+  lp_monomial_construct(ctx, &mono);
+  for (size_t i = 0; i < p->size; ++i) {
+    lp_monomial_set_coefficient(ctx, &mono, &p->monomials[i].coefficient);
+    lp_monomial_push(&mono, var, p->monomials[i].degree);
+    lp_polynomial_add_monomial(poly, &mono);
+    lp_monomial_clear(ctx, &mono);
+  }
+  lp_monomial_destruct(&mono);
+
+  return poly;
+}