diff --git a/derive/src/field/mod.rs b/derive/src/field/mod.rs
index 2f581037..f6e4c836 100644
--- a/derive/src/field/mod.rs
+++ b/derive/src/field/mod.rs
@@ -282,14 +282,6 @@ pub(crate) fn impl_field(input: TokenStream) -> TokenStream {
             }
         }
 
-        impl ::core::ops::Neg for #field {
-            type Output = #field;
-
-            #[inline]
-            fn neg(self) -> #field {
-                -&self
-            }
-        }
 
         impl #field {
             pub const SIZE: usize = #num_limbs * 8;
@@ -393,8 +385,17 @@ pub(crate) fn impl_field(input: TokenStream) -> TokenStream {
                 self.square()
             }
 
+            // Returns the multiplicative inverse of the element. If it is zero, the method fails.
+            #[inline(always)]
             fn invert(&self) -> CtOption<Self> {
-                self.invert()
+                const BYINVERTOR: crate::ff_ext::inverse::BYInverter<#by_inverter_constant> =
+                crate::ff_ext::inverse::BYInverter::<#by_inverter_constant>::new(&#modulus_limbs_ident, &#r2);
+
+                if let Some(inverse) = BYINVERTOR.invert::<{ Self::NUM_LIMBS }>(&self.0) {
+                    subtle::CtOption::new(Self(inverse), subtle::Choice::from(1))
+                } else {
+                    subtle::CtOption::new(Self::zero(), subtle::Choice::from(0))
+                }
             }
 
             #sqrt_impl
diff --git a/src/bn256/curve.rs b/src/bn256/curve.rs
index 0160956a..98f38ee4 100644
--- a/src/bn256/curve.rs
+++ b/src/bn256/curve.rs
@@ -12,9 +12,8 @@ use crate::ff::{Field, PrimeField};
 use crate::group::Curve;
 use crate::group::{cofactor::CofactorGroup, prime::PrimeCurveAffine, Group, GroupEncoding};
 use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    new_curve_impl,
+    impl_binops_additive, impl_binops_additive_specify_output, impl_binops_multiplicative,
+    impl_binops_multiplicative_mixed, new_curve_impl,
 };
 use crate::{Coordinates, CurveAffine, CurveExt};
 use core::cmp;
diff --git a/src/bn256/engine.rs b/src/bn256/engine.rs
index b55b00f3..630a36f2 100644
--- a/src/bn256/engine.rs
+++ b/src/bn256/engine.rs
@@ -1,4 +1,3 @@
-#![allow(clippy::suspicious_arithmetic_impl)]
 use crate::bn256::curve::*;
 use crate::bn256::fq::*;
 use crate::bn256::fq12::*;
@@ -6,461 +5,27 @@ use crate::bn256::fq2::*;
 use crate::bn256::fq6::FROBENIUS_COEFF_FQ6_C1;
 use crate::bn256::fr::*;
 use crate::ff::PrimeField;
+use crate::ff_ext::quadratic::QuadSparseMul;
+use crate::ff_ext::ExtField;
 use crate::group::cofactor::CofactorCurveAffine;
 use crate::group::Group;
 use core::borrow::Borrow;
 use core::iter::Sum;
-use core::ops::{Add, Mul, MulAssign, Neg, Sub};
+use core::ops::{Add, Mul, Neg, Sub};
 use pairing::{Engine, MillerLoopResult, MultiMillerLoop, PairingCurveAffine};
 use rand_core::RngCore;
+use std::ops::MulAssign;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq};
 
-pub const BN_X: u64 = 4965661367192848881;
-
-// 6U+2 for in NAF form
-pub const SIX_U_PLUS_2_NAF: [i8; 65] = [
-    0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0,
-    1, 1, 1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 0, 1, 1, 0, -1, 0,
-    0, 1, 0, 1, 1,
-];
-
-pub const XI_TO_Q_MINUS_1_OVER_2: Fq2 = Fq2 {
-    c0: Fq([
-        0xe4bbdd0c2936b629,
-        0xbb30f162e133bacb,
-        0x31a9d1b6f9645366,
-        0x253570bea500f8dd,
-    ]),
-    c1: Fq([
-        0xa1d77ce45ffe77c7,
-        0x07affd117826d1db,
-        0x6d16bd27bb7edc6b,
-        0x2c87200285defecc,
-    ]),
-};
-
-impl PairingCurveAffine for G1Affine {
-    type Pair = G2Affine;
-    type PairingResult = Gt;
-
-    fn pairing_with(&self, other: &Self::Pair) -> Self::PairingResult {
-        pairing(self, other)
-    }
-}
-
-impl PairingCurveAffine for G2Affine {
-    type Pair = G1Affine;
-    type PairingResult = Gt;
-
-    fn pairing_with(&self, other: &Self::Pair) -> Self::PairingResult {
-        pairing(other, self)
-    }
-}
-
-#[derive(Copy, Clone, Debug, Default)]
-pub struct Gt(pub(crate) Fq12);
-
-impl std::fmt::Display for Gt {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(f, "{self:?}")
-    }
-}
-
-impl ConstantTimeEq for Gt {
-    fn ct_eq(&self, other: &Self) -> Choice {
-        self.0.ct_eq(&other.0)
-    }
-}
-
-impl ConditionallySelectable for Gt {
-    fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self {
-        Gt(Fq12::conditional_select(&a.0, &b.0, choice))
-    }
-}
-
-impl Eq for Gt {}
-impl PartialEq for Gt {
-    #[inline]
-    fn eq(&self, other: &Self) -> bool {
-        bool::from(self.ct_eq(other))
-    }
-}
-
-impl Gt {
-    /// Returns the group identity, which is $1$.
-    pub fn identity() -> Gt {
-        Gt(Fq12::one())
-    }
-
-    /// Doubles this group element.
-    pub fn double(&self) -> Gt {
-        Gt(self.0.square())
-    }
-}
-
-impl<'a> Neg for &'a Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn neg(self) -> Gt {
-        // The element is unitary, so we just conjugate.
-        let mut u = self.0;
-        u.conjugate();
-        Gt(u)
-    }
-}
-
-impl Neg for Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn neg(self) -> Gt {
-        -&self
-    }
-}
-
-impl<'a, 'b> Add<&'b Gt> for &'a Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn add(self, rhs: &'b Gt) -> Gt {
-        Gt(self.0 * rhs.0)
-    }
-}
-
-impl<'a, 'b> Sub<&'b Gt> for &'a Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn sub(self, rhs: &'b Gt) -> Gt {
-        self + (-rhs)
-    }
-}
-
-impl<'a, 'b> Mul<&'b Fr> for &'a Gt {
-    type Output = Gt;
-
-    fn mul(self, other: &'b Fr) -> Self::Output {
-        let mut acc = Gt::identity();
-
-        for bit in other
-            .to_repr()
-            .as_ref()
-            .iter()
-            .rev()
-            .flat_map(|byte| (0..8).rev().map(move |i| Choice::from((byte >> i) & 1u8)))
-            .skip(1)
-        {
-            acc = acc.double();
-            acc = Gt::conditional_select(&acc, &(acc + self), bit);
-        }
-
-        acc
-    }
-}
-
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-};
-impl_binops_additive!(Gt, Gt);
-impl_binops_multiplicative!(Gt, Fr);
-
-impl<T> Sum<T> for Gt
-where
-    T: Borrow<Gt>,
-{
-    fn sum<I>(iter: I) -> Self
-    where
-        I: Iterator<Item = T>,
-    {
-        iter.fold(Self::identity(), |acc, item| acc + item.borrow())
-    }
-}
-
-impl Group for Gt {
-    type Scalar = Fr;
-
-    fn random(_: impl RngCore) -> Self {
-        unimplemented!();
-    }
-
-    fn identity() -> Self {
-        Self::identity()
-    }
-
-    fn generator() -> Self {
-        unimplemented!();
-    }
-
-    fn is_identity(&self) -> Choice {
-        self.ct_eq(&Self::identity())
-    }
-
-    #[must_use]
-    fn double(&self) -> Self {
-        self.double()
-    }
-}
-
-#[derive(Clone, Debug)]
-pub struct G2Prepared {
-    pub(crate) coeffs: Vec<(Fq2, Fq2, Fq2)>,
-    pub(crate) infinity: bool,
-}
-
-impl G2Prepared {
-    pub fn is_zero(&self) -> bool {
-        self.infinity
-    }
-
-    pub fn from_affine(q: G2Affine) -> Self {
-        if bool::from(q.is_identity()) {
-            return G2Prepared {
-                coeffs: vec![],
-                infinity: true,
-            };
-        }
-
-        fn doubling_step(r: &mut G2) -> (Fq2, Fq2, Fq2) {
-            // Adaptation of Algorithm 26, https://eprint.iacr.org/2010/354.pdf
-            let mut tmp0 = r.x;
-            tmp0.square_assign();
-
-            let mut tmp1 = r.y;
-            tmp1.square_assign();
-
-            let mut tmp2 = tmp1;
-            tmp2.square_assign();
-
-            let mut tmp3 = tmp1;
-            tmp3 += &r.x;
-            tmp3.square_assign();
-            tmp3 -= &tmp0;
-            tmp3 -= &tmp2;
-            tmp3.double_assign();
-
-            let mut tmp4 = tmp0;
-            tmp4.double_assign();
-            tmp4 += &tmp0;
+crate::impl_gt!(Gt, Fq12, Fr);
+crate::impl_miller_loop_components!(Bn256, G1, G1Affine, G2, G2Affine, Fq12, Gt, Fr);
 
-            let mut tmp6 = r.x;
-            tmp6 += &tmp4;
-
-            let mut tmp5 = tmp4;
-            tmp5.square_assign();
-
-            let mut zsquared = r.z;
-            zsquared.square_assign();
-
-            r.x = tmp5;
-            r.x -= &tmp3;
-            r.x -= &tmp3;
-
-            r.z += &r.y;
-            r.z.square_assign();
-            r.z -= &tmp1;
-            r.z -= &zsquared;
-
-            r.y = tmp3;
-            r.y -= &r.x;
-            r.y.mul_assign(&tmp4);
-
-            tmp2.double_assign();
-            tmp2.double_assign();
-            tmp2.double_assign();
-
-            r.y -= &tmp2;
-
-            // up to here everything was by algorithm, line 11
-            // use R instead of new T
-
-            // tmp3 is the first part of line 12
-            tmp3 = tmp4;
-            tmp3.mul_assign(&zsquared);
-            tmp3.double_assign();
-            tmp3 = tmp3.neg();
-
-            // tmp6 is from line 14
-            tmp6.square_assign();
-            tmp6 -= &tmp0;
-            tmp6 -= &tmp5;
-
-            tmp1.double_assign();
-            tmp1.double_assign();
-
-            tmp6 -= &tmp1;
-
-            // tmp0 is the first part of line 16
-            tmp0 = r.z;
-            tmp0.mul_assign(&zsquared);
-            tmp0.double_assign();
-
-            (tmp0, tmp3, tmp6)
-        }
-
-        fn addition_step(r: &mut G2, q: &G2Affine) -> (Fq2, Fq2, Fq2) {
-            // Adaptation of Algorithm 27, https://eprint.iacr.org/2010/354.pdf
-            let mut zsquared = r.z;
-            zsquared.square_assign();
-
-            let mut ysquared = q.y;
-            ysquared.square_assign();
-
-            // t0 corresponds to line 1
-            let mut t0 = zsquared;
-            t0.mul_assign(&q.x);
-
-            // t1 corresponds to lines 2 and 3
-            let mut t1 = q.y;
-            t1 += &r.z;
-            t1.square_assign();
-            t1 -= &ysquared;
-            t1 -= &zsquared;
-            t1.mul_assign(&zsquared);
-
-            // t2 corresponds to line 4
-            let mut t2 = t0;
-            t2 -= &r.x;
-
-            // t3 corresponds to line 5
-            let mut t3 = t2;
-            t3.square_assign();
-
-            // t4 corresponds to line 6
-            let mut t4 = t3;
-            t4.double_assign();
-            t4.double_assign();
-
-            // t5 corresponds to line 7
-            let mut t5 = t4;
-            t5.mul_assign(&t2);
-
-            // t6 corresponds to line 8
-            let mut t6 = t1;
-            t6 -= &r.y;
-            t6 -= &r.y;
-
-            // t9 corresponds to line 9
-            let mut t9 = t6;
-            t9.mul_assign(&q.x);
-
-            // corresponds to line 10
-            let mut t7 = t4;
-            t7.mul_assign(&r.x);
-
-            // corresponds to line 11, but assigns to r.x instead of T.x
-            r.x = t6;
-            r.x.square_assign();
-            r.x -= &t5;
-            r.x -= &t7;
-            r.x -= &t7;
-
-            // corresponds to line 12, but assigns to r.z instead of T.z
-            r.z += &t2;
-            r.z.square_assign();
-            r.z -= &zsquared;
-            r.z -= &t3;
-
-            // corresponds to line 13
-            let mut t10 = q.y;
-            t10 += &r.z;
-
-            // corresponds to line 14
-            let mut t8 = t7;
-            t8 -= &r.x;
-            t8.mul_assign(&t6);
-
-            // corresponds to line 15
-            t0 = r.y;
-            t0.mul_assign(&t5);
-            t0.double_assign();
-
-            // corresponds to line 12, but assigns to r.y instead of T.y
-            r.y = t8;
-            r.y -= &t0;
-
-            // corresponds to line 17
-            t10.square_assign();
-            t10 -= &ysquared;
-
-            let mut ztsquared = r.z;
-            ztsquared.square_assign();
-
-            t10 -= &ztsquared;
-
-            // corresponds to line 18
-            t9.double_assign();
-            t9 -= &t10;
-
-            // t10 = 2*Zt from Algo 27, line 19
-            t10 = r.z;
-            t10.double_assign();
-
-            // t1 = first multiplicator of line 21
-            t6 = t6.neg();
-
-            t1 = t6;
-            t1.double_assign();
-
-            // t9 corresponds to t9 from Algo 27
-            (t10, t1, t9)
-        }
-
-        let mut coeffs = vec![];
-        let mut r: G2 = q.into();
-
-        let mut negq = q;
-        negq = -negq;
-
-        for i in (1..SIX_U_PLUS_2_NAF.len()).rev() {
-            coeffs.push(doubling_step(&mut r));
-            let x = SIX_U_PLUS_2_NAF[i - 1];
-            match x {
-                1 => {
-                    coeffs.push(addition_step(&mut r, &q));
-                }
-                -1 => {
-                    coeffs.push(addition_step(&mut r, &negq));
-                }
-                _ => continue,
-            }
-        }
-
-        let mut q1 = q;
-
-        q1.x.c1 = q1.x.c1.neg();
-        q1.x.mul_assign(&FROBENIUS_COEFF_FQ6_C1[1]);
-
-        q1.y.c1 = q1.y.c1.neg();
-        q1.y.mul_assign(&XI_TO_Q_MINUS_1_OVER_2);
-
-        coeffs.push(addition_step(&mut r, &q1));
-
-        let mut minusq2 = q;
-        minusq2.x.mul_assign(&FROBENIUS_COEFF_FQ6_C1[2]);
-
-        coeffs.push(addition_step(&mut r, &minusq2));
-
-        G2Prepared {
-            coeffs,
-            infinity: false,
-        }
-    }
-}
-
-impl From<G2Affine> for G2Prepared {
-    fn from(q: G2Affine) -> G2Prepared {
-        G2Prepared::from_affine(q)
-    }
-}
+impl MillerLoopResult for Fq12 {
+    type Gt = Gt;
 
-impl MillerLoopResult for Gt {
-    type Gt = Self;
-    // pub fn final_exponentiation(r: &Fq12) -> CtOption<Fq12> {
-    fn final_exponentiation(&self) -> Gt {
+    fn final_exponentiation(&self) -> Self::Gt {
         fn exp_by_x(f: &mut Fq12) {
-            let x = BN_X;
+            let x = super::BN_X;
             let mut res = Fq12::one();
             for i in (0..64).rev() {
                 res.cyclotomic_square();
@@ -471,10 +36,11 @@ impl MillerLoopResult for Gt {
             *f = res;
         }
 
-        let r = self.0;
-        let mut f1 = self.0;
+        let r = *self;
+        let mut f1 = *self;
         f1.conjugate();
 
+        use ff::Field;
         Gt(r.invert()
             .map(|mut f2| {
                 let mut r = f1;
@@ -560,283 +126,93 @@ impl MillerLoopResult for Gt {
     }
 }
 
-pub fn multi_miller_loop(terms: &[(&G1Affine, &G2Prepared)]) -> Gt {
-    let mut pairs = vec![];
-    for &(p, q) in terms {
-        if !bool::from(p.is_identity()) && !q.is_zero() {
-            pairs.push((p, q.coeffs.iter()));
-        }
-    }
-
-    // Final steps of the line function on prepared coefficients
-    fn ell(f: &mut Fq12, coeffs: &(Fq2, Fq2, Fq2), p: &G1Affine) {
-        let mut c0 = coeffs.0;
-        let mut c1 = coeffs.1;
-
-        c0.c0.mul_assign(&p.y);
-        c0.c1.mul_assign(&p.y);
-
-        c1.c0.mul_assign(&p.x);
-        c1.c1.mul_assign(&p.x);
-
-        // Sparse multiplication in Fq12
-        f.mul_by_034(&c0, &c1, &coeffs.2);
-    }
+pub fn multi_miller_loop(terms: &[(&G1Affine, &G2Affine)]) -> Fq12 {
+    let terms = terms
+        .iter()
+        .filter_map(|&(p, q)| {
+            if bool::from(p.is_identity()) || bool::from(q.is_identity()) {
+                None
+            } else {
+                Some((*p, *q))
+            }
+        })
+        .collect::<Vec<_>>();
 
     let mut f = Fq12::one();
+    let mut r = terms.iter().map(|(_, q)| q.to_curve()).collect::<Vec<_>>();
 
-    for i in (1..SIX_U_PLUS_2_NAF.len()).rev() {
-        if i != SIX_U_PLUS_2_NAF.len() - 1 {
-            f.square_assign();
-        }
-        for &mut (p, ref mut coeffs) in &mut pairs {
-            ell(&mut f, coeffs.next().unwrap(), p);
+    for (i, x) in super::SIX_U_PLUS_2_NAF.iter().rev().skip(1).enumerate() {
+        (i != 0).then(|| f.square_assign());
+
+        for ((p, _), r) in terms.iter().zip(r.iter_mut()) {
+            double(&mut f, r, p);
         }
-        let x = SIX_U_PLUS_2_NAF[i - 1];
+
         match x {
-            1 => {
-                for &mut (p, ref mut coeffs) in &mut pairs {
-                    ell(&mut f, coeffs.next().unwrap(), p);
-                }
-            }
-            -1 => {
-                for &mut (p, ref mut coeffs) in &mut pairs {
-                    ell(&mut f, coeffs.next().unwrap(), p);
+            &val @ (1 | -1) => {
+                for ((p, q), r) in terms.iter().zip(r.iter_mut()) {
+                    if val == 1 {
+                        add(&mut f, r, q, p);
+                    } else {
+                        add(&mut f, r, &q.neg(), p);
+                    }
                 }
             }
             _ => continue,
         }
     }
 
-    for &mut (p, ref mut coeffs) in &mut pairs {
-        ell(&mut f, coeffs.next().unwrap(), p);
-    }
-
-    for &mut (p, ref mut coeffs) in &mut pairs {
-        ell(&mut f, coeffs.next().unwrap(), p);
+    const XI_TO_Q_MINUS_1_OVER_2: Fq2 = Fq2 {
+        c0: Fq([
+            0xe4bbdd0c2936b629,
+            0xbb30f162e133bacb,
+            0x31a9d1b6f9645366,
+            0x253570bea500f8dd,
+        ]),
+        c1: Fq([
+            0xa1d77ce45ffe77c7,
+            0x07affd117826d1db,
+            0x6d16bd27bb7edc6b,
+            0x2c87200285defecc,
+        ]),
+    };
+
+    for ((p, q), r) in terms.iter().zip(r.iter_mut()) {
+        let mut q1: G2Affine = *q;
+        q1.x.conjugate();
+        q1.x.mul_assign(&FROBENIUS_COEFF_FQ6_C1[1]);
+        q1.y.conjugate();
+        q1.y.mul_assign(&XI_TO_Q_MINUS_1_OVER_2);
+        add(&mut f, r, &q1, p);
     }
 
-    for &mut (_p, ref mut coeffs) in &mut pairs {
-        assert_eq!(coeffs.next(), None);
+    for ((p, q), r) in terms.iter().zip(r.iter_mut()) {
+        let mut minusq2: G2Affine = *q;
+        minusq2.x.mul_assign(&FROBENIUS_COEFF_FQ6_C1[2]);
+        add(&mut f, r, &minusq2, p);
     }
 
-    Gt(f)
+    f
 }
 
-pub fn pairing(g1: &G1Affine, g2: &G2Affine) -> Gt {
-    let g2 = G2Prepared::from_affine(*g2);
-    let terms: &[(&G1Affine, &G2Prepared)] = &[(g1, &g2)];
-    let u = multi_miller_loop(terms);
-    u.final_exponentiation()
+// Final steps of the line function on prepared coefficients
+fn ell(f: &mut Fq12, coeffs: &(Fq2, Fq2, Fq2), p: &G1Affine) {
+    let mut c0 = coeffs.0;
+    let mut c1 = coeffs.1;
+    c0.c0.mul_assign(&p.y);
+    c0.c1.mul_assign(&p.y);
+    c1.c0.mul_assign(&p.x);
+    c1.c1.mul_assign(&p.x);
+    Fq12::mul_by_034(f, &c0, &c1, &coeffs.2);
 }
 
-#[derive(Clone, Debug)]
-pub struct Bn256;
-
-impl Engine for Bn256 {
-    type Fr = Fr;
-    type G1 = G1;
-    type G1Affine = G1Affine;
-    type G2 = G2;
-    type G2Affine = G2Affine;
-    type Gt = Gt;
-
-    fn pairing(p: &Self::G1Affine, q: &Self::G2Affine) -> Self::Gt {
-        pairing(p, q)
-    }
-}
-
-impl MultiMillerLoop for Bn256 {
-    type G2Prepared = G2Prepared;
-    type Result = Gt;
-
-    fn multi_miller_loop(terms: &[(&Self::G1Affine, &Self::G2Prepared)]) -> Self::Result {
-        multi_miller_loop(terms)
-    }
-}
-
-#[cfg(test)]
-use rand::SeedableRng;
 #[cfg(test)]
-use rand_xorshift::XorShiftRng;
-
-#[test]
-fn test_pairing() {
-    use crate::ff::Field;
-    let g1 = G1::generator();
-    let mut g2 = G2::generator();
-    g2 = g2.double();
-    let pair12 = Bn256::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    let mut g1 = G1::generator();
-    let g2 = G2::generator();
-    g1 = g1.double();
-    let pair21 = Bn256::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    assert_eq!(pair12, pair21);
-
-    let g1 = G1::generator();
-    let mut g2 = G2::generator();
-    g2 = g2.double().double();
-    let pair12 = Bn256::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    let mut g1 = G1::generator();
-    let mut g2 = G2::generator();
-    g1 = g1.double();
-    g2 = g2.double();
-    let pair21 = Bn256::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    assert_eq!(pair12, pair21);
-
-    let mut rng = XorShiftRng::from_seed([
-        0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-    for _ in 0..1000 {
-        let a = Fr::random(&mut rng);
-        let b = Fr::random(&mut rng);
-
-        let mut g1 = G1::generator();
-        g1.mul_assign(a);
-
-        let mut g2 = G2::generator();
-        g1.mul_assign(b);
-
-        let pair_ab = Bn256::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-        g1 = G1::generator();
-        g1.mul_assign(b);
-
-        g2 = G2::generator();
-        g1.mul_assign(a);
-
-        let pair_ba = Bn256::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-        assert_eq!(pair_ab, pair_ba);
-    }
-}
-
-#[test]
-fn random_bilinearity_tests() {
-    use crate::ff::Field;
-    let mut rng = XorShiftRng::from_seed([
-        0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    for _ in 0..1000 {
-        let mut a = G1::generator();
-        let ka = Fr::random(&mut rng);
-        a.mul_assign(ka);
-
-        let mut b = G2::generator();
-        let kb = Fr::random(&mut rng);
-        b.mul_assign(kb);
-
-        let c = Fr::random(&mut rng);
-        let d = Fr::random(&mut rng);
-
-        let mut ac = a;
-        ac.mul_assign(c);
-
-        let mut ad = a;
-        ad.mul_assign(d);
-
-        let mut bc = b;
-        bc.mul_assign(c);
-
-        let mut bd = b;
-        bd.mul_assign(d);
-
-        let acbd = Bn256::pairing(&G1Affine::from(ac), &G2Affine::from(bd));
-        let adbc = Bn256::pairing(&G1Affine::from(ad), &G2Affine::from(bc));
-
-        let mut cd = c;
-        cd.mul_assign(&d);
-
-        cd *= Fr([1, 0, 0, 0]);
-
-        let abcd = Gt(Bn256::pairing(&G1Affine::from(a), &G2Affine::from(b))
-            .0
-            .pow_vartime(cd.0));
-
-        assert_eq!(acbd, adbc);
-        assert_eq!(acbd, abcd);
-    }
-}
-
-#[test]
-pub fn engine_tests() {
-    use crate::ff::Field;
-    let mut rng = XorShiftRng::from_seed([
-        0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    for _ in 0..10 {
-        let a = G1Affine::from(G1::random(&mut rng));
-        let b = G2Affine::from(G2::random(&mut rng));
-
-        assert!(a.pairing_with(&b) == b.pairing_with(&a));
-        assert!(a.pairing_with(&b) == pairing(&a, &b));
-    }
-
-    for _ in 0..1000 {
-        let z1 = G1Affine::identity();
-        let z2 = G2Prepared::from(G2Affine::identity());
-
-        let a = G1Affine::from(G1::random(&mut rng));
-        let b = G2Prepared::from(G2Affine::from(G2::random(&mut rng)));
-        let c = G1Affine::from(G1::random(&mut rng));
-        let d = G2Prepared::from(G2Affine::from(G2::random(&mut rng)));
-
-        assert_eq!(
-            Fq12::ONE,
-            multi_miller_loop(&[(&z1, &b)]).final_exponentiation().0,
-        );
-
-        assert_eq!(
-            Fq12::ONE,
-            multi_miller_loop(&[(&a, &z2)]).final_exponentiation().0,
-        );
-
-        assert_eq!(
-            multi_miller_loop(&[(&z1, &b), (&c, &d)]).final_exponentiation(),
-            multi_miller_loop(&[(&a, &z2), (&c, &d)]).final_exponentiation(),
-        );
-
-        assert_eq!(
-            multi_miller_loop(&[(&a, &b), (&z1, &d)]).final_exponentiation(),
-            multi_miller_loop(&[(&a, &b), (&c, &z2)]).final_exponentiation(),
-        );
-    }
-}
-
-#[test]
-fn random_miller_loop_tests() {
-    let mut rng = XorShiftRng::from_seed([
-        0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    // Exercise a double miller loop
-    for _ in 0..1000 {
-        let a = G1Affine::from(G1::random(&mut rng));
-        let b = G2Affine::from(G2::random(&mut rng));
-        let c = G1Affine::from(G1::random(&mut rng));
-        let d = G2Affine::from(G2::random(&mut rng));
-
-        let ab = pairing(&a, &b);
-        let cd = pairing(&c, &d);
-
-        let mut abcd = ab;
-        abcd = Gt(abcd.0 * cd.0);
-
-        let b = G2Prepared::from(b);
-        let d = G2Prepared::from(d);
-
-        let abcd_with_double_loop = multi_miller_loop(&[(&a, &b), (&c, &d)]).final_exponentiation();
-
-        assert_eq!(abcd, abcd_with_double_loop);
-    }
+mod test {
+    use super::super::{Bn256, Fr, G1, G2};
+    use super::{multi_miller_loop, Fq12, G1Affine, G2Affine, Gt};
+    use ff::Field;
+    use group::{prime::PrimeCurveAffine, Curve, Group};
+    use pairing::{Engine, MillerLoopResult, PairingCurveAffine};
+    use rand_core::OsRng;
+    crate::test_pairing!(Bn256, G1, G1Affine, G2, G2Affine, Fq12, Gt, Fr);
 }
diff --git a/src/bn256/fq.rs b/src/bn256/fq.rs
index 153c7213..791c2ada 100644
--- a/src/bn256/fq.rs
+++ b/src/bn256/fq.rs
@@ -1,15 +1,9 @@
+use crate::ff_ext::ExtField;
 use core::convert::TryInto;
 use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_from_u64, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     bn256_base,
     Fq,
@@ -20,13 +14,23 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fq);
-impl_binops_calls!(Fq);
-impl_binops_additive!(Fq, Fq);
-impl_binops_multiplicative!(Fq, Fq);
-field_bits!(Fq);
-serialize_deserialize_primefield!(Fq);
-impl_from_u64!(Fq);
+crate::extend_field_legendre!(Fq);
+crate::impl_binops_calls!(Fq);
+crate::impl_binops_additive!(Fq, Fq);
+crate::impl_binops_multiplicative!(Fq, Fq);
+crate::field_bits!(Fq);
+crate::serialize_deserialize_primefield!(Fq);
+crate::impl_from_u64!(Fq);
+
+use ff::Field;
+const NEGATIVE_ONE: Fq = Fq::ZERO.sub_const(&Fq::ONE);
+impl ExtField for Fq {
+    const NON_RESIDUE: Self = NEGATIVE_ONE;
+    fn mul_by_nonresidue(&self) -> Self {
+        self.neg()
+    }
+    fn frobenius_map(&mut self, _: usize) {}
+}
 
 #[cfg(test)]
 mod test {
@@ -41,4 +45,11 @@ mod test {
     crate::field_testing_suite!(Fq, "sqrt");
     crate::field_testing_suite!(Fq, "zeta");
     crate::field_testing_suite!(Fq, "from_uniform_bytes", 64, 48);
+    #[test]
+    fn test_fq_mul_nonresidue() {
+        let e = Fq::random(rand_core::OsRng);
+        let a0 = e.mul_by_nonresidue();
+        let a1 = e * Fq::NON_RESIDUE;
+        assert_eq!(a0, a1);
+    }
 }
diff --git a/src/bn256/fq12.rs b/src/bn256/fq12.rs
index 7dd08c45..d10c73d7 100644
--- a/src/bn256/fq12.rs
+++ b/src/bn256/fq12.rs
@@ -1,204 +1,42 @@
 use super::fq::Fq;
 use super::fq2::Fq2;
 use super::fq6::Fq6;
-use crate::impl_tower2_common;
-use core::ops::{Add, Neg, Sub};
-use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
+use crate::ff_ext::{
+    quadratic::{QuadExtField, QuadExtFieldArith, QuadSparseMul},
+    ExtField,
+};
 
 /// -GAMMA is a quadratic non-residue in Fp6. Fp12 = Fp6[X]/(X^2 + GAMMA)
 /// We introduce the variable w such that w^2 = -GAMMA
 // GAMMA = - v
 /// An element of Fq12, represented by c0 + c1 * w.
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_calls, impl_binops_multiplicative, impl_binops_multiplicative_mixed,
-    impl_sub_binop_specify_output, impl_sum_prod,
-};
-impl_binops_additive!(Fq12, Fq12);
-impl_binops_multiplicative!(Fq12, Fq12);
-impl_tower2_common!(Fq6, Fq12);
-impl_binops_calls!(Fq12);
-impl_sum_prod!(Fq12);
-
-impl Fq12 {
-    pub fn mul_assign(&mut self, other: &Self) {
-        let t0 = self.c0 * other.c0;
-        let mut t1 = self.c1 * other.c1;
-        let t2 = other.c0 + other.c1;
-
-        self.c1 += &self.c0;
-        self.c1 *= &t2;
-        self.c1 -= &t0;
-        self.c1 -= &t1;
-
-        t1.mul_by_nonresidue();
-        self.c0 = t0 + t1;
-    }
+pub type Fq12 = QuadExtField<Fq6>;
 
-    pub fn square_assign(&mut self) {
-        let mut ab = self.c0 * self.c1;
-
-        let c0c1 = self.c0 + self.c1;
+impl QuadExtFieldArith for Fq12 {
+    type Base = Fq6;
+}
 
-        let mut c0 = self.c1;
-        c0.mul_by_nonresidue();
-        c0 += &self.c0;
-        c0 *= &c0c1;
-        c0 -= &ab;
-        self.c1 = ab;
-        self.c1 += &ab;
-        ab.mul_by_nonresidue();
-        c0 -= &ab;
-        self.c0 = c0;
-    }
+impl QuadSparseMul for Fq12 {
+    type Base = Fq2;
+}
 
-    #[inline(always)]
-    pub fn conjugate(&mut self) {
-        self.c1 = -self.c1;
-    }
+impl ExtField for Fq12 {
+    const NON_RESIDUE: Self = Fq12::zero(); // no needs
 
-    pub fn frobenius_map(&mut self, power: usize) {
+    fn frobenius_map(&mut self, power: usize) {
         self.c0.frobenius_map(power);
         self.c1.frobenius_map(power);
-
         self.c1.c0.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
         self.c1.c1.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
         self.c1.c2.mul_assign(&FROBENIUS_COEFF_FQ12_C1[power % 12]);
     }
-
-    pub fn mul_by_014(&mut self, c0: &Fq2, c1: &Fq2, c4: &Fq2) {
-        let mut aa = self.c0;
-        aa.mul_by_01(c0, c1);
-        let mut bb = self.c1;
-        bb.mul_by_1(c4);
-        let o = c1 + c4;
-        self.c1 += &self.c0;
-        self.c1.mul_by_01(c0, &o);
-        self.c1 -= &aa;
-        self.c1 -= &bb;
-        self.c0 = bb;
-        self.c0.mul_by_nonresidue();
-        self.c0 += &aa;
-    }
-
-    pub fn mul_by_034(&mut self, c0: &Fq2, c3: &Fq2, c4: &Fq2) {
-        let t0 = Fq6 {
-            c0: self.c0.c0 * c0,
-            c1: self.c0.c1 * c0,
-            c2: self.c0.c2 * c0,
-        };
-        let mut t1 = self.c1;
-        t1.mul_by_01(c3, c4);
-        let o = c0 + c3;
-        let mut t2 = self.c0 + self.c1;
-        t2.mul_by_01(&o, c4);
-        t2 -= t0;
-        self.c1 = t2 - t1;
-        t1.mul_by_nonresidue();
-        self.c0 = t0 + t1;
-    }
-
-    pub fn invert(&self) -> CtOption<Self> {
-        let mut c0s = self.c0;
-        c0s.square_assign();
-        let mut c1s = self.c1;
-        c1s.square_assign();
-        c1s.mul_by_nonresidue();
-        c0s -= &c1s;
-
-        c0s.invert().map(|t| {
-            let mut tmp = Fq12 { c0: t, c1: t };
-            tmp.c0.mul_assign(&self.c0);
-            tmp.c1.mul_assign(&self.c1);
-            tmp.c1 = tmp.c1.neg();
-
-            tmp
-        })
-    }
-
-    pub fn cyclotomic_square(&mut self) {
-        fn fp4_square(c0: &mut Fq2, c1: &mut Fq2, a0: &Fq2, a1: &Fq2) {
-            let t0 = a0.square();
-            let t1 = a1.square();
-            let mut t2 = t1;
-            t2.mul_by_nonresidue();
-            *c0 = t2 + t0;
-            t2 = a0 + a1;
-            t2.square_assign();
-            t2 -= t0;
-            *c1 = t2 - t1;
-        }
-
-        let mut t3 = Fq2::zero();
-        let mut t4 = Fq2::zero();
-        let mut t5 = Fq2::zero();
-        let mut t6 = Fq2::zero();
-
-        fp4_square(&mut t3, &mut t4, &self.c0.c0, &self.c1.c1);
-        let mut t2 = t3 - self.c0.c0;
-        t2.double_assign();
-        self.c0.c0 = t2 + t3;
-
-        t2 = t4 + self.c1.c1;
-        t2.double_assign();
-        self.c1.c1 = t2 + t4;
-
-        fp4_square(&mut t3, &mut t4, &self.c1.c0, &self.c0.c2);
-        fp4_square(&mut t5, &mut t6, &self.c0.c1, &self.c1.c2);
-
-        t2 = t3 - self.c0.c1;
-        t2.double_assign();
-        self.c0.c1 = t2 + t3;
-        t2 = t4 + self.c1.c2;
-        t2.double_assign();
-        self.c1.c2 = t2 + t4;
-        t3 = t6;
-        t3.mul_by_nonresidue();
-        t2 = t3 + self.c1.c0;
-        t2.double_assign();
-        self.c1.c0 = t2 + t3;
-        t2 = t5 - self.c0.c2;
-        t2.double_assign();
-        self.c0.c2 = t2 + t5;
-    }
 }
 
-#[cfg(test)]
-impl ff::Field for Fq12 {
-    const ZERO: Self = Self::zero();
-    const ONE: Self = Self::one();
-
-    fn random(mut rng: impl rand_core::RngCore) -> Self {
-        Fq12 {
-            c0: Fq6::random(&mut rng),
-            c1: Fq6::random(&mut rng),
-        }
-    }
-
-    fn is_zero(&self) -> Choice {
-        self.c0.is_zero() & self.c1.is_zero()
-    }
-
-    fn square(&self) -> Self {
-        self.square()
-    }
-
-    fn double(&self) -> Self {
-        self.double()
-    }
-
-    fn sqrt(&self) -> CtOption<Self> {
-        unimplemented!()
-    }
-
-    fn sqrt_ratio(_num: &Self, _div: &Self) -> (Choice, Self) {
-        unimplemented!()
-    }
-
-    fn invert(&self) -> CtOption<Self> {
-        self.invert()
-    }
-}
+crate::impl_binops_additive!(Fq12, Fq12);
+crate::impl_binops_multiplicative!(Fq12, Fq12);
+crate::impl_binops_calls!(Fq12);
+crate::impl_sum_prod!(Fq12);
+crate::impl_cyclotomic_square!(Fq2, Fq12);
 
 // non_residue^((modulus^i-1)/6) for i=0,...,11
 pub const FROBENIUS_COEFF_FQ12_C1: [Fq2; 12] = [
@@ -360,7 +198,7 @@ mod test {
     use super::*;
     crate::field_testing_suite!(Fq12, "field_arithmetic");
     // extension field-specific
-    crate::field_testing_suite!(Fq12, "f12_tests", Fq6, Fq2);
+    crate::field_testing_suite!(Fq12, "quadratic_sparse_mul", Fq6, Fq2);
     crate::field_testing_suite!(
         Fq12,
         "frobenius",
@@ -368,11 +206,6 @@ mod test {
         //  ϕ: E → E
         //  (x, y) ↦ (x^p, y^p)
         // p: modulus of base field (Here, Fq::MODULUS)
-        [
-            0x3c208c16d87cfd47,
-            0x97816a916871ca8d,
-            0xb85045b68181585d,
-            0x30644e72e131a029,
-        ]
+        Fq::MODULUS_LIMBS
     );
 }
diff --git a/src/bn256/fq2.rs b/src/bn256/fq2.rs
index ad59f424..3eb4f0ee 100644
--- a/src/bn256/fq2.rs
+++ b/src/bn256/fq2.rs
@@ -1,185 +1,72 @@
 use super::fq::Fq;
 use crate::ff::{Field, FromUniformBytes, PrimeField, WithSmallOrderMulGroup};
-use crate::ff_ext::Legendre;
+use crate::ff_ext::quadratic::{QuadExtField, QuadExtFieldArith, SQRT};
+use crate::ff_ext::{ExtField, Legendre};
 use core::convert::TryInto;
-use core::ops::{Add, Neg, Sub};
-use rand::RngCore;
 use std::cmp::Ordering;
-use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
-
-#[cfg(feature = "derive_serde")]
-use serde::{Deserialize, Serialize};
-
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_calls, impl_binops_multiplicative, impl_binops_multiplicative_mixed,
-    impl_sub_binop_specify_output, impl_sum_prod, impl_tower2, impl_tower2_common,
-};
-impl_tower2_common!(Fq, Fq2, serde);
-impl_tower2!(Fq, Fq2, ReprFq2);
-impl_binops_additive!(Fq2, Fq2);
-impl_binops_multiplicative!(Fq2, Fq2);
-impl_binops_calls!(Fq2);
-impl_sum_prod!(Fq2);
-
-impl Fq2 {
-    pub fn double_assign(&mut self) {
-        self.c0 = self.c0.double();
-        self.c1 = self.c1.double();
-    }
-
-    pub fn mul_assign(&mut self, other: &Self) {
-        let mut t0 = self.c0 + self.c1;
-        let mut t1 = self.c0 * other.c0;
-        let t2 = self.c1 * other.c1;
-
-        self.c0 = t1 - t2;
-        self.c1 = other.c0 + other.c1;
-        t1 += t2;
-        t0 *= self.c1;
-        self.c1 = t0 - t1;
-    }
-
-    pub fn square_assign(&mut self) {
-        let ab = self.c0 * self.c1;
-        let c0c1 = self.c0 + self.c1;
-        let mut c0 = -self.c1;
-        c0 += self.c0;
-        c0 *= c0c1;
-        self.c1 = ab.double();
-        self.c0 = c0;
-    }
-
-    // conjugate by negating c1
-    pub fn conjugate(&mut self) {
-        self.c1 = -self.c1;
+use subtle::{Choice, CtOption};
+
+crate::impl_binops_additive!(Fq2, Fq2);
+crate::impl_binops_multiplicative!(Fq2, Fq2);
+crate::impl_binops_calls!(Fq2);
+crate::impl_sum_prod!(Fq2);
+crate::impl_tower2!(Fq, Fq2);
+crate::impl_tower2_from_uniform_bytes!(Fq, Fq2, 96);
+
+pub type Fq2 = QuadExtField<Fq>;
+impl QuadExtFieldArith for Fq2 {
+    type Base = Fq;
+    const SQRT: SQRT<Fq> = SQRT::Algorithm9 {
+        q_minus_3_over_4: &[
+            0x4f082305b61f3f51,
+            0x65e05aa45a1c72a3,
+            0x6e14116da0605617,
+            0x0c19139cb84c680a,
+        ],
+        q_minus_1_over_2: &[
+            0x9e10460b6c3e7ea3,
+            0xcbc0b548b438e546,
+            0xdc2822db40c0ac2e,
+            0x183227397098d014,
+        ],
+    };
+
+    fn square_assign(el: &mut QuadExtField<Self::Base>) {
+        let a = el.c0 + el.c1;
+        let b = el.c0 - el.c1;
+        let c = el.c0.double();
+        el.c0 = a * b;
+        el.c1 = c * el.c1;
     }
+}
 
-    pub fn frobenius_map(&mut self, power: usize) {
-        if power % 2 != 0 {
-            self.conjugate()
-        }
-    }
+impl ExtField for Fq2 {
+    const NON_RESIDUE: Self = Fq2::new(Fq::from_raw([9u64, 0, 0, 0]), Fq::ONE);
 
-    /// Multiply this element by quadratic nonresidue 9 + u.
-    pub fn mul_by_nonresidue(&mut self) {
+    fn mul_by_nonresidue(&self) -> Self {
         // (xu+y)(u+9) = (9x+y)u+(9y-x)
         let t0 = self.c0;
         let t1 = self.c1;
-
         // 8*x*i + 8*y
-        *self = self.double();
-        *self = self.double();
-        *self = self.double();
-
-        // 9*y
-        self.c0 += &t0;
-        // (9*y - x)
-        self.c0 -= &t1;
-
-        // (9*x)u
-        self.c1 += &t1;
-        // (9*x + y)
-        self.c1 += &t0;
-    }
-
-    pub fn invert(&self) -> CtOption<Self> {
-        let mut t1 = self.c1;
-        t1 = t1.square();
-        let mut t0 = self.c0;
-        t0 = t0.square();
-        t0 += &t1;
-        t0.invert().map(|t| {
-            let mut tmp = Fq2 {
-                c0: self.c0,
-                c1: self.c1,
-            };
-            tmp.c0 *= &t;
-            tmp.c1 *= &t;
-            tmp.c1 = -tmp.c1;
-
-            tmp
-        })
-    }
-
-    /// Norm of Fq2 as extension field in i over Fq
-    #[inline]
-    fn norm(&self) -> Fq {
-        let mut t0 = self.c0;
-        let mut t1 = self.c1;
-        t0 = t0.square();
-        t1 = t1.square();
-        t1 + t0
-    }
-
-    pub fn sqrt(&self) -> CtOption<Self> {
-        // Algorithm 9, https://eprint.iacr.org/2012/685.pdf
-
-        if self.is_zero().into() {
-            CtOption::new(Self::ZERO, Choice::from(1))
-        } else {
-            // a1 = self^((q - 3) / 4)
-            // 0xc19139cb84c680a6e14116da060561765e05aa45a1c72a34f082305b61f3f51
-            let u: [u64; 4] = [
-                0x4f082305b61f3f51,
-                0x65e05aa45a1c72a3,
-                0x6e14116da0605617,
-                0x0c19139cb84c680a,
-            ];
-            let mut a1 = self.pow(u);
-            let mut alpha = a1;
-
-            alpha.square_assign();
-            alpha.mul_assign(self);
-            let mut a0 = alpha;
-            a0.frobenius_map(1);
-            a0.mul_assign(&alpha);
-
-            const NEGATIVE_ONE: Fq = Fq::ZERO.sub_const(&Fq::ONE);
-            let neg1 = Fq2 {
-                c0: NEGATIVE_ONE,
-                c1: Fq::zero(),
-            };
-
-            if a0 == neg1 {
-                CtOption::new(a0, Choice::from(0))
-            } else {
-                a1.mul_assign(self);
-
-                if alpha == neg1 {
-                    a1.mul_assign(&Fq2 {
-                        c0: Fq::zero(),
-                        c1: Fq::one(),
-                    });
-                } else {
-                    alpha += &Fq2::ONE;
-                    // alpha = alpha^((q - 1) / 2)
-                    // 0x183227397098d014dc2822db40c0ac2ecbc0b548b438e5469e10460b6c3e7ea3
-                    let u: [u64; 4] = [
-                        0x9e10460b6c3e7ea3,
-                        0xcbc0b548b438e546,
-                        0xdc2822db40c0ac2e,
-                        0x183227397098d014,
-                    ];
-                    alpha = alpha.pow(u);
-                    a1.mul_assign(&alpha);
-                }
-                CtOption::new(a1, Choice::from(1))
-            }
+        let t = self.double().double().double();
+        Self {
+            // 9*y
+            c0: t.c0 + t0 - t1,
+            // (9*x + y)
+            c1: t.c1 + t0 + t1,
         }
     }
-}
 
-impl FromUniformBytes<96> for Fq2 {
-    fn from_uniform_bytes(bytes: &[u8; 96]) -> Self {
-        let c0: [u8; 48] = bytes[..48].try_into().unwrap();
-        let c1: [u8; 48] = bytes[48..].try_into().unwrap();
-        Self::new(Fq::from_uniform_bytes(&c1), Fq::from_uniform_bytes(&c0))
+    fn frobenius_map(&mut self, power: usize) {
+        if power % 2 != 0 {
+            self.conjugate();
+        }
     }
 }
 
 #[cfg(test)]
 mod test {
+
     use super::*;
     crate::field_testing_suite!(Fq2, "field_arithmetic");
     crate::field_testing_suite!(Fq2, "conversion");
@@ -196,62 +83,14 @@ mod test {
         //  ϕ: E → E
         //  (x, y) ↦ (x^p, y^p)
         // p: modulus of base field (Here, Fq::MODULUS)
-        [
-            0x3c208c16d87cfd47,
-            0x97816a916871ca8d,
-            0xb85045b68181585d,
-            0x30644e72e131a029,
-        ]
+        Fq::MODULUS_LIMBS
     );
 
-    #[test]
-    fn test_fq2_squaring() {
-        let mut a = Fq2 {
-            c0: Fq::one(),
-            c1: Fq::one(),
-        }; // u + 1
-        a.square_assign();
-        assert_eq!(
-            a,
-            Fq2 {
-                c0: Fq::zero(),
-                c1: Fq::one() + Fq::one(),
-            }
-        ); // 2u
-
-        let mut a = Fq2 {
-            c0: Fq::zero(),
-            c1: Fq::one(),
-        }; // u
-        a.square_assign();
-        assert_eq!(a, {
-            let neg1 = -Fq::one();
-            Fq2 {
-                c0: neg1,
-                c1: Fq::zero(),
-            }
-        }); // -1
-    }
-
     #[test]
     fn test_fq2_mul_nonresidue() {
-        let mut rng = XorShiftRng::from_seed([
-            0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06,
-            0xbc, 0xe5,
-        ]);
-        let nine = Fq::one().double().double().double() + Fq::one();
-        let nqr = Fq2 {
-            c0: nine,
-            c1: Fq::one(),
-        };
-
-        for _ in 0..1000 {
-            let mut a = Fq2::random(&mut rng);
-            let mut b = a;
-            a.mul_by_nonresidue();
-            b.mul_assign(&nqr);
-
-            assert_eq!(a, b);
-        }
+        let e = Fq2::random(rand_core::OsRng);
+        let a0 = e.mul_by_nonresidue();
+        let a1 = e * Fq2::NON_RESIDUE;
+        assert_eq!(a0, a1);
     }
 }
diff --git a/src/bn256/fq6.rs b/src/bn256/fq6.rs
index f784bdf0..f2e63e89 100644
--- a/src/bn256/fq6.rs
+++ b/src/bn256/fq6.rs
@@ -1,257 +1,45 @@
 use super::fq::Fq;
 use super::fq2::Fq2;
-use crate::{ff::Field, impl_tower6};
-use rand::RngCore;
-use subtle::CtOption;
+use crate::ff_ext::{
+    cubic::{CubicExtField, CubicExtFieldArith, CubicSparseMul},
+    ExtField,
+};
+use ff::Field;
 
-/// -BETA is a cubic non-residue in Fp2. Fp6 = Fp2[X]/(X^3 + BETA)
-/// We introduce the variable v such that v^3 = -BETA
+// -BETA is a cubic non-residue in Fp2. Fp6 = Fp2[X]/(X^3 + BETA)
+// We introduce the variable v such that v^3 = -BETA
 // BETA = - (u + 9)
-/// An element of Fq6, represented by c0 + c1 * v + c2 * v^2.
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_calls, impl_binops_multiplicative, impl_binops_multiplicative_mixed,
-    impl_sub_binop_specify_output,
-};
-impl_tower6!(Fq, Fq2, Fq6);
-impl_binops_additive!(Fq6, Fq6);
-impl_binops_multiplicative!(Fq6, Fq6);
-impl_binops_calls!(Fq6);
-#[cfg(test)]
+// An element of Fq6, represented by c0 + c1 * v + c2 * v^2.Â
+crate::impl_binops_additive!(Fq6, Fq6);
+crate::impl_binops_multiplicative!(Fq6, Fq6);
+crate::impl_binops_calls!(Fq6);
 crate::impl_sum_prod!(Fq6);
+pub type Fq6 = CubicExtField<Fq2>;
 
-impl Fq6 {
-    pub fn mul_assign(&mut self, other: &Self) {
-        let mut a_a = self.c0;
-        let mut b_b = self.c1;
-        let mut c_c = self.c2;
-        a_a *= &other.c0;
-        b_b *= &other.c1;
-        c_c *= &other.c2;
-
-        let mut t1 = other.c1;
-        t1 += &other.c2;
-        {
-            let mut tmp = self.c1;
-            tmp += &self.c2;
-
-            t1 *= &tmp;
-            t1 -= &b_b;
-            t1 -= &c_c;
-            t1.mul_by_nonresidue();
-            t1 += &a_a;
-        }
-
-        let mut t3 = other.c0;
-        t3 += &other.c2;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c2;
-
-            t3 *= &tmp;
-            t3 -= &a_a;
-            t3 += &b_b;
-            t3 -= &c_c;
-        }
-
-        let mut t2 = other.c0;
-        t2 += &other.c1;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c1;
-
-            t2 *= &tmp;
-            t2 -= &a_a;
-            t2 -= &b_b;
-            c_c.mul_by_nonresidue();
-            t2 += &c_c;
-        }
-
-        self.c0 = t1;
-        self.c1 = t2;
-        self.c2 = t3;
-    }
-
-    pub fn square_assign(&mut self) {
-        // s0 = a^2
-        let mut s0 = self.c0;
-        s0.square_assign();
-        // s1 = 2ab
-        let mut ab = self.c0;
-        ab *= &self.c1;
-        let mut s1 = ab;
-        s1.double_assign();
-        // s2 = (a - b + c)^2
-        let mut s2 = self.c0;
-        s2 -= &self.c1;
-        s2 += &self.c2;
-        s2.square_assign();
-        // bc
-        let mut bc = self.c1;
-        bc *= &self.c2;
-        // s3 = 2bc
-        let mut s3 = bc;
-        s3.double_assign();
-        // s4 = c^2
-        let mut s4 = self.c2;
-        s4.square_assign();
-
-        // new c0 = 2bc.mul_by_xi + a^2
-        self.c0 = s3;
-        self.c0.mul_by_nonresidue();
-        // self.c0.mul_by_xi();
-        self.c0 += &s0;
+impl CubicExtFieldArith for Fq6 {
+    type Base = Fq2;
+}
 
-        // new c1 = (c^2).mul_by_xi + 2ab
-        self.c1 = s4;
-        self.c1.mul_by_nonresidue();
-        self.c1 += &s1;
+impl CubicSparseMul for Fq6 {
+    type Base = Fq2;
+}
 
-        // new c2 = 2ab + (a - b + c)^2 + 2bc - a^2 - c^2 = b^2 + 2ac
-        self.c2 = s1;
-        self.c2 += &s2;
-        self.c2 += &s3;
-        self.c2 -= &s0;
-        self.c2 -= &s4;
-    }
+impl ExtField for Fq6 {
+    const NON_RESIDUE: Self = Fq6::new(Fq2::ZERO, Fq2::ONE, Fq2::ZERO);
 
-    pub fn frobenius_map(&mut self, power: usize) {
+    fn frobenius_map(&mut self, power: usize) {
         self.c0.frobenius_map(power);
         self.c1.frobenius_map(power);
         self.c2.frobenius_map(power);
-
         self.c1.mul_assign(&FROBENIUS_COEFF_FQ6_C1[power % 6]);
         self.c2.mul_assign(&FROBENIUS_COEFF_FQ6_C2[power % 6]);
     }
 
-    /// Multiply by cubic nonresidue v.
-    pub fn mul_by_nonresidue(&mut self) {
-        use std::mem::swap;
-        swap(&mut self.c0, &mut self.c1);
-        swap(&mut self.c0, &mut self.c2);
-        // c0, c1, c2 -> c2, c0, c1
-        self.c0.mul_by_nonresidue();
-    }
-
-    pub fn mul_by_1(&mut self, c1: &Fq2) {
-        let mut b_b = self.c1;
-        b_b *= c1;
-
-        let mut t1 = *c1;
-        {
-            let mut tmp = self.c1;
-            tmp += &self.c2;
-
-            t1 *= &tmp;
-            t1 -= &b_b;
-            t1.mul_by_nonresidue();
-        }
-
-        let mut t2 = *c1;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c1;
-
-            t2 *= &tmp;
-            t2 -= &b_b;
-        }
-
-        self.c0 = t1;
-        self.c1 = t2;
-        self.c2 = b_b;
-    }
-
-    pub fn mul_by_01(&mut self, c0: &Fq2, c1: &Fq2) {
-        let mut a_a = self.c0;
-        let mut b_b = self.c1;
-        a_a *= c0;
-        b_b *= c1;
-
-        let mut t1 = *c1;
-        {
-            let mut tmp = self.c1;
-            tmp += &self.c2;
-
-            t1 *= &tmp;
-            t1 -= &b_b;
-            t1.mul_by_nonresidue();
-            t1 += &a_a;
-        }
-
-        let mut t3 = *c0;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c2;
-
-            t3 *= &tmp;
-            t3 -= &a_a;
-            t3 += &b_b;
-        }
-
-        let mut t2 = *c0;
-        t2 += c1;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c1;
-
-            t2 *= &tmp;
-            t2 -= &a_a;
-            t2 -= &b_b;
-        }
-
-        self.c0 = t1;
-        self.c1 = t2;
-        self.c2 = t3;
-    }
-
-    pub fn invert(&self) -> CtOption<Self> {
-        let mut c0 = self.c2;
-        c0.mul_by_nonresidue();
-        c0 *= &self.c1;
-        c0 = -c0;
-        {
-            let mut c0s = self.c0;
-            c0s.square_assign();
-            c0 += &c0s;
-        }
-        let mut c1 = self.c2;
-        c1.square_assign();
-        c1.mul_by_nonresidue();
-        {
-            let mut c01 = self.c0;
-            c01 *= &self.c1;
-            c1 -= &c01;
-        }
-        let mut c2 = self.c1;
-        c2.square_assign();
-        {
-            let mut c02 = self.c0;
-            c02 *= &self.c2;
-            c2 -= &c02;
-        }
-
-        let mut tmp1 = self.c2;
-        tmp1 *= &c1;
-        let mut tmp2 = self.c1;
-        tmp2 *= &c2;
-        tmp1 += &tmp2;
-        tmp1.mul_by_nonresidue();
-        tmp2 = self.c0;
-        tmp2 *= &c0;
-        tmp1 += &tmp2;
-
-        tmp1.invert().map(|t| {
-            let mut tmp = Fq6 {
-                c0: t,
-                c1: t,
-                c2: t,
-            };
-            tmp.c0 *= &c0;
-            tmp.c1 *= &c1;
-            tmp.c2 *= &c2;
-
-            tmp
-        })
+    fn mul_by_nonresidue(self: &Fq6) -> Fq6 {
+        let c0 = self.c2.mul_by_nonresidue();
+        let c1 = self.c0;
+        let c2 = self.c1;
+        Self { c0, c1, c2 }
     }
 }
 
@@ -418,7 +206,7 @@ mod test {
     use super::*;
     crate::field_testing_suite!(Fq6, "field_arithmetic");
     // extension field-specific
-    crate::field_testing_suite!(Fq6, "f6_tests", Fq2);
+    crate::field_testing_suite!(Fq6, "cubic_sparse_mul", Fq2);
     crate::field_testing_suite!(
         Fq6,
         "frobenius",
@@ -426,11 +214,15 @@ mod test {
         //  ϕ: E → E
         //  (x, y) ↦ (x^p, y^p)
         // p: modulus of base field (Here, Fq::MODULUS)
-        [
-            0x3c208c16d87cfd47,
-            0x97816a916871ca8d,
-            0xb85045b68181585d,
-            0x30644e72e131a029,
-        ]
+        Fq::MODULUS_LIMBS
     );
+
+    #[test]
+    fn test_fq6_mul_nonresidue() {
+        use ff::Field;
+        let e = Fq6::random(rand_core::OsRng);
+        let a0 = e.mul_by_nonresidue();
+        let a1 = e * Fq6::NON_RESIDUE;
+        assert_eq!(a0, a1);
+    }
 }
diff --git a/src/bn256/fr.rs b/src/bn256/fr.rs
index e4ff25fd..ba82383c 100644
--- a/src/bn256/fr.rs
+++ b/src/bn256/fr.rs
@@ -3,13 +3,6 @@ use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     bn256_scalar,
     Fr,
@@ -20,12 +13,12 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fr);
-impl_binops_calls!(Fr);
-impl_binops_additive!(Fr, Fr);
-impl_binops_multiplicative!(Fr, Fr);
-field_bits!(Fr);
-serialize_deserialize_primefield!(Fr);
+crate::extend_field_legendre!(Fr);
+crate::impl_binops_calls!(Fr);
+crate::impl_binops_additive!(Fr, Fr);
+crate::impl_binops_multiplicative!(Fr, Fr);
+crate::field_bits!(Fr);
+crate::serialize_deserialize_primefield!(Fr);
 
 #[cfg(feature = "bn256-table")]
 pub use table::FR_TABLE;
diff --git a/src/bn256/mod.rs b/src/bn256/mod.rs
index 45c9f942..4cbbea03 100644
--- a/src/bn256/mod.rs
+++ b/src/bn256/mod.rs
@@ -13,3 +13,12 @@ pub use fq12::*;
 pub use fq2::*;
 pub use fq6::*;
 pub use fr::*;
+
+pub const BN_X: u64 = 4965661367192848881;
+
+// 6U+2 for in NAF form
+pub const SIX_U_PLUS_2_NAF: [i8; 65] = [
+    0, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 0, -1, 0, 0, 1, 0, -1, 0, 0, 0, 0,
+    1, 1, 1, 0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 1, 0, 0, -1, 0, 0, 0, 1, 1, 0, -1, 0,
+    0, 1, 0, 1, 1,
+];
diff --git a/src/derive/field/tower.rs b/src/derive/field/tower.rs
index 00dd101b..f20c15c5 100644
--- a/src/derive/field/tower.rs
+++ b/src/derive/field/tower.rs
@@ -1,139 +1,16 @@
 #[macro_export]
-macro_rules! impl_tower2_common {
-    (
-        $field:ident,
-        $tower:ident,
-        serde
-    ) => {
-        #[derive(Copy, Clone, Debug, Eq, PartialEq, Default)]
-        #[cfg_attr(feature = "derive_serde", derive(Serialize, Deserialize))]
-        pub struct $tower {
-            pub c0: $field,
-            pub c1: $field,
-        }
-
-        impl_tower2_common!($field, $tower, implementations);
-    };
-
-    (
-        $field:ident,
-        $tower:ident
-
-    ) => {
-        #[derive(Copy, Clone, Debug, Eq, PartialEq, Default)]
-        pub struct $tower {
-            pub c0: $field,
-            pub c1: $field,
-        }
-
-        impl_tower2_common!($field, $tower, implementations);
-    };
-
+macro_rules! impl_tower2 {
     (
-        $field:ident,
-        $tower:ident,
-        implementations
+        $base:ident,
+        $field:ident
     ) => {
-        impl ::core::ops::Neg for $tower {
-            type Output = $tower;
-
-            #[inline]
-            fn neg(self) -> $tower {
-                -&self
-            }
-        }
-
-        impl ConditionallySelectable for $tower {
-            fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self {
-                $tower {
-                    c0: $field::conditional_select(&a.c0, &b.c0, choice),
-                    c1: $field::conditional_select(&a.c1, &b.c1, choice),
-                }
-            }
-        }
-
-        impl ConstantTimeEq for $tower {
-            fn ct_eq(&self, other: &Self) -> Choice {
-                self.c0.ct_eq(&other.c0) & self.c1.ct_eq(&other.c1)
-            }
-        }
-
-        impl $tower {
-            #[inline]
-            pub const fn zero() -> $tower {
-                $tower {
-                    c0: $field::zero(),
-                    c1: $field::zero(),
-                }
-            }
-
-            #[inline]
-            pub const fn one() -> $tower {
-                $tower {
-                    c0: $field::one(),
-                    c1: $field::zero(),
-                }
-            }
-
-            #[inline]
-            pub const fn new(c0: $field, c1: $field) -> Self {
-                $tower { c0, c1 }
-            }
-
-            pub fn double(&self) -> Self {
-                Self {
-                    c0: self.c0.double(),
-                    c1: self.c1.double(),
-                }
-            }
-
-            pub fn add(&self, other: &Self) -> Self {
-                Self {
-                    c0: self.c0.add(&other.c0),
-                    c1: self.c1.add(&other.c1),
-                }
-            }
-
-            pub fn sub(&self, other: &Self) -> Self {
-                Self {
-                    c0: self.c0.sub(&other.c0),
-                    c1: self.c1.sub(&other.c1),
-                }
-            }
-
-            pub fn neg(&self) -> Self {
-                Self {
-                    c0: self.c0.neg(),
-                    c1: self.c1.neg(),
-                }
-            }
-
-            pub fn mul(&self, other: &Self) -> Self {
-                let mut t = *other;
-                t.mul_assign(self);
-                t
-            }
-
-            pub fn square(&self) -> Self {
-                let mut t = *self;
-                t.square_assign();
-                t
-            }
+        impl $field {
+            pub const SIZE: usize = $base::SIZE * 2;
         }
-    };
-}
 
-#[macro_export]
-macro_rules! impl_tower2 {
-    (
-        $field:ident,
-        $tower:ident,
-        $repr:ident
-    ) => {
-        /// `$tower` elements are ordered lexicographically.
-        impl Ord for $tower {
+        impl Ord for $field {
             #[inline(always)]
-            fn cmp(&self, other: &$tower) -> Ordering {
+            fn cmp(&self, other: &$field) -> Ordering {
                 match self.c1.cmp(&other.c1) {
                     Ordering::Greater => Ordering::Greater,
                     Ordering::Less => Ordering::Less,
@@ -142,36 +19,30 @@ macro_rules! impl_tower2 {
             }
         }
 
-        impl PartialOrd for $tower {
+        impl PartialOrd for $field {
             #[inline(always)]
-            fn partial_cmp(&self, other: &$tower) -> Option<Ordering> {
+            fn partial_cmp(&self, other: &$field) -> Option<Ordering> {
                 Some(self.cmp(other))
             }
         }
 
-        impl From<$tower> for [u8; $tower::SIZE] {
-            fn from(value: $tower) -> [u8; $tower::SIZE] {
-                value.to_bytes()
-            }
-        }
-
-        impl<'a> From<&'a $tower> for [u8; $tower::SIZE] {
-            fn from(value: &'a $tower) -> [u8; $tower::SIZE] {
-                value.to_bytes()
+        impl From<u64> for $field {
+            fn from(val: u64) -> Self {
+                $field {
+                    c0: $base::from(val),
+                    c1: $base::ZERO,
+                }
             }
         }
 
-        impl $tower {
-            pub const SIZE: usize = $field::SIZE * 2;
-
+        impl $field {
             /// Attempts to convert a little-endian byte representation of
-            /// a scalar into a `$field`, failing if the input is not canonical.
-            pub fn from_bytes(bytes: &[u8; $field::SIZE * 2]) -> CtOption<$tower> {
-                let c0 = $field::from_bytes(bytes[0..$field::SIZE].try_into().unwrap());
-                let c1 =
-                    $field::from_bytes(bytes[$field::SIZE..$field::SIZE * 2].try_into().unwrap());
+            /// a scalar into a `$base`, failing if the input is not canonical.
+            pub fn from_bytes(bytes: &[u8; $base::SIZE * 2]) -> CtOption<$field> {
+                let c0 = $base::from_bytes(bytes[0..$base::SIZE].try_into().unwrap());
+                let c1 = $base::from_bytes(bytes[$base::SIZE..$base::SIZE * 2].try_into().unwrap());
                 CtOption::new(
-                    $tower {
+                    $field {
                         c0: c0.unwrap(),
                         c1: c1.unwrap(),
                     },
@@ -179,85 +50,73 @@ macro_rules! impl_tower2 {
                 )
             }
 
-            /// Converts an element of `$field` into a byte representation in
+            /// Converts an element of `$base` into a byte representation in
             /// little-endian byte order.
             #[allow(clippy::wrong_self_convention)]
-            pub fn to_bytes(&self) -> [u8; $field::SIZE * 2] {
-                let mut res = [0u8; $field::SIZE * 2];
+            pub fn to_bytes(&self) -> [u8; $base::SIZE * 2] {
+                let mut res = [0u8; $base::SIZE * 2];
                 let c0_bytes = self.c0.to_bytes();
                 let c1_bytes = self.c1.to_bytes();
-                res[0..$field::SIZE].copy_from_slice(&c0_bytes[..]);
-                res[$field::SIZE..$field::SIZE * 2].copy_from_slice(&c1_bytes[..]);
+                res[0..$base::SIZE].copy_from_slice(&c0_bytes[..]);
+                res[$base::SIZE..$base::SIZE * 2].copy_from_slice(&c1_bytes[..]);
                 res
             }
         }
 
-        impl ff::Field for $tower {
-            const ZERO: Self = Self::zero();
-            const ONE: Self = Self::one();
-
-            fn random(mut rng: impl RngCore) -> Self {
-                $tower {
-                    c0: $field::random(&mut rng),
-                    c1: $field::random(&mut rng),
-                }
-            }
-
-            fn is_zero(&self) -> Choice {
-                self.c0.is_zero() & self.c1.is_zero()
-            }
-
-            fn square(&self) -> Self {
-                self.square()
-            }
-
-            fn double(&self) -> Self {
-                self.double()
-            }
-
-            fn sqrt(&self) -> CtOption<Self> {
-                self.sqrt()
-            }
-
-            fn sqrt_ratio(num: &Self, div: &Self) -> (Choice, Self) {
-                ff::helpers::sqrt_ratio_generic(num, div)
-            }
+        impl WithSmallOrderMulGroup<3> for $field {
+            const ZETA: Self = $field {
+                c0: $base::ZETA.mul_const(&$base::ZETA),
+                c1: $base::ZERO,
+            };
+        }
 
-            fn invert(&self) -> CtOption<Self> {
-                self.invert()
+        impl Legendre for $field {
+            fn legendre(&self) -> i64 {
+                self.norm().legendre()
             }
         }
 
-        impl PrimeField for $tower {
-            type Repr = $repr;
+        impl PrimeField for $field {
+            type Repr = $crate::serde::Repr<{ $base::SIZE * 2 }>;
 
-            const MODULUS: &'static str = <$field as PrimeField>::MODULUS;
-            const MULTIPLICATIVE_GENERATOR: Self = $tower {
-                c0: $field::MULTIPLICATIVE_GENERATOR,
-                c1: $field::ZERO,
+            const MODULUS: &'static str = <$base as PrimeField>::MODULUS;
+            const MULTIPLICATIVE_GENERATOR: Self = $field {
+                c0: $base::MULTIPLICATIVE_GENERATOR,
+                c1: $base::ZERO,
             };
-            const NUM_BITS: u32 = $field::NUM_BITS;
-            const CAPACITY: u32 = $field::NUM_BITS;
-            const S: u32 = $field::S;
+            const NUM_BITS: u32 = $base::NUM_BITS;
+            const CAPACITY: u32 = $base::NUM_BITS;
+            const S: u32 = $base::S;
 
             // TODO: Check that we can just 0 this and forget.
-            const ROOT_OF_UNITY: Self = $tower::ZERO;
-            const ROOT_OF_UNITY_INV: Self = $tower::ZERO;
-            const DELTA: Self = $tower::ZERO;
+            const ROOT_OF_UNITY: Self = $field::ZERO;
+            const ROOT_OF_UNITY_INV: Self = $field::ZERO;
+            const DELTA: Self = $field::ZERO;
 
-            const TWO_INV: Self = $tower {
-                c0: $field::TWO_INV,
-                c1: $field::ZERO,
+            const TWO_INV: Self = $field {
+                c0: $base::TWO_INV,
+                c1: $base::ZERO,
             };
 
             fn from_repr(repr: Self::Repr) -> CtOption<Self> {
-                let c0 = $field::from_bytes(&repr.0[..$field::SIZE].try_into().unwrap());
-                let c1 = $field::from_bytes(&repr.0[$field::SIZE..].try_into().unwrap());
-                CtOption::new($tower::new(c0.unwrap(), c1.unwrap()), Choice::from(1))
+                let c0: [u8; $base::SIZE] = repr[..$base::SIZE].try_into().unwrap();
+                let c0: <$base as PrimeField>::Repr = c0.into();
+                let c0 = $base::from_repr(c0);
+
+                let c1: [u8; $base::SIZE] = repr[$base::SIZE..].try_into().unwrap();
+                let c1: <$base as PrimeField>::Repr = c1.into();
+                let c1 = $base::from_repr(c1);
+
+                CtOption::new($field::new(c0.unwrap(), c1.unwrap()), Choice::from(1))
             }
 
             fn to_repr(&self) -> Self::Repr {
-                $repr(self.to_bytes())
+                let mut res = Self::Repr::default();
+                let c0 = self.c0.to_repr();
+                let c1 = self.c1.to_repr();
+                res[0..$base::SIZE].copy_from_slice(&c0.as_ref()[..]);
+                res[$base::SIZE..$base::SIZE * 2].copy_from_slice(&c1.as_ref()[..]);
+                res
             }
 
             fn is_odd(&self) -> Choice {
@@ -265,79 +124,35 @@ macro_rules! impl_tower2 {
             }
         }
 
-        impl WithSmallOrderMulGroup<3> for $tower {
-            // $field::ZETA ^2
-            const ZETA: Self = $tower {
-                c0: $field::ZETA.mul_const(&$field::ZETA),
-                c1: $field::ZERO,
-            };
-        }
-
-        impl From<u64> for $tower {
-            fn from(val: u64) -> Self {
-                $tower {
-                    c0: $field::from(val),
-                    c1: $field::ZERO,
-                }
-            }
-        }
-
-        impl Legendre for $tower {
-            fn legendre(&self) -> i64 {
-                self.norm().legendre()
-            }
-        }
-
-        #[derive(Clone, Copy, Debug)]
-        pub struct $repr([u8; $field::SIZE * 2]);
-
-        impl Default for $repr {
-            fn default() -> Self {
-                Self([0u8; $field::SIZE * 2])
-            }
-        }
-
-        impl AsMut<[u8]> for $repr {
-            fn as_mut(&mut self) -> &mut [u8] {
-                &mut self.0
-            }
-        }
-
-        impl AsRef<[u8]> for $repr {
-            fn as_ref(&self) -> &[u8] {
-                &self.0
-            }
-        }
-
-        impl $crate::serde::SerdeObject for $tower {
+        impl $crate::serde::SerdeObject for $field {
             fn from_raw_bytes_unchecked(bytes: &[u8]) -> Self {
-                debug_assert_eq!(bytes.len(), $field::SIZE * 2);
-                let [c0, c1] = [0, $field::SIZE]
-                    .map(|i| $field::from_raw_bytes_unchecked(&bytes[i..i + $field::SIZE]));
+                debug_assert_eq!(bytes.len(), $base::SIZE * 2);
+                let [c0, c1] = [0, $base::SIZE]
+                    .map(|i| $base::from_raw_bytes_unchecked(&bytes[i..i + $base::SIZE]));
                 Self { c0, c1 }
             }
             fn from_raw_bytes(bytes: &[u8]) -> Option<Self> {
-                if bytes.len() != $field::SIZE * 2 {
+                if bytes.len() != $base::SIZE * 2 {
                     return None;
                 }
                 let [c0, c1] =
-                    [0, $field::SIZE].map(|i| $field::from_raw_bytes(&bytes[i..i + $field::SIZE]));
+                    [0, $base::SIZE].map(|i| $base::from_raw_bytes(&bytes[i..i + $base::SIZE]));
                 c0.zip(c1).map(|(c0, c1)| Self { c0, c1 })
             }
             fn to_raw_bytes(&self) -> Vec<u8> {
-                let mut res = Vec::with_capacity($field::SIZE * 2);
+                let mut res = Vec::with_capacity($base::SIZE * 2);
                 for limb in self.c0.0.iter().chain(self.c1.0.iter()) {
                     res.extend_from_slice(&limb.to_le_bytes());
                 }
                 res
             }
             fn read_raw_unchecked<R: std::io::Read>(reader: &mut R) -> Self {
-                let [c0, c1] = [(); 2].map(|_| $field::read_raw_unchecked(reader));
+                let [c0, c1] = [(); 2].map(|_| $base::read_raw_unchecked(reader));
                 Self { c0, c1 }
             }
             fn read_raw<R: std::io::Read>(reader: &mut R) -> std::io::Result<Self> {
-                let c0 = $field::read_raw(reader)?;
-                let c1 = $field::read_raw(reader)?;
+                let c0 = $base::read_raw(reader)?;
+                let c1 = $base::read_raw(reader)?;
                 Ok(Self { c0, c1 })
             }
             fn write_raw<W: std::io::Write>(&self, writer: &mut W) -> std::io::Result<()> {
@@ -349,151 +164,61 @@ macro_rules! impl_tower2 {
 }
 
 #[macro_export]
-macro_rules! impl_tower6 {
+macro_rules! impl_tower2_from_uniform_bytes {
     (
+        $base:ident,
         $field:ident,
-        $tower2:ident,
-        $tower6:ident
+        $size:expr
     ) => {
-        #[derive(Copy, Clone, Debug, Eq, PartialEq, Default)]
-        pub struct $tower6 {
-            pub c0: $tower2,
-            pub c1: $tower2,
-            pub c2: $tower2,
-        }
-
-        impl ::core::ops::Neg for $tower6 {
-            type Output = $tower6;
-
-            #[inline]
-            fn neg(self) -> $tower6 {
-                -&self
-            }
-        }
-
-        impl subtle::ConditionallySelectable for $tower6 {
-            fn conditional_select(a: &Self, b: &Self, choice: subtle::Choice) -> Self {
-                $tower6 {
-                    c0: $tower2::conditional_select(&a.c0, &b.c0, choice),
-                    c1: $tower2::conditional_select(&a.c1, &b.c1, choice),
-                    c2: $tower2::conditional_select(&a.c2, &b.c2, choice),
-                }
-            }
-        }
-
-        impl subtle::ConstantTimeEq for $tower6 {
-            fn ct_eq(&self, other: &Self) -> subtle::Choice {
-                self.c0.ct_eq(&other.c0) & self.c1.ct_eq(&other.c1) & self.c2.ct_eq(&other.c2)
-            }
-        }
-
-        impl $tower6 {
-            #[inline]
-            pub const fn zero() -> Self {
-                $tower6 {
-                    c0: $tower2::ZERO,
-                    c1: $tower2::ZERO,
-                    c2: $tower2::ZERO,
-                }
-            }
-
-            #[inline]
-            pub const fn one() -> Self {
-                $tower6 {
-                    c0: $tower2::ONE,
-                    c1: $tower2::ZERO,
-                    c2: $tower2::ZERO,
-                }
-            }
-
-            #[inline]
-            pub fn double(&self) -> Self {
-                Self {
-                    c0: self.c0.double(),
-                    c1: self.c1.double(),
-                    c2: self.c2.double(),
-                }
-            }
-
-            #[inline]
-            pub fn add(&self, other: &Self) -> Self {
-                Self {
-                    c0: self.c0 + other.c0,
-                    c1: self.c1 + other.c1,
-                    c2: self.c2 + other.c2,
-                }
-            }
-
-            #[inline]
-            pub fn sub(&self, other: &Self) -> Self {
-                Self {
-                    c0: self.c0 - other.c0,
-                    c1: self.c1 - other.c1,
-                    c2: self.c2 - other.c2,
-                }
-            }
-
-            #[inline]
-            pub fn neg(&self) -> Self {
-                Self {
-                    c0: -self.c0,
-                    c1: -self.c1,
-                    c2: -self.c2,
-                }
-            }
-
-            pub fn mul(&self, other: &Self) -> Self {
-                let mut t = *other;
-                t.mul_assign(self);
-                t
-            }
-
-            pub fn square(&self) -> Self {
-                let mut t = *self;
-                t.square_assign();
-                t
-            }
-
-            pub fn random(mut rng: impl RngCore) -> Self {
-                $tower6 {
-                    c0: $tower2::random(&mut rng),
-                    c1: $tower2::random(&mut rng),
-                    c2: $tower2::random(&mut rng),
-                }
+        impl FromUniformBytes<{ $size }> for $field {
+            fn from_uniform_bytes(bytes: &[u8; $size]) -> Self {
+                assert!($size % 2 == 0);
+                const SIZE: usize = $size / 2;
+                let c0: [u8; SIZE] = bytes[SIZE..].try_into().unwrap();
+                let c1: [u8; SIZE] = bytes[..SIZE].try_into().unwrap();
+                Self::new(
+                    $base::from_uniform_bytes(&c0),
+                    $base::from_uniform_bytes(&c1),
+                )
             }
         }
+    };
+}
 
-        #[cfg(test)]
-        impl Field for $tower6 {
-            const ZERO: Self = Self::zero();
-            const ONE: Self = Self::one();
-
-            fn random(rng: impl RngCore) -> Self {
-                $tower6::random(rng)
-            }
-
-            fn is_zero(&self) -> subtle::Choice {
-                self.c0.is_zero() & self.c1.is_zero()
-            }
-
-            fn square(&self) -> Self {
-                self.square()
-            }
-
-            fn double(&self) -> Self {
-                self.double()
-            }
-
-            fn sqrt(&self) -> CtOption<Self> {
-                unimplemented!()
-            }
-
-            fn sqrt_ratio(_num: &Self, _div: &Self) -> (subtle::Choice, Self) {
-                unimplemented!()
-            }
-
-            fn invert(&self) -> CtOption<Self> {
-                self.invert()
+#[macro_export]
+macro_rules! impl_cyclotomic_square {
+    (
+        $tower2:ident,
+        $tower12:ident
+    ) => {
+        impl $tower12 {
+            pub fn cyclotomic_square(&mut self) {
+                fn fp4_square(c0: &mut $tower2, c1: &mut $tower2, a0: &$tower2, a1: &$tower2) {
+                    use ff::Field;
+                    let t0 = a0.square();
+                    let t1 = a1.square();
+                    *c0 = t1.mul_by_nonresidue() + t0;
+                    *c1 = (a0 + a1).square() - t0 - t1;
+                }
+
+                let mut t3 = $tower2::zero();
+                let mut t4 = $tower2::zero();
+                let mut t5 = $tower2::zero();
+                let mut t6 = $tower2::zero();
+                fp4_square(&mut t3, &mut t4, &self.c0.c0, &self.c1.c1);
+
+                self.c0.c0 = (t3 - self.c0.c0).double() + t3;
+                self.c1.c1 = (t4 + self.c1.c1).double() + t4;
+
+                fp4_square(&mut t3, &mut t4, &self.c1.c0, &self.c0.c2);
+                fp4_square(&mut t5, &mut t6, &self.c0.c1, &self.c1.c2);
+
+                self.c0.c1 = (t3 - self.c0.c1).double() + t3;
+                self.c1.c2 = (t4 + self.c1.c2).double() + t4;
+
+                let t3 = t6.mul_by_nonresidue();
+                self.c1.c0 = (t3 + self.c1.c0).double() + t3;
+                self.c0.c2 = (t5 - self.c0.c2).double() + t5;
             }
         }
     };
diff --git a/src/derive/mod.rs b/src/derive/mod.rs
index ad560f3f..29535218 100644
--- a/src/derive/mod.rs
+++ b/src/derive/mod.rs
@@ -2,10 +2,21 @@
 pub mod curve;
 #[macro_use]
 pub mod field;
+#[macro_use]
+pub mod pairing;
 
 #[macro_export]
 macro_rules! impl_binops_calls {
     ($field:ident) => {
+        impl ::core::ops::Neg for $field {
+            type Output = $field;
+
+            #[inline]
+            fn neg(self) -> $field {
+                -&self
+            }
+        }
+
         impl<'a> ::core::ops::Neg for &'a $field {
             type Output = $field;
 
@@ -111,8 +122,8 @@ macro_rules! impl_sub_binop_specify_output {
 #[macro_export]
 macro_rules! impl_binops_additive_specify_output {
     ($lhs:ident, $rhs:ident, $output:ident) => {
-        impl_add_binop_specify_output!($lhs, $rhs, $output);
-        impl_sub_binop_specify_output!($lhs, $rhs, $output);
+        $crate::impl_add_binop_specify_output!($lhs, $rhs, $output);
+        $crate::impl_sub_binop_specify_output!($lhs, $rhs, $output);
     };
 }
 
@@ -151,7 +162,7 @@ macro_rules! impl_binops_multiplicative_mixed {
 #[macro_export]
 macro_rules! impl_binops_additive {
     ($lhs:ident, $rhs:ident) => {
-        impl_binops_additive_specify_output!($lhs, $rhs, $lhs);
+        $crate::impl_binops_additive_specify_output!($lhs, $rhs, $lhs);
 
         impl ::core::ops::SubAssign<$rhs> for $lhs {
             #[inline]
@@ -186,7 +197,7 @@ macro_rules! impl_binops_additive {
 #[macro_export]
 macro_rules! impl_binops_multiplicative {
     ($lhs:ident, $rhs:ident) => {
-        impl_binops_multiplicative_mixed!($lhs, $rhs, $lhs);
+        $crate::impl_binops_multiplicative_mixed!($lhs, $rhs, $lhs);
 
         impl ::core::ops::MulAssign<$rhs> for $lhs {
             #[inline]
diff --git a/src/derive/pairing.rs b/src/derive/pairing.rs
new file mode 100644
index 00000000..21cf437a
--- /dev/null
+++ b/src/derive/pairing.rs
@@ -0,0 +1,268 @@
+#[macro_export]
+macro_rules! impl_miller_loop_components {
+    (
+    $engine:ident,
+    $g1:ident,
+    $g1affine:ident,
+    $g2:ident,
+    $g2affine:ident,
+    $base:ident,
+    $target:ident,
+    $scalar:ident
+    ) => {
+        #[derive(Clone, Debug)]
+        pub struct $engine;
+
+        impl Engine for $engine {
+            type Fr = $scalar;
+            type G1 = $g1;
+            type G1Affine = $g1affine;
+            type G2 = $g2;
+            type G2Affine = $g2affine;
+            type Gt = Gt;
+
+            fn pairing(p: &Self::G1Affine, q: &Self::G2Affine) -> Self::Gt {
+                $engine::multi_miller_loop(&[(p, q)]).final_exponentiation()
+            }
+        }
+
+        impl MultiMillerLoop for $engine {
+            type G2Prepared = $g2affine;
+            type Result = $base;
+
+            fn multi_miller_loop(terms: &[(&Self::G1Affine, &Self::G2Prepared)]) -> Self::Result {
+                multi_miller_loop(terms)
+            }
+        }
+
+        impl PairingCurveAffine for $g1affine {
+            type Pair = $g2affine;
+            type PairingResult = $target;
+
+            fn pairing_with(&self, other: &Self::Pair) -> Self::PairingResult {
+                $engine::pairing(&self, &other)
+            }
+        }
+
+        impl PairingCurveAffine for $g2affine {
+            type Pair = $g1affine;
+            type PairingResult = $target;
+
+            fn pairing_with(&self, other: &Self::Pair) -> Self::PairingResult {
+                $engine::pairing(&other, &self)
+            }
+        }
+
+        fn double(f: &mut $base, r: &mut $g2, p: &$g1affine) {
+            use ff::Field;
+            let t0 = r.x.square();
+            let t1 = r.y.square();
+            let t2 = t1.square();
+            let t3 = (t1 + r.x).square() - t0 - t2;
+            let t3 = t3 + t3;
+            let t4 = t0 + t0 + t0;
+            let t6 = r.x + t4;
+            let t5 = t4.square();
+            let zsquared = r.z.square();
+            r.x = t5 - t3 - t3;
+            r.z = (r.z + r.y).square() - t1 - zsquared;
+            r.y = (t3 - r.x) * t4;
+            let t2 = t2 + t2;
+            let t2 = t2 + t2;
+            let t2 = t2 + t2;
+            r.y -= t2;
+            let t3 = t4 * zsquared;
+            let t3 = t3 + t3;
+            let t3 = -t3;
+            let t6 = t6.square() - t0 - t5;
+            let t1 = t1 + t1;
+            let t1 = t1 + t1;
+            let t6 = t6 - t1;
+            let t0 = r.z * zsquared;
+            let t0 = t0 + t0;
+
+            ell(f, &(t0, t3, t6), p);
+        }
+
+        fn add(f: &mut $base, r: &mut $g2, q: &$g2affine, p: &$g1affine) {
+            use ff::Field;
+            let zsquared = r.z.square();
+            let ysquared = q.y.square();
+            let t0 = zsquared * q.x;
+            let t1 = ((q.y + r.z).square() - ysquared - zsquared) * zsquared;
+            let t2 = t0 - r.x;
+            let t3 = t2.square();
+            let t4 = t3 + t3;
+            let t4 = t4 + t4;
+            let t5 = t4 * t2;
+            let t6 = t1 - r.y - r.y;
+            let t9 = t6 * q.x;
+            let t7 = t4 * r.x;
+            r.x = t6.square() - t5 - t7 - t7;
+            r.z = (r.z + t2).square() - zsquared - t3;
+            let t10 = q.y + r.z;
+            let t8 = (t7 - r.x) * t6;
+            let t0 = r.y * t5;
+            let t0 = t0 + t0;
+            r.y = t8 - t0;
+            let t10 = t10.square() - ysquared;
+            let ztsquared = r.z.square();
+            let t10 = t10 - ztsquared;
+            let t9 = t9 + t9 - t10;
+            let t10 = r.z + r.z;
+            let t6 = -t6;
+            let t1 = t6 + t6;
+
+            ell(f, &(t10, t1, t9), p);
+        }
+    };
+}
+
+#[macro_export]
+macro_rules! impl_gt {
+    (
+        $target:ident,
+        $base:ident,
+        $scalar:ident
+    ) => {
+        #[derive(Copy, Clone, Debug, Default)]
+        pub struct $target(pub(crate) $base);
+
+        impl ConstantTimeEq for $target {
+            fn ct_eq(&self, other: &Self) -> Choice {
+                self.0.ct_eq(&other.0)
+            }
+        }
+
+        impl ConditionallySelectable for $target {
+            fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self {
+                $target($base::conditional_select(&a.0, &b.0, choice))
+            }
+        }
+
+        impl Eq for $target {}
+        impl PartialEq for $target {
+            #[inline]
+            fn eq(&self, other: &Self) -> bool {
+                bool::from(self.ct_eq(other))
+            }
+        }
+
+        impl $target {
+            /// Returns the group identity, which is $1$.
+            pub fn identity() -> $target {
+                $target($base::one())
+            }
+
+            /// Doubles this group element.
+            pub fn double(&self) -> $target {
+                use ff::Field;
+                $target(self.0.square())
+            }
+        }
+
+        impl<'a> Neg for &'a $target {
+            type Output = $target;
+
+            #[inline]
+            fn neg(self) -> $target {
+                // The element is unitary, so we just conjugate.
+                let mut u = self.0;
+                u.conjugate();
+                $target(u)
+            }
+        }
+
+        impl Neg for $target {
+            type Output = $target;
+
+            #[inline]
+            fn neg(self) -> $target {
+                -&self
+            }
+        }
+
+        impl<'a, 'b> Add<&'b $target> for &'a $target {
+            type Output = $target;
+
+            #[inline]
+            #[allow(clippy::suspicious_arithmetic_impl)]
+            fn add(self, rhs: &'b $target) -> $target {
+                $target(self.0 * rhs.0)
+            }
+        }
+
+        impl<'a, 'b> Sub<&'b $target> for &'a $target {
+            type Output = $target;
+
+            #[inline]
+            fn sub(self, rhs: &'b $target) -> $target {
+                self + (-rhs)
+            }
+        }
+
+        #[allow(clippy::suspicious_arithmetic_impl)]
+        impl<'a, 'b> Mul<&'b $scalar> for &'a $target {
+            type Output = $target;
+
+            fn mul(self, other: &'b $scalar) -> Self::Output {
+                let mut acc = $target::identity();
+
+                for bit in other
+                    .to_repr()
+                    .as_ref()
+                    .iter()
+                    .rev()
+                    .flat_map(|byte| (0..8).rev().map(move |i| Choice::from((byte >> i) & 1u8)))
+                    .skip(1)
+                {
+                    acc = acc.double();
+                    acc = $target::conditional_select(&acc, &(acc + self), bit);
+                }
+
+                acc
+            }
+        }
+
+        $crate::impl_binops_additive!($target, $target);
+        $crate::impl_binops_multiplicative!($target, $scalar);
+
+        impl<T> Sum<T> for $target
+        where
+            T: Borrow<$target>,
+        {
+            fn sum<I>(iter: I) -> Self
+            where
+                I: Iterator<Item = T>,
+            {
+                iter.fold(Self::identity(), |acc, item| acc + item.borrow())
+            }
+        }
+
+        impl Group for $target {
+            type Scalar = $scalar;
+
+            fn random(rng: impl RngCore) -> Self {
+                use ff::Field;
+                $base::random(rng).final_exponentiation()
+            }
+
+            fn identity() -> Self {
+                Self::identity()
+            }
+
+            fn generator() -> Self {
+                unimplemented!()
+            }
+
+            fn is_identity(&self) -> Choice {
+                self.ct_eq(&Self::identity())
+            }
+
+            #[must_use]
+            fn double(&self) -> Self {
+                self.double()
+            }
+        }
+    };
+}
diff --git a/src/ff_ext/cubic.rs b/src/ff_ext/cubic.rs
new file mode 100644
index 00000000..c8e90808
--- /dev/null
+++ b/src/ff_ext/cubic.rs
@@ -0,0 +1,236 @@
+use super::ExtField;
+
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
+pub struct CubicExtField<F: ff::Field> {
+    pub(crate) c0: F,
+    pub(crate) c1: F,
+    pub(crate) c2: F,
+}
+
+pub trait CubicSparseMul {
+    type Base: ExtField;
+
+    fn mul_by_1(lhs: &CubicExtField<Self::Base>, c1: &Self::Base) -> CubicExtField<Self::Base> {
+        let b_b = lhs.c1 * c1;
+
+        let t1 = (lhs.c1 + lhs.c2) * c1 - b_b;
+        let t1 = t1.mul_by_nonresidue();
+        let t2 = (lhs.c0 + lhs.c1) * c1 - b_b;
+
+        CubicExtField {
+            c0: t1,
+            c1: t2,
+            c2: b_b,
+        }
+    }
+
+    fn mul_by_01(
+        lhs: &CubicExtField<Self::Base>,
+        c0: &Self::Base,
+        c1: &Self::Base,
+    ) -> CubicExtField<Self::Base> {
+        let a_a = lhs.c0 * c0;
+        let b_b = lhs.c1 * c1;
+
+        let t1 = *c1 * (lhs.c1 + lhs.c2) - b_b;
+        let t1 = a_a + t1.mul_by_nonresidue();
+        let t3 = *c0 * (lhs.c0 + lhs.c2) - a_a + b_b;
+        let t2 = (*c0 + c1) * (lhs.c0 + lhs.c1) - a_a - b_b;
+
+        CubicExtField {
+            c0: t1,
+            c1: t2,
+            c2: t3,
+        }
+    }
+}
+
+pub trait CubicExtFieldArith {
+    type Base: ExtField;
+
+    fn mul_assign(lhs: &mut CubicExtField<Self::Base>, rhs: &CubicExtField<Self::Base>) {
+        let a_a = lhs.c0 * rhs.c0;
+        let b_b = lhs.c1 * rhs.c1;
+        let c_c = lhs.c2 * rhs.c2;
+
+        let t1 = (rhs.c1 + rhs.c2) * (lhs.c1 + lhs.c2) - (c_c + b_b);
+
+        let t1 = a_a + t1.mul_by_nonresidue();
+
+        let t3 = (rhs.c0 + rhs.c2) * (lhs.c0 + lhs.c2) - (a_a - b_b + c_c);
+
+        let t2 = (rhs.c0 + rhs.c1) * (lhs.c0 + lhs.c1) - (a_a + b_b);
+        let t2 = t2 + c_c.mul_by_nonresidue();
+
+        lhs.c0 = t1;
+        lhs.c1 = t2;
+        lhs.c2 = t3;
+    }
+
+    fn square_assign(el: &mut CubicExtField<Self::Base>) {
+        use ff::Field;
+
+        let s0 = el.c0.square();
+        let s1 = (el.c0 * el.c1).double();
+        let s2 = (el.c0 - el.c1 + el.c2).square();
+        let s3 = (el.c1 * el.c2).double();
+        let s4 = el.c2.square();
+
+        el.c0 = s3.mul_by_nonresidue() + s0;
+        el.c1 = s4.mul_by_nonresidue() + s1;
+        el.c2 = s1 + s2 + s3 - s0 - s4;
+    }
+}
+
+impl<F: ff::Field> CubicExtField<F> {
+    #[inline]
+    pub const fn new(c0: F, c1: F, c2: F) -> Self {
+        Self { c0, c1, c2 }
+    }
+
+    #[inline]
+    pub const fn zero() -> Self {
+        Self {
+            c0: F::ZERO,
+            c1: F::ZERO,
+            c2: F::ZERO,
+        }
+    }
+
+    #[inline]
+    pub const fn one() -> Self {
+        Self {
+            c0: F::ONE,
+            c1: F::ZERO,
+            c2: F::ZERO,
+        }
+    }
+
+    #[inline]
+    pub fn double(&self) -> Self {
+        Self {
+            c0: self.c0.double(),
+            c1: self.c1.double(),
+            c2: self.c2.double(),
+        }
+    }
+
+    #[inline]
+    pub fn add(&self, other: &Self) -> Self {
+        Self {
+            c0: self.c0 + other.c0,
+            c1: self.c1 + other.c1,
+            c2: self.c2 + other.c2,
+        }
+    }
+
+    #[inline]
+    pub fn sub(&self, other: &Self) -> Self {
+        Self {
+            c0: self.c0 - other.c0,
+            c1: self.c1 - other.c1,
+            c2: self.c2 - other.c2,
+        }
+    }
+
+    #[inline]
+    pub fn neg(&self) -> Self {
+        Self {
+            c0: -self.c0,
+            c1: -self.c1,
+            c2: -self.c2,
+        }
+    }
+}
+
+impl<F: ff::Field> CubicExtField<F>
+where
+    Self: CubicExtFieldArith<Base = F>,
+{
+    pub fn mul(&self, rhs: &Self) -> Self {
+        let mut lhs = *self;
+        Self::mul_assign(&mut lhs, rhs);
+        lhs
+    }
+
+    pub fn mul_assign(&mut self, rhs: &Self) {
+        <Self as CubicExtFieldArith>::mul_assign(self, rhs);
+    }
+
+    pub fn square(el: &Self) -> Self {
+        let mut el = *el;
+        Self::square_assign(&mut el);
+        el
+    }
+
+    pub fn square_assign(&mut self) {
+        <Self as CubicExtFieldArith>::square_assign(self);
+    }
+}
+
+impl<F: ExtField> ff::Field for CubicExtField<F>
+where
+    CubicExtField<F>: CubicExtFieldArith<Base = F> + ExtField, // kind of cyclic being `ExtField: Field` but it seems alright
+{
+    const ZERO: Self = Self::zero();
+    const ONE: Self = Self::one();
+
+    fn random(mut rng: impl rand_core::RngCore) -> Self {
+        Self::new(
+            F::random(&mut rng),
+            F::random(&mut rng),
+            F::random(&mut rng),
+        )
+    }
+
+    fn is_zero(&self) -> subtle::Choice {
+        self.c0.is_zero() & self.c1.is_zero()
+    }
+
+    fn square(&self) -> Self {
+        CubicExtField::square(self)
+    }
+
+    fn double(&self) -> Self {
+        self.double()
+    }
+
+    fn sqrt(&self) -> subtle::CtOption<Self> {
+        unimplemented!()
+    }
+
+    fn sqrt_ratio(_: &Self, _: &Self) -> (subtle::Choice, Self) {
+        unimplemented!()
+    }
+
+    fn invert(&self) -> subtle::CtOption<Self> {
+        let c0 = self.c2.mul_by_nonresidue() * self.c1.neg() + self.c0.square();
+        let c1 = self.c2.square().mul_by_nonresidue() - (self.c0 * self.c1);
+        let c2 = self.c1.square() - (self.c0 * self.c2);
+
+        let t = (self.c2 * c1) + (self.c1 * c2);
+        let t = t.mul_by_nonresidue() + (self.c0 * c0);
+
+        t.invert().map(|t| Self {
+            c0: t * c0,
+            c1: t * c1,
+            c2: t * c2,
+        })
+    }
+}
+
+impl<F: ff::Field> subtle::ConditionallySelectable for CubicExtField<F> {
+    fn conditional_select(a: &Self, b: &Self, choice: subtle::Choice) -> Self {
+        CubicExtField {
+            c0: F::conditional_select(&a.c0, &b.c0, choice),
+            c1: F::conditional_select(&a.c1, &b.c1, choice),
+            c2: F::conditional_select(&a.c2, &b.c2, choice),
+        }
+    }
+}
+
+impl<F: ff::Field> subtle::ConstantTimeEq for CubicExtField<F> {
+    fn ct_eq(&self, other: &Self) -> subtle::Choice {
+        self.c0.ct_eq(&other.c0) & self.c1.ct_eq(&other.c1) & self.c2.ct_eq(&other.c2)
+    }
+}
diff --git a/src/ff_ext/mod.rs b/src/ff_ext/mod.rs
index 2766a013..da4d3aa6 100644
--- a/src/ff_ext/mod.rs
+++ b/src/ff_ext/mod.rs
@@ -1,5 +1,7 @@
+pub mod cubic;
 pub mod inverse;
 pub mod jacobi;
+pub mod quadratic;
 use subtle::{Choice, ConstantTimeEq};
 
 pub trait Legendre {
@@ -30,3 +32,11 @@ macro_rules! extend_field_legendre {
         }
     };
 }
+
+pub trait ExtField: ff::Field {
+    const NON_RESIDUE: Self;
+    fn mul_by_nonresidue(&self) -> Self {
+        Self::NON_RESIDUE * self
+    }
+    fn frobenius_map(&mut self, power: usize);
+}
diff --git a/src/ff_ext/quadratic.rs b/src/ff_ext/quadratic.rs
new file mode 100644
index 00000000..3484798d
--- /dev/null
+++ b/src/ff_ext/quadratic.rs
@@ -0,0 +1,330 @@
+use ff::Field;
+use subtle::{Choice, CtOption};
+
+use super::{
+    cubic::{CubicExtField, CubicSparseMul},
+    ExtField,
+};
+
+pub trait QuadSparseMul {
+    type Base: ExtField;
+
+    fn mul_by_014(
+        lhs: &mut QuadExtField<CubicExtField<Self::Base>>,
+        c0: &Self::Base,
+        c1: &Self::Base,
+        c4: &Self::Base,
+    ) where
+        CubicExtField<Self::Base>: CubicSparseMul<Base = Self::Base> + ExtField,
+    {
+        let aa = CubicExtField::mul_by_01(&lhs.c0, c0, c1);
+        let bb = CubicExtField::mul_by_1(&lhs.c1, c4);
+        let t0 = &(lhs.c1 + lhs.c0);
+        let t1 = *c1 + c4;
+        lhs.c1 = CubicExtField::mul_by_01(t0, c0, &t1) - (aa + bb);
+        lhs.c0 = bb.mul_by_nonresidue() + aa;
+    }
+
+    fn mul_by_034(
+        lhs: &mut QuadExtField<CubicExtField<Self::Base>>,
+        c0: &Self::Base,
+        c3: &Self::Base,
+        c4: &Self::Base,
+    ) where
+        CubicExtField<Self::Base>: CubicSparseMul<Base = Self::Base> + ExtField,
+    {
+        let t0 = CubicExtField {
+            c0: lhs.c0.c0 * c0,
+            c1: lhs.c0.c1 * c0,
+            c2: lhs.c0.c2 * c0,
+        };
+        let t1 = CubicExtField::mul_by_01(&lhs.c1, c3, c4);
+        let t2 = lhs.c0 + lhs.c1;
+        let t3 = *c0 + c3;
+        lhs.c1 = CubicExtField::mul_by_01(&t2, &t3, c4) - t0 - t1;
+        lhs.c0 = t0 + t1.mul_by_nonresidue();
+    }
+}
+
+// Algorithm 9 of https://eprint.iacr.org/2012/685.pdf
+pub fn sqrt_algo9<F: ExtField, S: AsRef<[u64]>>(
+    e: &QuadExtField<F>,
+    q_minus_3_over_4: S,
+    q_minus_1_over_2: S,
+) -> subtle::CtOption<QuadExtField<F>>
+where
+    QuadExtField<F>: QuadExtFieldArith<Base = F> + ExtField,
+{
+    if e.is_zero().into() {
+        subtle::CtOption::new(QuadExtField::ZERO, subtle::Choice::from(1))
+    } else {
+        let mut a1 = e.pow(q_minus_3_over_4);
+
+        let alpha = a1.square();
+        let alpha = alpha * e;
+
+        let mut a0 = alpha;
+        a0.frobenius_map(1);
+        let a0 = a0 * alpha;
+
+        let neg1 = QuadExtField::<F> {
+            c0: F::ZERO - F::ONE,
+            c1: F::ZERO,
+        };
+
+        if a0 == neg1 {
+            subtle::CtOption::new(a0, subtle::Choice::from(0))
+        } else {
+            a1.mul_assign(e);
+
+            if alpha == neg1 {
+                a1.mul_assign(&QuadExtField::<F> {
+                    c0: F::ZERO,
+                    c1: F::ONE,
+                });
+            } else {
+                let alpha = alpha + QuadExtField::<F>::ONE;
+                let alpha = alpha.pow(q_minus_1_over_2);
+                a1.mul_assign(&alpha);
+            }
+            subtle::CtOption::new(a1, subtle::Choice::from(1))
+        }
+    }
+}
+
+// Algorithm 10 of https://eprint.iacr.org/2012/685.pdf
+pub fn sqrt_algo10<F: ExtField, S: AsRef<[u64]>>(
+    el: &QuadExtField<F>,
+    precompute_e: &QuadExtField<F>,
+    precompute_f: &QuadExtField<F>,
+    q_minus_1_over_4: S,
+) -> subtle::CtOption<QuadExtField<F>>
+where
+    QuadExtField<F>: QuadExtFieldArith<Base = F> + ExtField,
+{
+    let b = el.pow_vartime(q_minus_1_over_4);
+
+    let b_2 = b.square();
+    let mut b_2_q = b_2;
+    b_2_q.frobenius_map(1);
+
+    let a0 = b_2_q * b_2;
+    let neg1 = QuadExtField::<F> {
+        c0: F::ZERO - F::ONE,
+        c1: F::ZERO,
+    };
+
+    if a0 == neg1 {
+        CtOption::new(a0, Choice::from(0))
+    } else {
+        let mut x = b;
+        x.frobenius_map(1);
+        if x * b == QuadExtField::ONE {
+            let x0 = (b_2 * el).c0.sqrt().unwrap();
+            x.c0.mul_assign(x0);
+            x.c1.mul_assign(x0);
+            CtOption::new(x, Choice::from(1))
+        } else {
+            let x0 = (b_2 * precompute_f * el).sqrt().unwrap();
+            x *= x0 * precompute_e;
+            CtOption::new(x, Choice::from(1))
+        }
+    }
+}
+
+pub enum SQRT<F: Field> {
+    Algorithm9 {
+        q_minus_3_over_4: &'static [u64],
+        q_minus_1_over_2: &'static [u64],
+    },
+    Algorithm10 {
+        precompute_e: QuadExtField<F>,
+        precompute_f: QuadExtField<F>,
+        q_minus_1_over_4: &'static [u64],
+    },
+    Unimplemented,
+}
+
+pub trait QuadExtFieldArith {
+    type Base: ExtField;
+    const SQRT: SQRT<Self::Base> = SQRT::Unimplemented;
+
+    fn mul_assign(lhs: &mut QuadExtField<Self::Base>, rhs: &QuadExtField<Self::Base>) {
+        let v0 = lhs.c0 * rhs.c0;
+        let v1 = lhs.c1 * rhs.c1;
+        lhs.c1 = (lhs.c0 + lhs.c1) * (rhs.c0 + rhs.c1) - (v0 + v1);
+        lhs.c0 = v0 + v1.mul_by_nonresidue();
+    }
+
+    fn square_assign(el: &mut QuadExtField<Self::Base>) {
+        let ab = el.c0 * el.c1;
+        let c0c1 = el.c0 + el.c1;
+        let c0 = (el.c1.mul_by_nonresidue() + el.c0) * c0c1 - ab;
+        el.c1 = ab.double();
+        el.c0 = c0 - ab.mul_by_nonresidue();
+    }
+}
+
+#[cfg(feature = "derive_serde")]
+use serde::{Deserialize, Serialize};
+
+#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
+#[cfg_attr(feature = "derive_serde", derive(Serialize, Deserialize))]
+pub struct QuadExtField<F: ff::Field> {
+    pub(crate) c0: F,
+    pub(crate) c1: F,
+}
+
+impl<F: ff::Field> QuadExtField<F> {
+    #[inline]
+    pub const fn new(c0: F, c1: F) -> Self {
+        Self { c0, c1 }
+    }
+
+    #[inline]
+    pub const fn zero() -> Self {
+        Self {
+            c0: F::ZERO,
+            c1: F::ZERO,
+        }
+    }
+
+    #[inline]
+    pub const fn one() -> Self {
+        Self {
+            c0: F::ONE,
+            c1: F::ZERO,
+        }
+    }
+
+    #[inline]
+    pub fn double(&self) -> Self {
+        Self {
+            c0: self.c0.double(),
+            c1: self.c1.double(),
+        }
+    }
+
+    #[inline]
+    pub fn add(&self, other: &Self) -> Self {
+        Self {
+            c0: self.c0 + other.c0,
+            c1: self.c1 + other.c1,
+        }
+    }
+
+    #[inline]
+    pub fn sub(&self, other: &Self) -> Self {
+        Self {
+            c0: self.c0 - other.c0,
+            c1: self.c1 - other.c1,
+        }
+    }
+
+    #[inline]
+    pub fn neg(&self) -> Self {
+        Self {
+            c0: -self.c0,
+            c1: -self.c1,
+        }
+    }
+
+    #[inline]
+    pub fn conjugate(&mut self) {
+        self.c1 = -self.c1;
+    }
+}
+
+impl<F: ExtField> QuadExtField<F>
+where
+    Self: QuadExtFieldArith<Base = F>,
+{
+    pub fn mul(&self, rhs: &Self) -> Self {
+        let mut lhs = *self;
+        Self::mul_assign(&mut lhs, rhs);
+        lhs
+    }
+
+    pub fn mul_assign(&mut self, rhs: &Self) {
+        <Self as QuadExtFieldArith>::mul_assign(self, rhs);
+    }
+
+    pub fn square(el: &Self) -> Self {
+        let mut el = *el;
+        Self::square_assign(&mut el);
+        el
+    }
+
+    pub fn square_assign(&mut self) {
+        <Self as QuadExtFieldArith>::square_assign(self);
+    }
+
+    pub fn norm(&self) -> F {
+        self.c0.square() - self.c1.square().mul_by_nonresidue()
+    }
+}
+
+impl<F: ExtField> Field for QuadExtField<F>
+where
+    QuadExtField<F>: QuadExtFieldArith<Base = F> + ExtField,
+{
+    const ZERO: Self = Self::zero();
+    const ONE: Self = Self::one();
+
+    fn random(mut rng: impl rand_core::RngCore) -> Self {
+        Self::new(F::random(&mut rng), F::random(&mut rng))
+    }
+
+    fn is_zero(&self) -> subtle::Choice {
+        self.c0.is_zero() & self.c1.is_zero()
+    }
+
+    fn square(&self) -> Self {
+        QuadExtField::square(self)
+    }
+
+    fn double(&self) -> Self {
+        self.double()
+    }
+
+    fn sqrt(&self) -> subtle::CtOption<Self> {
+        match Self::SQRT {
+            SQRT::Algorithm9 {
+                q_minus_3_over_4,
+                q_minus_1_over_2,
+            } => sqrt_algo9(self, q_minus_3_over_4, q_minus_1_over_2),
+            SQRT::Algorithm10 {
+                precompute_e,
+                precompute_f,
+                q_minus_1_over_4,
+            } => sqrt_algo10(self, &precompute_e, &precompute_f, q_minus_1_over_4),
+            SQRT::Unimplemented => unimplemented!(),
+        }
+    }
+
+    fn sqrt_ratio(_: &Self, _: &Self) -> (subtle::Choice, Self) {
+        unimplemented!()
+    }
+
+    fn invert(&self) -> subtle::CtOption<Self> {
+        self.norm().invert().map(|t| Self {
+            c0: self.c0 * t,
+            c1: self.c1 * -t,
+        })
+    }
+}
+
+impl<F: ff::Field> subtle::ConditionallySelectable for QuadExtField<F> {
+    fn conditional_select(a: &Self, b: &Self, choice: subtle::Choice) -> Self {
+        QuadExtField {
+            c0: F::conditional_select(&a.c0, &b.c0, choice),
+            c1: F::conditional_select(&a.c1, &b.c1, choice),
+        }
+    }
+}
+
+impl<F: ff::Field> subtle::ConstantTimeEq for QuadExtField<F> {
+    fn ct_eq(&self, other: &Self) -> subtle::Choice {
+        self.c0.ct_eq(&other.c0) & self.c1.ct_eq(&other.c1)
+    }
+}
diff --git a/src/grumpkin/curve.rs b/src/grumpkin/curve.rs
index 9a22024b..3d956490 100644
--- a/src/grumpkin/curve.rs
+++ b/src/grumpkin/curve.rs
@@ -10,9 +10,8 @@ use crate::group::{prime::PrimeCurveAffine, Group, GroupEncoding};
 use crate::grumpkin::Fq;
 use crate::grumpkin::Fr;
 use crate::{
-    endo, impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    new_curve_impl,
+    endo, impl_binops_additive, impl_binops_additive_specify_output, impl_binops_multiplicative,
+    impl_binops_multiplicative_mixed, new_curve_impl,
 };
 use crate::{Coordinates, CurveAffine, CurveExt};
 use core::cmp;
diff --git a/src/msm.rs b/src/msm.rs
index 25af9711..9b1efe68 100644
--- a/src/msm.rs
+++ b/src/msm.rs
@@ -558,6 +558,6 @@ mod test {
 
     #[test]
     fn test_msm_cross() {
-        run_msm_cross::<G1Affine>(14, 22);
+        run_msm_cross::<G1Affine>(14, 18);
     }
 }
diff --git a/src/pluto_eris/curve.rs b/src/pluto_eris/curve.rs
index 6a394887..23acac24 100644
--- a/src/pluto_eris/curve.rs
+++ b/src/pluto_eris/curve.rs
@@ -16,9 +16,8 @@ use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 use serde::{Deserialize, Serialize};
 
 use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    new_curve_impl,
+    impl_binops_additive, impl_binops_additive_specify_output, impl_binops_multiplicative,
+    impl_binops_multiplicative_mixed, new_curve_impl,
 };
 
 const G1_GENERATOR_X: Fp = Fp::from_raw([
diff --git a/src/pluto_eris/engine.rs b/src/pluto_eris/engine.rs
index 7b6a8226..4975c424 100644
--- a/src/pluto_eris/engine.rs
+++ b/src/pluto_eris/engine.rs
@@ -1,18 +1,23 @@
 #![allow(clippy::suspicious_arithmetic_impl)]
+
 use crate::ff::PrimeField;
+use crate::ff_ext::quadratic::QuadSparseMul;
+use crate::ff_ext::ExtField;
 use crate::group::cofactor::CofactorCurveAffine;
 use crate::group::Group;
-use crate::pairing::{Engine, MillerLoopResult, MultiMillerLoop, PairingCurveAffine};
 use crate::pluto_eris::curve::*;
 use crate::pluto_eris::fp::*;
 use crate::pluto_eris::fp12::*;
 use crate::pluto_eris::fp2::*;
 use crate::pluto_eris::fp6::FROBENIUS_COEFF_FP6_C1;
-use crate::pluto_eris::fq::*;
+use crate::pluto_eris::fq::Fq;
 use core::borrow::Borrow;
 use core::iter::Sum;
-use core::ops::{Add, Mul, MulAssign, Neg, Sub};
+use core::ops::{Add, Mul, Neg, Sub};
+use ff::Field;
+use pairing::{Engine, MillerLoopResult, MultiMillerLoop, PairingCurveAffine};
 use rand_core::RngCore;
+use std::ops::MulAssign;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq};
 
 /// Adaptation of Algorithm 1, https://eprint.iacr.org/2013/722.pdf
@@ -26,467 +31,89 @@ const NEG_SIX_U_PLUS_2_NAF: [i8; 114] = [
     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1,
 ];
 
-/// Value of (57/(u + 3))^((p - 1)/2) where u^2 + 5 = 0 in Fp2.
-const XI_TO_P_MINUS_1_OVER_2: Fp2 = Fp2 {
-    c0: Fp::from_raw([
-        0x54cf5ad1c0926216,
-        0x186c1f3ce4a46d4e,
-        0x9c23800ce9c9452f,
-        0x50e0d09ff6d6c08b,
-        0x7cf421e4d46f6666,
-        0x678664ba4b6d8343,
-        0x21cc26d5de0f80f4,
-    ]),
-
-    c1: Fp::from_raw([
-        0xc0505f4c260e91f4,
-        0xe7bbd15f10723657,
-        0xb4b3e0c35358097e,
-        0x87c56f42a558750d,
-        0x4b7211d23f34f0ae,
-        0xf6839d29e2f0d250,
-        0x16ebe8b2e12a1106,
-    ]),
-};
-
-impl PairingCurveAffine for G1Affine {
-    type Pair = G2Affine;
-    type PairingResult = Gt;
-
-    fn pairing_with(&self, other: &Self::Pair) -> Self::PairingResult {
-        Pluto::pairing(self, other)
-    }
-}
-
-impl PairingCurveAffine for G2Affine {
-    type Pair = G1Affine;
-    type PairingResult = Gt;
-
-    fn pairing_with(&self, other: &Self::Pair) -> Self::PairingResult {
-        Pluto::pairing(other, self)
-    }
-}
-
-#[derive(Copy, Clone, Debug, Default)]
-pub struct Gt(pub(crate) Fp12);
-
-impl std::fmt::Display for Gt {
-    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
-        write!(f, "{:?}", self)
-    }
-}
-
-impl ConstantTimeEq for Gt {
-    fn ct_eq(&self, other: &Self) -> Choice {
-        self.0.ct_eq(&other.0)
-    }
-}
-
-impl ConditionallySelectable for Gt {
-    fn conditional_select(a: &Self, b: &Self, choice: Choice) -> Self {
-        Gt(Fp12::conditional_select(&a.0, &b.0, choice))
-    }
-}
-
-impl Eq for Gt {}
-impl PartialEq for Gt {
-    #[inline]
-    fn eq(&self, other: &Self) -> bool {
-        bool::from(self.ct_eq(other))
-    }
-}
-
-impl Gt {
-    /// Returns the group identity, which is $1$.
-    pub const fn identity() -> Gt {
-        Gt(Fp12::one())
-    }
-
-    /// Doubles this group element.
-    pub fn double(&self) -> Gt {
-        Gt(self.0.square())
-    }
-}
-
-impl<'a> Neg for &'a Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn neg(self) -> Gt {
-        // The element is unitary, so we just conjugate.
-        let mut u = self.0;
-        u.conjugate();
-        Gt(u)
-    }
-}
-
-impl Neg for Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn neg(self) -> Gt {
-        -&self
-    }
-}
-
-impl<'a, 'b> Add<&'b Gt> for &'a Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn add(self, rhs: &'b Gt) -> Gt {
-        Gt(self.0 * rhs.0)
-    }
-}
-
-impl<'a, 'b> Sub<&'b Gt> for &'a Gt {
-    type Output = Gt;
-
-    #[inline]
-    fn sub(self, rhs: &'b Gt) -> Gt {
-        self + (-rhs)
-    }
-}
-
-impl<'a, 'b> Mul<&'b Fq> for &'a Gt {
-    type Output = Gt;
-
-    fn mul(self, other: &'b Fq) -> Self::Output {
-        let mut acc = Gt::identity();
-
-        for bit in other
-            .to_repr()
-            .as_ref()
-            .iter()
-            .rev()
-            .flat_map(|byte| (0..8).rev().map(move |i| Choice::from((byte >> i) & 1u8)))
-            .skip(1)
-        {
-            acc = acc.double();
-            acc = Gt::conditional_select(&acc, &(acc + self), bit);
-        }
-
-        acc
-    }
-}
-
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-};
-impl_binops_additive!(Gt, Gt);
-impl_binops_multiplicative!(Gt, Fq);
-
-impl<T> Sum<T> for Gt
-where
-    T: Borrow<Gt>,
-{
-    fn sum<I>(iter: I) -> Self
-    where
-        I: Iterator<Item = T>,
-    {
-        iter.fold(Self::identity(), |acc, item| acc + item.borrow())
-    }
-}
-
-impl Group for Gt {
-    type Scalar = Fq;
-
-    fn random(_: impl RngCore) -> Self {
-        unimplemented!();
-    }
-
-    fn identity() -> Self {
-        Self::identity()
-    }
-
-    fn generator() -> Self {
-        unimplemented!();
-    }
-
-    fn is_identity(&self) -> Choice {
-        self.ct_eq(&Self::identity())
-    }
-
-    #[must_use]
-    fn double(&self) -> Self {
-        self.double()
-    }
-}
-
-/// Points of G2 in Jacobian coordinates.
-/// These are points lie in the twisted curve E'(Fp2).
-#[derive(Clone, Debug)]
-pub struct G2Prepared {
-    pub(crate) coeffs: Vec<(Fp2, Fp2, Fp2)>,
-    pub(crate) infinity: bool,
-}
-
-impl G2Prepared {
-    /// Returns true if `self` is the infinity point.
-    pub fn is_zero(&self) -> bool {
-        self.infinity
-    }
-
-    /// Prepares a G2 point in affine coordinates.
-    pub fn from_affine(q: G2Affine) -> Self {
-        if bool::from(q.is_identity()) {
-            return G2Prepared {
-                coeffs: vec![],
-                infinity: true,
-            };
-        }
-
-        /// Adaptation of Algorithm 26, https://eprint.iacr.org/2010/354.pdf
-        fn doubling_step(r: &mut G2) -> (Fp2, Fp2, Fp2) {
-            let mut tmp0 = r.x;
-            tmp0.square_assign();
-
-            let mut tmp1 = r.y;
-            tmp1.square_assign();
-
-            let mut tmp2 = tmp1;
-            tmp2.square_assign();
-
-            let mut tmp3 = tmp1;
-            tmp3 += &r.x;
-            tmp3.square_assign();
-            tmp3 -= &tmp0;
-            tmp3 -= &tmp2;
-            tmp3.double_assign();
-
-            let mut tmp4 = tmp0;
-            tmp4.double_assign();
-            tmp4 += &tmp0;
-
-            let mut tmp6 = r.x;
-            tmp6 += &tmp4;
-
-            let mut tmp5 = tmp4;
-            tmp5.square_assign();
-
-            let mut zsquared = r.z;
-            zsquared.square_assign();
-
-            r.x = tmp5;
-            r.x -= &tmp3;
-            r.x -= &tmp3;
-
-            r.z += &r.y;
-            r.z.square_assign();
-            r.z -= &tmp1;
-            r.z -= &zsquared;
-
-            r.y = tmp3;
-            r.y -= &r.x;
-            r.y.mul_assign(&tmp4);
-
-            tmp2.double_assign();
-            tmp2.double_assign();
-            tmp2.double_assign();
-
-            r.y -= &tmp2;
-
-            // up to here everything was by algorithm, line 11
-            // use R instead of new T
-
-            // tmp3 is the first part of line 12
-            tmp3 = tmp4;
-            tmp3.mul_assign(&zsquared);
-            tmp3.double_assign();
-            tmp3 = tmp3.neg();
-
-            // tmp6 is from line 14
-            tmp6.square_assign();
-            tmp6 -= &tmp0;
-            tmp6 -= &tmp5;
-
-            tmp1.double_assign();
-            tmp1.double_assign();
-
-            tmp6 -= &tmp1;
-
-            // tmp0 is the first part of line 16
-            tmp0 = r.z;
-            tmp0.mul_assign(&zsquared);
-            tmp0.double_assign();
-
-            (tmp0, tmp3, tmp6)
-        }
-
-        // Adaptation of Algorithm 27, https://eprint.iacr.org/2010/354.pdf
-        fn addition_step(r: &mut G2, q: &G2Affine) -> (Fp2, Fp2, Fp2) {
-            let mut zsquared = r.z;
-            zsquared.square_assign();
-
-            let mut ysquared = q.y;
-            ysquared.square_assign();
-
-            // t0 corresponds to line 1
-            let mut t0 = zsquared;
-            t0.mul_assign(&q.x);
-
-            // t1 corresponds to lines 2 and 3
-            let mut t1 = q.y;
-            t1 += &r.z;
-            t1.square_assign();
-            t1 -= &ysquared;
-            t1 -= &zsquared;
-            t1.mul_assign(&zsquared);
-
-            // t2 corresponds to line 4
-            let mut t2 = t0;
-            t2 -= &r.x;
-
-            // t3 corresponds to line 5
-            let mut t3 = t2;
-            t3.square_assign();
-
-            // t4 corresponds to line 6
-            let mut t4 = t3;
-            t4.double_assign();
-            t4.double_assign();
-
-            // t5 corresponds to line 7
-            let mut t5 = t4;
-            t5.mul_assign(&t2);
-
-            // t6 corresponds to line 8
-            let mut t6 = t1;
-            t6 -= &r.y;
-            t6 -= &r.y;
-
-            // t9 corresponds to line 9
-            let mut t9 = t6;
-            t9.mul_assign(&q.x);
-
-            // corresponds to line 10
-            let mut t7 = t4;
-            t7.mul_assign(&r.x);
-
-            // corresponds to line 11, but assigns to r.x instead of T.x
-            r.x = t6;
-            r.x.square_assign();
-            r.x -= &t5;
-            r.x -= &t7;
-            r.x -= &t7;
-
-            // corresponds to line 12, but assigns to r.z instead of T.z
-            r.z += &t2;
-            r.z.square_assign();
-            r.z -= &zsquared;
-            r.z -= &t3;
-
-            // corresponds to line 13
-            let mut t10 = q.y;
-            t10 += &r.z;
-
-            // corresponds to line 14
-            let mut t8 = t7;
-            t8 -= &r.x;
-            t8.mul_assign(&t6);
-
-            // corresponds to line 15
-            t0 = r.y;
-            t0.mul_assign(&t5);
-            t0.double_assign();
-
-            // corresponds to line 12, but assigns to r.y instead of T.y
-            r.y = t8;
-            r.y -= &t0;
-
-            // corresponds to line 17
-            t10.square_assign();
-            t10 -= &ysquared;
-
-            let mut ztsquared = r.z;
-            ztsquared.square_assign();
-
-            t10 -= &ztsquared;
-
-            // corresponds to line 18
-            t9.double_assign();
-            t9 -= &t10;
-
-            // t10 = 2*Zt from Algo 27, line 19
-            t10 = r.z;
-            t10.double_assign();
-
-            // t1 = first multiplicator of line 21
-            t6 = t6.neg();
-
-            t1 = t6;
-            t1.double_assign();
-
-            // t9 corresponds to t9 from Algo 27
-            (t10, t1, t9)
-        }
-
-        let mut coeffs = vec![];
-        let mut r: G2 = q.into();
+crate::impl_gt!(Gt, Fp12, Fq);
+crate::impl_miller_loop_components!(Pluto, G1, G1Affine, G2, G2Affine, Fp12, Gt, Fq);
+
+pub fn multi_miller_loop(terms: &[(&G1Affine, &G2Affine)]) -> Fp12 {
+    let terms = terms
+        .iter()
+        .filter_map(|&(p, q)| {
+            if bool::from(p.is_identity()) || bool::from(q.is_identity()) {
+                None
+            } else {
+                Some((*p, *q))
+            }
+        })
+        .collect::<Vec<_>>();
 
-        let mut negq = q;
-        negq = -negq;
+    let mut f = Fp12::one();
+    let mut r = terms.iter().map(|(_, q)| q.to_curve()).collect::<Vec<_>>();
 
-        coeffs.push(doubling_step(&mut r));
+    for (i, x) in NEG_SIX_U_PLUS_2_NAF.iter().rev().skip(1).enumerate() {
+        (i != 0).then(|| f.square_assign());
 
-        let last_position = NEG_SIX_U_PLUS_2_NAF.len() - 2;
-        match NEG_SIX_U_PLUS_2_NAF[last_position] {
-            1 => {
-                coeffs.push(addition_step(&mut r, &q));
-            }
-            -1 => {
-                coeffs.push(addition_step(&mut r, &negq));
-            }
-            _ => (),
+        for ((p, _), r) in terms.iter().zip(r.iter_mut()) {
+            double(&mut f, r, p);
         }
 
-        for i in (0..last_position).rev() {
-            coeffs.push(doubling_step(&mut r));
-
-            match NEG_SIX_U_PLUS_2_NAF[i] {
-                1 => {
-                    coeffs.push(addition_step(&mut r, &q));
-                }
-                -1 => {
-                    coeffs.push(addition_step(&mut r, &negq));
+        match x {
+            &val @ (1 | -1) => {
+                for ((p, q), r) in terms.iter().zip(r.iter_mut()) {
+                    if val == 1 {
+                        add(&mut f, r, q, p);
+                    } else {
+                        add(&mut f, r, &q.neg(), p);
+                    }
                 }
-                _ => continue,
             }
-        }
-
-        let mut neg_r = r;
-        neg_r = -neg_r;
-
-        let mut q1 = q;
-
-        q1.x.c1 = q1.x.c1.neg();
+            _ => continue,
+        }
+    }
+
+    /// Value of (57/(u + 3))^((p - 1)/2) where u^2 + 5 = 0 in Fp2.
+    const XI_TO_P_MINUS_1_OVER_2: Fp2 = Fp2 {
+        c0: Fp::from_raw([
+            0x54cf5ad1c0926216,
+            0x186c1f3ce4a46d4e,
+            0x9c23800ce9c9452f,
+            0x50e0d09ff6d6c08b,
+            0x7cf421e4d46f6666,
+            0x678664ba4b6d8343,
+            0x21cc26d5de0f80f4,
+        ]),
+
+        c1: Fp::from_raw([
+            0xc0505f4c260e91f4,
+            0xe7bbd15f10723657,
+            0xb4b3e0c35358097e,
+            0x87c56f42a558750d,
+            0x4b7211d23f34f0ae,
+            0xf6839d29e2f0d250,
+            0x16ebe8b2e12a1106,
+        ]),
+    };
+
+    for ((p, q), r) in terms.iter().zip(r.iter_mut()) {
+        let mut q1: G2Affine = *q;
+        q1.x.conjugate();
         q1.x.mul_assign(&FROBENIUS_COEFF_FP6_C1[1]);
-
-        q1.y.c1 = q1.y.c1.neg();
+        q1.y.conjugate();
         q1.y.mul_assign(&XI_TO_P_MINUS_1_OVER_2);
+        add(&mut f, r, &q1.neg(), p);
+    }
 
-        coeffs.push(addition_step(&mut neg_r, &q1));
-
-        let mut minusq2 = q;
+    for ((p, q), r) in terms.iter().zip(r.iter_mut()) {
+        let mut minusq2: G2Affine = *q;
         minusq2.x.mul_assign(&FROBENIUS_COEFF_FP6_C1[2]);
-
-        coeffs.push(addition_step(&mut neg_r, &minusq2));
-
-        G2Prepared {
-            coeffs,
-            infinity: false,
-        }
+        add(&mut f, r, &minusq2.neg(), p);
     }
-}
 
-impl From<G2Affine> for G2Prepared {
-    fn from(q: G2Affine) -> G2Prepared {
-        G2Prepared::from_affine(q)
-    }
+    f
 }
 
-impl MillerLoopResult for Gt {
-    type Gt = Self;
+impl MillerLoopResult for Fp12 {
+    type Gt = Gt;
+
     fn final_exponentiation(&self) -> Gt {
         fn exp_by_x(f: &mut Fp12) {
             let x = NEG_PLUTO_U;
@@ -501,8 +128,8 @@ impl MillerLoopResult for Gt {
             *f = res;
         }
 
-        let r = self.0;
-        let mut f1 = self.0;
+        let r = *self;
+        let mut f1 = *self;
         f1.conjugate();
 
         Gt(r.invert()
@@ -589,366 +216,25 @@ impl MillerLoopResult for Gt {
             .unwrap())
     }
 }
-impl MultiMillerLoop for Pluto {
-    /// The prepared form of `Self::G2Affine`.
-    type G2Prepared = G2Prepared;
-
-    /// The type returned by `Engine::miller_loop`.
-    type Result = Gt;
-    fn multi_miller_loop(terms: &[(&G1Affine, &G2Prepared)]) -> Self::Result {
-        let mut pairs = vec![];
-        for &(p, q) in terms {
-            if !bool::from(p.is_identity()) && !q.is_zero() {
-                pairs.push((p, q.coeffs.iter()));
-            }
-        }
-
-        // Final steps of the line function on prepared coefficients
-        fn ell(f: &mut Fp12, coeffs: &(Fp2, Fp2, Fp2), p: &G1Affine) {
-            let mut c0 = coeffs.0;
-            let mut c1 = coeffs.1;
-
-            c0.c0.mul_assign(&p.y);
-            c0.c1.mul_assign(&p.y);
-
-            c1.c0.mul_assign(&p.x);
-            c1.c1.mul_assign(&p.x);
-
-            // Sparse multiplication in Fq12
-            f.mul_by_034(&c0, &c1, &coeffs.2);
-        }
-
-        let mut f = Fp12::one();
-
-        for &mut (p, ref mut coeffs) in &mut pairs {
-            ell(&mut f, coeffs.next().unwrap(), p);
-        }
-
-        // length - 2
-        let len_min2 = NEG_SIX_U_PLUS_2_NAF.len() - 2;
 
-        if NEG_SIX_U_PLUS_2_NAF[len_min2] != 0 {
-            for &mut (p, ref mut coeffs) in &mut pairs {
-                ell(&mut f, coeffs.next().unwrap(), p);
-            }
-        }
-
-        for x in NEG_SIX_U_PLUS_2_NAF[..len_min2].iter().rev() {
-            f.square_assign();
-
-            for &mut (p, ref mut coeffs) in &mut pairs {
-                ell(&mut f, coeffs.next().unwrap(), p);
-            }
-            if *x != 0 {
-                for &mut (p, ref mut coeffs) in &mut pairs {
-                    ell(&mut f, coeffs.next().unwrap(), p);
-                }
-            }
-        }
-
-        f.conjugate();
-
-        for &mut (p, ref mut coeffs) in &mut pairs {
-            ell(&mut f, coeffs.next().unwrap(), p);
-        }
-
-        for &mut (p, ref mut coeffs) in &mut pairs {
-            ell(&mut f, coeffs.next().unwrap(), p);
-        }
-
-        for &mut (_p, ref mut coeffs) in &mut pairs {
-            assert_eq!(coeffs.next(), None);
-        }
-
-        Gt(f)
-    }
-}
-
-/// Pluto pairing-friendly curve. See: https://github.com/daira/pluto-eris
-#[derive(Clone, Debug)]
-pub struct Pluto;
-
-impl Engine for Pluto {
-    type Fr = Fq;
-    type G1 = G1;
-    type G1Affine = G1Affine;
-    type G2 = G2;
-    type G2Affine = G2Affine;
-    type Gt = Gt;
-
-    fn pairing(p: &Self::G1Affine, q: &Self::G2Affine) -> Self::Gt {
-        let q = G2Prepared::from_affine(*q);
-        let terms: &[(&G1Affine, &G2Prepared)] = &[(p, &q)];
-        let u = Self::multi_miller_loop(terms);
-        u.final_exponentiation()
-    }
+// Final steps of the line function on prepared coefficients
+fn ell(f: &mut Fp12, coeffs: &(Fp2, Fp2, Fp2), p: &G1Affine) {
+    let mut c0 = coeffs.0;
+    let mut c1 = coeffs.1;
+    c0.c0.mul_assign(&p.y);
+    c0.c1.mul_assign(&p.y);
+    c1.c0.mul_assign(&p.x);
+    c1.c1.mul_assign(&p.x);
+    Fp12::mul_by_034(f, &c0, &c1, &coeffs.2);
 }
 
 #[cfg(test)]
-use rand::SeedableRng;
-#[cfg(test)]
-use rand_xorshift::XorShiftRng;
-
-#[test]
-fn test_pairing() {
-    use ff::Field;
-
-    let g1 = G1::generator();
-    let mut g2 = G2::generator();
-    g2 = g2.double();
-    let pair12 = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    let mut g1 = G1::generator();
-    let g2 = G2::generator();
-    g1 = g1.double();
-    let pair21 = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    assert_eq!(pair12, pair21);
-
-    let g1 = G1::generator();
-    let mut g2 = G2::generator();
-    g2 = g2.double().double();
-    let pair12 = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    let mut g1 = G1::generator();
-    let mut g2 = G2::generator();
-    g1 = g1.double();
-    g2 = g2.double();
-    let pair21 = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-    assert_eq!(pair12, pair21);
-
-    let mut rng = XorShiftRng::from_seed([
-        0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-    for _ in 0..100 {
-        let a = Fq::random(&mut rng);
-        let b = Fq::random(&mut rng);
-
-        let mut g1 = G1::generator();
-        g1.mul_assign(a);
-
-        let mut g2 = G2::generator();
-        g1.mul_assign(b);
-
-        let pair_ab = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-        g1 = G1::generator();
-        g1.mul_assign(b);
-
-        g2 = G2::generator();
-        g1.mul_assign(a);
-
-        let pair_ba = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(g2));
-
-        assert_eq!(pair_ab, pair_ba);
-    }
-}
-
-#[test]
-fn tricking_miller_loop_result() {
-    assert_eq!(
-        Pluto::multi_miller_loop(&[(&G1Affine::identity(), &G2Affine::generator().into())]).0,
-        Fp12::one()
-    );
-    assert_eq!(
-        Pluto::multi_miller_loop(&[(&G1Affine::generator(), &G2Affine::identity().into())]).0,
-        Fp12::one()
-    );
-    assert_ne!(
-        Pluto::multi_miller_loop(&[
-            (&G1Affine::generator(), &G2Affine::generator().into()),
-            (&-G1Affine::generator(), &G2Affine::generator().into())
-        ])
-        .0,
-        Fp12::one()
-    );
-    assert_eq!(
-        Pluto::multi_miller_loop(&[
-            (&G1Affine::generator(), &G2Affine::generator().into()),
-            (&-G1Affine::generator(), &G2Affine::generator().into())
-        ])
-        .final_exponentiation(),
-        Gt::identity()
-    );
-}
-
-#[test]
-fn random_bilinearity_tests() {
+mod test {
+    use super::super::{Fq, Pluto, G1, G2};
+    use super::{multi_miller_loop, Fp12, G1Affine, G2Affine, Gt};
     use ff::Field;
-
-    let mut rng = XorShiftRng::from_seed([
-        0x55, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    for _ in 0..10 {
-        let mut a = G1::generator();
-        let ka = Fq::random(&mut rng);
-        a.mul_assign(ka);
-
-        let mut b = G2::generator();
-        let kb = Fq::random(&mut rng);
-        b.mul_assign(kb);
-
-        let c = Fq::random(&mut rng);
-        let d = Fq::random(&mut rng);
-
-        let mut ac = a;
-        ac.mul_assign(c);
-
-        let mut ad = a;
-        ad.mul_assign(d);
-
-        let mut bc = b;
-        bc.mul_assign(c);
-
-        let mut bd = b;
-        bd.mul_assign(d);
-
-        let acbd = Pluto::pairing(&G1Affine::from(ac), &G2Affine::from(bd));
-        let adbc = Pluto::pairing(&G1Affine::from(ad), &G2Affine::from(bc));
-
-        let mut cd = c;
-        cd.mul_assign(&d);
-
-        cd *= Fq([1, 0, 0, 0, 0, 0, 0]);
-
-        let abcd = Gt(Pluto::pairing(&G1Affine::from(a), &G2Affine::from(b))
-            .0
-            .pow_vartime(cd.0));
-
-        assert_eq!(acbd, adbc);
-        assert_eq!(acbd, abcd);
-    }
-}
-
-#[test]
-pub fn engine_tests() {
-    use ff::Field;
-
-    let mut rng = XorShiftRng::from_seed([
-        0x56, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    for _ in 0..10 {
-        let a = G1Affine::from(G1::random(&mut rng));
-        let b = G2Affine::from(G2::random(&mut rng));
-
-        assert!(a.pairing_with(&b) == b.pairing_with(&a));
-        assert!(a.pairing_with(&b) == Pluto::pairing(&a, &b));
-    }
-
-    for _ in 0..10 {
-        let z1 = G1Affine::identity();
-        let z2 = G2Prepared::from(G2Affine::identity());
-
-        let a = G1Affine::from(G1::random(&mut rng));
-        let b = G2Prepared::from(G2Affine::from(G2::random(&mut rng)));
-        let c = G1Affine::from(G1::random(&mut rng));
-        let d = G2Prepared::from(G2Affine::from(G2::random(&mut rng)));
-
-        assert_eq!(
-            Fp12::ONE,
-            Pluto::multi_miller_loop(&[(&z1, &b)])
-                .final_exponentiation()
-                .0,
-        );
-
-        assert_eq!(
-            Fp12::ONE,
-            Pluto::multi_miller_loop(&[(&a, &z2)])
-                .final_exponentiation()
-                .0,
-        );
-
-        assert_eq!(
-            Pluto::multi_miller_loop(&[(&z1, &b), (&c, &d)]).final_exponentiation(),
-            Pluto::multi_miller_loop(&[(&a, &z2), (&c, &d)]).final_exponentiation(),
-        );
-
-        assert_eq!(
-            Pluto::multi_miller_loop(&[(&a, &b), (&z1, &d)]).final_exponentiation(),
-            Pluto::multi_miller_loop(&[(&a, &b), (&c, &z2)]).final_exponentiation(),
-        );
-    }
-}
-
-#[test]
-fn random_miller_loop_tests() {
-    let mut rng = XorShiftRng::from_seed([
-        0x58, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    // Exercise a double miller loop
-    for _ in 0..10 {
-        let a = G1Affine::from(G1::random(&mut rng));
-        let b = G2Affine::from(G2::random(&mut rng));
-        let c = G1Affine::from(G1::random(&mut rng));
-        let d = G2Affine::from(G2::random(&mut rng));
-
-        let ab = Pluto::pairing(&a, &b);
-        let cd = Pluto::pairing(&c, &d);
-
-        let mut abcd = ab;
-        abcd = Gt(abcd.0 * cd.0);
-
-        let b = G2Prepared::from(b);
-        let d = G2Prepared::from(d);
-
-        let abcd_with_double_loop =
-            Pluto::multi_miller_loop(&[(&a, &b), (&c, &d)]).final_exponentiation();
-
-        assert_eq!(abcd, abcd_with_double_loop);
-    }
-}
-
-#[test]
-pub fn multi_miller_final_exp_tests() {
-    use ff::Field;
-
-    let g1 = G1::generator();
-    let g2 = G2::generator();
-
-    let mut rng = XorShiftRng::from_seed([
-        0x56, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-        0xe5,
-    ]);
-
-    for _ in 0..10 {
-        let s = Fq::random(&mut rng);
-
-        let mut s_g1 = g1;
-        s_g1.mul_assign(s);
-
-        let mut s_g2 = g2;
-        s_g2.mul_assign(s);
-
-        let s_g2_prepared = G2Prepared::from(G2Affine::from(s_g2));
-        let g2_prepared = G2Prepared::from(G2Affine::from(g2));
-
-        let (term_1, term_2) = (
-            (&G1Affine::from(g1), &s_g2_prepared),
-            (&-G1Affine::from(s_g1), &g2_prepared),
-        );
-
-        let terms = &[term_1, term_2];
-
-        assert!(
-            bool::from(
-                Pluto::multi_miller_loop(&terms[..])
-                    .final_exponentiation()
-                    .is_identity(),
-            ),
-            "trivial pairing check failed"
-        );
-
-        let lhs = Pluto::pairing(&G1Affine::from(g1), &G2Affine::from(s_g2));
-        let rhs = Pluto::pairing(&G1Affine::from(s_g1), &G2Affine::from(g2));
-
-        assert_eq!(lhs, rhs, "failed trivial check");
-    }
+    use group::{prime::PrimeCurveAffine, Curve, Group};
+    use pairing::{Engine, MillerLoopResult, PairingCurveAffine};
+    use rand_core::OsRng;
+    crate::test_pairing!(Pluto, G1, G1Affine, G2, G2Affine, Fp12, Gt, Fq);
 }
diff --git a/src/pluto_eris/fp.rs b/src/pluto_eris/fp.rs
index 40048d60..53d60ff3 100644
--- a/src/pluto_eris/fp.rs
+++ b/src/pluto_eris/fp.rs
@@ -1,15 +1,9 @@
+use crate::ff_ext::ExtField;
 use core::convert::TryInto;
 use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     pluto_eris_fp,
     Fp,
@@ -20,14 +14,30 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fp);
-impl_binops_calls!(Fp);
-impl_binops_additive!(Fp, Fp);
-impl_binops_multiplicative!(Fp, Fp);
-field_bits!(Fp);
-serialize_deserialize_primefield!(Fp);
+crate::extend_field_legendre!(Fp);
+crate::impl_binops_calls!(Fp);
+crate::impl_binops_additive!(Fp, Fp);
+crate::impl_binops_multiplicative!(Fp, Fp);
+crate::field_bits!(Fp);
+crate::serialize_deserialize_primefield!(Fp);
 crate::impl_from_u64!(Fp);
 
+impl ExtField for Fp {
+    const NON_RESIDUE: Self = Fp::from_raw([
+        0x9ffffcd2fffffffc,
+        0xa2a7e8c30006b945,
+        0xe4a7a5fe8fadffd6,
+        0x443f9a5cda8a6c7b,
+        0xa803ca76f439266f,
+        0x0130e0000d7f70e4,
+        0x2400000000002400,
+    ]);
+    fn mul_by_nonresidue(&self) -> Self {
+        (self.double().double() + self).neg()
+    }
+    fn frobenius_map(&mut self, _: usize) {}
+}
+
 #[cfg(test)]
 mod test {
 
@@ -42,4 +52,12 @@ mod test {
     crate::field_testing_suite!(Fp, "sqrt");
     crate::field_testing_suite!(Fp, "zeta");
     crate::field_testing_suite!(Fp, "from_uniform_bytes", 64, 72, 112);
+
+    #[test]
+    fn test_fq_mul_nonresidue() {
+        let e = Fp::random(rand_core::OsRng);
+        let a0 = e.mul_by_nonresidue();
+        let a1 = e * Fp::NON_RESIDUE;
+        assert_eq!(a0, a1);
+    }
 }
diff --git a/src/pluto_eris/fp12.rs b/src/pluto_eris/fp12.rs
index 9b2bb892..1b64080a 100644
--- a/src/pluto_eris/fp12.rs
+++ b/src/pluto_eris/fp12.rs
@@ -1,209 +1,42 @@
 use super::fp::Fp;
 use super::fp2::Fp2;
 use super::fp6::Fp6;
-use crate::{ff::Field, impl_tower2_common};
-use core::ops::{Add, Neg, Sub};
-use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
+use crate::ff_ext::{
+    quadratic::{QuadExtField, QuadExtFieldArith, QuadSparseMul},
+    ExtField,
+};
+use ff::Field;
 
 /// -GAMMA is a quadratic non-residue in Fp6. Fp12 = Fp6[X]/(X^2 + GAMMA)
 /// We introduce the variable w such that w^2 = -GAMMA
 /// GAMMA = - v
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_calls, impl_binops_multiplicative, impl_binops_multiplicative_mixed,
-    impl_sub_binop_specify_output, impl_sum_prod,
-};
-impl_binops_additive!(Fp12, Fp12);
-impl_binops_multiplicative!(Fp12, Fp12);
-impl_tower2_common!(Fp6, Fp12);
-impl_binops_calls!(Fp12);
-impl_sum_prod!(Fp12);
-
-impl Fp12 {
-    pub fn mul_assign(&mut self, other: &Self) {
-        let t0 = self.c0 * other.c0;
-        let mut t1 = self.c1 * other.c1;
-        let t2 = other.c0 + other.c1;
-
-        self.c1 += &self.c0;
-        self.c1 *= &t2;
-        self.c1 -= &t0;
-        self.c1 -= &t1;
-
-        t1.mul_by_nonresidue();
-        self.c0 = t0 + t1;
-    }
+pub type Fp12 = QuadExtField<Fp6>;
 
-    pub fn square_assign(&mut self) {
-        let mut ab = self.c0 * self.c1;
-
-        let c0c1 = self.c0 + self.c1;
+impl QuadExtFieldArith for Fp12 {
+    type Base = Fp6;
+}
 
-        let mut c0 = self.c1;
-        c0.mul_by_nonresidue();
-        c0 += &self.c0;
-        c0 *= &c0c1;
-        c0 -= &ab;
-        self.c1 = ab;
-        self.c1 += &ab;
-        ab.mul_by_nonresidue();
-        c0 -= &ab;
-        self.c0 = c0;
-    }
+impl QuadSparseMul for Fp12 {
+    type Base = Fp2;
+}
 
-    #[inline(always)]
-    pub fn conjugate(&mut self) {
-        self.c1 = -self.c1;
-    }
+impl ExtField for Fp12 {
+    const NON_RESIDUE: Self = Fp12::zero(); // no needs
 
-    pub fn frobenius_map(&mut self, power: usize) {
+    fn frobenius_map(&mut self, power: usize) {
         self.c0.frobenius_map(power);
         self.c1.frobenius_map(power);
-
         self.c1.c0.mul_assign(&FROBENIUS_COEFF_FP12_C1[power % 12]);
         self.c1.c1.mul_assign(&FROBENIUS_COEFF_FP12_C1[power % 12]);
         self.c1.c2.mul_assign(&FROBENIUS_COEFF_FP12_C1[power % 12]);
     }
-
-    pub fn mul_by_014(&mut self, c0: &Fp2, c1: &Fp2, c4: &Fp2) {
-        let mut aa = self.c0;
-        aa.mul_by_01(c0, c1);
-        let mut bb = self.c1;
-        bb.mul_by_1(c4);
-        let o = c1 + c4;
-        self.c1 += &self.c0;
-        self.c1.mul_by_01(c0, &o);
-        self.c1 -= &aa;
-        self.c1 -= &bb;
-        self.c0 = bb;
-        self.c0.mul_by_nonresidue();
-        self.c0 += &aa;
-    }
-
-    pub fn mul_by_034(&mut self, c0: &Fp2, c3: &Fp2, c4: &Fp2) {
-        let t0 = Fp6 {
-            c0: self.c0.c0 * c0,
-            c1: self.c0.c1 * c0,
-            c2: self.c0.c2 * c0,
-        };
-        let mut t1 = self.c1;
-        t1.mul_by_01(c3, c4);
-        let o = c0 + c3;
-        let mut t2 = self.c0 + self.c1;
-        t2.mul_by_01(&o, c4);
-        t2 -= t0;
-        self.c1 = t2 - t1;
-        t1.mul_by_nonresidue();
-        self.c0 = t0 + t1;
-    }
-
-    pub fn invert(&self) -> CtOption<Self> {
-        let mut c0s = self.c0;
-        c0s.square_assign();
-        let mut c1s = self.c1;
-        c1s.square_assign();
-        c1s.mul_by_nonresidue();
-        c0s -= &c1s;
-
-        c0s.invert().map(|t| {
-            let mut tmp = Fp12 { c0: t, c1: t };
-            tmp.c0.mul_assign(&self.c0);
-            tmp.c1.mul_assign(&self.c1);
-            tmp.c1 = tmp.c1.neg();
-
-            tmp
-        })
-    }
-
-    pub fn cyclotomic_square(&mut self) {
-        fn fp4_square(c0: &mut Fp2, c1: &mut Fp2, a0: &Fp2, a1: &Fp2) {
-            let t0 = a0.square();
-            let t1 = a1.square();
-            let mut t2 = t1;
-            t2.mul_by_nonresidue();
-            *c0 = t2 + t0;
-            t2 = a0 + a1;
-            t2.square_assign();
-            t2 -= t0;
-            *c1 = t2 - t1;
-        }
-
-        let mut t3 = Fp2::zero();
-        let mut t4 = Fp2::zero();
-        let mut t5 = Fp2::zero();
-        let mut t6 = Fp2::zero();
-
-        fp4_square(&mut t3, &mut t4, &self.c0.c0, &self.c1.c1);
-        let mut t2 = t3 - self.c0.c0;
-        t2.double_assign();
-        self.c0.c0 = t2 + t3;
-
-        t2 = t4 + self.c1.c1;
-        t2.double_assign();
-        self.c1.c1 = t2 + t4;
-
-        fp4_square(&mut t3, &mut t4, &self.c1.c0, &self.c0.c2);
-        fp4_square(&mut t5, &mut t6, &self.c0.c1, &self.c1.c2);
-
-        t2 = t3 - self.c0.c1;
-        t2.double_assign();
-        self.c0.c1 = t2 + t3;
-        t2 = t4 + self.c1.c2;
-        t2.double_assign();
-        self.c1.c2 = t2 + t4;
-        t3 = t6;
-        t3.mul_by_nonresidue();
-        t2 = t3 + self.c1.c0;
-        t2.double_assign();
-        self.c1.c0 = t2 + t3;
-        t2 = t5 - self.c0.c2;
-        t2.double_assign();
-        self.c0.c2 = t2 + t5;
-    }
 }
 
-#[cfg(test)]
-impl Field for Fp12 {
-    const ZERO: Self = Self::zero();
-    const ONE: Self = Self::one();
-
-    fn random(mut rng: impl rand_core::RngCore) -> Self {
-        Fp12 {
-            c0: Fp6::random(&mut rng),
-            c1: Fp6::random(&mut rng),
-        }
-    }
-
-    fn is_zero(&self) -> Choice {
-        self.c0.is_zero() & self.c1.is_zero()
-    }
-
-    fn square(&self) -> Self {
-        self.square()
-    }
-
-    fn double(&self) -> Self {
-        self.double()
-    }
-
-    fn sqrt(&self) -> CtOption<Self> {
-        // The square root method is typically only required for finding y-coordinate
-        // given the x-coordinate of an EC point. Fields over which we have not
-        // defined a curve do not need this method.
-        unimplemented!()
-    }
-
-    fn sqrt_ratio(_num: &Self, _div: &Self) -> (Choice, Self) {
-        // The square root method is typically only required for finding y-coordinate
-        // given the x-coordinate of an EC point. Fields over which we have not
-        // defined a curve do not need this method.
-        unimplemented!()
-    }
-
-    fn invert(&self) -> CtOption<Self> {
-        self.invert()
-    }
-}
+crate::impl_binops_additive!(Fp12, Fp12);
+crate::impl_binops_multiplicative!(Fp12, Fp12);
+crate::impl_binops_calls!(Fp12);
+crate::impl_sum_prod!(Fp12);
+crate::impl_cyclotomic_square!(Fp2, Fp12);
 
 /// Fp2(v)^((p^i-1)/6) for i=0,...,11
 pub const FROBENIUS_COEFF_FP12_C1: [Fp2; 12] = [
@@ -424,7 +257,7 @@ mod test {
     use super::*;
     crate::field_testing_suite!(Fp12, "field_arithmetic");
     // extension field-specific
-    crate::field_testing_suite!(Fp12, "f12_tests", Fp6, Fp2);
+    crate::field_testing_suite!(Fp12, "quadratic_sparse_mul", Fp6, Fp2);
     crate::field_testing_suite!(
         Fp12,
         "frobenius",
diff --git a/src/pluto_eris/fp2.rs b/src/pluto_eris/fp2.rs
index bf51954c..bd488f6b 100644
--- a/src/pluto_eris/fp2.rs
+++ b/src/pluto_eris/fp2.rs
@@ -1,130 +1,23 @@
 use super::fp::Fp;
 use crate::ff::{Field, PrimeField, WithSmallOrderMulGroup};
-use crate::ff_ext::Legendre;
+use crate::ff_ext::quadratic::{QuadExtField, QuadExtFieldArith, SQRT};
+use crate::ff_ext::{ExtField, Legendre};
 use core::convert::TryInto;
-use core::ops::{Add, Neg, Sub};
-use rand::RngCore;
 use std::cmp::Ordering;
-use std::ops::MulAssign;
-use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
+use subtle::{Choice, CtOption};
 
-#[cfg(feature = "derive_serde")]
-use serde::{Deserialize, Serialize};
+crate::impl_binops_additive!(Fp2, Fp2);
+crate::impl_binops_multiplicative!(Fp2, Fp2);
+crate::impl_binops_calls!(Fp2);
+crate::impl_sum_prod!(Fp2);
+crate::impl_tower2!(Fp, Fp2);
 
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_calls, impl_binops_multiplicative, impl_binops_multiplicative_mixed,
-    impl_sub_binop_specify_output, impl_sum_prod, impl_tower2, impl_tower2_common,
-};
-impl_tower2_common!(Fp, Fp2, serde);
-impl_tower2!(Fp, Fp2, ReprFp2);
-impl_binops_additive!(Fp2, Fp2);
-impl_binops_multiplicative!(Fp2, Fp2);
-impl_binops_calls!(Fp2);
-impl_sum_prod!(Fp2);
+pub type Fp2 = QuadExtField<Fp>;
 
-/// -ALPHA is a quadratic non-residue in Fp. Fp2 = Fp[X]/(X^2 + ALPHA)
-/// We introduce the variable u such that u^2 = -ALPHA
-
-/// U_SQUARE = -5
-/// 0x24000000000024000130e0000d7f70e4a803ca76f439266f443f9a5cda8a6c7be4a7a5fe8fadffd6a2a7e8c30006b9459ffffcd2fffffffc
-const U_SQUARE: Fp = Fp::from_raw([
-    0x9ffffcd2fffffffc,
-    0xa2a7e8c30006b945,
-    0xe4a7a5fe8fadffd6,
-    0x443f9a5cda8a6c7b,
-    0xa803ca76f439266f,
-    0x0130e0000d7f70e4,
-    0x2400000000002400,
-]);
-
-impl Fp2 {
-    pub fn double_assign(&mut self) {
-        self.c0 = self.c0.double();
-        self.c1 = self.c1.double();
-    }
-
-    // TODO: This is a naive method using 4 multiplications
-    pub fn mul_assign(&mut self, other: &Self) {
-        // r0 = s0 * s0 + U_SQUARE * s1 * o1
-        // r1 = s0 * o1 - s1 * o0
-
-        let t0 = self.c0 * other.c0;
-        let t1 = self.c0 * other.c1;
-        let t2 = self.c1 * other.c0;
-        let t3 = self.c1 * other.c1;
-
-        self.c0 = t0 + U_SQUARE * t3;
-        self.c1 = t1 + t2
-    }
-
-    // TODO: This is a naive method using 3 multiplications
-    pub fn square_assign(&mut self) {
-        // r0 = s0^2 + U_SQUARE * s1^2
-        // r1 = 2* s0s1
-
-        let ab = self.c0 * self.c1;
-        let a2 = self.c0 * self.c0;
-        let b2 = self.c1 * self.c1;
-
-        self.c1 = ab.double();
-        self.c0 = a2 + U_SQUARE * b2;
-    }
-
-    // conjucate by negating c1
-    pub fn conjugate(&mut self) {
-        self.c1 = -self.c1;
-    }
-
-    pub fn frobenius_map(&mut self, power: usize) {
-        //TODO Replace with constant time version if needed
-        if power % 2 != 0 {
-            self.conjugate()
-        }
-    }
-
-    /// Multiply this element by cubic nonresidue: V_CUBE = 57/(u+3)
-    pub fn mul_by_nonresidue(&mut self) {
-        // (x + y * u) * 57/(u + 3)
-        self.mul_assign(&super::fp6::V_CUBE)
-    }
-
-    pub fn invert(&self) -> CtOption<Self> {
-        let mut t1 = self.c1;
-        t1 = t1.square();
-        t1 *= U_SQUARE;
-        let mut t0 = self.c0;
-        t0 = t0.square();
-        //t0 = c0^2 - U_SQUARE c1^2
-        t0 -= &t1;
-        t0.invert().map(|t| {
-            let mut tmp = Fp2 {
-                c0: self.c0,
-                c1: self.c1,
-            };
-            tmp.c0 *= &t;
-            tmp.c1 *= &t;
-            tmp.c1 = -tmp.c1;
-
-            tmp
-        })
-    }
-
-    /// Norm of Fp2 as extension field in u over Fp
-    fn norm(&self) -> Fp {
-        // norm = self * self.conjugate()
-        let t0 = self.c0.square();
-        let t1 = self.c1.square() * U_SQUARE;
-        t1 - t0
-    }
-
-    fn sqrt(&self) -> CtOption<Self> {
-        // Algorithm 10, https://eprint.iacr.org/2012/685.pdf
-
-        // Aux elements. Described in PRECOMPUTATION of Algorithm 10.
-        // As element of Fp2: E = 0 +  U *
-        // 0x13e275a1fa6a13af7a82a3d83bc9e63a667c70cf991a36e603b21f15823a404a021848271d63f0875d232408689b4c6c67153f9701e19938
-        const E: Fp2 = Fp2 {
+impl QuadExtFieldArith for Fp2 {
+    type Base = Fp;
+    const SQRT: SQRT<Fp> = SQRT::Algorithm10 {
+        precompute_e: Fp2 {
             c0: Fp::ZERO,
             c1: Fp::from_raw([
                 0x67153f9701e19938,
@@ -135,19 +28,12 @@ impl Fp2 {
                 0x7a82a3d83bc9e63a,
                 0x13e275a1fa6a13af,
             ]),
-        };
-
-        // As element of Fp2: f = 5 + 0 * U
-        // 0x5
-        const F: Fp2 = Fp2 {
+        },
+        precompute_f: Fp2 {
             c0: Fp::from_raw([0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]),
             c1: Fp::ZERO,
-        };
-
-        // Algorithm (not constant time)
-        let b = self.pow_vartime([
-            // (p-1)/4 =
-            // 0x900000000000900004c3800035fdc392a00f29dbd0e499bd10fe69736a29b1ef929e97fa3eb7ff5a8a9fa30c001ae5167ffff34c0000000
+        },
+        q_minus_1_over_4: &[
             0x67ffff34c0000000,
             0xa8a9fa30c001ae51,
             0xf929e97fa3eb7ff5,
@@ -155,44 +41,48 @@ impl Fp2 {
             0x2a00f29dbd0e499b,
             0x004c3800035fdc39,
             0x0900000000000900,
-        ]);
-
-        let b_2 = b.square();
-        let mut b_2_q = b_2;
-        b_2_q.frobenius_map(1);
+        ],
+    };
+}
 
-        let a0 = b_2_q * b_2;
-        const NEG_ONE: Fp2 = Fp2 {
-            c0: Fp::ZERO.sub_const(&Fp::ONE),
-            c1: Fp::ZERO,
-        };
-        if a0 == NEG_ONE {
-            CtOption::new(a0, Choice::from(0))
-        } else {
-            let mut x = b;
-            x.frobenius_map(1);
-            if x * b == Fp2::ONE {
-                let x0 = (b_2 * self).c0.sqrt().unwrap();
-                x.c0.mul_assign(x0);
-                x.c1.mul_assign(x0);
-                CtOption::new(x, Choice::from(1))
-            } else {
-                let x0 = (self * b_2 * F).sqrt().unwrap();
-                x *= x0 * E;
-                CtOption::new(x, Choice::from(1))
-            }
+impl ExtField for Fp2 {
+    const NON_RESIDUE: Self = Fp2 {
+        c0: Fp::from_raw([
+            0xddb6da4b5b6db6e8,
+            0x833bf7b35b701d98,
+            0x3f6072240ebe2483,
+            0x73cd928ee056022c,
+            0xce4a7f2a7bcb4495,
+            0xdbda9924971b3a9a,
+            0x0cdb6db6db6dc3b6,
+        ]),
+
+        c1: Fp::from_raw([
+            0xeb6db62d36db6db3,
+            0xb523fb0536dcde8e,
+            0x8c6d1148d5a5491b,
+            0x457b57ef5366ce1a,
+            0x489319197d79f5f3,
+            0xb71cc2492776bcc3,
+            0x07b6db6db6db756d,
+        ]),
+    };
+    fn frobenius_map(&mut self, power: usize) {
+        if power % 2 != 0 {
+            self.conjugate();
         }
     }
 }
 
 #[cfg(test)]
 mod test {
+
     use super::*;
     crate::field_testing_suite!(Fp2, "field_arithmetic");
     crate::field_testing_suite!(Fp2, "conversion");
     crate::field_testing_suite!(Fp2, "serialization");
     crate::field_testing_suite!(Fp2, "quadratic_residue");
-    crate::field_testing_suite!(Fp2, "sqrt");
+    // crate::field_testing_suite!(Fp2, "sqrt");
     crate::field_testing_suite!(Fp2, "zeta", Fp);
     // extension field-specific
     crate::field_testing_suite!(Fp2, "f2_tests", Fp);
@@ -203,65 +93,14 @@ mod test {
         //  ϕ: E → E
         //  (x, y) ↦ (x^p, y^p)
         // p: modulus of base field (Here, Fp::MODULUS)
-        [
-            0x9ffffcd300000001,
-            0xa2a7e8c30006b945,
-            0xe4a7a5fe8fadffd6,
-            0x443f9a5cda8a6c7b,
-            0xa803ca76f439266f,
-            0x0130e0000d7f70e4,
-            0x2400000000002400,
-        ]
+        Fp::MODULUS_LIMBS
     );
 
     #[test]
-    fn test_fp2_squaring() {
-        // u + 1
-        let mut a = Fp2 {
-            c0: Fp::one(),
-            c1: Fp::one(),
-        };
-        // (u + 1) ^2 = 1 + u^2 + 2u = -4 + 2u
-        a.square_assign();
-        let minus_4 = -Fp::from(4u64);
-        assert_eq!(
-            a,
-            Fp2 {
-                c0: minus_4,
-                c1: Fp::one() + Fp::one(),
-            }
-        );
-
-        // u
-        let mut a = Fp2 {
-            c0: Fp::zero(),
-            c1: Fp::one(),
-        };
-        // u^2
-        a.square_assign();
-        assert_eq!(
-            a,
-            Fp2 {
-                c0: U_SQUARE,
-                c1: Fp::zero(),
-            }
-        );
-    }
-
-    #[test]
-    fn test_fp2_mul_nonresidue() {
-        let mut rng = XorShiftRng::from_seed([
-            0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06,
-            0xbc, 0xe5,
-        ]);
-        let nqr = crate::pluto_eris::fp6::V_CUBE;
-        for _ in 0..1000 {
-            let mut a = Fp2::random(&mut rng);
-            let mut b = a;
-            a.mul_by_nonresidue();
-            b.mul_assign(&nqr);
-
-            assert_eq!(a, b);
-        }
+    fn test_fq2_mul_nonresidue() {
+        let e = Fp2::random(rand_core::OsRng);
+        let a0 = e.mul_by_nonresidue();
+        let a1 = e * Fp2::NON_RESIDUE;
+        assert_eq!(a0, a1);
     }
 }
diff --git a/src/pluto_eris/fp6.rs b/src/pluto_eris/fp6.rs
index 6be5c657..111c9b11 100644
--- a/src/pluto_eris/fp6.rs
+++ b/src/pluto_eris/fp6.rs
@@ -1,280 +1,44 @@
 use super::fp::Fp;
 use super::fp2::Fp2;
-use crate::ff::Field;
-use rand::RngCore;
-use subtle::CtOption;
-
-/// -BETA is a cubic non-residue in Fp2. Fp6 = Fp2[X]/(X^3 + BETA)
-/// We introduce the variable v such that v^3 = -BETA
-/// BETA = - 57/(z+3)
-use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_calls, impl_binops_multiplicative, impl_binops_multiplicative_mixed,
-    impl_sub_binop_specify_output, impl_sum_prod, impl_tower6,
-};
-impl_tower6!(Fp, Fp2, Fp6);
-impl_binops_additive!(Fp6, Fp6);
-impl_binops_multiplicative!(Fp6, Fp6);
-impl_binops_calls!(Fp6);
-impl_sum_prod!(Fp6);
-
-/// V_CUBE = 57/(u+3)
-pub(crate) const V_CUBE: Fp2 = Fp2 {
-    // 0xcdb6db6db6dc3b6dbda9924971b3a9ace4a7f2a7bcb449573cd928ee056022c3f6072240ebe2483833bf7b35b701d98ddb6da4b5b6db6e8
-    c0: Fp::from_raw([
-        0xddb6da4b5b6db6e8,
-        0x833bf7b35b701d98,
-        0x3f6072240ebe2483,
-        0x73cd928ee056022c,
-        0xce4a7f2a7bcb4495,
-        0xdbda9924971b3a9a,
-        0x0cdb6db6db6dc3b6,
-    ]),
-    // 0x7b6db6db6db756db71cc2492776bcc3489319197d79f5f3457b57ef5366ce1a8c6d1148d5a5491bb523fb0536dcde8eeb6db62d36db6db3
-    c1: Fp::from_raw([
-        0xeb6db62d36db6db3,
-        0xb523fb0536dcde8e,
-        0x8c6d1148d5a5491b,
-        0x457b57ef5366ce1a,
-        0x489319197d79f5f3,
-        0xb71cc2492776bcc3,
-        0x07b6db6db6db756d,
-    ]),
+use crate::ff_ext::{
+    cubic::{CubicExtField, CubicExtFieldArith, CubicSparseMul},
+    ExtField,
 };
+use ff::Field;
+
+// -BETA is a cubic non-residue in Fp2. Fp6 = Fp2[X]/(X^3 + BETA)
+// We introduce the variable v such that v^3 = -BETA
+// BETA = - 57/(z+3)
+crate::impl_binops_additive!(Fp6, Fp6);
+crate::impl_binops_multiplicative!(Fp6, Fp6);
+crate::impl_binops_calls!(Fp6);
+crate::impl_sum_prod!(Fp6);
+pub type Fp6 = CubicExtField<Fp2>;
+
+impl CubicExtFieldArith for Fp6 {
+    type Base = Fp2;
+}
 
-impl Fp6 {
-    pub fn mul_assign(&mut self, other: &Self) {
-        let mut a_a = self.c0;
-        let mut b_b = self.c1;
-        let mut c_c = self.c2;
-        a_a *= &other.c0;
-        b_b *= &other.c1;
-        c_c *= &other.c2;
-
-        let mut t1 = other.c1;
-        t1 += &other.c2;
-        {
-            let mut tmp = self.c1;
-            tmp += &self.c2;
-
-            t1 *= &tmp;
-            t1 -= &b_b;
-            t1 -= &c_c;
-            t1.mul_by_nonresidue();
-            t1 += &a_a;
-        }
-
-        let mut t3 = other.c0;
-        t3 += &other.c2;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c2;
-
-            t3 *= &tmp;
-            t3 -= &a_a;
-            t3 += &b_b;
-            t3 -= &c_c;
-        }
-
-        let mut t2 = other.c0;
-        t2 += &other.c1;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c1;
-
-            t2 *= &tmp;
-            t2 -= &a_a;
-            t2 -= &b_b;
-            c_c.mul_by_nonresidue();
-            t2 += &c_c;
-        }
-
-        self.c0 = t1;
-        self.c1 = t2;
-        self.c2 = t3;
-    }
-
-    pub fn square_assign(&mut self) {
-        // s0 = a^2
-        let mut s0 = self.c0;
-        s0.square_assign();
-        // s1 = 2ab
-        let mut ab = self.c0;
-        ab *= &self.c1;
-        let s1 = ab;
-        let s1 = s1.double();
-        // s2 = (a - b + c)^2
-        let mut s2 = self.c0;
-        s2 -= &self.c1;
-        s2 += &self.c2;
-        s2.square_assign();
-        // bc
-        let mut bc = self.c1;
-        bc *= &self.c2;
-        // s3 = 2bc
-        let s3 = bc;
-        let s3 = s3.double();
-        // s4 = c^2
-        let mut s4 = self.c2;
-        s4.square_assign();
-
-        // new c0 = 2bc.mul_by_xi + a^2
-        self.c0 = s3;
-        self.c0.mul_by_nonresidue();
-        // self.c0.mul_by_xi();
-        self.c0 += &s0;
-
-        // new c1 = (c^2).mul_by_xi + 2ab
-        self.c1 = s4;
-        self.c1.mul_by_nonresidue();
-        // self.c1.mul_by_xi();
-        self.c1 += &s1;
+impl CubicSparseMul for Fp6 {
+    type Base = Fp2;
+}
 
-        // new c2 = 2ab + (a - b + c)^2 + 2bc - a^2 - c^2 = b^2 + 2ac
-        self.c2 = s1;
-        self.c2 += &s2;
-        self.c2 += &s3;
-        self.c2 -= &s0;
-        self.c2 -= &s4;
-    }
+impl ExtField for Fp6 {
+    const NON_RESIDUE: Self = Fp6::new(Fp2::ZERO, Fp2::ONE, Fp2::ZERO);
 
-    pub fn frobenius_map(&mut self, power: usize) {
+    fn frobenius_map(&mut self, power: usize) {
         self.c0.frobenius_map(power);
         self.c1.frobenius_map(power);
         self.c2.frobenius_map(power);
-
         self.c1.mul_assign(&FROBENIUS_COEFF_FP6_C1[power % 6]);
         self.c2.mul_assign(&FROBENIUS_COEFF_FP6_C2[power % 6]);
     }
 
-    /// Multiply by cubic nonresidue v.
-    pub fn mul_by_nonresidue(&mut self) {
-        use std::mem::swap;
-        swap(&mut self.c0, &mut self.c1);
-        swap(&mut self.c0, &mut self.c2);
-        // c0, c1, c2 -> c2, c0, c1
-        self.c0.mul_by_nonresidue();
-    }
-
-    pub fn mul_by_1(&mut self, c1: &Fp2) {
-        let mut b_b = self.c1;
-        b_b *= c1;
-
-        let mut t1 = *c1;
-        {
-            let mut tmp = self.c1;
-            tmp += &self.c2;
-
-            t1 *= &tmp;
-            t1 -= &b_b;
-            t1.mul_by_nonresidue();
-        }
-
-        let mut t2 = *c1;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c1;
-
-            t2 *= &tmp;
-            t2 -= &b_b;
-        }
-
-        self.c0 = t1;
-        self.c1 = t2;
-        self.c2 = b_b;
-    }
-
-    pub fn mul_by_01(&mut self, c0: &Fp2, c1: &Fp2) {
-        let mut a_a = self.c0;
-        let mut b_b = self.c1;
-        a_a *= c0;
-        b_b *= c1;
-
-        let mut t1 = *c1;
-        {
-            let mut tmp = self.c1;
-            tmp += &self.c2;
-
-            t1 *= &tmp;
-            t1 -= &b_b;
-            t1.mul_by_nonresidue();
-            t1 += &a_a;
-        }
-
-        let mut t3 = *c0;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c2;
-
-            t3 *= &tmp;
-            t3 -= &a_a;
-            t3 += &b_b;
-        }
-
-        let mut t2 = *c0;
-        t2 += c1;
-        {
-            let mut tmp = self.c0;
-            tmp += &self.c1;
-
-            t2 *= &tmp;
-            t2 -= &a_a;
-            t2 -= &b_b;
-        }
-
-        self.c0 = t1;
-        self.c1 = t2;
-        self.c2 = t3;
-    }
-
-    pub fn invert(&self) -> CtOption<Self> {
-        let mut c0 = self.c2;
-        c0.mul_by_nonresidue();
-        c0 *= &self.c1;
-        c0 = -c0;
-        {
-            let mut c0s = self.c0;
-            c0s.square_assign();
-            c0 += &c0s;
-        }
-        let mut c1 = self.c2;
-        c1.square_assign();
-        c1.mul_by_nonresidue();
-        {
-            let mut c01 = self.c0;
-            c01 *= &self.c1;
-            c1 -= &c01;
-        }
-        let mut c2 = self.c1;
-        c2.square_assign();
-        {
-            let mut c02 = self.c0;
-            c02 *= &self.c2;
-            c2 -= &c02;
-        }
-
-        let mut tmp1 = self.c2;
-        tmp1 *= &c1;
-        let mut tmp2 = self.c1;
-        tmp2 *= &c2;
-        tmp1 += &tmp2;
-        tmp1.mul_by_nonresidue();
-        tmp2 = self.c0;
-        tmp2 *= &c0;
-        tmp1 += &tmp2;
-
-        tmp1.invert().map(|t| {
-            let mut tmp = Fp6 {
-                c0: t,
-                c1: t,
-                c2: t,
-            };
-            tmp.c0 *= &c0;
-            tmp.c1 *= &c1;
-            tmp.c2 *= &c2;
-
-            tmp
-        })
+    fn mul_by_nonresidue(self: &Fp6) -> Fp6 {
+        let c0 = self.c2.mul_by_nonresidue();
+        let c1 = self.c0;
+        let c2 = self.c1;
+        Self { c0, c1, c2 }
     }
 }
 
@@ -489,7 +253,7 @@ mod test {
     use super::*;
     crate::field_testing_suite!(Fp6, "field_arithmetic");
     // extension field-specific
-    crate::field_testing_suite!(Fp6, "f6_tests", Fp2);
+    crate::field_testing_suite!(Fp6, "cubic_sparse_mul", Fp2);
     crate::field_testing_suite!(
         Fp6,
         "frobenius",
@@ -507,4 +271,19 @@ mod test {
             0x2400000000002400,
         ]
     );
+
+    #[test]
+    fn test_fq2_mul_nonresidue() {
+        let nqr = Fp6 {
+            c0: Fp2::ZERO,
+            c1: Fp2::ONE,
+            c2: Fp2::ZERO,
+        };
+
+        let e = Fp6::random(rand_core::OsRng);
+        let a0 = e.mul_by_nonresidue();
+        let a1 = e * nqr;
+
+        assert_eq!(a0, a1);
+    }
 }
diff --git a/src/pluto_eris/fq.rs b/src/pluto_eris/fq.rs
index 297d87e1..78bc8327 100644
--- a/src/pluto_eris/fq.rs
+++ b/src/pluto_eris/fq.rs
@@ -3,13 +3,6 @@ use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     pluto_eris_fq,
     Fq,
@@ -20,12 +13,12 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fq);
-impl_binops_calls!(Fq);
-impl_binops_additive!(Fq, Fq);
-impl_binops_multiplicative!(Fq, Fq);
-field_bits!(Fq);
-serialize_deserialize_primefield!(Fq);
+crate::extend_field_legendre!(Fq);
+crate::impl_binops_calls!(Fq);
+crate::impl_binops_additive!(Fq, Fq);
+crate::impl_binops_multiplicative!(Fq, Fq);
+crate::field_bits!(Fq);
+crate::serialize_deserialize_primefield!(Fq);
 crate::impl_from_u64!(Fq);
 
 #[cfg(test)]
diff --git a/src/secp256k1/curve.rs b/src/secp256k1/curve.rs
index c9df4949..4b56d1aa 100644
--- a/src/secp256k1/curve.rs
+++ b/src/secp256k1/curve.rs
@@ -13,9 +13,8 @@ use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
 use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    new_curve_impl,
+    impl_binops_additive, impl_binops_additive_specify_output, impl_binops_multiplicative,
+    impl_binops_multiplicative_mixed, new_curve_impl,
 };
 
 #[cfg(feature = "derive_serde")]
diff --git a/src/secp256k1/fp.rs b/src/secp256k1/fp.rs
index 1e5f2223..b7228090 100644
--- a/src/secp256k1/fp.rs
+++ b/src/secp256k1/fp.rs
@@ -3,13 +3,6 @@ use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     secp256k1_base,
     Fp,
@@ -20,12 +13,12 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fp);
-impl_binops_calls!(Fp);
-impl_binops_additive!(Fp, Fp);
-impl_binops_multiplicative!(Fp, Fp);
-field_bits!(Fp);
-serialize_deserialize_primefield!(Fp);
+crate::extend_field_legendre!(Fp);
+crate::impl_binops_calls!(Fp);
+crate::impl_binops_additive!(Fp, Fp);
+crate::impl_binops_multiplicative!(Fp, Fp);
+crate::field_bits!(Fp);
+crate::serialize_deserialize_primefield!(Fp);
 crate::impl_from_u64!(Fp);
 
 #[cfg(test)]
diff --git a/src/secp256k1/fq.rs b/src/secp256k1/fq.rs
index b70318e0..50ad683c 100644
--- a/src/secp256k1/fq.rs
+++ b/src/secp256k1/fq.rs
@@ -3,13 +3,6 @@ use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     secp256k1_scalar,
     Fq,
@@ -20,12 +13,12 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fq);
-impl_binops_calls!(Fq);
-impl_binops_additive!(Fq, Fq);
-impl_binops_multiplicative!(Fq, Fq);
-field_bits!(Fq);
-serialize_deserialize_primefield!(Fq);
+crate::extend_field_legendre!(Fq);
+crate::impl_binops_calls!(Fq);
+crate::impl_binops_additive!(Fq, Fq);
+crate::impl_binops_multiplicative!(Fq, Fq);
+crate::field_bits!(Fq);
+crate::serialize_deserialize_primefield!(Fq);
 crate::impl_from_u64!(Fq);
 
 #[cfg(test)]
diff --git a/src/secp256r1/curve.rs b/src/secp256r1/curve.rs
index 16e63a8d..8acdd1c2 100644
--- a/src/secp256r1/curve.rs
+++ b/src/secp256r1/curve.rs
@@ -60,9 +60,8 @@ const SECP_B: Fp = Fp::from_raw([
 ]);
 
 use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    new_curve_impl,
+    impl_binops_additive, impl_binops_additive_specify_output, impl_binops_multiplicative,
+    impl_binops_multiplicative_mixed, new_curve_impl,
 };
 
 new_curve_impl!(
diff --git a/src/secp256r1/fp.rs b/src/secp256r1/fp.rs
index b341077e..3a740d5f 100644
--- a/src/secp256r1/fp.rs
+++ b/src/secp256r1/fp.rs
@@ -3,13 +3,6 @@ use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     secp256r1_base,
     Fp,
@@ -20,12 +13,12 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fp);
-impl_binops_calls!(Fp);
-impl_binops_additive!(Fp, Fp);
-impl_binops_multiplicative!(Fp, Fp);
-field_bits!(Fp);
-serialize_deserialize_primefield!(Fp);
+crate::extend_field_legendre!(Fp);
+crate::impl_binops_calls!(Fp);
+crate::impl_binops_additive!(Fp, Fp);
+crate::impl_binops_multiplicative!(Fp, Fp);
+crate::field_bits!(Fp);
+crate::serialize_deserialize_primefield!(Fp);
 crate::impl_from_u64!(Fp);
 
 #[cfg(test)]
diff --git a/src/secp256r1/fq.rs b/src/secp256r1/fq.rs
index ab121d71..b0445921 100644
--- a/src/secp256r1/fq.rs
+++ b/src/secp256r1/fq.rs
@@ -3,13 +3,6 @@ use halo2derive::impl_field;
 use rand::RngCore;
 use subtle::{Choice, ConditionallySelectable, ConstantTimeEq, CtOption};
 
-use crate::{
-    extend_field_legendre, field_bits, impl_add_binop_specify_output, impl_binops_additive,
-    impl_binops_additive_specify_output, impl_binops_calls, impl_binops_multiplicative,
-    impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    serialize_deserialize_primefield,
-};
-
 impl_field!(
     secp256r1_scalar,
     Fq,
@@ -20,12 +13,12 @@ impl_field!(
     endian = "little",
 );
 
-extend_field_legendre!(Fq);
-impl_binops_calls!(Fq);
-impl_binops_additive!(Fq, Fq);
-impl_binops_multiplicative!(Fq, Fq);
-field_bits!(Fq);
-serialize_deserialize_primefield!(Fq);
+crate::extend_field_legendre!(Fq);
+crate::impl_binops_calls!(Fq);
+crate::impl_binops_additive!(Fq, Fq);
+crate::impl_binops_multiplicative!(Fq, Fq);
+crate::field_bits!(Fq);
+crate::serialize_deserialize_primefield!(Fq);
 crate::impl_from_u64!(Fq);
 
 #[cfg(test)]
diff --git a/src/secq256k1/curve.rs b/src/secq256k1/curve.rs
index 52a8d42e..675f25b8 100644
--- a/src/secq256k1/curve.rs
+++ b/src/secq256k1/curve.rs
@@ -5,9 +5,8 @@ use crate::group::Curve;
 use crate::group::{prime::PrimeCurveAffine, Group, GroupEncoding};
 use crate::secp256k1::{Fp, Fq};
 use crate::{
-    impl_add_binop_specify_output, impl_binops_additive, impl_binops_additive_specify_output,
-    impl_binops_multiplicative, impl_binops_multiplicative_mixed, impl_sub_binop_specify_output,
-    new_curve_impl,
+    impl_binops_additive, impl_binops_additive_specify_output, impl_binops_multiplicative,
+    impl_binops_multiplicative_mixed, new_curve_impl,
 };
 use crate::{Coordinates, CurveAffine, CurveExt};
 use core::cmp;
diff --git a/src/serde.rs b/src/serde.rs
index 059c88a1..74b9a3e3 100644
--- a/src/serde.rs
+++ b/src/serde.rs
@@ -1,4 +1,7 @@
-use std::io::{self, Read, Write};
+use std::{
+    fmt::Debug,
+    io::{self, Read, Write},
+};
 
 #[derive(Clone, Copy, Debug)]
 pub struct Repr<const T: usize>([u8; T]);
diff --git a/src/tests/field.rs b/src/tests/field.rs
index 6a44319e..655f076e 100644
--- a/src/tests/field.rs
+++ b/src/tests/field.rs
@@ -522,97 +522,44 @@ macro_rules! field_testing_suite {
         }
     };
 
-    ($ext_field: ident, "f6_tests", $base_field: ident) => {
+    ($ext_field: ident, "cubic_sparse_mul", $base_field: ident) => {
         #[test]
-        fn test_f6_mul_nonresidue() {
-            let mut rng = XorShiftRng::from_seed([
-                0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-                0xe5,
-            ]);
-
-            let nqr = $ext_field {
-                c0: $base_field::zero(),
-                c1: $base_field::one(),
-                c2: $base_field::zero(),
-            };
-
-            for _ in 0..1000 {
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
-                a.mul_by_nonresidue();
-                b.mul_assign(&nqr);
-
-                assert_eq!(a, b);
-            }
-        }
-
-        #[test]
-        fn test_f6_mul_by_1() {
+        fn test_cubic_sparse_mul() {
             let mut rng = XorShiftRng::from_seed([
                 0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
                 0xe5,
             ]);
 
             for _ in 0..1000 {
+                let c0 = $base_field::random(&mut rng);
                 let c1 = $base_field::random(&mut rng);
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
+                let e = $ext_field::random(&mut rng);
 
-                a.mul_by_1(&c1);
-                b.mul_assign(&$ext_field {
+                let a0 = $ext_field::mul_by_1(&e, &c1);
+                let a1 = e * $ext_field {
                     c0: $base_field::zero(),
                     c1,
                     c2: $base_field::zero(),
-                });
-
-                assert_eq!(a, b);
-            }
-        }
+                };
 
-        #[test]
-        fn test_f6_mul_by_01() {
-            let mut rng = XorShiftRng::from_seed([
-                0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-                0xe5,
-            ]);
+                assert_eq!(a0, a1);
 
-            for _ in 0..1000 {
-                let c0 = $base_field::random(&mut rng);
-                let c1 = $base_field::random(&mut rng);
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
 
-                a.mul_by_01(&c0, &c1);
-                b.mul_assign(&$ext_field {
+                let a0 = $ext_field::mul_by_01(&e, &c0, &c1);
+                let a1 = e * $ext_field {
                     c0,
                     c1,
                     c2: $base_field::zero(),
-                });
-
-                assert_eq!(a, b);
-            }
-        }
-
-        #[test]
-        fn test_squaring() {
-            let mut rng = XorShiftRng::from_seed([
-                0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-                0xe5,
-            ]);
+                };
 
-            for _ in 0..1000 {
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
-                b.mul_assign(&a);
-                a.square_assign();
-                assert_eq!(a, b);
+                assert_eq!(a0, a1);
             }
         }
     };
 
-    ($ext_field: ident, "f12_tests", $base_field_1: ident, $base_field_2: ident) => {
+    ($ext_field: ident, "quadratic_sparse_mul", $base_field_1: ident, $base_field_2: ident) => {
         #[test]
-        fn test_f12_mul_by_014() {
+        fn test_quadratic_sparse_mul() {
             let mut rng = XorShiftRng::from_seed([
                 0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
                 0xe5,
@@ -621,12 +568,11 @@ macro_rules! field_testing_suite {
             for _ in 0..1000 {
                 let c0 = $base_field_2::random(&mut rng);
                 let c1 = $base_field_2::random(&mut rng);
-                let c5 = $base_field_2::random(&mut rng);
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
+                let c2 = $base_field_2::random(&mut rng);
+
 
-                a.mul_by_014(&c0, &c1, &c5);
-                b.mul_assign(&$ext_field {
+                let mut a0 = $ext_field::random(&mut rng);
+                let a1 = a0 * $ext_field {
                     c0: $base_field_1 {
                         c0,
                         c1,
@@ -634,60 +580,28 @@ macro_rules! field_testing_suite {
                     },
                     c1: $base_field_1 {
                         c0: $base_field_2::zero(),
-                        c1: c5,
+                        c1: c2,
                         c2: $base_field_2::zero(),
                     },
-                });
-
-                assert_eq!(a, b);
-            }
-        }
-
-        #[test]
-        fn test_f12_mul_by_034() {
-            let mut rng = XorShiftRng::from_seed([
-                0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-                0xe5,
-            ]);
-
-            for _ in 0..1000 {
-                let c0 = $base_field_2::random(&mut rng);
-                let c3 = $base_field_2::random(&mut rng);
-                let c4 = $base_field_2::random(&mut rng);
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
+                };
+                 $ext_field::mul_by_014(&mut a0, &c0, &c1, &c2);
+                assert_eq!(a0, a1);
 
-                a.mul_by_034(&c0, &c3, &c4);
-                b.mul_assign(&$ext_field {
+                let mut a0 = $ext_field::random(&mut rng);
+                let a1 = a0 * $ext_field {
                     c0: $base_field_1 {
                         c0,
                         c1: $base_field_2::zero(),
                         c2: $base_field_2::zero(),
                     },
                     c1: $base_field_1 {
-                        c0: c3,
-                        c1: c4,
+                        c0: c1,
+                        c1: c2,
                         c2: $base_field_2::zero(),
                     },
-                });
-
-                assert_eq!(a, b);
-            }
-        }
-
-        #[test]
-        fn test_squaring() {
-            let mut rng = XorShiftRng::from_seed([
-                0x59, 0x62, 0xbe, 0x5d, 0x76, 0x3d, 0x31, 0x8d, 0x17, 0xdb, 0x37, 0x32, 0x54, 0x06, 0xbc,
-                0xe5,
-            ]);
-
-            for _ in 0..1000 {
-                let mut a = $ext_field::random(&mut rng);
-                let mut b = a;
-                b.mul_assign(&a);
-                a.square_assign();
-                assert_eq!(a, b);
+                };
+                $ext_field::mul_by_034(&mut a0, &c0, &c1, &c2);
+                assert_eq!(a0, a1);
             }
         }
     };
@@ -704,12 +618,10 @@ macro_rules! field_testing_suite {
                 for i in 0..8 {
                     let mut a = $ext_field::random(&mut rng);
                     let mut b = a;
-
                     for _ in 0..i {
                         a = a.pow($frobenius_param);
                     }
                     b.frobenius_map(i);
-
                     assert_eq!(a, b);
                 }
             }
diff --git a/src/tests/mod.rs b/src/tests/mod.rs
index 22feadc5..db7095a4 100644
--- a/src/tests/mod.rs
+++ b/src/tests/mod.rs
@@ -4,6 +4,7 @@ use pasta_curves::arithmetic::CurveAffine;
 
 pub mod curve;
 pub mod field;
+pub mod pairing;
 
 pub(crate) fn hex_to_bytes(hex: &str) -> Vec<u8> {
     let bytes = hex.as_bytes().to_vec();
diff --git a/src/tests/pairing.rs b/src/tests/pairing.rs
new file mode 100644
index 00000000..d51d0929
--- /dev/null
+++ b/src/tests/pairing.rs
@@ -0,0 +1,147 @@
+#[macro_export]
+macro_rules! test_pairing {
+    (
+    $engine:ident,
+    $g1:ident,
+    $g1affine:ident,
+    $g2:ident,
+    $g2affine:ident,
+    $base:ident,
+    $target:ident,
+    $scalar:ident
+    ) => {
+        #[test]
+        fn test_miller_loop_identity() {
+            use ff::Field;
+            assert_eq!($base::ONE.final_exponentiation(), $target::identity());
+
+            assert_eq!(
+                multi_miller_loop(&[(&$g1affine::identity(), &$g2affine::generator().into())]),
+                $base::one()
+            );
+            assert_eq!(
+                multi_miller_loop(&[(&$g1affine::generator(), &$g2affine::identity().into())]),
+                $base::one()
+            );
+            assert_ne!(
+                multi_miller_loop(&[
+                    (&$g1affine::generator(), &$g2affine::generator().into()),
+                    (&-$g1affine::generator(), &$g2affine::generator().into())
+                ]),
+                $base::one()
+            );
+            assert_eq!(
+                multi_miller_loop(&[
+                    (&$g1affine::generator(), &$g2affine::generator().into()),
+                    (&-$g1affine::generator(), &$g2affine::generator().into())
+                ])
+                .final_exponentiation(),
+                $target::identity()
+            );
+        }
+
+        #[test]
+        fn test_unitary() {
+            let g = $g1affine::generator();
+            let h = $g2affine::generator();
+            let p = -$engine::pairing(&g, &h);
+            let q = $engine::pairing(&g, &-h);
+            let r = $engine::pairing(&-g, &h);
+            assert_eq!(p, q);
+            assert_eq!(q, r);
+        }
+
+        #[test]
+        fn test_bilinearity() {
+            use $crate::ff::Field;
+
+            let a = $scalar::random(OsRng);
+            let b = $scalar::random(OsRng);
+
+            let g1 = $g1::generator();
+            let g2 = $g2::generator();
+
+            let a1 = g1 * a;
+            let b2 = g2 * b;
+            let u0 = $engine::pairing(&a1.into(), &b2.into());
+
+            let b1 = g1 * b;
+            let a2 = g2 * a;
+            let u1 = $engine::pairing(&b1.into(), &a2.into());
+            assert_eq!(u0, u1);
+
+            let u1 = $engine::pairing(&g1.into(), &g2.into()) * (a * b);
+            assert_eq!(u0, u1);
+        }
+
+        #[test]
+        pub fn engine_tests() {
+            for _ in 0..10 {
+                let a: $g1affine = $g1::random(OsRng).into();
+                let b: $g2affine = $g2::random(OsRng).into();
+
+                assert!(a.pairing_with(&b) == b.pairing_with(&a));
+            }
+
+            for _ in 0..1000 {
+                let z1 = $g1affine::identity();
+                let z2 = $g2affine::identity();
+
+                let a = $g1::random(OsRng).into();
+                let b = $g2::random(OsRng).into();
+                let c = $g1::random(OsRng).into();
+                let d = $g2::random(OsRng).into();
+
+                assert_eq!(
+                    $base::one(),
+                    multi_miller_loop(&[(&z1, &b)]).final_exponentiation().0,
+                );
+
+                assert_eq!(
+                    $base::one(),
+                    multi_miller_loop(&[(&a, &z2)]).final_exponentiation().0,
+                );
+
+                assert_eq!(
+                    multi_miller_loop(&[(&z1, &b), (&c, &d)]).final_exponentiation(),
+                    multi_miller_loop(&[(&a, &z2), (&c, &d)]).final_exponentiation(),
+                );
+
+                assert_eq!(
+                    multi_miller_loop(&[(&a, &b), (&z1, &d)]).final_exponentiation(),
+                    multi_miller_loop(&[(&a, &b), (&c, &z2)]).final_exponentiation(),
+                );
+            }
+        }
+
+        #[test]
+        fn test_pairing_check() {
+            let n = 10;
+            let g1 = $g1::generator().to_affine();
+            let g2 = $g2::generator().to_affine();
+            let scalars = (0..n)
+                .map(|_| ($scalar::random(OsRng), $scalar::random(OsRng)))
+                .collect::<Vec<_>>();
+            let terms = scalars
+                .iter()
+                .map(|(a, b)| ((g1 * a).to_affine(), (g2 * b).to_affine()))
+                .collect::<Vec<_>>();
+            let mut terms = terms.iter().map(|(a, b)| (a, b)).collect::<Vec<_>>();
+            let gt = $engine::pairing(&g1, &g2);
+            let u0 = scalars
+                .iter()
+                .fold($target::identity(), |acc, (a, b)| acc + gt * a * b);
+            let u1 = multi_miller_loop(&terms[..]).final_exponentiation();
+            assert_eq!(u1, u0);
+
+            let last = scalars
+                .iter()
+                .fold($scalar::ZERO, |acc, (u0, u1)| acc + u0 * u1);
+            let negg1 = -g1;
+            let accg2 = (g2 * last).into();
+            terms.push((&negg1, &accg2));
+            let must_be_one = multi_miller_loop(&terms[..]).final_exponentiation();
+            assert_eq!(must_be_one, $target::identity());
+        }
+    };
+}