diff --git a/Cargo.toml b/Cargo.toml
index d551a24..39061a4 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -11,6 +11,8 @@ license-file = "LICENSE"
 keywords = ["zkSNARKs", "cryptography", "proofs"]
 
 [dependencies]
+bellperson =  { git = "https://github.com/filecoin-project/bellperson", branch = "nova-minimal" }
+ff = "0.11.0"
 merlin = "2.0.0"
 rand = "0.8.4"
 digest = "0.8.1"
diff --git a/src/bellperson/mod.rs b/src/bellperson/mod.rs
new file mode 100644
index 0000000..73f248d
--- /dev/null
+++ b/src/bellperson/mod.rs
@@ -0,0 +1,7 @@
+//! Support for generating R1CS from [Bellperson].
+//!
+//! [Bellperson]: https://github.com/filecoin-project/bellperson
+
+pub mod prover;
+pub mod r1cs;
+pub mod shape_cs;
diff --git a/src/bellperson/prover.rs b/src/bellperson/prover.rs
new file mode 100644
index 0000000..3f79835
--- /dev/null
+++ b/src/bellperson/prover.rs
@@ -0,0 +1,216 @@
+//! Support for generating R1CS witness using bellperson.
+
+use crate::traits::Group;
+use ff::PrimeField;
+
+use bellperson::{
+  multiexp::DensityTracker, ConstraintSystem, Index, LinearCombination, SynthesisError, Variable,
+};
+
+/// A `ConstraintSystem` which calculates witness values for a concrete instance of an R1CS circuit.
+pub struct ProvingAssignment<G: Group>
+where
+  G::Scalar: PrimeField,
+{
+  // Density of queries
+  a_aux_density: DensityTracker,
+  b_input_density: DensityTracker,
+  b_aux_density: DensityTracker,
+
+  // Evaluations of A, B, C polynomials
+  a: Vec<G::Scalar>,
+  b: Vec<G::Scalar>,
+  c: Vec<G::Scalar>,
+
+  // Assignments of variables
+  pub(crate) input_assignment: Vec<G::Scalar>,
+  pub(crate) aux_assignment: Vec<G::Scalar>,
+}
+use std::fmt;
+
+impl<G: Group> fmt::Debug for ProvingAssignment<G>
+where
+  G::Scalar: PrimeField,
+{
+  fn fmt(&self, fmt: &mut fmt::Formatter) -> fmt::Result {
+    fmt
+      .debug_struct("ProvingAssignment")
+      .field("a_aux_density", &self.a_aux_density)
+      .field("b_input_density", &self.b_input_density)
+      .field("b_aux_density", &self.b_aux_density)
+      .field(
+        "a",
+        &self
+          .a
+          .iter()
+          .map(|v| format!("Fr({:?})", v))
+          .collect::<Vec<_>>(),
+      )
+      .field(
+        "b",
+        &self
+          .b
+          .iter()
+          .map(|v| format!("Fr({:?})", v))
+          .collect::<Vec<_>>(),
+      )
+      .field(
+        "c",
+        &self
+          .c
+          .iter()
+          .map(|v| format!("Fr({:?})", v))
+          .collect::<Vec<_>>(),
+      )
+      .field("input_assignment", &self.input_assignment)
+      .field("aux_assignment", &self.aux_assignment)
+      .finish()
+  }
+}
+
+impl<G: Group> PartialEq for ProvingAssignment<G>
+where
+  G::Scalar: PrimeField,
+{
+  fn eq(&self, other: &ProvingAssignment<G>) -> bool {
+    self.a_aux_density == other.a_aux_density
+      && self.b_input_density == other.b_input_density
+      && self.b_aux_density == other.b_aux_density
+      && self.a == other.a
+      && self.b == other.b
+      && self.c == other.c
+      && self.input_assignment == other.input_assignment
+      && self.aux_assignment == other.aux_assignment
+  }
+}
+
+impl<G: Group> ConstraintSystem<G::Scalar> for ProvingAssignment<G>
+where
+  G::Scalar: PrimeField,
+{
+  type Root = Self;
+
+  fn new() -> Self {
+    Self {
+      a_aux_density: DensityTracker::new(),
+      b_input_density: DensityTracker::new(),
+      b_aux_density: DensityTracker::new(),
+      a: vec![],
+      b: vec![],
+      c: vec![],
+      input_assignment: vec![],
+      aux_assignment: vec![],
+    }
+  }
+
+  fn alloc<F, A, AR>(&mut self, _: A, f: F) -> Result<Variable, SynthesisError>
+  where
+    F: FnOnce() -> Result<G::Scalar, SynthesisError>,
+    A: FnOnce() -> AR,
+    AR: Into<String>,
+  {
+    self.aux_assignment.push(f()?);
+    self.a_aux_density.add_element();
+    self.b_aux_density.add_element();
+
+    Ok(Variable(Index::Aux(self.aux_assignment.len() - 1)))
+  }
+
+  fn alloc_input<F, A, AR>(&mut self, _: A, f: F) -> Result<Variable, SynthesisError>
+  where
+    F: FnOnce() -> Result<G::Scalar, SynthesisError>,
+    A: FnOnce() -> AR,
+    AR: Into<String>,
+  {
+    self.input_assignment.push(f()?);
+    self.b_input_density.add_element();
+
+    Ok(Variable(Index::Input(self.input_assignment.len() - 1)))
+  }
+
+  fn enforce<A, AR, LA, LB, LC>(&mut self, _: A, a: LA, b: LB, c: LC)
+  where
+    A: FnOnce() -> AR,
+    AR: Into<String>,
+    LA: FnOnce(LinearCombination<G::Scalar>) -> LinearCombination<G::Scalar>,
+    LB: FnOnce(LinearCombination<G::Scalar>) -> LinearCombination<G::Scalar>,
+    LC: FnOnce(LinearCombination<G::Scalar>) -> LinearCombination<G::Scalar>,
+  {
+    let a = a(LinearCombination::zero());
+    let b = b(LinearCombination::zero());
+    let c = c(LinearCombination::zero());
+
+    let input_assignment = &self.input_assignment;
+    let aux_assignment = &self.aux_assignment;
+    let a_aux_density = &mut self.a_aux_density;
+    let b_input_density = &mut self.b_input_density;
+    let b_aux_density = &mut self.b_aux_density;
+
+    let a_res = a.eval(
+      // Inputs have full density in the A query
+      // because there are constraints of the
+      // form x * 0 = 0 for each input.
+      None,
+      Some(a_aux_density),
+      input_assignment,
+      aux_assignment,
+    );
+
+    let b_res = b.eval(
+      Some(b_input_density),
+      Some(b_aux_density),
+      input_assignment,
+      aux_assignment,
+    );
+
+    let c_res = c.eval(
+      // There is no C polynomial query,
+      // though there is an (beta)A + (alpha)B + C
+      // query for all aux variables.
+      // However, that query has full density.
+      None,
+      None,
+      input_assignment,
+      aux_assignment,
+    );
+
+    self.a.push(a_res);
+    self.b.push(b_res);
+    self.c.push(c_res);
+  }
+
+  fn push_namespace<NR, N>(&mut self, _: N)
+  where
+    NR: Into<String>,
+    N: FnOnce() -> NR,
+  {
+    // Do nothing; we don't care about namespaces in this context.
+  }
+
+  fn pop_namespace(&mut self) {
+    // Do nothing; we don't care about namespaces in this context.
+  }
+
+  fn get_root(&mut self) -> &mut Self::Root {
+    self
+  }
+
+  fn is_extensible() -> bool {
+    true
+  }
+
+  fn extend(&mut self, other: Self) {
+    self.a_aux_density.extend(other.a_aux_density, false);
+    self.b_input_density.extend(other.b_input_density, true);
+    self.b_aux_density.extend(other.b_aux_density, false);
+
+    self.a.extend(other.a);
+    self.b.extend(other.b);
+    self.c.extend(other.c);
+
+    self.input_assignment
+            // Skip first input, which must have been a temporarily allocated one variable.
+            .extend(&other.input_assignment[1..]);
+    self.aux_assignment.extend(other.aux_assignment);
+  }
+}
diff --git a/src/bellperson/r1cs.rs b/src/bellperson/r1cs.rs
new file mode 100644
index 0000000..fb0d3ee
--- /dev/null
+++ b/src/bellperson/r1cs.rs
@@ -0,0 +1,136 @@
+//! Support for generating R1CS using bellperson.
+
+#![allow(non_snake_case)]
+
+use super::prover::ProvingAssignment;
+use super::shape_cs::ShapeCS;
+use bellperson::{Index, LinearCombination};
+
+use ff::PrimeField;
+
+use crate::{
+  errors::NovaError,
+  r1cs::{R1CSGens, R1CSInstance, R1CSShape, R1CSWitness},
+  traits::Group,
+};
+
+/// `NovaWitness` provide a method for acquiring an `R1CSInstance` and `R1CSWitness` from implementers.
+pub trait NovaWitness<G: Group> {
+    /// Return an instance and witness, given a shape and gens.
+  fn r1cs_instance_and_witness(
+    &self,
+    shape: &R1CSShape<G>,
+    gens: &R1CSGens<G>,
+  ) -> Result<(R1CSInstance<G>, R1CSWitness<G>), NovaError>;
+}
+
+/// `NovaShape` provides methods for acquiring `R1CSShape` and `R1CSGens` from implementers.
+pub trait NovaShape<G: Group> {
+    /// Return an appropriate `R1CSShape` struct.
+  fn r1cs_shape(&self) -> R1CSShape<G>;
+    /// Return an appropriate `R1CSGens` struct.
+  fn r1cs_gens(&self) -> R1CSGens<G>;
+}
+
+impl<G: Group> NovaWitness<G> for ProvingAssignment<G>
+where
+  G::Scalar: PrimeField,
+{
+  fn r1cs_instance_and_witness(
+    &self,
+    shape: &R1CSShape<G>,
+    gens: &R1CSGens<G>,
+  ) -> Result<(R1CSInstance<G>, R1CSWitness<G>), NovaError> {
+    let W = R1CSWitness::<G>::new(shape, &self.aux_assignment)?;
+    let X = &self.input_assignment;
+
+    let comm_W = W.commit(gens);
+
+    let instance = R1CSInstance::<G>::new(shape, &comm_W, X)?;
+
+    Ok((instance, W))
+  }
+}
+
+impl<G: Group> NovaShape<G> for ShapeCS<G>
+where
+  G::Scalar: PrimeField,
+{
+  fn r1cs_shape(&self) -> R1CSShape<G> {
+    let mut A: Vec<(usize, usize, G::Scalar)> = Vec::new();
+    let mut B: Vec<(usize, usize, G::Scalar)> = Vec::new();
+    let mut C: Vec<(usize, usize, G::Scalar)> = Vec::new();
+
+    let mut num_cons_added = 0;
+    let mut X = (&mut A, &mut B, &mut C, &mut num_cons_added);
+
+    let num_inputs = self.num_inputs();
+    let num_constraints = self.num_constraints();
+    let num_vars = self.num_aux();
+
+    for constraint in self.constraints.iter() {
+      add_constraint(
+        &mut X,
+        num_vars,
+        &constraint.0,
+        &constraint.1,
+        &constraint.2,
+      );
+    }
+
+    assert_eq!(num_cons_added, num_constraints);
+
+    let S: R1CSShape<G> = {
+      // Don't count One as an input for shape's purposes.
+      let res = R1CSShape::new(num_constraints, num_vars, num_inputs - 1, &A, &B, &C);
+      res.unwrap()
+    };
+
+    S
+  }
+
+  fn r1cs_gens(&self) -> R1CSGens<G> {
+    R1CSGens::<G>::new(self.num_constraints(), self.num_aux())
+  }
+}
+
+fn add_constraint<S: PrimeField>(
+  X: &mut (
+    &mut Vec<(usize, usize, S)>,
+    &mut Vec<(usize, usize, S)>,
+    &mut Vec<(usize, usize, S)>,
+    &mut usize,
+  ),
+  num_vars: usize,
+  a_lc: &LinearCombination<S>,
+  b_lc: &LinearCombination<S>,
+  c_lc: &LinearCombination<S>,
+) {
+  let (A, B, C, nn) = X;
+  let n = **nn;
+  let one = S::one();
+
+  let add_constraint_component = |index: Index, coeff, V: &mut Vec<_>| {
+    match index {
+      Index::Input(idx) => {
+        // Inputs come last, with input 0, reprsenting 'one',
+        // at position num_vars within the witness vector.
+        let i = idx + num_vars;
+        V.push((n, i, one * coeff))
+      }
+      Index::Aux(idx) => V.push((n, idx, one * coeff)),
+    }
+  };
+
+  for (index, coeff) in a_lc.iter() {
+    add_constraint_component(index.0, coeff, A);
+  }
+  for (index, coeff) in b_lc.iter() {
+    add_constraint_component(index.0, coeff, B)
+  }
+  for (index, coeff) in c_lc.iter() {
+    add_constraint_component(index.0, coeff, C)
+  }
+
+  **nn += 1;
+}
diff --git a/src/bellperson/shape_cs.rs b/src/bellperson/shape_cs.rs
new file mode 100644
index 0000000..2e07a07
--- /dev/null
+++ b/src/bellperson/shape_cs.rs
@@ -0,0 +1,328 @@
+//! Support for generating R1CS shape using bellperson.
+
+use std::cmp::Ordering;
+use std::collections::{BTreeMap, HashMap};
+
+use crate::traits::{Group, PrimeField as PF};
+use ff::{Field, PrimeField};
+
+use bellperson::{ConstraintSystem, Index, LinearCombination, SynthesisError, Variable};
+
+#[derive(Clone, Copy)]
+struct OrderedVariable(Variable);
+
+#[derive(Debug)]
+enum NamedObject {
+  Constraint(usize),
+  Var(Variable),
+  Namespace,
+}
+
+impl Eq for OrderedVariable {}
+impl PartialEq for OrderedVariable {
+  fn eq(&self, other: &OrderedVariable) -> bool {
+    match (self.0.get_unchecked(), other.0.get_unchecked()) {
+      (Index::Input(ref a), Index::Input(ref b)) => a == b,
+      (Index::Aux(ref a), Index::Aux(ref b)) => a == b,
+      _ => false,
+    }
+  }
+}
+impl PartialOrd for OrderedVariable {
+  fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
+    Some(self.cmp(other))
+  }
+}
+impl Ord for OrderedVariable {
+  fn cmp(&self, other: &Self) -> Ordering {
+    match (self.0.get_unchecked(), other.0.get_unchecked()) {
+      (Index::Input(ref a), Index::Input(ref b)) => a.cmp(b),
+      (Index::Aux(ref a), Index::Aux(ref b)) => a.cmp(b),
+      (Index::Input(_), Index::Aux(_)) => Ordering::Less,
+      (Index::Aux(_), Index::Input(_)) => Ordering::Greater,
+    }
+  }
+}
+
+#[allow(clippy::upper_case_acronyms)]
+/// TODO: document
+pub struct ShapeCS<G: Group>
+where
+  G::Scalar: PrimeField + Field,
+{
+  named_objects: HashMap<String, NamedObject>,
+  current_namespace: Vec<String>,
+  #[allow(clippy::type_complexity)]
+  /// TODO: document
+  pub constraints: Vec<(
+    LinearCombination<G::Scalar>,
+    LinearCombination<G::Scalar>,
+    LinearCombination<G::Scalar>,
+    String,
+  )>,
+  inputs: Vec<String>,
+  aux: Vec<String>,
+}
+
+fn proc_lc<Scalar: PrimeField>(
+  terms: &LinearCombination<Scalar>,
+) -> BTreeMap<OrderedVariable, Scalar> {
+  let mut map = BTreeMap::new();
+  for (var, &coeff) in terms.iter() {
+    map
+      .entry(OrderedVariable(var))
+      .or_insert_with(Scalar::zero)
+      .add_assign(&coeff);
+  }
+
+  // Remove terms that have a zero coefficient to normalize
+  let mut to_remove = vec![];
+  for (var, coeff) in map.iter() {
+    if coeff.is_zero().into() {
+      to_remove.push(*var)
+    }
+  }
+
+  for var in to_remove {
+    map.remove(&var);
+  }
+
+  map
+}
+
+impl<G: Group> ShapeCS<G>
+where
+  G::Scalar: PrimeField,
+{
+  /// TODO: document
+  pub fn new() -> Self {
+    ShapeCS::default()
+  }
+
+  // pub fn constraints(
+  //     &self,
+  // ) -> Vec<(
+  //     LinearCombination<G::Scalar>,
+  //     LinearCombination<G::Scalar>,
+  //     LinearCombination<G::Scalar>,
+  //     String,
+  // )> {
+  //     self.constraints.clone()
+  // }
+
+  /// TODO: document
+  pub fn num_constraints(&self) -> usize {
+    self.constraints.len()
+  }
+
+  /// TODO: document
+  pub fn num_inputs(&self) -> usize {
+    self.inputs.len()
+  }
+
+  /// TODO: document
+  pub fn num_aux(&self) -> usize {
+    self.aux.len()
+  }
+
+  /// TODO: document
+  pub fn pretty_print_list(&self) -> Vec<String> {
+    let mut result = Vec::new();
+
+    for input in &self.inputs {
+      result.push(format!("INPUT {}", input));
+    }
+    for aux in &self.aux {
+      result.push(format!("AUX {}", aux));
+    }
+
+    for &(ref _a, ref _b, ref _c, ref name) in &self.constraints {
+      result.push(name.to_string());
+    }
+
+    result
+  }
+
+  /// TODO: document
+  pub fn pretty_print(&self) -> String {
+    let mut s = String::new();
+
+    for input in &self.inputs {
+      s.push_str(&format!("INPUT {}\n", &input))
+    }
+
+    let negone = -<G::Scalar as PF>::one();
+
+    let powers_of_two = (0..G::Scalar::NUM_BITS)
+      .map(|i| G::Scalar::from(2u64).pow_vartime(&[u64::from(i)]))
+      .collect::<Vec<_>>();
+
+    let pp = |s: &mut String, lc: &LinearCombination<G::Scalar>| {
+      s.push('(');
+      let mut is_first = true;
+      for (var, coeff) in proc_lc::<G::Scalar>(&lc) {
+        if coeff == negone {
+          s.push_str(" - ")
+        } else if !is_first {
+          s.push_str(" + ")
+        }
+        is_first = false;
+
+        if coeff != <G::Scalar as PF>::one() && coeff != negone {
+          for (i, x) in powers_of_two.iter().enumerate() {
+            if x == &coeff {
+              s.push_str(&format!("2^{} . ", i));
+              break;
+            }
+          }
+
+          s.push_str(&format!("{:?} . ", coeff))
+        }
+
+        match var.0.get_unchecked() {
+          Index::Input(i) => {
+            s.push_str(&format!("`I{}`", &self.inputs[i]));
+          }
+          Index::Aux(i) => {
+            s.push_str(&format!("`A{}`", &self.aux[i]));
+          }
+        }
+      }
+      if is_first {
+        // Nothing was visited, print 0.
+        s.push('0');
+      }
+      s.push(')');
+    };
+
+    for &(ref a, ref b, ref c, ref name) in &self.constraints {
+      s.push('\n');
+
+      s.push_str(&format!("{}: ", name));
+      pp(&mut s, a);
+      s.push_str(" * ");
+      pp(&mut s, b);
+      s.push_str(" = ");
+      pp(&mut s, c);
+    }
+
+    s.push('\n');
+
+    s
+  }
+
+  /// TODO: document
+  fn set_named_obj(&mut self, path: String, to: NamedObject) {
+    if self.named_objects.contains_key(&path) {
+      panic!("tried to create object at existing path: {}", path);
+    }
+
+    self.named_objects.insert(path, to);
+  }
+}
+
+impl<G: Group> Default for ShapeCS<G>
+where
+  G::Scalar: PrimeField,
+{
+  fn default() -> Self {
+    let mut map = HashMap::new();
+    map.insert("ONE".into(), NamedObject::Var(ShapeCS::<G>::one()));
+    ShapeCS {
+      named_objects: map,
+      current_namespace: vec![],
+      constraints: vec![],
+      inputs: vec![String::from("ONE")],
+      aux: vec![],
+    }
+  }
+}
+
+impl<G: Group> ConstraintSystem<G::Scalar> for ShapeCS<G>
+where
+  G::Scalar: PrimeField,
+{
+  type Root = Self;
+
+  fn alloc<F, A, AR>(&mut self, annotation: A, _f: F) -> Result<Variable, SynthesisError>
+  where
+    F: FnOnce() -> Result<G::Scalar, SynthesisError>,
+    A: FnOnce() -> AR,
+    AR: Into<String>,
+  {
+    let path = compute_path(&self.current_namespace, &annotation().into());
+    self.aux.push(path);
+
+    Ok(Variable::new_unchecked(Index::Aux(self.aux.len() - 1)))
+  }
+
+  fn alloc_input<F, A, AR>(&mut self, annotation: A, _f: F) -> Result<Variable, SynthesisError>
+  where
+    F: FnOnce() -> Result<G::Scalar, SynthesisError>,
+    A: FnOnce() -> AR,
+    AR: Into<String>,
+  {
+    let path = compute_path(&self.current_namespace, &annotation().into());
+    self.inputs.push(path);
+
+    Ok(Variable::new_unchecked(Index::Input(self.inputs.len() - 1)))
+  }
+
+  fn enforce<A, AR, LA, LB, LC>(&mut self, annotation: A, a: LA, b: LB, c: LC)
+  where
+    A: FnOnce() -> AR,
+    AR: Into<String>,
+    LA: FnOnce(LinearCombination<G::Scalar>) -> LinearCombination<G::Scalar>,
+    LB: FnOnce(LinearCombination<G::Scalar>) -> LinearCombination<G::Scalar>,
+    LC: FnOnce(LinearCombination<G::Scalar>) -> LinearCombination<G::Scalar>,
+  {
+    let path = compute_path(&self.current_namespace, &annotation().into());
+    let index = self.constraints.len();
+    self.set_named_obj(path.clone(), NamedObject::Constraint(index));
+
+    let a = a(LinearCombination::zero());
+    let b = b(LinearCombination::zero());
+    let c = c(LinearCombination::zero());
+
+    self.constraints.push((a, b, c, path));
+  }
+
+  fn push_namespace<NR, N>(&mut self, name_fn: N)
+  where
+    NR: Into<String>,
+    N: FnOnce() -> NR,
+  {
+    let name = name_fn().into();
+    let path = compute_path(&self.current_namespace, &name);
+    self.set_named_obj(path, NamedObject::Namespace);
+    self.current_namespace.push(name);
+  }
+
+  fn pop_namespace(&mut self) {
+    assert!(self.current_namespace.pop().is_some());
+  }
+
+  fn get_root(&mut self) -> &mut Self::Root {
+    self
+  }
+}
+
+fn compute_path(ns: &[String], this: &str) -> String {
+  if this.chars().any(|a| a == '/') {
+    panic!("'/' is not allowed in names");
+  }
+
+  let mut name = String::new();
+
+  let mut needs_separation = false;
+  for ns in ns.iter().chain(Some(this.to_string()).iter()) {
+    if needs_separation {
+      name += "/";
+    }
+
+    name += ns;
+    needs_separation = true;
+  }
+
+  name
+}
diff --git a/src/lib.rs b/src/lib.rs
index 2fcd318..0bf82a3 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -5,6 +5,8 @@
 #![deny(missing_docs)]
 
 mod commitments;
+
+pub mod bellperson;
 pub mod errors;
 pub mod pasta;
 pub mod r1cs;