From a3b97e1d16f92ea947df0a2f01e1070bed5feb5c Mon Sep 17 00:00:00 2001
From: Ben Ruijl <benruyl@gmail.com>
Date: Sun, 24 Nov 2024 10:38:38 +0100
Subject: [PATCH] Allow negative exponents in polynomials

- Improve speed of polynomial multiplication
- Support LaTeX output for series
---
 src/api/python.rs                             |  69 ++
 src/domains/algebraic_number.rs               |   8 +-
 src/domains/factorized_rational_polynomial.rs |  54 +-
 src/domains/rational_polynomial.rs            |  67 +-
 src/parser.rs                                 |   4 +-
 src/poly.rs                                   | 559 ++++++++--
 src/poly/evaluate.rs                          |  10 +-
 src/poly/factor.rs                            |  24 +-
 src/poly/gcd.rs                               |  31 +-
 src/poly/polynomial.rs                        | 997 ++++++++++--------
 src/poly/series.rs                            |  78 +-
 src/poly/univariate.rs                        |  17 +-
 src/solve.rs                                  |   6 +-
 symbolica.pyi                                 |  26 +
 14 files changed, 1296 insertions(+), 654 deletions(-)

diff --git a/src/api/python.rs b/src/api/python.rs
index 829bb76..36bceb7 100644
--- a/src/api/python.rs
+++ b/src/api/python.rs
@@ -4941,6 +4941,75 @@ impl PythonSeries {
         Ok(format!("{}", self.series))
     }
 
+    /// Convert the series into a LaTeX string.
+    pub fn to_latex(&self) -> PyResult<String> {
+        Ok(format!(
+            "$${}$$",
+            self.series
+                .format_string(&PrintOptions::latex(), PrintState::new())
+        ))
+    }
+
+    /// Convert the expression into a human-readable string, with tunable settings.
+    ///
+    /// Examples
+    /// --------
+    /// >>> a = Expression.parse('128378127123 z^(2/3)*w^2/x/y + y^4 + z^34 + x^(x+2)+3/5+f(x,x^2)')
+    /// >>> print(a.format(number_thousands_separator='_', multiplication_operator=' '))
+    #[pyo3(signature =
+        (terms_on_new_line = false,
+            color_top_level_sum = true,
+            color_builtin_symbols = true,
+            print_finite_field = true,
+            symmetric_representation_for_finite_field = false,
+            explicit_rational_polynomial = false,
+            number_thousands_separator = None,
+            multiplication_operator = '*',
+            double_star_for_exponentiation = false,
+            square_brackets_for_function = false,
+            num_exp_as_superscript = true,
+            latex = false,
+            precision = None)
+        )]
+    pub fn format(
+        &self,
+        terms_on_new_line: bool,
+        color_top_level_sum: bool,
+        color_builtin_symbols: bool,
+        print_finite_field: bool,
+        symmetric_representation_for_finite_field: bool,
+        explicit_rational_polynomial: bool,
+        number_thousands_separator: Option<char>,
+        multiplication_operator: char,
+        double_star_for_exponentiation: bool,
+        square_brackets_for_function: bool,
+        num_exp_as_superscript: bool,
+        latex: bool,
+        precision: Option<usize>,
+    ) -> PyResult<String> {
+        Ok(format!(
+            "{}",
+            self.series.format_string(
+                &PrintOptions {
+                    terms_on_new_line,
+                    color_top_level_sum,
+                    color_builtin_symbols,
+                    print_finite_field,
+                    symmetric_representation_for_finite_field,
+                    explicit_rational_polynomial,
+                    number_thousands_separator,
+                    multiplication_operator,
+                    double_star_for_exponentiation,
+                    square_brackets_for_function,
+                    num_exp_as_superscript,
+                    latex,
+                    precision,
+                },
+                PrintState::new()
+            )
+        ))
+    }
+
     pub fn sin(&self) -> PyResult<Self> {
         Ok(Self {
             series: self
diff --git a/src/domains/algebraic_number.rs b/src/domains/algebraic_number.rs
index 00e913d..a13116c 100644
--- a/src/domains/algebraic_number.rs
+++ b/src/domains/algebraic_number.rs
@@ -6,8 +6,8 @@ use crate::{
     coefficient::ConvertToRing,
     combinatorics::CombinationIterator,
     poly::{
-        factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial, Exponent,
-        Variable,
+        factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial,
+        PositiveExponent, Variable,
     },
 };
 
@@ -577,7 +577,9 @@ impl<R: Field + PolynomialGCD<u16>> AlgebraicExtension<R> {
     }
 }
 
-impl<R: Field + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<AlgebraicExtension<R>, E> {
+impl<R: Field + PolynomialGCD<E>, E: PositiveExponent>
+    MultivariatePolynomial<AlgebraicExtension<R>, E>
+{
     /// Get the norm of a non-constant square-free polynomial `f` in the algebraic number field.
     pub fn norm(&self) -> MultivariatePolynomial<R, E> {
         self.norm_impl().3
diff --git a/src/domains/factorized_rational_polynomial.rs b/src/domains/factorized_rational_polynomial.rs
index 5a4e04a..102fb5a 100644
--- a/src/domains/factorized_rational_polynomial.rs
+++ b/src/domains/factorized_rational_polynomial.rs
@@ -9,8 +9,8 @@ use std::{
 
 use crate::{
     poly::{
-        factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial, Exponent,
-        Variable,
+        factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial,
+        PositiveExponent, Variable,
     },
     printer::{PrintOptions, PrintState},
 };
@@ -23,13 +23,13 @@ use super::{
 };
 
 #[derive(Clone, PartialEq, Eq, Hash, Debug)]
-pub struct FactorizedRationalPolynomialField<R: Ring, E: Exponent> {
+pub struct FactorizedRationalPolynomialField<R: Ring, E: PositiveExponent> {
     ring: R,
     var_map: Arc<Vec<Variable>>,
     _phantom_exp: PhantomData<E>,
 }
 
-impl<R: Ring, E: Exponent> FactorizedRationalPolynomialField<R, E> {
+impl<R: Ring, E: PositiveExponent> FactorizedRationalPolynomialField<R, E> {
     pub fn new(
         coeff_ring: R,
         var_map: Arc<Vec<Variable>>,
@@ -52,7 +52,7 @@ impl<R: Ring, E: Exponent> FactorizedRationalPolynomialField<R, E> {
     }
 }
 
-pub trait FromNumeratorAndFactorizedDenominator<R: Ring, OR: Ring, E: Exponent> {
+pub trait FromNumeratorAndFactorizedDenominator<R: Ring, OR: Ring, E: PositiveExponent> {
     /// Construct a rational polynomial from a numerator and a factorized denominator.
     /// An empty denominator means a denominator of 1.
     fn from_num_den(
@@ -64,21 +64,21 @@ pub trait FromNumeratorAndFactorizedDenominator<R: Ring, OR: Ring, E: Exponent>
 }
 
 #[derive(Clone, PartialEq, Eq, Hash, Debug)]
-pub struct FactorizedRationalPolynomial<R: Ring, E: Exponent> {
+pub struct FactorizedRationalPolynomial<R: Ring, E: PositiveExponent> {
     pub numerator: MultivariatePolynomial<R, E>,
     pub numer_coeff: R::Element,
     pub denom_coeff: R::Element,
     pub denominators: Vec<(MultivariatePolynomial<R, E>, usize)>, // TODO: sort factors?
 }
 
-impl<R: Ring, E: Exponent> InternalOrdering for FactorizedRationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> InternalOrdering for FactorizedRationalPolynomial<R, E> {
     /// An ordering of rational polynomials that has no intuitive meaning.
     fn internal_cmp(&self, _other: &Self) -> Ordering {
         todo!()
     }
 }
 
-impl<R: Ring, E: Exponent> FactorizedRationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> FactorizedRationalPolynomial<R, E> {
     pub fn new(field: &R, var_map: Arc<Vec<Variable>>) -> FactorizedRationalPolynomial<R, E> {
         let num = MultivariatePolynomial::new(field, None, var_map);
 
@@ -147,7 +147,7 @@ impl<R: Ring, E: Exponent> FactorizedRationalPolynomial<R, E> {
     }
 }
 
-impl<R: Ring, E: Exponent> SelfRing for FactorizedRationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> SelfRing for FactorizedRationalPolynomial<R, E> {
     fn is_zero(&self) -> bool {
         self.is_zero()
     }
@@ -385,7 +385,7 @@ impl<R: Ring, E: Exponent> SelfRing for FactorizedRationalPolynomial<R, E> {
     }
 }
 
-impl<E: Exponent> FromNumeratorAndFactorizedDenominator<RationalField, IntegerRing, E>
+impl<E: PositiveExponent> FromNumeratorAndFactorizedDenominator<RationalField, IntegerRing, E>
     for FactorizedRationalPolynomial<IntegerRing, E>
 {
     fn from_num_den(
@@ -428,7 +428,7 @@ impl<E: Exponent> FromNumeratorAndFactorizedDenominator<RationalField, IntegerRi
     }
 }
 
-impl<E: Exponent> FromNumeratorAndFactorizedDenominator<IntegerRing, IntegerRing, E>
+impl<E: PositiveExponent> FromNumeratorAndFactorizedDenominator<IntegerRing, IntegerRing, E>
     for FactorizedRationalPolynomial<IntegerRing, E>
 {
     fn from_num_den(
@@ -534,7 +534,7 @@ impl<E: Exponent> FromNumeratorAndFactorizedDenominator<IntegerRing, IntegerRing
     }
 }
 
-impl<UField: FiniteFieldWorkspace, E: Exponent>
+impl<UField: FiniteFieldWorkspace, E: PositiveExponent>
     FromNumeratorAndFactorizedDenominator<FiniteField<UField>, FiniteField<UField>, E>
     for FactorizedRationalPolynomial<FiniteField<UField>, E>
 where
@@ -618,7 +618,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> FactorizedRationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> FactorizedRationalPolynomial<R, E>
 where
     Self: FromNumeratorAndFactorizedDenominator<R, R, E>,
     MultivariatePolynomial<R, E>: Factorize,
@@ -644,7 +644,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> FactorizedRationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> FactorizedRationalPolynomial<R, E>
 where
     Self: FromNumeratorAndFactorizedDenominator<R, R, E>,
 {
@@ -727,20 +727,20 @@ where
     }
 }
 
-impl<R: Ring, E: Exponent> Display for FactorizedRationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> Display for FactorizedRationalPolynomial<R, E> {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         self.format(&PrintOptions::from_fmt(f), PrintState::from_fmt(f), f)
             .map(|_| ())
     }
 }
 
-impl<R: Ring, E: Exponent> Display for FactorizedRationalPolynomialField<R, E> {
+impl<R: Ring, E: PositiveExponent> Display for FactorizedRationalPolynomialField<R, E> {
     fn fmt(&self, _: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         Ok(())
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Ring
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Ring
     for FactorizedRationalPolynomialField<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
@@ -864,7 +864,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> EuclideanDomain
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> EuclideanDomain
     for FactorizedRationalPolynomialField<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
@@ -888,7 +888,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Field
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Field
     for FactorizedRationalPolynomialField<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
@@ -907,7 +907,7 @@ where
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E> + PolynomialGCD<E>, E: Exponent>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E> + PolynomialGCD<E>, E: PositiveExponent>
     Add<&'a FactorizedRationalPolynomial<R, E>> for &'b FactorizedRationalPolynomial<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
@@ -1014,7 +1014,8 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Sub for FactorizedRationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Sub
+    for FactorizedRationalPolynomial<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
 {
@@ -1025,7 +1026,7 @@ where
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
     Sub<&'a FactorizedRationalPolynomial<R, E>> for &'b FactorizedRationalPolynomial<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
@@ -1037,7 +1038,8 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Neg for FactorizedRationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Neg
+    for FactorizedRationalPolynomial<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
 {
@@ -1052,7 +1054,7 @@ where
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
     Mul<&'a FactorizedRationalPolynomial<R, E>> for &'b FactorizedRationalPolynomial<R, E>
 {
     type Output = FactorizedRationalPolynomial<R, E>;
@@ -1145,7 +1147,7 @@ impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent>
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
     Div<&'a FactorizedRationalPolynomial<R, E>> for &'b FactorizedRationalPolynomial<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
@@ -1256,7 +1258,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> FactorizedRationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> FactorizedRationalPolynomial<R, E>
 where
     FactorizedRationalPolynomial<R, E>: FromNumeratorAndFactorizedDenominator<R, R, E>,
     MultivariatePolynomial<R, E>: Factorize,
diff --git a/src/domains/rational_polynomial.rs b/src/domains/rational_polynomial.rs
index ee95ee1..2c7e91b 100644
--- a/src/domains/rational_polynomial.rs
+++ b/src/domains/rational_polynomial.rs
@@ -12,7 +12,7 @@ use ahash::HashMap;
 use crate::{
     poly::{
         factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial,
-        univariate::UnivariatePolynomial, Exponent, Variable,
+        univariate::UnivariatePolynomial, PositiveExponent, Variable,
     },
     printer::{PrintOptions, PrintState},
 };
@@ -25,12 +25,12 @@ use super::{
 };
 
 #[derive(Clone, PartialEq, Eq, Hash, Debug)]
-pub struct RationalPolynomialField<R: Ring, E: Exponent> {
+pub struct RationalPolynomialField<R: Ring, E: PositiveExponent> {
     ring: R,
     _phantom_exp: PhantomData<E>,
 }
 
-impl<R: Ring, E: Exponent> RationalPolynomialField<R, E> {
+impl<R: Ring, E: PositiveExponent> RationalPolynomialField<R, E> {
     pub fn new(coeff_ring: R) -> RationalPolynomialField<R, E> {
         RationalPolynomialField {
             ring: coeff_ring,
@@ -46,7 +46,7 @@ impl<R: Ring, E: Exponent> RationalPolynomialField<R, E> {
     }
 }
 
-pub trait FromNumeratorAndDenominator<R: Ring, OR: Ring, E: Exponent> {
+pub trait FromNumeratorAndDenominator<R: Ring, OR: Ring, E: PositiveExponent> {
     fn from_num_den(
         num: MultivariatePolynomial<R, E>,
         den: MultivariatePolynomial<R, E>,
@@ -56,12 +56,12 @@ pub trait FromNumeratorAndDenominator<R: Ring, OR: Ring, E: Exponent> {
 }
 
 #[derive(Clone, PartialEq, Eq, Hash, Debug)]
-pub struct RationalPolynomial<R: Ring, E: Exponent> {
+pub struct RationalPolynomial<R: Ring, E: PositiveExponent> {
     pub numerator: MultivariatePolynomial<R, E>,
     pub denominator: MultivariatePolynomial<R, E>,
 }
 
-impl<R: Ring, E: Exponent> InternalOrdering for RationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> InternalOrdering for RationalPolynomial<R, E> {
     /// An ordering of rational polynomials that has no intuitive meaning.
     fn internal_cmp(&self, other: &Self) -> Ordering {
         self.numerator
@@ -81,7 +81,7 @@ impl<R: Ring, E: Exponent> InternalOrdering for RationalPolynomial<R, E> {
     }
 }
 
-impl<R: Ring, E: Exponent> From<MultivariatePolynomial<R, E>> for RationalPolynomial<R, E>
+impl<R: Ring, E: PositiveExponent> From<MultivariatePolynomial<R, E>> for RationalPolynomial<R, E>
 where
     Self: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -92,7 +92,7 @@ where
     }
 }
 
-impl<R: Ring, E: Exponent> RationalPolynomial<R, E>
+impl<R: Ring, E: PositiveExponent> RationalPolynomial<R, E>
 where
     Self: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -120,7 +120,7 @@ where
     }
 }
 
-impl<R: Ring, E: Exponent> RationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> RationalPolynomial<R, E> {
     pub fn new(field: &R, var_map: Arc<Vec<Variable>>) -> RationalPolynomial<R, E> {
         let num = MultivariatePolynomial::new(field, None, var_map);
         let den = num.one();
@@ -165,7 +165,7 @@ impl<R: Ring, E: Exponent> RationalPolynomial<R, E> {
     }
 }
 
-impl<R: Ring, E: Exponent> SelfRing for RationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> SelfRing for RationalPolynomial<R, E> {
     fn is_zero(&self) -> bool {
         self.is_zero()
     }
@@ -244,7 +244,7 @@ impl<R: Ring, E: Exponent> SelfRing for RationalPolynomial<R, E> {
     }
 }
 
-impl<E: Exponent> FromNumeratorAndDenominator<RationalField, IntegerRing, E>
+impl<E: PositiveExponent> FromNumeratorAndDenominator<RationalField, IntegerRing, E>
     for RationalPolynomial<IntegerRing, E>
 {
     fn from_num_den(
@@ -293,7 +293,7 @@ impl<E: Exponent> FromNumeratorAndDenominator<RationalField, IntegerRing, E>
     }
 }
 
-impl<E: Exponent> FromNumeratorAndDenominator<IntegerRing, IntegerRing, E>
+impl<E: PositiveExponent> FromNumeratorAndDenominator<IntegerRing, IntegerRing, E>
     for RationalPolynomial<IntegerRing, E>
 {
     fn from_num_den(
@@ -333,7 +333,7 @@ impl<E: Exponent> FromNumeratorAndDenominator<IntegerRing, IntegerRing, E>
     }
 }
 
-impl<UField: FiniteFieldWorkspace, E: Exponent>
+impl<UField: FiniteFieldWorkspace, E: PositiveExponent>
     FromNumeratorAndDenominator<FiniteField<UField>, FiniteField<UField>, E>
     for RationalPolynomial<FiniteField<UField>, E>
 where
@@ -378,7 +378,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> RationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> RationalPolynomial<R, E>
 where
     Self: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -550,20 +550,21 @@ where
     }
 }
 
-impl<R: Ring, E: Exponent> Display for RationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> Display for RationalPolynomial<R, E> {
     fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         self.format(&PrintOptions::from_fmt(f), PrintState::from_fmt(f), f)
             .map(|_| ())
     }
 }
 
-impl<R: Ring, E: Exponent> Display for RationalPolynomialField<R, E> {
+impl<R: Ring, E: PositiveExponent> Display for RationalPolynomialField<R, E> {
     fn fmt(&self, _: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
         Ok(())
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Ring for RationalPolynomialField<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Ring
+    for RationalPolynomialField<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -684,7 +685,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> EuclideanDomain
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> EuclideanDomain
     for RationalPolynomialField<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
@@ -705,7 +706,8 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Field for RationalPolynomialField<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Field
+    for RationalPolynomialField<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -722,7 +724,7 @@ where
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E> + PolynomialGCD<E>, E: Exponent>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E> + PolynomialGCD<E>, E: PositiveExponent>
     Add<&'a RationalPolynomial<R, E>> for &'b RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
@@ -774,7 +776,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Sub for RationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Sub for RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -785,8 +787,8 @@ where
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Sub<&'a RationalPolynomial<R, E>>
-    for &'b RationalPolynomial<R, E>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
+    Sub<&'a RationalPolynomial<R, E>> for &'b RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -797,7 +799,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Neg for RationalPolynomial<R, E> {
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Neg for RationalPolynomial<R, E> {
     type Output = Self;
     fn neg(self) -> Self::Output {
         RationalPolynomial {
@@ -807,8 +809,8 @@ impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Neg for RationalPolynom
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Mul<&'a RationalPolynomial<R, E>>
-    for &'b RationalPolynomial<R, E>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
+    Mul<&'a RationalPolynomial<R, E>> for &'b RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -851,8 +853,8 @@ where
     }
 }
 
-impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Div<&'a RationalPolynomial<R, E>>
-    for &'b RationalPolynomial<R, E>
+impl<'a, 'b, R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
+    Div<&'a RationalPolynomial<R, E>> for &'b RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -864,7 +866,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> RationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -897,7 +899,8 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> Derivable for RationalPolynomialField<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> Derivable
+    for RationalPolynomialField<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
 {
@@ -910,7 +913,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> RationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
     MultivariatePolynomial<R, E>: Factorize,
@@ -999,7 +1002,7 @@ where
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> RationalPolynomial<R, E>
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> RationalPolynomial<R, E>
 where
     RationalPolynomial<R, E>: FromNumeratorAndDenominator<R, R, E>,
     MultivariatePolynomial<R, E>: Factorize,
diff --git a/src/parser.rs b/src/parser.rs
index 63b53a3..322e223 100644
--- a/src/parser.rs
+++ b/src/parser.rs
@@ -10,7 +10,7 @@ use crate::{
     atom::Atom,
     coefficient::{Coefficient, ConvertToRing},
     domains::{float::Float, integer::Integer, Ring},
-    poly::{polynomial::MultivariatePolynomial, Exponent, Variable},
+    poly::{polynomial::MultivariatePolynomial, PositiveExponent, Variable},
     state::{State, Workspace},
     LicenseManager,
 };
@@ -1085,7 +1085,7 @@ impl Token {
 
     /// A special routine that can parse a polynomial written in expanded form,
     /// where the coefficient comes first.
-    pub fn parse_polynomial<'a, R: Ring + ConvertToRing, E: Exponent>(
+    pub fn parse_polynomial<'a, R: Ring + ConvertToRing, E: PositiveExponent>(
         mut input: &'a [u8],
         var_map: &Arc<Vec<Variable>>,
         var_name_map: &[SmartString<LazyCompact>],
diff --git a/src/poly.rs b/src/poly.rs
index 1e8968e..a373908 100644
--- a/src/poly.rs
+++ b/src/poly.rs
@@ -27,7 +27,7 @@ use crate::domains::factorized_rational_polynomial::{
 };
 use crate::domains::integer::Integer;
 use crate::domains::rational_polynomial::{FromNumeratorAndDenominator, RationalPolynomial};
-use crate::domains::{EuclideanDomain, Ring, RingPrinter};
+use crate::domains::{EuclideanDomain, Ring, SelfRing};
 use crate::parser::{Operator, Token};
 use crate::printer::{PrintOptions, PrintState};
 use crate::state::{State, Workspace};
@@ -55,13 +55,14 @@ pub trait Exponent:
     + Copy
     + PartialEq
     + Eq
+    + TryFrom<i32>
 {
     fn zero() -> Self;
     fn one() -> Self;
-    /// Convert the exponent to `u32`. This is always possible, as `u32` is the largest supported exponent type.
-    fn to_u32(&self) -> u32;
-    /// Convert from `u32`. This function may panic if the exponent is too large.
-    fn from_u32(n: u32) -> Self;
+    /// Convert the exponent to `i32`. This is always possible, as `i32` is the largest supported exponent type.
+    fn to_i32(&self) -> i32;
+    /// Convert from `i32`. This function may panic if the exponent is too large.
+    fn from_i32(n: i32) -> Self;
     fn is_zero(&self) -> bool;
     fn checked_add(&self, other: &Self) -> Option<Self>;
     fn gcd(&self, other: &Self) -> Self;
@@ -93,13 +94,16 @@ impl Exponent for u32 {
     }
 
     #[inline]
-    fn to_u32(&self) -> u32 {
-        *self
+    fn to_i32(&self) -> i32 {
+        *self as i32
     }
 
     #[inline]
-    fn from_u32(n: u32) -> Self {
-        n
+    fn from_i32(n: i32) -> Self {
+        if n < 0 {
+            panic!("Exponent {} is negative", n);
+        }
+        n as u32
     }
 
     #[inline]
@@ -109,7 +113,7 @@ impl Exponent for u32 {
 
     #[inline]
     fn checked_add(&self, other: &Self) -> Option<Self> {
-        u32::checked_add(*self, *other)
+        i32::checked_add(*self as i32, *other as i32).map(|x| x as u32)
     }
 
     #[inline]
@@ -150,6 +154,77 @@ impl Exponent for u32 {
     }
 }
 
+impl Exponent for i32 {
+    #[inline]
+    fn zero() -> Self {
+        0
+    }
+
+    #[inline]
+    fn one() -> Self {
+        1
+    }
+
+    #[inline]
+    fn to_i32(&self) -> i32 {
+        *self
+    }
+
+    #[inline]
+    fn from_i32(n: i32) -> Self {
+        n
+    }
+
+    #[inline]
+    fn is_zero(&self) -> bool {
+        *self == 0
+    }
+
+    #[inline]
+    fn checked_add(&self, other: &Self) -> Option<Self> {
+        i32::checked_add(*self, *other)
+    }
+
+    #[inline]
+    fn gcd(&self, other: &Self) -> Self {
+        utils::gcd_signed(*self as i64, *other as i64) as Self
+    }
+
+    // Pack a list of positive exponents.
+    fn pack(list: &[Self]) -> u64 {
+        let mut num: u64 = 0;
+        for x in list.iter().rev() {
+            num = (num << 8) + (*x as u8 as u64);
+        }
+        num.swap_bytes()
+    }
+
+    fn unpack(mut n: u64, out: &mut [Self]) {
+        n = n.swap_bytes();
+        let s = unsafe { std::slice::from_raw_parts(&n as *const u64 as *const u8, out.len()) };
+        for (o, ss) in out.iter_mut().zip(s) {
+            *o = *ss as i32;
+        }
+    }
+
+    // Pack a list of positive exponents.
+    fn pack_u16(list: &[Self]) -> u64 {
+        let mut num: u64 = 0;
+        for x in list.iter().rev() {
+            num = (num << 16) + ((*x as u16).to_be() as u64);
+        }
+        num.swap_bytes()
+    }
+
+    fn unpack_u16(mut n: u64, out: &mut [Self]) {
+        n = n.swap_bytes();
+        let s = unsafe { std::slice::from_raw_parts(&n as *const u64 as *const u16, out.len()) };
+        for (o, ss) in out.iter_mut().zip(s) {
+            *o = ss.swap_bytes() as i32;
+        }
+    }
+}
+
 impl Exponent for u16 {
     #[inline]
     fn zero() -> Self {
@@ -162,13 +237,13 @@ impl Exponent for u16 {
     }
 
     #[inline]
-    fn to_u32(&self) -> u32 {
-        *self as u32
+    fn to_i32(&self) -> i32 {
+        *self as i32
     }
 
     #[inline]
-    fn from_u32(n: u32) -> Self {
-        if n <= u16::MAX as u32 {
+    fn from_i32(n: i32) -> Self {
+        if n >= 0 && n <= u16::MAX as i32 {
             n as u16
         } else {
             panic!("Exponent {} too large for u16", n);
@@ -223,6 +298,81 @@ impl Exponent for u16 {
     }
 }
 
+impl Exponent for i16 {
+    #[inline]
+    fn zero() -> Self {
+        0
+    }
+
+    #[inline]
+    fn one() -> Self {
+        1
+    }
+
+    #[inline]
+    fn to_i32(&self) -> i32 {
+        *self as i32
+    }
+
+    #[inline]
+    fn from_i32(n: i32) -> Self {
+        if n >= i16::MIN as i32 && n <= i16::MAX as i32 {
+            n as i16
+        } else {
+            panic!("Exponent {} too large for i16", n);
+        }
+    }
+
+    #[inline]
+    fn is_zero(&self) -> bool {
+        *self == 0
+    }
+
+    #[inline]
+    fn checked_add(&self, other: &Self) -> Option<Self> {
+        i16::checked_add(*self, *other)
+    }
+
+    #[inline]
+    fn gcd(&self, other: &Self) -> Self {
+        utils::gcd_signed(*self as i64, *other as i64) as Self
+    }
+
+    // Pack a list of positive exponents.
+    fn pack(list: &[Self]) -> u64 {
+        let mut num: u64 = 0;
+        for x in list.iter().rev() {
+            num = (num << 8) + (*x as u8 as u64);
+        }
+        num.swap_bytes()
+    }
+
+    fn unpack(mut n: u64, out: &mut [Self]) {
+        n = n.swap_bytes();
+        let s = unsafe { std::slice::from_raw_parts(&n as *const u64 as *const u8, out.len()) };
+        for (o, ss) in out.iter_mut().zip(s) {
+            *o = *ss as i16;
+        }
+    }
+
+    // Pack a list of positive exponents.
+    fn pack_u16(list: &[Self]) -> u64 {
+        let mut num: u64 = 0;
+        for x in list.iter().rev() {
+            num = (num << 16) + ((*x as u16).to_be() as u64);
+        }
+        num.swap_bytes()
+    }
+
+    fn unpack_u16(mut n: u64, out: &mut [Self]) {
+        n = n.swap_bytes();
+        let s = unsafe { std::slice::from_raw_parts(&n as *const u64 as *const u16, out.len()) };
+        for (o, ss) in out.iter_mut().zip(s) {
+            *o = ss.swap_bytes() as i16;
+        }
+    }
+}
+
 /// An exponent limited to 255 for efficiency
 impl Exponent for u8 {
     #[inline]
@@ -236,13 +386,13 @@ impl Exponent for u8 {
     }
 
     #[inline]
-    fn to_u32(&self) -> u32 {
-        *self as u32
+    fn to_i32(&self) -> i32 {
+        *self as i32
     }
 
     #[inline]
-    fn from_u32(n: u32) -> Self {
-        if n <= u8::MAX as u32 {
+    fn from_i32(n: i32) -> Self {
+        if n >= 0 && n <= u8::MAX as i32 {
             n as u8
         } else {
             panic!("Exponent {} too large for u8", n);
@@ -295,6 +445,150 @@ impl Exponent for u8 {
     }
 }
 
+impl Exponent for i8 {
+    #[inline]
+    fn zero() -> Self {
+        0
+    }
+
+    #[inline]
+    fn one() -> Self {
+        1
+    }
+
+    #[inline]
+    fn to_i32(&self) -> i32 {
+        *self as i32
+    }
+
+    #[inline]
+    fn from_i32(n: i32) -> Self {
+        if n >= i8::MIN as i32 && n <= i8::MAX as i32 {
+            n as i8
+        } else {
+            panic!("Exponent {} too large for i8", n);
+        }
+    }
+
+    #[inline]
+    fn is_zero(&self) -> bool {
+        *self == 0
+    }
+
+    #[inline]
+    fn checked_add(&self, other: &Self) -> Option<Self> {
+        i8::checked_add(*self, *other)
+    }
+
+    #[inline]
+    fn gcd(&self, other: &Self) -> Self {
+        utils::gcd_signed(*self as i64, *other as i64) as Self
+    }
+
+    // Pack a list of positive exponents.
+    fn pack(list: &[Self]) -> u64 {
+        let mut num: u64 = 0;
+        for x in list.iter().rev() {
+            num = (num << 8) + (*x as u8 as u64);
+        }
+        num.swap_bytes()
+    }
+
+    fn unpack(mut n: u64, out: &mut [Self]) {
+        n = n.swap_bytes();
+        let s = unsafe { std::slice::from_raw_parts(&n as *const u64 as *const u8, out.len()) };
+        for (o, ss) in out.iter_mut().zip(s) {
+            *o = *ss as i8;
+        }
+    }
+
+    // Pack a list of positive exponents.
+    fn pack_u16(list: &[Self]) -> u64 {
+        let mut num: u64 = 0;
+        for x in list.iter().rev() {
+            num = (num << 16) + ((*x as u16).to_be() as u64);
+        }
+        num.swap_bytes()
+    }
+
+    fn unpack_u16(mut n: u64, out: &mut [Self]) {
+        n = n.swap_bytes();
+        let s = unsafe { std::slice::from_raw_parts(&n as *const u64 as *const u16, out.len()) };
+        for (o, ss) in out.iter_mut().zip(s) {
+            *o = ss.swap_bytes() as i8;
+        }
+    }
+}
+
+pub trait PositiveExponent: Exponent {
+    fn from_u32(n: u32) -> Self {
+        if n > i32::MAX as u32 {
+            panic!("Exponent {} too large for i32", n);
+        }
+        Self::from_i32(n as i32)
+    }
+    fn to_u32(&self) -> u32;
+}
+
+impl PositiveExponent for u8 {
+    #[inline]
+    fn to_u32(&self) -> u32 {
+        *self as u32
+    }
+}
+impl PositiveExponent for u16 {
+    #[inline]
+    fn to_u32(&self) -> u32 {
+        *self as u32
+    }
+}
+impl PositiveExponent for u32 {
+    #[inline]
+    fn to_u32(&self) -> u32 {
+        *self
+    }
+}
+
+macro_rules! to_positive {
+    ($neg: ty, $pos: ty) => {
+        impl<R: Ring> MultivariatePolynomial<R, $neg> {
+            /// Convert a polynomial with positive exponents to its unsigned type equivalent
+            /// by a safe and almost zero-cost cast.
+            ///
+            /// Panics if the polynomial has negative exponents.
+            pub fn to_positive(self) -> MultivariatePolynomial<R, $pos> {
+                if !self.is_polynomial() {
+                    panic!("Polynomial has negative exponent");
+                }
+
+                unsafe { std::mem::transmute_copy(&std::mem::ManuallyDrop::new(self)) }
+            }
+        }
+
+        impl<R: Ring> MultivariatePolynomial<R, $pos> {
+            /// Convert a polynomial with positive exponents to its signed type equivalent
+            /// by a safe and almost zero-cost cast.
+            ///
+            /// Panics if the polynomial has exponents that are too large.
+            pub fn to_signed(self) -> MultivariatePolynomial<R, $neg> {
+                if self
+                    .exponents
+                    .iter()
+                    .any(|x| x.to_i32() > <$neg>::MAX as i32)
+                {
+                    panic!("Polynomial has exponents that are too large");
+                }
+
+                unsafe { std::mem::transmute_copy(&std::mem::ManuallyDrop::new(self)) }
+            }
+        }
+    };
+}
+
+to_positive!(i8, u8);
+to_positive!(i16, u16);
+to_positive!(i32, u32);
+
 /// A well-order of monomials.
 pub trait MonomialOrder: Clone {
     fn cmp<E: Exponent>(a: &[E], b: &[E]) -> Ordering;
@@ -344,7 +638,7 @@ impl MonomialOrder for LexOrder {
 /// A polynomial variable. It is either a (global) symbol
 /// a temporary variable (for internal use), an array entry,
 /// a function or any other non-polynomial part.
-#[derive(Clone, Hash, PartialEq, Eq, Debug)]
+#[derive(Clone, Hash, PartialEq, Eq, PartialOrd, Ord, Debug)]
 pub enum Variable {
     Symbol(Symbol),
     Temporary(usize), // a temporary variable, for internal use
@@ -394,19 +688,11 @@ impl Variable {
         }
     }
 
-    pub fn to_string_with_state(&self, state: PrintState) -> String {
+    fn format_string(&self, opts: &PrintOptions, state: PrintState) -> String {
         match self {
             Variable::Symbol(v) => State::get_name(*v).to_string(),
             Variable::Temporary(t) => format!("_TMP_{}", *t),
-            Variable::Function(_, a) | Variable::Other(a) => format!(
-                "{}",
-                RingPrinter {
-                    element: a.as_ref(),
-                    ring: &AtomField::new(),
-                    opts: PrintOptions::default(),
-                    state,
-                }
-            ),
+            Variable::Function(_, a) | Variable::Other(a) => a.format_string(opts, state),
         }
     }
 
@@ -462,6 +748,18 @@ impl Atom {
         self.as_view().to_polynomial(field, var_map)
     }
 
+    /// Convert the atom to a polynomial in specific variables.
+    /// All other parts will be collected into the coefficient, which
+    /// is a general expression.
+    ///
+    /// This routine does not perform expansions.
+    pub fn to_polynomial_in_vars<E: Exponent>(
+        &self,
+        var_map: &Arc<Vec<Variable>>,
+    ) -> MultivariatePolynomial<AtomField, E> {
+        self.as_view().to_polynomial_in_vars(var_map)
+    }
+
     /// Convert the atom to a rational polynomial, optionally in the variable ordering
     /// specified by `var_map`. If new variables are encountered, they are
     /// added to the variable map. Similarly, non-rational polynomial parts are automatically
@@ -469,7 +767,7 @@ impl Atom {
     pub fn to_rational_polynomial<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -491,7 +789,7 @@ impl Atom {
     pub fn to_factorized_rational_polynomial<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -660,11 +958,11 @@ impl<'a> AtomView<'a> {
                     match exp {
                         AtomView::Num(n) => match n.get_coeff_view() {
                             CoefficientView::Natural(r, _) => {
-                                exponents[var_index] += E::from_u32(r as u32)
+                                exponents[var_index] += E::from_i32(r as i32)
                             }
                             CoefficientView::Large(r) => {
                                 exponents[var_index] +=
-                                    E::from_u32(r.to_rat().numerator_ref().to_i64().unwrap() as u32)
+                                    E::from_i32(r.to_rat().numerator_ref().to_i64().unwrap() as i32)
                             }
                             _ => unreachable!(),
                         },
@@ -749,8 +1047,39 @@ impl<'a> AtomView<'a> {
                 if let AtomView::Num(n) = exp {
                     let num_n = n.get_coeff_view();
                     if let CoefficientView::Natural(nn, nd) = num_n {
-                        if nd == 1 && nn > 0 && nn < u32::MAX as i64 {
-                            return base.to_polynomial_impl(field, var_map).pow(nn as usize);
+                        if nd == 1 {
+                            if nn > 0 && nn < i32::MAX as i64 {
+                                return base.to_polynomial_impl(field, var_map).pow(nn as usize);
+                            } else if nn < 0 && nn > i32::MIN as i64 {
+                                // allow x^-2 as a term if supported by the exponent
+                                if let Ok(e) = (nn as i32).try_into() {
+                                    if let AtomView::Var(v) = base {
+                                        let s = Variable::Symbol(v.get_symbol());
+                                        if let Some(id) = var_map.iter().position(|v| v == &s) {
+                                            let mut exp = vec![E::zero(); var_map.len()];
+                                            exp[id] = e;
+                                            return MultivariatePolynomial::new(
+                                                field,
+                                                None,
+                                                var_map.clone(),
+                                            )
+                                            .monomial(field.one(), exp);
+                                        } else {
+                                            let mut var_map = var_map.as_ref().clone();
+                                            var_map.push(s);
+                                            let mut exp = vec![E::zero(); var_map.len()];
+                                            exp[var_map.len() - 1] = e;
+
+                                            return MultivariatePolynomial::new(
+                                                field,
+                                                None,
+                                                Arc::new(var_map),
+                                            )
+                                            .monomial(field.one(), exp);
+                                        }
+                                    }
+                                }
+                            }
                         }
                     }
                 }
@@ -829,6 +1158,22 @@ impl<'a> AtomView<'a> {
     pub fn to_polynomial_in_vars<E: Exponent>(
         &self,
         var_map: &Arc<Vec<Variable>>,
+    ) -> MultivariatePolynomial<AtomField, E> {
+        let poly = MultivariatePolynomial::<_, E>::new(&AtomField::new(), None, var_map.clone());
+        let mut cache = Vec::new();
+        self.to_polynomial_in_vars_impl(var_map, &poly, &mut cache)
+    }
+
+    /// Convert the atom to a polynomial in specific variables.
+    /// All other parts will be collected into the coefficient, which
+    /// is a general expression.
+    ///
+    /// This routine does not perform expansions.
+    fn to_polynomial_in_vars_impl<E: Exponent>(
+        &self,
+        var_map: &Arc<Vec<Variable>>,
+        poly: &MultivariatePolynomial<AtomField, E>,
+        cache: &mut Vec<MultivariatePolynomial<AtomField, E>>,
     ) -> MultivariatePolynomial<AtomField, E> {
         let field = AtomField::new();
         // see if the current term can be cast into a polynomial using a fast routine
@@ -837,9 +1182,7 @@ impl<'a> AtomView<'a> {
         }
 
         match self {
-            AtomView::Num(_) | AtomView::Var(_) => {
-                MultivariatePolynomial::new(&field, None, var_map.clone()).constant(self.to_owned())
-            }
+            AtomView::Num(_) | AtomView::Var(_) => poly.constant(self.to_owned()),
             AtomView::Pow(p) => {
                 let (base, exp) = p.get_base_exp();
 
@@ -847,7 +1190,12 @@ impl<'a> AtomView<'a> {
                     let num_n = n.get_coeff_view();
                     if let CoefficientView::Natural(nn, nd) = num_n {
                         if nd == 1 && nn > 0 && nn < u32::MAX as i64 {
-                            return base.to_polynomial_in_vars(var_map).pow(nn as usize);
+                            let mut b = base.to_polynomial_in_vars_impl(var_map, poly, cache);
+                            let r = b.pow(nn as usize);
+
+                            b.clear();
+                            cache.push(b);
+                            return r;
                         }
                     }
                 }
@@ -858,11 +1206,13 @@ impl<'a> AtomView<'a> {
                 }) {
                     let mut exp = vec![E::zero(); var_map.len()];
                     exp[id] = E::one();
-                    MultivariatePolynomial::new(&field, None, var_map.clone())
-                        .monomial(field.one(), exp)
+                    cache
+                        .pop()
+                        .unwrap_or_else(|| poly.monomial(field.one(), exp))
                 } else {
-                    MultivariatePolynomial::new(&field, None, var_map.clone())
-                        .constant(self.to_owned())
+                    cache
+                        .pop()
+                        .unwrap_or_else(|| poly.constant(self.to_owned()))
                 }
             }
             AtomView::Fun(_) => {
@@ -872,29 +1222,30 @@ impl<'a> AtomView<'a> {
                 }) {
                     let mut exp = vec![E::zero(); var_map.len()];
                     exp[id] = E::one();
-                    MultivariatePolynomial::new(&field, None, var_map.clone())
-                        .monomial(field.one(), exp)
+                    poly.monomial(field.one(), exp)
                 } else {
-                    MultivariatePolynomial::new(&field, None, var_map.clone())
-                        .constant(self.to_owned())
+                    poly.constant(self.to_owned())
                 }
             }
             AtomView::Mul(m) => {
-                let mut r = MultivariatePolynomial::new(&field, None, var_map.clone())
-                    .constant(field.one());
+                let mut r = cache.pop().unwrap_or_else(|| poly.one());
                 for arg in m {
-                    let mut arg_r = arg.to_polynomial_in_vars(&r.variables);
-                    r.unify_variables(&mut arg_r);
+                    let mut arg_r = arg.to_polynomial_in_vars_impl(&r.variables, poly, cache);
                     r = &r * &arg_r;
+
+                    arg_r.clear();
+                    cache.push(arg_r);
                 }
                 r
             }
             AtomView::Add(a) => {
-                let mut r = MultivariatePolynomial::new(&field, None, var_map.clone());
+                let mut r = cache.pop().unwrap_or_else(|| poly.zero());
                 for arg in a {
-                    let mut arg_r = arg.to_polynomial_in_vars(&r.variables);
-                    r.unify_variables(&mut arg_r);
+                    let mut arg_r = arg.to_polynomial_in_vars_impl(&r.variables, poly, cache);
                     r = &r + &arg_r;
+
+                    arg_r.clear();
+                    cache.push(arg_r);
                 }
                 r
             }
@@ -908,7 +1259,7 @@ impl<'a> AtomView<'a> {
     pub fn to_rational_polynomial<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -929,7 +1280,7 @@ impl<'a> AtomView<'a> {
     fn to_rational_polynomial_impl<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -1051,7 +1402,7 @@ impl<'a> AtomView<'a> {
     pub fn to_factorized_rational_polynomial<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -1073,7 +1424,7 @@ impl<'a> AtomView<'a> {
     pub fn to_factorized_rational_polynomial_impl<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -1210,33 +1561,27 @@ impl<E: Exponent, O: MonomialOrder> MultivariatePolynomial<AtomField, E, O> {
         let add = out.to_add();
 
         let mut mul_h = workspace.new_atom();
-        let mut var_h = workspace.new_atom();
         let mut num_h = workspace.new_atom();
         let mut pow_h = workspace.new_atom();
 
+        let vars: Vec<_> = self.variables.iter().map(|v| v.to_atom()).collect();
+
+        let mut sorted_vars = (0..vars.len()).collect::<Vec<_>>();
+        sorted_vars.sort_by_key(|&i| vars[i].clone());
+
         for monomial in self {
             let mul = mul_h.to_mul();
 
-            for (var_id, &pow) in self.variables.iter().zip(monomial.exponents) {
-                if pow > E::zero() {
-                    match var_id {
-                        Variable::Symbol(v) => {
-                            var_h.to_var(*v);
-                        }
-                        Variable::Temporary(_) => {
-                            unreachable!("Temporary variable in expression")
-                        }
-                        Variable::Function(_, a) | Variable::Other(a) => {
-                            var_h.set_from_view(&a.as_view());
-                        }
-                    }
-
-                    if pow > E::one() {
-                        num_h.to_num((pow.to_u32() as i64).into());
-                        pow_h.to_pow(var_h.as_view(), num_h.as_view());
+            for i in &sorted_vars {
+                let var = &vars[*i];
+                let pow = monomial.exponents[*i];
+                if pow != E::zero() {
+                    if pow != E::one() {
+                        num_h.to_num((pow.to_i32() as i64).into());
+                        pow_h.to_pow(var.as_view(), num_h.as_view());
                         mul.extend(pow_h.as_view());
                     } else {
-                        mul.extend(var_h.as_view());
+                        mul.extend(var.as_view());
                     }
                 }
             }
@@ -1284,34 +1629,38 @@ impl<R: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<R, E, O> {
         let add = out.to_add();
 
         let mut mul_h = workspace.new_atom();
-        let mut var_h = workspace.new_atom();
         let mut num_h = workspace.new_atom();
         let mut pow_h = workspace.new_atom();
 
+        let vars: Vec<_> = self
+            .variables
+            .iter()
+            .map(|v| {
+                if let Variable::Temporary(_) = v {
+                    let a = map.get(v).expect("Variable missing from map");
+                    a.to_owned()
+                } else {
+                    v.to_atom()
+                }
+            })
+            .collect();
+
+        let mut sorted_vars = (0..vars.len()).collect::<Vec<_>>();
+        sorted_vars.sort_by_key(|&i| vars[i].clone());
+
         for monomial in self {
             let mul = mul_h.to_mul();
 
-            for (var_id, &pow) in self.variables.iter().zip(monomial.exponents) {
-                if pow > E::zero() {
-                    match var_id {
-                        Variable::Symbol(v) => {
-                            var_h.to_var(*v);
-                        }
-                        Variable::Temporary(_) => {
-                            let a = map.get(var_id).expect("Variable missing from map");
-                            var_h.set_from_view(a);
-                        }
-                        Variable::Function(_, a) | Variable::Other(a) => {
-                            var_h.set_from_view(&a.as_view());
-                        }
-                    }
-
-                    if pow > E::one() {
-                        num_h.to_num((pow.to_u32() as i64).into());
-                        pow_h.to_pow(var_h.as_view(), num_h.as_view());
+            for i in &sorted_vars {
+                let var = &vars[*i];
+                let pow = monomial.exponents[*i];
+                if pow != E::zero() {
+                    if pow != E::one() {
+                        num_h.to_num((pow.to_i32() as i64).into());
+                        pow_h.to_pow(var.as_view(), num_h.as_view());
                         mul.extend(pow_h.as_view());
                     } else {
-                        mul.extend(var_h.as_view());
+                        mul.extend(var.as_view());
                     }
                 }
             }
@@ -1364,7 +1713,7 @@ impl<R: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<R, E, O> {
             let mul = mul_h.to_mul();
 
             for (var_id, &pow) in self.variables.iter().zip(monomial.exponents) {
-                if pow > E::zero() {
+                if pow != E::zero() {
                     match var_id {
                         Variable::Symbol(v) => {
                             var_h.to_var(*v);
@@ -1377,8 +1726,8 @@ impl<R: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<R, E, O> {
                         }
                     }
 
-                    if pow > E::one() {
-                        num_h.to_num((pow.to_u32() as i64).into());
+                    if pow != E::one() {
+                        num_h.to_num((pow.to_i32() as i64).into());
                         pow_h.to_pow(var_h.as_view(), num_h.as_view());
                         mul.extend(pow_h.as_view());
                     } else {
@@ -1398,7 +1747,7 @@ impl<R: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<R, E, O> {
     }
 }
 
-impl<R: Ring, E: Exponent> RationalPolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> RationalPolynomial<R, E> {
     pub fn to_expression(&self) -> Atom
     where
         R::Element: Into<Coefficient>,
@@ -1500,8 +1849,8 @@ impl Token {
 
                     match &args[1] {
                         Token::Number(n) => {
-                            if let Ok(x) = n.parse::<u32>() {
-                                exponents[var_index] += E::from_u32(x);
+                            if let Ok(x) = n.parse::<i32>() {
+                                exponents[var_index] += E::from_i32(x);
                             } else {
                                 Err("Invalid exponent")?
                             };
@@ -1596,7 +1945,7 @@ impl Token {
     pub fn to_rational_polynomial<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + ConvertToRing + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
@@ -1704,7 +2053,7 @@ impl Token {
     pub fn to_factorized_rational_polynomial<
         R: EuclideanDomain + ConvertToRing,
         RO: EuclideanDomain + ConvertToRing + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     >(
         &self,
         field: &R,
diff --git a/src/poly/evaluate.rs b/src/poly/evaluate.rs
index 4ba813a..dd7f15d 100644
--- a/src/poly/evaluate.rs
+++ b/src/poly/evaluate.rs
@@ -23,7 +23,7 @@ use crate::{
     state::State,
 };
 
-use super::{polynomial::MultivariatePolynomial, Exponent};
+use super::{polynomial::MultivariatePolynomial, PositiveExponent};
 
 /// A borrowed version of a Horner node, suitable as a key in a
 /// hashmap. It uses precomputed hashes for the complete node
@@ -315,7 +315,7 @@ where
     }
 }
 
-impl<E: Exponent> MultivariatePolynomial<RationalField, E> {
+impl<E: PositiveExponent> MultivariatePolynomial<RationalField, E> {
     /// Write the polynomial in a Horner scheme with the variable ordering
     /// defined in `order`.
     pub fn to_horner_scheme(&self, order: &[usize]) -> HornerScheme<RationalField> {
@@ -464,7 +464,7 @@ impl<E: Exponent> MultivariatePolynomial<RationalField, E> {
         let mut h = AHasher::default();
         h.write_u8(0);
         var.hash(&mut h);
-        (min_pow.to_u32() as usize).hash(&mut h);
+        (min_pow.to_i32() as usize).hash(&mut h);
 
         let pow_hash = h.finish(); // hash var^pow
 
@@ -513,7 +513,7 @@ impl<E: Exponent> MultivariatePolynomial<RationalField, E> {
 
         HornerScheme::Node(HornerNode {
             var,
-            pow: min_pow.to_u32() as usize,
+            pow: min_pow.to_i32() as usize,
             gcd,
             hash: (pow_hash, pow_content_hash, full_hash),
             content_rest: boxed_children,
@@ -545,7 +545,7 @@ impl<E: Exponent> MultivariatePolynomial<RationalField, E> {
 }
 
 impl HornerScheme<RationalField> {
-    pub fn optimize_multiple<E: Exponent>(
+    pub fn optimize_multiple<E: PositiveExponent>(
         polys: &[&MultivariatePolynomial<RationalField, E>],
         num_tries: usize,
     ) -> (Vec<HornerScheme<RationalField>>, usize, Vec<usize>) {
diff --git a/src/poly/factor.rs b/src/poly/factor.rs
index 3edbdd5..aaf9de8 100644
--- a/src/poly/factor.rs
+++ b/src/poly/factor.rs
@@ -20,7 +20,7 @@ use crate::{
     utils,
 };
 
-use super::{gcd::PolynomialGCD, polynomial::MultivariatePolynomial, Exponent, LexOrder};
+use super::{gcd::PolynomialGCD, polynomial::MultivariatePolynomial, LexOrder, PositiveExponent};
 
 pub trait Factorize: Sized {
     /// Perform a square-free factorization.
@@ -32,7 +32,9 @@ pub trait Factorize: Sized {
     fn is_irreducible(&self) -> bool;
 }
 
-impl<F: EuclideanDomain + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
+impl<F: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
+    MultivariatePolynomial<F, E, LexOrder>
+{
     /// Find factors that do not contain all variables.
     pub fn factor_separable(&self) -> Vec<Self> {
         let mut stripped = self.clone();
@@ -205,7 +207,7 @@ impl<F: EuclideanDomain + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<
     }
 }
 
-impl<E: Exponent> Factorize for MultivariatePolynomial<IntegerRing, E, LexOrder> {
+impl<E: PositiveExponent> Factorize for MultivariatePolynomial<IntegerRing, E, LexOrder> {
     fn square_free_factorization(&self) -> Vec<(Self, usize)> {
         if self.is_zero() {
             return vec![];
@@ -343,7 +345,7 @@ impl<E: Exponent> Factorize for MultivariatePolynomial<IntegerRing, E, LexOrder>
     }
 }
 
-impl<E: Exponent> Factorize for MultivariatePolynomial<RationalField, E, LexOrder> {
+impl<E: PositiveExponent> Factorize for MultivariatePolynomial<RationalField, E, LexOrder> {
     fn square_free_factorization(&self) -> Vec<(Self, usize)> {
         let c = self.content();
 
@@ -411,7 +413,7 @@ impl<E: Exponent> Factorize for MultivariatePolynomial<RationalField, E, LexOrde
     }
 }
 
-impl<E: Exponent> Factorize
+impl<E: PositiveExponent> Factorize
     for MultivariatePolynomial<AlgebraicExtension<RationalField>, E, LexOrder>
 {
     fn square_free_factorization(&self) -> Vec<(Self, usize)> {
@@ -493,7 +495,7 @@ impl<E: Exponent> Factorize
 impl<
         UField: FiniteFieldWorkspace,
         F: GaloisField<Base = FiniteField<UField>> + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     > Factorize for MultivariatePolynomial<F, E, LexOrder>
 where
     FiniteField<UField>: Field + FiniteFieldCore<UField> + PolynomialGCD<u16>,
@@ -648,7 +650,7 @@ where
 impl<
         UField: FiniteFieldWorkspace,
         F: GaloisField<Base = FiniteField<UField>> + PolynomialGCD<E>,
-        E: Exponent,
+        E: PositiveExponent,
     > MultivariatePolynomial<F, E, LexOrder>
 where
     FiniteField<UField>: Field + FiniteFieldCore<UField> + PolynomialGCD<u16>,
@@ -1650,7 +1652,7 @@ where
     }
 }
 
-impl<F: Field, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
+impl<F: Field, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
     fn multivariate_diophantine(
         univariate_deltas: &[Self],
         univariate_factors: &mut [Self],
@@ -1987,7 +1989,7 @@ impl<F: Field, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
     }
 }
 
-impl<E: Exponent> MultivariatePolynomial<IntegerRing, E, LexOrder> {
+impl<E: PositiveExponent> MultivariatePolynomial<IntegerRing, E, LexOrder> {
     /// Hensel lift a solution of `self = u * w mod p` to `self = u * w mod max_p`
     /// where `max_p` is a power of `p`.
     ///
@@ -2159,7 +2161,7 @@ impl<E: Exponent> MultivariatePolynomial<IntegerRing, E, LexOrder> {
         }
 
         let bound = self.coefficient_bound();
-        let p: Integer = (field.get_prime().to_u32() as i64).into();
+        let p: Integer = (field.get_prime() as i64).into();
         let mut max_p = p.clone();
         while max_p < bound {
             max_p = &max_p * &p;
@@ -3361,7 +3363,7 @@ impl<E: Exponent> MultivariatePolynomial<IntegerRing, E, LexOrder> {
     }
 }
 
-impl<E: Exponent> MultivariatePolynomial<FiniteField<Integer>, E, LexOrder> {
+impl<E: PositiveExponent> MultivariatePolynomial<FiniteField<Integer>, E, LexOrder> {
     /// Compute a univariate diophantine equation in `Z_p^k` by Newton iteration.
     fn get_univariate_factors_and_deltas(
         factors: &[Self],
diff --git a/src/poly/gcd.rs b/src/poly/gcd.rs
index b57bbcb..acd156c 100755
--- a/src/poly/gcd.rs
+++ b/src/poly/gcd.rs
@@ -18,7 +18,7 @@ use crate::poly::INLINED_EXPONENTS;
 use crate::tensors::matrix::{Matrix, MatrixError};
 
 use super::polynomial::MultivariatePolynomial;
-use super::Exponent;
+use super::PositiveExponent;
 
 // 100 large u32 primes starting from the 203213901st prime number
 pub const LARGE_U32_PRIMES: [u32; 100] = [
@@ -122,7 +122,7 @@ enum GCDError {
     BadCurrentImage,
 }
 
-impl<R: Ring, E: Exponent> MultivariatePolynomial<R, E> {
+impl<R: Ring, E: PositiveExponent> MultivariatePolynomial<R, E> {
     /// Evaluation of the exponents by filling in the variables
     #[inline(always)]
     fn evaluate_exponents(
@@ -184,7 +184,7 @@ impl<R: Ring, E: Exponent> MultivariatePolynomial<R, E> {
     }
 }
 
-impl<F: Field, E: Exponent> MultivariatePolynomial<F, E> {
+impl<F: Field, E: PositiveExponent> MultivariatePolynomial<F, E> {
     /// Compute the univariate GCD using Euclid's algorithm. The result is normalized to 1.
     pub fn univariate_gcd(&self, b: &Self) -> Self {
         if self.is_zero() {
@@ -1077,7 +1077,7 @@ impl<F: Field, E: Exponent> MultivariatePolynomial<F, E> {
     }
 }
 
-impl<F: Field + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<F, E> {
+impl<F: Field + PolynomialGCD<E>, E: PositiveExponent> MultivariatePolynomial<F, E> {
     /// Compute the gcd shape of two polynomials in a finite field by filling in random
     /// numbers.
     #[instrument(level = "debug", skip_all)]
@@ -1376,7 +1376,7 @@ impl<F: Field + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<F, E> {
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<R, E> {
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> MultivariatePolynomial<R, E> {
     /// Get the content of a multivariate polynomial viewed as a
     /// univariate polynomial in `x`.
     pub fn univariate_content(&self, x: usize) -> MultivariatePolynomial<R, E> {
@@ -1628,7 +1628,7 @@ impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<
 
         /// Undo simplifications made to the input polynomials and normalize the gcd.
         #[inline(always)]
-        fn rescale_gcd<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent>(
+        fn rescale_gcd<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>(
             mut g: MultivariatePolynomial<R, E>,
             shared_degree: &[E],
             base_degree: &[Option<E>],
@@ -1852,11 +1852,11 @@ pub enum HeuristicGCDError {
     BadReconstruction,
 }
 
-impl<E: Exponent> MultivariatePolynomial<IntegerRing, E> {
+impl<E: PositiveExponent> MultivariatePolynomial<IntegerRing, E> {
     /// Perform a heuristic GCD algorithm.
     #[instrument(level = "debug", skip_all)]
     pub fn heuristic_gcd(&self, b: &Self) -> Result<(Self, Self, Self), HeuristicGCDError> {
-        fn interpolate<E: Exponent>(
+        fn interpolate<E: PositiveExponent>(
             mut gamma: MultivariatePolynomial<IntegerRing, E>,
             var: usize,
             xi: &Integer,
@@ -2391,7 +2391,7 @@ impl<E: Exponent> MultivariatePolynomial<IntegerRing, E> {
 }
 
 /// Polynomial GCD functions for a certain coefficient type `Self`.
-pub trait PolynomialGCD<E: Exponent>: Ring {
+pub trait PolynomialGCD<E: PositiveExponent>: Ring {
     fn heuristic_gcd(
         a: &MultivariatePolynomial<Self, E>,
         b: &MultivariatePolynomial<Self, E>,
@@ -2417,7 +2417,7 @@ pub trait PolynomialGCD<E: Exponent>: Ring {
     fn normalize(a: MultivariatePolynomial<Self, E>) -> MultivariatePolynomial<Self, E>;
 }
 
-impl<E: Exponent> PolynomialGCD<E> for IntegerRing {
+impl<E: PositiveExponent> PolynomialGCD<E> for IntegerRing {
     fn heuristic_gcd(
         a: &MultivariatePolynomial<Self, E>,
         b: &MultivariatePolynomial<Self, E>,
@@ -2530,7 +2530,7 @@ impl<E: Exponent> PolynomialGCD<E> for IntegerRing {
     }
 }
 
-impl<E: Exponent> PolynomialGCD<E> for RationalField {
+impl<E: PositiveExponent> PolynomialGCD<E> for RationalField {
     fn heuristic_gcd(
         _a: &MultivariatePolynomial<Self, E>,
         _b: &MultivariatePolynomial<Self, E>,
@@ -2605,8 +2605,11 @@ impl<E: Exponent> PolynomialGCD<E> for RationalField {
     }
 }
 
-impl<UField: FiniteFieldWorkspace, F: GaloisField<Base = FiniteField<UField>>, E: Exponent>
-    PolynomialGCD<E> for F
+impl<
+        UField: FiniteFieldWorkspace,
+        F: GaloisField<Base = FiniteField<UField>>,
+        E: PositiveExponent,
+    > PolynomialGCD<E> for F
 where
     FiniteField<UField>: FiniteFieldCore<UField>,
     <FiniteField<UField> as Ring>::Element: Copy,
@@ -2669,7 +2672,7 @@ where
     }
 }
 
-impl<E: Exponent> PolynomialGCD<E> for AlgebraicExtension<RationalField> {
+impl<E: PositiveExponent> PolynomialGCD<E> for AlgebraicExtension<RationalField> {
     fn heuristic_gcd(
         _a: &MultivariatePolynomial<Self, E>,
         _b: &MultivariatePolynomial<Self, E>,
diff --git a/src/poly/polynomial.rs b/src/poly/polynomial.rs
index 275a6b3..4af53e4 100755
--- a/src/poly/polynomial.rs
+++ b/src/poly/polynomial.rs
@@ -1,5 +1,5 @@
 use ahash::{HashMap, HashMapExt};
-use std::cell::Cell;
+use std::cell::{Cell, UnsafeCell};
 use std::cmp::{Ordering, Reverse};
 use std::collections::{BTreeMap, BinaryHeap};
 use std::fmt::Display;
@@ -16,7 +16,7 @@ use crate::printer::{PrintOptions, PrintState};
 
 use super::gcd::PolynomialGCD;
 use super::univariate::UnivariatePolynomial;
-use super::{Exponent, LexOrder, MonomialOrder, Variable, INLINED_EXPONENTS};
+use super::{Exponent, LexOrder, MonomialOrder, PositiveExponent, Variable, INLINED_EXPONENTS};
 use smallvec::{smallvec, SmallVec};
 
 const MAX_DENSE_MUL_BUFFER_SIZE: usize = 1 << 24;
@@ -160,7 +160,9 @@ impl<R: Ring, E: Exponent> Ring for PolynomialRing<R, E> {
     }
 }
 
-impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> EuclideanDomain for PolynomialRing<R, E> {
+impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> EuclideanDomain
+    for PolynomialRing<R, E>
+{
     fn rem(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
         a.rem(b)
     }
@@ -175,6 +177,7 @@ impl<R: EuclideanDomain + PolynomialGCD<E>, E: Exponent> EuclideanDomain for Pol
 }
 
 /// Multivariate polynomial with a sparse degree and dense variable representation.
+/// Negative exponents are supported, if they are allowed by the exponent type.
 #[derive(Clone)]
 pub struct MultivariatePolynomial<F: Ring, E: Exponent = u16, O: MonomialOrder = LexOrder> {
     // Data format: the i-th monomial is stored as coefficients[i] and
@@ -429,7 +432,7 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
         self.exponents.chunks_mut(nvars)
     }
 
-    /// Returns the number of variables in the polynomial.
+    /// Reset the polynomial to 0.
     #[inline]
     pub fn clear(&mut self) {
         self.coefficients.clear();
@@ -446,7 +449,7 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
         self.variables.as_ref()
     }
 
-    /// Renaname a variable.
+    /// Rename a variable.
     pub fn rename_variable(&mut self, old: &Variable, new: &Variable) {
         if let Some(pos) = self.variables.iter().position(|v| v == old) {
             let mut new_vars = self.variables.as_ref().clone();
@@ -493,6 +496,26 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
         self.variables = Arc::new(new_var_map);
         self.exponents = newexp;
 
+        // check if term ordering remains unchanged
+        if new_var_pos_other.windows(2).all(|w| w[0] <= w[1]) {
+            let mut newexp = vec![E::zero(); self.nvars() * other.nterms()];
+
+            if other.nvars() > 0 {
+                for (d, t) in newexp
+                    .chunks_mut(self.nvars())
+                    .zip(other.exponents.chunks(other.nvars()))
+                {
+                    for (var, e) in t.iter().enumerate() {
+                        d[new_var_pos_other[var]] = *e;
+                    }
+                }
+            }
+
+            other.variables = self.variables.clone();
+            other.exponents = newexp;
+            return;
+        }
+
         // reconstruct 'other' with correct monomial ordering
         let mut newother = Self::new(&other.ring, other.nterms().into(), self.variables.clone());
         let mut newexp = vec![E::zero(); self.nvars()];
@@ -688,25 +711,6 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
         let i = l * self.nvars();
         self.exponents.splice(i..i, exponents.iter().cloned());
     }
-
-    /// Take the derivative of the polynomial w.r.t the variable `var`.
-    pub fn derivative(&self, var: usize) -> Self {
-        debug_assert!(var < self.nvars());
-
-        let mut res = self.zero_with_capacity(self.nterms());
-
-        let mut exp = vec![E::zero(); self.nvars()];
-        for x in self {
-            if x.exponents[var] > E::zero() {
-                exp.copy_from_slice(x.exponents);
-                let pow = exp[var].to_u32() as u64;
-                exp[var] = exp[var] - E::one();
-                res.append_monomial(self.ring.mul(x.coefficient, &self.ring.nth(pow)), &exp);
-            }
-        }
-
-        res
-    }
 }
 
 impl<F: Ring, E: Exponent, O: MonomialOrder> SelfRing for MultivariatePolynomial<F, E, O> {
@@ -761,10 +765,13 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> SelfRing for MultivariatePolynomial
             .as_ref()
             .iter()
             .map(|v| {
-                v.to_string_with_state(PrintState {
-                    in_exp: true,
-                    ..state
-                })
+                v.format_string(
+                    opts,
+                    PrintState {
+                        in_exp: true,
+                        ..state
+                    },
+                )
             })
             .collect();
 
@@ -787,7 +794,7 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> SelfRing for MultivariatePolynomial
 
                 f.write_str(var_id)?;
 
-                if e.to_u32() != 1 {
+                if e.to_i32() != 1 {
                     if opts.latex {
                         write!(f, "^{{{}}}", e)?;
                     } else if opts.double_star_for_exponentiation {
@@ -1158,8 +1165,8 @@ impl<'a, F: Ring, E: Exponent> Mul<&'a MultivariatePolynomial<F, E, LexOrder>>
     }
 }
 
-impl<'a, 'b, F: EuclideanDomain, E: Exponent> Div<&'a MultivariatePolynomial<F, E, LexOrder>>
-    for &'b MultivariatePolynomial<F, E, LexOrder>
+impl<'a, 'b, F: EuclideanDomain, E: PositiveExponent>
+    Div<&'a MultivariatePolynomial<F, E, LexOrder>> for &'b MultivariatePolynomial<F, E, LexOrder>
 {
     type Output = MultivariatePolynomial<F, E, LexOrder>;
 
@@ -1169,7 +1176,7 @@ impl<'a, 'b, F: EuclideanDomain, E: Exponent> Div<&'a MultivariatePolynomial<F,
     }
 }
 
-impl<'a, F: EuclideanDomain, E: Exponent> Div<&'a MultivariatePolynomial<F, E, LexOrder>>
+impl<'a, F: EuclideanDomain, E: PositiveExponent> Div<&'a MultivariatePolynomial<F, E, LexOrder>>
     for MultivariatePolynomial<F, E, LexOrder>
 {
     type Output = MultivariatePolynomial<F, E, LexOrder>;
@@ -1345,7 +1352,335 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
     }
 }
 
+impl<F: Ring, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
+    /// Take the derivative of the polynomial w.r.t the variable `var`.
+    pub fn derivative(&self, var: usize) -> Self {
+        debug_assert!(var < self.nvars());
+
+        let mut res = self.zero_with_capacity(self.nterms());
+
+        let mut exp = vec![E::zero(); self.nvars()];
+        for x in self {
+            if x.exponents[var] > E::zero() {
+                exp.copy_from_slice(x.exponents);
+                let pow = exp[var].to_i32() as u64;
+                exp[var] = exp[var] - E::one();
+                res.append_monomial(self.ring.mul(x.coefficient, &self.ring.nth(pow)), &exp);
+            }
+        }
+
+        res
+    }
+
+    /// Replace a variable `n` in the polynomial by an element from
+    /// the ring `v`.
+    pub fn replace(&self, n: usize, v: &F::Element) -> MultivariatePolynomial<F, E, LexOrder> {
+        if (n + 1..self.nvars()).all(|i| self.degree(i) == E::zero()) {
+            return self.replace_last(n, v);
+        }
+
+        let mut res = self.zero_with_capacity(self.nterms());
+        let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
+
+        // TODO: cache power taking?
+        for t in self {
+            if t.exponents[n] == E::zero() {
+                res.append_monomial(t.coefficient.clone(), t.exponents);
+                continue;
+            }
+
+            let c = self.ring.mul(
+                t.coefficient,
+                &self.ring.pow(v, t.exponents[n].to_i32() as u64),
+            );
+
+            e.copy_from_slice(t.exponents);
+            e[n] = E::zero();
+            res.append_monomial(c, &e);
+        }
+
+        res
+    }
+
+    /// Replace the last variable `n` in the polynomial by an element from
+    /// the ring `v`.
+    pub fn replace_last(&self, n: usize, v: &F::Element) -> MultivariatePolynomial<F, E, LexOrder> {
+        let mut res = self.zero_with_capacity(self.nterms());
+        let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
+
+        // TODO: cache power taking?
+        for t in self {
+            if t.exponents[n] == E::zero() {
+                res.append_monomial(t.coefficient.clone(), t.exponents);
+                continue;
+            }
+
+            let c = self.ring.mul(
+                t.coefficient,
+                &self.ring.pow(v, t.exponents[n].to_i32() as u64),
+            );
+
+            if F::is_zero(&c) {
+                continue;
+            }
+
+            e.copy_from_slice(t.exponents);
+            e[n] = E::zero();
+
+            if res.is_zero() || res.last_exponents() != e.as_slice() {
+                res.coefficients.push(c);
+                res.exponents.extend_from_slice(&e);
+            } else {
+                let l = res.coefficients.last_mut().unwrap();
+                self.ring.add_assign(l, &c);
+
+                if F::is_zero(l) {
+                    res.coefficients.pop();
+                    res.exponents.truncate(res.exponents.len() - self.nvars());
+                }
+            }
+        }
+
+        res
+    }
+
+    /// Replace a variable `n` in the polynomial by an element from
+    /// the ring `v`.
+    pub fn replace_all(&self, r: &[F::Element]) -> F::Element {
+        let mut res = self.ring.zero();
+
+        // TODO: cache power taking?
+        for t in self {
+            let mut c = t.coefficient.clone();
+
+            for (i, v) in r.iter().zip(t.exponents) {
+                if v != &E::zero() {
+                    self.ring
+                        .mul_assign(&mut c, &self.ring.pow(i, v.to_i32() as u64));
+                }
+            }
+
+            self.ring.add_assign(&mut res, &c);
+        }
+
+        res
+    }
+
+    /// Replace a variable `n` in the polynomial by a polynomial `v`.
+    pub fn replace_with_poly(&self, n: usize, v: &Self) -> Self {
+        assert_eq!(self.variables, v.variables);
+
+        if v.is_constant() {
+            return self.replace(n, &v.lcoeff());
+        }
+
+        let mut res = self.zero_with_capacity(self.nterms());
+        let mut exp = vec![E::zero(); self.nvars()];
+        for t in self {
+            if t.exponents[n] == E::zero() {
+                res.append_monomial(t.coefficient.clone(), &t.exponents[..self.nvars()]);
+                continue;
+            }
+
+            exp.copy_from_slice(t.exponents);
+            exp[n] = E::zero();
+
+            // TODO: cache v^e
+            res = res
+                + (&v.pow(t.exponents[n].to_i32() as usize)
+                    * &self.monomial(t.coefficient.clone(), exp.clone()));
+        }
+        res
+    }
+
+    /// Replace all variables except `v` in the polynomial by elements from
+    /// the ring.
+    pub fn replace_all_except(
+        &self,
+        v: usize,
+        r: &[(usize, F::Element)],
+        cache: &mut [Vec<F::Element>],
+    ) -> MultivariatePolynomial<F, E, LexOrder> {
+        let mut tm: HashMap<E, F::Element> = HashMap::new();
+
+        for t in self {
+            let mut c = t.coefficient.clone();
+            for (n, vv) in r {
+                let p = t.exponents[*n].to_i32() as usize;
+                if p > 0 {
+                    if p < cache[*n].len() {
+                        if F::is_zero(&cache[*n][p]) {
+                            cache[*n][p] = self.ring.pow(vv, p as u64);
+                        }
+
+                        self.ring.mul_assign(&mut c, &cache[*n][p]);
+                    } else {
+                        self.ring.mul_assign(&mut c, &self.ring.pow(vv, p as u64));
+                    }
+                }
+            }
+
+            tm.entry(t.exponents[v])
+                .and_modify(|e| self.ring.add_assign(e, &c))
+                .or_insert(c);
+        }
+
+        let mut res = self.zero();
+        let mut e = vec![E::zero(); self.nvars()];
+        for (k, c) in tm {
+            e[v] = k;
+            res.append_monomial(c, &e);
+            e[v] = E::zero();
+        }
+
+        res
+    }
+
+    /// Shift a variable `var` to `var+shift`.
+    pub fn shift_var(&self, var: usize, shift: &F::Element) -> Self {
+        let d = self.degree(var).to_i32() as usize;
+
+        let y_poly = self.to_univariate_polynomial_list(var);
+
+        let mut v = vec![self.zero(); d + 1];
+        for (x_poly, p) in y_poly {
+            v[p.to_i32() as usize] = x_poly;
+        }
+
+        for k in 0..d {
+            for j in (k..d).rev() {
+                v[j] = &v[j] + &v[j + 1].clone().mul_coeff(shift.clone());
+            }
+        }
+
+        let mut poly = self.zero();
+        for (i, mut v) in v.into_iter().enumerate() {
+            for x in v.exponents.chunks_mut(self.nvars()) {
+                x[var] = E::from_i32(i as i32);
+            }
+
+            for m in &v {
+                poly.append_monomial(m.coefficient.clone(), m.exponents);
+            }
+        }
+
+        poly
+    }
+
+    /// Synthetic division for univariate polynomials, where `div` is monic.
+    // TODO: create UnivariatePolynomial?
+    pub fn quot_rem_univariate_monic(
+        &self,
+        div: &MultivariatePolynomial<F, E, LexOrder>,
+    ) -> (
+        MultivariatePolynomial<F, E, LexOrder>,
+        MultivariatePolynomial<F, E, LexOrder>,
+    ) {
+        debug_assert_eq!(div.lcoeff(), self.ring.one());
+        if self.is_zero() {
+            return (self.clone(), self.clone());
+        }
+
+        let mut dividendpos = self.nterms() - 1; // work from the back
+
+        let mut q = self.zero_with_capacity(self.nterms());
+        let mut r = self.zero();
+
+        // determine the variable
+        let mut var = 0;
+        for (i, x) in self.last_exponents().iter().enumerate() {
+            if !x.is_zero() {
+                var = i;
+                break;
+            }
+        }
+
+        let m = div.ldegree_max();
+        let mut pow = self.ldegree_max();
+
+        loop {
+            // find the power in the dividend if it exists
+            let mut coeff = loop {
+                if self.exponents(dividendpos)[var] == pow {
+                    break self.coefficients[dividendpos].clone();
+                }
+                if dividendpos == 0 || self.exponents(dividendpos)[var] < pow {
+                    break self.ring.zero();
+                }
+                dividendpos -= 1;
+            };
+
+            let mut qindex = 0; // starting from highest
+            let mut bindex = 0; // starting from lowest
+            while bindex < div.nterms() && qindex < q.nterms() {
+                while bindex + 1 < div.nterms()
+                    && div.exponents(bindex)[var] + q.exponents(qindex)[var] < pow
+                {
+                    bindex += 1;
+                }
+
+                if div.exponents(bindex)[var] + q.exponents(qindex)[var] == pow {
+                    self.ring.sub_mul_assign(
+                        &mut coeff,
+                        &div.coefficients[bindex],
+                        &q.coefficients[qindex],
+                    );
+                }
+
+                qindex += 1;
+            }
+
+            if !F::is_zero(&coeff) {
+                // can the division be performed? if not, add to rest
+                // TODO: refactor
+                let (quot, div) = if pow >= m {
+                    (coeff, true)
+                } else {
+                    (coeff, false)
+                };
+
+                if div {
+                    let nterms = q.nterms();
+                    let nvars = q.nvars();
+                    q.coefficients.push(quot);
+                    q.exponents.resize((nterms + 1) * nvars, E::zero());
+                    q.exponents[nterms * nvars + var] = pow - m;
+                } else {
+                    let nterms = r.nterms();
+                    let nvars = r.nvars();
+                    r.coefficients.push(quot);
+                    r.exponents.resize((nterms + 1) * nvars, E::zero());
+                    r.exponents[nterms * nvars + var] = pow;
+                }
+            }
+
+            if pow.is_zero() {
+                break;
+            }
+
+            pow = pow - E::one();
+        }
+
+        q.reverse();
+        r.reverse();
+
+        #[cfg(debug_assertions)]
+        {
+            if !(&q * div + r.clone() - self.clone()).is_zero() {
+                panic!("Division failed: ({})/({}): q={}, r={}", self, div, q, r);
+            }
+        }
+
+        (q, r)
+    }
+}
+
 impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
+    /// Check if all exponents are positive.
+    pub fn is_polynomial(&self) -> bool {
+        self.is_zero() || self.exponents.iter().all(|e| *e >= E::zero())
+    }
+
     /// Get the leading coefficient under a given variable ordering.
     /// This operation is O(n) if the variables are out of order.
     pub fn lcoeff_varorder(&self, vars: &[usize]) -> F::Element {
@@ -1532,205 +1867,41 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
     /// The map can also be reversed, by setting `inverse` to `true`.
     pub fn rearrange(
         &self,
-        order: &[usize],
-        inverse: bool,
-    ) -> MultivariatePolynomial<F, E, LexOrder> {
-        self.rearrange_impl(order, inverse, true)
-    }
-
-    /// Change the order of the variables in the polynomial, using `order`.
-    /// The order may contain `None`, to signal unmapped indices. This operation
-    /// allows the polynomial to grow in size.
-    ///
-    /// Note that the polynomial `var_map` is not updated.
-    pub fn rearrange_with_growth(
-        &self,
-        order: &[Option<usize>],
-    ) -> MultivariatePolynomial<F, E, LexOrder> {
-        let mut new_exp = vec![E::zero(); self.nterms() * order.len()];
-        for (e, er) in new_exp.chunks_mut(order.len()).zip(self.exponents_iter()) {
-            for x in 0..order.len() {
-                if let Some(v) = order[x] {
-                    e[x] = er[v];
-                }
-            }
-        }
-
-        let mut indices: Vec<usize> = (0..self.nterms()).collect();
-        indices.sort_unstable_by_key(|&i| &new_exp[i * order.len()..(i + 1) * order.len()]);
-
-        let mut res =
-            MultivariatePolynomial::new(&self.ring, self.nterms().into(), self.variables.clone());
-
-        for i in indices {
-            res.append_monomial(
-                self.coefficients[i].clone(),
-                &new_exp[i * order.len()..(i + 1) * order.len()],
-            );
-        }
-
-        res
-    }
-
-    /// Replace a variable `n` in the polynomial by an element from
-    /// the ring `v`.
-    pub fn replace(&self, n: usize, v: &F::Element) -> MultivariatePolynomial<F, E, LexOrder> {
-        if (n + 1..self.nvars()).all(|i| self.degree(i) == E::zero()) {
-            return self.replace_last(n, v);
-        }
-
-        let mut res = self.zero_with_capacity(self.nterms());
-        let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
-
-        // TODO: cache power taking?
-        for t in self {
-            if t.exponents[n] == E::zero() {
-                res.append_monomial(t.coefficient.clone(), t.exponents);
-                continue;
-            }
-
-            let c = self.ring.mul(
-                t.coefficient,
-                &self.ring.pow(v, t.exponents[n].to_u32() as u64),
-            );
-
-            e.copy_from_slice(t.exponents);
-            e[n] = E::zero();
-            res.append_monomial(c, &e);
-        }
-
-        res
-    }
-
-    /// Replace the last variable `n` in the polynomial by an element from
-    /// the ring `v`.
-    pub fn replace_last(&self, n: usize, v: &F::Element) -> MultivariatePolynomial<F, E, LexOrder> {
-        let mut res = self.zero_with_capacity(self.nterms());
-        let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
-
-        // TODO: cache power taking?
-        for t in self {
-            if t.exponents[n] == E::zero() {
-                res.append_monomial(t.coefficient.clone(), t.exponents);
-                continue;
-            }
-
-            let c = self.ring.mul(
-                t.coefficient,
-                &self.ring.pow(v, t.exponents[n].to_u32() as u64),
-            );
-
-            if F::is_zero(&c) {
-                continue;
-            }
-
-            e.copy_from_slice(t.exponents);
-            e[n] = E::zero();
-
-            if res.is_zero() || res.last_exponents() != e.as_slice() {
-                res.coefficients.push(c);
-                res.exponents.extend_from_slice(&e);
-            } else {
-                let l = res.coefficients.last_mut().unwrap();
-                self.ring.add_assign(l, &c);
-
-                if F::is_zero(l) {
-                    res.coefficients.pop();
-                    res.exponents.truncate(res.exponents.len() - self.nvars());
-                }
-            }
-        }
-
-        res
-    }
-
-    /// Replace a variable `n` in the polynomial by an element from
-    /// the ring `v`.
-    pub fn replace_all(&self, r: &[F::Element]) -> F::Element {
-        let mut res = self.ring.zero();
-
-        // TODO: cache power taking?
-        for t in self {
-            let mut c = t.coefficient.clone();
-
-            for (i, v) in r.iter().zip(t.exponents) {
-                if v != &E::zero() {
-                    self.ring
-                        .mul_assign(&mut c, &self.ring.pow(i, v.to_u32() as u64));
-                }
-            }
-
-            self.ring.add_assign(&mut res, &c);
-        }
-
-        res
-    }
-
-    /// Replace a variable `n` in the polynomial by a polynomial `v`.
-    pub fn replace_with_poly(&self, n: usize, v: &Self) -> Self {
-        assert_eq!(self.variables, v.variables);
-
-        if v.is_constant() {
-            return self.replace(n, &v.lcoeff());
-        }
-
-        let mut res = self.zero_with_capacity(self.nterms());
-        let mut exp = vec![E::zero(); self.nvars()];
-        for t in self {
-            if t.exponents[n] == E::zero() {
-                res.append_monomial(t.coefficient.clone(), &t.exponents[..self.nvars()]);
-                continue;
-            }
-
-            exp.copy_from_slice(t.exponents);
-            exp[n] = E::zero();
-
-            // TODO: cache v^e
-            res = res
-                + (&v.pow(t.exponents[n].to_u32() as usize)
-                    * &self.monomial(t.coefficient.clone(), exp.clone()));
-        }
-        res
+        order: &[usize],
+        inverse: bool,
+    ) -> MultivariatePolynomial<F, E, LexOrder> {
+        self.rearrange_impl(order, inverse, true)
     }
 
-    /// Replace all variables except `v` in the polynomial by elements from
-    /// the ring.
-    pub fn replace_all_except(
+    /// Change the order of the variables in the polynomial, using `order`.
+    /// The order may contain `None`, to signal unmapped indices. This operation
+    /// allows the polynomial to grow in size.
+    ///
+    /// Note that the polynomial `var_map` is not updated.
+    pub fn rearrange_with_growth(
         &self,
-        v: usize,
-        r: &[(usize, F::Element)],
-        cache: &mut [Vec<F::Element>],
+        order: &[Option<usize>],
     ) -> MultivariatePolynomial<F, E, LexOrder> {
-        let mut tm: HashMap<E, F::Element> = HashMap::new();
-
-        for t in self {
-            let mut c = t.coefficient.clone();
-            for (n, vv) in r {
-                let p = t.exponents[*n].to_u32() as usize;
-                if p > 0 {
-                    if p < cache[*n].len() {
-                        if F::is_zero(&cache[*n][p]) {
-                            cache[*n][p] = self.ring.pow(vv, p as u64);
-                        }
-
-                        self.ring.mul_assign(&mut c, &cache[*n][p]);
-                    } else {
-                        self.ring.mul_assign(&mut c, &self.ring.pow(vv, p as u64));
-                    }
+        let mut new_exp = vec![E::zero(); self.nterms() * order.len()];
+        for (e, er) in new_exp.chunks_mut(order.len()).zip(self.exponents_iter()) {
+            for x in 0..order.len() {
+                if let Some(v) = order[x] {
+                    e[x] = er[v];
                 }
             }
-
-            tm.entry(t.exponents[v])
-                .and_modify(|e| self.ring.add_assign(e, &c))
-                .or_insert(c);
         }
 
-        let mut res = self.zero();
-        let mut e = vec![E::zero(); self.nvars()];
-        for (k, c) in tm {
-            e[v] = k;
-            res.append_monomial(c, &e);
-            e[v] = E::zero();
+        let mut indices: Vec<usize> = (0..self.nterms()).collect();
+        indices.sort_unstable_by_key(|&i| &new_exp[i * order.len()..(i + 1) * order.len()]);
+
+        let mut res =
+            MultivariatePolynomial::new(&self.ring, self.nterms().into(), self.variables.clone());
+
+        for i in indices {
+            res.append_monomial(
+                self.coefficients[i].clone(),
+                &new_exp[i * order.len()..(i + 1) * order.len()],
+            );
         }
 
         res
@@ -1774,9 +1945,13 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
             return p;
         }
 
-        p.coefficients = vec![self.zero(); c.last().unwrap().1.to_u32() as usize + 1];
+        p.coefficients = vec![self.zero(); c.last().unwrap().1.to_i32() as usize + 1];
         for (q, e) in c {
-            p.coefficients[e.to_u32() as usize] = q;
+            if e < E::zero() {
+                panic!("Negative exponent in univariate conversion");
+            }
+
+            p.coefficients[e.to_i32() as usize] = q;
         }
 
         p
@@ -1790,9 +1965,13 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
             return p;
         }
 
-        p.coefficients = vec![p.ring.zero(); self.degree(var).to_u32() as usize + 1];
+        p.coefficients = vec![p.ring.zero(); self.degree(var).to_i32() as usize + 1];
         for (q, e) in self.coefficients.iter().zip(self.exponents_iter()) {
-            p.coefficients[e[var].to_u32() as usize] = q.clone();
+            if e[var] < E::zero() {
+                panic!("Negative exponent in univariate conversion");
+            }
+
+            p.coefficients[e[var].to_i32() as usize] = q.clone();
         }
 
         p
@@ -1809,22 +1988,26 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         }
 
         // get maximum degree for variable x
+        let mut mindeg = E::zero();
         let mut maxdeg = E::zero();
         for t in 0..self.nterms() {
             let d = self.exponents(t)[x];
             if d > maxdeg {
                 maxdeg = d;
             }
+            if d < mindeg {
+                mindeg = d;
+            }
         }
 
         // construct the coefficient per power of x
         let mut result = vec![];
         let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
-        for d in 0..maxdeg.to_u32() + 1 {
+        for d in mindeg.to_i32()..maxdeg.to_i32() + 1 {
             // TODO: add bounds estimate
             let mut a = self.zero();
             for t in 0..self.nterms() {
-                if self.exponents(t)[x].to_u32() == d {
+                if self.exponents(t)[x].to_i32() == d {
                     for (i, ee) in self.exponents(t).iter().enumerate() {
                         e[i] = *ee;
                     }
@@ -1834,7 +2017,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
             }
 
             if !a.is_zero() {
-                result.push((a, E::from_u32(d)));
+                result.push((a, E::from_i32(d)));
             }
         }
 
@@ -1923,7 +2106,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
 
             c2.exponents = c
                 .exponents_iter()
-                .map(|x| x[var_index].to_u32() as u16)
+                .map(|x| x[var_index].to_i32() as u16)
                 .collect();
             c2.coefficients = c.coefficients;
 
@@ -1936,9 +2119,9 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         if self.is_constant() {
             if let Some(m) = max_pow {
                 if let Some(var) = rhs.last_exponents().iter().position(|e| *e != E::zero()) {
-                    if rhs.degree(var).to_u32() > m as u32 {
+                    if rhs.degree(var).to_i32() > m as i32 {
                         return rhs
-                            .mod_var(var, E::from_u32(m as u32 + 1))
+                            .mod_var(var, E::from_i32(m as i32 + 1))
                             .mul_coeff(self.lcoeff());
                     }
                 }
@@ -1949,9 +2132,9 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         if rhs.is_constant() {
             if let Some(m) = max_pow {
                 if let Some(var) = self.last_exponents().iter().position(|e| *e != E::zero()) {
-                    if self.degree(var).to_u32() > m as u32 {
+                    if self.degree(var).to_i32() > m as i32 {
                         return self
-                            .mod_var(var, E::from_u32(m as u32 + 1))
+                            .mod_var(var, E::from_i32(m as i32 + 1))
                             .mul_coeff(rhs.lcoeff());
                     }
                 }
@@ -1967,7 +2150,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
 
         let d1 = self.degree(var);
         let d2 = rhs.degree(var);
-        let mut max = (d1.to_u32() + d2.to_u32()) as usize;
+        let mut max = (d1.to_i32() + d2.to_i32()) as usize;
         if let Some(m) = max_pow {
             max = max.min(m);
         }
@@ -1976,7 +2159,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
 
         for x in self {
             for y in rhs {
-                let pos = x.exponents[var].to_u32() + y.exponents[var].to_u32();
+                let pos = x.exponents[var].to_i32() + y.exponents[var].to_i32();
                 if pos as usize > max {
                     continue;
                 }
@@ -1990,7 +2173,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         let mut res = self.zero_with_capacity(coeffs.len());
         for (p, c) in coeffs.into_iter().enumerate() {
             if !F::is_zero(&c) {
-                exp[var] = E::from_u32(p as u32);
+                exp[var] = E::from_i32(p as i32);
                 res.append_monomial(c, &exp);
             }
         }
@@ -2001,9 +2184,13 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         &self,
         rhs: &MultivariatePolynomial<F, E, LexOrder>,
     ) -> Option<MultivariatePolynomial<F, E, LexOrder>> {
+        if !self.is_polynomial() || !rhs.is_polynomial() {
+            return None;
+        }
+
         let max_degs_rev = (0..self.nvars())
             .rev()
-            .map(|i| 1 + self.degree(i).to_u32() as usize + rhs.degree(i).to_u32() as usize)
+            .map(|i| 1 + self.degree(i).to_i32() as usize + rhs.degree(i).to_i32() as usize)
             .collect::<Vec<_>>();
 
         if max_degs_rev.iter().filter(|x| **x > 1).count() == 1 {
@@ -2034,10 +2221,10 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         #[inline(always)]
         fn to_uni_var<E: Exponent>(s: &[E], max_degs_rev: &[usize]) -> u32 {
             let mut shift = 1;
-            let mut res = s.last().unwrap().to_u32();
+            let mut res = s.last().unwrap().to_i32() as u32;
             for (ee, &x) in s.iter().rev().skip(1).zip(max_degs_rev) {
-                shift = shift.to_u32() * x as u32;
-                res += ee.to_u32() * shift;
+                shift = shift * x as u32;
+                res += ee.to_i32() as u32 * shift;
             }
             res
         }
@@ -2045,7 +2232,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         #[inline(always)]
         fn from_uni_var<E: Exponent>(mut p: u32, max_degs_rev: &[usize], exp: &mut [E]) {
             for (ee, &x) in exp.iter_mut().rev().zip(max_degs_rev) {
-                *ee = E::from_u32(p % x as u32);
+                *ee = E::from_i32((p % x as u32) as i32);
                 p /= x as u32;
             }
         }
@@ -2069,7 +2256,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
 
             for (c1, e1) in self.coefficients.iter().zip(&uni_exp_self) {
                 for (c2, e2) in rhs.coefficients.iter().zip(&uni_exp_rhs) {
-                    let pos = e1.to_u32() as usize + e2.to_u32() as usize;
+                    let pos = *e1 as usize + *e2 as usize;
                     self.ring.add_mul_assign(&mut coeffs[pos], c1, c2);
                 }
             }
@@ -2093,7 +2280,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
 
             for (c1, e1) in self.coefficients.iter().zip(&uni_exp_self) {
                 for (c2, e2) in rhs.coefficients.iter().zip(&uni_exp_rhs) {
-                    let pos = e1.to_u32() as usize + e2.to_u32() as usize;
+                    let pos = *e1 as usize + *e2 as usize;
                     if coeff_index[pos] == 0 {
                         coeffs.push(self.ring.mul(c1, c2));
                         coeff_index[pos] = coeffs.len() as u32;
@@ -2124,17 +2311,6 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         }
     }
 
-    /// Multiplication for multivariate polynomials using a custom variation of the heap method
-    /// described in "Sparse polynomial division using a heap" by Monagan, Pearce (2011) and using
-    /// the sorting described in "Sparse Polynomial Powering Using Heaps".
-    /// It uses a heap to obtain the next monomial of the result in an ordered fashion.
-    /// Additionally, this method uses a hashmap with the monomial exponent as a key and a vector of all pairs
-    /// of indices in `self` and `other` that have that monomial exponent when multiplied together.
-    /// When a multiplication of two monomials is considered, its indices are added to the hashmap,
-    /// but they are only added to the heap if the monomial exponent is new. As a result, the heap
-    /// only has new monomials, and by taking (and removing) the corresponding entry from the hashmap, all
-    /// monomials that have that exponent can be summed. Then, new monomials combinations are added that
-    /// should be considered next as they are smaller than the current monomial.
     fn heap_mul(
         &self,
         rhs: &MultivariatePolynomial<F, E, LexOrder>,
@@ -2146,12 +2322,14 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         }
 
         let degree_sum: Vec<_> = (0..self.nvars())
-            .map(|i| self.degree(i).to_u32() as usize + rhs.degree(i).to_u32() as usize)
+            .map(|i| self.degree(i).to_i32() as i64 + rhs.degree(i).to_i32() as i64)
             .collect();
 
         // use a special routine if the exponents can be packed into a u64
         let mut pack_u8 = true;
         if self.nvars() <= 8
+            && self.is_polynomial()
+            && rhs.is_polynomial()
             && degree_sum.iter().all(|deg| {
                 if *deg > 255 {
                     pack_u8 = false;
@@ -2163,24 +2341,85 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
             return self.heap_mul_packed_exp(rhs, pack_u8);
         }
 
+        let mut monomials = Vec::with_capacity(self.nterms() * self.nvars());
+        monomials.extend(
+            self.exponents(0)
+                .iter()
+                .zip(rhs.exponents(0))
+                .map(|(e1, e2)| *e1 + *e2),
+        );
+
+        let monomials = UnsafeCell::new((self.nvars(), monomials));
+
+        /// In order to prevent allocations of the exponents, store them in a single
+        /// append-only vector and use a key to index into it. For performance,
+        /// we use an unsafe cell.
+        #[derive(Clone, Copy)]
+        struct Key<'a, E: Exponent> {
+            index: usize,
+            monomials: &'a UnsafeCell<(usize, Vec<E>)>,
+        }
+
+        impl<'a, E: Exponent> PartialEq for Key<'a, E> {
+            #[inline(always)]
+            fn eq(&self, other: &Self) -> bool {
+                unsafe {
+                    let b1 = &*self.monomials.get();
+                    b1.1.get_unchecked(self.index..self.index + b1.0)
+                        == b1.1.get_unchecked(other.index..other.index + b1.0)
+                }
+            }
+        }
+
+        impl<'a, E: Exponent> Eq for Key<'a, E> {}
+
+        impl<'a, E: Exponent> PartialOrd for Key<'a, E> {
+            #[inline(always)]
+            fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+                Some(self.cmp(other))
+            }
+        }
+
+        impl<'a, E: Exponent> Ord for Key<'a, E> {
+            #[inline(always)]
+            fn cmp(&self, other: &Self) -> Ordering {
+                unsafe {
+                    let b1 = &*self.monomials.get();
+                    b1.1.get_unchecked(self.index..self.index + b1.0)
+                        .cmp(&b1.1.get_unchecked(other.index..other.index + b1.0))
+                }
+            }
+        }
+
+        impl<'a, E: Exponent> std::hash::Hash for Key<'_, E> {
+            #[inline(always)]
+            fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
+                unsafe {
+                    let b = &*self.monomials.get();
+                    b.1.get_unchecked(self.index..self.index + b.0).hash(state);
+                }
+            }
+        }
+
         let mut res = self.zero_with_capacity(self.nterms().max(rhs.nterms()));
 
-        let mut cache: BTreeMap<Vec<E>, Vec<(usize, usize)>> = BTreeMap::new();
+        let mut cache: HashMap<_, Vec<(usize, usize)>> = HashMap::new();
         let mut q_cache: Vec<Vec<(usize, usize)>> = vec![];
 
         // create a min-heap since our polynomials are sorted smallest to largest
-        let mut h: BinaryHeap<Reverse<Vec<E>>> = BinaryHeap::with_capacity(self.nterms());
-
-        let monom: Vec<E> = self
-            .exponents(0)
-            .iter()
-            .zip(rhs.exponents(0))
-            .map(|(e1, e2)| *e1 + *e2)
-            .collect();
-        cache.insert(monom.clone(), vec![(0, 0)]);
-        h.push(Reverse(monom));
-
-        let mut m_cache: Vec<E> = vec![E::zero(); self.nvars()];
+        let mut h: BinaryHeap<Reverse<_>> = BinaryHeap::with_capacity(self.nterms());
+
+        cache.insert(
+            Key {
+                index: 0,
+                monomials: &monomials,
+            },
+            vec![(0, 0)],
+        );
+        h.push(Reverse(Key {
+            index: 0,
+            monomials: &monomials,
+        }));
 
         // i=merged_index[j] signifies that self[i]*other[j] has been merged
         let mut merged_index = vec![0; rhs.nterms()];
@@ -2205,23 +2444,31 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
                 merged_index[j] = i + 1;
 
                 if i + 1 < self.nterms() && (j == 0 || merged_index[j - 1] > i + 1) {
-                    for ((m, e1), e2) in m_cache
-                        .iter_mut()
-                        .zip(self.exponents(i + 1))
-                        .zip(rhs.exponents(j))
-                    {
-                        *m = *e1 + *e2;
-                    }
+                    let m = unsafe {
+                        let b = &mut *monomials.get();
+                        let index = b.1.len();
+                        b.1.extend(
+                            self.exponents(i + 1)
+                                .iter()
+                                .zip(rhs.exponents(j))
+                                .map(|(e1, e2)| *e1 + *e2),
+                        );
 
-                    if let Some(e) = cache.get_mut(&m_cache) {
+                        Key {
+                            index,
+                            monomials: &monomials,
+                        }
+                    };
+
+                    if let Some(e) = cache.get_mut(&m) {
                         e.push((i + 1, j));
                     } else {
-                        h.push(Reverse(m_cache.clone())); // only add when new
+                        h.push(Reverse(m)); // only add when new
                         if let Some(mut qq) = q_cache.pop() {
                             qq.push((i + 1, j));
-                            cache.insert(m_cache.clone(), qq);
+                            cache.insert(m, qq);
                         } else {
-                            cache.insert(m_cache.clone(), vec![(i + 1, j)]);
+                            cache.insert(m, vec![(i + 1, j)]);
                         }
                     }
                 } else {
@@ -2229,24 +2476,32 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
                 }
 
                 if j + 1 < rhs.nterms() && !in_heap[j + 1] {
-                    for ((m, e1), e2) in m_cache
-                        .iter_mut()
-                        .zip(self.exponents(i))
-                        .zip(rhs.exponents(j + 1))
-                    {
-                        *m = *e1 + *e2;
-                    }
+                    let m = unsafe {
+                        let b = &mut *monomials.get();
+                        let index = b.1.len();
+                        b.1.extend(
+                            self.exponents(i)
+                                .iter()
+                                .zip(rhs.exponents(j + 1))
+                                .map(|(e1, e2)| *e1 + *e2),
+                        );
 
-                    if let Some(e) = cache.get_mut(&m_cache) {
+                        Key {
+                            index,
+                            monomials: &monomials,
+                        }
+                    };
+
+                    if let Some(e) = cache.get_mut(&m) {
                         e.push((i, j + 1));
                     } else {
-                        h.push(Reverse(m_cache.clone())); // only add when new
+                        h.push(Reverse(m)); // only add when new
 
                         if let Some(mut qq) = q_cache.pop() {
                             qq.push((i, j + 1));
-                            cache.insert(m_cache.clone(), qq);
+                            cache.insert(m, qq);
                         } else {
-                            cache.insert(m_cache.clone(), vec![(i, j + 1)]);
+                            cache.insert(m, vec![(i, j + 1)]);
                         }
                     }
 
@@ -2258,9 +2513,15 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
 
             if !F::is_zero(&coefficient) {
                 res.coefficients.push(coefficient);
-                res.exponents.extend_from_slice(&cur_mon.0);
+
+                unsafe {
+                    let b = &*monomials.get();
+                    res.exponents
+                        .extend_from_slice(&b.1[cur_mon.0.index..cur_mon.0.index + b.0]);
+                }
             }
         }
+
         res
     }
 
@@ -2368,147 +2629,9 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         }
         res
     }
-
-    /// Synthetic division for univariate polynomials, where `div` is monic.
-    // TODO: create UnivariatePolynomial?
-    pub fn quot_rem_univariate_monic(
-        &self,
-        div: &MultivariatePolynomial<F, E, LexOrder>,
-    ) -> (
-        MultivariatePolynomial<F, E, LexOrder>,
-        MultivariatePolynomial<F, E, LexOrder>,
-    ) {
-        debug_assert_eq!(div.lcoeff(), self.ring.one());
-        if self.is_zero() {
-            return (self.clone(), self.clone());
-        }
-
-        let mut dividendpos = self.nterms() - 1; // work from the back
-
-        let mut q = self.zero_with_capacity(self.nterms());
-        let mut r = self.zero();
-
-        // determine the variable
-        let mut var = 0;
-        for (i, x) in self.last_exponents().iter().enumerate() {
-            if !x.is_zero() {
-                var = i;
-                break;
-            }
-        }
-
-        let m = div.ldegree_max();
-        let mut pow = self.ldegree_max();
-
-        loop {
-            // find the power in the dividend if it exists
-            let mut coeff = loop {
-                if self.exponents(dividendpos)[var] == pow {
-                    break self.coefficients[dividendpos].clone();
-                }
-                if dividendpos == 0 || self.exponents(dividendpos)[var] < pow {
-                    break self.ring.zero();
-                }
-                dividendpos -= 1;
-            };
-
-            let mut qindex = 0; // starting from highest
-            let mut bindex = 0; // starting from lowest
-            while bindex < div.nterms() && qindex < q.nterms() {
-                while bindex + 1 < div.nterms()
-                    && div.exponents(bindex)[var] + q.exponents(qindex)[var] < pow
-                {
-                    bindex += 1;
-                }
-
-                if div.exponents(bindex)[var] + q.exponents(qindex)[var] == pow {
-                    self.ring.sub_mul_assign(
-                        &mut coeff,
-                        &div.coefficients[bindex],
-                        &q.coefficients[qindex],
-                    );
-                }
-
-                qindex += 1;
-            }
-
-            if !F::is_zero(&coeff) {
-                // can the division be performed? if not, add to rest
-                // TODO: refactor
-                let (quot, div) = if pow >= m {
-                    (coeff, true)
-                } else {
-                    (coeff, false)
-                };
-
-                if div {
-                    let nterms = q.nterms();
-                    let nvars = q.nvars();
-                    q.coefficients.push(quot);
-                    q.exponents.resize((nterms + 1) * nvars, E::zero());
-                    q.exponents[nterms * nvars + var] = pow - m;
-                } else {
-                    let nterms = r.nterms();
-                    let nvars = r.nvars();
-                    r.coefficients.push(quot);
-                    r.exponents.resize((nterms + 1) * nvars, E::zero());
-                    r.exponents[nterms * nvars + var] = pow;
-                }
-            }
-
-            if pow.is_zero() {
-                break;
-            }
-
-            pow = pow - E::one();
-        }
-
-        q.reverse();
-        r.reverse();
-
-        #[cfg(debug_assertions)]
-        {
-            if !(&q * div + r.clone() - self.clone()).is_zero() {
-                panic!("Division failed: ({})/({}): q={}, r={}", self, div, q, r);
-            }
-        }
-
-        (q, r)
-    }
-
-    /// Shift a variable `var` to `var+shift`.
-    pub fn shift_var(&self, var: usize, shift: &F::Element) -> Self {
-        let d = self.degree(var).to_u32() as usize;
-
-        let y_poly = self.to_univariate_polynomial_list(var);
-
-        let mut v = vec![self.zero(); d + 1];
-        for (x_poly, p) in y_poly {
-            v[p.to_u32() as usize] = x_poly;
-        }
-
-        for k in 0..d {
-            for j in (k..d).rev() {
-                v[j] = &v[j] + &v[j + 1].clone().mul_coeff(shift.clone());
-            }
-        }
-
-        let mut poly = self.zero();
-        for (i, mut v) in v.into_iter().enumerate() {
-            for x in v.exponents.chunks_mut(self.nvars()) {
-                x[var] = E::from_u32(i as u32);
-            }
-
-            for m in &v {
-                poly.append_monomial(m.coefficient.clone(), m.exponents);
-            }
-        }
-
-        poly
-    }
 }
 
-impl<F: EuclideanDomain, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
+impl<F: EuclideanDomain, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
     /// Get the content from the coefficients.
     pub fn content(&self) -> F::Element {
         if self.coefficients.is_empty() {
@@ -2541,7 +2664,9 @@ impl<F: EuclideanDomain, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
         let c = self.content();
         self.div_coeff(&c)
     }
+}
 
+impl<F: EuclideanDomain, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
     pub fn divides(
         &self,
         div: &MultivariatePolynomial<F, E, LexOrder>,
@@ -3254,7 +3379,9 @@ impl<F: Field, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
             self
         }
     }
+}
 
+impl<F: Field, E: PositiveExponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
     /// Integrate the polynomial w.r.t the variable `var`,
     /// producing the antiderivative with zero constant.
     pub fn integrate(&self, var: usize) -> Self {
@@ -3277,7 +3404,7 @@ impl<F: Field, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
     }
 }
 
-impl<F: Field, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
+impl<F: Field, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
     /// Optimized division routine for univariate polynomials over a field, which
     /// makes the divisor monic first.
     pub fn quot_rem_univariate(
@@ -3502,7 +3629,7 @@ impl<F: Field, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
     }
 }
 
-impl<R: Ring, E: Exponent> Derivable for PolynomialRing<R, E> {
+impl<R: Ring, E: PositiveExponent> Derivable for PolynomialRing<R, E> {
     fn derivative(
         &self,
         p: &MultivariatePolynomial<R, E>,
@@ -3542,7 +3669,7 @@ impl<R: Ring, E: Exponent> MultivariatePolynomial<AlgebraicExtension<R>, E> {
         for t in self {
             exp[..self.nvars()].copy_from_slice(t.exponents);
             for t2 in &t.coefficient.poly {
-                exp[var_index] = E::from_u32(t2.exponents[0].to_u32());
+                exp[var_index] = E::from_i32(t2.exponents[0].to_i32());
                 poly.append_monomial(t2.coefficient.clone(), &exp);
             }
         }
diff --git a/src/poly/series.rs b/src/poly/series.rs
index 6bbafea..f23f17d 100644
--- a/src/poly/series.rs
+++ b/src/poly/series.rs
@@ -12,7 +12,7 @@ use crate::{
         atom::AtomField,
         integer::Integer,
         rational::{Rational, Q},
-        EuclideanDomain, InternalOrdering, Ring,
+        EuclideanDomain, InternalOrdering, Ring, SelfRing,
     },
     printer::{PrintOptions, PrintState},
     state::State,
@@ -465,20 +465,40 @@ impl<F: Ring> Series<F> {
 
         self
     }
+}
+
+impl<F: Ring> SelfRing for Series<F> {
+    fn is_zero(&self) -> bool {
+        self.is_zero()
+    }
 
-    pub fn format<W: std::fmt::Write>(
+    fn is_one(&self) -> bool {
+        self.is_one()
+    }
+
+    fn format<W: std::fmt::Write>(
         &self,
         opts: &PrintOptions,
         mut state: PrintState,
         f: &mut W,
     ) -> Result<bool, std::fmt::Error> {
-        let v = self.variable.to_string_with_state(PrintState {
-            in_exp: true,
-            ..state
-        });
+        let v = self.variable.format_string(
+            opts,
+            PrintState {
+                in_exp: true,
+                ..state
+            },
+        );
 
         if self.coefficients.is_empty() {
-            write!(f, "𝒪({}^{})", v, self.absolute_order())?;
+            let o = self.absolute_order();
+            if opts.latex {
+                write!(f, "\\mathcal{{O}}\\left({}^{{{}}}\\right)", v, o)?;
+            } else {
+                write!(f, "𝒪({}^", v)?;
+                Q.format(&o, opts, state.step(false, false, true), f)?;
+                f.write_char(')')?;
+            }
             return Ok(false);
         }
 
@@ -493,6 +513,7 @@ impl<F: Ring> Series<F> {
             state.in_exp = false;
             f.write_str("(")?;
         }
+        let in_product = state.in_product;
 
         for (e, c) in self.coefficients.iter().enumerate() {
             if F::is_zero(c) {
@@ -501,6 +522,7 @@ impl<F: Ring> Series<F> {
 
             let e = self.get_exponent(e);
 
+            state.in_product = in_product || !e.is_zero();
             state.suppress_one = !e.is_zero();
             let suppressed_one = self.field.format(
                 c,
@@ -517,18 +539,30 @@ impl<F: Ring> Series<F> {
                 write!(f, "{}", v)?;
             } else if !e.is_zero() {
                 write!(f, "{}^", v)?;
+                state.suppress_one = false;
+
+                if opts.latex {
+                    f.write_char('{')?;
+                }
+
                 Q.format(&e, opts, state.step(false, false, true), f)?;
+
+                if opts.latex {
+                    f.write_char('}')?;
+                }
             }
 
             state.in_sum = true;
-            state.in_product = true;
         }
 
         let o = self.absolute_order();
-        if o.is_integer() {
-            write!(f, "+𝒪({}^{})", v, o)?;
+
+        if opts.latex {
+            write!(f, "+\\mathcal{{O}}\\left({}^{{{}}}\\right)", v, o)?;
         } else {
-            write!(f, "+𝒪({}^({}))", v, o)?;
+            write!(f, "+𝒪({}^", v)?;
+            Q.format(&o, opts, state.step(false, false, true), f)?;
+            f.write_char(')')?;
         }
 
         if add_paren {
@@ -539,6 +573,28 @@ impl<F: Ring> Series<F> {
     }
 }
 
+impl<F: Ring> InternalOrdering for Series<F> {
+    fn internal_cmp(&self, other: &Self) -> Ordering {
+        if self.variable != other.variable {
+            return self.variable.cmp(&other.variable);
+        }
+
+        if self.shift != other.shift {
+            return self.shift.cmp(&other.shift);
+        }
+
+        if self.ramification != other.ramification {
+            return self.ramification.cmp(&other.ramification);
+        }
+
+        if self.order != other.order {
+            return self.order.cmp(&other.order);
+        }
+
+        self.coefficients.internal_cmp(&other.coefficients)
+    }
+}
+
 impl<F: Ring> PartialEq for Series<F> {
     #[inline]
     fn eq(&self, other: &Self) -> bool {
diff --git a/src/poly/univariate.rs b/src/poly/univariate.rs
index aabd824..fe49da2 100644
--- a/src/poly/univariate.rs
+++ b/src/poly/univariate.rs
@@ -17,7 +17,7 @@ use crate::{
 use super::{
     factor::Factorize,
     polynomial::{MultivariatePolynomial, PolynomialRing},
-    Exponent, Variable,
+    PositiveExponent, Variable,
 };
 
 #[derive(Clone, PartialEq, Eq, Hash, Debug)]
@@ -457,7 +457,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     }
 
     /// Convert from a univariate polynomial to a multivariate polynomial.
-    pub fn to_multivariate<E: Exponent>(self) -> MultivariatePolynomial<F, E> {
+    pub fn to_multivariate<E: PositiveExponent>(self) -> MultivariatePolynomial<F, E> {
         let mut res = MultivariatePolynomial::new(
             &self.ring,
             self.degree().into(),
@@ -538,10 +538,13 @@ impl<F: Ring> SelfRing for UnivariatePolynomial<F> {
             f.write_str("(")?;
         }
 
-        let v = self.variable.to_string_with_state(PrintState {
-            in_exp: true,
-            ..state
-        });
+        let v = self.variable.format_string(
+            opts,
+            PrintState {
+                in_exp: true,
+                ..state
+            },
+        );
 
         for (e, c) in self.coefficients.iter().enumerate() {
             state.suppress_one = e > 0;
@@ -1730,7 +1733,7 @@ impl<F: Field> UnivariatePolynomial<F> {
     }
 }
 
-impl<R: Ring, E: Exponent> UnivariatePolynomial<PolynomialRing<R, E>> {
+impl<R: Ring, E: PositiveExponent> UnivariatePolynomial<PolynomialRing<R, E>> {
     /// Convert a univariate polynomial of multivariate polynomials to a multivariate polynomial.
     pub fn flatten(self) -> MultivariatePolynomial<R, E> {
         if self.is_zero() {
diff --git a/src/solve.rs b/src/solve.rs
index 9b97f81..d31850d 100644
--- a/src/solve.rs
+++ b/src/solve.rs
@@ -10,7 +10,7 @@ use crate::{
         InternalOrdering,
     },
     evaluate::FunctionMap,
-    poly::{Exponent, Variable},
+    poly::{PositiveExponent, Variable},
     tensors::matrix::Matrix,
 };
 
@@ -41,7 +41,7 @@ impl Atom {
 
     /// Solve a system that is linear in `vars`, if possible.
     /// Each expression in `system` is understood to yield 0.
-    pub fn solve_linear_system<E: Exponent>(
+    pub fn solve_linear_system<E: PositiveExponent>(
         system: &[AtomView],
         vars: &[Symbol],
     ) -> Result<Vec<Atom>, String> {
@@ -189,7 +189,7 @@ impl<'a> AtomView<'a> {
 
     /// Solve a system that is linear in `vars`, if possible.
     /// Each expression in `system` is understood to yield 0.
-    pub fn solve_linear_system<E: Exponent>(
+    pub fn solve_linear_system<E: PositiveExponent>(
         system: &[AtomView],
         vars: &[Symbol],
     ) -> Result<Vec<Atom>, String> {
diff --git a/symbolica.pyi b/symbolica.pyi
index 4d1bc54..4472867 100644
--- a/symbolica.pyi
+++ b/symbolica.pyi
@@ -2057,6 +2057,32 @@ class Series:
     >>> print(s)
     """
 
+    def __str__(self) -> str:
+        """Print the series in a human-readable format."""
+
+    def to_latex(self) -> str:
+        """Convert the series into a LaTeX string."""
+
+    def format(
+        self,
+        terms_on_new_line: bool = False,
+        color_top_level_sum: bool = True,
+        color_builtin_symbols: bool = True,
+        print_finite_field: bool = True,
+        symmetric_representation_for_finite_field: bool = False,
+        explicit_rational_polynomial: bool = False,
+        number_thousands_separator: Optional[str] = None,
+        multiplication_operator: str = "*",
+        double_star_for_exponentiation: bool = False,
+        square_brackets_for_function: bool = False,
+        num_exp_as_superscript: bool = True,
+        latex: bool = False,
+        precision: Optional[int] = None,
+    ) -> str:
+        """
+        Convert the series into a human-readable string.
+        """
+
     def __add__(self, other: Series | Expression) -> Series:
         """Add another series or expression to this series, returning the result."""