diff --git a/crates/math/src/field/fields/binary/README.md b/crates/math/src/field/fields/binary/README.md
index e6d5bcef4..768634d41 100644
--- a/crates/math/src/field/fields/binary/README.md
+++ b/crates/math/src/field/fields/binary/README.md
@@ -33,12 +33,14 @@ $$\begin{aligned}
 
 **Level 2:** At level 2, we define the field extension $GF(2^{2^2})$. In this case the elements are polynomials with binary coefficients and two variables, $x_0$ and $x_1$. The first one keeps satisfying ${x_0}^2 = x_0 + 1$ and in addition the second one satisfies ${x_1}^2 = x_1 \cdot x_0 + 1$. This means that the polynomials are lineal in each variable. Therefore, this field extension has $2^{2^2}$ elements:
 
-$\begin{array}{llll}
-0000 = 0  & 0100 = x_1  & 1000 = x_1x_0 & 1100 = x_1x_0 + x_1 \\
-0001 = 1  & 0101 = x_1 + 1  & 1001 = x_1x_0 + 1  & 1101 = x_1x_0 + x_1 + 1 \\
-0010 = x_0  & 0110 = x_1 + x_0  & 1010 = x_1x_0 + x_0  & 1110 = x_1x_0 + x_1 + x_0 \\
-0011 = x_0 + 1  & 0111 = x_1 + x_0 + 1  & 1011 = x_1x_0 + x_0 + 1  & 1111 = x_1x_0 + x_1 + x_0 + 1
-\end{array}$
+$$
+\begin{array}{llll}
+0000 = 0 & 0100 = x_1 & 1000 = x_1x_0 & 1100 = x_1x_0 + x_1 \\
+0001 = 1 & 0101 = x_1 + 1 & 1001 = x_1x_0 + 1 & 1101 = x_1x_0 + x_1 + 1 \\
+0010 = x_0 & 0110 = x_1 + x_0 & 1010 = x_1x_0 + x_0 & 1110 = x_1x_0 + x_1 + x_0 \\
+0011 = x_0 + 1 & 0111 = x_1 + x_0 + 1 & 1011 = x_1x_0 + x_0 + 1 & 1111 = x_1x_0 + x_1 + x_0 + 1
+\end{array}
+$$ 
 
 **Level 3:** At level 3, we define $GF(2^{2^3})$ in the same way. This time the polynomials have three variables $x_0$, $x_1$ and $x_2$. The first two variables satisfy the equations mentioned before and in addition the last one satisfies ${x_2}^2 = x_2 \cdot x_1 + 1$. This field extension has $2^{2^3}$ elements.
 
diff --git a/crates/math/src/field/fields/binary/field.rs b/crates/math/src/field/fields/binary/field.rs
index 6f9562e93..4637653ce 100644
--- a/crates/math/src/field/fields/binary/field.rs
+++ b/crates/math/src/field/fields/binary/field.rs
@@ -96,6 +96,57 @@ impl TowerFieldElement {
         format!("{:0width$b}", self.value, width = self.num_bits())
     }
 
+    /// Returns polynomial representation with variables
+    #[cfg(feature = "std")]
+    pub fn to_polynomial_string(&self) -> String {
+        if self.is_zero() {
+            return "0".to_string();
+        }
+
+        debug_assert!(
+            self.value < (1 << self.num_bits()),
+            "Invalid value for level"
+        );
+
+        let mut terms = Vec::new();
+        let num_possible_monomials = self.num_bits();
+        let num_vars = self.num_level();
+
+        // Collect terms in reverse order to get higher degree terms first
+        for i in (0..num_possible_monomials).rev() {
+            if (self.value >> i) & 1 == 1 {
+                let term = match i {
+                    0 => "1".to_string(),
+                    1 => "x_0".to_string(),
+                    2 => "x_1".to_string(),
+                    3 => "x_1 x_0".to_string(),
+                    4 => "x_2".to_string(),
+                    5 => "x_2 x_0".to_string(),
+                    6 => "x_2 x_1".to_string(),
+                    7 => "x_2 x_1 x_0".to_string(),
+                    _ => {
+                        let mut monomial_str = String::new();
+                        let mut first = true;
+                        // Iterate in reverse to get higher indices first
+                        for j in (0..num_vars).rev() {
+                            if ((i >> j) & 1) == 1 {
+                                if !first {
+                                    monomial_str.push('x');
+                                }
+                                monomial_str.push_str(&format!("_{}", j));
+                                first = false;
+                            }
+                        }
+                        monomial_str
+                    }
+                };
+                terms.push(term);
+            }
+        }
+
+        terms.join(" + ")
+    }
+
     /// Splits element into high and low parts.
     /// For example, if a = xy + y + x, then a = (x + 1)y + x and
     /// therefore, a_hi = x + 1 and a_lo = x.
@@ -469,6 +520,29 @@ mod tests {
         assert_eq!(a * b, c);
     }
 
+    #[test]
+    fn mul_in_level_3() {
+        let a = TowerFieldElement::new(0b01000010, 3);
+        let b = TowerFieldElement::new(0b00100101, 3);
+        let result = a * b;
+
+        println!(
+            "a = {} (binary) = {} (polynomial)",
+            a.to_binary_string(),
+            a.to_polynomial_string()
+        );
+        println!(
+            "b = {} (binary) = {} (polynomial)",
+            b.to_binary_string(),
+            b.to_polynomial_string()
+        );
+        println!(
+            "result = {} (binary) = {} (polynomial)",
+            result.to_binary_string(),
+            result.to_polynomial_string()
+        );
+    }
+
     #[test]
     fn test_correct_level_mul() {
         let a = TowerFieldElement::new(0b1111, 5);