diff --git a/Cargo.toml b/Cargo.toml
index fe0a0ee..56709c9 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -11,8 +11,8 @@ license-file = "LICENSE"
 keywords = ["zkSNARKs", "cryptography", "proofs"]
 
 [dependencies]
-bellperson =  "0.19.0"
-ff = "0.11.0"
+bellperson =  "0.19.1"
+ff = {version = "0.11.0", features = ["derive"]}
 merlin = "2.0.0"
 rand = "0.8.4"
 digest = "0.8.1"
@@ -22,3 +22,7 @@ rand_core = { version = "0.5", default-features = false }
 itertools = "0.9.0"
 subtle = "2.4"
 pasta_curves = "0.3.0"
+neptune = "6"
+generic-array = "0.14.4"
+bellperson-nonnative = "0.1.0"
+rug = { version = "1.10", default-features = false, features = ["integer", "serde", "rand"] }
diff --git a/src/circuit.rs b/src/circuit.rs
new file mode 100644
index 0000000..24e7125
--- /dev/null
+++ b/src/circuit.rs
@@ -0,0 +1,805 @@
+//! There are two Verification Circuits. Each of them is over a Pasta curve but
+//! only one of them executes the next step of the computation by applying the inner function F.
+//! There are also two running relaxed r1cs instances.
+//!
+//! When we build a circuit we denote u1 the running relaxed r1cs instance of
+//! the circuit and u2 the running relaxed r1cs instance of the other circuit.
+//! The circuit takes as input two hashes h1 and h2.
+//! If the circuit applies the inner function F, then
+//! h1 = H(params = H(shape, gens), u2, i, z0, zi) and h2 = H(u1, i)
+//! otherwise
+//! h1 = H(u2, i) and h2 = H(params = H(shape, gens), u1, i, z0, zi)
+
+use super::commitments::Commitment;
+use super::gadgets::{
+  ecc_circuit::AllocatedPoint,
+  utils::{
+    alloc_bignat_constant, alloc_num_equals, alloc_one, alloc_zero, conditionally_select,
+    conditionally_select_bignat, le_bits_to_num,
+  },
+};
+use super::poseidon::NovaPoseidonConstants;
+use super::poseidon::PoseidonROGadget;
+use super::r1cs::RelaxedR1CSInstance;
+use super::traits::{Group, InnerCircuit, PrimeField};
+use bellperson::{
+  gadgets::{boolean::Boolean, num::AllocatedNum, Assignment},
+  Circuit, ConstraintSystem, SynthesisError,
+};
+use bellperson_nonnative::{
+  mp::bignat::BigNat,
+  util::{convert::f_to_nat, num::Num},
+};
+use ff::PrimeFieldBits;
+
+#[derive(Debug, Clone)]
+pub struct VerificationCircuitParams {
+  limb_width: usize,
+  n_limbs: usize,
+}
+
+impl VerificationCircuitParams {
+  #[allow(dead_code)]
+  pub fn new(limb_width: usize, n_limbs: usize) -> Self {
+    Self {
+      limb_width,
+      n_limbs,
+    }
+  }
+}
+
+pub struct VerificationCircuitInputs<G>
+where
+  G: Group,
+{
+  h1: G::Base,
+  h2: G::Scalar,
+  u2: RelaxedR1CSInstance<G>,
+  i: G::Base,
+  z0: G::Base,
+  zi: G::Base,
+  params: G::Base, // Hash(Shape of u2, Gens for u2). Needed for computing the challenge.
+  T: Commitment<G>,
+  w: Commitment<G>, // The commitment to the witness of the fresh r1cs instance
+}
+
+impl<G> VerificationCircuitInputs<G>
+where
+  G: Group,
+{
+  /// Create new inputs/witness for the verification circuit
+  #[allow(dead_code, clippy::too_many_arguments)]
+  pub fn new(
+    h1: G::Base,
+    u2: RelaxedR1CSInstance<G>,
+    i: G::Base,
+    z0: G::Base,
+    zi: G::Base,
+    h2: G::Scalar,
+    params: G::Base,
+    T: Commitment<G>,
+    w: Commitment<G>,
+  ) -> Self {
+    Self {
+      h1,
+      u2,
+      i,
+      z0,
+      zi,
+      h2,
+      params,
+      T,
+      w,
+    }
+  }
+}
+
+/// Circuit that encodes only the folding verifier
+pub struct VerificationCircuit<G, IC>
+where
+  G: Group,
+  <G as Group>::Base: ff::PrimeField,
+  IC: InnerCircuit<G::Base>,
+{
+  params: VerificationCircuitParams,
+  inputs: Option<VerificationCircuitInputs<G>>,
+  inner_circuit: Option<IC>, // The function that is applied for each step. may be None.
+  poseidon_constants: NovaPoseidonConstants<G::Base>,
+}
+
+impl<G, IC> VerificationCircuit<G, IC>
+where
+  G: Group,
+  <G as Group>::Base: ff::PrimeField,
+  IC: InnerCircuit<G::Base>,
+{
+  /// Create a new verification circuit for the input relaxed r1cs instances
+  #[allow(dead_code)]
+  pub fn new(
+    params: VerificationCircuitParams,
+    inputs: Option<VerificationCircuitInputs<G>>,
+    inner_circuit: Option<IC>,
+    poseidon_constants: NovaPoseidonConstants<G::Base>,
+  ) -> Self
+  where
+    <G as Group>::Base: ff::PrimeField,
+  {
+    Self {
+      params,
+      inputs,
+      inner_circuit,
+      poseidon_constants,
+    }
+  }
+}
+
+impl<G, IC> Circuit<<G as Group>::Base> for VerificationCircuit<G, IC>
+where
+  G: Group,
+  <G as Group>::Base: ff::PrimeField + PrimeField + PrimeFieldBits,
+  <G as Group>::Scalar: PrimeFieldBits,
+  IC: InnerCircuit<G::Base>,
+{
+  fn synthesize<CS: ConstraintSystem<<G as Group>::Base>>(
+    self,
+    cs: &mut CS,
+  ) -> Result<(), SynthesisError> {
+    /***********************************************************************/
+    // This circuit does not modify h2 but it outputs it.
+    // Allocate it and output it.
+    /***********************************************************************/
+
+    // Allocate h2 as a big number with 8 limbs
+    let h2_bn = BigNat::alloc_from_nat(
+      cs.namespace(|| "allocate h2"),
+      || Ok(f_to_nat(&self.inputs.get()?.h2)),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    let _ = h2_bn.inputize(cs.namespace(|| "Output 1"))?;
+
+    /***********************************************************************/
+    // Allocate h1
+    /***********************************************************************/
+
+    let h1 = AllocatedNum::alloc(cs.namespace(|| "allocate h1"), || Ok(self.inputs.get()?.h1))?;
+    let h1_bn = BigNat::from_num(
+      cs.namespace(|| "allocate h1_bn"),
+      Num::from(h1.clone()),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    /***********************************************************************/
+    // Allocate u2 by allocating W_r, E_r, u_r, X_r
+    /***********************************************************************/
+
+    // W_r = (x, y, infinity)
+    let W_r_x = AllocatedNum::alloc(cs.namespace(|| "W_r.x"), || {
+      Ok(self.inputs.get()?.u2.comm_W.comm.to_coordinates().0)
+    })?;
+    let W_r_y = AllocatedNum::alloc(cs.namespace(|| "W_r.y"), || {
+      Ok(self.inputs.get()?.u2.comm_W.comm.to_coordinates().1)
+    })?;
+    let W_r_inf = AllocatedNum::alloc(cs.namespace(|| "W_r.inf"), || {
+      let infty = self.inputs.get()?.u2.comm_W.comm.to_coordinates().2;
+      if infty {
+        Ok(G::Base::one())
+      } else {
+        Ok(G::Base::zero())
+      }
+    })?;
+
+    let W_r = AllocatedPoint::new(W_r_x.clone(), W_r_y.clone(), W_r_inf.clone());
+    let _ = W_r.check_is_infinity(cs.namespace(|| "W_r check is_infinity"))?;
+
+    // E_r = (x, y, infinity)
+    let E_r_x = AllocatedNum::alloc(cs.namespace(|| "E_r.x"), || {
+      Ok(self.inputs.get()?.u2.comm_E.comm.to_coordinates().0)
+    })?;
+    let E_r_y = AllocatedNum::alloc(cs.namespace(|| "E_r.y"), || {
+      Ok(self.inputs.get()?.u2.comm_E.comm.to_coordinates().1)
+    })?;
+    let E_r_inf = AllocatedNum::alloc(cs.namespace(|| "E_r.inf"), || {
+      let infty = self.inputs.get()?.u2.comm_E.comm.to_coordinates().2;
+      if infty {
+        Ok(G::Base::one())
+      } else {
+        Ok(G::Base::zero())
+      }
+    })?;
+
+    let E_r = AllocatedPoint::new(E_r_x.clone(), E_r_y.clone(), E_r_inf.clone());
+    let _ = E_r.check_is_infinity(cs.namespace(|| "E_r check is_infinity"));
+
+    // u_r << |G::Base| despite the fact that u_r is a scalar.
+    // So we parse all of its bytes as a G::Base element
+    let u_r = AllocatedNum::alloc(cs.namespace(|| "u_r"), || {
+      let u_bits = self.inputs.get()?.u2.u.to_le_bits();
+      let mut mult = G::Base::one();
+      let mut u = G::Base::zero();
+      for bit in u_bits {
+        if bit {
+          u += mult;
+        }
+        mult = mult + mult;
+      }
+      Ok(u)
+    })?;
+
+    // The running X is two items! the running h1 and the running h2
+    let Xr0 = BigNat::alloc_from_nat(
+      cs.namespace(|| "allocate X_r[0]"),
+      || Ok(f_to_nat(&self.inputs.get()?.u2.X[0])),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    // Analyze Xr0 as limbs to use later
+    let Xr0_bn = Xr0
+      .as_limbs::<CS>()
+      .iter()
+      .enumerate()
+      .map(|(i, limb)| {
+        limb
+          .as_sapling_allocated_num(cs.namespace(|| format!("convert limb {} of X_r[0] to num", i)))
+      })
+      .collect::<Result<Vec<AllocatedNum<G::Base>>, _>>()?;
+
+    let Xr1 = BigNat::alloc_from_nat(
+      cs.namespace(|| "allocate X_r[1]"),
+      || Ok(f_to_nat(&self.inputs.get()?.u2.X[1])),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    // Analyze Xr1 as limbs to use later
+    let Xr1_bn = Xr1
+      .as_limbs::<CS>()
+      .iter()
+      .enumerate()
+      .map(|(i, limb)| {
+        limb
+          .as_sapling_allocated_num(cs.namespace(|| format!("convert limb {} of X_r[1] to num", i)))
+      })
+      .collect::<Result<Vec<AllocatedNum<G::Base>>, _>>()?;
+
+    /***********************************************************************/
+    // Allocate i
+    /***********************************************************************/
+
+    let i = AllocatedNum::alloc(cs.namespace(|| "i"), || Ok(self.inputs.get()?.i))?;
+
+    /***********************************************************************/
+    // Allocate T
+    /***********************************************************************/
+
+    // T = (x, y, infinity)
+    let T_x = AllocatedNum::alloc(cs.namespace(|| "T.x"), || {
+      Ok(self.inputs.get()?.T.comm.to_coordinates().0)
+    })?;
+    let T_y = AllocatedNum::alloc(cs.namespace(|| "T.y"), || {
+      Ok(self.inputs.get()?.T.comm.to_coordinates().1)
+    })?;
+    let T_inf = AllocatedNum::alloc(cs.namespace(|| "T.inf"), || {
+      let infty = self.inputs.get()?.T.comm.to_coordinates().2;
+      if infty {
+        Ok(G::Base::one())
+      } else {
+        Ok(G::Base::zero())
+      }
+    })?;
+
+    let T = AllocatedPoint::new(T_x.clone(), T_y.clone(), T_inf.clone());
+    let _ = T.check_is_infinity(cs.namespace(|| "T check is_infinity"));
+
+    /***********************************************************************/
+    // Allocate params
+    /***********************************************************************/
+
+    let params = AllocatedNum::alloc(cs.namespace(|| "params"), || Ok(self.inputs.get()?.params))?;
+
+    /***********************************************************************/
+    // Allocate W
+    /***********************************************************************/
+
+    // T = (x, y, infinity)
+    let W_x = AllocatedNum::alloc(cs.namespace(|| "W.x"), || {
+      Ok(self.inputs.get()?.w.comm.to_coordinates().0)
+    })?;
+    let W_y = AllocatedNum::alloc(cs.namespace(|| "W.y"), || {
+      Ok(self.inputs.get()?.w.comm.to_coordinates().1)
+    })?;
+    let W_inf = AllocatedNum::alloc(cs.namespace(|| "w.inf"), || {
+      let infty = self.inputs.get()?.w.comm.to_coordinates().2;
+      if infty {
+        Ok(G::Base::one())
+      } else {
+        Ok(G::Base::zero())
+      }
+    })?;
+
+    let W = AllocatedPoint::new(W_x.clone(), W_y.clone(), W_inf.clone());
+    let _ = W.check_is_infinity(cs.namespace(|| "W check is_infinity"));
+
+    /***********************************************************************/
+    // U2' = default if i == 0, otherwise NIFS.V(pp, u_new, U, T)
+    /***********************************************************************/
+
+    //Allocate 0 and 1
+    let zero = alloc_zero(cs.namespace(|| "zero"))?;
+    // Hack: We just do this because the number of inputs must be even!!
+    zero.inputize(cs.namespace(|| "allocate zero as input"))?;
+    let one = alloc_one(cs.namespace(|| "one"))?;
+
+    // Compute default values of U2':
+    let zero_commitment = AllocatedPoint::new(zero.clone(), zero.clone(), one);
+
+    //W_default and E_default are a commitment to zero
+    let W_default = zero_commitment.clone();
+    let E_default = zero_commitment;
+
+    // u_default = 0
+    let u_default = zero.clone();
+
+    // X_default = 0
+    let X0_default = BigNat::alloc_from_nat(
+      cs.namespace(|| "allocate x_default[0]"),
+      || Ok(f_to_nat(&G::Scalar::zero())),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    let X1_default = BigNat::alloc_from_nat(
+      cs.namespace(|| "allocate x_default[1]"),
+      || Ok(f_to_nat(&G::Scalar::zero())),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    // Compute r:
+
+    let mut ro: PoseidonROGadget<G::Base> = PoseidonROGadget::new(self.poseidon_constants.clone());
+
+    ro.absorb(h1.clone());
+    // absorb each of the limbs of h2
+    // TODO: Check if it is more efficient to treat h2 as allocNum
+    for (i, limb) in h2_bn.as_limbs::<CS>().iter().enumerate() {
+      let limb_num = limb
+        .as_sapling_allocated_num(cs.namespace(|| format!("convert limb {} of h2 to num", i)))?;
+      ro.absorb(limb_num);
+    }
+    ro.absorb(W_x);
+    ro.absorb(W_y);
+    ro.absorb(W_inf);
+    ro.absorb(T_x);
+    ro.absorb(T_y);
+    ro.absorb(T_inf);
+    // absorb each of the limbs of X_r[0]
+    for limb in Xr0_bn.clone().into_iter() {
+      ro.absorb(limb);
+    }
+
+    // absorb each of the limbs of X_r[1]
+    for limb in Xr1_bn.clone().into_iter() {
+      ro.absorb(limb);
+    }
+
+    let r_bits = ro.get_challenge(cs.namespace(|| "r bits"))?;
+    let r = le_bits_to_num(cs.namespace(|| "r"), r_bits.clone())?;
+
+    // W_fold = W_r + r * W
+    let rW = W.scalar_mul(cs.namespace(|| "r * W"), r_bits.clone())?;
+    let W_fold = W_r.add(cs.namespace(|| "W_r + r * W"), &rW)?;
+
+    // E_fold = E_r + r * T
+    let rT = T.scalar_mul(cs.namespace(|| "r * T"), r_bits)?;
+    let E_fold = E_r.add(cs.namespace(|| "E_r + r * T"), &rT)?;
+
+    // u_fold = u_r + r
+    let u_fold = AllocatedNum::alloc(cs.namespace(|| "u_fold"), || {
+      Ok(*u_r.get_value().get()? + r.get_value().get()?)
+    })?;
+    cs.enforce(
+      || "Check u_fold",
+      |lc| lc,
+      |lc| lc,
+      |lc| lc + u_fold.get_variable() - u_r.get_variable() - r.get_variable(),
+    );
+
+    // Fold the IO:
+    // Analyze r into limbs
+    let r_bn = BigNat::from_num(
+      cs.namespace(|| "allocate r_bn"),
+      Num::from(r.clone()),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    // Allocate the order of the non-native field as a constant
+    let m_bn = alloc_bignat_constant(
+      cs.namespace(|| "alloc m"),
+      &G::Scalar::get_order(),
+      self.params.limb_width,
+      self.params.n_limbs,
+    )?;
+
+    // First the fold h1 with X_r[0];
+    let (_, r_0) = h1_bn.mult_mod(cs.namespace(|| "r*h1"), &r_bn, &m_bn)?;
+    // add X_r[0]
+    let r_new_0 = Xr0.add::<CS>(&r_0)?;
+    // Now reduce
+    let Xr0_fold = r_new_0.red_mod(cs.namespace(|| "reduce folded X_r[0]"), &m_bn)?;
+
+    // First the fold h2 with X_r[1];
+    let (_, r_1) = h2_bn.mult_mod(cs.namespace(|| "r*h2"), &r_bn, &m_bn)?;
+    // add X_r[1]
+    let r_new_1 = Xr1.add::<CS>(&r_1)?;
+    // Now reduce
+    let Xr1_fold = r_new_1.red_mod(cs.namespace(|| "reduce folded X_r[1]"), &m_bn)?;
+
+    // Now select the default values if i == 0 otherwise the fold values
+    let base_case = Boolean::from(alloc_num_equals(
+      cs.namespace(|| "Check if base case"),
+      i.clone(),
+      zero,
+    )?);
+
+    let W_new = AllocatedPoint::conditionally_select(
+      cs.namespace(|| "W_new"),
+      &W_default,
+      &W_fold,
+      &base_case,
+    )?;
+
+    let E_new = AllocatedPoint::conditionally_select(
+      cs.namespace(|| "E_new"),
+      &E_default,
+      &E_fold,
+      &base_case,
+    )?;
+
+    let u_new = conditionally_select(cs.namespace(|| "u_new"), &u_default, &u_fold, &base_case)?;
+
+    let Xr0_new = conditionally_select_bignat(
+      cs.namespace(|| "X_r_new[0]"),
+      &X0_default,
+      &Xr0_fold,
+      &base_case,
+    )?;
+
+    // Analyze Xr0_new as limbs to use later
+    let Xr0_new_bn = Xr0_new
+      .as_limbs::<CS>()
+      .iter()
+      .enumerate()
+      .map(|(i, limb)| {
+        limb.as_sapling_allocated_num(
+          cs.namespace(|| format!("convert limb {} of X_r_new[0] to num", i)),
+        )
+      })
+      .collect::<Result<Vec<AllocatedNum<G::Base>>, _>>()?;
+
+    let Xr1_new = conditionally_select_bignat(
+      cs.namespace(|| "X_r_new[1]"),
+      &X1_default,
+      &Xr1_fold,
+      &base_case,
+    )?;
+
+    // Analyze Xr1_new as limbs to use later
+    let Xr1_new_bn = Xr1_new
+      .as_limbs::<CS>()
+      .iter()
+      .enumerate()
+      .map(|(i, limb)| {
+        limb.as_sapling_allocated_num(
+          cs.namespace(|| format!("convert limb {} of X_r_new[1] to num", i)),
+        )
+      })
+      .collect::<Result<Vec<AllocatedNum<G::Base>>, _>>()?;
+
+    /***********************************************************************/
+    //Compute i + 1
+    /***********************************************************************/
+
+    let next_i = AllocatedNum::alloc(cs.namespace(|| "i + 1"), || {
+      Ok(*i.get_value().get()? + G::Base::one())
+    })?;
+
+    cs.enforce(
+      || "check i + 1",
+      |lc| lc,
+      |lc| lc,
+      |lc| lc + next_i.get_variable() - CS::one() - i.get_variable(),
+    );
+
+    if self.inner_circuit.is_some() {
+      /***********************************************************************/
+      //Allocate z0
+      /***********************************************************************/
+
+      let z_0 = AllocatedNum::alloc(cs.namespace(|| "z0"), || Ok(self.inputs.get()?.z0))?;
+
+      /***********************************************************************/
+      //Allocate zi
+      /***********************************************************************/
+
+      let z_i = AllocatedNum::alloc(cs.namespace(|| "zi"), || Ok(self.inputs.get()?.zi))?;
+
+      /***********************************************************************/
+      //Check that if i == 0, z0 = zi, that is (i == 0) AND (z0 != zi) = false
+      /***********************************************************************/
+
+      let z0_is_zi = Boolean::from(alloc_num_equals(
+        cs.namespace(|| "z0 = zi"),
+        z_0.clone(),
+        z_i.clone(),
+      )?);
+
+      cs.enforce(
+        || "i == 0 and z0 != zi = false",
+        |_| base_case.lc(CS::one(), G::Base::one()),
+        |_| z0_is_zi.not().lc(CS::one(), G::Base::one()),
+        |lc| lc,
+      );
+
+      /***********************************************************************/
+      // Check that h1 = Hash(params, u2,i,z0,zi)
+      /***********************************************************************/
+
+      let mut h1_hash: PoseidonROGadget<G::Base> =
+        PoseidonROGadget::new(self.poseidon_constants.clone());
+
+      h1_hash.absorb(params.clone());
+      h1_hash.absorb(W_r_x);
+      h1_hash.absorb(W_r_y);
+      h1_hash.absorb(W_r_inf);
+      h1_hash.absorb(E_r_x);
+      h1_hash.absorb(E_r_y);
+      h1_hash.absorb(E_r_inf);
+      h1_hash.absorb(u_r.clone());
+
+      // absorb each of the limbs of X_r[0]
+      for limb in Xr0_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      // absorb each of the limbs of X_r[1]
+      for limb in Xr1_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      h1_hash.absorb(i.clone());
+      h1_hash.absorb(z_0.clone());
+      h1_hash.absorb(z_i.clone());
+
+      let hash_bits = h1_hash.get_challenge(cs.namespace(|| "Input hash"))?;
+      let hash = le_bits_to_num(cs.namespace(|| "bits to hash"), hash_bits)?;
+
+      cs.enforce(
+        || "check h1",
+        |lc| lc,
+        |lc| lc,
+        |lc| lc + h1.get_variable() - hash.get_variable(),
+      );
+
+      /***********************************************************************/
+      // Compute z_{i+1}
+      /***********************************************************************/
+
+      let z_next = self
+        .inner_circuit
+        .unwrap()
+        .synthesize(&mut cs.namespace(|| "F"), z_i)?;
+
+      /***********************************************************************/
+      // Compute the new hash H(params, u2_new, i+1, z0, z_{i+1})
+      /***********************************************************************/
+
+      h1_hash.flush_state();
+      h1_hash.absorb(params);
+      h1_hash.absorb(W_new.x.clone());
+      h1_hash.absorb(W_new.y.clone());
+      h1_hash.absorb(W_new.is_infinity);
+      h1_hash.absorb(E_new.x.clone());
+      h1_hash.absorb(E_new.y.clone());
+      h1_hash.absorb(E_new.is_infinity);
+      h1_hash.absorb(u_new);
+
+      // absorb each of the limbs of X_r_new[0]
+      for limb in Xr0_new_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      // absorb each of the limbs of X_r_new[1]
+      for limb in Xr1_new_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      h1_hash.absorb(next_i.clone());
+      h1_hash.absorb(z_0);
+      h1_hash.absorb(z_next);
+      let h1_new_bits = h1_hash.get_challenge(cs.namespace(|| "h1_new bits"))?;
+      let h1_new = le_bits_to_num(cs.namespace(|| "h1_new"), h1_new_bits)?;
+      let _ = h1_new.inputize(cs.namespace(|| "output h1_new"))?;
+    } else {
+      /***********************************************************************/
+      // Check that h1 = Hash(params, u2, i)
+      /***********************************************************************/
+
+      let mut h1_hash: PoseidonROGadget<G::Base> = PoseidonROGadget::new(self.poseidon_constants);
+
+      h1_hash.absorb(params.clone());
+      h1_hash.absorb(W_r_x);
+      h1_hash.absorb(W_r_y);
+      h1_hash.absorb(W_r_inf);
+      h1_hash.absorb(E_r_x);
+      h1_hash.absorb(E_r_y);
+      h1_hash.absorb(E_r_inf);
+      h1_hash.absorb(u_r);
+      h1_hash.absorb(i.clone());
+
+      //absorb each of the limbs of X_r[0]
+      for limb in Xr0_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      //absorb each of the limbs of X_r[1]
+      for limb in Xr1_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      let hash_bits = h1_hash.get_challenge(cs.namespace(|| "Input hash"))?;
+      let hash = le_bits_to_num(cs.namespace(|| "bits to hash"), hash_bits)?;
+
+      cs.enforce(
+        || "check h1",
+        |lc| lc,
+        |lc| lc,
+        |lc| lc + h1.get_variable() - hash.get_variable(),
+      );
+
+      /***********************************************************************/
+      // Compute the new hash H(params, u2')
+      /***********************************************************************/
+
+      h1_hash.flush_state();
+      h1_hash.absorb(params);
+      h1_hash.absorb(W_new.x.clone());
+      h1_hash.absorb(W_new.y.clone());
+      h1_hash.absorb(W_new.is_infinity);
+      h1_hash.absorb(E_new.x.clone());
+      h1_hash.absorb(E_new.y.clone());
+      h1_hash.absorb(E_new.is_infinity);
+      h1_hash.absorb(u_new);
+      h1_hash.absorb(next_i.clone());
+
+      // absorb each of the limbs of X_r_new[0]
+      for limb in Xr0_new_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      // absorb each of the limbs of X_r_new[1]
+      for limb in Xr1_new_bn.into_iter() {
+        h1_hash.absorb(limb);
+      }
+
+      let h1_new_bits = h1_hash.get_challenge(cs.namespace(|| "h1_new bits"))?;
+      let h1_new = le_bits_to_num(cs.namespace(|| "h1_new"), h1_new_bits)?;
+      let _ = h1_new.inputize(cs.namespace(|| "output h1_new"))?;
+    }
+
+    Ok(())
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+  use crate::bellperson::shape_cs::ShapeCS;
+  use crate::bellperson::solver::SatisfyingAssignment;
+  type G1 = pasta_curves::pallas::Point;
+  type G2 = pasta_curves::vesta::Point;
+  use crate::bellperson::r1cs::{NovaShape, NovaWitness};
+  use crate::commitments::CommitTrait;
+  use std::marker::PhantomData;
+
+  struct TestCircuit<F>
+  where
+    F: PrimeField + ff::PrimeField,
+  {
+    _p: PhantomData<F>,
+  }
+
+  impl<F> InnerCircuit<F> for TestCircuit<F>
+  where
+    F: PrimeField + ff::PrimeField,
+  {
+    fn synthesize<CS: ConstraintSystem<F>>(
+      &self,
+      _cs: &mut CS,
+      z: AllocatedNum<F>,
+    ) -> Result<AllocatedNum<F>, SynthesisError> {
+      Ok(z)
+    }
+  }
+
+  #[test]
+  fn test_verification_circuit() {
+    // We experiment with 8 limbs of 32 bits each
+    let params = VerificationCircuitParams::new(32, 8);
+    // The first circuit that verifies G2
+    let poseidon_constants1: NovaPoseidonConstants<<G2 as Group>::Base> =
+      NovaPoseidonConstants::new();
+    let circuit1: VerificationCircuit<G2, TestCircuit<<G2 as Group>::Base>> =
+      VerificationCircuit::new(
+        params.clone(),
+        None,
+        Some(TestCircuit {
+          _p: Default::default(),
+        }),
+        poseidon_constants1.clone(),
+      );
+    // First create the shape
+    let mut cs: ShapeCS<G1> = ShapeCS::new();
+    let _ = circuit1.synthesize(&mut cs);
+    let shape1 = cs.r1cs_shape();
+    let gens1 = cs.r1cs_gens();
+    println!(
+      "Circuit1 -> Number of constraints: {}",
+      cs.num_constraints()
+    );
+
+    // The second circuit that verifies G1
+    let poseidon_constants2: NovaPoseidonConstants<<G1 as Group>::Base> =
+      NovaPoseidonConstants::new();
+    let circuit2: VerificationCircuit<G1, TestCircuit<<G1 as Group>::Base>> =
+      VerificationCircuit::new(params.clone(), None, None, poseidon_constants2);
+    // First create the shape
+    let mut cs: ShapeCS<G2> = ShapeCS::new();
+    let _ = circuit2.synthesize(&mut cs);
+    let shape2 = cs.r1cs_shape();
+    let gens2 = cs.r1cs_gens();
+    println!(
+      "Circuit2 -> Number of constraints: {}",
+      cs.num_constraints()
+    );
+
+    //TODO: We need to hardwire default hash or give it as input
+    let default_hash = <<G2 as Group>::Base as ff::PrimeField>::from_str_vartime(
+      "332553638888022689042501686561503049809",
+    )
+    .unwrap();
+    let T = vec![<G2 as Group>::Scalar::zero()].commit(&gens2.gens_E);
+    let w = vec![<G2 as Group>::Scalar::zero()].commit(&gens2.gens_E);
+    // Now get an assignment
+    let mut cs: SatisfyingAssignment<G1> = SatisfyingAssignment::new();
+    let inputs: VerificationCircuitInputs<G2> = VerificationCircuitInputs::new(
+      default_hash,
+      RelaxedR1CSInstance::default(&gens2, &shape2),
+      <<G2 as Group>::Base as PrimeField>::zero(), // TODO: provide real inputs
+      <<G2 as Group>::Base as PrimeField>::zero(), // TODO: provide real inputs
+      <<G2 as Group>::Base as PrimeField>::zero(), // TODO: provide real inputs
+      <<G2 as Group>::Scalar as PrimeField>::zero(), // TODO: provide real inputs
+      <<G2 as Group>::Base as PrimeField>::zero(), // TODO: provide real inputs
+      T,                                           // TODO: provide real inputs
+      w,
+    );
+    let circuit: VerificationCircuit<G2, TestCircuit<<G2 as Group>::Base>> =
+      VerificationCircuit::new(
+        params,
+        Some(inputs),
+        Some(TestCircuit {
+          _p: Default::default(),
+        }),
+        poseidon_constants1,
+      );
+    let _ = circuit.synthesize(&mut cs);
+    let (inst, witness) = cs.r1cs_instance_and_witness(&shape1, &gens1).unwrap();
+
+    // Make sure that this is satisfiable
+    assert!(shape1.is_sat(&gens1, &inst, &witness).is_ok());
+  }
+}
diff --git a/src/commitments.rs b/src/commitments.rs
index 8d64f2b..9ac322c 100644
--- a/src/commitments.rs
+++ b/src/commitments.rs
@@ -14,7 +14,7 @@ pub struct CommitGens<G: Group> {
 
 #[derive(Clone, Copy, Debug, PartialEq, Eq)]
 pub struct Commitment<G: Group> {
-  comm: G,
+  pub(crate) comm: G,
 }
 
 #[derive(Clone, Debug, PartialEq, Eq)]
@@ -64,9 +64,9 @@ pub trait CommitTrait<G: Group> {
 
 impl<G: Group> CommitTrait<G> for [G::Scalar] {
   fn commit(&self, gens: &CommitGens<G>) -> Commitment<G> {
-    assert_eq!(gens.gens.len(), self.len());
+    assert!(gens.gens.len() >= self.len());
     Commitment {
-      comm: G::vartime_multiscalar_mul(self, &gens.gens),
+      comm: G::vartime_multiscalar_mul(self, &gens.gens[..self.len()]),
     }
   }
 }
diff --git a/src/gadgets/ecc.rs b/src/gadgets/ecc.rs
new file mode 100644
index 0000000..613442b
--- /dev/null
+++ b/src/gadgets/ecc.rs
@@ -0,0 +1,219 @@
+#![allow(non_snake_case)]
+use ff::{PrimeField, PrimeFieldBits};
+use rand::rngs::OsRng;
+use std::marker::PhantomData;
+
+#[derive(Debug, Clone)]
+pub struct Point<Fp, Fq>
+where
+  Fp: PrimeField,
+  Fq: PrimeField + PrimeFieldBits,
+{
+  pub(crate) x: Fp, //TODO: Make this not public
+  pub(crate) y: Fp,
+  is_infinity: bool,
+  _p: PhantomData<Fq>,
+}
+
+impl<Fp, Fq> Point<Fp, Fq>
+where
+  Fp: PrimeField,
+  Fq: PrimeField + PrimeFieldBits,
+{
+  #[allow(dead_code)]
+  pub fn new(x: Fp, y: Fp, is_infinity: bool) -> Self {
+    Self {
+      x,
+      y,
+      is_infinity,
+      _p: Default::default(),
+    }
+  }
+
+  #[allow(dead_code)]
+  pub fn random_vartime() -> Self {
+    loop {
+      let x = Fp::random(&mut OsRng);
+      let y = (x * x * x + Fp::one() + Fp::one() + Fp::one() + Fp::one() + Fp::one()).sqrt();
+      if y.is_some().unwrap_u8() == 1 {
+        return Self {
+          x,
+          y: y.unwrap(),
+          is_infinity: false,
+          _p: Default::default(),
+        };
+      }
+    }
+  }
+
+  pub fn add(&self, other: &Point<Fp, Fq>) -> Self {
+    if self.is_infinity {
+      return other.clone();
+    }
+
+    if other.is_infinity {
+      return self.clone();
+    }
+
+    let lambda = (other.y - self.y) * (other.x - self.x).invert().unwrap();
+    let x = lambda * lambda - self.x - other.x;
+    let y = lambda * (self.x - x) - self.y;
+    Self {
+      x,
+      y,
+      is_infinity: false,
+      _p: Default::default(),
+    }
+  }
+
+  pub fn double(&self) -> Self {
+    if self.is_infinity {
+      return Self {
+        x: Fp::zero(),
+        y: Fp::zero(),
+        is_infinity: true,
+        _p: Default::default(),
+      };
+    }
+
+    let lambda = (Fp::one() + Fp::one() + Fp::one())
+      * self.x
+      * self.x
+      * ((Fp::one() + Fp::one()) * self.y).invert().unwrap();
+    let x = lambda * lambda - self.x - self.x;
+    let y = lambda * (self.x - x) - self.y;
+    Self {
+      x,
+      y,
+      is_infinity: false,
+      _p: Default::default(),
+    }
+  }
+
+  #[allow(dead_code)]
+  pub fn scalar_mul_mont(&self, scalar: &Fq) -> Self {
+    let mut R0 = Self {
+      x: Fp::zero(),
+      y: Fp::zero(),
+      is_infinity: true,
+      _p: Default::default(),
+    };
+
+    let mut R1 = self.clone();
+    let bits = scalar.to_le_bits();
+    for i in (0..bits.len()).rev() {
+      if bits[i] {
+        R0 = R0.add(&R1);
+        R1 = R1.double();
+      } else {
+        R1 = R0.add(&R1);
+        R0 = R0.double();
+      }
+    }
+    R0
+  }
+
+  #[allow(dead_code)]
+  pub fn scalar_mul(&self, scalar: &Fq) -> Self {
+    let mut res = Self {
+      x: Fp::zero(),
+      y: Fp::zero(),
+      is_infinity: true,
+      _p: Default::default(),
+    };
+
+    let bits = scalar.to_le_bits();
+    for i in (0..bits.len()).rev() {
+      res = res.double();
+      if bits[i] {
+        res = self.add(&res);
+      }
+    }
+    res
+  }
+}
+
+#[cfg(test)]
+#[allow(clippy::too_many_arguments)]
+mod fp {
+  use ff::PrimeField;
+
+  #[derive(PrimeField)]
+  #[PrimeFieldModulus = "28948022309329048855892746252171976963363056481941560715954676764349967630337"]
+  #[PrimeFieldGenerator = "5"]
+  #[PrimeFieldReprEndianness = "little"]
+  pub struct Fp([u64; 4]);
+}
+
+#[cfg(test)]
+#[allow(clippy::too_many_arguments)]
+mod fq {
+  use ff::PrimeField;
+
+  #[derive(PrimeField)]
+  #[PrimeFieldModulus = "28948022309329048855892746252171976963363056481941647379679742748393362948097"]
+  #[PrimeFieldGenerator = "5"]
+  #[PrimeFieldReprEndianness = "little"]
+  pub struct Fq([u64; 4]);
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+  use super::{fp::Fp, fq::Fq};
+  use ff::Field;
+  use pasta_curves::arithmetic::CurveAffine;
+  use pasta_curves::group::Curve;
+  use pasta_curves::EpAffine;
+  use std::ops::Mul;
+
+  #[test]
+  fn test_ecc_ops() {
+    // perform some curve arithmetic
+    let a = Point::<Fp, Fq>::random_vartime();
+    let b = Point::<Fp, Fq>::random_vartime();
+    let c = a.add(&b);
+    let d = a.double();
+    let s = Fq::random(&mut OsRng);
+    let e = a.scalar_mul(&s);
+
+    // perform the same computation by translating to pasta_curve types
+    let a_pasta = EpAffine::from_xy(
+      pasta_curves::Fp::from_repr(a.x.to_repr().0).unwrap(),
+      pasta_curves::Fp::from_repr(a.y.to_repr().0).unwrap(),
+    )
+    .unwrap();
+    let b_pasta = EpAffine::from_xy(
+      pasta_curves::Fp::from_repr(b.x.to_repr().0).unwrap(),
+      pasta_curves::Fp::from_repr(b.y.to_repr().0).unwrap(),
+    )
+    .unwrap();
+    let c_pasta = (a_pasta + b_pasta).to_affine();
+    let d_pasta = (a_pasta + a_pasta).to_affine();
+    let e_pasta = a_pasta
+      .mul(pasta_curves::Fq::from_repr(s.to_repr().0).unwrap())
+      .to_affine();
+
+    // transform c, d, and e into pasta_curve types
+    let c_pasta_2 = EpAffine::from_xy(
+      pasta_curves::Fp::from_repr(c.x.to_repr().0).unwrap(),
+      pasta_curves::Fp::from_repr(c.y.to_repr().0).unwrap(),
+    )
+    .unwrap();
+    let d_pasta_2 = EpAffine::from_xy(
+      pasta_curves::Fp::from_repr(d.x.to_repr().0).unwrap(),
+      pasta_curves::Fp::from_repr(d.y.to_repr().0).unwrap(),
+    )
+    .unwrap();
+    let e_pasta_2 = EpAffine::from_xy(
+      pasta_curves::Fp::from_repr(e.x.to_repr().0).unwrap(),
+      pasta_curves::Fp::from_repr(e.y.to_repr().0).unwrap(),
+    )
+    .unwrap();
+
+    // check that we have the same outputs
+    assert_eq!(c_pasta, c_pasta_2);
+    assert_eq!(d_pasta, d_pasta_2);
+    assert_eq!(e_pasta, e_pasta_2);
+  }
+}
diff --git a/src/gadgets/ecc_circuit.rs b/src/gadgets/ecc_circuit.rs
new file mode 100644
index 0000000..3d13417
--- /dev/null
+++ b/src/gadgets/ecc_circuit.rs
@@ -0,0 +1,552 @@
+#![allow(non_snake_case)]
+use crate::gadgets::utils::{
+  alloc_one, alloc_zero, conditionally_select, conditionally_select2, select_one_or, select_zero_or,
+};
+use bellperson::{
+  gadgets::{
+    boolean::{AllocatedBit, Boolean},
+    num::AllocatedNum,
+    Assignment,
+  },
+  ConstraintSystem, SynthesisError,
+};
+use ff::PrimeField;
+use rand::rngs::OsRng;
+
+#[derive(Clone)]
+pub struct AllocatedPoint<Fp>
+where
+  Fp: PrimeField,
+{
+  pub(crate) x: AllocatedNum<Fp>,
+  pub(crate) y: AllocatedNum<Fp>,
+  pub(crate) is_infinity: AllocatedNum<Fp>,
+}
+
+impl<Fp> AllocatedPoint<Fp>
+where
+  Fp: PrimeField,
+{
+  // Creates a new allocated point from allocated nums.
+  pub fn new(x: AllocatedNum<Fp>, y: AllocatedNum<Fp>, is_infinity: AllocatedNum<Fp>) -> Self {
+    Self { x, y, is_infinity }
+  }
+
+  // Check that is infinity is 0/1
+  #[allow(dead_code)]
+  pub fn check_is_infinity<CS: ConstraintSystem<Fp>>(
+    &self,
+    mut cs: CS,
+  ) -> Result<(), SynthesisError> {
+    // Check that is_infinity * ( 1 - is_infinity ) = 0
+    cs.enforce(
+      || "is_infinity is bit",
+      |lc| lc + self.is_infinity.get_variable(),
+      |lc| lc + CS::one() - self.is_infinity.get_variable(),
+      |lc| lc,
+    );
+    Ok(())
+  }
+
+  #[allow(dead_code)]
+  // Allocate a random point. Only used for testing
+  pub fn random_vartime<CS: ConstraintSystem<Fp>>(mut cs: CS) -> Result<Self, SynthesisError> {
+    loop {
+      let x = Fp::random(&mut OsRng);
+      let y = (x * x * x + Fp::one() + Fp::one() + Fp::one() + Fp::one() + Fp::one()).sqrt();
+      if y.is_some().unwrap_u8() == 1 {
+        let x_alloc = AllocatedNum::alloc(cs.namespace(|| "x"), || Ok(x))?;
+        let y_alloc = AllocatedNum::alloc(cs.namespace(|| "y"), || Ok(y.unwrap()))?;
+        let is_infinity = alloc_zero(cs.namespace(|| "Is Infinity"))?;
+        return Ok(Self::new(x_alloc, y_alloc, is_infinity));
+      }
+    }
+  }
+
+  // Make the point io
+  #[allow(dead_code)]
+  pub fn inputize<CS: ConstraintSystem<Fp>>(&self, mut cs: CS) -> Result<(), SynthesisError> {
+    let _ = self.x.inputize(cs.namespace(|| "Input point.x"));
+    let _ = self.y.inputize(cs.namespace(|| "Input point.y"));
+    let _ = self
+      .is_infinity
+      .inputize(cs.namespace(|| "Input point.is_infinity"));
+    Ok(())
+  }
+
+  // Adds other point to this point and returns the result
+  // Assumes that both other.is_infinity and this.is_infinty are bits
+  pub fn add<CS: ConstraintSystem<Fp>>(
+    &self,
+    mut cs: CS,
+    other: &AllocatedPoint<Fp>,
+  ) -> Result<Self, SynthesisError> {
+    // Allocate the boolean variables that check if either of the points is infinity
+
+    //************************************************************************/
+    // lambda = (other.y - self.y) * (other.x - self.x).invert().unwrap();
+    //************************************************************************/
+    // First compute (other.x - self.x).inverse()
+    // If either self or other are 1 then compute bogus values
+
+    // x_diff = other != inf && self != inf ? (other.x - self.x) : 1
+    let x_diff_actual = AllocatedNum::alloc(cs.namespace(|| "actual x diff"), || {
+      Ok(*other.x.get_value().get()? - *self.x.get_value().get()?)
+    })?;
+    cs.enforce(
+      || "actual x_diff is correct",
+      |lc| lc + other.x.get_variable() - self.x.get_variable(),
+      |lc| lc + CS::one(),
+      |lc| lc + x_diff_actual.get_variable(),
+    );
+
+    // Compute self.is_infinity OR other.is_infinity
+    let at_least_one_inf = AllocatedNum::alloc(cs.namespace(|| "at least one inf"), || {
+      Ok(*self.is_infinity.get_value().get()? * *other.is_infinity.get_value().get()?)
+    })?;
+    cs.enforce(
+      || "at least one inf = self.is_infinity * other.is_infinity",
+      |lc| lc + self.is_infinity.get_variable(),
+      |lc| lc + other.is_infinity.get_variable(),
+      |lc| lc + at_least_one_inf.get_variable(),
+    );
+
+    // x_diff = 1 if either self.is_infinity or other.is_infinity else x_diff_actual
+    let x_diff = select_one_or(
+      cs.namespace(|| "Compute x_diff"),
+      &x_diff_actual,
+      &at_least_one_inf,
+    )?;
+
+    let x_diff_inv = AllocatedNum::alloc(cs.namespace(|| "x diff inverse"), || {
+      if *at_least_one_inf.get_value().get()? == Fp::one() {
+        // Set to default
+        Ok(Fp::one())
+      } else {
+        // Set to the actual inverse
+        let inv = (*other.x.get_value().get()? - *self.x.get_value().get()?).invert();
+        if inv.is_some().unwrap_u8() == 1 {
+          Ok(inv.unwrap())
+        } else {
+          Err(SynthesisError::DivisionByZero)
+        }
+      }
+    })?;
+
+    cs.enforce(
+      || "Check inverse",
+      |lc| lc + x_diff.get_variable(),
+      |lc| lc + x_diff_inv.get_variable(),
+      |lc| lc + CS::one(),
+    );
+
+    let lambda = AllocatedNum::alloc(cs.namespace(|| "lambda"), || {
+      Ok(
+        (*other.y.get_value().get()? - *self.y.get_value().get()?)
+          * x_diff_inv.get_value().get()?,
+      )
+    })?;
+    cs.enforce(
+      || "Check that lambda is correct",
+      |lc| lc + other.y.get_variable() - self.y.get_variable(),
+      |lc| lc + x_diff_inv.get_variable(),
+      |lc| lc + lambda.get_variable(),
+    );
+
+    //************************************************************************/
+    // x = lambda * lambda - self.x - other.x;
+    //************************************************************************/
+    let x = AllocatedNum::alloc(cs.namespace(|| "x"), || {
+      Ok(
+        *lambda.get_value().get()? * lambda.get_value().get()?
+          - *self.x.get_value().get()?
+          - *other.x.get_value().get()?,
+      )
+    })?;
+    cs.enforce(
+      || "check that x is correct",
+      |lc| lc + lambda.get_variable(),
+      |lc| lc + lambda.get_variable(),
+      |lc| lc + x.get_variable() + self.x.get_variable() + other.x.get_variable(),
+    );
+
+    //************************************************************************/
+    // y = lambda * (self.x - x) - self.y;
+    //************************************************************************/
+    let y = AllocatedNum::alloc(cs.namespace(|| "y"), || {
+      Ok(
+        *lambda.get_value().get()? * (*self.x.get_value().get()? - *x.get_value().get()?)
+          - *self.y.get_value().get()?,
+      )
+    })?;
+
+    cs.enforce(
+      || "Check that y is correct",
+      |lc| lc + lambda.get_variable(),
+      |lc| lc + self.x.get_variable() - x.get_variable(),
+      |lc| lc + y.get_variable() + self.y.get_variable(),
+    );
+
+    let is_infinity = AllocatedNum::alloc(cs.namespace(|| "is infinity"), || Ok(Fp::zero()))?;
+
+    //************************************************************************/
+    // We only return the computed x, y if neither of the points is infinity.
+    // if self.is_infinity return other.clone()
+    // elif other.is_infinity return self.clone()
+    // Otherwise return the computed points.
+    //************************************************************************/
+    // Now compute the output x
+    let inner_x = conditionally_select2(
+      cs.namespace(|| "final x: inner if"),
+      &self.x,
+      &x,
+      &other.is_infinity,
+    )?;
+    let final_x = conditionally_select2(
+      cs.namespace(|| "final x: outer if"),
+      &other.x,
+      &inner_x,
+      &self.is_infinity,
+    )?;
+
+    // The output y
+    let inner_y = conditionally_select2(
+      cs.namespace(|| "final y: inner if"),
+      &self.y,
+      &y,
+      &other.is_infinity,
+    )?;
+    let final_y = conditionally_select2(
+      cs.namespace(|| "final y: outer if"),
+      &other.y,
+      &inner_y,
+      &self.is_infinity,
+    )?;
+
+    // The output is_infinity
+    let inner_is_infinity = conditionally_select2(
+      cs.namespace(|| "final is infinity: inner if"),
+      &self.is_infinity,
+      &is_infinity,
+      &other.is_infinity,
+    )?;
+    let final_is_infinity = conditionally_select2(
+      cs.namespace(|| "final is infinity: outer if"),
+      &other.is_infinity,
+      &inner_is_infinity,
+      &self.is_infinity,
+    )?;
+    Ok(Self::new(final_x, final_y, final_is_infinity))
+  }
+
+  pub fn double<CS: ConstraintSystem<Fp>>(&self, mut cs: CS) -> Result<Self, SynthesisError> {
+    //*************************************************************/
+    // lambda = (Fp::one() + Fp::one() + Fp::one())
+    //  * self.x
+    //  * self.x
+    //  * ((Fp::one() + Fp::one()) * self.y).invert().unwrap();
+    /*************************************************************/
+
+    // Compute tmp = (Fp::one() + Fp::one())* self.y ? self != inf : 1
+    let tmp_actual = AllocatedNum::alloc(cs.namespace(|| "tmp_actual"), || {
+      Ok(*self.y.get_value().get()? + *self.y.get_value().get()?)
+    })?;
+    cs.enforce(
+      || "check tmp_actual",
+      |lc| lc + CS::one() + CS::one(),
+      |lc| lc + self.y.get_variable(),
+      |lc| lc + tmp_actual.get_variable(),
+    );
+
+    let tmp = select_one_or(cs.namespace(|| "tmp"), &tmp_actual, &self.is_infinity)?;
+
+    // Compute inv = tmp.invert
+    let tmp_inv = AllocatedNum::alloc(cs.namespace(|| "tmp inverse"), || {
+      if *self.is_infinity.get_value().get()? == Fp::one() {
+        // Return default value 1
+        Ok(Fp::one())
+      } else {
+        // Return the actual inverse
+        let inv = (*tmp.get_value().get()?).invert();
+        if inv.is_some().unwrap_u8() == 1 {
+          Ok(inv.unwrap())
+        } else {
+          Err(SynthesisError::DivisionByZero)
+        }
+      }
+    })?;
+    cs.enforce(
+      || "Check inverse",
+      |lc| lc + tmp.get_variable(),
+      |lc| lc + tmp_inv.get_variable(),
+      |lc| lc + CS::one(),
+    );
+
+    // Now compute lambda as (Fp::one() + Fp::one + Fp::one()) * self.x * self.x * tmp_inv
+    let prod_1 = AllocatedNum::alloc(cs.namespace(|| "alloc prod 1"), || {
+      Ok(*tmp_inv.get_value().get()? * self.x.get_value().get()?)
+    })?;
+    cs.enforce(
+      || "Check prod 1",
+      |lc| lc + self.x.get_variable(),
+      |lc| lc + tmp_inv.get_variable(),
+      |lc| lc + prod_1.get_variable(),
+    );
+
+    let prod_2 = AllocatedNum::alloc(cs.namespace(|| "alloc prod 2"), || {
+      Ok(*prod_1.get_value().get()? * self.x.get_value().get()?)
+    })?;
+    cs.enforce(
+      || "Check prod 2",
+      |lc| lc + self.x.get_variable(),
+      |lc| lc + prod_1.get_variable(),
+      |lc| lc + prod_2.get_variable(),
+    );
+
+    let lambda = AllocatedNum::alloc(cs.namespace(|| "lambda"), || {
+      Ok(*prod_2.get_value().get()? * (Fp::one() + Fp::one() + Fp::one()))
+    })?;
+    cs.enforce(
+      || "Check lambda",
+      |lc| lc + CS::one() + CS::one() + CS::one(),
+      |lc| lc + prod_2.get_variable(),
+      |lc| lc + lambda.get_variable(),
+    );
+
+    /*************************************************************/
+    //          x = lambda * lambda - self.x - self.x;
+    /*************************************************************/
+
+    let x = AllocatedNum::alloc(cs.namespace(|| "x"), || {
+      Ok(
+        ((*lambda.get_value().get()?) * (*lambda.get_value().get()?))
+          - *self.x.get_value().get()?
+          - self.x.get_value().get()?,
+      )
+    })?;
+    cs.enforce(
+      || "Check x",
+      |lc| lc + lambda.get_variable(),
+      |lc| lc + lambda.get_variable(),
+      |lc| lc + x.get_variable() + self.x.get_variable() + self.x.get_variable(),
+    );
+
+    /*************************************************************/
+    //        y = lambda * (self.x - x) - self.y;
+    /*************************************************************/
+
+    let y = AllocatedNum::alloc(cs.namespace(|| "y"), || {
+      Ok(
+        (*lambda.get_value().get()?) * (*self.x.get_value().get()? - x.get_value().get()?)
+          - self.y.get_value().get()?,
+      )
+    })?;
+    cs.enforce(
+      || "Check y",
+      |lc| lc + lambda.get_variable(),
+      |lc| lc + self.x.get_variable() - x.get_variable(),
+      |lc| lc + y.get_variable() + self.y.get_variable(),
+    );
+
+    /*************************************************************/
+    // Only return the computed x and y if the point is not infinity
+    /*************************************************************/
+
+    // x
+    let final_x = select_zero_or(cs.namespace(|| "final x"), &x, &self.is_infinity)?;
+
+    // y
+    let final_y = select_zero_or(cs.namespace(|| "final y"), &y, &self.is_infinity)?;
+
+    // is_infinity
+    let final_is_infinity = self.is_infinity.clone();
+
+    Ok(Self::new(final_x, final_y, final_is_infinity))
+  }
+
+  #[allow(dead_code)]
+  pub fn scalar_mul_mont<CS: ConstraintSystem<Fp>>(
+    &self,
+    mut cs: CS,
+    scalar: Vec<AllocatedBit>,
+  ) -> Result<Self, SynthesisError> {
+    /*************************************************************/
+    // Initialize RO = Self {
+    //  x: Fp::zero(),
+    //  y: Fp::zero(),
+    //  is_infinity: true,
+    //  _p: Default::default(),
+    //};
+    /*************************************************************/
+
+    let zero = alloc_zero(cs.namespace(|| "Allocate zero"))?;
+    let one = alloc_one(cs.namespace(|| "Allocate one"))?;
+    let mut R0 = Self::new(zero.clone(), zero, one);
+
+    /*************************************************************/
+    // Initialize R1 and the bits of the scalar
+    /*************************************************************/
+
+    let mut R1 = self.clone();
+
+    for i in (0..scalar.len()).rev() {
+      /*************************************************************/
+      //if bits[i] {
+      //  R0 = R0.add(&R1);
+      //  R1 = R1.double();
+      //} else {
+      //  R0 = R0.double();
+      //  R1 = R0.add(&R1);
+      //}
+      /*************************************************************/
+
+      let R0_and_R1 = R0.add(cs.namespace(|| format!("{}: R0 + R1", i)), &R1)?;
+      let R0_double = R0.double(cs.namespace(|| format!("{}: 2 * R0", i)))?;
+      let R1_double = R1.double(cs.namespace(|| format!("{}: 2 * R1", i)))?;
+
+      R0 = Self::conditionally_select(
+        cs.namespace(|| format!("{}: Update R0", i)),
+        &R0_and_R1,
+        &R0_double,
+        &Boolean::from(scalar[i].clone()),
+      )?;
+
+      R1 = Self::conditionally_select(
+        cs.namespace(|| format!("{}: Update R1", i)),
+        &R1_double,
+        &R0_and_R1,
+        &Boolean::from(scalar[i].clone()),
+      )?;
+    }
+    Ok(R0)
+  }
+
+  #[allow(dead_code)]
+  pub fn scalar_mul<CS: ConstraintSystem<Fp>>(
+    &self,
+    mut cs: CS,
+    scalar: Vec<AllocatedBit>,
+  ) -> Result<Self, SynthesisError> {
+    /*************************************************************/
+    // Initialize res = Self {
+    //  x: Fp::zero(),
+    //  y: Fp::zero(),
+    //  is_infinity: true,
+    //  _p: Default::default(),
+    //};
+    /*************************************************************/
+
+    let zero = alloc_zero(cs.namespace(|| "Allocate zero"))?;
+    let one = alloc_one(cs.namespace(|| "Allocate one"))?;
+    let mut res = Self::new(zero.clone(), zero, one);
+
+    for i in (0..scalar.len()).rev() {
+      /*************************************************************/
+      //  res = res.double();
+      /*************************************************************/
+
+      res = res.double(cs.namespace(|| format!("{}: double", i)))?;
+
+      /*************************************************************/
+      //  if scalar[i] {
+      //    res = self.add(&res);
+      //  }
+      /*************************************************************/
+
+      let self_and_res = self.add(cs.namespace(|| format!("{}: add", i)), &res)?;
+      res = Self::conditionally_select(
+        cs.namespace(|| format!("{}: Update res", i)),
+        &self_and_res,
+        &res,
+        &Boolean::from(scalar[i].clone()),
+      )?;
+    }
+    Ok(res)
+  }
+
+  /// If condition outputs a otherwise outputs b
+  pub fn conditionally_select<CS: ConstraintSystem<Fp>>(
+    mut cs: CS,
+    a: &Self,
+    b: &Self,
+    condition: &Boolean,
+  ) -> Result<Self, SynthesisError> {
+    let x = conditionally_select(cs.namespace(|| "select x"), &a.x, &b.x, condition)?;
+
+    let y = conditionally_select(cs.namespace(|| "select y"), &a.y, &b.y, condition)?;
+
+    let is_infinity = conditionally_select(
+      cs.namespace(|| "select is_infinity"),
+      &a.is_infinity,
+      &b.is_infinity,
+      condition,
+    )?;
+
+    Ok(Self::new(x, y, is_infinity))
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+  use crate::bellperson::shape_cs::ShapeCS;
+  use crate::bellperson::solver::SatisfyingAssignment;
+  type G = pasta_curves::pallas::Point;
+  type Fp = pasta_curves::pallas::Scalar;
+  type Fq = pasta_curves::vesta::Scalar;
+  use crate::bellperson::r1cs::{NovaShape, NovaWitness};
+  use crate::gadgets::ecc::Point;
+  use ff::PrimeFieldBits;
+
+  fn synthesize_smul<Fp, Fq, CS>(mut cs: CS) -> (AllocatedPoint<Fp>, AllocatedPoint<Fp>, Fq)
+  where
+    Fp: PrimeField,
+    Fq: PrimeField + PrimeFieldBits,
+    CS: ConstraintSystem<Fp>,
+  {
+    let a = AllocatedPoint::<Fp>::random_vartime(cs.namespace(|| "a")).unwrap();
+    let _ = a.inputize(cs.namespace(|| "inputize a")).unwrap();
+    let s = Fq::random(&mut OsRng);
+    // Allocate random bits and only keep 128 bits
+    let bits: Vec<AllocatedBit> = s
+      .to_le_bits()
+      .into_iter()
+      .enumerate()
+      .map(|(i, bit)| AllocatedBit::alloc(cs.namespace(|| format!("bit {}", i)), Some(bit)))
+      .collect::<Result<Vec<AllocatedBit>, SynthesisError>>()
+      .unwrap();
+    let e = a.scalar_mul(cs.namespace(|| "Scalar Mul"), bits).unwrap();
+    let _ = e.inputize(cs.namespace(|| "inputize e")).unwrap();
+    (a, e, s)
+  }
+
+  #[test]
+  fn test_ecc_circuit_ops() {
+    // First create the shape
+    let mut cs: ShapeCS<G> = ShapeCS::new();
+    let _ = synthesize_smul::<Fp, Fq, _>(cs.namespace(|| "synthesize"));
+    println!("Number of constraints: {}", cs.num_constraints());
+    let shape = cs.r1cs_shape();
+    let gens = cs.r1cs_gens();
+
+    // Then the satisfying assignment
+    let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+    let (a, e, s) = synthesize_smul::<Fp, Fq, _>(cs.namespace(|| "synthesize"));
+    let (inst, witness) = cs.r1cs_instance_and_witness(&shape, &gens).unwrap();
+
+    let a_p: Point<Fp, Fq> = Point::new(
+      a.x.get_value().unwrap(),
+      a.y.get_value().unwrap(),
+      a.is_infinity.get_value().unwrap() == Fp::one(),
+    );
+    let e_p: Point<Fp, Fq> = Point::new(
+      e.x.get_value().unwrap(),
+      e.y.get_value().unwrap(),
+      e.is_infinity.get_value().unwrap() == Fp::one(),
+    );
+    let e_new = a_p.scalar_mul(&s);
+    assert!(e_p.x == e_new.x && e_p.y == e_new.y);
+    // Make sure that this is satisfiable
+    assert!(shape.is_sat(&gens, &inst, &witness).is_ok());
+  }
+}
diff --git a/src/gadgets/mod.rs b/src/gadgets/mod.rs
new file mode 100644
index 0000000..18bad97
--- /dev/null
+++ b/src/gadgets/mod.rs
@@ -0,0 +1,3 @@
+mod ecc;
+pub mod ecc_circuit;
+pub mod utils;
diff --git a/src/gadgets/utils.rs b/src/gadgets/utils.rs
new file mode 100644
index 0000000..54a742c
--- /dev/null
+++ b/src/gadgets/utils.rs
@@ -0,0 +1,338 @@
+use bellperson::{
+  gadgets::{
+    boolean::{AllocatedBit, Boolean},
+    num::AllocatedNum,
+    Assignment,
+  },
+  ConstraintSystem, LinearCombination, SynthesisError,
+};
+use bellperson_nonnative::mp::bignat::{nat_to_limbs, BigNat};
+use ff::{PrimeField, PrimeFieldBits};
+use rug::Integer;
+
+/// Gets as input the little indian representation of a number and spits out the number
+#[allow(dead_code)]
+pub fn le_bits_to_num<Scalar, CS>(
+  mut cs: CS,
+  bits: Vec<AllocatedBit>,
+) -> Result<AllocatedNum<Scalar>, SynthesisError>
+where
+  Scalar: PrimeField + PrimeFieldBits,
+  CS: ConstraintSystem<Scalar>,
+{
+  // We loop over the input bits and construct the constraint
+  // and the field element that corresponds to the result
+  let mut lc = LinearCombination::zero();
+  let mut coeff = Scalar::one();
+  let mut fe = Some(Scalar::zero());
+  for bit in bits.iter() {
+    lc = lc + (coeff, bit.get_variable());
+    fe = bit.get_value().map(|val| {
+      if val {
+        fe.unwrap() + coeff
+      } else {
+        fe.unwrap()
+      }
+    });
+    coeff = coeff.double();
+  }
+  let num = AllocatedNum::alloc(cs.namespace(|| "Field element"), || {
+    fe.ok_or(SynthesisError::AssignmentMissing)
+  })?;
+  lc = lc - num.get_variable();
+  cs.enforce(|| "compute number from bits", |lc| lc, |lc| lc, |_| lc);
+  Ok(num)
+}
+
+/// Allocate a variable that is set to zero
+pub fn alloc_zero<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+) -> Result<AllocatedNum<F>, SynthesisError> {
+  let zero = AllocatedNum::alloc(cs.namespace(|| "alloc"), || Ok(F::zero()))?;
+  cs.enforce(
+    || "check zero is valid",
+    |lc| lc,
+    |lc| lc,
+    |lc| lc + zero.get_variable(),
+  );
+  Ok(zero)
+}
+
+/// Allocate a variable that is set to one
+pub fn alloc_one<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+) -> Result<AllocatedNum<F>, SynthesisError> {
+  let one = AllocatedNum::alloc(cs.namespace(|| "alloc"), || Ok(F::one()))?;
+  cs.enforce(
+    || "check one is valid",
+    |lc| lc + CS::one(),
+    |lc| lc + CS::one(),
+    |lc| lc + one.get_variable(),
+  );
+
+  Ok(one)
+}
+
+/// Allocate bignat a constant
+pub fn alloc_bignat_constant<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  val: &Integer,
+  limb_width: usize,
+  n_limbs: usize,
+) -> Result<BigNat<F>, SynthesisError> {
+  let limbs = nat_to_limbs(val, limb_width, n_limbs).unwrap();
+  let bignat = BigNat::alloc_from_limbs(
+    cs.namespace(|| "alloc bignat"),
+    || Ok(limbs.clone()),
+    None,
+    limb_width,
+    n_limbs,
+  )?;
+  // Now enforce that the limbs are all equal to the constants
+  #[allow(clippy::needless_range_loop)]
+  for i in 0..n_limbs {
+    cs.enforce(
+      || format!("check limb {}", i),
+      |lc| lc + &bignat.limbs[i],
+      |lc| lc + CS::one(),
+      |lc| lc + (limbs[i], CS::one()),
+    );
+  }
+  Ok(bignat)
+}
+
+/// Check that two numbers are equal and return a bit
+pub fn alloc_num_equals<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  a: AllocatedNum<F>,
+  b: AllocatedNum<F>,
+) -> Result<AllocatedBit, SynthesisError> {
+  // Allocate and constrain `r`: result boolean bit.
+  // It equals `true` if `a` equals `b`, `false` otherwise
+
+  let r_value = match (a.get_value(), b.get_value()) {
+    (Some(a), Some(b)) => Some(a == b),
+    _ => None,
+  };
+
+  let r = AllocatedBit::alloc(cs.namespace(|| "r"), r_value)?;
+
+  let delta = AllocatedNum::alloc(cs.namespace(|| "delta"), || {
+    let a_value = *a.get_value().get()?;
+    let b_value = *b.get_value().get()?;
+
+    let mut delta = a_value;
+    delta.sub_assign(&b_value);
+
+    Ok(delta)
+  })?;
+
+  cs.enforce(
+    || "delta = (a - b)",
+    |lc| lc + a.get_variable() - b.get_variable(),
+    |lc| lc + CS::one(),
+    |lc| lc + delta.get_variable(),
+  );
+
+  let delta_inv = AllocatedNum::alloc(cs.namespace(|| "delta_inv"), || {
+    let delta = *delta.get_value().get()?;
+
+    if delta.is_zero().unwrap_u8() == 1 {
+      Ok(F::one()) // we can return any number here, it doesn't matter
+    } else {
+      Ok(delta.invert().unwrap())
+    }
+  })?;
+
+  // Allocate `t = delta * delta_inv`
+  // If `delta` is non-zero (a != b), `t` will equal 1
+  // If `delta` is zero (a == b), `t` cannot equal 1
+
+  let t = AllocatedNum::alloc(cs.namespace(|| "t"), || {
+    let mut tmp = *delta.get_value().get()?;
+    tmp.mul_assign(&(*delta_inv.get_value().get()?));
+
+    Ok(tmp)
+  })?;
+
+  // Constrain allocation:
+  // t = (a - b) * delta_inv
+  cs.enforce(
+    || "t = (a - b) * delta_inv",
+    |lc| lc + a.get_variable() - b.get_variable(),
+    |lc| lc + delta_inv.get_variable(),
+    |lc| lc + t.get_variable(),
+  );
+
+  // Constrain:
+  // (a - b) * (t - 1) == 0
+  // This enforces that correct `delta_inv` was provided,
+  // and thus `t` is 1 if `(a - b)` is non zero (a != b )
+  cs.enforce(
+    || "(a - b) * (t - 1) == 0",
+    |lc| lc + a.get_variable() - b.get_variable(),
+    |lc| lc + t.get_variable() - CS::one(),
+    |lc| lc,
+  );
+
+  // Constrain:
+  // (a - b) * r == 0
+  // This enforces that `r` is zero if `(a - b)` is non-zero (a != b)
+  cs.enforce(
+    || "(a - b) * r == 0",
+    |lc| lc + a.get_variable() - b.get_variable(),
+    |lc| lc + r.get_variable(),
+    |lc| lc,
+  );
+
+  // Constrain:
+  // (t - 1) * (r - 1) == 0
+  // This enforces that `r` is one if `t` is not one (a == b)
+  cs.enforce(
+    || "(t - 1) * (r - 1) == 0",
+    |lc| lc + t.get_variable() - CS::one(),
+    |lc| lc + r.get_variable() - CS::one(),
+    |lc| lc,
+  );
+
+  Ok(r)
+}
+
+/// If condition return a otherwise b
+pub fn conditionally_select<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  a: &AllocatedNum<F>,
+  b: &AllocatedNum<F>,
+  condition: &Boolean,
+) -> Result<AllocatedNum<F>, SynthesisError> {
+  let c = AllocatedNum::alloc(cs.namespace(|| "conditional select result"), || {
+    if *condition.get_value().get()? {
+      Ok(*a.get_value().get()?)
+    } else {
+      Ok(*b.get_value().get()?)
+    }
+  })?;
+
+  // a * condition + b*(1-condition) = c ->
+  // a * condition - b*condition = c - b
+
+  cs.enforce(
+    || "conditional select constraint",
+    |lc| lc + a.get_variable() - b.get_variable(),
+    |_| condition.lc(CS::one(), F::one()),
+    |lc| lc + c.get_variable() - b.get_variable(),
+  );
+
+  Ok(c)
+}
+
+/// If condition return a otherwise b where a and b are BigNats
+pub fn conditionally_select_bignat<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  a: &BigNat<F>,
+  b: &BigNat<F>,
+  condition: &Boolean,
+) -> Result<BigNat<F>, SynthesisError> {
+  assert!(a.limbs.len() == b.limbs.len());
+  let c = BigNat::alloc_from_nat(
+    cs.namespace(|| "conditional select result"),
+    || {
+      if *condition.get_value().get()? {
+        Ok(a.value.get()?.clone())
+      } else {
+        Ok(b.value.get()?.clone())
+      }
+    },
+    a.params.limb_width,
+    a.params.n_limbs,
+  )?;
+
+  // a * condition + b*(1-condition) = c ->
+  // a * condition - b*condition = c - b
+  for i in 0..c.limbs.len() {
+    cs.enforce(
+      || format!("conditional select constraint {}", i),
+      |lc| lc + &a.limbs[i] - &b.limbs[i],
+      |_| condition.lc(CS::one(), F::one()),
+      |lc| lc + &c.limbs[i] - &b.limbs[i],
+    );
+  }
+  Ok(c)
+}
+
+/// Same as the above but Condition is an AllocatedNum that needs to be
+/// 0 or 1. 1 => True, 0 => False
+pub fn conditionally_select2<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  a: &AllocatedNum<F>,
+  b: &AllocatedNum<F>,
+  condition: &AllocatedNum<F>,
+) -> Result<AllocatedNum<F>, SynthesisError> {
+  let c = AllocatedNum::alloc(cs.namespace(|| "conditional select result"), || {
+    if *condition.get_value().get()? == F::one() {
+      Ok(*a.get_value().get()?)
+    } else {
+      Ok(*b.get_value().get()?)
+    }
+  })?;
+
+  // a * condition + b*(1-condition) = c ->
+  // a * condition - b*condition = c - b
+
+  cs.enforce(
+    || "conditional select constraint",
+    |lc| lc + a.get_variable() - b.get_variable(),
+    |lc| lc + condition.get_variable(),
+    |lc| lc + c.get_variable() - b.get_variable(),
+  );
+
+  Ok(c)
+}
+
+/// If condition set to 0 otherwise a
+pub fn select_zero_or<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  a: &AllocatedNum<F>,
+  condition: &AllocatedNum<F>,
+) -> Result<AllocatedNum<F>, SynthesisError> {
+  let c = AllocatedNum::alloc(cs.namespace(|| "conditional select result"), || {
+    if *condition.get_value().get()? == F::one() {
+      Ok(F::zero())
+    } else {
+      Ok(*a.get_value().get()?)
+    }
+  })?;
+
+  // a * (1 - condition) = c
+  cs.enforce(
+    || "conditional select constraint",
+    |lc| lc + a.get_variable(),
+    |lc| lc + CS::one() - condition.get_variable(),
+    |lc| lc + c.get_variable(),
+  );
+
+  Ok(c)
+}
+
+/// If condition set to 1 otherwise a
+pub fn select_one_or<F: PrimeField, CS: ConstraintSystem<F>>(
+  mut cs: CS,
+  a: &AllocatedNum<F>,
+  condition: &AllocatedNum<F>,
+) -> Result<AllocatedNum<F>, SynthesisError> {
+  let c = AllocatedNum::alloc(cs.namespace(|| "conditional select result"), || {
+    if *condition.get_value().get()? == F::one() {
+      Ok(F::one())
+    } else {
+      Ok(*a.get_value().get()?)
+    }
+  })?;
+
+  cs.enforce(
+    || "conditional select constraint",
+    |lc| lc + CS::one() - a.get_variable(),
+    |lc| lc + condition.get_variable(),
+    |lc| lc + c.get_variable() - a.get_variable(),
+  );
+  Ok(c)
+}
diff --git a/src/lib.rs b/src/lib.rs
index 185715b..d556d0a 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -6,8 +6,11 @@
 mod commitments;
 
 pub mod bellperson;
+mod circuit;
 pub mod errors;
+mod gadgets;
 pub mod pasta;
+mod poseidon;
 pub mod r1cs;
 pub mod traits;
 
diff --git a/src/pasta.rs b/src/pasta.rs
index be60449..dd5b98e 100644
--- a/src/pasta.rs
+++ b/src/pasta.rs
@@ -1,16 +1,36 @@
-//! This module implements the Nova traits for pallas::Point and pallas::Scalar.
+//! This module implements the Nova traits for pallas::Point, pallas::Scalar, vesta::Point, vesta::Scalar.
 use crate::traits::{ChallengeTrait, CompressedGroup, Group, PrimeField};
 use merlin::Transcript;
-use pasta_curves::arithmetic::{CurveExt, FieldExt, Group as Grp};
-use pasta_curves::group::GroupEncoding;
-use pasta_curves::{self, pallas, Ep, Fq};
+use pasta_curves::arithmetic::{CurveAffine, CurveExt, FieldExt, Group as Grp};
+use pasta_curves::group::{Curve, GroupEncoding};
+use pasta_curves::{self, pallas, vesta, Ep, Eq, Fp, Fq};
 use rand::{CryptoRng, RngCore};
+use rug::Integer;
 use std::borrow::Borrow;
 use std::ops::Mul;
 
+//////////////////////////////////////Pallas///////////////////////////////////////////////
+
+/// A wrapper for compressed group elements that come from the pallas curve
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+pub struct PallasCompressedElementWrapper {
+  repr: [u8; 32],
+}
+
+impl PallasCompressedElementWrapper {
+  /// Wraps repr into the wrapper
+  pub fn new(repr: [u8; 32]) -> Self {
+    Self { repr }
+  }
+}
+
+unsafe impl Send for PallasCompressedElementWrapper {}
+unsafe impl Sync for PallasCompressedElementWrapper {}
+
 impl Group for pallas::Point {
+  type Base = pallas::Base;
   type Scalar = pallas::Scalar;
-  type CompressedGroupElement = <pallas::Point as GroupEncoding>::Repr;
+  type CompressedGroupElement = PallasCompressedElementWrapper;
 
   fn vartime_multiscalar_mul<I, J>(scalars: I, points: J) -> Self
   where
@@ -29,7 +49,7 @@ impl Group for pallas::Point {
   }
 
   fn compress(&self) -> Self::CompressedGroupElement {
-    self.to_bytes()
+    PallasCompressedElementWrapper::new(self.to_bytes())
   }
 
   fn from_uniform_bytes(bytes: &[u8]) -> Option<Self> {
@@ -43,6 +63,15 @@ impl Group for pallas::Point {
       Some(hash(&arr))
     }
   }
+
+  fn to_coordinates(&self) -> (Self::Base, Self::Base, bool) {
+    let coordinates = self.to_affine().coordinates();
+    if coordinates.is_some().unwrap_u8() == 1 {
+      (*coordinates.unwrap().x(), *coordinates.unwrap().y(), false)
+    } else {
+      (Self::Base::zero(), Self::Base::zero(), true)
+    }
+  }
 }
 
 impl PrimeField for pallas::Scalar {
@@ -65,6 +94,14 @@ impl PrimeField for pallas::Scalar {
   fn random(rng: &mut (impl RngCore + CryptoRng)) -> Self {
     <Fq as ff::Field>::random(rng)
   }
+
+  fn get_order() -> Integer {
+    Integer::from_str_radix(
+      "40000000000000000000000000000000224698fc0994a8dd8c46eb2100000001",
+      16,
+    )
+    .unwrap()
+  }
 }
 
 impl ChallengeTrait for pallas::Scalar {
@@ -75,12 +112,127 @@ impl ChallengeTrait for pallas::Scalar {
   }
 }
 
-impl CompressedGroup for <pallas::Point as GroupEncoding>::Repr {
+impl CompressedGroup for PallasCompressedElementWrapper {
   type GroupElement = pallas::Point;
-  fn decompress(&self) -> Option<<Self as CompressedGroup>::GroupElement> {
-    Some(Ep::from_bytes(self).unwrap())
+
+  fn decompress(&self) -> Option<pallas::Point> {
+    Some(Ep::from_bytes(&self.repr).unwrap())
+  }
+  fn as_bytes(&self) -> &[u8] {
+    &self.repr
+  }
+}
+
+//////////////////////////////////////Vesta////////////////////////////////////////////////
+
+/// A wrapper for compressed group elements that come from the vesta curve
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+pub struct VestaCompressedElementWrapper {
+  repr: [u8; 32],
+}
+
+impl VestaCompressedElementWrapper {
+  /// Wraps repr into the wrapper
+  pub fn new(repr: [u8; 32]) -> Self {
+    Self { repr }
+  }
+}
+
+unsafe impl Send for VestaCompressedElementWrapper {}
+unsafe impl Sync for VestaCompressedElementWrapper {}
+
+impl Group for vesta::Point {
+  type Base = vesta::Base;
+  type Scalar = vesta::Scalar;
+  type CompressedGroupElement = VestaCompressedElementWrapper;
+
+  fn vartime_multiscalar_mul<I, J>(scalars: I, points: J) -> Self
+  where
+    I: IntoIterator,
+    I::Item: Borrow<Self::Scalar>,
+    J: IntoIterator,
+    J::Item: Borrow<Self>,
+    Self: Clone,
+  {
+    // Unoptimized.
+    scalars
+      .into_iter()
+      .zip(points)
+      .map(|(scalar, point)| (*point.borrow()).mul(*scalar.borrow()))
+      .fold(Eq::group_zero(), |acc, x| acc + x)
+  }
+
+  fn compress(&self) -> Self::CompressedGroupElement {
+    VestaCompressedElementWrapper::new(self.to_bytes())
+  }
+
+  fn from_uniform_bytes(bytes: &[u8]) -> Option<Self> {
+    if bytes.len() != 64 {
+      None
+    } else {
+      let mut arr = [0; 32];
+      arr.copy_from_slice(&bytes[0..32]);
+
+      let hash = Eq::hash_to_curve("from_uniform_bytes");
+      Some(hash(&arr))
+    }
+  }
+
+  fn to_coordinates(&self) -> (Self::Base, Self::Base, bool) {
+    let coordinates = self.to_affine().coordinates();
+    if coordinates.is_some().unwrap_u8() == 1 {
+      (*coordinates.unwrap().x(), *coordinates.unwrap().y(), false)
+    } else {
+      (Self::Base::zero(), Self::Base::zero(), true)
+    }
+  }
+}
+
+impl PrimeField for vesta::Scalar {
+  fn zero() -> Self {
+    Fp::zero()
+  }
+  fn one() -> Self {
+    Fp::one()
+  }
+  fn from_bytes_mod_order_wide(bytes: &[u8]) -> Option<Self> {
+    if bytes.len() != 64 {
+      None
+    } else {
+      let mut arr = [0; 64];
+      arr.copy_from_slice(&bytes[0..64]);
+      Some(Fp::from_bytes_wide(&arr))
+    }
+  }
+
+  fn random(rng: &mut (impl RngCore + CryptoRng)) -> Self {
+    <Fp as ff::Field>::random(rng)
+  }
+
+  fn get_order() -> Integer {
+    Integer::from_str_radix(
+      "40000000000000000000000000000000224698fc094cf91b992d30ed00000001",
+      16,
+    )
+    .unwrap()
+  }
+}
+
+impl ChallengeTrait for vesta::Scalar {
+  fn challenge(label: &'static [u8], transcript: &mut Transcript) -> Self {
+    let mut buf = [0u8; 64];
+    transcript.challenge_bytes(label, &mut buf);
+    vesta::Scalar::from_bytes_mod_order_wide(&buf).unwrap()
+  }
+}
+
+impl CompressedGroup for VestaCompressedElementWrapper {
+  type GroupElement = vesta::Point;
+
+  fn decompress(&self) -> Option<vesta::Point> {
+    Some(Eq::from_bytes(&self.repr).unwrap())
   }
   fn as_bytes(&self) -> &[u8] {
-    self
+    &self.repr
   }
 }
diff --git a/src/poseidon.rs b/src/poseidon.rs
new file mode 100644
index 0000000..48ff399
--- /dev/null
+++ b/src/poseidon.rs
@@ -0,0 +1,220 @@
+//! Poseidon Constants and Poseidon-based RO used in Nova
+use bellperson::{
+  gadgets::{
+    boolean::{AllocatedBit, Boolean},
+    num::AllocatedNum,
+  },
+  ConstraintSystem, SynthesisError,
+};
+use ff::{PrimeField, PrimeFieldBits};
+use generic_array::typenum::{U25, U27, U31};
+use neptune::{
+  circuit::poseidon_hash,
+  poseidon::{Poseidon, PoseidonConstants},
+};
+
+#[cfg(test)]
+use neptune::Strength;
+
+/// All Poseidon Constants that are used in Nova
+#[derive(Clone)]
+pub struct NovaPoseidonConstants<F>
+where
+  F: PrimeField,
+{
+  pub(crate) constants25: PoseidonConstants<F, U25>,
+  pub(crate) constants27: PoseidonConstants<F, U27>,
+  pub(crate) constants31: PoseidonConstants<F, U31>,
+}
+
+#[cfg(test)]
+impl<F> NovaPoseidonConstants<F>
+where
+  F: PrimeField,
+{
+  /// Generate Poseidon constants for the arities that Nova uses
+  pub fn new() -> Self {
+    let constants25 = PoseidonConstants::<F, U25>::new_with_strength(Strength::Strengthened);
+    let constants27 = PoseidonConstants::<F, U27>::new_with_strength(Strength::Strengthened);
+    let constants31 = PoseidonConstants::<F, U31>::new_with_strength(Strength::Strengthened);
+    Self {
+      constants25,
+      constants27,
+      constants31,
+    }
+  }
+}
+
+/// A Poseidon-based RO to use outside circuits
+pub struct PoseidonRO<Scalar>
+where
+  Scalar: PrimeField + PrimeFieldBits,
+{
+  // Internal State
+  state: Vec<Scalar>,
+  // Constants for Poseidon
+  constants: NovaPoseidonConstants<Scalar>,
+}
+
+impl<Scalar> PoseidonRO<Scalar>
+where
+  Scalar: PrimeField + PrimeFieldBits,
+{
+  #[allow(dead_code)]
+  pub fn new(constants: NovaPoseidonConstants<Scalar>) -> Self {
+    Self {
+      state: Vec::new(),
+      constants,
+    }
+  }
+
+  #[allow(dead_code)]
+  /// Flush the state of the RO
+  pub fn flush_state(&mut self) {
+    self.state = Vec::new();
+  }
+
+  /// Absorb a new number into the state of the oracle
+  #[allow(dead_code)]
+  pub fn absorb(&mut self, e: Scalar) {
+    self.state.push(e);
+  }
+
+  /// Compute a challenge by hashing the current state
+  #[allow(dead_code)]
+  pub fn get_challenge(&mut self) -> Scalar {
+    let hash = match self.state.len() {
+      25 => {
+        Poseidon::<Scalar, U25>::new_with_preimage(&self.state, &self.constants.constants25).hash()
+      }
+      27 => {
+        Poseidon::<Scalar, U27>::new_with_preimage(&self.state, &self.constants.constants27).hash()
+      }
+      31 => {
+        Poseidon::<Scalar, U31>::new_with_preimage(&self.state, &self.constants.constants31).hash()
+      }
+      _ => {
+        panic!("Number of elements in the RO state does not match any of the arities used in Nova")
+      }
+    };
+    // Only keep 128 bits
+    let bits = hash.to_le_bits();
+    let mut res = Scalar::zero();
+    let mut coeff = Scalar::one();
+    for bit in bits[0..128].into_iter() {
+      if *bit {
+        res += coeff;
+      }
+      coeff += coeff;
+    }
+    res
+  }
+}
+
+/// A Poseidon-based RO gadget to use inside the verifier circuit.
+pub struct PoseidonROGadget<Scalar>
+where
+  Scalar: PrimeField + PrimeFieldBits,
+{
+  // Internal state
+  state: Vec<AllocatedNum<Scalar>>,
+  constants: NovaPoseidonConstants<Scalar>,
+}
+
+impl<Scalar> PoseidonROGadget<Scalar>
+where
+  Scalar: PrimeField + PrimeFieldBits,
+{
+  /// Initialize the internal state and set the poseidon constants
+  #[allow(dead_code)]
+  pub fn new(constants: NovaPoseidonConstants<Scalar>) -> Self {
+    Self {
+      state: Vec::new(),
+      constants,
+    }
+  }
+
+  /// Flush the state of the RO
+  pub fn flush_state(&mut self) {
+    self.state = Vec::new();
+  }
+
+  /// Absorb a new number into the state of the oracle
+  #[allow(dead_code)]
+  pub fn absorb(&mut self, e: AllocatedNum<Scalar>) {
+    self.state.push(e);
+  }
+
+  /// Compute a challenge by hashing the current state
+  #[allow(dead_code)]
+  pub fn get_challenge<CS>(&mut self, mut cs: CS) -> Result<Vec<AllocatedBit>, SynthesisError>
+  where
+    CS: ConstraintSystem<Scalar>,
+  {
+    let out = match self.state.len() {
+      25 => poseidon_hash(
+        cs.namespace(|| "Poseidon hash"),
+        self.state.clone(),
+        &self.constants.constants25,
+      )?,
+      27 => poseidon_hash(
+        cs.namespace(|| "Poseidon hash"),
+        self.state.clone(),
+        &self.constants.constants27,
+      )?,
+      31 => poseidon_hash(
+        cs.namespace(|| "Poseidon hash"),
+        self.state.clone(),
+        &self.constants.constants31,
+      )?,
+      _ => {
+        panic!("Number of elements in the RO state does not match any of the arities used in Nova")
+      }
+    };
+    // Only keep 128 bits
+    let bits: Vec<AllocatedBit> = out
+      .to_bits_le_strict(cs.namespace(|| "poseidon hash to boolean"))?
+      .iter()
+      .map(|boolean| match boolean {
+        Boolean::Is(ref x) => x.clone(),
+        _ => panic!("Wrong type of input. We should have never reached there"),
+      })
+      .collect();
+    Ok(bits[..128].into())
+  }
+}
+
+#[cfg(test)]
+mod tests {
+  use super::*;
+  type S = pasta_curves::pallas::Scalar;
+  type G = pasta_curves::pallas::Point;
+  use crate::bellperson::solver::SatisfyingAssignment;
+  use crate::gadgets::utils::le_bits_to_num;
+  use crate::traits::PrimeField;
+  use rand::rngs::OsRng;
+
+  #[test]
+  fn test_poseidon_ro() {
+    // Check that the number computed inside the circuit is equal to the number computed outside the circuit
+    let mut csprng: OsRng = OsRng;
+    let constants = NovaPoseidonConstants::new();
+    let mut ro: PoseidonRO<S> = PoseidonRO::new(constants.clone());
+    let mut ro_gadget: PoseidonROGadget<S> = PoseidonROGadget::new(constants);
+    let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+    for i in 0..31 {
+      let num = S::random(&mut csprng);
+      ro.absorb(num);
+      let num_gadget =
+        AllocatedNum::alloc(cs.namespace(|| format!("data {}", i)), || Ok(num)).unwrap();
+      let _ = num_gadget
+        .inputize(&mut cs.namespace(|| format!("input {}", i)))
+        .unwrap();
+      ro_gadget.absorb(num_gadget);
+    }
+    let num = ro.get_challenge();
+    let num2_bits = ro_gadget.get_challenge(&mut cs).unwrap();
+    let num2 = le_bits_to_num(&mut cs, num2_bits).unwrap();
+    assert_eq!(num, num2.get_value().unwrap());
+  }
+}
diff --git a/src/r1cs.rs b/src/r1cs.rs
index 115b2e6..fed9356 100644
--- a/src/r1cs.rs
+++ b/src/r1cs.rs
@@ -9,8 +9,8 @@ use rayon::prelude::*;
 /// Public parameters for a given R1CS
 #[derive(Debug)]
 pub struct R1CSGens<G: Group> {
-  gens_W: CommitGens<G>,
-  gens_E: CommitGens<G>,
+  pub(crate) gens_W: CommitGens<G>, // TODO: avoid pub(crate)
+  pub(crate) gens_E: CommitGens<G>,
 }
 
 /// A type that holds the shape of the R1CS matrices
@@ -47,10 +47,10 @@ pub struct RelaxedR1CSWitness<G: Group> {
 /// A type that holds a Relaxed R1CS instance
 #[derive(Clone, Debug, PartialEq, Eq)]
 pub struct RelaxedR1CSInstance<G: Group> {
-  comm_W: Commitment<G>,
-  comm_E: Commitment<G>,
-  X: Vec<G::Scalar>,
-  u: G::Scalar,
+  pub(crate) comm_W: Commitment<G>,
+  pub(crate) comm_E: Commitment<G>,
+  pub(crate) X: Vec<G::Scalar>,
+  pub(crate) u: G::Scalar,
   Y_last: Vec<G::Scalar>, // output of the last instance that was folded
   counter: usize,         // the number of folds thus far
 }
diff --git a/src/traits.rs b/src/traits.rs
index 864200e..d98c9d8 100644
--- a/src/traits.rs
+++ b/src/traits.rs
@@ -1,9 +1,11 @@
 //! This module defines various traits required by the users of the library to implement.
+use bellperson::{gadgets::num::AllocatedNum, ConstraintSystem, SynthesisError};
 use core::borrow::Borrow;
 use core::fmt::Debug;
 use core::ops::{Add, AddAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 use merlin::Transcript;
 use rand::{CryptoRng, RngCore};
+use rug::Integer;
 
 /// Represents an element of a prime field
 pub trait PrimeField:
@@ -40,6 +42,9 @@ pub trait PrimeField:
 
   /// returns an uniformly random element from the finite field
   fn random(rng: &mut (impl RngCore + CryptoRng)) -> Self;
+
+  /// Get prime field order as a rug::Integer
+  fn get_order() -> Integer;
 }
 
 /// Represents an element of a group
@@ -54,6 +59,9 @@ pub trait Group:
   + ScalarMul<<Self as Group>::Scalar>
   + ScalarMulOwned<<Self as Group>::Scalar>
 {
+  /// A type representing an element of the base field of the group
+  type Base: PrimeField;
+
   /// A type representing an element of the scalar field of the group
   type Scalar: PrimeField + ChallengeTrait;
 
@@ -75,6 +83,9 @@ pub trait Group:
   /// Attempts to create a group element from a sequence of bytes,
   /// failing with a `None` if the supplied bytes do not encode the group element
   fn from_uniform_bytes(bytes: &[u8]) -> Option<Self>;
+
+  /// Returns the affine coordinates (x, y, infinty) for the point
+  fn to_coordinates(&self) -> (Self::Base, Self::Base, bool);
 }
 
 /// Represents a compressed version of a group element
@@ -119,3 +130,13 @@ impl<T, Rhs, Output> ScalarMul<Rhs, Output> for T where T: Mul<Rhs, Output = Out
 /// A helper trait for references implementing group scalar multiplication.
 pub trait ScalarMulOwned<Rhs, Output = Self>: for<'r> ScalarMul<&'r Rhs, Output> {}
 impl<T, Rhs, Output> ScalarMulOwned<Rhs, Output> for T where T: for<'r> ScalarMul<&'r Rhs, Output> {}
+
+///A helper trait for the inner circuit F
+pub trait InnerCircuit<F: PrimeField + ff::PrimeField> {
+  ///Sythesize the circuit for a computation step and return variable that corresponds to z_{i+1}
+  fn synthesize<CS: ConstraintSystem<F>>(
+    &self,
+    cs: &mut CS,
+    z: AllocatedNum<F>,
+  ) -> Result<AllocatedNum<F>, SynthesisError>;
+}
diff --git a/tests/bit.rs b/tests/bit.rs
new file mode 100644
index 0000000..7c29305
--- /dev/null
+++ b/tests/bit.rs
@@ -0,0 +1,50 @@
+use bellperson::{gadgets::num::AllocatedNum, ConstraintSystem, SynthesisError};
+use ff::PrimeField;
+use nova_snark::bellperson::{
+  r1cs::{NovaShape, NovaWitness},
+  shape_cs::ShapeCS,
+  solver::SatisfyingAssignment,
+};
+
+fn synthesize_alloc_bit<Fr: PrimeField, CS: ConstraintSystem<Fr>>(
+  cs: &mut CS,
+) -> Result<(), SynthesisError> {
+  //get two bits as input and check that they are indeed bits
+  let a = AllocatedNum::alloc(cs.namespace(|| "a"), || Ok(Fr::one()))?;
+  let _ = a.inputize(cs.namespace(|| "a is input"));
+  cs.enforce(
+    || "check a is 0 or 1",
+    |lc| lc + CS::one() - a.get_variable(),
+    |lc| lc + a.get_variable(),
+    |lc| lc,
+  );
+  let b = AllocatedNum::alloc(cs.namespace(|| "b"), || Ok(Fr::one()))?;
+  let _ = b.inputize(cs.namespace(|| "b is input"));
+  cs.enforce(
+    || "check b is 0 or 1",
+    |lc| lc + CS::one() - b.get_variable(),
+    |lc| lc + b.get_variable(),
+    |lc| lc,
+  );
+  Ok(())
+}
+
+#[test]
+fn test_alloc_bit() {
+  type G = pasta_curves::pallas::Point;
+
+  //First create the shape
+  let mut cs: ShapeCS<G> = ShapeCS::new();
+  let _ = synthesize_alloc_bit(&mut cs);
+  let shape = cs.r1cs_shape();
+  let gens = cs.r1cs_gens();
+  println!("Mult mod constraint no: {}", cs.num_constraints());
+
+  //Now get the assignment
+  let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+  let _ = synthesize_alloc_bit(&mut cs);
+  let (inst, witness) = cs.r1cs_instance_and_witness(&shape, &gens).unwrap();
+
+  //Make sure that this is satisfiable
+  assert!(shape.is_sat(&gens, &inst, &witness).is_ok());
+}
diff --git a/tests/nonnative.rs b/tests/nonnative.rs
new file mode 100644
index 0000000..007839b
--- /dev/null
+++ b/tests/nonnative.rs
@@ -0,0 +1,307 @@
+use bellperson::{ConstraintSystem, SynthesisError};
+use bellperson_nonnative::{
+  mp::bignat::BigNat,
+  util::{convert::nat_to_f, num::Num},
+};
+use ff::PrimeField;
+use nova_snark::bellperson::{
+  r1cs::{NovaShape, NovaWitness},
+  shape_cs::ShapeCS,
+  solver::SatisfyingAssignment,
+};
+use rug::Integer;
+
+fn synthesize_is_equal<Fr: PrimeField, CS: ConstraintSystem<Fr>>(
+  cs: &mut CS,
+  a_val: &Integer,
+  limb_width: usize,
+  n_limbs: usize,
+) -> Result<(), SynthesisError> {
+  let a1 = BigNat::alloc_from_nat(
+    cs.namespace(|| "alloc a2"),
+    || Ok(a_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let _ = a1.inputize(cs.namespace(|| "make a input"));
+
+  let a_num = Num::alloc(cs.namespace(|| "alloc a num"), || {
+    Ok(nat_to_f(a_val).unwrap())
+  })?;
+  let a2 = BigNat::from_num(
+    cs.namespace(|| "allocate a1_limbs"),
+    a_num,
+    limb_width,
+    n_limbs,
+  )?;
+
+  let _ = a1.equal_when_carried(cs.namespace(|| "check equal"), &a2)?;
+  Ok(())
+}
+
+#[allow(clippy::too_many_arguments)]
+fn synthesize_mult_mod<Fr: PrimeField, CS: ConstraintSystem<Fr>>(
+  cs: &mut CS,
+  a_val: &Integer,
+  b_val: &Integer,
+  m_val: &Integer,
+  q_val: &Integer,
+  r_val: &Integer,
+  limb_width: usize,
+  n_limbs: usize,
+) -> Result<(), SynthesisError> {
+  let a_num = Num::alloc(cs.namespace(|| "alloc a num"), || {
+    Ok(nat_to_f(a_val).unwrap())
+  })?;
+  let m = BigNat::alloc_from_nat(
+    cs.namespace(|| "m"),
+    || Ok(m_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let _ = m.inputize(cs.namespace(|| "modulus m"))?;
+
+  let a = BigNat::from_num(
+    cs.namespace(|| "allocate a_limbs"),
+    a_num,
+    limb_width,
+    n_limbs,
+  )?;
+  let b = BigNat::alloc_from_nat(
+    cs.namespace(|| "b"),
+    || Ok(b_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let q = BigNat::alloc_from_nat(
+    cs.namespace(|| "q"),
+    || Ok(q_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let r = BigNat::alloc_from_nat(
+    cs.namespace(|| "r"),
+    || Ok(r_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let (qa, ra) = a.mult_mod(cs.namespace(|| "prod"), &b, &m)?;
+  qa.equal(cs.namespace(|| "qcheck"), &q)?;
+  ra.equal(cs.namespace(|| "rcheck"), &r)?;
+  Ok(())
+}
+
+fn synthesize_add<Fr: PrimeField, CS: ConstraintSystem<Fr>>(
+  cs: &mut CS,
+  a_val: &Integer,
+  b_val: &Integer,
+  c_val: &Integer,
+  limb_width: usize,
+  n_limbs: usize,
+) -> Result<(), SynthesisError> {
+  let a = BigNat::alloc_from_nat(
+    cs.namespace(|| "a"),
+    || Ok(a_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let _ = a.inputize(cs.namespace(|| "input a"))?;
+  let b = BigNat::alloc_from_nat(
+    cs.namespace(|| "b"),
+    || Ok(b_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let _ = b.inputize(cs.namespace(|| "input b"))?;
+  let c = BigNat::alloc_from_nat(
+    cs.namespace(|| "c"),
+    || Ok(c_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let ca = a.add::<CS>(&b)?;
+  ca.equal(cs.namespace(|| "ccheck"), &c)?;
+  Ok(())
+}
+
+fn synthesize_add_mod<Fr: PrimeField, CS: ConstraintSystem<Fr>>(
+  cs: &mut CS,
+  a_val: &Integer,
+  b_val: &Integer,
+  c_val: &Integer,
+  m_val: &Integer,
+  limb_width: usize,
+  n_limbs: usize,
+) -> Result<(), SynthesisError> {
+  let a = BigNat::alloc_from_nat(
+    cs.namespace(|| "a"),
+    || Ok(a_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let _ = a.inputize(cs.namespace(|| "input a"))?;
+  let b = BigNat::alloc_from_nat(
+    cs.namespace(|| "b"),
+    || Ok(b_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let _ = b.inputize(cs.namespace(|| "input b"))?;
+  let c = BigNat::alloc_from_nat(
+    cs.namespace(|| "c"),
+    || Ok(c_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let m = BigNat::alloc_from_nat(
+    cs.namespace(|| "m"),
+    || Ok(m_val.clone()),
+    limb_width,
+    n_limbs,
+  )?;
+  let d = a.add::<CS>(&b)?;
+  let ca = d.red_mod(cs.namespace(|| "reduce"), &m)?;
+  ca.equal(cs.namespace(|| "ccheck"), &c)?;
+  Ok(())
+}
+
+#[test]
+fn test_mult_mod() {
+  type G = pasta_curves::pallas::Point;
+
+  //Set the inputs
+  let a_val = Integer::from_str_radix(
+    "11572336752428856981970994795408771577024165681374400871001196932361466228192",
+    10,
+  )
+  .unwrap();
+  let b_val = Integer::from_str_radix(
+    "87673389408848523602668121701204553693362841135953267897017930941776218798802",
+    10,
+  )
+  .unwrap();
+  let m_val = Integer::from_str_radix(
+    "40000000000000000000000000000000224698fc094cf91b992d30ed00000001",
+    16,
+  )
+  .unwrap();
+  let q_val = Integer::from_str_radix(
+    "35048542371029440058224000662033175648615707461806414787901284501179083518342",
+    10,
+  )
+  .unwrap();
+  let r_val = Integer::from_str_radix(
+    "26362617993085418618858432307761590013874563896298265114483698919121453084730",
+    10,
+  )
+  .unwrap();
+
+  //First create the shape
+  let mut cs: ShapeCS<G> = ShapeCS::new();
+  let _ = synthesize_mult_mod(&mut cs, &a_val, &b_val, &m_val, &q_val, &r_val, 32, 8);
+  let shape = cs.r1cs_shape();
+  let gens = cs.r1cs_gens();
+  println!("Mult mod constraint no: {}", cs.num_constraints());
+
+  //Now get the assignment
+  let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+  let _ = synthesize_mult_mod(&mut cs, &a_val, &b_val, &m_val, &q_val, &r_val, 32, 8);
+  let (inst, witness) = cs.r1cs_instance_and_witness(&shape, &gens).unwrap();
+
+  //Make sure that this is satisfiable
+  assert!(shape.is_sat(&gens, &inst, &witness).is_ok());
+}
+
+#[test]
+fn test_add() {
+  type G = pasta_curves::pallas::Point;
+
+  //Set the inputs
+  let a_val = Integer::from_str_radix(
+    "11572336752428856981970994795408771577024165681374400871001196932361466228192",
+    10,
+  )
+  .unwrap();
+  let b_val = Integer::from_str_radix("1", 10).unwrap();
+  let c_val = Integer::from_str_radix(
+    "11572336752428856981970994795408771577024165681374400871001196932361466228193",
+    10,
+  )
+  .unwrap();
+
+  //First create the shape
+  let mut cs: ShapeCS<G> = ShapeCS::new();
+  let _ = synthesize_add(&mut cs, &a_val, &b_val, &c_val, 32, 8);
+  let shape = cs.r1cs_shape();
+  let gens = cs.r1cs_gens();
+  println!("Add mod constraint no: {}", cs.num_constraints());
+
+  //Now get the assignment
+  let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+  let _ = synthesize_add(&mut cs, &a_val, &b_val, &c_val, 32, 8);
+  let (inst, witness) = cs.r1cs_instance_and_witness(&shape, &gens).unwrap();
+
+  //Make sure that this is satisfiable
+  assert!(shape.is_sat(&gens, &inst, &witness).is_ok());
+}
+
+#[test]
+fn test_add_mod() {
+  type G = pasta_curves::pallas::Point;
+
+  //Set the inputs
+  let a_val = Integer::from_str_radix(
+    "11572336752428856981970994795408771577024165681374400871001196932361466228192",
+    10,
+  )
+  .unwrap();
+  let b_val = Integer::from_str_radix("1", 10).unwrap();
+  let c_val = Integer::from_str_radix(
+    "11572336752428856981970994795408771577024165681374400871001196932361466228193",
+    10,
+  )
+  .unwrap();
+  let m_val = Integer::from_str_radix(
+    "40000000000000000000000000000000224698fc094cf91b992d30ed00000001",
+    16,
+  )
+  .unwrap();
+
+  //First create the shape
+  let mut cs: ShapeCS<G> = ShapeCS::new();
+  let _ = synthesize_add_mod(&mut cs, &a_val, &b_val, &c_val, &m_val, 32, 8);
+  let shape = cs.r1cs_shape();
+  let gens = cs.r1cs_gens();
+  println!("Add mod constraint no: {}", cs.num_constraints());
+
+  //Now get the assignment
+  let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+  let _ = synthesize_add_mod(&mut cs, &a_val, &b_val, &c_val, &m_val, 32, 8);
+  let (inst, witness) = cs.r1cs_instance_and_witness(&shape, &gens).unwrap();
+
+  //Make sure that this is satisfiable
+  assert!(shape.is_sat(&gens, &inst, &witness).is_ok());
+}
+
+#[test]
+fn test_equal() {
+  type G = pasta_curves::pallas::Point;
+
+  //Set the inputs
+  let a_val = Integer::from_str_radix("1157233675242885698197099479540877", 10).unwrap();
+
+  //First create the shape
+  let mut cs: ShapeCS<G> = ShapeCS::new();
+  let _ = synthesize_is_equal(&mut cs, &a_val, 32, 8);
+  let shape = cs.r1cs_shape();
+  let gens = cs.r1cs_gens();
+  println!("Equal constraint no: {}", cs.num_constraints());
+
+  //Now get the assignment
+  let mut cs: SatisfyingAssignment<G> = SatisfyingAssignment::new();
+  let _ = synthesize_is_equal(&mut cs, &a_val, 32, 8);
+  let (inst, witness) = cs.r1cs_instance_and_witness(&shape, &gens).unwrap();
+
+  //Make sure that this is satisfiable
+  assert!(shape.is_sat(&gens, &inst, &witness).is_ok());
+}