diff --git a/src/groth/domain.rs b/src/groth/domain.rs
new file mode 100644
index 0000000..e48ea59
--- /dev/null
+++ b/src/groth/domain.rs
@@ -0,0 +1,192 @@
+use curves::{Engine, Field, SnarkField, PrimeField};
+
+pub struct EvaluationDomain<E: Engine> {
+    pub m: u64,
+    exp: u64,
+    omega: E::Fr,
+    omegainv: E::Fr,
+    geninv: E::Fr,
+    minv: E::Fr
+}
+
+impl<E: Engine> EvaluationDomain<E> {
+    pub fn new(e: &E, needed: u64) -> Self {
+        if needed > 268435456 {
+            panic!("circuit depths larger than 2^28 are not supported");
+        }
+
+        let mut m = 1;
+        let mut exp = 0;
+        while m < needed {
+            m *= 2;
+            exp += 1;
+
+            assert!(exp < E::Fr::s(e));
+        }
+
+        let mut omega = E::Fr::root_of_unity(e);
+        for _ in exp..E::Fr::s(e) {
+            omega.square(e);
+        }
+
+        EvaluationDomain {
+            m: m,
+            exp: exp,
+            omega: omega,
+            omegainv: omega.inverse(e).unwrap(),
+            geninv: E::Fr::multiplicative_generator(e).inverse(e).unwrap(),
+            minv: E::Fr::from_u64(e, m).inverse(e).unwrap()
+        }
+    }
+
+    pub fn z(&self, e: &E, tau: &E::Fr) -> E::Fr {
+        let mut tmp = tau.pow(e, &[self.m]);
+        tmp.sub_assign(e, &E::Fr::one(e));
+
+        tmp
+    }
+
+    pub fn ifft(&self, e: &E, v: &mut [E::Fr])
+    {
+        assert!(v.len() == self.m as usize);
+        self._fft(e, v, &self.omegainv);
+        for v in v {
+            v.mul_assign(e, &self.minv);
+        }
+    }
+
+    fn mul_coset(&self, e: &E, v: &mut [E::Fr], g: &E::Fr)
+    {
+        let mut u = *g;
+        for v in v.iter_mut().skip(1) {
+            v.mul_assign(e, &u);
+            u.mul_assign(e, g);
+        }
+    }
+
+    pub fn coset_fft(&self, e: &E, v: &mut [E::Fr])
+    {
+        self.mul_coset(e, v, &E::Fr::multiplicative_generator(e));
+        self.fft(e, v);
+    }
+
+    pub fn icoset_fft(&self, e: &E, v: &mut [E::Fr])
+    {
+        self.ifft(e, v);
+        self.mul_coset(e, v, &self.geninv);
+    }
+
+    pub fn divide_by_z_on_coset(&self, e: &E, v: &mut [E::Fr])
+    {
+        let i = self.z(e, &E::Fr::multiplicative_generator(e)).inverse(e).unwrap();
+        for v in v {
+            v.mul_assign(e, &i);
+        }
+    }
+
+    pub fn fft(&self, e: &E, a: &mut [E::Fr])
+    {
+        self._fft(e, a, &self.omega);
+    }
+
+    fn _fft(&self, e: &E, a: &mut [E::Fr], omega: &E::Fr)
+    {
+        fn bitreverse(mut n: usize, l: u64) -> usize {
+            let mut r = 0;
+            for _ in 0..l {
+                r = (r << 1) | (n & 1);
+                n >>= 1;
+            }
+            r
+        }
+
+        for k in 0..(self.m as usize) {
+            let rk = bitreverse(k, self.exp);
+            if k < rk {
+                let tmp1 = a[rk];
+                let tmp2 = a[k];
+                a[rk] = tmp2;
+                a[k] = tmp1;
+            }
+        }
+
+        let mut m = 1;
+        for _ in 0..self.exp {
+            let w_m = omega.pow(e, &[(self.m / (2*m)) as u64]);
+
+            let mut k = 0;
+            while k < self.m {
+                let mut w = E::Fr::one(e);
+                for j in 0..m {
+                    let mut t = w;
+                    t.mul_assign(e, &a[(k+j+m) as usize]);
+                    let mut tmp = a[(k+j) as usize];
+                    tmp.sub_assign(e, &t);
+                    a[(k+j+m) as usize] = tmp;
+                    a[(k+j) as usize].add_assign(e, &t);
+                    w.mul_assign(e, &w_m);
+                }
+
+                k += 2*m;
+            }
+
+            m *= 2;
+        }
+    }
+}
+
+// Test multiplying various (low degree) polynomials together and
+// comparing with naive evaluations.
+#[test]
+fn polynomial_arith() {
+    use curves::*;
+    use curves::bls381::Bls381;
+    use rand;
+
+    fn test_mul<E: Engine, R: rand::Rng>(e: &E, rng: &mut R)
+    {
+        for coeffs_a in 1..70 {
+            for coeffs_b in 1..70 {
+                let final_degree = coeffs_a + coeffs_b - 1;
+
+                let domain = EvaluationDomain::new(e, final_degree as u64);
+                let mut a: Vec<_> = (0..coeffs_a).map(|_| E::Fr::random(e, rng)).collect();
+                let mut b: Vec<_> = (0..coeffs_b).map(|_| E::Fr::random(e, rng)).collect();
+
+                // naive evaluation
+                let mut naive = vec![E::Fr::zero(); domain.m as usize];
+                for (i1, a) in a.iter().enumerate() {
+                    for (i2, b) in b.iter().enumerate() {
+                        let mut prod = *a;
+                        prod.mul_assign(e, b);
+                        naive[i1 + i2].add_assign(e, &prod);
+                    }
+                }
+
+                a.resize(domain.m as usize, E::Fr::zero());
+                b.resize(domain.m as usize, E::Fr::zero());
+                let mut c = vec![];
+                c.resize(domain.m as usize, E::Fr::zero());
+
+                domain.fft(e, &mut a);
+                domain.fft(e, &mut b);
+
+                for ((a, b), c) in a.iter().zip(b.iter()).zip(c.iter_mut()) {
+                    *c = *a;
+                    c.mul_assign(e, b);
+                }
+
+                domain.ifft(e, &mut c);
+
+                for (naive, fft) in naive.iter().zip(c.iter()) {
+                    assert_eq!(naive, fft);
+                }
+            }
+        }
+    }
+
+    let e = &Bls381::new();
+    let rng = &mut rand::thread_rng();
+
+    test_mul(e, rng);
+}
diff --git a/src/groth/mod.rs b/src/groth/mod.rs
new file mode 100644
index 0000000..18d9506
--- /dev/null
+++ b/src/groth/mod.rs
@@ -0,0 +1,515 @@
+use curves::*;
+use super::*;
+
+mod domain;
+
+pub struct ProvingKey<E: Engine> {
+    a_inputs: Vec<<E::G1 as Group<E>>::Affine>,
+    b1_inputs: Vec<<E::G1 as Group<E>>::Affine>,
+    b2_inputs: Vec<<E::G2 as Group<E>>::Affine>,
+    a_aux: Vec<<E::G1 as Group<E>>::Affine>,
+    b1_aux: Vec<<E::G1 as Group<E>>::Affine>,
+    b2_aux: Vec<<E::G2 as Group<E>>::Affine>,
+    h: Vec<<E::G1 as Group<E>>::Affine>,
+    l: Vec<<E::G1 as Group<E>>::Affine>,
+    alpha_g1: <E::G1 as Group<E>>::Affine,
+    beta_g1: <E::G1 as Group<E>>::Affine,
+    beta_g2: <E::G2 as Group<E>>::Affine,
+    delta_g1: <E::G1 as Group<E>>::Affine,
+    delta_g2: <E::G2 as Group<E>>::Affine
+}
+
+pub struct VerifyingKey<E: Engine> {
+    alpha_g1: <E::G1 as Group<E>>::Affine,
+    beta_g2: <E::G2 as Group<E>>::Affine,
+    gamma_g2: <E::G2 as Group<E>>::Affine,
+    delta_g2: <E::G2 as Group<E>>::Affine,
+    ic: Vec<<E::G1 as Group<E>>::Affine>
+}
+
+pub struct PreparedVerifyingKey<E: Engine> {
+    alpha_g1_beta_g2: E::Fqk,
+    neg_gamma_g2: <E::G2 as Group<E>>::Prepared,
+    neg_delta_g2: <E::G2 as Group<E>>::Prepared,
+    ic: Vec<<E::G1 as Group<E>>::Affine>
+}
+
+pub struct Proof<E: Engine> {
+    a: E::G1,
+    b: E::G2,
+    c: E::G1
+}
+
+pub fn keypair<E: Engine, C: Circuit<E>>(
+    e: &E,
+    circuit: C,
+    tau: &E::Fr,
+    alpha: &E::Fr,
+    beta: &E::Fr,
+    gamma: &E::Fr,
+    delta: &E::Fr
+) -> (ProvingKey<E>, VerifyingKey<E>)
+{
+    struct KeypairAssembly<E: Engine> {
+        num_inputs: usize,
+        num_aux: usize,
+        num_constraints: usize,
+        at_inputs: Vec<Vec<(E::Fr, usize)>>,
+        bt_inputs: Vec<Vec<(E::Fr, usize)>>,
+        ct_inputs: Vec<Vec<(E::Fr, usize)>>,
+        at_aux: Vec<Vec<(E::Fr, usize)>>,
+        bt_aux: Vec<Vec<(E::Fr, usize)>>,
+        ct_aux: Vec<Vec<(E::Fr, usize)>>
+    }
+
+    impl<E: Engine> PublicConstraintSystem<E> for KeypairAssembly<E> {
+        fn alloc_input(&mut self, _: E::Fr) -> Variable {
+            let index = self.num_inputs;
+            self.num_inputs += 1;
+
+            self.at_inputs.push(vec![]);
+            self.bt_inputs.push(vec![]);
+            self.ct_inputs.push(vec![]);
+
+            Variable(Index::Input(index))
+        }
+    }
+
+    impl<E: Engine> ConstraintSystem<E> for KeypairAssembly<E> {
+        fn alloc(&mut self, _: E::Fr) -> Variable {
+            let index = self.num_aux;
+            self.num_aux += 1;
+
+            self.at_aux.push(vec![]);
+            self.bt_aux.push(vec![]);
+            self.ct_aux.push(vec![]);
+
+            Variable(Index::Aux(index))
+        }
+
+        fn enforce(
+            &mut self,
+            a: LinearCombination<E>,
+            b: LinearCombination<E>,
+            c: LinearCombination<E>
+        )
+        {
+            fn qap_eval<E: Engine>(
+                l: LinearCombination<E>,
+                inputs: &mut [Vec<(E::Fr, usize)>],
+                aux: &mut [Vec<(E::Fr, usize)>],
+                this_constraint: usize
+            )
+            {
+                for (index, coeff) in l.0 {
+                    match index {
+                        Index::Input(id) => inputs[id].push((coeff, this_constraint)),
+                        Index::Aux(id) => aux[id].push((coeff, this_constraint))
+                    }
+                }
+            }
+
+            qap_eval(a, &mut self.at_inputs, &mut self.at_aux, self.num_constraints);
+            qap_eval(b, &mut self.bt_inputs, &mut self.bt_aux, self.num_constraints);
+            qap_eval(c, &mut self.ct_inputs, &mut self.ct_aux, self.num_constraints);
+
+            self.num_constraints += 1;
+        }
+    }
+
+    let mut assembly = KeypairAssembly {
+        num_inputs: 0,
+        num_aux: 0,
+        num_constraints: 0,
+        at_inputs: vec![],
+        bt_inputs: vec![],
+        ct_inputs: vec![],
+        at_aux: vec![],
+        bt_aux: vec![],
+        ct_aux: vec![]
+    };
+
+    assembly.alloc_input(E::Fr::one(e));
+
+    circuit.synthesize(e, &mut assembly).synthesize(e, &mut assembly);
+
+    // Input consistency constraints: x * 0 = 0
+    for i in 0..assembly.num_inputs {
+        assembly.enforce(LinearCombination::zero(e).add(E::Fr::one(e), Variable(Index::Input(i))),
+                         LinearCombination::zero(e),
+                         LinearCombination::zero(e));
+    }
+
+    let domain = domain::EvaluationDomain::new(e, assembly.num_constraints as u64);
+
+    let mut u = Vec::with_capacity(domain.m as usize);
+    {
+        let mut acc = E::Fr::one(e);
+        for _ in 0..domain.m {
+            u.push(acc);
+            acc.mul_assign(e, tau);
+        }
+    }
+
+    let gamma_inverse = gamma.inverse(e).unwrap();
+    let delta_inverse = delta.inverse(e).unwrap();
+
+    let g1_table;
+    let h;
+    {
+        let mut powers_of_tau = u.clone();
+        powers_of_tau.truncate((domain.m - 1) as usize);
+
+        let mut coeff = delta_inverse;
+        coeff.mul_assign(e, &domain.z(e, tau));
+        for h in &mut powers_of_tau {
+            h.mul_assign(e, &coeff);
+        }
+
+        g1_table = E::G1::one(e).optimal_window_batch(e,
+            (domain.m - 1) as usize + (assembly.num_inputs + assembly.num_aux) * 3
+        );
+
+        h = e.batch_baseexp(&g1_table, powers_of_tau);
+    }
+
+    domain.ifft(e, &mut u);
+
+    fn eval<E: Engine>(
+        e: &E,
+        u: &[E::Fr],
+        alpha: &E::Fr,
+        beta: &E::Fr,
+        inv: &E::Fr,
+        at_in: Vec<Vec<(E::Fr, usize)>>,
+        bt_in: Vec<Vec<(E::Fr, usize)>>,
+        ct_in: Vec<Vec<(E::Fr, usize)>>
+    ) -> (Vec<E::Fr>, Vec<E::Fr>, Vec<E::Fr>)
+    {
+        assert_eq!(at_in.len(), bt_in.len());
+        assert_eq!(bt_in.len(), ct_in.len());
+
+        fn eval_at_tau<E: Engine>(
+            e: &E,
+            val: Vec<(E::Fr, usize)>,
+            u: &[E::Fr]
+        ) -> E::Fr
+        {
+            let mut acc = E::Fr::zero();
+
+            for (coeff, index) in val {
+                let mut n = u[index];
+                n.mul_assign(e, &coeff);
+                acc.add_assign(e, &n);
+            }
+
+            acc
+        }
+
+        let mut a_out = Vec::with_capacity(at_in.len());
+        let mut b_out = Vec::with_capacity(at_in.len());
+        let mut l_out = Vec::with_capacity(at_in.len());
+
+        for ((a, b), c) in at_in.into_iter().zip(bt_in.into_iter()).zip(ct_in.into_iter()) {
+            let a = eval_at_tau(e, a, u);
+            let b = eval_at_tau(e, b, u);
+
+            let mut t0 = a;
+            t0.mul_assign(e, beta);
+
+            let mut t1 = b;
+            t1.mul_assign(e, alpha);
+
+            t0.add_assign(e, &t1);
+            t0.add_assign(e, &eval_at_tau(e, c, u));
+            t0.mul_assign(e, inv);
+
+            a_out.push(a);
+            b_out.push(b);
+            l_out.push(t0);
+        }
+
+        (a_out, b_out, l_out)
+    }
+
+    let (a_inputs, b_inputs, ic_coeffs) = eval(e, &u, alpha, beta, &gamma_inverse, assembly.at_inputs, assembly.bt_inputs, assembly.ct_inputs);
+    let a_inputs = e.batch_baseexp(&g1_table, a_inputs);
+    let b1_inputs = e.batch_baseexp(&g1_table, &b_inputs);
+    let ic_coeffs = e.batch_baseexp(&g1_table, ic_coeffs);
+    let (a_aux, b_aux, l_coeffs) = eval(e, &u, alpha, beta, &delta_inverse, assembly.at_aux, assembly.bt_aux, assembly.ct_aux);
+    let a_aux = e.batch_baseexp(&g1_table, a_aux);
+    let b1_aux = e.batch_baseexp(&g1_table, &b_aux);
+    let l_coeffs = e.batch_baseexp(&g1_table, l_coeffs);
+
+    drop(g1_table);
+
+    let g2_table = E::G2::one(e).optimal_window_batch(e,
+        (assembly.num_inputs + assembly.num_aux)
+    );
+
+    let b2_inputs = e.batch_baseexp(&g2_table, b_inputs);
+    let b2_aux = e.batch_baseexp(&g2_table, b_aux);
+
+    let mut alpha_g1 = E::G1::one(e);
+    alpha_g1.mul_assign(e, alpha);
+    let mut beta_g1 = E::G1::one(e);
+    beta_g1.mul_assign(e, beta);
+    let mut beta_g2 = E::G2::one(e);
+    beta_g2.mul_assign(e, beta);
+    let mut gamma_g2 = E::G2::one(e);
+    gamma_g2.mul_assign(e, gamma);
+    let mut delta_g1 = E::G1::one(e);
+    delta_g1.mul_assign(e, delta);
+    let mut delta_g2 = E::G2::one(e);
+    delta_g2.mul_assign(e, delta);
+
+    (
+        ProvingKey {
+            a_inputs: a_inputs,
+            b1_inputs: b1_inputs,
+            b2_inputs: b2_inputs,
+            a_aux: a_aux,
+            b1_aux: b1_aux,
+            b2_aux: b2_aux,
+            h: h,
+            l: l_coeffs,
+            delta_g1: delta_g1.to_affine(e),
+            delta_g2: delta_g2.to_affine(e),
+            alpha_g1: alpha_g1.to_affine(e),
+            beta_g1: beta_g1.to_affine(e),
+            beta_g2: beta_g2.to_affine(e)
+        },
+        VerifyingKey {
+            alpha_g1: alpha_g1.to_affine(e),
+            beta_g2: beta_g2.to_affine(e),
+            gamma_g2: gamma_g2.to_affine(e),
+            delta_g2: delta_g2.to_affine(e),
+            ic: ic_coeffs
+        }
+    )
+}
+
+pub fn prepare_verifying_key<E: Engine>(
+    e: &E,
+    vk: &VerifyingKey<E>
+) -> PreparedVerifyingKey<E>
+{
+    let mut gamma = vk.gamma_g2;
+    gamma.negate(e);
+    let mut delta = vk.delta_g2;
+    delta.negate(e);
+
+    PreparedVerifyingKey {
+        alpha_g1_beta_g2: e.pairing(&vk.alpha_g1, &vk.beta_g2),
+        neg_gamma_g2: gamma.prepare(e),
+        neg_delta_g2: delta.prepare(e),
+        ic: vk.ic.clone()
+    }
+}
+
+pub fn verify<E: Engine, C: Witness<E>, F: FnOnce(&mut ConstraintSystem<E>) -> C>(
+    e: &E,
+    circuit: F,
+    proof: &Proof<E>,
+    pvk: &PreparedVerifyingKey<E>
+) -> bool
+{
+    struct VerifierInput<'a, E: Engine + 'a> {
+        e: &'a E,
+        acc: E::G1,
+        ic: &'a [<E::G1 as Group<E>>::Affine],
+        insufficient_inputs: bool,
+        num_inputs: usize,
+        num_aux: usize
+    }
+
+    impl<'a, E: Engine> PublicConstraintSystem<E> for VerifierInput<'a, E> {
+        fn alloc_input(&mut self, value: E::Fr) -> Variable {
+            if self.ic.len() == 0 {
+                self.insufficient_inputs = true;
+            } else {
+                self.acc.add_assign(self.e, &self.ic[0].mul(self.e, &value));
+                self.ic = &self.ic[1..];
+            }
+
+            let index = self.num_inputs;
+            self.num_inputs += 1;
+
+            Variable(Index::Input(index))
+        }
+    }
+
+    impl<'a, E: Engine> ConstraintSystem<E> for VerifierInput<'a, E> {
+        fn alloc(&mut self, _: E::Fr) -> Variable {
+            let index = self.num_aux;
+            self.num_aux += 1;
+
+            Variable(Index::Aux(index))
+        }
+
+        fn enforce(
+            &mut self,
+            _: LinearCombination<E>,
+            _: LinearCombination<E>,
+            _: LinearCombination<E>
+        )
+        {
+            // Do nothing; we don't care about the constraint system
+            // in this context.
+        }
+    }
+
+    let mut witness = VerifierInput {
+        e: e,
+        acc: pvk.ic[0].to_jacobian(e),
+        ic: &pvk.ic[1..],
+        insufficient_inputs: false,
+        num_inputs: 1,
+        num_aux: 0
+    };
+
+    circuit(&mut witness).synthesize(e, &mut witness);
+
+    if witness.ic.len() != 0 || witness.insufficient_inputs {
+        return false;
+    }
+
+    e.final_exponentiation(
+        &e.miller_loop([
+            (&proof.a.prepare(e), &proof.b.prepare(e)),
+            (&witness.acc.prepare(e), &pvk.neg_gamma_g2),
+            (&proof.c.prepare(e), &pvk.neg_delta_g2)
+        ].into_iter())
+    ) == pvk.alpha_g1_beta_g2
+}
+
+pub fn prove<E: Engine, C: Circuit<E>>(
+    e: &E,
+    circuit: C,
+    r: &E::Fr,
+    s: &E::Fr,
+    pk: &ProvingKey<E>
+) -> Result<Proof<E>, ()>
+{
+    struct ProvingAssignment<'a, E: Engine + 'a> {
+        e: &'a E,
+        // Evaluations of A, B, C polynomials
+        a: Vec<E::Fr>,
+        b: Vec<E::Fr>,
+        c: Vec<E::Fr>,
+        // Assignments of variables
+        input_assignment: Vec<E::Fr>,
+        aux_assignment: Vec<E::Fr>
+    }
+
+    impl<'a, E: Engine> PublicConstraintSystem<E> for ProvingAssignment<'a, E> {
+        fn alloc_input(&mut self, value: E::Fr) -> Variable {
+            self.input_assignment.push(value);
+
+            Variable(Index::Input(self.input_assignment.len() - 1))
+        }
+    }
+
+    impl<'a, E: Engine> ConstraintSystem<E> for ProvingAssignment<'a, E> {
+        fn alloc(&mut self, value: E::Fr) -> Variable {
+            self.aux_assignment.push(value);
+
+            Variable(Index::Aux(self.aux_assignment.len() - 1))
+        }
+
+        fn enforce(
+            &mut self,
+            a: LinearCombination<E>,
+            b: LinearCombination<E>,
+            c: LinearCombination<E>
+        )
+        {
+            self.a.push(a.evaluate(self.e, &self.input_assignment, &self.aux_assignment));
+            self.b.push(b.evaluate(self.e, &self.input_assignment, &self.aux_assignment));
+            self.c.push(c.evaluate(self.e, &self.input_assignment, &self.aux_assignment));
+        }
+    }
+
+    let mut prover = ProvingAssignment {
+        e: e,
+        a: vec![],
+        b: vec![],
+        c: vec![],
+        input_assignment: vec![],
+        aux_assignment: vec![]
+    };
+
+    prover.alloc_input(E::Fr::one(e));
+
+    circuit.synthesize(e, &mut prover).synthesize(e, &mut prover);
+
+    // Input consistency constraints: x * 0 = 0
+    for i in 0..prover.input_assignment.len() {
+        prover.enforce(LinearCombination::zero(e).add(E::Fr::one(e), Variable(Index::Input(i))),
+                       LinearCombination::zero(e),
+                       LinearCombination::zero(e));
+    }
+
+    // Perform FFTs
+    let h = {
+        let domain = domain::EvaluationDomain::new(e, prover.a.len() as u64);
+        prover.a.resize(domain.m as usize, E::Fr::zero());
+        prover.b.resize(domain.m as usize, E::Fr::zero());
+        prover.c.resize(domain.m as usize, E::Fr::zero());
+        domain.ifft(e, &mut prover.a);
+        domain.coset_fft(e, &mut prover.a);
+        domain.ifft(e, &mut prover.b);
+        domain.coset_fft(e, &mut prover.b);
+        domain.ifft(e, &mut prover.c);
+        domain.coset_fft(e, &mut prover.c);
+
+        let mut h = prover.a;
+        for (h, b) in h.iter_mut().zip(prover.b.into_iter()) {
+            h.mul_assign(e, &b);
+        }
+        for (h, c) in h.iter_mut().zip(prover.c.into_iter()) {
+            h.sub_assign(e, &c);
+        }
+        domain.divide_by_z_on_coset(e, &mut h);
+        domain.icoset_fft(e, &mut h);
+
+        e.multiexp(&pk.h, &h[0..(domain.m-1) as usize])?
+    };
+
+    // Construct proof
+    let mut g_a = pk.delta_g1.mul(e, r);
+    g_a.add_assign(e, &pk.alpha_g1.to_jacobian(e));
+    let mut g_b = pk.delta_g2.mul(e, s);
+    g_b.add_assign(e, &pk.beta_g2.to_jacobian(e));
+    let mut g_c;
+    {
+        let mut rs = *r;
+        rs.mul_assign(e, s);
+        g_c = pk.delta_g1.mul(e, &rs);
+        g_c.add_assign(e, &pk.alpha_g1.mul(e, s));
+        g_c.add_assign(e, &pk.beta_g1.mul(e, r));
+    }
+    let mut a_answer: E::G1 = e.multiexp(&pk.a_inputs, &prover.input_assignment)?;
+    a_answer.add_assign(e, &e.multiexp(&pk.a_aux, &prover.aux_assignment)?);
+    g_a.add_assign(e, &a_answer);
+    a_answer.mul_assign(e, s);
+    g_c.add_assign(e, &a_answer);
+    let mut b1_answer: E::G1 = e.multiexp(&pk.b1_inputs, &prover.input_assignment)?;
+    b1_answer.add_assign(e, &e.multiexp(&pk.b1_aux, &prover.aux_assignment)?);
+    let mut b2_answer: E::G2 = e.multiexp(&pk.b2_inputs, &prover.input_assignment)?;
+    b2_answer.add_assign(e, &e.multiexp(&pk.b2_aux, &prover.aux_assignment)?);
+    g_b.add_assign(e, &b2_answer);
+    b1_answer.mul_assign(e, r);
+    g_c.add_assign(e, &b1_answer);
+    g_c.add_assign(e, &h);
+    g_c.add_assign(e, &e.multiexp(&pk.l, &prover.aux_assignment)?);
+
+    Ok(Proof {
+        a: g_a,
+        b: g_b,
+        c: g_c
+    })
+}
+
+#[cfg(test)]
+mod tests;
diff --git a/src/groth/tests/mod.rs b/src/groth/tests/mod.rs
new file mode 100644
index 0000000..e486c5f
--- /dev/null
+++ b/src/groth/tests/mod.rs
@@ -0,0 +1,207 @@
+use super::*;
+use rand::{Rng, thread_rng};
+
+struct RootCircuit<E: Engine> {
+    root: E::Fr
+}
+
+impl<E: Engine> Circuit<E> for RootCircuit<E> {
+    type WitnessMap = RootWitness<E>;
+
+    fn synthesize<CS: ConstraintSystem<E>>(self,
+                                           e: &E,
+                                           cs: &mut CS)
+                                           -> Self::WitnessMap
+    {
+        let root_var = cs.alloc(self.root);
+
+        let mut cur = root_var;
+        let mut cur_val = self.root;
+
+        for _ in 0..99 {
+            cur_val.mul_assign(e, &self.root);
+            let new = cs.alloc(cur_val);
+
+            cs.enforce(
+                LinearCombination::zero(e) + (E::Fr::from_str(e, "3").unwrap(), cur),
+                LinearCombination::zero(e) + (E::Fr::from_str(e, "4").unwrap(), root_var),
+                LinearCombination::zero(e) + (E::Fr::from_str(e, "12").unwrap(), new),
+            );
+
+            cur = new;
+        }
+
+        RootWitness {
+            num: cur_val,
+            num_var: cur
+        }
+    }
+}
+
+struct RootWitness<E: Engine> {
+    num: E::Fr,
+    num_var: Variable
+}
+
+impl<E: Engine> Witness<E> for RootWitness<E> {
+    fn synthesize<CS: PublicConstraintSystem<E>>(
+        self,
+        e: &E,
+        cs: &mut CS
+    )
+    {
+        let result_input = cs.alloc_input(self.num);
+        cs.enforce(
+            LinearCombination::zero(e) + result_input,
+            LinearCombination::one(e),
+            LinearCombination::zero(e) + self.num_var
+        );
+    }
+}
+
+fn test_snark_system<E: Engine, R: Rng>(
+    e: &E,
+    rng: &mut R
+)
+{
+    let tau = E::Fr::random(e, rng);
+    let alpha = E::Fr::random(e, rng);
+    let beta = E::Fr::random(e, rng);
+    let gamma = E::Fr::random(e, rng);
+    let delta = E::Fr::random(e, rng);
+
+    // create keypair
+    let (pk, vk) = {
+        let c = RootCircuit {
+            root: E::Fr::zero()
+        };
+
+        keypair(e, c, &tau, &alpha, &beta, &gamma, &delta)
+    };
+
+    // construct proof
+    let proof = {
+        let r = E::Fr::random(e, rng);
+        let s = E::Fr::random(e, rng);
+
+        let c = RootCircuit {
+            root: E::Fr::from_str(e, "2").unwrap()
+        };
+
+        prove(e, c, &r, &s, &pk).unwrap()
+    };
+
+    // prepare verifying key
+    let pvk = prepare_verifying_key(e, &vk);
+
+    // verify proof
+    assert!(verify(e, |cs| {
+        RootWitness {
+            num: E::Fr::from_str(e, "1267650600228229401496703205376").unwrap(),
+            num_var: cs.alloc(E::Fr::one(e))
+        }
+    }, &proof, &pvk));
+
+    // verify invalid proof
+    assert!(!verify(e, |cs| {
+        RootWitness {
+            num: E::Fr::from_str(e, "1267650600228229401496703205375").unwrap(),
+            num_var: cs.alloc(E::Fr::one(e))
+        }
+    }, &proof, &pvk));
+
+    // simulate a groth proof with trapdoors
+    // ----------------
+    // 99: a1 * a0 = l*
+    // 100: a0 * 0 = 0
+    // 101: a1 * 0 = 0
+    // ---
+    // u_0(tau) = tau^100
+    // u_1(tau) = tau^99 + tau^101
+    // v_0(tau) = tau^99
+    // v_1(tau) = 0
+    // w_0(tau) = 0
+    // w_1(tau) = 0
+    // ---
+
+    let mut lagrange_coeffs: Vec<E::Fr> = (0..128).map(|i| tau.pow(e, &[i])).collect();
+
+    let d = domain::EvaluationDomain::new(e, 128);
+    d.ifft(e, &mut lagrange_coeffs);
+
+    let a = E::Fr::random(e, rng);
+    let b = E::Fr::random(e, rng);
+
+    let mut c = a;
+    c.mul_assign(e, &b);
+
+    let mut alphabeta = alpha;
+    alphabeta.mul_assign(e, &beta);
+    c.sub_assign(e, &alphabeta);
+
+    let mut ic = E::Fr::zero();
+    {
+        let mut ic_i_beta = lagrange_coeffs[100];
+        ic_i_beta.mul_assign(e, &beta);
+
+        let mut ic_i_alpha = lagrange_coeffs[99];
+        ic_i_alpha.mul_assign(e, &alpha);
+
+        ic_i_beta.add_assign(e, &ic_i_alpha);
+
+        ic.add_assign(e, &ic_i_beta);
+    }
+    {
+        let mut ic_i_beta = lagrange_coeffs[99];
+        ic_i_beta.add_assign(e, &lagrange_coeffs[101]);
+        ic_i_beta.mul_assign(e, &beta);
+
+        ic_i_beta.mul_assign(e, &E::Fr::from_str(e, "100").unwrap());
+
+        ic.add_assign(e, &ic_i_beta);
+    }
+
+    c.sub_assign(e, &ic);
+    c.mul_assign(e, &delta.inverse(e).unwrap());
+
+    let mut a_g = E::G1::one(e);
+    a_g.mul_assign(e, &a);
+
+    let mut b_g = E::G2::one(e);
+    b_g.mul_assign(e, &b);
+
+    let mut c_g = E::G1::one(e);
+    c_g.mul_assign(e, &c);
+
+    let fake_proof = Proof {
+        a: a_g,
+        b: b_g,
+        c: c_g
+    };
+
+    // verify fake proof
+    assert!(verify(e, |cs| {
+        RootWitness {
+            num: E::Fr::from_str(e, "100").unwrap(),
+            num_var: cs.alloc(E::Fr::one(e))
+        }
+    }, &fake_proof, &pvk));
+
+    // verify fake proof with wrong input
+    assert!(!verify(e, |cs| {
+        RootWitness {
+            num: E::Fr::from_str(e, "101").unwrap(),
+            num_var: cs.alloc(E::Fr::one(e))
+        }
+    }, &fake_proof, &pvk));
+}
+
+#[test]
+fn groth_with_bls381() {
+    use curves::bls381::Bls381;
+
+    let e = &Bls381::new();
+    let rng = &mut thread_rng();
+
+    test_snark_system(e, rng);
+}
diff --git a/src/lib.rs b/src/lib.rs
index 102d2ff..d72b212 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -1,8 +1,151 @@
 #![feature(i128_type)]
 
 extern crate rand;
-extern crate rayon;
+extern crate num_cpus;
+extern crate crossbeam;
 extern crate byteorder;
 extern crate serde;
 
 pub mod curves;
+pub mod groth;
+
+use std::collections::HashMap;
+use std::ops;
+
+use curves::{Engine, Field};
+
+#[derive(Copy, Clone)]
+pub struct Variable(Index);
+
+impl Variable {
+    pub fn one() -> Self {
+        Variable(Index::Input(0))
+    }
+}
+
+#[derive(Clone, Copy, PartialEq, Eq, Hash)]
+enum Index {
+    Input(usize),
+    Aux(usize)
+}
+
+pub struct LinearCombination<'a, E: Engine + 'a>(HashMap<Index, E::Fr>, &'a E);
+
+impl<'a, E: Engine + 'a> ops::Add<Variable> for LinearCombination<'a, E> {
+    type Output = LinearCombination<'a, E>;
+
+    fn add(self, other: Variable) -> LinearCombination<'a, E> {
+        let one = E::Fr::one(self.1);
+
+        self.add(one, other)
+    }
+}
+
+impl<'a, E: Engine + 'a> ops::Add<(E::Fr, Variable)> for LinearCombination<'a, E> {
+    type Output = LinearCombination<'a, E>;
+
+    fn add(self, (coeff, var): (E::Fr, Variable)) -> LinearCombination<'a, E> {
+        self.add(coeff, var)
+    }
+}
+
+impl<'a, E: Engine + 'a> ops::Sub<Variable> for LinearCombination<'a, E> {
+    type Output = LinearCombination<'a, E>;
+
+    fn sub(self, other: Variable) -> LinearCombination<'a, E> {
+        let one = E::Fr::one(self.1);
+
+        self.sub(one, other)
+    }
+}
+
+impl<'a, E: Engine + 'a> ops::Sub<(E::Fr, Variable)> for LinearCombination<'a, E> {
+    type Output = LinearCombination<'a, E>;
+
+    fn sub(self, (coeff, var): (E::Fr, Variable)) -> LinearCombination<'a, E> {
+        self.sub(coeff, var)
+    }
+}
+
+impl<'a, E: Engine> LinearCombination<'a, E> {
+    pub fn zero(e: &'a E) -> LinearCombination<'a, E> {
+        LinearCombination(HashMap::new(), e)
+    }
+
+    pub fn one(e: &'a E) -> LinearCombination<'a, E> {
+        LinearCombination::zero(e).add(E::Fr::one(e), Variable::one())
+    }
+
+    pub fn add(mut self, coeff: E::Fr, var: Variable) -> Self
+    {
+        self.0.entry(var.0)
+              .or_insert(E::Fr::zero())
+              .add_assign(self.1, &coeff);
+
+        self
+    }
+
+    pub fn sub(self, mut coeff: E::Fr, var: Variable) -> Self
+    {
+        coeff.negate(self.1);
+
+        self.add(coeff, var)
+    }
+
+    fn evaluate(
+        &self,
+        e: &E,
+        input_assignment: &[E::Fr],
+        aux_assignment: &[E::Fr]
+    ) -> E::Fr
+    {
+        let mut acc = E::Fr::zero();
+        for (index, coeff) in self.0.iter() {
+            let mut n = *coeff;
+            match index {
+                &Index::Input(id) => {
+                    n.mul_assign(e, &input_assignment[id]);
+                },
+                &Index::Aux(id) => {
+                    n.mul_assign(e, &aux_assignment[id]);
+                }
+            }
+            acc.add_assign(e, &n);
+        }
+
+        acc
+    }
+}
+
+pub trait Circuit<E: Engine> {
+    type WitnessMap: Witness<E>;
+
+    /// Synthesize the circuit into a rank-1 quadratic constraint system
+    #[must_use]
+    fn synthesize<CS: ConstraintSystem<E>>(self, engine: &E, cs: &mut CS) -> Self::WitnessMap;
+}
+
+pub trait Witness<E: Engine> {
+    /// Synthesize the circuit, except with additional access to public input
+    /// variables
+    fn synthesize<CS: PublicConstraintSystem<E>>(self, engine: &E, cs: &mut CS);
+}
+
+pub trait PublicConstraintSystem<E: Engine>: ConstraintSystem<E> {
+    /// Allocate a public input that the verifier knows.
+    fn alloc_input(&mut self, value: E::Fr) -> Variable;
+}
+
+pub trait ConstraintSystem<E: Engine> {
+    /// Allocate a private variable in the constraint system, setting it to
+    /// the provided value.
+    fn alloc(&mut self, value: E::Fr) -> Variable;
+    
+    /// Enforce that `A` * `B` = `C`.
+    fn enforce(
+        &mut self,
+        a: LinearCombination<E>,
+        b: LinearCombination<E>,
+        c: LinearCombination<E>
+    );
+}