diff --git a/src/field/lagrange.rs b/src/field/interpolation.rs
similarity index 81%
rename from src/field/lagrange.rs
rename to src/field/interpolation.rs
index 8911feb0..968eeca2 100644
--- a/src/field/lagrange.rs
+++ b/src/field/interpolation.rs
@@ -65,11 +65,23 @@ pub fn barycentric_weights<F: Field>(points: &[(F, F)]) -> Vec<F> {
     )
 }
 
+/// Interpolate the linear polynomial passing through `points` on `x`.
+pub fn interpolate2<F: Field>(points: [(F, F); 2], x: F) -> F {
+    // a0 -> a1
+    // b0 -> b1
+    // x  -> a1 + (x-a0)*(b1-a1)/(b0-a0)
+    let (a0, a1) = points[0];
+    let (b0, b1) = points[1];
+    assert_ne!(a0, b0);
+    a1 + (x - a0) * (b1 - a1) / (b0 - a0)
+}
+
 #[cfg(test)]
 mod tests {
+    use super::*;
     use crate::field::crandall_field::CrandallField;
+    use crate::field::extension_field::quartic::QuarticCrandallField;
     use crate::field::field::Field;
-    use crate::field::lagrange::interpolant;
     use crate::polynomial::polynomial::PolynomialCoeffs;
 
     #[test]
@@ -120,4 +132,18 @@ mod tests {
     fn eval_naive<F: Field>(coeffs: &PolynomialCoeffs<F>, domain: &[F]) -> Vec<(F, F)> {
         domain.iter().map(|&x| (x, coeffs.eval(x))).collect()
     }
+
+    #[test]
+    fn test_interpolate2() {
+        type F = QuarticCrandallField;
+        let points = [(F::rand(), F::rand()), (F::rand(), F::rand())];
+        let x = F::rand();
+
+        let ev0 = interpolant(&points).eval(x);
+        let ev1 = interpolate(&points, x, &barycentric_weights(&points));
+        let ev2 = interpolate2(points, x);
+
+        assert_eq!(ev0, ev1);
+        assert_eq!(ev0, ev2);
+    }
 }
diff --git a/src/field/mod.rs b/src/field/mod.rs
index 179fb10d..15efe280 100644
--- a/src/field/mod.rs
+++ b/src/field/mod.rs
@@ -3,7 +3,7 @@ pub mod crandall_field;
 pub mod extension_field;
 pub mod fft;
 pub mod field;
-pub(crate) mod lagrange;
+pub(crate) mod interpolation;
 
 #[cfg(test)]
 mod field_testing;
diff --git a/src/gadgets/interpolation.rs b/src/gadgets/interpolation.rs
index a94eb582..40b91c1e 100644
--- a/src/gadgets/interpolation.rs
+++ b/src/gadgets/interpolation.rs
@@ -1,9 +1,10 @@
+use std::marker::PhantomData;
+
 use crate::circuit_builder::CircuitBuilder;
 use crate::field::extension_field::target::ExtensionTarget;
 use crate::field::extension_field::Extendable;
 use crate::gates::interpolation::InterpolationGate;
 use crate::target::Target;
-use std::marker::PhantomData;
 
 impl<F: Extendable<D>, const D: usize> CircuitBuilder<F, D> {
     /// Interpolate two points. No need for an `InterpolationGate` since the coefficients
@@ -56,15 +57,16 @@ impl<F: Extendable<D>, const D: usize> CircuitBuilder<F, D> {
 
 #[cfg(test)]
 mod tests {
+    use std::convert::TryInto;
+
     use super::*;
     use crate::circuit_data::CircuitConfig;
     use crate::field::crandall_field::CrandallField;
     use crate::field::extension_field::quartic::QuarticCrandallField;
     use crate::field::extension_field::FieldExtension;
     use crate::field::field::Field;
-    use crate::field::lagrange::{interpolant, interpolate};
+    use crate::field::interpolation::{interpolant, interpolate};
     use crate::witness::PartialWitness;
-    use std::convert::TryInto;
 
     #[test]
     fn test_interpolate() {
diff --git a/src/gates/interpolation.rs b/src/gates/interpolation.rs
index ac2ca49f..ccf8d57d 100644
--- a/src/gates/interpolation.rs
+++ b/src/gates/interpolation.rs
@@ -6,7 +6,7 @@ use crate::circuit_builder::CircuitBuilder;
 use crate::field::extension_field::algebra::PolynomialCoeffsAlgebra;
 use crate::field::extension_field::target::ExtensionTarget;
 use crate::field::extension_field::{Extendable, FieldExtension};
-use crate::field::lagrange::interpolant;
+use crate::field::interpolation::interpolant;
 use crate::gadgets::polynomial::PolynomialCoeffsExtAlgebraTarget;
 use crate::gates::gate::{Gate, GateRef};
 use crate::generator::{SimpleGenerator, WitnessGenerator};
diff --git a/src/polynomial/division.rs b/src/polynomial/division.rs
index b74fd00f..50e1f8a6 100644
--- a/src/polynomial/division.rs
+++ b/src/polynomial/division.rs
@@ -26,7 +26,7 @@ impl<F: Field> PolynomialCoeffs<F> {
                 .to_vec()
                 .into();
             let mut q = rev_q.rev();
-            let mut qb = &q * b;
+            let qb = &q * b;
             let mut r = self - &qb;
             q.trim();
             r.trim();
@@ -59,8 +59,7 @@ impl<F: Field> PolynomialCoeffs<F> {
                 quotient.coeffs[cur_q_degree] = cur_q_coeff;
 
                 for (i, &div_coeff) in b.coeffs.iter().enumerate() {
-                    remainder.coeffs[cur_q_degree + i] =
-                        remainder.coeffs[cur_q_degree + i] - (cur_q_coeff * div_coeff);
+                    remainder.coeffs[cur_q_degree + i] -= cur_q_coeff * div_coeff;
                 }
                 remainder.trim();
             }
@@ -97,7 +96,7 @@ impl<F: Field> PolynomialCoeffs<F> {
         let denominators = (0..a_eval.len())
             .map(|i| {
                 if i != 0 {
-                    root_pow = root_pow * root_n;
+                    root_pow *= root_n;
                 }
                 denominator_g * root_pow - F::ONE
             })
@@ -125,8 +124,25 @@ impl<F: Field> PolynomialCoeffs<F> {
         p
     }
 
+    /// Let `self=p(X)`, this returns `(p(X)-p(z))/(X-z)` and `p(z)`.
+    /// See https://en.wikipedia.org/wiki/Horner%27s_method
+    pub(crate) fn divide_by_linear(&self, z: F) -> (PolynomialCoeffs<F>, F) {
+        let mut bs = self
+            .coeffs
+            .iter()
+            .rev()
+            .scan(F::ZERO, |acc, &c| {
+                *acc = *acc * z + c;
+                Some(*acc)
+            })
+            .collect::<Vec<_>>();
+        let ev = bs.pop().unwrap_or(F::ZERO);
+        bs.reverse();
+        (Self { coeffs: bs }, ev)
+    }
+
     /// Computes the inverse of `self` modulo `x^n`.
-    pub(crate) fn inv_mod_xn(&self, n: usize) -> Self {
+    pub fn inv_mod_xn(&self, n: usize) -> Self {
         assert!(self.coeffs[0].is_nonzero(), "Inverse doesn't exist.");
 
         let h = if self.len() < n {
@@ -166,7 +182,10 @@ impl<F: Field> PolynomialCoeffs<F> {
 
 #[cfg(test)]
 mod tests {
+    use std::time::Instant;
+
     use crate::field::crandall_field::CrandallField;
+    use crate::field::extension_field::quartic::QuarticCrandallField;
     use crate::field::field::Field;
     use crate::polynomial::polynomial::PolynomialCoeffs;
 
@@ -199,4 +218,49 @@ mod tests {
         let computed_q = a.divide_by_z_h(4);
         assert_eq!(computed_q, q);
     }
+
+    #[test]
+    #[ignore]
+    fn test_division_by_linear() {
+        type F = QuarticCrandallField;
+        let n = 1_000_000;
+        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 = QuarticCrandallField;
+        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());
+    }
 }
diff --git a/src/polynomial/polynomial.rs b/src/polynomial/polynomial.rs
index 9f605051..aefcc0c6 100644
--- a/src/polynomial/polynomial.rs
+++ b/src/polynomial/polynomial.rs
@@ -1,6 +1,6 @@
 use std::cmp::max;
 use std::iter::Sum;
-use std::ops::{Add, Mul, Sub};
+use std::ops::{Add, AddAssign, Mul, Sub, SubAssign};
 
 use crate::field::extension_field::Extendable;
 use crate::field::fft::{fft, ifft};
@@ -243,6 +243,26 @@ impl<F: Field> Sub for &PolynomialCoeffs<F> {
     }
 }
 
+impl<F: Field> AddAssign for PolynomialCoeffs<F> {
+    fn add_assign(&mut self, rhs: Self) {
+        let len = max(self.len(), rhs.len());
+        self.coeffs.resize(len, F::ZERO);
+        for (l, r) in self.coeffs.iter_mut().zip(rhs.coeffs) {
+            *l += r;
+        }
+    }
+}
+
+impl<F: Field> SubAssign for PolynomialCoeffs<F> {
+    fn sub_assign(&mut self, rhs: Self) {
+        let len = max(self.len(), rhs.len());
+        self.coeffs.resize(len, F::ZERO);
+        for (l, r) in self.coeffs.iter_mut().zip(rhs.coeffs) {
+            *l -= r;
+        }
+    }
+}
+
 impl<F: Field> Mul<F> for &PolynomialCoeffs<F> {
     type Output = PolynomialCoeffs<F>;