diff --git a/CHANGELOG.md b/CHANGELOG.md
index 6d275092..70aab304 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Fixed
 
 - Fix inconsistency in gate ordering of arithmetic verifier key [#797]
+- Fix leading coefficients might be zero [#796]
 
 ## [0.18.0] - 2023-12-13
 
@@ -537,6 +538,7 @@ is necessary since `rkyv/validation` was required as a bound.
 
 <!-- ISSUES -->
 [#797]: https://github.com/dusk-network/plonk/issues/797
+[#796]: https://github.com/dusk-network/plonk/issues/796
 [#784]: https://github.com/dusk-network/plonk/issues/784
 [#773]: https://github.com/dusk-network/plonk/issues/773
 [#774]: https://github.com/dusk-network/plonk/issues/774
diff --git a/src/fft/polynomial.rs b/src/fft/polynomial.rs
index f08bc831..4276482c 100644
--- a/src/fft/polynomial.rs
+++ b/src/fft/polynomial.rs
@@ -147,7 +147,10 @@ impl Polynomial {
             .map(BlsScalar::from_slice)
             .collect::<Result<Vec<BlsScalar>, dusk_bytes::Error>>()?;
 
-        Ok(Polynomial { coeffs })
+        let mut p = Polynomial { coeffs };
+        p.truncate_leading_zeros();
+
+        Ok(p)
     }
 
     /// Returns an iterator over the polynomial coefficients.
@@ -213,8 +216,9 @@ impl<'a> AddAssign<&'a Polynomial> for Polynomial {
             for (a, b) in self.coeffs.iter_mut().zip(&other.coeffs) {
                 *a += b
             }
-            self.truncate_leading_zeros();
         }
+        // If the leading coefficient ends up being zero, pop it off.
+        self.truncate_leading_zeros();
     }
 }
 
@@ -235,8 +239,9 @@ impl<'a> AddAssign<(BlsScalar, &'a Polynomial)> for Polynomial {
             for (a, b) in self.coeffs.iter_mut().zip(&other.coeffs) {
                 *a += &(f * b);
             }
-            self.truncate_leading_zeros();
         }
+        // If the leading coefficient ends up being zero, pop it off.
+        self.truncate_leading_zeros();
     }
 }
 
@@ -303,9 +308,9 @@ impl<'a> SubAssign<&'a Polynomial> for Polynomial {
             for (a, b) in self.coeffs.iter_mut().zip(&other.coeffs) {
                 *a -= b
             }
-            // If the leading coefficient ends up being zero, pop it off.
-            self.truncate_leading_zeros();
         }
+        // If the leading coefficient ends up being zero, pop it off.
+        self.truncate_leading_zeros();
     }
 }
 
@@ -313,7 +318,14 @@ impl Polynomial {
     #[allow(dead_code)]
     #[inline]
     fn leading_coefficient(&self) -> Option<&BlsScalar> {
-        self.last()
+        // skip leading coefficients that are zero
+        let len = self.len();
+        for i in 0..len {
+            if self[len - 1 - i] != BlsScalar::zero() {
+                return Some(&self[len - i]);
+            }
+        }
+        None
     }
 
     #[allow(dead_code)]
@@ -458,6 +470,7 @@ mod test {
         ]);
         assert_eq!(quotient, expected_quotient);
     }
+
     #[test]
     fn test_ruffini_zero() {
         // Tests the two situations where zero can be added to Ruffini:
@@ -465,15 +478,15 @@ mod test {
         // (2) Zero as the constant term for the polynomial you are dividing by
         // In the first case, we should get zero as the quotient
         // In the second case, this is the same as dividing by X
+
         // (1)
-        //
         // Zero polynomial
         let zero = Polynomial::zero();
         // Quotient is invariant under any argument we pass
         let quotient = zero.ruffini(-BlsScalar::from(2));
         assert_eq!(quotient, Polynomial::zero());
+
         // (2)
-        //
         // X^2 + X
         let p = Polynomial::from_coefficients_vec(vec![
             BlsScalar::zero(),
@@ -489,4 +502,46 @@ mod test {
         ]);
         assert_eq!(quotient, expected_quotient);
     }
+
+    #[test]
+    fn test_add_assign() {
+        let mut p1 = Polynomial::from_coefficients_vec(vec![
+            BlsScalar::from(21),
+            BlsScalar::zero(),
+            BlsScalar::from(1),
+        ]);
+        let p2 = Polynomial::from_coefficients_vec(vec![
+            BlsScalar::from(21),
+            BlsScalar::zero(),
+            -BlsScalar::from(1),
+        ]);
+
+        p1 += &p2;
+
+        assert_eq!(
+            p1,
+            Polynomial::from_coefficients_vec(vec![BlsScalar::from(42)])
+        );
+    }
+
+    #[test]
+    fn test_sub_assign() {
+        let mut p1 = Polynomial::from_coefficients_vec(vec![
+            BlsScalar::from(21),
+            BlsScalar::zero(),
+            BlsScalar::from(1),
+        ]);
+        let p2 = Polynomial::from_coefficients_vec(vec![
+            -BlsScalar::from(21),
+            BlsScalar::zero(),
+            BlsScalar::from(1),
+        ]);
+
+        p1 -= &p2;
+
+        assert_eq!(
+            p1,
+            Polynomial::from_coefficients_vec(vec![BlsScalar::from(42)])
+        );
+    }
 }