diff --git a/field/src/polynomial/division.rs b/field/src/polynomial/division.rs
index 4f3cafae..d761ab50 100644
--- a/field/src/polynomial/division.rs
+++ b/field/src/polynomial/division.rs
@@ -67,9 +67,9 @@ impl<F: Field> PolynomialCoeffs<F> {
         }
     }
 
-    /// Let `self=p(X)`, this returns `(p(X)-p(z))/(X-z)` and `p(z)`.
+    /// Let `self=p(X)`, this returns `(p(X)-p(z))/(X-z)`.
     /// See https://en.wikipedia.org/wiki/Horner%27s_method
-    pub fn divide_by_linear(&self, z: F) -> (PolynomialCoeffs<F>, F) {
+    pub fn divide_by_linear(&self, z: F) -> PolynomialCoeffs<F> {
         let mut bs = self
             .coeffs
             .iter()
@@ -79,9 +79,9 @@ impl<F: Field> PolynomialCoeffs<F> {
                 Some(*acc)
             })
             .collect::<Vec<_>>();
-        let ev = bs.pop().unwrap_or(F::ZERO);
+        bs.pop();
         bs.reverse();
-        (Self { coeffs: bs }, ev)
+        Self { coeffs: bs }
     }
 
     /// Computes the inverse of `self` modulo `x^n`.
@@ -125,7 +125,7 @@ impl<F: Field> PolynomialCoeffs<F> {
 
 #[cfg(test)]
 mod tests {
-    use std::time::Instant;
+    use rand::{thread_rng, Rng};
 
     use crate::extension_field::quartic::QuarticExtension;
     use crate::field_types::Field;
@@ -133,47 +133,17 @@ mod tests {
     use crate::polynomial::PolynomialCoeffs;
 
     #[test]
-    #[ignore]
     fn test_division_by_linear() {
         type F = QuarticExtension<GoldilocksField>;
-        let n = 1_000_000;
+        let n = thread_rng().gen_range(1..1000);
         let poly = PolynomialCoeffs::new(F::rand_vec(n));
         let z = F::rand();
         let ev = poly.eval(z);
 
-        let timer = Instant::now();
-        let (_quotient, ev2) = poly.div_rem(&PolynomialCoeffs::new(vec![-z, F::ONE]));
-        println!("{:.3}s for usual", timer.elapsed().as_secs_f32());
-        assert_eq!(ev2.trimmed().coeffs, vec![ev]);
-
-        let timer = Instant::now();
-        let (quotient, ev3) = poly.div_rem_long_division(&PolynomialCoeffs::new(vec![-z, F::ONE]));
-        println!("{:.3}s for long division", timer.elapsed().as_secs_f32());
-        assert_eq!(ev3.trimmed().coeffs, vec![ev]);
-
-        let timer = Instant::now();
-        let horn = poly.divide_by_linear(z);
-        println!("{:.3}s for Horner", timer.elapsed().as_secs_f32());
-        assert_eq!((quotient, ev), horn);
-    }
-
-    #[test]
-    #[ignore]
-    fn test_division_by_quadratic() {
-        type F = QuarticExtension<GoldilocksField>;
-        let n = 1_000_000;
-        let poly = PolynomialCoeffs::new(F::rand_vec(n));
-        let quad = PolynomialCoeffs::new(F::rand_vec(2));
-
-        let timer = Instant::now();
-        let (quotient0, rem0) = poly.div_rem(&quad);
-        println!("{:.3}s for usual", timer.elapsed().as_secs_f32());
-
-        let timer = Instant::now();
-        let (quotient1, rem1) = poly.div_rem_long_division(&quad);
-        println!("{:.3}s for long division", timer.elapsed().as_secs_f32());
-
-        assert_eq!(quotient0.trimmed(), quotient1.trimmed());
-        assert_eq!(rem0.trimmed(), rem1.trimmed());
+        let quotient = poly.divide_by_linear(z);
+        assert_eq!(
+            poly,
+            &(&quotient * &vec![-z, F::ONE].into()) + &vec![ev].into() // `quotient * (X-z) + ev`
+        );
     }
 }
diff --git a/plonky2/src/fri/oracle.rs b/plonky2/src/fri/oracle.rs
index f1dc84eb..c705e125 100644
--- a/plonky2/src/fri/oracle.rs
+++ b/plonky2/src/fri/oracle.rs
@@ -150,11 +150,10 @@ impl<F: RichField + Extendable<D>, C: GenericConfig<D, F = F>, const D: usize>
                 &format!("reduce batch of {} polynomials", polynomials.len()),
                 alpha.reduce_polys_base(polys_coeff)
             );
-            let quotient = Self::compute_quotient([*point], composition_poly);
+            let quotient = composition_poly.divide_by_linear(*point);
             alpha.shift_poly(&mut final_poly);
             final_poly += quotient;
         }
-        final_poly.trim();
         // Multiply the final polynomial by `X`, so that `final_poly` has the maximum degree for
         // which the LDT will pass. See github.com/mir-protocol/plonky2/pull/436 for details.
         final_poly.coeffs.insert(0, F::Extension::ZERO);
@@ -180,28 +179,4 @@ impl<F: RichField + Extendable<D>, C: GenericConfig<D, F = F>, const D: usize>
 
         fri_proof
     }
-
-    /// Given `points=(x_i)`, `evals=(y_i)` and `poly=P` with `P(x_i)=y_i`, computes the polynomial
-    /// `Q=(P-I)/Z` where `I` interpolates `(x_i, y_i)` and `Z` is the vanishing polynomial on `(x_i)`.
-    fn compute_quotient<const N: usize>(
-        points: [F::Extension; N],
-        poly: PolynomialCoeffs<F::Extension>,
-    ) -> PolynomialCoeffs<F::Extension> {
-        let quotient = if N == 1 {
-            poly.divide_by_linear(points[0]).0
-        } else if N == 2 {
-            // The denominator is `(X - p0)(X - p1) = p0 p1 - (p0 + p1) X + X^2`.
-            let denominator = vec![
-                points[0] * points[1],
-                -points[0] - points[1],
-                F::Extension::ONE,
-            ]
-            .into();
-            poly.div_rem_long_division(&denominator).0 // Could also use `divide_by_linear` twice.
-        } else {
-            unreachable!("This shouldn't happen. Plonk should open polynomials at 1 or 2 points.")
-        };
-
-        quotient.padded(quotient.degree_plus_one().next_power_of_two())
-    }
 }