diff --git a/include/polyxx/context.h b/include/polyxx/context.h
index a1b21106..b378c004 100644
--- a/include/polyxx/context.h
+++ b/include/polyxx/context.h
@@ -19,7 +19,18 @@ namespace poly {
     deleting_unique_ptr<lp_polynomial_context_t> mPolynomialContext;
 
    public:
+    /** Constructs a new (empty) context */
     Context();
+
+    /** Wraps the lp context. */
+    explicit Context(lp_polynomial_context_t* ctx);
+
+    /** Copy constructor. */
+    Context(const Context& other);
+
+    /** No assignment as this would mess with internal reference counting. */
+    Context& operator=(const Context& other) = delete;
+
     /** Get a non-const pointer to the internal lp_variable_db_t.
      * Handle with care!
      */
diff --git a/include/polyxx/dyadic_interval.h b/include/polyxx/dyadic_interval.h
index 6836e048..08961ee6 100644
--- a/include/polyxx/dyadic_interval.h
+++ b/include/polyxx/dyadic_interval.h
@@ -33,12 +33,14 @@ namespace poly {
     /** Construct an open interval. */
     DyadicInterval(const Integer& a, const Integer& b);
     /** Construct an interval from the given bounds. */
-    DyadicInterval(const Integer& a, bool a_open, const Integer& b,
-                   bool b_open);
+    DyadicInterval(const Integer& a, bool a_open,
+                   const Integer& b, bool b_open);
+    /** Construct a point interval. */
+    explicit DyadicInterval(long i);
     /** Construct an open interval. */
-    DyadicInterval(long a, long b);
+    explicit DyadicInterval(long a, long b);
     /** Construct an interval from the given bounds. */
-    DyadicInterval(long a, bool a_open, long b, bool b_open);
+    explicit DyadicInterval(long a, bool a_open, long b, bool b_open);
     /** Copy from a DyadicInterval. */
     DyadicInterval(const DyadicInterval& i);
     /** Move from a DyadicInterval. */
@@ -58,11 +60,11 @@ namespace poly {
 
     /** Collapse this interval to a single point. */
     void collapse(const DyadicRational& dr);
-    /** The the lower bound. */
+    /** The lower bound. */
     void set_lower(const DyadicRational& dr, bool open);
-    /** The the upper bound. */
+    /** The upper bound. */
     void set_upper(const DyadicRational& dr, bool open);
-    /** The this interval by 2^n. */
+    /** The interval by 2^n. */
     void scale(int n);
   };
 
diff --git a/include/polyxx/dyadic_rational.h b/include/polyxx/dyadic_rational.h
index a8be408d..8ef68a29 100644
--- a/include/polyxx/dyadic_rational.h
+++ b/include/polyxx/dyadic_rational.h
@@ -31,6 +31,8 @@ namespace poly {
     explicit DyadicRational(double d);
     /** Construct from an int. */
     DyadicRational(int i);
+    /** Construct from a long. */
+    DyadicRational(long i);
 
     /** Construct from an internal lp_dyadic_rational_t pointer. */
     explicit DyadicRational(const lp_dyadic_rational_t* dr);
diff --git a/include/polyxx/integer.h b/include/polyxx/integer.h
index b4e22d9d..a3f8c3e4 100644
--- a/include/polyxx/integer.h
+++ b/include/polyxx/integer.h
@@ -22,9 +22,9 @@ namespace poly {
     /** Constructs zero. */
     Integer();
     /** Constructs from an int. */
-    explicit Integer(int i);
+    Integer(int i);
     /** Constructs from a long. */
-    explicit Integer(long i);
+    Integer(long i);
     /** Constructs from a long into the given ring. */
     Integer(const IntegerRing& ir, long i);
     /** Constructs from a string. */
@@ -122,16 +122,28 @@ namespace poly {
 
   /** Compare two integers. */
   bool operator==(const Integer& lhs, const Integer& rhs);
+  bool operator==(const Integer& lhs, long rhs);
+  bool operator==(long lhs, const Integer& rhs);
   /** Compare two integers. */
   bool operator!=(const Integer& lhs, const Integer& rhs);
+  bool operator!=(const Integer& lhs, long rhs);
+  bool operator!=(long lhs, const Integer& rhs);
   /** Compare two integers according to the lexicographic ordering on (lower bound,upper bound). */
   bool operator<(const Integer& lhs, const Integer& rhs);
+  bool operator<(const Integer& lhs, long rhs);
+  bool operator<(long lhs, const Integer& rhs);
   /** Compare two integers according to the lexicographic ordering on (lower bound,upper bound). */
   bool operator<=(const Integer& lhs, const Integer& rhs);
+  bool operator<=(const Integer& lhs, long rhs);
+  bool operator<=(long lhs, const Integer& rhs);
   /** Compare two integers according to the lexicographic ordering on (lower bound,upper bound). */
   bool operator>(const Integer& lhs, const Integer& rhs);
+  bool operator>(const Integer& lhs, long rhs);
+  bool operator>(long lhs, const Integer& rhs);
   /** Compare two integers according to the lexicographic ordering on (lower bound,upper bound). */
   bool operator>=(const Integer& lhs, const Integer& rhs);
+  bool operator>=(const Integer& lhs, long rhs);
+  bool operator>=(long lhs, const Integer& rhs);
 
   /** Compare two integers over the given ring. */
   int compare(const IntegerRing& ir, const Integer& lhs, const Integer& rhs);
diff --git a/include/polyxx/polynomial.h b/include/polyxx/polynomial.h
index fd4f8e68..9e2ada82 100644
--- a/include/polyxx/polynomial.h
+++ b/include/polyxx/polynomial.h
@@ -33,26 +33,27 @@ namespace poly {
     /** Construct a zero polynomial. */
     Polynomial();
 
+    /** Construct a constant polynomial from an internal context pointer. */
+    Polynomial(const lp_polynomial_context_t* c, const Variable& v);
     /** Construct from a variable and a custom context. */
-    Polynomial(const Context& c, Variable v);
+    Polynomial(const Context& c, const Variable& v);
     /** Construct from a variable. */
-    Polynomial(Variable v);
+    Polynomial(const Variable& v);
 
+    /** Construct a constant polynomial from an internal context pointer. */
+    Polynomial(const lp_polynomial_context_t* c, const Integer& i, const Variable& v, unsigned n);
     /** Construct i * v^n from a custom context. */
-    Polynomial(const Context& c, const Integer& i, Variable v, unsigned n);
+    Polynomial(const Context& c, const Integer& i, const Variable& v, unsigned n);
     /** Construct i * v^n. */
-    Polynomial(const Integer& i, Variable v, unsigned n);
+    Polynomial(const Integer& i, const Variable& v, unsigned n);
 
+    /** Construct a constant polynomial from an internal context pointer. */
+    Polynomial(const lp_polynomial_context_t* c, const Integer& i);
     /** Construct from an integer and a custom context. */
     Polynomial(const Context& c, const Integer& i);
     /** Construct from an integer. */
     Polynomial(const Integer& i);
 
-    /** Construct from an integer and a custom context. */
-    Polynomial(const Context& c, long i);
-    /** Construct from an integer. */
-    Polynomial(long i);
-
     /** Copy from a Polynomial. */
     Polynomial(const Polynomial& p);
     /** Move from a Polynomial. */
@@ -132,16 +133,28 @@ namespace poly {
 
   /** Compare polynomials. */
   bool operator==(const Polynomial& lhs, const Polynomial& rhs);
+  bool operator==(const Integer& lhs, const Polynomial& rhs);
+  bool operator==(const Polynomial& lhs, const Integer& rhs);
   /** Compare polynomials. */
   bool operator!=(const Polynomial& lhs, const Polynomial& rhs);
+  bool operator!=(const Integer& lhs, const Polynomial& rhs);
+  bool operator!=(const Polynomial& lhs, const Integer& rhs);
   /** Compare polynomials. */
   bool operator<(const Polynomial& lhs, const Polynomial& rhs);
+  bool operator<(const Integer& lhs, const Polynomial& rhs);
+  bool operator<(const Polynomial& lhs, const Integer& rhs);
   /** Compare polynomials. */
   bool operator<=(const Polynomial& lhs, const Polynomial& rhs);
+  bool operator<=(const Integer& lhs, const Polynomial& rhs);
+  bool operator<=(const Polynomial& lhs, const Integer& rhs);
   /** Compare polynomials. */
   bool operator>(const Polynomial& lhs, const Polynomial& rhs);
+  bool operator>(const Integer& lhs, const Polynomial& rhs);
+  bool operator>(const Polynomial& lhs, const Integer& rhs);
   /** Compare polynomials. */
   bool operator>=(const Polynomial& lhs, const Polynomial& rhs);
+  bool operator>=(const Integer& lhs, const Polynomial& rhs);
+  bool operator>=(const Polynomial& lhs, const Integer& rhs);
 
   /** Add two polynomials. */
   Polynomial operator+(const Polynomial& lhs, const Polynomial& rhs);
@@ -151,6 +164,8 @@ namespace poly {
   Polynomial operator+(const Integer& lhs, const Polynomial& rhs);
   /** Add and assign two polynomials. */
   Polynomial& operator+=(Polynomial& lhs, const Polynomial& rhs);
+  /** Add and assign a polynomial with an integer. */
+  Polynomial& operator+=(Polynomial& lhs, const Integer& rhs);
   /** Compute lhs += rhs1 * rhs2. */
   Polynomial& add_mul(Polynomial& lhs, const Polynomial& rhs1, const Polynomial& rhs2);
 
@@ -164,6 +179,8 @@ namespace poly {
   Polynomial operator-(const Integer& lhs, const Polynomial& rhs);
   /** Subtract and assign two polynomials. */
   Polynomial& operator-=(Polynomial& lhs, const Polynomial& rhs);
+  /** Subtract and assigns a polynomial with an integer. */
+  Polynomial& operator-=(Polynomial& lhs, const Integer& rhs);
   /** Compute lhs -= rhs1 * rhs2. */
   Polynomial& sub_mul(Polynomial& lhs, const Polynomial& rhs1, const Polynomial& rhs2);
 
@@ -175,6 +192,8 @@ namespace poly {
   Polynomial operator*(const Integer& lhs, const Polynomial& rhs);
   /** Multiply and assign two polynomials. */
   Polynomial& operator*=(Polynomial& lhs, const Polynomial& rhs);
+  /** Multiply and assign a polynomial with an integer */
+  Polynomial& operator*=(Polynomial& lhs, const Integer& rhs);
 
   /** Multiply with x^n where x is the main variable. */
   Polynomial shl(const Polynomial& lhs, unsigned n);
diff --git a/include/polyxx/variable.h b/include/polyxx/variable.h
index 614994ce..9c2e97df 100644
--- a/include/polyxx/variable.h
+++ b/include/polyxx/variable.h
@@ -9,6 +9,8 @@ namespace poly {
 
   /**
    * Implements a wrapper for lp_variable_t.
+   * WARNING: Variable does not store its context, i.e. the variable database
+   * it belongs to. Only use it with a Polynomial / Variable of the same context!
    */
   class Variable {
     /** The actual variable. */
@@ -18,11 +20,11 @@ namespace poly {
     /** Construct with a null variable. */
     Variable();
     /** Construct from a lp_variable_t. */
-    Variable(lp_variable_t var);
+    explicit Variable(lp_variable_t var);
     /** Construct a new variable with the given name in the specified context. */
-    Variable(const Context& c, const char* name);
+    explicit Variable(const Context& c, const char* name);
     /** Construct a new variable with the given name in the default context. */
-    Variable(const char* name);
+    explicit Variable(const char* name);
 
     /** Get the internal lp_variable_t. Note that it's only a type alias for
      * long.
diff --git a/include/variable_db.h b/include/variable_db.h
index 57dff087..dd2f026d 100644
--- a/include/variable_db.h
+++ b/include/variable_db.h
@@ -46,6 +46,9 @@ int lp_variable_db_print(const lp_variable_db_t* var_db, FILE* out);
 /** Get the name of the variable */
 const char* lp_variable_db_get_name(const lp_variable_db_t* var_db, lp_variable_t var);
 
+/** Checks if var is a valid index */
+int lp_variable_db_is_valid(const lp_variable_db_t* var_db, lp_variable_t var);
+
 #ifdef __cplusplus
 } /* close extern "C" { */
 #endif
diff --git a/src/polyxx/context.cpp b/src/polyxx/context.cpp
index bad8945a..2c65dd59 100644
--- a/src/polyxx/context.cpp
+++ b/src/polyxx/context.cpp
@@ -16,6 +16,25 @@ namespace poly {
           lp_polynomial_context_detach(ptr);
         });
   }
+
+  Context::Context(lp_polynomial_context_t* ctx) {
+    mVariableDB = deleting_unique_ptr<lp_variable_db_t>(
+      ctx->var_db,
+      [](lp_variable_db_t *ptr) { lp_variable_db_detach(ptr); });
+    mVariableOrder = deleting_unique_ptr<lp_variable_order_t>(
+      ctx->var_order,
+      [](lp_variable_order_t *ptr) { lp_variable_order_detach(ptr); });
+    mPolynomialContext = deleting_unique_ptr<lp_polynomial_context_t>(
+      ctx,
+      [](lp_polynomial_context_t *ptr) { lp_polynomial_context_detach(ptr); });
+
+    lp_variable_db_attach(ctx->var_db);
+    lp_variable_order_attach(ctx->var_order);
+    lp_polynomial_context_attach(ctx);
+  }
+
+  Context::Context(const Context& other) : Context(other.get_polynomial_context()) {}
+
   lp_variable_db_t* Context::get_variable_db() const {
     return const_cast<lp_variable_db_t*>(mVariableDB.get());
   }
diff --git a/src/polyxx/dyadic_interval.cpp b/src/polyxx/dyadic_interval.cpp
index c2d27a4a..f3b06428 100644
--- a/src/polyxx/dyadic_interval.cpp
+++ b/src/polyxx/dyadic_interval.cpp
@@ -27,6 +27,8 @@ namespace poly {
     lp_dyadic_interval_construct_from_integer(get_internal(), a.get_internal(),
                                               a_open, b.get_internal(), b_open);
   }
+  DyadicInterval::DyadicInterval(long i)
+      : DyadicInterval(DyadicRational(i)) {}
   DyadicInterval::DyadicInterval(long a, long b)
       : DyadicInterval(a, true, b, true) {}
   DyadicInterval::DyadicInterval(long a, bool a_open, long b, bool b_open) {
diff --git a/src/polyxx/dyadic_rational.cpp b/src/polyxx/dyadic_rational.cpp
index c80273a3..e77873b8 100644
--- a/src/polyxx/dyadic_rational.cpp
+++ b/src/polyxx/dyadic_rational.cpp
@@ -25,6 +25,7 @@ namespace poly {
     lp_dyadic_rational_construct_from_double(&mDRat, d);
   }
   DyadicRational::DyadicRational(int i) : DyadicRational(i, 0) {}
+  DyadicRational::DyadicRational(long i) : DyadicRational(i, 0) {}
 
   DyadicRational::DyadicRational(const lp_dyadic_rational_t* dr) {
     lp_dyadic_rational_construct_copy(&mDRat, dr);
diff --git a/src/polyxx/integer.cpp b/src/polyxx/integer.cpp
index 5067a2c5..ec33f561 100644
--- a/src/polyxx/integer.cpp
+++ b/src/polyxx/integer.cpp
@@ -80,21 +80,57 @@ namespace poly {
   bool operator==(const Integer& lhs, const Integer& rhs) {
     return compare(IntegerRing::Z, lhs, rhs) == 0;
   }
+  bool operator==(const Integer& lhs, long rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) == 0;
+  }
+  bool operator==(long lhs, const Integer& rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) == 0;
+  }
   bool operator!=(const Integer& lhs, const Integer& rhs) {
     return compare(IntegerRing::Z, lhs, rhs) != 0;
   }
+  bool operator!=(const Integer& lhs, long rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) != 0;
+  }
+  bool operator!=(long lhs, const Integer& rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) != 0;
+  }
   bool operator<(const Integer& lhs, const Integer& rhs) {
     return compare(IntegerRing::Z, lhs, rhs) < 0;
   }
+  bool operator<(const Integer& lhs, long rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) < 0;
+  }
+  bool operator<(long lhs, const Integer& rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) < 0;
+  }
   bool operator<=(const Integer& lhs, const Integer& rhs) {
     return compare(IntegerRing::Z, lhs, rhs) <= 0;
   }
+  bool operator<=(const Integer& lhs, long rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) <= 0;
+  }
+  bool operator<=(long lhs, const Integer& rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) <= 0;
+  }
   bool operator>(const Integer& lhs, const Integer& rhs) {
     return compare(IntegerRing::Z, lhs, rhs) > 0;
   }
+  bool operator>(const Integer& lhs, long rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) > 0;
+  }
+  bool operator>(long lhs, const Integer& rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) > 0;
+  }
   bool operator>=(const Integer& lhs, const Integer& rhs) {
     return compare(IntegerRing::Z, lhs, rhs) >= 0;
   }
+  bool operator>=(const Integer& lhs, long rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) >= 0;
+  }
+  bool operator>=(long lhs, const Integer& rhs) {
+    return compare(IntegerRing::Z, lhs, rhs) >= 0;
+  }
 
   int compare(const IntegerRing& ir, const Integer& lhs, const Integer& rhs) {
     return lp_integer_cmp(ir.get_internal(), lhs.get_internal(),
diff --git a/src/polyxx/polynomial.cpp b/src/polyxx/polynomial.cpp
index 87f33fd0..1f40fb62 100644
--- a/src/polyxx/polynomial.cpp
+++ b/src/polyxx/polynomial.cpp
@@ -36,25 +36,35 @@ namespace poly {
       : Polynomial(c.get_polynomial_context()) {}
   Polynomial::Polynomial() : Polynomial(Context::get_context()) {}
 
-  Polynomial::Polynomial(const Context& c, Variable v)
-      : Polynomial(c, Integer(1), v, 1) {}
-  Polynomial::Polynomial(Variable v) : Polynomial(Context::get_context(), v) {}
-  Polynomial::Polynomial(const Context& c, const Integer &i, Variable v, unsigned n)
+  Polynomial::Polynomial(const lp_polynomial_context_t* c, const Variable& v)
       : mPoly(lp_polynomial_alloc(), polynomial_deleter) {
-    lp_polynomial_construct_simple(get_internal(), c.get_polynomial_context(),
+    assert(lp_variable_db_is_valid(c->var_db, v.get_internal()));
+    lp_polynomial_construct_simple(get_internal(), c,
+                                   Integer(1).get_internal(),
+                                   v.get_internal(), 1);
+  }
+  Polynomial::Polynomial(const Context& c, const Variable& v)
+      : Polynomial(c.get_polynomial_context(), v) {}
+  Polynomial::Polynomial(const Variable& v) : Polynomial(Context::get_context(), v) {}
+
+  Polynomial::Polynomial(const lp_polynomial_context_t* c, const Integer& i, const Variable& v, unsigned n)
+      : mPoly(lp_polynomial_alloc(), polynomial_deleter) {
+    assert(lp_variable_db_is_valid(c->var_db, v.get_internal()));
+    lp_polynomial_construct_simple(get_internal(), c,
                                    i.get_internal(), v.get_internal(), n);
   }
-  Polynomial::Polynomial(const Integer& i, Variable v, unsigned n)
+  Polynomial::Polynomial(const Context& c, const Integer& i, const Variable& v, unsigned n)
+      : Polynomial(c.get_polynomial_context(), i, v, n) {}
+  Polynomial::Polynomial(const Integer& i, const Variable& v, unsigned n)
       : Polynomial(Context::get_context(), i, v, n) {}
-  Polynomial::Polynomial(const Context& c, const Integer& i)
+
+  Polynomial::Polynomial(const lp_polynomial_context_t* c, const Integer & i)
       : mPoly(lp_polynomial_alloc(), polynomial_deleter) {
-    lp_polynomial_construct_simple(get_internal(), c.get_polynomial_context(),
+    lp_polynomial_construct_simple(get_internal(), c,
                                    i.get_internal(), lp_variable_null, 0);
   }
+  Polynomial::Polynomial(const Context& c, const Integer& i) : Polynomial(c.get_polynomial_context(), i) {}
   Polynomial::Polynomial(const Integer& i) : Polynomial(Context::get_context(), i){};
-  Polynomial::Polynomial(const Context& c, long i)
-      : Polynomial(c, Integer(i)) {}
-  Polynomial::Polynomial(long i) : Polynomial(Context::get_context(), i){};
 
   Polynomial::Polynomial(const Polynomial& p)
       : mPoly(lp_polynomial_new_copy(p.get_internal()), polynomial_deleter) {}
@@ -107,7 +117,7 @@ namespace poly {
     return lp_polynomial_degree(p.get_internal());
   }
   Variable main_variable(const Polynomial& p) {
-    return lp_polynomial_top_variable(p.get_internal());
+    return Variable(lp_polynomial_top_variable(p.get_internal()));
   }
   Polynomial coefficient(const Polynomial& p, std::size_t k) {
     Polynomial res(detail::context(p));
@@ -164,21 +174,67 @@ namespace poly {
   bool operator==(const Polynomial& lhs, const Polynomial& rhs) {
     return lp_polynomial_eq(lhs.get_internal(), rhs.get_internal());
   }
+  bool operator==(const Polynomial& lhs, const Integer& rhs) {
+    Polynomial tmp(detail::context(lhs), rhs);
+    return lp_polynomial_eq(lhs.get_internal(), tmp.get_internal());
+  }
+  bool operator==(const Integer& lhs, const Polynomial& rhs) {
+    Polynomial tmp(detail::context(rhs), lhs);
+    return lp_polynomial_eq(tmp.get_internal(), rhs.get_internal());
+  }
   bool operator!=(const Polynomial& lhs, const Polynomial& rhs) {
     return !(lhs == rhs);
   }
+  bool operator!=(const Polynomial& lhs, const Integer& rhs) {
+    return !(lhs == rhs);
+  }
+  bool operator!=(const Integer& lhs, const Polynomial& rhs) {
+    return !(lhs == rhs);
+  }
   bool operator<(const Polynomial& lhs, const Polynomial& rhs) {
     return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) < 0;
   }
+  bool operator<(const Polynomial& lhs, const Integer& rhs) {
+    Polynomial tmp(detail::context(lhs), rhs);
+    return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) < 0;
+  }
+  bool operator<(const Integer& lhs, const Polynomial& rhs) {
+    Polynomial tmp(detail::context(rhs), lhs);
+    return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) < 0;
+  }
   bool operator<=(const Polynomial& lhs, const Polynomial& rhs) {
     return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) <= 0;
   }
+  bool operator<=(const Polynomial& lhs, const Integer& rhs) {
+    Polynomial tmp(detail::context(lhs), rhs);
+    return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) <= 0;
+  }
+  bool operator<=(const Integer& lhs, const Polynomial& rhs) {
+    Polynomial tmp(detail::context(rhs), lhs);
+    return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) <= 0;
+  }
   bool operator>(const Polynomial& lhs, const Polynomial& rhs) {
     return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) > 0;
   }
+  bool operator>(const Polynomial& lhs, const Integer& rhs) {
+    Polynomial tmp(detail::context(lhs), rhs);
+    return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) > 0;
+  }
+  bool operator>(const Integer& lhs, const Polynomial& rhs) {
+    Polynomial tmp(detail::context(rhs), lhs);
+    return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) > 0;
+  }
   bool operator>=(const Polynomial& lhs, const Polynomial& rhs) {
     return lp_polynomial_cmp(lhs.get_internal(), rhs.get_internal()) >= 0;
   }
+  bool operator>=(const Polynomial& lhs, const Integer& rhs) {
+    Polynomial tmp(detail::context(lhs), rhs);
+    return lp_polynomial_cmp(lhs.get_internal(), tmp.get_internal()) >= 0;
+  }
+  bool operator>=(const Integer& lhs, const Polynomial& rhs) {
+    Polynomial tmp(detail::context(rhs), lhs);
+    return lp_polynomial_cmp(tmp.get_internal(), rhs.get_internal()) >= 0;
+  }
 
   Polynomial operator+(const Polynomial& lhs, const Polynomial& rhs) {
     Polynomial res(detail::context(lhs, rhs));
@@ -187,23 +243,30 @@ namespace poly {
     return res;
   }
   Polynomial operator+(const Polynomial& lhs, const Integer& rhs) {
-    lp_monomial_t monomial;
-    lp_monomial_construct(detail::context(lhs), &monomial);
-    lp_monomial_set_coefficient(detail::context(lhs), &monomial,
-                                rhs.get_internal());
     Polynomial res(lhs);
-    lp_polynomial_add_monomial(res.get_internal(), &monomial);
-    lp_monomial_destruct(&monomial);
+    res += rhs;
     return res;
   }
   Polynomial operator+(const Integer& lhs, const Polynomial& rhs) {
-    return rhs + lhs;
+    Polynomial res(rhs);
+    res += lhs;
+    return res;
   }
   Polynomial& operator+=(Polynomial& lhs, const Polynomial& rhs) {
     lp_polynomial_add(lhs.get_internal(), lhs.get_internal(),
                       rhs.get_internal());
     return lhs;
   }
+  Polynomial& operator+=(Polynomial& lhs, const Integer& rhs) {
+    lp_monomial_t monomial;
+    lp_monomial_construct(detail::context(lhs), &monomial);
+    lp_monomial_set_coefficient(detail::context(lhs), &monomial,
+                                rhs.get_internal());
+    lp_polynomial_add_monomial(lhs.get_internal(), &monomial);
+    lp_monomial_destruct(&monomial);
+    return lhs;
+  }
+
   Polynomial& add_mul(Polynomial& lhs, const Polynomial& rhs1,
                       const Polynomial& rhs2) {
     lp_polynomial_add_mul(lhs.get_internal(), rhs1.get_internal(),
@@ -226,13 +289,18 @@ namespace poly {
     return lhs + (-rhs);
   }
   Polynomial operator-(const Integer& lhs, const Polynomial& rhs) {
-    return -rhs + lhs;
+    return (-rhs) + lhs;
   }
   Polynomial& operator-=(Polynomial& lhs, const Polynomial& rhs) {
     lp_polynomial_sub(lhs.get_internal(), lhs.get_internal(),
                       rhs.get_internal());
     return lhs;
   }
+  Polynomial& operator-=(Polynomial& lhs, const Integer& rhs) {
+    lhs += (-rhs);
+    return lhs;
+  }
+
   Polynomial& sub_mul(Polynomial& lhs, const Polynomial& rhs1,
                       const Polynomial& rhs2) {
     lp_polynomial_sub_mul(lhs.get_internal(), rhs1.get_internal(),
@@ -260,6 +328,11 @@ namespace poly {
                       rhs.get_internal());
     return lhs;
   }
+  Polynomial& operator*=(Polynomial& lhs, const Integer& rhs) {
+    lp_polynomial_mul_integer(lhs.get_internal(), lhs.get_internal(),
+                              rhs.get_internal());
+    return lhs;
+  }
 
   Polynomial shl(const Polynomial& lhs, unsigned exp) {
     Polynomial res(detail::context(lhs));
@@ -368,10 +441,9 @@ namespace poly {
   }
   Polynomial discriminant(const Polynomial& p) {
     if (degree(p) == 1) {
-      // Derivative is constant, making the resultant trivial (and resultant()
-      // does not cope with that)
-      // This creates the constant polynomial 1 in the context of p.
-      return Polynomial(detail::context(p)) + Integer(1);
+      // Derivative is constant, making the resultant trivial
+      // (and resultant() does not cope with that)
+      return Polynomial(detail::context(p), Integer(1));
     }
     return div(resultant(p, derivative(p)), leading_coefficient(p));
   }
diff --git a/src/polyxx/rational_interval.cpp b/src/polyxx/rational_interval.cpp
index 9983e69a..d365d919 100644
--- a/src/polyxx/rational_interval.cpp
+++ b/src/polyxx/rational_interval.cpp
@@ -1,6 +1,7 @@
 #include "polyxx/rational_interval.h"
 
 #include <cassert>
+#include <iostream>
 
 namespace poly {
 
diff --git a/src/variable/variable_db.c b/src/variable/variable_db.c
index 60f3df4c..eeb82327 100644
--- a/src/variable/variable_db.c
+++ b/src/variable/variable_db.c
@@ -127,3 +127,7 @@ const char* lp_variable_db_get_name(const lp_variable_db_t* var_db, lp_variable_
   assert(var < var_db->size);
   return var_db->variable_names[var];
 }
+
+int lp_variable_db_is_valid(const lp_variable_db_t* var_db, lp_variable_t var) {
+  return var < var_db->size;
+}
diff --git a/test/polyxx/test_context.cpp b/test/polyxx/test_context.cpp
index 112aa08b..c3538534 100644
--- a/test/polyxx/test_context.cpp
+++ b/test/polyxx/test_context.cpp
@@ -9,6 +9,16 @@ TEST_CASE("context::discriminant") {
   Context ctx;
   Variable x(ctx, "x");
   Polynomial p1(ctx, x);
-  Polynomial d = poly::discriminant(p1); // d built in the default context
-  Polynomial p2 = p1 + d; // assertion failure!
+  Polynomial d = discriminant(p1);
+  Polynomial p2 = p1 + d;
+  CHECK(1 == d);
+  CHECK(p2 == Polynomial(ctx, x) + 1);
+}
+
+TEST_CASE("context::integer") {
+  Context ctx;
+  Polynomial p(ctx);
+  Polynomial p1 = p + 1;
+  Polynomial p2 = p + Integer(1);
+  CHECK(p1 == p2);
 }