diff --git a/src/api/python.rs b/src/api/python.rs
index 31d9e1a..829bb76 100644
--- a/src/api/python.rs
+++ b/src/api/python.rs
@@ -2345,7 +2345,7 @@ impl PythonExpression {
     /// 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.pretty_str(number_thousands_separator='_', multiplication_operator=' '))
+    /// >>> print(a.format(number_thousands_separator='_', multiplication_operator=' '))
     #[pyo3(signature =
         (terms_on_new_line = false,
             color_top_level_sum = true,
@@ -2361,7 +2361,7 @@ impl PythonExpression {
             latex = false,
             precision = None)
         )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
@@ -5354,7 +5354,7 @@ impl PythonPolynomial {
     /// Examples
     /// --------
     /// >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-    /// >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+    /// >>> print(p.format(symmetric_representation_for_finite_field=True))
     #[pyo3(signature =
         (terms_on_new_line = false,
             color_top_level_sum = true,
@@ -5370,7 +5370,7 @@ impl PythonPolynomial {
             latex = false,
             precision = None)
         )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
@@ -6177,7 +6177,7 @@ impl PythonFiniteFieldPolynomial {
     /// Examples
     /// --------
     /// >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-    /// >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+    /// >>> print(p.format(symmetric_representation_for_finite_field=True))
     #[pyo3(signature =
         (terms_on_new_line = false,
             color_top_level_sum = true,
@@ -6193,7 +6193,7 @@ impl PythonFiniteFieldPolynomial {
             latex = false,
             precision = None)
         )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
@@ -6757,7 +6757,7 @@ impl PythonPrimeTwoPolynomial {
     /// Examples
     /// --------
     /// >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-    /// >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+    /// >>> print(p.format(symmetric_representation_for_finite_field=True))
     #[pyo3(signature =
         (terms_on_new_line = false,
             color_top_level_sum = true,
@@ -6773,7 +6773,7 @@ impl PythonPrimeTwoPolynomial {
             latex = false,
             precision = None)
         )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
@@ -7297,7 +7297,7 @@ impl PythonGaloisFieldPrimeTwoPolynomial {
     /// Examples
     /// --------
     /// >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-    /// >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+    /// >>> print(p.format(symmetric_representation_for_finite_field=True))
     #[pyo3(signature =
     (terms_on_new_line = false,
         color_top_level_sum = true,
@@ -7313,7 +7313,7 @@ impl PythonGaloisFieldPrimeTwoPolynomial {
         latex = false,
             precision = None)
     )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
@@ -7841,7 +7841,7 @@ impl PythonGaloisFieldPolynomial {
     /// Examples
     /// --------
     /// >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-    /// >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+    /// >>> print(p.format(symmetric_representation_for_finite_field=True))
     #[pyo3(signature =
         (terms_on_new_line = false,
             color_top_level_sum = true,
@@ -7857,7 +7857,7 @@ impl PythonGaloisFieldPolynomial {
             latex = false,
             precision = None)
         )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
@@ -8386,7 +8386,7 @@ impl PythonNumberFieldPolynomial {
     /// Examples
     /// --------
     /// >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-    /// >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+    /// >>> print(p.format(symmetric_representation_for_finite_field=True))
     #[pyo3(signature =
     (terms_on_new_line = false,
         color_top_level_sum = true,
@@ -8402,7 +8402,7 @@ impl PythonNumberFieldPolynomial {
         latex = false,
             precision = None)
     )]
-    pub fn pretty_str(
+    pub fn format(
         &self,
         terms_on_new_line: bool,
         color_top_level_sum: bool,
diff --git a/src/domains/factorized_rational_polynomial.rs b/src/domains/factorized_rational_polynomial.rs
index 090ed8e..5a4e04a 100644
--- a/src/domains/factorized_rational_polynomial.rs
+++ b/src/domains/factorized_rational_polynomial.rs
@@ -246,7 +246,7 @@ impl<R: Ring, E: Exponent> SelfRing for FactorizedRationalPolynomial<R, E> {
                 if self.numerator.is_one() {
                     return Ok(false);
                 }
-                f.write_char('*')?;
+                f.write_char(opts.multiplication_operator)?;
             }
 
             if write_par {
@@ -311,7 +311,7 @@ impl<R: Ring, E: Exponent> SelfRing for FactorizedRationalPolynomial<R, E> {
             return Ok(false);
         }
 
-        state.in_product |= has_numer_coeff && !self.numerator.is_one();
+        state.in_product = true;
         if has_numer_coeff || self.numerator.is_one() {
             self.numerator
                 .ring
@@ -319,7 +319,7 @@ impl<R: Ring, E: Exponent> SelfRing for FactorizedRationalPolynomial<R, E> {
             state.in_sum = false;
 
             if !self.numerator.is_one() {
-                f.write_char('*')?;
+                f.write_char(opts.multiplication_operator)?;
             }
         }
 
@@ -355,7 +355,7 @@ impl<R: Ring, E: Exponent> SelfRing for FactorizedRationalPolynomial<R, E> {
 
         for (i, (d, p)) in self.denominators.iter().enumerate() {
             if has_denom_coeff || i > 0 {
-                f.write_char('*')?;
+                f.write_char(opts.multiplication_operator)?;
             }
 
             d.format(opts, state.step(false, true, *p != 1), f)?;
diff --git a/src/domains/rational_polynomial.rs b/src/domains/rational_polynomial.rs
index b6dac38..ee95ee1 100644
--- a/src/domains/rational_polynomial.rs
+++ b/src/domains/rational_polynomial.rs
@@ -513,8 +513,8 @@ where
     ) -> RationalPolynomial<R, E> {
         if f.is_zero() {
             return RationalPolynomial {
-                numerator: MultivariatePolynomial::new_zero(&f.field.ring),
-                denominator: MultivariatePolynomial::new_one(&f.field.ring),
+                numerator: MultivariatePolynomial::new_zero(&f.ring.ring),
+                denominator: MultivariatePolynomial::new_one(&f.ring.ring),
             };
         }
 
@@ -535,7 +535,7 @@ where
             });
 
         let mut res =
-            RationalPolynomial::new(&f.field.ring, f.coefficients[0].get_variables().clone());
+            RationalPolynomial::new(&f.ring.ring, f.coefficients[0].get_variables().clone());
 
         let mut exp = vec![E::zero(); f.coefficients[0].get_variables().len()];
         exp[pos] = E::one();
@@ -1103,7 +1103,7 @@ where
 
                 d_exp[i - 1] = d_exp[i - 1].clone()
                     + s_cor
-                    + t_cor.derivative().div_coeff(&(t_cor.field.nth(i as u64)));
+                    + t_cor.derivative().div_coeff(&(t_cor.ring.nth(i as u64)));
 
                 let t_full = Self::from_univariate(t_cor);
 
@@ -1186,7 +1186,7 @@ where
                     .to_multivariate_polynomial_list(&[var, new_var], true);
 
                 let mut bivar_poly = MultivariatePolynomial::new(
-                    &p.field,
+                    &p.ring,
                     Some(ll.len()),
                     p.coefficients[0].get_variables().clone(),
                 );
@@ -1197,7 +1197,7 @@ where
                 // convert defining polynomial to a univariate polynomial in t with rational polynomial coefficients
                 let def_uni = sqf
                     .to_univariate(new_var)
-                    .map_coeff(|c| c.clone().into(), p.field.clone());
+                    .map_coeff(|c| c.clone().into(), p.ring.clone());
 
                 // write the polynomial in x and t as a polynomial in x with rational polynomial coefficients in t and
                 // all other variables and solve a diophantine equation
@@ -1227,7 +1227,7 @@ where
                 let monic = bivar_poly_scaled.rem(&def_biv);
 
                 // convert the result to a multivariate rational polynomial
-                let mut res = p.field.zero();
+                let mut res = p.ring.zero();
                 for t in &monic {
                     let mut exp = vec![E::zero(); p.lcoeff().numerator.nvars()];
                     exp.copy_from_slice(t.exponents);
@@ -1250,7 +1250,7 @@ where
 
                         let eval = monic.replace(new_var, &sol);
 
-                        let mut res = p.field.zero();
+                        let mut res = p.ring.zero();
                         for t in &eval {
                             let mut exp = vec![E::zero(); p.lcoeff().numerator.nvars()];
                             exp.copy_from_slice(t.exponents);
diff --git a/src/poly/polynomial.rs b/src/poly/polynomial.rs
index c4afe34..275a6b3 100755
--- a/src/poly/polynomial.rs
+++ b/src/poly/polynomial.rs
@@ -741,7 +741,8 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> SelfRing for MultivariatePolynomial
         let add_paren = self.nterms() > 1 && state.in_product
             || (state.in_exp
                 && (self.nterms() > 1
-                    || self.exponents(0).iter().filter(|e| **e > E::zero()).count() > 1));
+                    || self.exponents(0).iter().filter(|e| **e > E::zero()).count() > 1
+                    || !self.ring.is_one(&self.coefficients[0])));
 
         if add_paren {
             if state.in_sum {
@@ -781,7 +782,7 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> SelfRing for MultivariatePolynomial
                 if suppressed_one {
                     suppressed_one = false;
                 } else if !opts.latex {
-                    write!(f, "*")?;
+                    f.write_char(opts.multiplication_operator)?;
                 }
 
                 f.write_str(var_id)?;
@@ -1789,7 +1790,7 @@ impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
             return p;
         }
 
-        p.coefficients = vec![p.field.zero(); self.degree(var).to_u32() as usize + 1];
+        p.coefficients = vec![p.ring.zero(); self.degree(var).to_u32() as usize + 1];
         for (q, e) in self.coefficients.iter().zip(self.exponents_iter()) {
             p.coefficients[e[var].to_u32() as usize] = q.clone();
         }
diff --git a/src/poly/resultant.rs b/src/poly/resultant.rs
index 993aef7..72cb574 100644
--- a/src/poly/resultant.rs
+++ b/src/poly/resultant.rs
@@ -7,24 +7,24 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
     pub fn resultant_prs(&self, other: &Self) -> F::Element {
         if self.degree() < other.degree() {
             if self.degree() % 2 == 1 && other.degree() % 2 == 1 {
-                return self.field.neg(&other.resultant_prs(self));
+                return self.ring.neg(&other.resultant_prs(self));
             } else {
                 return other.resultant_prs(self);
             }
         }
 
         if other.is_constant() {
-            return self.field.pow(&other.get_constant(), self.degree() as u64);
+            return self.ring.pow(&other.get_constant(), self.degree() as u64);
         }
 
         let mut a = self.clone();
         let mut a_new = other.clone();
 
         let mut deg = a.degree() as u64 - a_new.degree() as u64;
-        let mut neg_lc = self.field.one(); //unused
+        let mut neg_lc = self.ring.one(); //unused
         let mut init = false;
-        let mut beta = self.field.pow(&self.field.neg(&self.field.one()), deg + 1);
-        let mut psi = self.field.neg(&self.field.one());
+        let mut beta = self.ring.pow(&self.ring.neg(&self.ring.one()), deg + 1);
+        let mut psi = self.ring.neg(&self.ring.one());
 
         let mut lcs = vec![(a.lcoeff(), a.degree() as u64)];
         while !a_new.is_constant() {
@@ -35,23 +35,23 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
                 } else if deg == 1 {
                     neg_lc.clone()
                 } else {
-                    let a = self.field.pow(&neg_lc, deg);
-                    let psi_old = self.field.pow(&psi, deg - 1);
-                    let (q, r) = self.field.quot_rem(&a, &psi_old);
+                    let a = self.ring.pow(&neg_lc, deg);
+                    let psi_old = self.ring.pow(&psi, deg - 1);
+                    let (q, r) = self.ring.quot_rem(&a, &psi_old);
                     debug_assert!(F::is_zero(&r));
                     q
                 };
 
                 deg = a.degree() as u64 - a_new.degree() as u64;
-                beta = self.field.mul(&neg_lc, &self.field.pow(&psi, deg));
+                beta = self.ring.mul(&neg_lc, &self.ring.pow(&psi, deg));
             } else {
                 init = true;
             }
 
-            neg_lc = self.field.neg(a_new.coefficients.last().unwrap());
+            neg_lc = self.ring.neg(a_new.coefficients.last().unwrap());
 
             let (_, mut r) = a
-                .mul_coeff(&self.field.pow(&neg_lc, deg + 1))
+                .mul_coeff(&self.ring.pow(&neg_lc, deg + 1))
                 .quot_rem(&a_new);
             if (deg + 1) % 2 == 1 {
                 r = -r;
@@ -65,12 +65,12 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
         lcs.push((a_new.lcoeff(), 0));
 
         if a_new.is_zero() {
-            return self.field.zero();
+            return self.ring.zero();
         }
 
         // compute the resultant from the PRS, using the fundamental theorem
-        let mut rho = self.field.one();
-        let mut den = self.field.one();
+        let mut rho = self.ring.one();
+        let mut den = self.ring.one();
         for k in 1..lcs.len() {
             let mut deg = lcs[k as usize - 1].1 as i64 - lcs[k as usize].1 as i64;
             for l in k..lcs.len() - 1 {
@@ -78,22 +78,22 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
             }
 
             if deg > 0 {
-                self.field
-                    .mul_assign(&mut rho, &self.field.pow(&lcs[k].0, deg as u64));
+                self.ring
+                    .mul_assign(&mut rho, &self.ring.pow(&lcs[k].0, deg as u64));
             } else if deg < 0 {
-                self.field
-                    .mul_assign(&mut den, &self.field.pow(&lcs[k].0, (-deg) as u64));
+                self.ring
+                    .mul_assign(&mut den, &self.ring.pow(&lcs[k].0, (-deg) as u64));
             }
         }
 
-        self.field.quot_rem(&rho, &den).0
+        self.ring.quot_rem(&rho, &den).0
     }
 
     /// Compute the resultant using a primitive polynomial remainder sequence.
     pub fn resultant_primitive(&self, other: &Self) -> F::Element {
         if self.degree() < other.degree() {
             if self.degree() % 2 == 1 && other.degree() % 2 == 1 {
-                return self.field.neg(&other.resultant_primitive(self));
+                return self.ring.neg(&other.resultant_primitive(self));
             } else {
                 return other.resultant_primitive(self);
             }
@@ -104,11 +104,11 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
 
         let mut v = vec![a.degree()];
         let mut c = vec![a.lcoeff()];
-        let mut ab = vec![(self.field.one(), self.field.one())];
+        let mut ab = vec![(self.ring.one(), self.ring.one())];
 
         while !a_new.is_constant() {
             let n = a.degree() as u64 + 1 - a_new.degree() as u64;
-            let alpha = self.field.pow(&a_new.lcoeff(), n);
+            let alpha = self.ring.pow(&a_new.lcoeff(), n);
 
             let (_, mut r) = a.clone().mul_coeff(&alpha).quot_rem(&a_new);
 
@@ -132,27 +132,27 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
             sign += w[0] * w[1];
         }
 
-        let mut res = self.field.pow(&r, *v.last().unwrap() as u64);
+        let mut res = self.ring.pow(&r, *v.last().unwrap() as u64);
         if sign % 2 == 1 {
-            res = self.field.neg(&res);
+            res = self.ring.neg(&res);
         };
 
         v.push(0);
-        let mut num = self.field.one();
-        let mut den = self.field.one();
+        let mut num = self.ring.one();
+        let mut den = self.ring.one();
         for i in 1..c.len() {
-            self.field.mul_assign(
+            self.ring.mul_assign(
                 &mut res,
-                &self.field.pow(&c[i], v[i - 1] as u64 - v[i + 1] as u64),
+                &self.ring.pow(&c[i], v[i - 1] as u64 - v[i + 1] as u64),
             );
 
-            self.field
-                .mul_assign(&mut num, &self.field.pow(&ab[i].1, v[i] as u64));
-            self.field
-                .mul_assign(&mut den, &self.field.pow(&ab[i].0, v[i] as u64));
+            self.ring
+                .mul_assign(&mut num, &self.ring.pow(&ab[i].1, v[i] as u64));
+            self.ring
+                .mul_assign(&mut den, &self.ring.pow(&ab[i].0, v[i] as u64));
         }
 
-        let (q, r) = self.field.quot_rem(&self.field.mul(&num, &res), &den);
+        let (q, r) = self.ring.quot_rem(&self.ring.mul(&num, &res), &den);
         assert!(F::is_zero(&r));
         q
     }
@@ -185,16 +185,16 @@ impl<F: Field> UnivariatePolynomial<F> {
             sign += w[0] * w[1];
         }
 
-        let mut res = self.field.pow(&r, *v.last().unwrap() as u64);
+        let mut res = self.ring.pow(&r, *v.last().unwrap() as u64);
         if sign % 2 == 1 {
-            res = self.field.neg(&res);
+            res = self.ring.neg(&res);
         };
 
         v.push(0);
         for i in 1..c.len() {
-            self.field.mul_assign(
+            self.ring.mul_assign(
                 &mut res,
-                &self.field.pow(&c[i], v[i - 1] as u64 - v[i + 1] as u64),
+                &self.ring.pow(&c[i], v[i - 1] as u64 - v[i + 1] as u64),
             );
         }
 
diff --git a/src/poly/univariate.rs b/src/poly/univariate.rs
index 3a72aaa..aabd824 100644
--- a/src/poly/univariate.rs
+++ b/src/poly/univariate.rs
@@ -36,7 +36,7 @@ impl<R: Ring> UnivariatePolynomialRing<R> {
 
     pub fn new_from_poly(poly: &UnivariatePolynomial<R>) -> UnivariatePolynomialRing<R> {
         UnivariatePolynomialRing {
-            ring: poly.field.clone(),
+            ring: poly.ring.clone(),
             variable: poly.variable.clone(),
         }
     }
@@ -157,7 +157,7 @@ impl<R: EuclideanDomain> EuclideanDomain for UnivariatePolynomialRing<R> {
 pub struct UnivariatePolynomial<F: Ring> {
     pub coefficients: Vec<F::Element>,
     pub variable: Arc<Variable>,
-    pub field: F,
+    pub ring: F,
 }
 
 impl<R: Ring> InternalOrdering for UnivariatePolynomial<R> {
@@ -201,7 +201,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     pub fn new(field: &F, cap: Option<usize>, variable: Arc<Variable>) -> Self {
         Self {
             coefficients: Vec::with_capacity(cap.unwrap_or(0)),
-            field: field.clone(),
+            ring: field.clone(),
             variable,
         }
     }
@@ -211,7 +211,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     pub fn zero(&self) -> Self {
         Self {
             coefficients: vec![],
-            field: self.field.clone(),
+            ring: self.ring.clone(),
             variable: self.variable.clone(),
         }
     }
@@ -222,7 +222,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     pub fn zero_with_capacity(&self, cap: usize) -> Self {
         Self {
             coefficients: Vec::with_capacity(cap),
-            field: self.field.clone(),
+            ring: self.ring.clone(),
             variable: self.variable.clone(),
         }
     }
@@ -237,7 +237,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
 
         Self {
             coefficients: vec![coeff],
-            field: self.field.clone(),
+            ring: self.ring.clone(),
             variable: self.variable.clone(),
         }
     }
@@ -246,8 +246,8 @@ impl<F: Ring> UnivariatePolynomial<F> {
     #[inline]
     pub fn one(&self) -> Self {
         Self {
-            coefficients: vec![self.field.one()],
-            field: self.field.clone(),
+            coefficients: vec![self.ring.one()],
+            ring: self.ring.clone(),
             variable: self.variable.clone(),
         }
     }
@@ -259,12 +259,12 @@ impl<F: Ring> UnivariatePolynomial<F> {
             return self.zero();
         }
 
-        let mut coefficients = vec![self.field.zero(); exponent + 1];
+        let mut coefficients = vec![self.ring.zero(); exponent + 1];
         coefficients[exponent] = coeff;
 
         Self {
             coefficients,
-            field: self.field.clone(),
+            ring: self.ring.clone(),
             variable: self.variable.clone(),
         }
     }
@@ -276,7 +276,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
 
     #[inline]
     pub fn is_one(&self) -> bool {
-        self.coefficients.len() == 1 && self.field.is_one(&self.coefficients[0])
+        self.coefficients.len() == 1 && self.ring.is_one(&self.coefficients[0])
     }
 
     /// Returns true if the polynomial is constant.
@@ -289,7 +289,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     #[inline]
     pub fn get_constant(&self) -> F::Element {
         if self.is_zero() {
-            return self.field.zero();
+            return self.ring.zero();
         }
 
         self.coefficients[0].clone()
@@ -309,7 +309,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     pub fn lcoeff(&self) -> F::Element {
         self.coefficients
             .last()
-            .unwrap_or(&self.field.zero())
+            .unwrap_or(&self.ring.zero())
             .clone()
     }
 
@@ -348,7 +348,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
         }
 
         let mut a = self.zero();
-        a.coefficients = vec![self.field.zero(); self.degree() + exp + 1];
+        a.coefficients = vec![self.ring.zero(); self.degree() + exp + 1];
 
         for (cn, c) in a.coefficients.iter_mut().skip(exp).zip(&self.coefficients) {
             *cn = c.clone();
@@ -368,7 +368,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
             return a;
         }
 
-        a.coefficients = vec![self.field.zero(); self.degree() - exp + 1];
+        a.coefficients = vec![self.ring.zero(); self.degree() - exp + 1];
 
         for (cn, c) in a
             .coefficients
@@ -384,7 +384,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
     pub fn mul_coeff(mut self, coeff: &F::Element) -> Self {
         for c in &mut self.coefficients {
             if !F::is_zero(c) {
-                self.field.mul_assign(c, coeff);
+                self.ring.mul_assign(c, coeff);
             }
         }
 
@@ -423,9 +423,9 @@ impl<F: Ring> UnivariatePolynomial<F> {
         let mut res = self.coefficients.last().unwrap().clone();
         for c in self.coefficients.iter().rev().skip(1) {
             if !F::is_zero(c) {
-                res = self.field.add(&self.field.mul(&res, x), c);
+                res = self.ring.add(&self.ring.mul(&res, x), c);
             } else {
-                self.field.mul_assign(&mut res, x);
+                self.ring.mul_assign(&mut res, x);
             }
         }
 
@@ -440,7 +440,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
 
         let mut res = self.zero();
         res.coefficients
-            .resize(self.coefficients.len() - 1, self.field.zero());
+            .resize(self.coefficients.len() - 1, self.ring.zero());
 
         for (p, (nc, oc)) in res
             .coefficients
@@ -449,7 +449,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
             .enumerate()
         {
             if !F::is_zero(oc) {
-                *nc = self.field.mul(oc, &self.field.nth(p as u64 + 1));
+                *nc = self.ring.mul(oc, &self.ring.nth(p as u64 + 1));
             }
         }
 
@@ -459,7 +459,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> {
         let mut res = MultivariatePolynomial::new(
-            &self.field,
+            &self.ring,
             self.degree().into(),
             Arc::new(vec![self.variable.as_ref().clone()]),
         );
@@ -480,7 +480,7 @@ impl<F: Ring> UnivariatePolynomial<F> {
         for k in 0..d {
             for j in (k..d).rev() {
                 let (s, c) = poly.coefficients.split_at_mut(j + 1);
-                self.field.add_mul_assign(&mut s[j], &c[0], shift);
+                self.ring.add_mul_assign(&mut s[j], &c[0], shift);
             }
         }
 
@@ -511,14 +511,21 @@ impl<F: Ring> SelfRing for UnivariatePolynomial<F> {
                 f.write_char('0')?;
                 return Ok(false);
             } else {
-                return self.field.format(&self.coefficients[0], opts, state, f);
+                return self.ring.format(&self.coefficients[0], opts, state, f);
             }
         }
 
         let non_zero = self.coefficients.iter().filter(|c| !F::is_zero(c)).count();
 
         let add_paren = non_zero > 1 && state.in_product
-            || (state.in_exp && (non_zero > 1 || F::is_zero(&self.coefficients[0])));
+            || (state.in_exp
+                && (non_zero > 1
+                    || self
+                        .coefficients
+                        .iter()
+                        .filter(|c| !self.ring.is_one(c))
+                        .count()
+                        > 0));
 
         if add_paren {
             if state.in_sum {
@@ -543,7 +550,7 @@ impl<F: Ring> SelfRing for UnivariatePolynomial<F> {
                 continue;
             }
 
-            let suppressed_one = self.field.format(
+            let suppressed_one = self.ring.format(
                 c,
                 opts,
                 state.step(state.in_sum, state.in_product, false),
@@ -551,7 +558,7 @@ impl<F: Ring> SelfRing for UnivariatePolynomial<F> {
             )?;
 
             if !suppressed_one && e > 0 {
-                f.write_str("*")?;
+                f.write_char(opts.multiplication_operator)?;
             }
 
             if e == 1 {
@@ -585,7 +592,7 @@ impl UnivariatePolynomial<RationalField> {
 
         let stripped = self.map_coeff(
             |coeff| {
-                let coeff = self.field.div(coeff, &c);
+                let coeff = self.ring.div(coeff, &c);
                 debug_assert!(coeff.is_integer());
                 coeff.numerator()
             },
@@ -1059,7 +1066,7 @@ impl<R: Real + SingleFloat + std::hash::Hash + Eq + PartialOrd + InternalOrderin
     /// Get an upper bound on the norm of all (complex) roots.
     pub fn get_root_upper_bound(&self) -> R {
         if self.is_zero() {
-            return self.field.zero().re;
+            return self.ring.zero().re;
         }
 
         let last = self.coefficients.last().unwrap();
@@ -1071,13 +1078,13 @@ impl<R: Real + SingleFloat + std::hash::Hash + Eq + PartialOrd + InternalOrderin
             }
         }
 
-        max + self.field.one().re
+        max + self.ring.one().re
     }
 
     /// Get a lower bound on the norm of all (complex) roots.
     pub fn get_root_lower_bound(&self) -> R {
         if self.is_zero() {
-            return self.field.zero().re;
+            return self.ring.zero().re;
         }
 
         let last = &self.coefficients[0];
@@ -1089,7 +1096,7 @@ impl<R: Real + SingleFloat + std::hash::Hash + Eq + PartialOrd + InternalOrderin
             }
         }
 
-        self.field.one().re / (max + self.field.one().re)
+        self.ring.one().re / (max + self.ring.one().re)
     }
 
     /// Compute all complex roots of the polynomial using Aberth's method.
@@ -1105,11 +1112,11 @@ impl<R: Real + SingleFloat + std::hash::Hash + Eq + PartialOrd + InternalOrderin
         if self.get_constant().is_zero() {
             match self.div_exp(1).roots(max_iterations, tolerance) {
                 Ok(mut roots) => {
-                    roots.push(self.field.zero());
+                    roots.push(self.ring.zero());
                     return Ok(roots);
                 }
                 Err(mut roots) => {
-                    roots.push(self.field.zero());
+                    roots.push(self.ring.zero());
                     return Err(roots);
                 }
             }
@@ -1216,7 +1223,7 @@ impl<F: Ring> Add for UnivariatePolynomial<F> {
     type Output = Self;
 
     fn add(mut self, mut other: Self) -> Self::Output {
-        assert_eq!(self.field, other.field);
+        assert_eq!(self.ring, other.ring);
 
         if self.variable != other.variable {
             panic!("Cannot multiply polynomials with different variables");
@@ -1234,7 +1241,7 @@ impl<F: Ring> Add for UnivariatePolynomial<F> {
         }
 
         for (i, c) in other.coefficients.iter().enumerate() {
-            self.field.add_assign(&mut self.coefficients[i], c);
+            self.ring.add_assign(&mut self.coefficients[i], c);
         }
 
         self.truncate();
@@ -1272,7 +1279,7 @@ impl<F: Ring> Neg for UnivariatePolynomial<F> {
     fn neg(mut self) -> Self::Output {
         // Negate coefficients of all terms.
         for c in &mut self.coefficients {
-            *c = self.field.neg(c);
+            *c = self.ring.neg(c);
         }
         self
     }
@@ -1283,7 +1290,7 @@ impl<'a, 'b, F: Ring> Mul<&'a UnivariatePolynomial<F>> for &'b UnivariatePolynom
 
     #[inline]
     fn mul(self, rhs: &'a UnivariatePolynomial<F>) -> Self::Output {
-        assert_eq!(self.field, rhs.field);
+        assert_eq!(self.ring, rhs.ring);
 
         if self.is_zero() || rhs.is_zero() {
             return self.zero();
@@ -1299,7 +1306,7 @@ impl<'a, 'b, F: Ring> Mul<&'a UnivariatePolynomial<F>> for &'b UnivariatePolynom
         if n == 0 {
             let mut r = rhs.clone();
             for c in &mut r.coefficients {
-                self.field.mul_assign(c, &self.coefficients[0]);
+                self.ring.mul_assign(c, &self.coefficients[0]);
             }
             return r;
         }
@@ -1307,13 +1314,13 @@ impl<'a, 'b, F: Ring> Mul<&'a UnivariatePolynomial<F>> for &'b UnivariatePolynom
         if m == 0 {
             let mut r = self.clone();
             for c in &mut r.coefficients {
-                self.field.mul_assign(c, &rhs.coefficients[0]);
+                self.ring.mul_assign(c, &rhs.coefficients[0]);
             }
             return r;
         }
 
         let mut res = self.zero();
-        res.coefficients = vec![self.field.zero(); n + m + 1];
+        res.coefficients = vec![self.ring.zero(); n + m + 1];
 
         for (e1, c1) in self.coefficients.iter().enumerate() {
             if F::is_zero(c1) {
@@ -1322,7 +1329,7 @@ impl<'a, 'b, F: Ring> Mul<&'a UnivariatePolynomial<F>> for &'b UnivariatePolynom
 
             for (e2, c2) in rhs.coefficients.iter().enumerate() {
                 if !F::is_zero(c2) {
-                    self.field
+                    self.ring
                         .add_mul_assign(&mut res.coefficients[e1 + e2], c1, c2);
                 }
             }
@@ -1363,16 +1370,16 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
     /// Get the content from the coefficients.
     pub fn content(&self) -> F::Element {
         if self.coefficients.is_empty() {
-            return self.field.zero();
+            return self.ring.zero();
         }
         let mut c = self.coefficients.first().unwrap().clone();
         for cc in self.coefficients.iter().skip(1) {
             // early return if possible (not possible for rationals)
-            if F::one_is_gcd_unit() && self.field.is_one(&c) {
+            if F::one_is_gcd_unit() && self.ring.is_one(&c) {
                 break;
             }
 
-            c = self.field.gcd(&c, cc);
+            c = self.ring.gcd(&c, cc);
         }
         c
     }
@@ -1380,7 +1387,7 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
     /// Divide every coefficient with `other`.
     pub fn div_coeff(mut self, other: &F::Element) -> Self {
         for c in &mut self.coefficients {
-            let (quot, rem) = self.field.quot_rem(c, other);
+            let (quot, rem) = self.ring.quot_rem(c, other);
             debug_assert!(F::is_zero(&rem));
             *c = quot;
         }
@@ -1407,7 +1414,7 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
         }
 
         // check if the leading coefficients divide
-        if !F::is_zero(&self.field.rem(&self.lcoeff(), &div.lcoeff())) {
+        if !F::is_zero(&self.ring.rem(&self.lcoeff(), &div.lcoeff())) {
             return None;
         }
 
@@ -1415,33 +1422,33 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
             return None;
         }
 
-        if self.field.characteristic().is_zero() {
+        if self.ring.characteristic().is_zero() {
             // test division of constant term (evaluation at x_i = 0)
             let c = div.get_constant();
             if !F::is_zero(&c)
-                && !self.field.is_one(&c)
-                && !F::is_zero(&self.field.rem(&self.get_constant(), &c))
+                && !self.ring.is_one(&c)
+                && !F::is_zero(&self.ring.rem(&self.get_constant(), &c))
             {
                 return None;
             }
 
             // test division at x_i = 1
-            let mut num = self.field.zero();
+            let mut num = self.ring.zero();
             for c in &self.coefficients {
                 if !F::is_zero(c) {
-                    self.field.add_assign(&mut num, c);
+                    self.ring.add_assign(&mut num, c);
                 }
             }
-            let mut den = self.field.zero();
+            let mut den = self.ring.zero();
             for c in &div.coefficients {
                 if !F::is_zero(c) {
-                    self.field.add_assign(&mut den, c);
+                    self.ring.add_assign(&mut den, c);
                 }
             }
 
             if !F::is_zero(&den)
-                && !self.field.is_one(&den)
-                && !F::is_zero(&self.field.rem(&num, &den))
+                && !self.ring.is_one(&den)
+                && !F::is_zero(&self.ring.rem(&num, &den))
             {
                 return None;
             }
@@ -1481,14 +1488,12 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
         }
 
         let mut q = self.zero();
-        q.coefficients = vec![self.field.zero(); n + 1 - m];
+        q.coefficients = vec![self.ring.zero(); n + 1 - m];
 
         let mut r = self.clone();
 
         while n >= m {
-            let (qq, rr) = self
-                .field
-                .quot_rem(&r.coefficients[n], &div.coefficients[m]);
+            let (qq, rr) = self.ring.quot_rem(&r.coefficients[n], &div.coefficients[m]);
             if !F::is_zero(&rr) {
                 return (self.zero(), r);
             }
@@ -1534,7 +1539,7 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
 
         let mut res = self.zero();
         res.coefficients
-            .resize(self.coefficients.len() + 1, self.field.zero());
+            .resize(self.coefficients.len() + 1, self.ring.zero());
 
         for (p, (nc, oc)) in res
             .coefficients
@@ -1544,7 +1549,7 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
             .enumerate()
         {
             if !F::is_zero(oc) {
-                let (q, r) = self.field.quot_rem(oc, &self.field.nth(p as u64 + 1));
+                let (q, r) = self.ring.quot_rem(oc, &self.ring.nth(p as u64 + 1));
                 if !F::is_zero(&r) {
                     panic!("Could not compute integral since there is a remainder in the division of the exponent number.");
                 }
@@ -1560,8 +1565,8 @@ impl<F: Field> UnivariatePolynomial<F> {
     /// Make the polynomial monic, i.e., make the leading coefficient `1` by
     /// multiplying all monomials with `1/lcoeff`.
     pub fn make_monic(self) -> Self {
-        if self.lcoeff() != self.field.one() {
-            let ci = self.field.inv(&self.lcoeff());
+        if self.lcoeff() != self.ring.one() {
+            let ci = self.ring.inv(&self.lcoeff());
             self.mul_coeff(&ci)
         } else {
             self
@@ -1599,10 +1604,10 @@ impl<F: Field> UnivariatePolynomial<F> {
 
         let mut r0 = self.clone().make_monic();
         let mut r1 = other.clone().make_monic();
-        let mut s0 = self.constant(self.field.inv(&self.lcoeff()));
+        let mut s0 = self.constant(self.ring.inv(&self.lcoeff()));
         let mut s1 = self.zero();
         let mut t0 = self.zero();
-        let mut t1 = self.constant(self.field.inv(&other.lcoeff()));
+        let mut t1 = self.constant(self.ring.inv(&other.lcoeff()));
 
         while !r1.is_zero() {
             let (q, r) = r0.quot_rem(&r1);
@@ -1610,7 +1615,7 @@ impl<F: Field> UnivariatePolynomial<F> {
                 return (r1, s1, t1);
             }
 
-            let a = self.field.inv(&r.lcoeff());
+            let a = self.ring.inv(&r.lcoeff());
             (r1, r0) = (r.mul_coeff(&a), r1);
             (s1, s0) = ((s0 - &q * &s1).mul_coeff(&a), s1);
             (t1, t0) = ((t0 - q * &t1).mul_coeff(&a), t1);
@@ -1682,7 +1687,7 @@ impl<F: Field> UnivariatePolynomial<F> {
         // normalize the gcd
         let l = d.coefficients.last().unwrap().clone();
         for x in &mut d.coefficients {
-            self.field.div_assign(x, &l);
+            self.ring.div_assign(x, &l);
         }
 
         d
@@ -1705,15 +1710,15 @@ impl<F: Field> UnivariatePolynomial<F> {
         let mut n = self.degree();
         let m = div.degree();
 
-        let u = self.field.inv(&div.coefficients[m]);
+        let u = self.ring.inv(&div.coefficients[m]);
 
         let mut q = self.zero();
-        q.coefficients = vec![self.field.zero(); n - m + 1];
+        q.coefficients = vec![self.ring.zero(); n - m + 1];
 
         let mut r = self.clone();
 
         while n >= m {
-            let qq = self.field.mul(&r.coefficients[n], &u);
+            let qq = self.ring.mul(&r.coefficients[n], &u);
             r = r - div.mul_exp(n - m).mul_coeff(&qq);
             q.coefficients[n - m] = qq;
             n = r.degree();
@@ -1729,7 +1734,7 @@ impl<R: Ring, E: Exponent> 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() {
-            return self.field.zero();
+            return self.ring.zero();
         }
 
         let Some(pos) = self.coefficients[0]
@@ -1742,7 +1747,7 @@ impl<R: Ring, E: Exponent> UnivariatePolynomial<PolynomialRing<R, E>> {
 
         let n_vars = self.coefficients[0].get_vars().len();
         let mut res = MultivariatePolynomial::new(
-            &self.field.ring,
+            &self.ring.ring,
             self.degree().into(),
             self.coefficients[0].get_vars().clone(),
         );
diff --git a/symbolica.pyi b/symbolica.pyi
index 8228f52..4d1bc54 100644
--- a/symbolica.pyi
+++ b/symbolica.pyi
@@ -376,7 +376,7 @@ class Expression:
     def get_byte_size(self) -> int:
         """ Get the number of bytes that this expression takes up in memory."""
 
-    def pretty_str(
+    def format(
         self,
         terms_on_new_line: bool = False,
         color_top_level_sum: bool = True,
@@ -390,6 +390,7 @@ class Expression:
         square_brackets_for_function: bool = False,
         num_exp_as_superscript: bool = True,
         latex: bool = False,
+        precision: Optional[int] = None,
     ) -> str:
         """
         Convert the expression into a human-readable string, with tunable settings.
@@ -397,7 +398,7 @@ class Expression:
         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.pretty_str(number_thousands_separator='_', multiplication_operator=' '))
+        >>> print(a.format(number_thousands_separator='_', multiplication_operator=' '))
 
         Yields `z³⁴+x^(x+2)+y⁴+f(x,x²)+128_378_127_123 z^(2/3) w² x⁻¹ y⁻¹+3/5`.
         """
@@ -2219,7 +2220,7 @@ class Polynomial:
     def to_latex(self) -> str:
         """Convert the polynomial into a LaTeX string."""
 
-    def pretty_str(
+    def format(
         self,
         terms_on_new_line: bool = False,
         color_top_level_sum: bool = True,
@@ -2233,6 +2234,7 @@ class Polynomial:
         square_brackets_for_function: bool = False,
         num_exp_as_superscript: bool = True,
         latex: bool = False,
+        precision: Optional[int] = None,
     ) -> str:
         """
         Convert the polynomial into a human-readable string, with tunable settings.
@@ -2240,7 +2242,7 @@ class Polynomial:
         Examples
         --------
         >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-        >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+        >>> print(p.format(symmetric_representation_for_finite_field=True))
 
         Yields `z³⁴+x^(x+2)+y⁴+f(x,x²)+128_378_127_123 z^(2/3) w² x⁻¹ y⁻¹+3/5`.
         """
@@ -2488,7 +2490,7 @@ class IntegerPolynomial:
     def to_latex(self) -> str:
         """Convert the polynomial into a LaTeX string."""
 
-    def pretty_str(
+    def format(
         self,
         terms_on_new_line: bool = False,
         color_top_level_sum: bool = True,
@@ -2502,6 +2504,7 @@ class IntegerPolynomial:
         square_brackets_for_function: bool = False,
         num_exp_as_superscript: bool = True,
         latex: bool = False,
+        precision: Optional[int] = None,
     ) -> str:
         """
         Convert the polynomial into a human-readable string, with tunable settings.
@@ -2509,7 +2512,7 @@ class IntegerPolynomial:
         Examples
         --------
         >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-        >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+        >>> print(p.format(symmetric_representation_for_finite_field=True))
 
         Yields `z³⁴+x^(x+2)+y⁴+f(x,x²)+128_378_127_123 z^(2/3) w² x⁻¹ y⁻¹+3/5`.
         """
@@ -2649,7 +2652,7 @@ class NumberFieldPolynomial:
     def to_latex(self) -> str:
         """Convert the polynomial into a LaTeX string."""
 
-    def pretty_str(
+    def format(
         self,
         terms_on_new_line: bool = False,
         color_top_level_sum: bool = True,
@@ -2663,6 +2666,7 @@ class NumberFieldPolynomial:
         square_brackets_for_function: bool = False,
         num_exp_as_superscript: bool = True,
         latex: bool = False,
+        precision: Optional[int] = None,
     ) -> str:
         """
         Convert the polynomial into a human-readable string, with tunable settings.
@@ -2670,7 +2674,7 @@ class NumberFieldPolynomial:
         Examples
         --------
         >>> p = FiniteFieldNumberFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-        >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+        >>> print(p.format(symmetric_representation_for_finite_field=True))
 
         Yields `z³⁴+x^(x+2)+y⁴+f(x,x²)+128_378_127_123 z^(2/3) w² x⁻¹ y⁻¹+3/5`.
         """
@@ -2859,7 +2863,7 @@ class FiniteFieldPolynomial:
     def to_latex(self) -> str:
         """Convert the polynomial into a LaTeX string."""
 
-    def pretty_str(
+    def format(
         self,
         terms_on_new_line: bool = False,
         color_top_level_sum: bool = True,
@@ -2873,6 +2877,7 @@ class FiniteFieldPolynomial:
         square_brackets_for_function: bool = False,
         num_exp_as_superscript: bool = True,
         latex: bool = False,
+        precision: Optional[int] = None,
     ) -> str:
         """
         Convert the polynomial into a human-readable string, with tunable settings.
@@ -2880,7 +2885,7 @@ class FiniteFieldPolynomial:
         Examples
         --------
         >>> p = FiniteFieldPolynomial.parse("3*x^2+2*x+7*x^3", ['x'], 11)
-        >>> print(p.pretty_str(symmetric_representation_for_finite_field=True))
+        >>> print(p.format(symmetric_representation_for_finite_field=True))
 
         Yields `z³⁴+x^(x+2)+y⁴+f(x,x²)+128_378_127_123 z^(2/3) w² x⁻¹ y⁻¹+3/5`.
         """