@@ -36,8 +36,7 @@ pub c: Vec<F>,
 }
 
 impl<F: PrimeField> Arith<F> for CCS<F> {
-    /// check that a CCS structure is satisfied by a z vector. Only for testing.
-    fn check_relation(&self, z: &[F]) -> Result<(), Error> {
+    fn eval_relation(&self, z: &[F]) -> Result<Vec<F>, Error> {
         let mut result = vec![F::zero(); self.m];
 
         for i in 0..self.q {
@@ -57,14 +56,7 @@ // add it to the final vector
             result = vec_add(&result, &c_M_j_z)?;
         }
 
-        // make sure the final vector is all zeroes
-        for e in result {
-            if !e.is_zero() {
-                return Err(Error::NotSatisfied);
-            }
-        }
-
-        Ok(())
+        Ok(result)
     }
 
     fn params_to_le_bytes(&self) -> Vec<u8> {
@@ -113,7 +105,10 @@
 #[cfg(test)]
 pub mod tests {
     use super::*;
-    use crate::arith::r1cs::tests::{get_test_r1cs, get_test_z as r1cs_get_test_z};
+    use crate::{
+        arith::r1cs::tests::{get_test_r1cs, get_test_z as r1cs_get_test_z},
+        utils::vec::is_zero_vec,
+    };
     use ark_pallas::Fr;
 
     pub fn get_test_ccs<F: PrimeField>() -> CCS<F> {
@@ -124,9 +119,22 @@ pub fn get_test_z<F: PrimeField>(input: usize) -> Vec<F> {
         r1cs_get_test_z(input)
     }
 
+    #[test]
+    fn test_eval_ccs_relation() {
+        let ccs = get_test_ccs::<Fr>();
+        let mut z = get_test_z(3);
+
+        let f_w = ccs.eval_relation(&z).unwrap();
+        assert!(is_zero_vec(&f_w));
+
+        z[1] = Fr::from(111);
+        let f_w = ccs.eval_relation(&z).unwrap();
+        assert!(!is_zero_vec(&f_w));
+    }
+
     /// Test that a basic CCS relation can be satisfied
     #[test]
-    fn test_ccs_relation() {
+    fn test_check_ccs_relation() {
         let ccs = get_test_ccs::<Fr>();
         let z = get_test_z(3);
 

@@ -6,8 +6,21 @@ pub mod ccs;
 pub mod r1cs;
 
 pub trait Arith<F: PrimeField> {
-    /// Checks that the given Arith structure is satisfied by a z vector. Used only for testing.
-    fn check_relation(&self, z: &[F]) -> Result<(), Error>;
+    /// Evaluate the given Arith structure at `z`, a vector of assignments, and
+    /// return the evaluation.
+    fn eval_relation(&self, z: &[F]) -> Result<Vec<F>, Error>;
+
+    /// Checks that the given Arith structure is satisfied by a z vector, i.e.,
+    /// if the evaluation is a zero vector
+    ///
+    /// Used only for testing.
+    fn check_relation(&self, z: &[F]) -> Result<(), Error> {
+        if self.eval_relation(z)?.iter().all(|f| f.is_zero()) {
+            Ok(())
+        } else {
+            Err(Error::NotSatisfied)
+        }
+    }
 
     /// Returns the bytes that represent the parameters, that is, the matrices sizes, the amount of
     /// public inputs, etc, without the matrices/polynomials values.

@@ -1,15 +1,13 @@
 use crate::commitment::CommitmentScheme;
-use crate::folding::nova::{CommittedInstance, Witness};
 use crate::RngCore;
-use ark_crypto_primitives::sponge::Absorb;
-use ark_ec::{CurveGroup, Group};
+use ark_ec::CurveGroup;
 use ark_ff::PrimeField;
 use ark_relations::r1cs::ConstraintSystem;
 use ark_serialize::{CanonicalDeserialize, CanonicalSerialize};
 use ark_std::rand::Rng;
 
 use super::Arith;
-use crate::utils::vec::{hadamard, mat_vec_mul, vec_add, vec_scalar_mul, vec_sub, SparseMatrix};
+use crate::utils::vec::{hadamard, mat_vec_mul, vec_scalar_mul, vec_sub, SparseMatrix};
 use crate::Error;
 
 #[derive(Debug, Clone, Eq, PartialEq, CanonicalSerialize, CanonicalDeserialize)]
@@ -21,16 +19,24 @@ pub C: SparseMatrix<F>,
 }
 
 impl<F: PrimeField> Arith<F> for R1CS<F> {
-    /// check that a R1CS structure is satisfied by a z vector. Only for testing.
-    fn check_relation(&self, z: &[F]) -> Result<(), Error> {
+    fn eval_relation(&self, z: &[F]) -> Result<Vec<F>, Error> {
+        if z.len() != self.A.n_cols {
+            return Err(Error::NotSameLength(
+                "z.len()".to_string(),
+                z.len(),
+                "number of variables in R1CS".to_string(),
+                self.A.n_cols,
+            ));
+        }
+
         let Az = mat_vec_mul(&self.A, z)?;
         let Bz = mat_vec_mul(&self.B, z)?;
         let Cz = mat_vec_mul(&self.C, z)?;
+        // Multiply Cz by z[0] (u) here, allowing this method to be reused for
+        // both relaxed and unrelaxed R1CS.
+        let uCz = vec_scalar_mul(&Cz, &z[0]);
         let AzBz = hadamard(&Az, &Bz)?;
-        if AzBz != Cz {
-            return Err(Error::NotSatisfied);
-        }
-        Ok(())
+        vec_sub(&AzBz, &uCz)
     }
 
     fn params_to_le_bytes(&self) -> Vec<u8> {
@@ -65,55 +71,50 @@ /// returns a tuple containing (w, x) (witness and public inputs respectively)
     pub fn split_z(&self, z: &[F]) -> (Vec<F>, Vec<F>) {
         (z[self.l + 1..].to_vec(), z[1..self.l + 1].to_vec())
     }
+}
 
-    /// converts the R1CS instance into a RelaxedR1CS as described in
-    /// [Nova](https://eprint.iacr.org/2021/370.pdf) section 4.1.
-    pub fn relax(self) -> RelaxedR1CS<F> {
-        RelaxedR1CS::<F> {
-            l: self.l,
-            E: vec![F::zero(); self.A.n_rows],
-            A: self.A,
-            B: self.B,
-            C: self.C,
-            u: F::one(),
-        }
-    }
-}
+pub trait RelaxedR1CS<C: CurveGroup, W, U>: Arith<C::ScalarField> {
+    /// returns a dummy running instance (Witness and CommittedInstance) for the current R1CS structure
+    fn dummy_running_instance(&self) -> (W, U);
+
+    /// returns a dummy incoming instance (Witness and CommittedInstance) for the current R1CS structure
+    fn dummy_incoming_instance(&self) -> (W, U);
+
+    /// checks if the given instance is relaxed
+    fn is_relaxed(w: &W, u: &U) -> bool;
+
+    /// extracts `z`, the vector of variables, from the given Witness and CommittedInstance
+    fn extract_z(w: &W, u: &U) -> Vec<C::ScalarField>;
 
-#[derive(Debug, Clone, Eq, PartialEq)]
-pub struct RelaxedR1CS<F: PrimeField> {
-    pub l: usize, // io len
-    pub A: SparseMatrix<F>,
-    pub B: SparseMatrix<F>,
-    pub C: SparseMatrix<F>,
-    pub u: F,
-    pub E: Vec<F>,
-}
+    /// checks if the computed error terms correspond to the actual one in `w`
+    /// or `u`
+    fn check_error_terms(w: &W, u: &U, e: Vec<C::ScalarField>) -> Result<(), Error>;
 
-impl<F: PrimeField> RelaxedR1CS<F> {
-    /// check that a RelaxedR1CS structure is satisfied by a z vector. Only for testing.
-    pub fn check_relation(&self, z: &[F]) -> Result<(), Error> {
-        let Az = mat_vec_mul(&self.A, z)?;
-        let Bz = mat_vec_mul(&self.B, z)?;
-        let Cz = mat_vec_mul(&self.C, z)?;
-        let uCz = vec_scalar_mul(&Cz, &self.u);
-        let uCzE = vec_add(&uCz, &self.E)?;
-        let AzBz = hadamard(&Az, &Bz)?;
-        if AzBz != uCzE {
-            return Err(Error::NotSatisfied);
+    /// checks the tight (unrelaxed) R1CS relation
+    fn check_tight_relation(&self, w: &W, u: &U) -> Result<(), Error> {
+        if Self::is_relaxed(w, u) {
+            return Err(Error::R1CSUnrelaxedFail);
         }
 
-        Ok(())
+        let z = Self::extract_z(w, u);
+        self.check_relation(&z)
+    }
+
+    /// checks the relaxed R1CS relation
+    fn check_relaxed_relation(&self, w: &W, u: &U) -> Result<(), Error> {
+        let z = Self::extract_z(w, u);
+        let e = self.eval_relation(&z)?;
+        Self::check_error_terms(w, u, e)
     }
 
     // Computes the E term, given A, B, C, z, u
     fn compute_E(
-        A: &SparseMatrix<F>,
-        B: &SparseMatrix<F>,
-        C: &SparseMatrix<F>,
-        z: &[F],
-        u: &F,
-    ) -> Result<Vec<F>, Error> {
+        A: &SparseMatrix<C::ScalarField>,
+        B: &SparseMatrix<C::ScalarField>,
+        C: &SparseMatrix<C::ScalarField>,
+        z: &[C::ScalarField],
+        u: &C::ScalarField,
+    ) -> Result<Vec<C::ScalarField>, Error> {
         let Az = mat_vec_mul(A, z)?;
         let Bz = mat_vec_mul(B, z)?;
         let AzBz = hadamard(&Az, &Bz)?;
@@ -123,66 +124,9 @@ let uCz = vec_scalar_mul(&Cz, u);
         vec_sub(&AzBz, &uCz)
     }
 
-    pub fn check_sampled_relaxed_r1cs(&self, u: F, E: &[F], z: &[F]) -> bool {
-        let sampled = RelaxedR1CS {
-            l: self.l,
-            A: self.A.clone(),
-            B: self.B.clone(),
-            C: self.C.clone(),
-            u,
-            E: E.to_vec(),
-        };
-        sampled.check_relation(z).is_ok()
-    }
-
-    // Implements sampling a (committed) RelaxedR1CS
-    // See construction 5 in https://eprint.iacr.org/2023/573.pdf
-    pub fn sample<C, CS>(
-        &self,
-        params: &CS::ProverParams,
-        mut rng: impl RngCore,
-    ) -> Result<(CommittedInstance<C>, Witness<C>), Error>
+    fn sample<CS>(&self, params: &CS::ProverParams, rng: impl RngCore) -> Result<(W, U), Error>
     where
-        C: CurveGroup,
-        C: CurveGroup<ScalarField = F>,
-        <C as Group>::ScalarField: Absorb,
-        CS: CommitmentScheme<C, true>,
-    {
-        let u = C::ScalarField::rand(&mut rng);
-        let rE = C::ScalarField::rand(&mut rng);
-        let rW = C::ScalarField::rand(&mut rng);
-
-        let W = (0..self.A.n_cols - self.l - 1)
-            .map(|_| F::rand(&mut rng))
-            .collect();
-        let x = (0..self.l).map(|_| F::rand(&mut rng)).collect::<Vec<F>>();
-        let mut z = vec![u];
-        z.extend(&x);
-        z.extend(&W);
-
-        let E = RelaxedR1CS::compute_E(&self.A, &self.B, &self.C, &z, &u)?;
-
-        debug_assert!(
-            z.len() == self.A.n_cols,
-            "Length of z is {}, while A has {} columns.",
-            z.len(),
-            self.A.n_cols
-        );
-
-        debug_assert!(
-            self.check_sampled_relaxed_r1cs(u, &E, &z),
-            "Sampled a non satisfiable relaxed R1CS, sampled u: {}, computed E: {:?}",
-            u,
-            E
-        );
-
-        let witness = Witness { E, rE, W, rW };
-        let mut cm_witness = witness.commit::<CS, true>(params, x)?;
-
-        // witness.commit() sets u to 1, we set it to the sampled u value
-        cm_witness.u = u;
-        Ok((cm_witness, witness))
-    }
+        CS: CommitmentScheme<C, true>;
 }
 
 /// extracts arkworks ConstraintSystem matrices into crate::utils::vec::SparseMatrix format as R1CS
@@ -229,9 +173,13 @@
 #[cfg(test)]
 pub mod tests {
     use super::*;
+    use crate::folding::nova::{CommittedInstance, Witness};
     use crate::{
         commitment::pedersen::Pedersen,
-        utils::vec::tests::{to_F_matrix, to_F_vec},
+        utils::vec::{
+            is_zero_vec,
+            tests::{to_F_matrix, to_F_vec},
+        },
     };
 
     use ark_pallas::{Fr, Projective};
@@ -242,9 +190,8 @@ let rng = rand::rngs::OsRng;
         let r1cs = get_test_r1cs::<Fr>();
         let (prover_params, _) = Pedersen::<Projective>::setup(rng, r1cs.A.n_rows).unwrap();
 
-        let relaxed_r1cs = r1cs.relax();
-        let sampled =
-            relaxed_r1cs.sample::<Projective, Pedersen<Projective, true>>(&prover_params, rng);
+        let sampled: Result<(Witness<Projective>, CommittedInstance<Projective>), _> =
+            r1cs.sample::<Pedersen<Projective, true>>(&prover_params, rng);
         assert!(sampled.is_ok());
     }
 
@@ -302,10 +249,23 @@ )
     }
 
     #[test]
-    fn test_check_relation() {
+    fn test_eval_r1cs_relation() {
+        let mut rng = ark_std::test_rng();
+        let r1cs = get_test_r1cs::<Fr>();
+        let mut z = get_test_z::<Fr>(rng.gen::<u16>() as usize);
+
+        let f_w = r1cs.eval_relation(&z).unwrap();
+        assert!(is_zero_vec(&f_w));
+
+        z[1] = Fr::from(111);
+        let f_w = r1cs.eval_relation(&z).unwrap();
+        assert!(!is_zero_vec(&f_w));
+    }
+
+    #[test]
+    fn test_check_r1cs_relation() {
         let r1cs = get_test_r1cs::<Fr>();
         let z = get_test_z(5);
         r1cs.check_relation(&z).unwrap();
-        r1cs.relax().check_relation(&z).unwrap();
     }
 }

@@ -360,6 +360,8 @@ transcript.squeeze_bits(NOVA_N_BITS_RO)
     }
 }
 
+/// `CycleFoldConfig` allows us to customize the behavior of CycleFold circuit
+/// according to the folding scheme we are working with.
 pub trait CycleFoldConfig {
     /// `N_INPUT_POINTS` specifies the number of input points that are folded in
     /// [`CycleFoldCircuit`] via random linear combinations.
@@ -465,7 +467,10 @@ // In Nova, this is: x == [r, p1, p2, p3] (wheere p3 is the p_folded).
         // In multifolding schemes such as HyperNova, this is:
         // computed_x = [r, p_0, p_1, p_2, ..., p_n, p_folded],
         // where each p_i is in fact p_i.to_constraint_field()
-        let r_fp = Boolean::le_bits_to_fp_var(&r_bits)?;
+        let r_fp = r_bits
+            .chunks(CFG::F::MODULUS_BIT_SIZE as usize - 1)
+            .map(Boolean::le_bits_to_fp_var)
+            .collect::<Result<Vec<_>, _>>()?;
         let points_aux: Vec<FpVar<CFG::F>> = points
             .iter()
             .map(|p_i| Ok(p_i.to_constraint_field()?[..2].to_vec()))
@@ -475,7 +480,7 @@ .flatten()
             .collect();
 
         let computed_x: Vec<FpVar<CFG::F>> = [
-            vec![r_fp],
+            r_fp,
             points_aux,
             p_folded.to_constraint_field()?[..2].to_vec(),
         ]

@@ -1,10 +1,12 @@
-use ark_ec::{AffineRepr, CurveGroup};
+use ark_ec::{short_weierstrass::SWFlags, AffineRepr, CurveGroup};
+use ark_ff::{Field, PrimeField};
 use ark_r1cs_std::{
     alloc::{AllocVar, AllocationMode},
     fields::fp::FpVar,
-    ToConstraintFieldGadget,
+    R1CSVar, ToConstraintFieldGadget,
 };
-use ark_relations::r1cs::{Namespace, SynthesisError};
+use ark_relations::r1cs::{ConstraintSystemRef, Namespace, SynthesisError};
+use ark_serialize::{CanonicalSerialize, CanonicalSerializeWithFlags};
 use ark_std::Zero;
 use core::borrow::Borrow;
 
@@ -45,6 +47,39 @@ })
     }
 }
 
+impl<C: CurveGroup> R1CSVar<C::ScalarField> for NonNativeAffineVar<C> {
+    type Value = C;
+
+    fn cs(&self) -> ConstraintSystemRef<C::ScalarField> {
+        self.x.cs().or(self.y.cs())
+    }
+
+    fn value(&self) -> Result<Self::Value, SynthesisError> {
+        debug_assert_eq!(C::BaseField::extension_degree(), 1);
+
+        let x = <C::BaseField as Field>::BasePrimeField::from_le_bytes_mod_order(
+            &self.x.value()?.to_bytes_le(),
+        );
+        let y = <C::BaseField as Field>::BasePrimeField::from_le_bytes_mod_order(
+            &self.y.value()?.to_bytes_le(),
+        );
+        let mut bytes = vec![];
+        x.serialize_uncompressed(&mut bytes).unwrap();
+        y.serialize_with_flags(
+            &mut bytes,
+            if x.is_zero() && y.is_zero() {
+                SWFlags::PointAtInfinity
+            } else if y <= -y {
+                SWFlags::YIsPositive
+            } else {
+                SWFlags::YIsNegative
+            },
+        )
+        .unwrap();
+        Ok(C::deserialize_uncompressed_unchecked(&bytes[..]).unwrap())
+    }
+}
+
 impl<C: CurveGroup> ToConstraintFieldGadget<C::ScalarField> for NonNativeAffineVar<C> {
     // Used for converting `NonNativeAffineVar` to a vector of `FpVar` with minimum length in
     // the circuit.
@@ -83,6 +118,10 @@ let x = NonNativeUintVar::inputize(*x);
         let y = NonNativeUintVar::inputize(*y);
         Ok((x, y))
     }
+
+    pub fn zero() -> Self {
+        Self::new_constant(ConstraintSystemRef::None, C::zero()).unwrap()
+    }
 }
 
 impl<C: CurveGroup> AbsorbNonNative<C::ScalarField> for C {
@@ -105,7 +144,6 @@ #[cfg(test)]
 mod tests {
     use super::*;
     use ark_pallas::{Fr, Projective};
-    use ark_r1cs_std::R1CSVar;
     use ark_relations::r1cs::ConstraintSystem;
     use ark_std::UniformRand;
 

@@ -475,30 +475,30 @@ const NU: usize,
 > where
     for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
 {
-    pub _c2: PhantomData<C2>,
-    pub _gc2: PhantomData<GC2>,
-    pub poseidon_config: PoseidonConfig<CF1<C1>>,
-    pub ccs: CCS<C1::ScalarField>, // CCS of the AugmentedFCircuit
-    pub pp_hash: Option<CF1<C1>>,
-    pub i: Option<CF1<C1>>,
-    pub i_usize: Option<usize>,
-    pub z_0: Option<Vec<C1::ScalarField>>,
-    pub z_i: Option<Vec<C1::ScalarField>>,
-    pub external_inputs: Option<Vec<C1::ScalarField>>,
-    pub U_i: Option<LCCCS<C1>>,
-    pub Us: Option<Vec<LCCCS<C1>>>, // other U_i's to be folded that are not the main running instance
-    pub u_i_C: Option<C1>,          // u_i.C
-    pub us: Option<Vec<CCCS<C1>>>, // other u_i's to be folded that are not the main incoming instance
-    pub U_i1_C: Option<C1>,        // U_{i+1}.C
-    pub F: FC,                     // F circuit
-    pub x: Option<CF1<C1>>,        // public input (u_{i+1}.x[0])
-    pub nimfs_proof: Option<NIMFSProof<C1>>,
+    pub(super) _c2: PhantomData<C2>,
+    pub(super) _gc2: PhantomData<GC2>,
+    pub(super) poseidon_config: PoseidonConfig<CF1<C1>>,
+    pub(super) ccs: CCS<C1::ScalarField>, // CCS of the AugmentedFCircuit
+    pub(super) pp_hash: Option<CF1<C1>>,
+    pub(super) i: Option<CF1<C1>>,
+    pub(super) i_usize: Option<usize>,
+    pub(super) z_0: Option<Vec<C1::ScalarField>>,
+    pub(super) z_i: Option<Vec<C1::ScalarField>>,
+    pub(super) external_inputs: Option<Vec<C1::ScalarField>>,
+    pub(super) U_i: Option<LCCCS<C1>>,
+    pub(super) Us: Option<Vec<LCCCS<C1>>>, // other U_i's to be folded that are not the main running instance
+    pub(super) u_i_C: Option<C1>,          // u_i.C
+    pub(super) us: Option<Vec<CCCS<C1>>>, // other u_i's to be folded that are not the main incoming instance
+    pub(super) U_i1_C: Option<C1>,        // U_{i+1}.C
+    pub(super) F: FC,                     // F circuit
+    pub(super) x: Option<CF1<C1>>,        // public input (u_{i+1}.x[0])
+    pub(super) nimfs_proof: Option<NIMFSProof<C1>>,
 
     // cyclefold verifier on C1
-    pub cf_u_i_cmW: Option<C2>, // input, cf_u_i.cmW
-    pub cf_U_i: Option<CycleFoldCommittedInstance<C2>>, // input, RelaxedR1CS CycleFold instance
-    pub cf_x: Option<CF1<C1>>,  // public input (cf_u_{i+1}.x[1])
-    pub cf_cmT: Option<C2>,
+    pub(super) cf_u_i_cmW: Option<C2>, // input, cf_u_i.cmW
+    pub(super) cf_U_i: Option<CycleFoldCommittedInstance<C2>>, // input, RelaxedR1CS CycleFold instance
+    pub(super) cf_x: Option<CF1<C1>>,                          // public input (cf_u_{i+1}.x[1])
+    pub(super) cf_cmT: Option<C2>,
 }
 
 impl<C1, C2, GC2, FC, const MU: usize, const NU: usize> AugmentedFCircuit<C1, C2, GC2, FC, MU, NU>
@@ -891,7 +891,7 @@ use super::*;
     use crate::{
         arith::{
             ccs::tests::{get_test_ccs, get_test_z},
-            r1cs::extract_w_x,
+            r1cs::{extract_w_x, RelaxedR1CS},
         },
         commitment::{pedersen::Pedersen, CommitmentScheme},
         folding::{
@@ -900,9 +900,8 @@ hypernova::{
                 utils::{compute_c, compute_sigmas_thetas},
                 HyperNovaCycleFoldCircuit,
             },
-            nova::traits::NovaR1CS,
         },
-        frontend::tests::CubicFCircuit,
+        frontend::utils::CubicFCircuit,
         transcript::poseidon::poseidon_canonical_config,
         utils::get_cm_coordinates,
     };
@@ -1216,7 +1215,7 @@ let u_dummy = CCCS::<Projective>::dummy(ccs.l);
         let (cf_W_dummy, cf_U_dummy): (
             CycleFoldWitness<Projective2>,
             CycleFoldCommittedInstance<Projective2>,
-        ) = cf_r1cs.dummy_instance();
+        ) = cf_r1cs.dummy_running_instance();
 
         // set the initial dummy instances
         let mut W_i = W_dummy.clone();
@@ -1455,9 +1454,7 @@ // check the new CCCS instance relation
             u_i.check_relation(&ccs, &w_i).unwrap();
 
             // check the CycleFold instance relation
-            cf_r1cs
-                .check_relaxed_instance_relation(&cf_W_i, &cf_U_i)
-                .unwrap();
+            cf_r1cs.check_relaxed_relation(&cf_W_i, &cf_U_i).unwrap();
 
             println!("augmented_f_circuit step {}: {:?}", i, start.elapsed());
         }

@@ -228,7 +228,7 @@
     use super::*;
     use crate::commitment::{kzg::KZG, pedersen::Pedersen};
     use crate::folding::hypernova::PreprocessorParam;
-    use crate::frontend::tests::CubicFCircuit;
+    use crate::frontend::utils::CubicFCircuit;
     use crate::transcript::poseidon::poseidon_canonical_config;
 
     #[test]

@@ -517,7 +517,7 @@
     use super::*;
     use crate::commitment::pedersen::Pedersen;
     use crate::folding::nova::PreprocessorParam;
-    use crate::frontend::tests::CubicFCircuit;
+    use crate::frontend::utils::CubicFCircuit;
     use crate::transcript::poseidon::poseidon_canonical_config;
     use crate::FoldingScheme;
 

@@ -21,7 +21,6 @@ use circuits::AugmentedFCircuit;
 use lcccs::LCCCS;
 use nimfs::NIMFS;
 
-use crate::commitment::CommitmentScheme;
 use crate::constants::NOVA_N_BITS_RO;
 use crate::folding::circuits::{
     cyclefold::{
@@ -30,10 +29,11 @@ CycleFoldWitness,
     },
     CF2,
 };
-use crate::folding::nova::{get_r1cs_from_cs, traits::NovaR1CS, PreprocessorParam};
+use crate::folding::nova::{get_r1cs_from_cs, PreprocessorParam};
 use crate::frontend::FCircuit;
 use crate::utils::{get_cm_coordinates, pp_hash};
 use crate::Error;
+use crate::{arith::r1cs::RelaxedR1CS, commitment::CommitmentScheme};
 use crate::{
     arith::{
         ccs::CCS,
@@ -42,7 +42,8 @@ },
     FoldingScheme, MultiFolding,
 };
 
-struct HyperNovaCycleFoldConfig<C: CurveGroup, const MU: usize, const NU: usize> {
+/// Configuration for HyperNova's CycleFold circuit
+pub struct HyperNovaCycleFoldConfig<C: CurveGroup, const MU: usize, const NU: usize> {
     _c: PhantomData<C>,
 }
 
@@ -55,7 +56,9 @@ type C = C;
     type F = C::BaseField;
 }
 
-type HyperNovaCycleFoldCircuit<C, GC, const MU: usize, const NU: usize> =
+/// CycleFold circuit for computing random linear combinations of group elements
+/// in HyperNova instances.
+pub type HyperNovaCycleFoldCircuit<C, GC, const MU: usize, const NU: usize> =
     CycleFoldCircuit<HyperNovaCycleFoldConfig<C, MU, NU>, GC>;
 
 /// Witness for the LCCCS & CCCS, containing the w vector, and the r_w used as randomness in the Pedersen commitment.
@@ -76,6 +79,7 @@ Witness::<F>::new(vec![F::zero(); ccs.n - ccs.l - 1])
     }
 }
 
+/// Proving parameters for HyperNova-based IVC
 #[derive(Debug, Clone)]
 pub struct ProverParams<C1, C2, CS1, CS2, const H: bool>
 where
@@ -84,13 +88,18 @@ C2: CurveGroup,
     CS1: CommitmentScheme<C1, H>,
     CS2: CommitmentScheme<C2, H>,
 {
+    /// Poseidon sponge configuration
     pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// Proving parameters of the underlying commitment scheme over C1
     pub cs_pp: CS1::ProverParams,
+    /// Proving parameters of the underlying commitment scheme over C2
     pub cf_cs_pp: CS2::ProverParams,
-    // if ccs is set, it will be used, if not, it will be computed at runtime
+    /// CCS of the Augmented Function circuit
+    /// If ccs is set, it will be used, if not, it will be computed at runtime
     pub ccs: Option<CCS<C1::ScalarField>>,
 }
 
+/// Verification parameters for HyperNova-based IVC
 #[derive(Debug, Clone)]
 pub struct VerifierParams<
     C1: CurveGroup,
@@ -99,10 +108,15 @@ CS1: CommitmentScheme<C1, H>,
     CS2: CommitmentScheme<C2, H>,
     const H: bool,
 > {
+    /// Poseidon sponge configuration
     pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// CCS of the Augmented step circuit
     pub ccs: CCS<C1::ScalarField>,
+    /// R1CS of the CycleFold circuit
     pub cf_r1cs: R1CS<C2::ScalarField>,
+    /// Verification parameters of the underlying commitment scheme over C1
     pub cs_vp: CS1::VerifierParams,
+    /// Verification parameters of the underlying commitment scheme over C2
     pub cf_cs_vp: CS2::VerifierParams,
 }
 
@@ -282,7 +296,7 @@ // prepare the initial dummy instances
         let U_i = LCCCS::<C1>::dummy(self.ccs.l, self.ccs.t, self.ccs.s);
         let mut u_i = CCCS::<C1>::dummy(self.ccs.l);
         let (_, cf_U_i): (CycleFoldWitness<C2>, CycleFoldCommittedInstance<C2>) =
-            self.cf_r1cs.dummy_instance();
+            self.cf_r1cs.dummy_running_instance();
 
         let sponge = PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config);
 
@@ -476,7 +490,7 @@ let U_dummy = LCCCS::<C1>::dummy(ccs.l, ccs.t, ccs.s);
         let w_dummy = W_dummy.clone();
         let mut u_dummy = CCCS::<C1>::dummy(ccs.l);
         let (cf_W_dummy, cf_U_dummy): (CycleFoldWitness<C2>, CycleFoldCommittedInstance<C2>) =
-            cf_r1cs.dummy_instance();
+            cf_r1cs.dummy_running_instance();
         u_dummy.x = vec![
             U_dummy.hash(
                 &sponge,
@@ -884,8 +898,7 @@ // check CCCS satisfiability
         u_i.check_relation(&vp.ccs, &w_i)?;
 
         // check CycleFold's RelaxedR1CS satisfiability
-        vp.cf_r1cs
-            .check_relaxed_instance_relation(&cf_W_i, &cf_U_i)?;
+        vp.cf_r1cs.check_relaxed_relation(&cf_W_i, &cf_U_i)?;
 
         Ok(())
     }
@@ -900,7 +913,7 @@ use ark_std::UniformRand;
 
     use super::*;
     use crate::commitment::pedersen::Pedersen;
-    use crate::frontend::tests::CubicFCircuit;
+    use crate::frontend::utils::CubicFCircuit;
     use crate::transcript::poseidon::poseidon_canonical_config;
 
     #[test]

@@ -237,31 +237,31 @@ FC: FCircuit<CF1<C1>>,
 > where
     for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
 {
-    pub _gc2: PhantomData<GC2>,
-    pub poseidon_config: PoseidonConfig<CF1<C1>>,
-    pub pp_hash: Option<CF1<C1>>,
-    pub i: Option<CF1<C1>>,
-    pub i_usize: Option<usize>,
-    pub z_0: Option<Vec<C1::ScalarField>>,
-    pub z_i: Option<Vec<C1::ScalarField>>,
-    pub external_inputs: Option<Vec<C1::ScalarField>>,
-    pub u_i_cmW: Option<C1>,
-    pub U_i: Option<CommittedInstance<C1>>,
-    pub U_i1_cmE: Option<C1>,
-    pub U_i1_cmW: Option<C1>,
-    pub cmT: Option<C1>,
-    pub F: FC,              // F circuit
-    pub x: Option<CF1<C1>>, // public input (u_{i+1}.x[0])
+    pub(super) _gc2: PhantomData<GC2>,
+    pub(super) poseidon_config: PoseidonConfig<CF1<C1>>,
+    pub(super) pp_hash: Option<CF1<C1>>,
+    pub(super) i: Option<CF1<C1>>,
+    pub(super) i_usize: Option<usize>,
+    pub(super) z_0: Option<Vec<C1::ScalarField>>,
+    pub(super) z_i: Option<Vec<C1::ScalarField>>,
+    pub(super) external_inputs: Option<Vec<C1::ScalarField>>,
+    pub(super) u_i_cmW: Option<C1>,
+    pub(super) U_i: Option<CommittedInstance<C1>>,
+    pub(super) U_i1_cmE: Option<C1>,
+    pub(super) U_i1_cmW: Option<C1>,
+    pub(super) cmT: Option<C1>,
+    pub(super) F: FC,              // F circuit
+    pub(super) x: Option<CF1<C1>>, // public input (u_{i+1}.x[0])
 
     // cyclefold verifier on C1
     // Here 'cf1, cf2' are for each of the CycleFold circuits, corresponding to the fold of cmW and
     // cmE respectively
-    pub cf1_u_i_cmW: Option<C2>,                        // input
-    pub cf2_u_i_cmW: Option<C2>,                        // input
-    pub cf_U_i: Option<CycleFoldCommittedInstance<C2>>, // input
-    pub cf1_cmT: Option<C2>,
-    pub cf2_cmT: Option<C2>,
-    pub cf_x: Option<CF1<C1>>, // public input (u_{i+1}.x[1])
+    pub(super) cf1_u_i_cmW: Option<C2>, // input
+    pub(super) cf2_u_i_cmW: Option<C2>, // input
+    pub(super) cf_U_i: Option<CycleFoldCommittedInstance<C2>>, // input
+    pub(super) cf1_cmT: Option<C2>,
+    pub(super) cf2_cmT: Option<C2>,
+    pub(super) cf_x: Option<CF1<C1>>, // public input (u_{i+1}.x[1])
 }
 
 impl<C1: CurveGroup, C2: CurveGroup, GC2: CurveVar<C2, CF2<C2>>, FC: FCircuit<CF1<C1>>>

@@ -342,7 +342,7 @@ use crate::commitment::pedersen::Pedersen;
     use crate::folding::nova::{
         PreprocessorParam, ProverParams as NovaProverParams, VerifierParams as NovaVerifierParams,
     };
-    use crate::frontend::tests::CubicFCircuit;
+    use crate::frontend::utils::CubicFCircuit;
     use crate::transcript::poseidon::poseidon_canonical_config;
 
     #[test]

@@ -577,6 +577,8 @@ }
 
 #[cfg(test)]
 pub mod tests {
+    use std::cmp::max;
+
     use ark_crypto_primitives::crh::{
         sha256::{
             constraints::{Sha256Gadget, UnitVar},
@@ -587,34 +589,61 @@ };
     use ark_pallas::{constraints::GVar, Fq, Fr, Projective};
     use ark_r1cs_std::bits::uint8::UInt8;
     use ark_relations::r1cs::ConstraintSystem;
-    use ark_std::{One, UniformRand};
+    use ark_std::{
+        rand::{thread_rng, Rng},
+        One, UniformRand,
+    };
     use ark_vesta::{constraints::GVar as GVar2, Projective as Projective2};
 
     use super::*;
     use crate::arith::{
         r1cs::{
+            extract_r1cs, extract_w_x,
             tests::{get_test_r1cs, get_test_z},
-            {extract_r1cs, extract_w_x},
+            RelaxedR1CS,
         },
         Arith,
     };
     use crate::commitment::pedersen::Pedersen;
     use crate::folding::nova::PreprocessorParam;
-    use crate::frontend::tests::{CubicFCircuit, CustomFCircuit, WrapperCircuit};
+    use crate::frontend::utils::{CubicFCircuit, CustomFCircuit, WrapperCircuit};
     use crate::transcript::poseidon::poseidon_canonical_config;
     use crate::FoldingScheme;
 
+    fn prepare_instances<C: CurveGroup, CS: CommitmentScheme<C>, R: Rng>(
+        mut rng: R,
+        r1cs: &R1CS<C::ScalarField>,
+        z: &[C::ScalarField],
+    ) -> (Witness<C>, CommittedInstance<C>)
+    where
+        C::ScalarField: Absorb,
+    {
+        let (w, x) = r1cs.split_z(z);
+
+        let (cs_pp, _) = CS::setup(&mut rng, max(w.len(), r1cs.A.n_rows)).unwrap();
+
+        let mut w = Witness::new::<false>(w, r1cs.A.n_rows, &mut rng);
+        w.E = r1cs.eval_relation(z).unwrap();
+        let mut u = w.commit::<CS, false>(&cs_pp, x).unwrap();
+        u.u = z[0];
+
+        (w, u)
+    }
+
     #[test]
     fn test_relaxed_r1cs_small_gadget_handcrafted() {
+        let rng = &mut thread_rng();
+
         let r1cs: R1CS<Fr> = get_test_r1cs();
-        let rel_r1cs = r1cs.clone().relax();
-        let z = get_test_z(3);
+        let mut z = get_test_z(3);
+        z[0] = Fr::rand(rng);
+        let (w, u) = prepare_instances::<_, Pedersen<Projective>, _>(rng, &r1cs, &z);
 
         let cs = ConstraintSystem::<Fr>::new_ref();
 
         let zVar = Vec::<FpVar<Fr>>::new_witness(cs.clone(), || Ok(z)).unwrap();
-        let EVar = Vec::<FpVar<Fr>>::new_witness(cs.clone(), || Ok(rel_r1cs.E)).unwrap();
-        let uVar = FpVar::<Fr>::new_witness(cs.clone(), || Ok(rel_r1cs.u)).unwrap();
+        let EVar = Vec::<FpVar<Fr>>::new_witness(cs.clone(), || Ok(w.E)).unwrap();
+        let uVar = FpVar::<Fr>::new_witness(cs.clone(), || Ok(u.u)).unwrap();
         let r1csVar = R1CSVar::<Fr, Fr, FpVar<Fr>>::new_witness(cs.clone(), || Ok(r1cs)).unwrap();
 
         RelaxedR1CSGadget::check_native(r1csVar, EVar, uVar, zVar).unwrap();
@@ -624,6 +653,8 @@
     // gets as input a circuit that implements the ConstraintSynthesizer trait, and that has been
     // initialized.
     fn test_relaxed_r1cs_gadget<CS: ConstraintSynthesizer<Fr>>(circuit: CS) {
+        let rng = &mut thread_rng();
+
         let cs = ConstraintSystem::<Fr>::new_ref();
 
         circuit.generate_constraints(cs.clone()).unwrap();
@@ -634,18 +665,19 @@ let cs = cs.into_inner().unwrap();
 
         let r1cs = extract_r1cs::<Fr>(&cs);
         let (w, x) = extract_w_x::<Fr>(&cs);
-        let z = [vec![Fr::one()], x, w].concat();
+        let mut z = [vec![Fr::one()], x, w].concat();
         r1cs.check_relation(&z).unwrap();
 
-        let relaxed_r1cs = r1cs.clone().relax();
-        relaxed_r1cs.check_relation(&z).unwrap();
+        z[0] = Fr::rand(rng);
+        let (w, u) = prepare_instances::<_, Pedersen<Projective>, _>(rng, &r1cs, &z);
+        r1cs.check_relaxed_relation(&w, &u).unwrap();
 
         // set new CS for the circuit that checks the RelaxedR1CS of our original circuit
         let cs = ConstraintSystem::<Fr>::new_ref();
         // prepare the inputs for our circuit
         let zVar = Vec::<FpVar<Fr>>::new_witness(cs.clone(), || Ok(z)).unwrap();
-        let EVar = Vec::<FpVar<Fr>>::new_witness(cs.clone(), || Ok(relaxed_r1cs.E)).unwrap();
-        let uVar = FpVar::<Fr>::new_witness(cs.clone(), || Ok(relaxed_r1cs.u)).unwrap();
+        let EVar = Vec::<FpVar<Fr>>::new_witness(cs.clone(), || Ok(w.E)).unwrap();
+        let uVar = FpVar::<Fr>::new_witness(cs.clone(), || Ok(u.u)).unwrap();
         let r1csVar = R1CSVar::<Fr, Fr, FpVar<Fr>>::new_witness(cs.clone(), || Ok(r1cs)).unwrap();
 
         RelaxedR1CSGadget::check_native(r1csVar, EVar, uVar, zVar).unwrap();
@@ -709,6 +741,8 @@ }
 
     #[test]
     fn test_relaxed_r1cs_nonnative_circuit() {
+        let rng = &mut thread_rng();
+
         let cs = ConstraintSystem::<Fq>::new_ref();
         // in practice we would use CycleFoldCircuit, but is a very big circuit (when computed
         // non-natively inside the RelaxedR1CS circuit), so in order to have a short test we use a
@@ -725,16 +759,15 @@ cs.finalize();
         let cs = cs.into_inner().unwrap();
         let r1cs = extract_r1cs::<Fq>(&cs);
         let (w, x) = extract_w_x::<Fq>(&cs);
-        let z = [vec![Fq::one()], x, w].concat();
+        let z = [vec![Fq::rand(rng)], x, w].concat();
 
-        let relaxed_r1cs = r1cs.clone().relax();
+        let (w, u) = prepare_instances::<_, Pedersen<Projective2>, _>(rng, &r1cs, &z);
 
         // natively
         let cs = ConstraintSystem::<Fq>::new_ref();
         let zVar = Vec::<FpVar<Fq>>::new_witness(cs.clone(), || Ok(z.clone())).unwrap();
-        let EVar =
-            Vec::<FpVar<Fq>>::new_witness(cs.clone(), || Ok(relaxed_r1cs.clone().E)).unwrap();
-        let uVar = FpVar::<Fq>::new_witness(cs.clone(), || Ok(relaxed_r1cs.u)).unwrap();
+        let EVar = Vec::<FpVar<Fq>>::new_witness(cs.clone(), || Ok(w.E.clone())).unwrap();
+        let uVar = FpVar::<Fq>::new_witness(cs.clone(), || Ok(u.u)).unwrap();
         let r1csVar =
             R1CSVar::<Fq, Fq, FpVar<Fq>>::new_witness(cs.clone(), || Ok(r1cs.clone())).unwrap();
         RelaxedR1CSGadget::check_native(r1csVar, EVar, uVar, zVar).unwrap();
@@ -742,8 +775,8 @@
         // non-natively
         let cs = ConstraintSystem::<Fr>::new_ref();
         let zVar = Vec::new_witness(cs.clone(), || Ok(z)).unwrap();
-        let EVar = Vec::new_witness(cs.clone(), || Ok(relaxed_r1cs.E)).unwrap();
-        let uVar = NonNativeUintVar::<Fr>::new_witness(cs.clone(), || Ok(relaxed_r1cs.u)).unwrap();
+        let EVar = Vec::new_witness(cs.clone(), || Ok(w.E)).unwrap();
+        let uVar = NonNativeUintVar::<Fr>::new_witness(cs.clone(), || Ok(u.u)).unwrap();
         let r1csVar =
             R1CSVar::<Fq, Fr, NonNativeUintVar<Fr>>::new_witness(cs.clone(), || Ok(r1cs)).unwrap();
         RelaxedR1CSGadget::check_nonnative(r1csVar, EVar, uVar, zVar).unwrap();

@@ -14,7 +14,6 @@ use ark_std::rand::RngCore;
 use ark_std::{One, UniformRand, Zero};
 use core::marker::PhantomData;
 
-use crate::commitment::CommitmentScheme;
 use crate::folding::circuits::cyclefold::{
     fold_cyclefold_circuit, CycleFoldCircuit, CycleFoldCommittedInstance, CycleFoldConfig,
     CycleFoldWitness,
@@ -25,6 +24,7 @@ use crate::transcript::{poseidon::poseidon_canonical_config, AbsorbNonNative, Transcript};
 use crate::utils::vec::is_zero_vec;
 use crate::Error;
 use crate::FoldingScheme;
+use crate::{arith::r1cs::RelaxedR1CS, commitment::CommitmentScheme};
 use crate::{
     arith::r1cs::{extract_r1cs, extract_w_x, R1CS},
     constants::NOVA_N_BITS_RO,
@@ -40,8 +40,8 @@ pub mod traits;
 pub mod zk;
 use circuits::{AugmentedFCircuit, ChallengeGadget};
 use nifs::NIFS;
-use traits::NovaR1CS;
 
+/// Configuration for Nova's CycleFold circuit
 pub struct NovaCycleFoldConfig<C: CurveGroup> {
     _c: PhantomData<C>,
 }
@@ -56,7 +56,9 @@ type C = C;
     type F = C::BaseField;
 }
 
-type NovaCycleFoldCircuit<C, GC> = CycleFoldCircuit<NovaCycleFoldConfig<C>, GC>;
+/// CycleFold circuit for computing random linear combinations of group elements
+/// in Nova instances.
+pub type NovaCycleFoldCircuit<C, GC> = CycleFoldCircuit<NovaCycleFoldConfig<C>, GC>;
 
 #[derive(Debug, Clone, Eq, PartialEq, CanonicalSerialize, CanonicalDeserialize)]
 pub struct CommittedInstance<C: CurveGroup> {
@@ -136,10 +138,7 @@ pub W: Vec<C::ScalarField>,
     pub rW: C::ScalarField,
 }
 
-impl<C: CurveGroup> Witness<C>
-where
-    <C as Group>::ScalarField: Absorb,
-{
+impl<C: CurveGroup> Witness<C> {
     pub fn new<const H: bool>(w: Vec<C::ScalarField>, e_len: usize, mut rng: impl RngCore) -> Self {
         let (rW, rE) = if H {
             (
@@ -227,6 +226,7 @@ }
     }
 }
 
+/// Proving parameters for Nova-based IVC
 #[derive(Debug, Clone)]
 pub struct ProverParams<C1, C2, CS1, CS2, const H: bool = false>
 where
@@ -235,8 +235,11 @@ C2: CurveGroup,
     CS1: CommitmentScheme<C1, H>,
     CS2: CommitmentScheme<C2, H>,
 {
+    /// Poseidon sponge configuration
     pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// Proving parameters of the underlying commitment scheme over C1
     pub cs_pp: CS1::ProverParams,
+    /// Proving parameters of the underlying commitment scheme over C2
     pub cf_cs_pp: CS2::ProverParams,
 }
 
@@ -302,6 +305,7 @@ })
     }
 }
 
+/// Verification parameters for Nova-based IVC
 #[derive(Debug, Clone)]
 pub struct VerifierParams<C1, C2, CS1, CS2, const H: bool = false>
 where
@@ -310,10 +314,15 @@ C2: CurveGroup,
     CS1: CommitmentScheme<C1, H>,
     CS2: CommitmentScheme<C2, H>,
 {
+    /// Poseidon sponge configuration
     pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// R1CS of the Augmented step circuit
     pub r1cs: R1CS<C1::ScalarField>,
+    /// R1CS of the CycleFold circuit
     pub cf_r1cs: R1CS<C2::ScalarField>,
+    /// Verification parameters of the underlying commitment scheme over C1
     pub cs_vp: CS1::VerifierParams,
+    /// Verification parameters of the underlying commitment scheme over C2
     pub cf_cs_vp: CS2::VerifierParams,
 }
 
@@ -553,8 +562,9 @@ // compute the public params hash
         let pp_hash = vp.pp_hash()?;
 
         // setup the dummy instances
-        let (w_dummy, u_dummy) = r1cs.dummy_instance();
-        let (cf_w_dummy, cf_u_dummy) = cf_r1cs.dummy_instance();
+        let (W_dummy, U_dummy) = r1cs.dummy_running_instance();
+        let (w_dummy, u_dummy) = r1cs.dummy_incoming_instance();
+        let (cf_W_dummy, cf_U_dummy) = cf_r1cs.dummy_running_instance();
 
         // W_dummy=W_0 is a 'dummy witness', all zeroes, but with the size corresponding to the
         // R1CS that we're working with.
@@ -572,13 +582,13 @@ pp_hash,
             i: C1::ScalarField::zero(),
             z_0: z_0.clone(),
             z_i: z_0,
-            w_i: w_dummy.clone(),
-            u_i: u_dummy.clone(),
-            W_i: w_dummy,
-            U_i: u_dummy,
+            w_i: w_dummy,
+            u_i: u_dummy,
+            W_i: W_dummy,
+            U_i: U_dummy,
             // cyclefold running instance
-            cf_W_i: cf_w_dummy.clone(),
-            cf_U_i: cf_u_dummy.clone(),
+            cf_W_i: cf_W_dummy,
+            cf_U_i: cf_U_dummy,
         })
     }
 
@@ -805,10 +815,10 @@ self.cf_U_i = cf_U_i1;
 
             #[cfg(test)]
             {
-                self.cf_r1cs.check_instance_relation(&_cfW_w_i, &cfW_u_i)?;
-                self.cf_r1cs.check_instance_relation(&_cfE_w_i, &cfE_u_i)?;
+                self.cf_r1cs.check_tight_relation(&_cfW_w_i, &cfW_u_i)?;
+                self.cf_r1cs.check_tight_relation(&_cfE_w_i, &cfE_u_i)?;
                 self.cf_r1cs
-                    .check_relaxed_instance_relation(&self.cf_W_i, &self.cf_U_i)?;
+                    .check_relaxed_relation(&self.cf_W_i, &self.cf_U_i)?;
             }
         }
 
@@ -840,9 +850,8 @@ self.U_i = U_i1;
 
         #[cfg(test)]
         {
-            self.r1cs.check_instance_relation(&self.w_i, &self.u_i)?;
-            self.r1cs
-                .check_relaxed_instance_relation(&self.W_i, &self.U_i)?;
+            self.r1cs.check_tight_relation(&self.w_i, &self.u_i)?;
+            self.r1cs.check_relaxed_relation(&self.W_i, &self.U_i)?;
         }
 
         Ok(())
@@ -908,19 +917,13 @@ if expected_cf_u_i_x != u_i.x[1] {
             return Err(Error::IVCVerificationFail);
         }
 
-        // check u_i.cmE==0, u_i.u==1 (=u_i is a un-relaxed instance)
-        if !u_i.cmE.is_zero() || !u_i.u.is_one() {
-            return Err(Error::IVCVerificationFail);
-        }
-
-        // check R1CS satisfiability
-        vp.r1cs.check_instance_relation(&w_i, &u_i)?;
+        // check R1CS satisfiability, which also enforces u_i.cmE==0, u_i.u==1
+        vp.r1cs.check_tight_relation(&w_i, &u_i)?;
         // check RelaxedR1CS satisfiability
-        vp.r1cs.check_relaxed_instance_relation(&W_i, &U_i)?;
+        vp.r1cs.check_relaxed_relation(&W_i, &U_i)?;
 
         // check CycleFold RelaxedR1CS satisfiability
-        vp.cf_r1cs
-            .check_relaxed_instance_relation(&cf_W_i, &cf_U_i)?;
+        vp.cf_r1cs.check_relaxed_relation(&cf_W_i, &cf_U_i)?;
 
         Ok(())
     }
@@ -1077,7 +1080,7 @@ use ark_grumpkin::{constraints::GVar as GVar2, Projective as Projective2};
 
     use super::*;
     use crate::commitment::pedersen::Pedersen;
-    use crate::frontend::tests::CubicFCircuit;
+    use crate::frontend::utils::CubicFCircuit;
     use crate::transcript::poseidon::poseidon_canonical_config;
 
     /// This test tests the Nova+CycleFold IVC, and by consequence it is also testing the

@@ -210,10 +210,12 @@ use ark_ff::{BigInteger, PrimeField};
     use ark_pallas::{Fr, Projective};
     use ark_std::{ops::Mul, test_rng, UniformRand};
 
-    use crate::arith::r1cs::tests::{get_test_r1cs, get_test_z};
+    use crate::arith::r1cs::{
+        tests::{get_test_r1cs, get_test_z},
+        RelaxedR1CS,
+    };
     use crate::commitment::pedersen::{Params as PedersenParams, Pedersen};
     use crate::folding::nova::circuits::ChallengeGadget;
-    use crate::folding::nova::traits::NovaR1CS;
     use crate::transcript::poseidon::poseidon_canonical_config;
 
     #[allow(clippy::type_complexity)]
@@ -316,8 +318,8 @@ let w_i = w_dummy.clone();
         let u_i = u_dummy.clone();
         let W_i = w_dummy.clone();
         let U_i = u_dummy.clone();
-        r1cs.check_relaxed_instance_relation(&w_i, &u_i).unwrap();
-        r1cs.check_relaxed_instance_relation(&W_i, &U_i).unwrap();
+        r1cs.check_relaxed_relation(&w_i, &u_i).unwrap();
+        r1cs.check_relaxed_relation(&W_i, &U_i).unwrap();
 
         let r_Fr = Fr::from(3_u32);
 
@@ -334,7 +336,7 @@ let (W_i1, U_i1) = NIFS::<Projective, Pedersen<Projective>>::fold_instances(
             r_Fr, &w_i, &u_i, &W_i, &U_i, &T, cmT,
         )
         .unwrap();
-        r1cs.check_relaxed_instance_relation(&W_i1, &U_i1).unwrap();
+        r1cs.check_relaxed_relation(&W_i1, &U_i1).unwrap();
     }
 
     // fold 2 instances into one
@@ -348,9 +350,9 @@ let ci3_v = NIFS::<Projective, Pedersen<Projective>>::verify(r, &ci1, &ci2, &cmT);
         assert_eq!(ci3_v, ci3);
 
         // check that relations hold for the 2 inputted instances and the folded one
-        r1cs.check_relaxed_instance_relation(&w1, &ci1).unwrap();
-        r1cs.check_relaxed_instance_relation(&w2, &ci2).unwrap();
-        r1cs.check_relaxed_instance_relation(&w3, &ci3).unwrap();
+        r1cs.check_relaxed_relation(&w1, &ci1).unwrap();
+        r1cs.check_relaxed_relation(&w2, &ci2).unwrap();
+        r1cs.check_relaxed_relation(&w3, &ci3).unwrap();
 
         // check that folded commitments from folded instance (ci) are equal to folding the
         // use folded rE, rW to commit w3
@@ -425,7 +427,7 @@ let mut running_committed_instance = running_instance_w
             .commit::<Pedersen<Projective>, false>(&pedersen_params, x)
             .unwrap();
 
-        r1cs.check_relaxed_instance_relation(&running_instance_w, &running_committed_instance)
+        r1cs.check_relaxed_relation(&running_instance_w, &running_committed_instance)
             .unwrap();
 
         let num_iters = 10;
@@ -438,11 +440,8 @@ Witness::<Projective>::new::<false>(w.clone(), r1cs.A.n_rows, test_rng());
             let incoming_committed_instance = incoming_instance_w
                 .commit::<Pedersen<Projective>, false>(&pedersen_params, x)
                 .unwrap();
-            r1cs.check_relaxed_instance_relation(
-                &incoming_instance_w,
-                &incoming_committed_instance,
-            )
-            .unwrap();
+            r1cs.check_relaxed_relation(&incoming_instance_w, &incoming_committed_instance)
+                .unwrap();
 
             let r = Fr::rand(&mut rng); // folding challenge would come from the RO
 
@@ -475,7 +474,7 @@ &incoming_committed_instance,
                 &cmT,
             );
 
-            r1cs.check_relaxed_instance_relation(&folded_w, &folded_committed_instance)
+            r1cs.check_relaxed_relation(&folded_w, &folded_committed_instance)
                 .unwrap();
 
             // set running_instance for next loop iteration

@@ -187,7 +187,7 @@
     use crate::{
         commitment::{kzg::KZG, pedersen::Pedersen},
         folding::nova::{Nova, PreprocessorParam},
-        frontend::{tests::CubicFCircuit, FCircuit},
+        frontend::{utils::CubicFCircuit, FCircuit},
         transcript::poseidon::poseidon_canonical_config,
         FoldingScheme,
     };

@@ -1,69 +1,90 @@
-use ark_crypto_primitives::sponge::Absorb;
-use ark_ec::{CurveGroup, Group};
-use ark_std::One;
+use ark_ec::CurveGroup;
+use ark_std::{rand::RngCore, One, UniformRand};
 
 use super::{CommittedInstance, Witness};
-use crate::arith::{r1cs::R1CS, Arith};
+use crate::arith::r1cs::{RelaxedR1CS, R1CS};
 use crate::Error;
 
-/// NovaR1CS extends R1CS methods with Nova specific methods
-pub trait NovaR1CS<C: CurveGroup> {
-    /// returns a dummy instance (Witness and CommittedInstance) for the current R1CS structure
-    fn dummy_instance(&self) -> (Witness<C>, CommittedInstance<C>);
-
-    /// checks the R1CS relation (un-relaxed) for the given Witness and CommittedInstance.
-    fn check_instance_relation(
-        &self,
-        W: &Witness<C>,
-        U: &CommittedInstance<C>,
-    ) -> Result<(), Error>;
-
-    /// checks the Relaxed R1CS relation (corresponding to the current R1CS) for the given Witness
-    /// and CommittedInstance.
-    fn check_relaxed_instance_relation(
-        &self,
-        W: &Witness<C>,
-        U: &CommittedInstance<C>,
-    ) -> Result<(), Error>;
-}
-
-impl<C: CurveGroup> NovaR1CS<C> for R1CS<C::ScalarField>
-where
-    <C as Group>::ScalarField: Absorb,
-    <C as ark_ec::CurveGroup>::BaseField: ark_ff::PrimeField,
-{
-    fn dummy_instance(&self) -> (Witness<C>, CommittedInstance<C>) {
+impl<C: CurveGroup> RelaxedR1CS<C, Witness<C>, CommittedInstance<C>> for R1CS<C::ScalarField> {
+    fn dummy_running_instance(&self) -> (Witness<C>, CommittedInstance<C>) {
         let w_len = self.A.n_cols - 1 - self.l;
         let w_dummy = Witness::<C>::dummy(w_len, self.A.n_rows);
         let u_dummy = CommittedInstance::<C>::dummy(self.l);
         (w_dummy, u_dummy)
     }
 
-    // notice that this method does not check the commitment correctness
-    fn check_instance_relation(
-        &self,
-        W: &Witness<C>,
-        U: &CommittedInstance<C>,
+    fn dummy_incoming_instance(&self) -> (Witness<C>, CommittedInstance<C>) {
+        self.dummy_running_instance()
+    }
+
+    fn is_relaxed(_w: &Witness<C>, u: &CommittedInstance<C>) -> bool {
+        u.cmE != C::zero() || u.u != C::ScalarField::one()
+    }
+
+    fn extract_z(w: &Witness<C>, u: &CommittedInstance<C>) -> Vec<C::ScalarField> {
+        [&[u.u][..], &u.x, &w.W].concat()
+    }
+
+    fn check_error_terms(
+        w: &Witness<C>,
+        _u: &CommittedInstance<C>,
+        e: Vec<C::ScalarField>,
     ) -> Result<(), Error> {
-        if U.cmE != C::zero() || U.u != C::ScalarField::one() {
-            return Err(Error::R1CSUnrelaxedFail);
+        if w.E == e {
+            Ok(())
+        } else {
+            Err(Error::NotSatisfied)
         }
-
-        let Z: Vec<C::ScalarField> = [vec![U.u], U.x.to_vec(), W.W.to_vec()].concat();
-        self.check_relation(&Z)
     }
 
-    // notice that this method does not check the commitment correctness
-    fn check_relaxed_instance_relation(
+    fn sample<CS>(
         &self,
-        W: &Witness<C>,
-        U: &CommittedInstance<C>,
-    ) -> Result<(), Error> {
-        let mut rel_r1cs = self.clone().relax();
-        rel_r1cs.u = U.u;
-        rel_r1cs.E = W.E.clone();
+        params: &CS::ProverParams,
+        mut rng: impl RngCore,
+    ) -> Result<(Witness<C>, CommittedInstance<C>), Error>
+    where
+        CS: crate::commitment::CommitmentScheme<C, true>,
+    {
+        // Implements sampling a (committed) RelaxedR1CS
+        // See construction 5 in https://eprint.iacr.org/2023/573.pdf
+        let u = C::ScalarField::rand(&mut rng);
+        let rE = C::ScalarField::rand(&mut rng);
+        let rW = C::ScalarField::rand(&mut rng);
+
+        let W = (0..self.A.n_cols - self.l - 1)
+            .map(|_| C::ScalarField::rand(&mut rng))
+            .collect();
+        let x = (0..self.l)
+            .map(|_| C::ScalarField::rand(&mut rng))
+            .collect::<Vec<C::ScalarField>>();
+        let mut z = vec![u];
+        z.extend(&x);
+        z.extend(&W);
+
+        let E = <Self as RelaxedR1CS<C, Witness<C>, CommittedInstance<C>>>::compute_E(
+            &self.A, &self.B, &self.C, &z, &u,
+        )?;
 
-        let Z: Vec<C::ScalarField> = [vec![U.u], U.x.to_vec(), W.W.to_vec()].concat();
-        rel_r1cs.check_relation(&Z)
+        debug_assert!(
+            z.len() == self.A.n_cols,
+            "Length of z is {}, while A has {} columns.",
+            z.len(),
+            self.A.n_cols
+        );
+
+        let witness = Witness { E, rE, W, rW };
+        let mut cm_witness = witness.commit::<CS, true>(params, x)?;
+
+        // witness.commit() sets u to 1, we set it to the sampled u value
+        cm_witness.u = u;
+
+        debug_assert!(
+            self.check_relaxed_relation(&witness, &cm_witness).is_ok(),
+            "Sampled a non satisfiable relaxed R1CS, sampled u: {}, computed E: {:?}",
+            u,
+            witness.E
+        );
+
+        Ok((witness, cm_witness))
     }
 }

@@ -30,7 +30,6 @@ /// instance (steps 2,3,4 from section D.4 of the [HyperNova](https://eprint.iacr.org/2023/573.pdf)
 /// paper).
 /// And the Use-case-2 would require a modified version of the Decider circuits.
 ///
-use crate::folding::nova::traits::NovaR1CS;
 use ark_crypto_primitives::sponge::CryptographicSponge;
 use ark_ff::{BigInteger, PrimeField};
 use ark_std::{One, Zero};
@@ -141,9 +140,8 @@
         // d. Store folding proof
         let pi = FoldingProof { cmT };
 
-        // 2. Sample a satisfying relaxed R1CS instance-witness pair (U_r, W_r)
-        let relaxed_instance = nova.r1cs.clone().relax();
-        let (U_r, W_r) = relaxed_instance.sample::<C1, CS1>(&nova.cs_pp, &mut rng)?;
+        // 2. Sample a satisfying relaxed R1CS instance-witness pair (W_r, U_r)
+        let (W_r, U_r) = nova.r1cs.sample::<CS1>(&nova.cs_pp, &mut rng)?;
 
         // 3. Fold the instance-witness pair (U_f, W_f) with (U_r, W_r)
         // a. Compute T
@@ -280,21 +278,10 @@ &proof.pi_prime.cmT,
         );
 
         // 5. Check that W^{\prime}_i is a satisfying witness
-        let mut z = vec![U_i_prime.u];
-        z.extend(&U_i_prime.x);
-        z.extend(&proof.W_i_prime.W);
-        let relaxed_r1cs = RelaxedR1CS {
-            l: r1cs.l,
-            A: r1cs.A.clone(),
-            B: r1cs.B.clone(),
-            C: r1cs.C.clone(),
-            u: U_i_prime.u,
-            E: proof.W_i_prime.E.clone(),
-        };
-        relaxed_r1cs.check_relation(&z)?;
+        r1cs.check_relaxed_relation(&proof.W_i_prime, &U_i_prime)?;
 
         // 6. Check that the cyclefold instance-witness pair satisfies the cyclefold relaxed r1cs
-        cf_r1cs.check_relaxed_instance_relation(&proof.cf_W_i, &proof.cf_U_i)?;
+        cf_r1cs.check_relaxed_relation(&proof.cf_W_i, &proof.cf_U_i)?;
 
         Ok(())
     }
@@ -305,7 +292,7 @@ pub mod tests {
     use super::*;
     use crate::commitment::pedersen::Pedersen;
     use crate::folding::nova::tests::test_ivc_opt;
-    use crate::frontend::tests::CubicFCircuit;
+    use crate::frontend::utils::CubicFCircuit;
     use crate::transcript::poseidon::poseidon_canonical_config;
     use ark_bn254::{Fr, G1Projective as Projective};
     use ark_grumpkin::{constraints::GVar as GVar2, Projective as Projective2};
@@ -380,11 +367,9 @@ poseidon_config.clone(),
             F_circuit,
             3,
         );
-        let (sampled_committed_instance, _) = nova
+        let (_, sampled_committed_instance) = nova
             .r1cs
-            .clone()
-            .relax()
-            .sample::<Projective, Pedersen<Projective, true>>(&nova.cs_pp, rng)
+            .sample::<Pedersen<Projective, true>>(&nova.cs_pp, rng)
             .unwrap();
 
         // proof verification fails with incorrect running instance
@@ -419,11 +404,9 @@ poseidon_config.clone(),
             F_circuit,
             3,
         );
-        let (_, sampled_committed_witness) = nova
+        let (sampled_committed_witness, _) = nova
             .r1cs
-            .clone()
-            .relax()
-            .sample::<Projective, Pedersen<Projective, true>>(&nova.cs_pp, rng)
+            .sample::<Pedersen<Projective, true>>(&nova.cs_pp, rng)
             .unwrap();
 
         // proof generation fails with incorrect running witness

@@ -1,26 +1,46 @@
-use ark_crypto_primitives::sponge::CryptographicSponge;
+use ark_crypto_primitives::sponge::{
+    constraints::CryptographicSpongeVar,
+    poseidon::{constraints::PoseidonSpongeVar, PoseidonConfig, PoseidonSponge},
+    Absorb, CryptographicSponge,
+};
 use ark_ec::CurveGroup;
+use ark_ff::PrimeField;
 use ark_poly::{univariate::DensePolynomial, EvaluationDomain, GeneralEvaluationDomain};
 use ark_r1cs_std::{
     alloc::AllocVar,
+    boolean::Boolean,
+    eq::EqGadget,
     fields::{fp::FpVar, FieldVar},
+    groups::{CurveVar, GroupOpsBounds},
     poly::polynomial::univariate::dense::DensePolynomialVar,
+    R1CSVar, ToBitsGadget, ToConstraintFieldGadget,
 };
-use ark_relations::r1cs::{ConstraintSystemRef, SynthesisError};
+use ark_relations::r1cs::{ConstraintSynthesizer, ConstraintSystemRef, SynthesisError};
+use ark_std::{fmt::Debug, marker::PhantomData, One, Zero};
 
 use super::{
     folding::lagrange_polys,
     utils::{all_powers_var, betas_star_var, exponential_powers_var},
-    CommittedInstanceVar,
+    CommittedInstance, CommittedInstanceVar, ProtoGalaxyCycleFoldConfig,
 };
 use crate::{
-    folding::circuits::nonnative::affine::NonNativeAffineVar, transcript::TranscriptVar,
+    folding::circuits::{
+        cyclefold::{
+            CycleFoldChallengeGadget, CycleFoldCommittedInstance, CycleFoldCommittedInstanceVar,
+            CycleFoldConfig, NIFSFullGadget,
+        },
+        nonnative::{affine::NonNativeAffineVar, uint::NonNativeUintVar},
+        CF1, CF2,
+    },
+    frontend::FCircuit,
+    transcript::{AbsorbNonNativeGadget, TranscriptVar},
     utils::gadgets::VectorGadget,
 };
 
 pub struct FoldingGadget {}
 
 impl FoldingGadget {
+    #[allow(clippy::type_complexity)]
     pub fn fold_committed_instance<C: CurveGroup, S: CryptographicSponge>(
         transcript: &mut impl TranscriptVar<C::ScalarField, S>,
         // running instance
@@ -30,9 +50,8 @@ vec_instances: &[CommittedInstanceVar<C>],
         // polys from P
         F_coeffs: Vec<FpVar<C::ScalarField>>,
         K_coeffs: Vec<FpVar<C::ScalarField>>,
-    ) -> Result<CommittedInstanceVar<C>, SynthesisError> {
+    ) -> Result<(CommittedInstanceVar<C>, Vec<FpVar<C::ScalarField>>), SynthesisError> {
         let t = instance.betas.len();
-        let n = F_coeffs.len();
 
         // absorb the committed instances
         transcript.absorb(instance)?;
@@ -44,7 +63,7 @@
         transcript.absorb(&F_coeffs)?;
 
         let alpha = transcript.get_challenge()?;
-        let alphas = all_powers_var(alpha.clone(), n);
+        let alphas = all_powers_var(alpha.clone(), t);
 
         // F(alpha) = e + \sum_t F_i * alpha^i
         let mut F_alpha = instance.e.clone();
@@ -88,34 +107,377 @@ .collect::<Result<Vec<_>, _>>()?;
 
         let e_star = F_alpha * &L_X_evals[0] + Z_X.evaluate(&gamma)? * K_X.evaluate(&gamma)?;
 
-        let mut u_star = &instance.u * &L_X_evals[0];
         let mut x_star = instance.x.mul_scalar(&L_X_evals[0])?;
         for i in 0..k {
-            u_star += &vec_instances[i].u * &L_X_evals[i + 1];
             x_star = x_star.add(&vec_instances[i].x.mul_scalar(&L_X_evals[i + 1])?)?;
         }
 
         // return the folded instance
-        Ok(CommittedInstanceVar {
-            betas: betas_star,
-            // phi will be computed in CycleFold
-            phi: NonNativeAffineVar::new_constant(ConstraintSystemRef::None, C::zero())?,
-            e: e_star,
-            u: u_star,
-            x: x_star,
-        })
+        Ok((
+            CommittedInstanceVar {
+                betas: betas_star,
+                // phi will be computed in CycleFold
+                phi: NonNativeAffineVar::new_constant(ConstraintSystemRef::None, C::zero())?,
+                e: e_star,
+                x: x_star,
+            },
+            L_X_evals,
+        ))
+    }
+}
+
+pub struct AugmentationGadget;
+
+impl AugmentationGadget {
+    #[allow(clippy::type_complexity)]
+    pub fn prepare_and_fold_primary<C: CurveGroup, S: CryptographicSponge>(
+        transcript: &mut impl TranscriptVar<CF1<C>, S>,
+        U: CommittedInstanceVar<C>,
+        u_phis: Vec<NonNativeAffineVar<C>>,
+        u_xs: Vec<Vec<FpVar<CF1<C>>>>,
+        new_U_phi: NonNativeAffineVar<C>,
+        F_coeffs: Vec<FpVar<CF1<C>>>,
+        K_coeffs: Vec<FpVar<CF1<C>>>,
+    ) -> Result<(CommittedInstanceVar<C>, Vec<FpVar<CF1<C>>>), SynthesisError> {
+        assert_eq!(u_phis.len(), u_xs.len());
+
+        // Prepare the incoming instances.
+        // For each instance `u`, we have `u.betas = []`, `u.e = 0`.
+        let us = u_phis
+            .into_iter()
+            .zip(u_xs)
+            .map(|(phi, x)| CommittedInstanceVar {
+                phi,
+                betas: vec![],
+                e: FpVar::zero(),
+                x,
+            })
+            .collect::<Vec<_>>();
+
+        // Fold the incoming instances `us` into the running instance `U`.
+        let (mut U, L_X_evals) =
+            FoldingGadget::fold_committed_instance(transcript, &U, &us, F_coeffs, K_coeffs)?;
+        // Notice that FoldingGadget::fold_committed_instance does not fold phi.
+        // We set `U.phi` to unconstrained witnesses `U_phi` here, whose
+        // correctness will be checked on the other curve.
+        U.phi = new_U_phi;
+
+        Ok((U, L_X_evals))
+    }
+
+    pub fn prepare_and_fold_cyclefold<
+        C1: CurveGroup<BaseField = C2::ScalarField, ScalarField = C2::BaseField>,
+        C2: CurveGroup,
+        GC2: CurveVar<C2, CF2<C2>> + ToConstraintFieldGadget<CF2<C2>>,
+        S: CryptographicSponge,
+    >(
+        transcript: &mut PoseidonSpongeVar<CF1<C1>>,
+        pp_hash: FpVar<CF1<C1>>,
+        mut cf_U: CycleFoldCommittedInstanceVar<C2, GC2>,
+        cf_u_cmWs: Vec<GC2>,
+        cf_u_xs: Vec<Vec<NonNativeUintVar<CF1<C1>>>>,
+        cf_cmTs: Vec<GC2>,
+    ) -> Result<CycleFoldCommittedInstanceVar<C2, GC2>, SynthesisError>
+    where
+        C2::BaseField: PrimeField + Absorb,
+        for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
+    {
+        assert_eq!(cf_u_cmWs.len(), cf_u_xs.len());
+        assert_eq!(cf_u_xs.len(), cf_cmTs.len());
+
+        // Fold the incoming CycleFold instances into the running CycleFold
+        // instance in a iterative way, since `NIFSFullGadget` only supports
+        // folding one incoming instance at a time.
+        for ((cmW, x), cmT) in cf_u_cmWs.into_iter().zip(cf_u_xs).zip(cf_cmTs) {
+            // Prepare the incoming CycleFold instance `cf_u` for the current
+            // iteration.
+            // For each CycleFold instance `cf_u`, we have `cf_u.cmE = 0`, and
+            // `cf_u.u = 1`.
+            let cf_u = CycleFoldCommittedInstanceVar {
+                cmE: GC2::zero(),
+                u: NonNativeUintVar::new_constant(ConstraintSystemRef::None, C1::BaseField::one())?,
+                cmW,
+                x,
+            };
+
+            let cf_r_bits = CycleFoldChallengeGadget::get_challenge_gadget(
+                transcript,
+                pp_hash.clone(),
+                cf_U.to_native_sponge_field_elements()?,
+                cf_u.clone(),
+                cmT.clone(),
+            )?;
+            // Fold the current incoming CycleFold instance `cf_u` into the
+            // running CycleFold instance `cf_U`.
+            cf_U = NIFSFullGadget::fold_committed_instance(cf_r_bits, cmT, cf_U, cf_u)?;
+        }
+
+        Ok(cf_U)
+    }
+}
+
+/// `AugmentedFCircuit` enhances the original step function `F`, so that it can
+/// be used in recursive arguments such as IVC.
+///
+/// The method for converting `F` to `AugmentedFCircuit` (`F'`) is defined in
+/// [Nova](https://eprint.iacr.org/2021/370.pdf), where `AugmentedFCircuit` not
+/// only invokes `F`, but also adds additional constraints for verifying the
+/// correct folding of primary instances (i.e., the instances over `C1`).
+/// In the paper, the primary instances are Nova's `CommittedInstance`, but we
+/// extend this method to support using ProtoGalaxy's `CommittedInstance` as
+/// primary instances.
+///
+/// Furthermore, to reduce circuit size over `C2`, we implement the constraints
+/// defined in [CycleFold](https://eprint.iacr.org/2023/1192.pdf). These extra
+/// constraints verify the correct folding of CycleFold instances.
+#[derive(Debug, Clone)]
+pub struct AugmentedFCircuit<
+    C1: CurveGroup,
+    C2: CurveGroup,
+    GC2: CurveVar<C2, CF2<C2>>,
+    FC: FCircuit<CF1<C1>>,
+> {
+    pub(super) _gc2: PhantomData<GC2>,
+    pub(super) poseidon_config: PoseidonConfig<CF1<C1>>,
+    pub(super) pp_hash: CF1<C1>,
+    pub(super) i: CF1<C1>,
+    pub(super) i_usize: usize,
+    pub(super) z_0: Vec<CF1<C1>>,
+    pub(super) z_i: Vec<CF1<C1>>,
+    pub(super) external_inputs: Vec<CF1<C1>>,
+    pub(super) F: FC, // F circuit
+    pub(super) u_i_phi: C1,
+    pub(super) U_i: CommittedInstance<C1>,
+    pub(super) U_i1_phi: C1,
+    pub(super) F_coeffs: Vec<CF1<C1>>,
+    pub(super) K_coeffs: Vec<CF1<C1>>,
+    pub(super) x: Option<CF1<C1>>, // public input (u_{i+1}.x[0])
+
+    pub(super) phi_stars: Vec<C1>,
+
+    pub(super) cf1_u_i_cmW: C2,                        // input
+    pub(super) cf2_u_i_cmW: C2,                        // input
+    pub(super) cf_U_i: CycleFoldCommittedInstance<C2>, // input
+    pub(super) cf1_cmT: C2,
+    pub(super) cf2_cmT: C2,
+    pub(super) cf_x: Option<CF1<C1>>, // public input (u_{i+1}.x[1])
+}
+
+impl<C1: CurveGroup, C2: CurveGroup, GC2: CurveVar<C2, CF2<C2>>, FC: FCircuit<CF1<C1>>>
+    AugmentedFCircuit<C1, C2, GC2, FC>
+where
+    for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
+{
+    pub fn empty(
+        poseidon_config: &PoseidonConfig<CF1<C1>>,
+        F_circuit: FC,
+        t: usize,
+        d: usize,
+        k: usize,
+    ) -> Self {
+        let u_dummy = CommittedInstance::dummy_running(2, t);
+        let cf_u_dummy =
+            CycleFoldCommittedInstance::dummy(ProtoGalaxyCycleFoldConfig::<C1>::IO_LEN);
+
+        Self {
+            _gc2: PhantomData,
+            poseidon_config: poseidon_config.clone(),
+            pp_hash: CF1::<C1>::zero(),
+            i: CF1::<C1>::zero(),
+            i_usize: 0,
+            z_0: vec![CF1::<C1>::zero(); F_circuit.state_len()],
+            z_i: vec![CF1::<C1>::zero(); F_circuit.state_len()],
+            external_inputs: vec![CF1::<C1>::zero(); F_circuit.external_inputs_len()],
+            u_i_phi: C1::zero(),
+            U_i: u_dummy,
+            U_i1_phi: C1::zero(),
+            F_coeffs: vec![CF1::<C1>::zero(); t],
+            K_coeffs: vec![CF1::<C1>::zero(); d * k + 1],
+            phi_stars: vec![C1::zero(); k],
+            F: F_circuit,
+            x: None,
+            // cyclefold values
+            cf1_u_i_cmW: C2::zero(),
+            cf2_u_i_cmW: C2::zero(),
+            cf_U_i: cf_u_dummy,
+            cf1_cmT: C2::zero(),
+            cf2_cmT: C2::zero(),
+            cf_x: None,
+        }
+    }
+}
+
+impl<C1, C2, GC2, FC> ConstraintSynthesizer<CF1<C1>> for AugmentedFCircuit<C1, C2, GC2, FC>
+where
+    C1: CurveGroup<BaseField = C2::ScalarField, ScalarField = C2::BaseField>,
+    C2: CurveGroup,
+    GC2: CurveVar<C2, CF2<C2>> + ToConstraintFieldGadget<CF2<C2>>,
+    FC: FCircuit<CF1<C1>>,
+    C2::BaseField: PrimeField + Absorb,
+    for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
+{
+    fn generate_constraints(self, cs: ConstraintSystemRef<CF1<C1>>) -> Result<(), SynthesisError> {
+        let pp_hash = FpVar::<CF1<C1>>::new_witness(cs.clone(), || Ok(self.pp_hash))?;
+        let i = FpVar::<CF1<C1>>::new_witness(cs.clone(), || Ok(self.i))?;
+        let z_0 = Vec::<FpVar<CF1<C1>>>::new_witness(cs.clone(), || Ok(self.z_0))?;
+        let z_i = Vec::<FpVar<CF1<C1>>>::new_witness(cs.clone(), || Ok(self.z_i))?;
+        let external_inputs =
+            Vec::<FpVar<CF1<C1>>>::new_witness(cs.clone(), || Ok(self.external_inputs))?;
+
+        let u_dummy = CommittedInstance::<C1>::dummy_running(2, self.U_i.betas.len());
+        let U_i = CommittedInstanceVar::<C1>::new_witness(cs.clone(), || Ok(self.U_i))?;
+        let u_i_phi = NonNativeAffineVar::new_witness(cs.clone(), || Ok(self.u_i_phi))?;
+        let U_i1_phi = NonNativeAffineVar::new_witness(cs.clone(), || Ok(self.U_i1_phi))?;
+        let phi_stars =
+            Vec::<NonNativeAffineVar<C1>>::new_witness(cs.clone(), || Ok(self.phi_stars))?;
+
+        let cf_u_dummy =
+            CycleFoldCommittedInstance::dummy(ProtoGalaxyCycleFoldConfig::<C1>::IO_LEN);
+        let cf_U_i =
+            CycleFoldCommittedInstanceVar::<C2, GC2>::new_witness(cs.clone(), || Ok(self.cf_U_i))?;
+        let cf1_cmT = GC2::new_witness(cs.clone(), || Ok(self.cf1_cmT))?;
+        let cf2_cmT = GC2::new_witness(cs.clone(), || Ok(self.cf2_cmT))?;
+
+        let F_coeffs = Vec::new_witness(cs.clone(), || Ok(self.F_coeffs))?;
+        let K_coeffs = Vec::new_witness(cs.clone(), || Ok(self.K_coeffs))?;
+
+        // `sponge` is for digest computation.
+        let sponge = PoseidonSpongeVar::<C1::ScalarField>::new(cs.clone(), &self.poseidon_config);
+        // `transcript` is for challenge generation.
+        let mut transcript = sponge.clone();
+
+        // get z_{i+1} from the F circuit
+        let i_usize = self.i_usize;
+        let z_i1 =
+            self.F
+                .generate_step_constraints(cs.clone(), i_usize, z_i.clone(), external_inputs)?;
+
+        let is_basecase = i.is_zero()?;
+
+        // Primary Part
+        // P.1. Compute u_i.x
+        // u_i.x[0] = H(i, z_0, z_i, U_i)
+        let (u_i_x, _) = U_i.clone().hash(
+            &sponge,
+            pp_hash.clone(),
+            i.clone(),
+            z_0.clone(),
+            z_i.clone(),
+        )?;
+        // u_i.x[1] = H(cf_U_i)
+        let (cf_u_i_x, _) = cf_U_i.clone().hash(&sponge, pp_hash.clone())?;
+
+        // P.2. Prepare incoming primary instances
+        // P.3. Fold incoming primary instances into the running instance
+        let (U_i1, r) = AugmentationGadget::prepare_and_fold_primary(
+            &mut transcript,
+            U_i.clone(),
+            vec![u_i_phi.clone()],
+            vec![vec![u_i_x, cf_u_i_x]],
+            U_i1_phi,
+            F_coeffs,
+            K_coeffs,
+        )?;
+
+        // P.4.a compute and check the first output of F'
+        // Base case: u_{i+1}.x[0] == H((i+1, z_0, z_{i+1}, U_{\bot})
+        // Non-base case: u_{i+1}.x[0] == H((i+1, z_0, z_{i+1}, U_{i+1})
+        let (u_i1_x, _) = U_i1.clone().hash(
+            &sponge,
+            pp_hash.clone(),
+            i + FpVar::<CF1<C1>>::one(),
+            z_0.clone(),
+            z_i1.clone(),
+        )?;
+        let (u_i1_x_base, _) = CommittedInstanceVar::new_constant(cs.clone(), u_dummy)?.hash(
+            &sponge,
+            pp_hash.clone(),
+            FpVar::<CF1<C1>>::one(),
+            z_0.clone(),
+            z_i1.clone(),
+        )?;
+        let x = FpVar::new_input(cs.clone(), || Ok(self.x.unwrap_or(u_i1_x_base.value()?)))?;
+        x.enforce_equal(&is_basecase.select(&u_i1_x_base, &u_i1_x)?)?;
+
+        // CycleFold part
+        // C.1. Compute cf1_u_i.x and cf2_u_i.x
+        let mut r0_bits = r[0].to_bits_le()?;
+        let mut r1_bits = r[1].to_bits_le()?;
+        r0_bits.resize(C1::ScalarField::MODULUS_BIT_SIZE as usize, Boolean::FALSE);
+        r1_bits.resize(C1::ScalarField::MODULUS_BIT_SIZE as usize, Boolean::FALSE);
+        let cf1_x = [
+            r0_bits
+                .chunks(C1::BaseField::MODULUS_BIT_SIZE as usize - 1)
+                .map(|bits| {
+                    let mut bits = bits.to_vec();
+                    bits.resize(C1::BaseField::MODULUS_BIT_SIZE as usize, Boolean::FALSE);
+                    NonNativeUintVar::from(&bits)
+                })
+                .collect::<Vec<_>>(),
+            vec![
+                NonNativeUintVar::new_constant(cs.clone(), C1::BaseField::zero())?,
+                NonNativeUintVar::new_constant(cs.clone(), C1::BaseField::zero())?,
+                U_i.phi.x.clone(),
+                U_i.phi.y.clone(),
+                phi_stars[0].x.clone(),
+                phi_stars[0].y.clone(),
+            ],
+        ]
+        .concat();
+        let cf2_x = [
+            r1_bits
+                .chunks(C1::BaseField::MODULUS_BIT_SIZE as usize - 1)
+                .map(|bits| {
+                    let mut bits = bits.to_vec();
+                    bits.resize(C1::BaseField::MODULUS_BIT_SIZE as usize, Boolean::FALSE);
+                    NonNativeUintVar::from(&bits)
+                })
+                .collect::<Vec<_>>(),
+            vec![
+                phi_stars[0].x.clone(),
+                phi_stars[0].y.clone(),
+                u_i_phi.x.clone(),
+                u_i_phi.y.clone(),
+                U_i1.phi.x.clone(),
+                U_i1.phi.y.clone(),
+            ],
+        ]
+        .concat();
+
+        // C.2. Prepare incoming CycleFold instances
+        // C.3. Fold incoming CycleFold instances into the running instance
+        let cf_U_i1 =
+            AugmentationGadget::prepare_and_fold_cyclefold::<C1, C2, GC2, PoseidonSponge<CF1<C1>>>(
+                &mut transcript,
+                pp_hash.clone(),
+                cf_U_i,
+                vec![
+                    GC2::new_witness(cs.clone(), || Ok(self.cf1_u_i_cmW))?,
+                    GC2::new_witness(cs.clone(), || Ok(self.cf2_u_i_cmW))?,
+                ],
+                vec![cf1_x, cf2_x],
+                vec![cf1_cmT, cf2_cmT],
+            )?;
+
+        // Back to Primary Part
+        // P.4.b compute and check the second output of F'
+        // Base case: u_{i+1}.x[1] == H(cf_U_{\bot})
+        // Non-base case: u_{i+1}.x[1] == H(cf_U_{i+1})
+        let (cf_u_i1_x, _) = cf_U_i1.clone().hash(&sponge, pp_hash.clone())?;
+        let (cf_u_i1_x_base, _) =
+            CycleFoldCommittedInstanceVar::<C2, GC2>::new_constant(cs.clone(), cf_u_dummy)?
+                .hash(&sponge, pp_hash.clone())?;
+        let cf_x = FpVar::new_input(cs.clone(), || {
+            Ok(self.cf_x.unwrap_or(cf_u_i1_x_base.value()?))
+        })?;
+        cf_x.enforce_equal(&is_basecase.select(&cf_u_i1_x_base, &cf_u_i1_x)?)?;
+
+        Ok(())
     }
 }
 
 #[cfg(test)]
 mod tests {
-    use ark_crypto_primitives::sponge::{
-        constraints::CryptographicSpongeVar,
-        poseidon::{constraints::PoseidonSpongeVar, PoseidonSponge},
-    };
-    use ark_pallas::{Fr, Projective};
-    use ark_r1cs_std::R1CSVar;
-    use ark_relations::r1cs::ConstraintSystem;
     use std::error::Error;
 
     use super::*;
@@ -125,8 +487,11 @@ folding::protogalaxy::folding::{tests::prepare_inputs, Folding},
         transcript::poseidon::poseidon_canonical_config,
     };
 
+    use ark_bn254::{Fr, G1Projective as Projective};
+    use ark_relations::r1cs::ConstraintSystem;
+
     #[test]
-    fn test_fold_gadget() -> Result<(), Box<dyn Error>> {
+    fn test_folding_gadget() -> Result<(), Box<dyn Error>> {
         let k = 7;
         let (witness, instance, witnesses, instances) = prepare_inputs(k);
         let r1cs = get_test_r1cs::<Fr>();
@@ -136,7 +501,7 @@ let poseidon_config = poseidon_canonical_config::<Fr>();
         let mut transcript_p = PoseidonSponge::new(&poseidon_config);
         let mut transcript_v = PoseidonSponge::new(&poseidon_config);
 
-        let (_, _, F_coeffs, K_coeffs) = Folding::<Projective>::prove(
+        let (_, _, F_coeffs, K_coeffs, _, _) = Folding::<Projective>::prove(
             &mut transcript_p,
             &r1cs,
             &instance,
@@ -147,7 +512,6 @@ )?;
 
         let folded_instance = Folding::<Projective>::verify(
             &mut transcript_v,
-            &r1cs,
             &instance,
             &instances,
             F_coeffs.clone(),
@@ -161,7 +525,7 @@ let instances_var = Vec::new_witness(cs.clone(), || Ok(instances))?;
         let F_coeffs_var = Vec::new_witness(cs.clone(), || Ok(F_coeffs))?;
         let K_coeffs_var = Vec::new_witness(cs.clone(), || Ok(K_coeffs))?;
 
-        let folded_instance_var = FoldingGadget::fold_committed_instance(
+        let (folded_instance_var, _) = FoldingGadget::fold_committed_instance(
             &mut transcript_var,
             &instance_var,
             &instances_var,
@@ -170,7 +534,6 @@ K_coeffs_var,
         )?;
         assert_eq!(folded_instance.betas, folded_instance_var.betas.value()?);
         assert_eq!(folded_instance.e, folded_instance_var.e.value()?);
-        assert_eq!(folded_instance.u, folded_instance_var.u.value()?);
         assert_eq!(folded_instance.x, folded_instance_var.x.value()?);
         assert!(cs.is_satisfied()?);
 

@@ -6,18 +6,19 @@ use ark_poly::{
     univariate::{DensePolynomial, SparsePolynomial},
     DenseUVPolynomial, EvaluationDomain, Evaluations, GeneralEvaluationDomain, Polynomial,
 };
-use ark_std::{cfg_into_iter, log2, Zero};
-use rayon::iter::{IntoParallelIterator, ParallelIterator};
+use ark_std::{cfg_into_iter, log2, One, Zero};
+use rayon::prelude::*;
 use std::marker::PhantomData;
 
-use super::utils::{all_powers, betas_star, exponential_powers};
+use super::utils::{all_powers, betas_star, exponential_powers, pow_i};
 use super::ProtoGalaxyError;
 use super::{CommittedInstance, Witness};
 
-use crate::arith::r1cs::R1CS;
+#[cfg(test)]
+use crate::arith::r1cs::RelaxedR1CS;
+use crate::arith::{r1cs::R1CS, Arith};
 use crate::transcript::Transcript;
 use crate::utils::vec::*;
-use crate::utils::virtual_polynomial::bit_decompose;
 use crate::Error;
 
 #[derive(Clone, Debug)]
@@ -48,6 +49,8 @@ CommittedInstance<C>,
             Witness<C::ScalarField>,
             Vec<C::ScalarField>, // F_X coeffs
             Vec<C::ScalarField>, // K_X coeffs
+            Vec<C::ScalarField>, // L_X evals
+            Vec<C>,              // phi_stars
         ),
         Error,
     > {
@@ -63,19 +66,25 @@ let d = 2; // for the moment hardcoded to 2 since it only supports R1CS
         let k = vec_instances.len();
         let t = instance.betas.len();
         let n = r1cs.A.n_cols;
+        let m = r1cs.A.n_rows;
 
-        let z = [vec![instance.u], instance.x.clone(), w.w.clone()].concat();
+        let z = [vec![C::ScalarField::one()], instance.x.clone(), w.w.clone()].concat();
 
         if z.len() != n {
             return Err(Error::NotSameLength(
                 "z.len()".to_string(),
                 z.len(),
-                "n".to_string(),
+                "number of variables in R1CS".to_string(), // hardcoded to R1CS
                 n,
             ));
         }
-        if log2(n) as usize != t {
-            return Err(Error::NotEqual);
+        if log2(m) as usize != t {
+            return Err(Error::NotSameLength(
+                "log2(number of constraints in R1CS)".to_string(),
+                log2(m) as usize,
+                "instance.betas.len()".to_string(),
+                t,
+            ));
         }
         if !(k + 1).is_power_of_two() {
             return Err(Error::ProtoGalaxy(ProtoGalaxyError::WrongNumInstances(k)));
@@ -88,13 +97,24 @@
         let delta = transcript.get_challenge();
         let deltas = exponential_powers(delta, t);
 
-        let f_z = eval_f(r1cs, &z)?;
+        let mut f_z = r1cs.eval_relation(&z)?;
+        if f_z.len() != m {
+            return Err(Error::NotSameLength(
+                "number of constraints in R1CS".to_string(),
+                m,
+                "f_z.len()".to_string(),
+                f_z.len(),
+            ));
+        }
+        f_z.resize(1 << t, C::ScalarField::zero());
 
         // F(X)
         let F_X: SparsePolynomial<C::ScalarField> =
             calc_f_from_btree(&f_z, &instance.betas, &deltas).expect("Error calculating F[x]");
         let F_X_dense = DensePolynomial::from(F_X.clone());
-        transcript.absorb(&F_X_dense.coeffs);
+        let mut F_coeffs = F_X_dense.coeffs;
+        F_coeffs.resize(t, C::ScalarField::zero());
+        transcript.absorb(&F_coeffs);
 
         let alpha = transcript.get_challenge();
 
@@ -107,16 +127,14 @@
         // sanity check: check that the new randomized instance (the original instance but with
         // 'refreshed' randomness) satisfies the relation.
         #[cfg(test)]
-        tests::check_instance(
-            r1cs,
+        r1cs.check_relaxed_relation(
+            w,
             &CommittedInstance {
                 phi: instance.phi,
                 betas: betas_star.clone(),
                 e: F_alpha,
-                u: instance.u,
                 x: instance.x.clone(),
             },
-            w,
         )?;
 
         let zs: Vec<Vec<C::ScalarField>> = std::iter::once(z.clone())
@@ -125,12 +143,12 @@ vec_w
                     .iter()
                     .zip(vec_instances)
                     .map(|(wj, uj)| {
-                        let zj = [vec![uj.u], uj.x.clone(), wj.w.clone()].concat();
+                        let zj = [vec![C::ScalarField::one()], uj.x.clone(), wj.w.clone()].concat();
                         if zj.len() != n {
                             return Err(Error::NotSameLength(
                                 "zj.len()".to_string(),
                                 zj.len(),
-                                "n".to_string(),
+                                "number of variables in R1CS".to_string(),
                                 n,
                             ));
                         }
@@ -153,26 +171,19 @@ for (hi, h) in G_domain.elements().enumerate() {
             // each iteration evaluates G(h)
             // inner = L_0(x) * z + \sum_k L_i(x) * z_j
             let mut inner: Vec<C::ScalarField> = vec![C::ScalarField::zero(); zs[0].len()];
-            for (i, z) in zs.iter().enumerate() {
+            for (z, L) in zs.iter().zip(&L_X) {
                 // Li_z_h = (Li(X)*zj)(h) = Li(h) * zj
-                let mut Liz_h: Vec<C::ScalarField> = vec![C::ScalarField::zero(); z.len()];
+                let Lh = L.evaluate(&h);
                 for (j, zj) in z.iter().enumerate() {
-                    Liz_h[j] = (&L_X[i] * *zj).evaluate(&h);
-                }
-
-                for j in 0..inner.len() {
-                    inner[j] += Liz_h[j];
+                    inner[j] += Lh * zj;
                 }
             }
-            let f_ev = eval_f(r1cs, &inner)?;
+            let f_ev = r1cs.eval_relation(&inner)?;
 
-            let mut Gsum = C::ScalarField::zero();
-            for (i, f_ev_i) in f_ev.iter().enumerate() {
-                let pow_i_betas = pow_i(i, &betas_star);
-                let curr = pow_i_betas * f_ev_i;
-                Gsum += curr;
-            }
-            G_evals[hi] = Gsum;
+            G_evals[hi] = cfg_into_iter!(f_ev)
+                .enumerate()
+                .map(|(i, f_ev_i)| pow_i(i, &betas_star) * f_ev_i)
+                .sum();
         }
         let G_X: DensePolynomial<C::ScalarField> =
             Evaluations::<C::ScalarField>::from_vec_and_domain(G_evals, G_domain).interpolate();
@@ -190,7 +201,9 @@ if !remainder.is_zero() {
             return Err(Error::ProtoGalaxy(ProtoGalaxyError::RemainderNotZero));
         }
 
-        transcript.absorb(&K_X.coeffs);
+        let mut K_coeffs = K_X.coeffs.clone();
+        K_coeffs.resize(d * k + 1, C::ScalarField::zero());
+        transcript.absorb(&K_coeffs);
 
         let gamma = transcript.get_challenge();
 
@@ -199,18 +212,19 @@ .iter()
             .take(k + 1)
             .map(|L| L.evaluate(&gamma))
             .collect::<Vec<_>>();
+
+        let mut phi_stars = vec![];
 
         let e_star = F_alpha * L_X_evals[0] + Z_X.evaluate(&gamma) * K_X.evaluate(&gamma);
         let mut w_star = vec_scalar_mul(&w.w, &L_X_evals[0]);
         let mut r_w_star = w.r_w * L_X_evals[0];
         let mut phi_star = instance.phi * L_X_evals[0];
-        let mut u_star = instance.u * L_X_evals[0];
         let mut x_star = vec_scalar_mul(&instance.x, &L_X_evals[0]);
         for i in 0..k {
             w_star = vec_add(&w_star, &vec_scalar_mul(&vec_w[i].w, &L_X_evals[i + 1]))?;
             r_w_star += vec_w[i].r_w * L_X_evals[i + 1];
+            phi_stars.push(phi_star); // Push before updating. We don't need the last one
             phi_star += vec_instances[i].phi * L_X_evals[i + 1];
-            u_star += vec_instances[i].u * L_X_evals[i + 1];
             x_star = vec_add(
                 &x_star,
                 &vec_scalar_mul(&vec_instances[i].x, &L_X_evals[i + 1]),
@@ -222,22 +236,22 @@ CommittedInstance {
                 betas: betas_star,
                 phi: phi_star,
                 e: e_star,
-                u: u_star,
                 x: x_star,
             },
             Witness {
                 w: w_star,
                 r_w: r_w_star,
             },
-            F_X_dense.coeffs,
-            K_X.coeffs,
+            F_coeffs,
+            K_coeffs,
+            L_X_evals,
+            phi_stars,
         ))
     }
 
     /// implements the non-interactive Verifier from the folding scheme described in section 4
     pub fn verify(
         transcript: &mut impl Transcript<C::ScalarField>,
-        r1cs: &R1CS<C::ScalarField>,
         // running instance
         instance: &CommittedInstance<C>,
         // incoming instances
@@ -247,7 +261,6 @@ F_coeffs: Vec<C::ScalarField>,
         K_coeffs: Vec<C::ScalarField>,
     ) -> Result<CommittedInstance<C>, Error> {
         let t = instance.betas.len();
-        let n = r1cs.A.n_cols;
 
         // absorb the committed instances
         transcript.absorb(instance);
@@ -259,7 +272,7 @@
         transcript.absorb(&F_coeffs);
 
         let alpha = transcript.get_challenge();
-        let alphas = all_powers(alpha, n);
+        let alphas = all_powers(alpha, t);
 
         // F(alpha) = e + \sum_t F_i * alpha^i
         let mut F_alpha = instance.e;
@@ -268,6 +281,8 @@ F_alpha += *F_i * alphas[i + 1];
         }
 
         let betas_star = betas_star(&instance.betas, &deltas, alpha);
+
+        transcript.absorb(&K_coeffs);
 
         let k = vec_instances.len();
         let H =
@@ -277,8 +292,6 @@ let Z_X: DensePolynomial<C::ScalarField> = H.vanishing_polynomial().into();
         let K_X: DensePolynomial<C::ScalarField> =
             DensePolynomial::<C::ScalarField>::from_coefficients_vec(K_coeffs);
 
-        transcript.absorb(&K_X.coeffs);
-
         let gamma = transcript.get_challenge();
 
         let L_X_evals = L_X
@@ -290,11 +303,9 @@
         let e_star = F_alpha * L_X_evals[0] + Z_X.evaluate(&gamma) * K_X.evaluate(&gamma);
 
         let mut phi_star = instance.phi * L_X_evals[0];
-        let mut u_star = instance.u * L_X_evals[0];
         let mut x_star = vec_scalar_mul(&instance.x, &L_X_evals[0]);
         for i in 0..k {
             phi_star += vec_instances[i].phi * L_X_evals[i + 1];
-            u_star += vec_instances[i].u * L_X_evals[i + 1];
             x_star = vec_add(
                 &x_star,
                 &vec_scalar_mul(&vec_instances[i].x, &L_X_evals[i + 1]),
@@ -306,28 +317,9 @@ Ok(CommittedInstance {
             betas: betas_star,
             phi: phi_star,
             e: e_star,
-            u: u_star,
             x: x_star,
         })
     }
-}
-
-// naive impl of pow_i for betas, assuming that betas=(b, b^2, b^4, ..., b^{2^{t-1}})
-fn pow_i<F: PrimeField>(i: usize, betas: &[F]) -> F {
-    // WIP check if makes more sense to do it with ifs instead of arithmetic
-
-    let n = 2_u64.pow(betas.len() as u32);
-    let b = bit_decompose(i as u64, n as usize);
-
-    let mut r: F = F::one();
-    for (j, beta_j) in betas.iter().enumerate() {
-        let mut b_j = F::zero();
-        if b[j] {
-            b_j = F::one();
-        }
-        r *= (F::one() - b_j) + b_j * beta_j;
-    }
-    r
 }
 
 /// calculates F[x] using the optimized binary-tree technique
@@ -394,16 +386,6 @@ }
     lagrange_polynomials
 }
 
-// f(w) in R1CS context. For the moment we use R1CS, in the future we will abstract this with a
-// trait
-fn eval_f<F: PrimeField>(r1cs: &R1CS<F>, z: &[F]) -> Result<Vec<F>, Error> {
-    let Az = mat_vec_mul(&r1cs.A, z)?;
-    let Bz = mat_vec_mul(&r1cs.B, z)?;
-    let Cz = mat_vec_mul(&r1cs.C, z)?;
-    let AzBz = hadamard(&Az, &Bz)?;
-    vec_sub(&AzBz, &Cz)
-}
-
 #[cfg(test)]
 pub mod tests {
     use super::*;
@@ -412,38 +394,10 @@ use ark_crypto_primitives::sponge::CryptographicSponge;
     use ark_pallas::{Fr, Projective};
     use ark_std::{rand::Rng, UniformRand};
 
-    use crate::arith::r1cs::tests::{get_test_r1cs, get_test_z, get_test_z_split};
+    use crate::arith::r1cs::tests::{get_test_r1cs, get_test_z_split};
     use crate::commitment::{pedersen::Pedersen, CommitmentScheme};
     use crate::transcript::poseidon::poseidon_canonical_config;
 
-    pub(crate) fn check_instance<C: CurveGroup>(
-        r1cs: &R1CS<C::ScalarField>,
-        instance: &CommittedInstance<C>,
-        w: &Witness<C::ScalarField>,
-    ) -> Result<(), Error> {
-        let z = [vec![instance.u], instance.x.clone(), w.w.clone()].concat();
-
-        if instance.betas.len() != log2(z.len()) as usize {
-            return Err(Error::NotSameLength(
-                "instance.betas.len()".to_string(),
-                instance.betas.len(),
-                "log2(z.len())".to_string(),
-                log2(z.len()) as usize,
-            ));
-        }
-
-        let f_z = eval_f(r1cs, &z)?; // f(z)
-
-        let mut r = C::ScalarField::zero();
-        for (i, f_z_i) in f_z.iter().enumerate() {
-            r += pow_i(i, &instance.betas) * f_z_i;
-        }
-        if instance.e == r {
-            return Ok(());
-        }
-        Err(Error::NotSatisfied)
-    }
-
     #[test]
     fn test_pow_i() {
         let mut rng = ark_std::test_rng();
@@ -459,76 +413,54 @@ assert_eq!(pow_i(i, &betas), not_betas[i]);
         }
     }
 
-    #[test]
-    fn test_eval_f() {
-        let mut rng = ark_std::test_rng();
-        let r1cs = get_test_r1cs::<Fr>();
-        let mut z = get_test_z::<Fr>(rng.gen::<u16>() as usize);
-
-        let f_w = eval_f(&r1cs, &z).unwrap();
-        assert!(is_zero_vec(&f_w));
-
-        z[1] = Fr::from(111);
-        let f_w = eval_f(&r1cs, &z).unwrap();
-        assert!(!is_zero_vec(&f_w));
-    }
-
     // k represents the number of instances to be fold, apart from the running instance
     #[allow(clippy::type_complexity)]
-    pub fn prepare_inputs(
+    pub fn prepare_inputs<C: CurveGroup>(
         k: usize,
     ) -> (
-        Witness<Fr>,
-        CommittedInstance<Projective>,
-        Vec<Witness<Fr>>,
-        Vec<CommittedInstance<Projective>>,
+        Witness<C::ScalarField>,
+        CommittedInstance<C>,
+        Vec<Witness<C::ScalarField>>,
+        Vec<CommittedInstance<C>>,
     ) {
         let mut rng = ark_std::test_rng();
 
-        let (u, x, w) = get_test_z_split::<Fr>(rng.gen::<u16>() as usize);
+        let (_, x, w) = get_test_z_split::<C::ScalarField>(rng.gen::<u16>() as usize);
 
-        let (pedersen_params, _) = Pedersen::<Projective>::setup(&mut rng, w.len()).unwrap();
+        let (pedersen_params, _) = Pedersen::<C>::setup(&mut rng, w.len()).unwrap();
 
-        let n = 1 + x.len() + w.len();
-        let t = log2(n) as usize;
+        let t = log2(get_test_r1cs::<C::ScalarField>().A.n_rows) as usize;
 
-        let beta = Fr::rand(&mut rng);
+        let beta = C::ScalarField::rand(&mut rng);
         let betas = exponential_powers(beta, t);
 
-        let witness = Witness::<Fr> {
+        let witness = Witness::<C::ScalarField> {
             w,
-            r_w: Fr::rand(&mut rng),
+            r_w: C::ScalarField::zero(),
         };
-        let phi = Pedersen::<Projective, true>::commit(&pedersen_params, &witness.w, &witness.r_w)
-            .unwrap();
-        let instance = CommittedInstance::<Projective> {
+        let phi = Pedersen::<C>::commit(&pedersen_params, &witness.w, &witness.r_w).unwrap();
+        let instance = CommittedInstance::<C> {
             phi,
             betas: betas.clone(),
-            e: Fr::zero(),
-            u,
+            e: C::ScalarField::zero(),
             x,
         };
         // same for the other instances
-        let mut witnesses: Vec<Witness<Fr>> = Vec::new();
-        let mut instances: Vec<CommittedInstance<Projective>> = Vec::new();
+        let mut witnesses: Vec<Witness<C::ScalarField>> = Vec::new();
+        let mut instances: Vec<CommittedInstance<C>> = Vec::new();
         #[allow(clippy::needless_range_loop)]
         for _ in 0..k {
-            let (u_i, x_i, w_i) = get_test_z_split::<Fr>(rng.gen::<u16>() as usize);
-            let witness_i = Witness::<Fr> {
+            let (_, x_i, w_i) = get_test_z_split::<C::ScalarField>(rng.gen::<u16>() as usize);
+            let witness_i = Witness::<C::ScalarField> {
                 w: w_i,
-                r_w: Fr::rand(&mut rng),
+                r_w: C::ScalarField::zero(),
             };
-            let phi_i = Pedersen::<Projective, true>::commit(
-                &pedersen_params,
-                &witness_i.w,
-                &witness_i.r_w,
-            )
-            .unwrap();
-            let instance_i = CommittedInstance::<Projective> {
+            let phi_i =
+                Pedersen::<C>::commit(&pedersen_params, &witness_i.w, &witness_i.r_w).unwrap();
+            let instance_i = CommittedInstance::<C> {
                 phi: phi_i,
                 betas: vec![],
-                e: Fr::zero(),
-                u: u_i,
+                e: C::ScalarField::zero(),
                 x: x_i,
             };
             witnesses.push(witness_i);
@@ -549,20 +481,20 @@ let poseidon_config = poseidon_canonical_config::<Fr>();
         let mut transcript_p = PoseidonSponge::<Fr>::new(&poseidon_config);
         let mut transcript_v = PoseidonSponge::<Fr>::new(&poseidon_config);
 
-        let (folded_instance, folded_witness, F_coeffs, K_coeffs) = Folding::<Projective>::prove(
-            &mut transcript_p,
-            &r1cs,
-            &instance,
-            &witness,
-            &instances,
-            &witnesses,
-        )
-        .unwrap();
+        let (folded_instance, folded_witness, F_coeffs, K_coeffs, _, _) =
+            Folding::<Projective>::prove(
+                &mut transcript_p,
+                &r1cs,
+                &instance,
+                &witness,
+                &instances,
+                &witnesses,
+            )
+            .unwrap();
 
         // verifier
         let folded_instance_v = Folding::<Projective>::verify(
             &mut transcript_v,
-            &r1cs,
             &instance,
             &instances,
             F_coeffs,
@@ -577,7 +509,8 @@ assert_eq!(folded_instance.e, folded_instance_v.e);
         assert!(!folded_instance.e.is_zero());
 
         // check that the folded instance satisfies the relation
-        check_instance(&r1cs, &folded_instance, &folded_witness).unwrap();
+        r1cs.check_relaxed_relation(&folded_witness, &folded_instance)
+            .unwrap();
     }
 
     #[test]
@@ -598,7 +531,7 @@ for _ in 0..num_iters {
             // generate the instances to be fold
             let (_, _, witnesses, instances) = prepare_inputs(k);
 
-            let (folded_instance, folded_witness, F_coeffs, K_coeffs) =
+            let (folded_instance, folded_witness, F_coeffs, K_coeffs, _, _) =
                 Folding::<Projective>::prove(
                     &mut transcript_p,
                     &r1cs,
@@ -612,7 +545,6 @@
             // verifier
             let folded_instance_v = Folding::<Projective>::verify(
                 &mut transcript_v,
-                &r1cs,
                 &running_instance,
                 &instances,
                 F_coeffs,
@@ -626,7 +558,8 @@
             assert!(!folded_instance.e.is_zero());
 
             // check that the folded instance satisfies the relation
-            check_instance(&r1cs, &folded_instance, &folded_witness).unwrap();
+            r1cs.check_relaxed_relation(&folded_witness, &folded_instance)
+                .unwrap();
 
             running_witness = folded_witness;
             running_instance = folded_instance;

@@ -1,37 +1,120 @@
-use std::borrow::Borrow;
-
 /// Implements the scheme described in [ProtoGalaxy](https://eprint.iacr.org/2023/1106.pdf)
-use ark_ec::CurveGroup;
-use ark_ff::PrimeField;
+use ark_crypto_primitives::sponge::{
+    constraints::{AbsorbGadget, CryptographicSpongeVar},
+    poseidon::{constraints::PoseidonSpongeVar, PoseidonConfig, PoseidonSponge},
+    Absorb, CryptographicSponge,
+};
+use ark_ec::{CurveGroup, Group};
+use ark_ff::{BigInteger, PrimeField};
 use ark_r1cs_std::{
     alloc::{AllocVar, AllocationMode},
     fields::fp::FpVar,
+    groups::{CurveVar, GroupOpsBounds},
+    R1CSVar, ToConstraintFieldGadget,
 };
-use ark_relations::r1cs::{Namespace, SynthesisError};
-use thiserror::Error;
+use ark_relations::r1cs::{
+    ConstraintSynthesizer, ConstraintSystem, ConstraintSystemRef, Namespace, SynthesisError,
+};
+use ark_std::{
+    borrow::Borrow, cmp::max, fmt::Debug, log2, marker::PhantomData, rand::RngCore, One, Zero,
+};
+use num_bigint::BigUint;
 
-use super::circuits::nonnative::affine::NonNativeAffineVar;
+use crate::{
+    arith::r1cs::{extract_r1cs, extract_w_x, RelaxedR1CS, R1CS},
+    commitment::CommitmentScheme,
+    folding::circuits::{
+        cyclefold::{
+            fold_cyclefold_circuit, CycleFoldCircuit, CycleFoldCommittedInstance, CycleFoldConfig,
+            CycleFoldWitness,
+        },
+        nonnative::affine::NonNativeAffineVar,
+        CF1, CF2,
+    },
+    frontend::{utils::DummyCircuit, FCircuit},
+    utils::{get_cm_coordinates, pp_hash},
+    Error, FoldingScheme,
+};
 
 pub mod circuits;
 pub mod folding;
 pub mod traits;
 pub(crate) mod utils;
 
+use circuits::AugmentedFCircuit;
+use folding::Folding;
+
+/// Configuration for ProtoGalaxy's CycleFold circuit
+pub struct ProtoGalaxyCycleFoldConfig<C: CurveGroup> {
+    _c: PhantomData<C>,
+}
+
+impl<C: CurveGroup> CycleFoldConfig for ProtoGalaxyCycleFoldConfig<C> {
+    const RANDOMNESS_BIT_LENGTH: usize = C::ScalarField::MODULUS_BIT_SIZE as usize;
+    const N_INPUT_POINTS: usize = 2;
+    type C = C;
+    type F = C::BaseField;
+}
+
+/// CycleFold circuit for computing random linear combinations of group elements
+/// in ProtoGalaxy instances.
+pub type ProtoGalaxyCycleFoldCircuit<C, GC> = CycleFoldCircuit<ProtoGalaxyCycleFoldConfig<C>, GC>;
+
 #[derive(Clone, Debug, PartialEq, Eq)]
 pub struct CommittedInstance<C: CurveGroup> {
     phi: C,
     betas: Vec<C::ScalarField>,
     e: C::ScalarField,
-    u: C::ScalarField,
     x: Vec<C::ScalarField>,
 }
 
+impl<C: CurveGroup> CommittedInstance<C> {
+    pub fn dummy_running(io_len: usize, t: usize) -> Self {
+        Self {
+            phi: C::zero(),
+            betas: vec![C::ScalarField::zero(); t],
+            e: C::ScalarField::zero(),
+            x: vec![C::ScalarField::zero(); io_len],
+        }
+    }
+
+    pub fn dummy_incoming(io_len: usize) -> Self {
+        Self::dummy_running(io_len, 0)
+    }
+}
+
+impl<C: CurveGroup> CommittedInstance<C>
+where
+    C::ScalarField: Absorb,
+    C::BaseField: PrimeField,
+{
+    /// hash implements the committed instance hash compatible with the gadget implemented in
+    /// CommittedInstanceVar.hash.
+    /// Returns `H(i, z_0, z_i, U_i)`, where `i` can be `i` but also `i+1`, and `U_i` is the
+    /// `CommittedInstance`.
+    pub fn hash(
+        &self,
+        sponge: &PoseidonSponge<C::ScalarField>,
+        pp_hash: C::ScalarField,
+        i: C::ScalarField,
+        z_0: Vec<C::ScalarField>,
+        z_i: Vec<C::ScalarField>,
+    ) -> C::ScalarField {
+        let mut sponge = sponge.clone();
+        sponge.absorb(&pp_hash);
+        sponge.absorb(&i);
+        sponge.absorb(&z_0);
+        sponge.absorb(&z_i);
+        sponge.absorb(&self);
+        sponge.squeeze_field_elements(1)[0]
+    }
+}
+
 #[derive(Clone, Debug)]
 pub struct CommittedInstanceVar<C: CurveGroup> {
     phi: NonNativeAffineVar<C>,
     betas: Vec<FpVar<C::ScalarField>>,
     e: FpVar<C::ScalarField>,
-    u: FpVar<C::ScalarField>,
     x: Vec<FpVar<C::ScalarField>>,
 }
 
@@ -50,20 +133,98 @@ Ok(Self {
                 phi: NonNativeAffineVar::new_variable(cs.clone(), || Ok(u.phi), mode)?,
                 betas: Vec::new_variable(cs.clone(), || Ok(u.betas.clone()), mode)?,
                 e: FpVar::new_variable(cs.clone(), || Ok(u.e), mode)?,
-                u: FpVar::new_variable(cs.clone(), || Ok(u.u), mode)?,
                 x: Vec::new_variable(cs.clone(), || Ok(u.x.clone()), mode)?,
             })
         })
     }
 }
 
+impl<C: CurveGroup> R1CSVar<C::ScalarField> for CommittedInstanceVar<C> {
+    type Value = CommittedInstance<C>;
+
+    fn cs(&self) -> ConstraintSystemRef<C::ScalarField> {
+        self.phi
+            .cs()
+            .or(self.betas.cs())
+            .or(self.e.cs())
+            .or(self.x.cs())
+    }
+
+    fn value(&self) -> Result<Self::Value, SynthesisError> {
+        Ok(CommittedInstance {
+            phi: self.phi.value()?,
+            betas: self
+                .betas
+                .iter()
+                .map(|v| v.value())
+                .collect::<Result<_, _>>()?,
+            e: self.e.value()?,
+            x: self.x.iter().map(|v| v.value()).collect::<Result<_, _>>()?,
+        })
+    }
+}
+
+impl<C: CurveGroup> CommittedInstanceVar<C>
+where
+    C::ScalarField: Absorb,
+    C::BaseField: PrimeField,
+{
+    /// hash implements the committed instance hash compatible with the native implementation from
+    /// CommittedInstance.hash.
+    /// Returns `H(i, z_0, z_i, U_i)`, where `i` can be `i` but also `i+1`, and `U` is the
+    /// `CommittedInstance`.
+    /// Additionally it returns the vector of the field elements from the self parameters, so they
+    /// can be reused in other gadgets avoiding recalculating (reconstraining) them.
+    #[allow(clippy::type_complexity)]
+    pub fn hash(
+        self,
+        sponge: &PoseidonSpongeVar<CF1<C>>,
+        pp_hash: FpVar<CF1<C>>,
+        i: FpVar<CF1<C>>,
+        z_0: Vec<FpVar<CF1<C>>>,
+        z_i: Vec<FpVar<CF1<C>>>,
+    ) -> Result<(FpVar<CF1<C>>, Vec<FpVar<CF1<C>>>), SynthesisError> {
+        let mut sponge = sponge.clone();
+        let U_vec = self.to_sponge_field_elements()?;
+        sponge.absorb(&pp_hash)?;
+        sponge.absorb(&i)?;
+        sponge.absorb(&z_0)?;
+        sponge.absorb(&z_i)?;
+        sponge.absorb(&U_vec)?;
+        Ok((sponge.squeeze_field_elements(1)?.pop().unwrap(), U_vec))
+    }
+}
+
 #[derive(Clone, Debug)]
 pub struct Witness<F: PrimeField> {
     w: Vec<F>,
     r_w: F,
 }
 
-#[derive(Debug, Error, PartialEq)]
+impl<F: PrimeField> Witness<F> {
+    pub fn new(w: Vec<F>) -> Self {
+        // note: at the current version, we don't use the blinding factors and we set them to 0
+        // always.
+        // Tracking issue: https://github.com/privacy-scaling-explorations/sonobe/issues/82
+        Self { w, r_w: F::zero() }
+    }
+
+    pub fn commit<CS: CommitmentScheme<C>, C: CurveGroup<ScalarField = F>>(
+        &self,
+        params: &CS::ProverParams,
+        x: Vec<F>,
+    ) -> Result<CommittedInstance<C>, crate::Error> {
+        let phi = CS::commit(params, &self.w, &self.r_w)?;
+        Ok(CommittedInstance {
+            phi,
+            x,
+            e: F::zero(),
+            betas: vec![],
+        })
+    }
+}
+
+#[derive(Debug, thiserror::Error, PartialEq)]
 pub enum ProtoGalaxyError {
     #[error("The remainder from G(X)-F(α)*L_0(X)) / Z(X) should be zero")]
     RemainderNotZero,
@@ -76,3 +237,821 @@ BTreeNotFull(usize),
     #[error("The lengths of β and δ do not equal: |β| = {0}, |δ|={0}")]
     WrongLenBetas(usize, usize),
 }
+
+/// Proving parameters for ProtoGalaxy-based IVC
+#[derive(Debug, Clone)]
+pub struct ProverParams<C1, C2, CS1, CS2>
+where
+    C1: CurveGroup,
+    C2: CurveGroup,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+{
+    /// Poseidon sponge configuration
+    pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// Proving parameters of the underlying commitment scheme over C1
+    pub cs_params: CS1::ProverParams,
+    /// Proving parameters of the underlying commitment scheme over C2
+    pub cf_cs_params: CS2::ProverParams,
+}
+
+/// Verification parameters for ProtoGalaxy-based IVC
+#[derive(Debug, Clone)]
+pub struct VerifierParams<C1, C2, CS1, CS2>
+where
+    C1: CurveGroup,
+    C2: CurveGroup,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+{
+    /// Poseidon sponge configuration
+    pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// R1CS of the Augmented step circuit
+    pub r1cs: R1CS<C1::ScalarField>,
+    /// R1CS of the CycleFold circuit
+    pub cf_r1cs: R1CS<C2::ScalarField>,
+    /// Verification parameters of the underlying commitment scheme over C1
+    pub cs_vp: CS1::VerifierParams,
+    /// Verification parameters of the underlying commitment scheme over C2
+    pub cf_cs_vp: CS2::VerifierParams,
+}
+
+impl<C1, C2, CS1, CS2> VerifierParams<C1, C2, CS1, CS2>
+where
+    C1: CurveGroup,
+    C2: CurveGroup,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+{
+    /// returns the hash of the public parameters of ProtoGalaxy
+    pub fn pp_hash(&self) -> Result<C1::ScalarField, Error> {
+        // TODO (@winderica): support hiding commitments in ProtoGalaxy.
+        // For now, `H` is set to false.
+        // Tracking issue: https://github.com/privacy-scaling-explorations/sonobe/issues/82
+        pp_hash::<C1, C2, CS1, CS2, false>(
+            &self.r1cs,
+            &self.cf_r1cs,
+            &self.cs_vp,
+            &self.cf_cs_vp,
+            &self.poseidon_config,
+        )
+    }
+}
+
+/// Implements ProtoGalaxy+CycleFold's IVC, described in [ProtoGalaxy] and
+/// [CycleFold], following the FoldingScheme trait
+///
+/// [ProtoGalaxy]: https://eprint.iacr.org/2023/1106.pdf
+/// [CycleFold]: https://eprint.iacr.org/2023/1192.pdf
+#[derive(Clone, Debug)]
+pub struct ProtoGalaxy<C1, GC1, C2, GC2, FC, CS1, CS2>
+where
+    C1: CurveGroup,
+    GC1: CurveVar<C1, CF2<C1>> + ToConstraintFieldGadget<CF2<C1>>,
+    C2: CurveGroup,
+    GC2: CurveVar<C2, CF2<C2>>,
+    FC: FCircuit<C1::ScalarField>,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+{
+    _gc1: PhantomData<GC1>,
+    _c2: PhantomData<C2>,
+    _gc2: PhantomData<GC2>,
+    /// R1CS of the Augmented Function circuit
+    pub r1cs: R1CS<C1::ScalarField>,
+    /// R1CS of the CycleFold circuit
+    pub cf_r1cs: R1CS<C2::ScalarField>,
+    pub poseidon_config: PoseidonConfig<C1::ScalarField>,
+    /// CommitmentScheme::ProverParams over C1
+    pub cs_params: CS1::ProverParams,
+    /// CycleFold CommitmentScheme::ProverParams, over C2
+    pub cf_cs_params: CS2::ProverParams,
+    /// F circuit, the circuit that is being folded
+    pub F: FC,
+    /// public params hash
+    pub pp_hash: C1::ScalarField,
+    pub i: C1::ScalarField,
+    /// initial state
+    pub z_0: Vec<C1::ScalarField>,
+    /// current i-th state
+    pub z_i: Vec<C1::ScalarField>,
+    /// ProtoGalaxy instances
+    pub w_i: Witness<C1::ScalarField>,
+    pub u_i: CommittedInstance<C1>,
+    pub W_i: Witness<C1::ScalarField>,
+    pub U_i: CommittedInstance<C1>,
+
+    /// CycleFold running instance
+    pub cf_W_i: CycleFoldWitness<C2>,
+    pub cf_U_i: CycleFoldCommittedInstance<C2>,
+}
+
+impl<C1, GC1, C2, GC2, FC, CS1, CS2> ProtoGalaxy<C1, GC1, C2, GC2, FC, CS1, CS2>
+where
+    C1: CurveGroup,
+    GC1: CurveVar<C1, CF2<C1>> + ToConstraintFieldGadget<CF2<C1>>,
+    C2: CurveGroup,
+    GC2: CurveVar<C2, CF2<C2>> + ToConstraintFieldGadget<CF2<C2>>,
+    FC: FCircuit<C1::ScalarField>,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+    <C1 as CurveGroup>::BaseField: PrimeField,
+    <C2 as CurveGroup>::BaseField: PrimeField,
+    C1::ScalarField: Absorb,
+    C2::ScalarField: Absorb,
+    C1: CurveGroup<BaseField = C2::ScalarField, ScalarField = C2::BaseField>,
+    for<'a> &'a GC1: GroupOpsBounds<'a, C1, GC1>,
+    for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
+{
+    /// This method computes the parameter `t` in ProtoGalaxy for folding `F'`,
+    /// the augmented circuit of `F`
+    fn compute_t(
+        poseidon_config: &PoseidonConfig<CF1<C1>>,
+        F: &FC,
+        d: usize,
+        k: usize,
+    ) -> Result<usize, SynthesisError> {
+        // In ProtoGalaxy, prover and verifier are parameterized by `t = log(n)`
+        // where `n` is the number of constraints in the circuit (known as the
+        // mapping `f` in the paper).
+        // For IVC, `f` is the augmented circuit `F'`, which not only includes
+        // the original computation `F`, but also the in-circuit verifier of
+        // ProtoGalaxy.
+        // Therefore, `t` depends on the size of `F'`, but the size of `F'` in
+        // turn depends on `t`.
+        // To address this circular dependency, we first find `t_lower_bound`,
+        // the lower bound of `t`. Then we incrementally increase `t` and build
+        // the circuit `F'` with `t` as ProtoGalaxy's parameter, until `t` is
+        // the smallest integer that equals the logarithm of the number of
+        // constraints.
+
+        // For `t_lower_bound`, we configure `F'` with `t = 1` and compute log2
+        // of the size of `F'`.
+        let state_len = F.state_len();
+        let external_inputs_len = F.external_inputs_len();
+
+        // `F'` includes `F` and `ProtoGalaxy.V`, where `F` might be costly.
+        // Observing that the cost of `F` is constant with respect to `t`, we
+        // separately compute `step_constraints`, the size of `F`.
+        // Later, we only need to re-run the rest of `F'` with updated `t` to
+        // get the size of `F'`.
+        let cs = ConstraintSystem::<C1::ScalarField>::new_ref();
+        F.generate_step_constraints(
+            cs.clone(),
+            0,
+            Vec::new_witness(cs.clone(), || Ok(vec![Zero::zero(); state_len]))?,
+            Vec::new_witness(cs.clone(), || Ok(vec![Zero::zero(); external_inputs_len]))?,
+        )?;
+        let step_constraints = cs.num_constraints();
+
+        // Create a dummy circuit with the same state length and external inputs
+        // length as `F`, which replaces `F` in the augmented circuit `F'`.
+        let dummy_circuit: DummyCircuit =
+            FCircuit::<C1::ScalarField>::new((state_len, external_inputs_len)).unwrap();
+
+        // Compute `augmentation_constraints`, the size of `F'` without `F`.
+        let cs = ConstraintSystem::<C1::ScalarField>::new_ref();
+        AugmentedFCircuit::<C1, C2, GC2, DummyCircuit>::empty(
+            poseidon_config,
+            dummy_circuit.clone(),
+            1,
+            d,
+            k,
+        )
+        .generate_constraints(cs.clone())?;
+        let augmentation_constraints = cs.num_constraints();
+
+        // The sum of `step_constraints` and `augmentation_constraints` is the
+        // size of `F'` with `t = 1`, and hence the actual `t` should have lower
+        // bound `log2(step_constraints + augmentation_constraints)`.
+        let t_lower_bound = log2(step_constraints + augmentation_constraints) as usize;
+        // Optimization: we in fact only need to try two values of `t`.
+        // This is because increasing `t` will only slightly affect the size of
+        // `F'` (more specifically, the size of `F'` will never be doubled).
+        // Thus, `t_lower_bound` (the log2 size of `F'` with `t = 1`) is very
+        // close to the actual `t` (either `t` or `t - 1`).
+        let t_upper_bound = t_lower_bound + 1;
+
+        for t in t_lower_bound..=t_upper_bound {
+            let cs = ConstraintSystem::<C1::ScalarField>::new_ref();
+            AugmentedFCircuit::<C1, C2, GC2, DummyCircuit>::empty(
+                poseidon_config,
+                dummy_circuit.clone(),
+                t,
+                d,
+                k,
+            )
+            .generate_constraints(cs.clone())?;
+            let augmentation_constraints = cs.num_constraints();
+            if t == log2(step_constraints + augmentation_constraints) as usize {
+                return Ok(t);
+            }
+        }
+        unreachable!()
+    }
+}
+
+impl<C1, GC1, C2, GC2, FC, CS1, CS2> FoldingScheme<C1, C2, FC>
+    for ProtoGalaxy<C1, GC1, C2, GC2, FC, CS1, CS2>
+where
+    C1: CurveGroup,
+    GC1: CurveVar<C1, CF2<C1>> + ToConstraintFieldGadget<CF2<C1>>,
+    C2: CurveGroup,
+    GC2: CurveVar<C2, CF2<C2>> + ToConstraintFieldGadget<CF2<C2>>,
+    FC: FCircuit<C1::ScalarField>,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+    <C1 as CurveGroup>::BaseField: PrimeField,
+    <C2 as CurveGroup>::BaseField: PrimeField,
+    C1::ScalarField: Absorb,
+    C2::ScalarField: Absorb,
+    C1: CurveGroup<BaseField = C2::ScalarField, ScalarField = C2::BaseField>,
+    for<'a> &'a GC1: GroupOpsBounds<'a, C1, GC1>,
+    for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
+{
+    type PreprocessorParam = (PoseidonConfig<CF1<C1>>, FC);
+    type ProverParam = ProverParams<C1, C2, CS1, CS2>;
+    type VerifierParam = VerifierParams<C1, C2, CS1, CS2>;
+    type RunningInstance = (CommittedInstance<C1>, Witness<C1::ScalarField>);
+    type IncomingInstance = (CommittedInstance<C1>, Witness<C1::ScalarField>);
+    type MultiCommittedInstanceWithWitness = (CommittedInstance<C1>, Witness<C1::ScalarField>);
+    type CFInstance = (CycleFoldCommittedInstance<C2>, CycleFoldWitness<C2>);
+
+    fn preprocess(
+        mut rng: impl RngCore,
+        (poseidon_config, F): &Self::PreprocessorParam,
+    ) -> Result<(Self::ProverParam, Self::VerifierParam), Error> {
+        // We fix `k`, the number of incoming instances, to 1, because
+        // multi-instances folding is not supported yet.
+        // TODO (@winderica): Support multi-instances folding and make `k` a
+        // constant generic parameter (as in HyperNova)
+        // Tracking issue: https://github.com/privacy-scaling-explorations/sonobe/issues/82
+        let k = 1;
+        // `d`, the degree of the constraint system, is set to 2, as we only
+        // support R1CS for now, whose highest degree is 2.
+        let d = 2;
+        let t = Self::compute_t(poseidon_config, F, d, k)?;
+
+        // prepare the circuit to obtain its R1CS
+        let cs = ConstraintSystem::<C1::ScalarField>::new_ref();
+        let cs2 = ConstraintSystem::<C1::BaseField>::new_ref();
+
+        let augmented_F_circuit =
+            AugmentedFCircuit::<C1, C2, GC2, FC>::empty(poseidon_config, F.clone(), t, d, k);
+        let cf_circuit = ProtoGalaxyCycleFoldCircuit::<C1, GC1>::empty();
+
+        augmented_F_circuit.generate_constraints(cs.clone())?;
+        cs.finalize();
+        let cs = cs.into_inner().ok_or(Error::NoInnerConstraintSystem)?;
+        let r1cs = extract_r1cs::<C1::ScalarField>(&cs);
+
+        cf_circuit.generate_constraints(cs2.clone())?;
+        cs2.finalize();
+        let cs2 = cs2.into_inner().ok_or(Error::NoInnerConstraintSystem)?;
+        let cf_r1cs = extract_r1cs::<C1::BaseField>(&cs2);
+
+        let (cs_pp, cs_vp) = CS1::setup(&mut rng, r1cs.A.n_rows)?;
+        let (cf_cs_pp, cf_cs_vp) = CS2::setup(&mut rng, max(cf_r1cs.A.n_rows, cf_r1cs.A.n_cols))?;
+
+        Ok((
+            Self::ProverParam {
+                poseidon_config: poseidon_config.clone(),
+                cs_params: cs_pp,
+                cf_cs_params: cf_cs_pp,
+            },
+            Self::VerifierParam {
+                poseidon_config: poseidon_config.clone(),
+                r1cs,
+                cf_r1cs,
+                cs_vp,
+                cf_cs_vp,
+            },
+        ))
+    }
+
+    /// Initializes the ProtoGalaxy+CycleFold's IVC for the given parameters and
+    /// initial state `z_0`.
+    fn init(
+        (pp, vp): &(Self::ProverParam, Self::VerifierParam),
+        F: FC,
+        z_0: Vec<C1::ScalarField>,
+    ) -> Result<Self, Error> {
+        // compute the public params hash
+        let pp_hash = vp.pp_hash()?;
+
+        // setup the dummy instances
+        let (W_dummy, U_dummy) = vp.r1cs.dummy_running_instance();
+        let (w_dummy, u_dummy) = vp.r1cs.dummy_incoming_instance();
+        let (cf_W_dummy, cf_U_dummy) = vp.cf_r1cs.dummy_running_instance();
+
+        // W_dummy=W_0 is a 'dummy witness', all zeroes, but with the size corresponding to the
+        // R1CS that we're working with.
+        Ok(Self {
+            _gc1: PhantomData,
+            _c2: PhantomData,
+            _gc2: PhantomData,
+            r1cs: vp.r1cs.clone(),
+            cf_r1cs: vp.cf_r1cs.clone(),
+            poseidon_config: pp.poseidon_config.clone(),
+            cs_params: pp.cs_params.clone(),
+            cf_cs_params: pp.cf_cs_params.clone(),
+            F,
+            pp_hash,
+            i: C1::ScalarField::zero(),
+            z_0: z_0.clone(),
+            z_i: z_0,
+            w_i: w_dummy,
+            u_i: u_dummy,
+            W_i: W_dummy,
+            U_i: U_dummy,
+            // cyclefold running instance
+            cf_W_i: cf_W_dummy,
+            cf_U_i: cf_U_dummy,
+        })
+    }
+
+    /// Implements IVC.P of ProtoGalaxy+CycleFold
+    fn prove_step(
+        &mut self,
+        mut rng: impl RngCore,
+        external_inputs: Vec<C1::ScalarField>,
+        _other_instances: Option<Self::MultiCommittedInstanceWithWitness>,
+    ) -> Result<(), Error> {
+        // Multi-instances folding is not supported yet.
+        if _other_instances.is_some() {
+            return Err(Error::NoMultiInstances);
+        }
+        // We fix `k`, the number of incoming instances, to 1, because
+        // multi-instances folding is not supported yet.
+        // TODO (@winderica): Support multi-instances folding and make `k` a
+        // constant generic parameter (as in HyperNova)
+        // Tracking issue: https://github.com/privacy-scaling-explorations/sonobe/issues/82
+        let k = 1;
+        // `d`, the degree of the constraint system, is set to 2, as we only
+        // support R1CS for now, whose highest degree is 2.
+        let d = 2;
+
+        // `sponge` is for digest computation.
+        let sponge = PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config);
+        // `transcript` is for challenge generation.
+        let mut transcript_prover = sponge.clone();
+
+        let mut augmented_F_circuit: AugmentedFCircuit<C1, C2, GC2, FC>;
+
+        if self.z_i.len() != self.F.state_len() {
+            return Err(Error::NotSameLength(
+                "z_i.len()".to_string(),
+                self.z_i.len(),
+                "F.state_len()".to_string(),
+                self.F.state_len(),
+            ));
+        }
+        if external_inputs.len() != self.F.external_inputs_len() {
+            return Err(Error::NotSameLength(
+                "F.external_inputs_len()".to_string(),
+                self.F.external_inputs_len(),
+                "external_inputs.len()".to_string(),
+                external_inputs.len(),
+            ));
+        }
+
+        let i_bn: BigUint = self.i.into();
+        let i_usize: usize = i_bn.try_into().map_err(|_| Error::MaxStep)?;
+
+        let z_i1 = self
+            .F
+            .step_native(i_usize, self.z_i.clone(), external_inputs.clone())?;
+
+        // folded instance output (public input, x)
+        // u_{i+1}.x[0] = H(i+1, z_0, z_{i+1}, U_{i+1})
+        let u_i1_x: C1::ScalarField;
+        // u_{i+1}.x[1] = H(cf_U_{i+1})
+        let cf_u_i1_x: C1::ScalarField;
+
+        if self.i.is_zero() {
+            // Take extra care of the base case
+            // `U_{i+1}` (i.e., `U_1`) is fixed to `U_dummy`, so we just use
+            // `self.U_i = U_0 = U_dummy`.
+            u_i1_x = self.U_i.hash(
+                &sponge,
+                self.pp_hash,
+                self.i + C1::ScalarField::one(),
+                self.z_0.clone(),
+                z_i1.clone(),
+            );
+            // `cf_U_{i+1}` (i.e., `cf_U_1`) is fixed to `cf_U_dummy`, so we
+            // just use `self.cf_U_i = cf_U_0 = cf_U_dummy`.
+            cf_u_i1_x = self.cf_U_i.hash_cyclefold(&sponge, self.pp_hash);
+
+            augmented_F_circuit = AugmentedFCircuit::empty(
+                &self.poseidon_config,
+                self.F.clone(),
+                self.U_i.betas.len(),
+                d,
+                k,
+            );
+            augmented_F_circuit.pp_hash = self.pp_hash;
+            augmented_F_circuit.z_0.clone_from(&self.z_0);
+            augmented_F_circuit.z_i.clone_from(&self.z_i);
+            augmented_F_circuit
+                .external_inputs
+                .clone_from(&external_inputs);
+
+            // There is no need to update `self.U_i` etc. as they are unchanged.
+        } else {
+            // Primary part:
+            // Compute `U_{i+1}` by folding `u_i` into `U_i`.
+            let (U_i1, W_i1, F_coeffs, K_coeffs, L_evals, phi_stars) = Folding::prove(
+                &mut transcript_prover,
+                &self.r1cs,
+                &self.U_i,
+                &self.W_i,
+                &[self.u_i.clone()],
+                &[self.w_i.clone()],
+            )?;
+
+            // CycleFold part:
+            // get the vector used as public inputs 'x' in the CycleFold circuit
+            let mut r0_bits = L_evals[0].into_bigint().to_bits_le();
+            let mut r1_bits = L_evals[1].into_bigint().to_bits_le();
+            r0_bits.resize(C1::ScalarField::MODULUS_BIT_SIZE as usize, false);
+            r1_bits.resize(C1::ScalarField::MODULUS_BIT_SIZE as usize, false);
+
+            // cyclefold circuit for enforcing:
+            // 0 + U_i.phi * L_evals[0] == phi_stars[0]
+            let cf1_u_i_x = [
+                r0_bits
+                    .chunks(C1::BaseField::MODULUS_BIT_SIZE as usize - 1)
+                    .map(BigInteger::from_bits_le)
+                    .map(C1::BaseField::from_bigint)
+                    .collect::<Option<Vec<_>>>()
+                    .unwrap(),
+                get_cm_coordinates(&C1::zero()),
+                get_cm_coordinates(&self.U_i.phi),
+                get_cm_coordinates(&phi_stars[0]),
+            ]
+            .concat();
+            let cf1_circuit = ProtoGalaxyCycleFoldCircuit::<C1, GC1> {
+                _gc: PhantomData,
+                r_bits: Some(r0_bits),
+                points: Some(vec![C1::zero(), self.U_i.phi]),
+                x: Some(cf1_u_i_x.clone()),
+            };
+
+            // cyclefold circuit for enforcing:
+            // phi_stars[0] + u_i.phi * L_evals[1] == U_i1.phi
+            // i.e., U_i.phi * L_evals[0] + u_i.phi * L_evals[1] == U_i1.phi
+            let cf2_u_i_x = [
+                r1_bits
+                    .chunks(C1::BaseField::MODULUS_BIT_SIZE as usize - 1)
+                    .map(BigInteger::from_bits_le)
+                    .map(C1::BaseField::from_bigint)
+                    .collect::<Option<Vec<_>>>()
+                    .unwrap(),
+                get_cm_coordinates(&phi_stars[0]),
+                get_cm_coordinates(&self.u_i.phi),
+                get_cm_coordinates(&U_i1.phi),
+            ]
+            .concat();
+            let cf2_circuit = ProtoGalaxyCycleFoldCircuit::<C1, GC1> {
+                _gc: PhantomData,
+                r_bits: Some(r1_bits),
+                points: Some(vec![phi_stars[0], self.u_i.phi]),
+                x: Some(cf2_u_i_x.clone()),
+            };
+
+            // fold self.cf_U_i + cf1_U -> folded running with cf1
+            let (_cf1_w_i, cf1_u_i, cf1_W_i1, cf1_U_i1, cf1_cmT, _) = self.fold_cyclefold_circuit(
+                &mut transcript_prover,
+                self.cf_W_i.clone(), // CycleFold running instance witness
+                self.cf_U_i.clone(), // CycleFold running instance
+                cf1_u_i_x,
+                cf1_circuit,
+                &mut rng,
+            )?;
+            // fold [the output from folding self.cf_U_i + cf1_U] + cf2_U = folded_running_with_cf1 + cf2
+            let (_cf2_w_i, cf2_u_i, cf_W_i1, cf_U_i1, cf2_cmT, _) = self.fold_cyclefold_circuit(
+                &mut transcript_prover,
+                cf1_W_i1,
+                cf1_U_i1.clone(),
+                cf2_u_i_x,
+                cf2_circuit,
+                &mut rng,
+            )?;
+
+            // Derive `u_{i+1}.x[0], u_{i+1}.x[1]` by hashing folded instances
+            u_i1_x = U_i1.hash(
+                &sponge,
+                self.pp_hash,
+                self.i + C1::ScalarField::one(),
+                self.z_0.clone(),
+                z_i1.clone(),
+            );
+            cf_u_i1_x = cf_U_i1.hash_cyclefold(&sponge, self.pp_hash);
+
+            augmented_F_circuit = AugmentedFCircuit {
+                _gc2: PhantomData,
+                poseidon_config: self.poseidon_config.clone(),
+                pp_hash: self.pp_hash,
+                i: self.i,
+                i_usize,
+                z_0: self.z_0.clone(),
+                z_i: self.z_i.clone(),
+                external_inputs: external_inputs.clone(),
+                u_i_phi: self.u_i.phi,
+                U_i: self.U_i.clone(),
+                U_i1_phi: U_i1.phi,
+                F_coeffs: F_coeffs.clone(),
+                K_coeffs: K_coeffs.clone(),
+                phi_stars,
+                F: self.F.clone(),
+                x: Some(u_i1_x),
+                // cyclefold values
+                cf1_u_i_cmW: cf1_u_i.cmW,
+                cf2_u_i_cmW: cf2_u_i.cmW,
+                cf_U_i: self.cf_U_i.clone(),
+                cf1_cmT,
+                cf2_cmT,
+                cf_x: Some(cf_u_i1_x),
+            };
+
+            #[cfg(test)]
+            {
+                let mut transcript_verifier = sponge.clone();
+                assert_eq!(
+                    Folding::verify(
+                        &mut transcript_verifier,
+                        &self.U_i,
+                        &[self.u_i.clone()],
+                        F_coeffs,
+                        K_coeffs
+                    )?,
+                    U_i1
+                );
+                self.cf_r1cs.check_tight_relation(&_cf1_w_i, &cf1_u_i)?;
+                self.cf_r1cs.check_tight_relation(&_cf2_w_i, &cf2_u_i)?;
+                self.cf_r1cs
+                    .check_relaxed_relation(&self.cf_W_i, &self.cf_U_i)?;
+            }
+
+            self.W_i = W_i1;
+            self.U_i = U_i1;
+            self.cf_W_i = cf_W_i1;
+            self.cf_U_i = cf_U_i1;
+        }
+
+        let cs = ConstraintSystem::<C1::ScalarField>::new_ref();
+
+        augmented_F_circuit.generate_constraints(cs.clone())?;
+
+        #[cfg(test)]
+        assert!(cs.is_satisfied().unwrap());
+
+        let cs = cs.into_inner().ok_or(Error::NoInnerConstraintSystem)?;
+        let (w_i1, x_i1) = extract_w_x::<C1::ScalarField>(&cs);
+        if x_i1[0] != u_i1_x || x_i1[1] != cf_u_i1_x {
+            return Err(Error::NotEqual);
+        }
+
+        #[cfg(test)]
+        if x_i1.len() != 2 {
+            return Err(Error::NotExpectedLength(x_i1.len(), 2));
+        }
+
+        // set values for next iteration
+        self.i += C1::ScalarField::one();
+        self.z_i = z_i1;
+        self.w_i = Witness::new(w_i1);
+        self.u_i = self.w_i.commit::<CS1, C1>(&self.cs_params, x_i1)?;
+
+        #[cfg(test)]
+        {
+            self.r1cs.check_tight_relation(&self.w_i, &self.u_i)?;
+            self.r1cs.check_relaxed_relation(&self.W_i, &self.U_i)?;
+        }
+
+        Ok(())
+    }
+
+    fn state(&self) -> Vec<C1::ScalarField> {
+        self.z_i.clone()
+    }
+    fn instances(
+        &self,
+    ) -> (
+        Self::RunningInstance,
+        Self::IncomingInstance,
+        Self::CFInstance,
+    ) {
+        (
+            (self.U_i.clone(), self.W_i.clone()),
+            (self.u_i.clone(), self.w_i.clone()),
+            (self.cf_U_i.clone(), self.cf_W_i.clone()),
+        )
+    }
+
+    /// Implements IVC.V of ProtoGalaxy+CycleFold
+    fn verify(
+        vp: Self::VerifierParam,
+        z_0: Vec<C1::ScalarField>, // initial state
+        z_i: Vec<C1::ScalarField>, // last state
+        num_steps: C1::ScalarField,
+        running_instance: Self::RunningInstance,
+        incoming_instance: Self::IncomingInstance,
+        cyclefold_instance: Self::CFInstance,
+    ) -> Result<(), Error> {
+        let sponge = PoseidonSponge::<C1::ScalarField>::new(&vp.poseidon_config);
+
+        let (U_i, W_i) = running_instance;
+        let (u_i, w_i) = incoming_instance;
+        let (cf_U_i, cf_W_i) = cyclefold_instance;
+
+        if u_i.x.len() != 2 || U_i.x.len() != 2 {
+            return Err(Error::IVCVerificationFail);
+        }
+
+        let pp_hash = vp.pp_hash()?;
+
+        // check that u_i's output points to the running instance
+        // u_i.X[0] == H(i, z_0, z_i, U_i)
+        let expected_u_i_x = U_i.hash(&sponge, pp_hash, num_steps, z_0, z_i.clone());
+        if expected_u_i_x != u_i.x[0] {
+            return Err(Error::IVCVerificationFail);
+        }
+        // u_i.X[1] == H(cf_U_i)
+        let expected_cf_u_i_x = cf_U_i.hash_cyclefold(&sponge, pp_hash);
+        if expected_cf_u_i_x != u_i.x[1] {
+            return Err(Error::IVCVerificationFail);
+        }
+
+        // check R1CS satisfiability
+        vp.r1cs.check_tight_relation(&w_i, &u_i)?;
+        // check RelaxedR1CS satisfiability
+        vp.r1cs.check_relaxed_relation(&W_i, &U_i)?;
+
+        // check CycleFold RelaxedR1CS satisfiability
+        vp.cf_r1cs.check_relaxed_relation(&cf_W_i, &cf_U_i)?;
+
+        Ok(())
+    }
+}
+
+impl<C1, GC1, C2, GC2, FC, CS1, CS2> ProtoGalaxy<C1, GC1, C2, GC2, FC, CS1, CS2>
+where
+    C1: CurveGroup,
+    GC1: CurveVar<C1, CF2<C1>> + ToConstraintFieldGadget<CF2<C1>>,
+    C2: CurveGroup,
+    GC2: CurveVar<C2, CF2<C2>> + ToConstraintFieldGadget<CF2<C2>>,
+    FC: FCircuit<C1::ScalarField>,
+    CS1: CommitmentScheme<C1>,
+    CS2: CommitmentScheme<C2>,
+    <C1 as CurveGroup>::BaseField: PrimeField,
+    <C2 as CurveGroup>::BaseField: PrimeField,
+    <C1 as Group>::ScalarField: Absorb,
+    <C2 as Group>::ScalarField: Absorb,
+    C1: CurveGroup<BaseField = C2::ScalarField, ScalarField = C2::BaseField>,
+    for<'a> &'a GC1: GroupOpsBounds<'a, C1, GC1>,
+    for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
+{
+    // folds the given cyclefold circuit and its instances
+    #[allow(clippy::type_complexity)]
+    fn fold_cyclefold_circuit(
+        &self,
+        transcript: &mut PoseidonSponge<C1::ScalarField>,
+        cf_W_i: CycleFoldWitness<C2>, // witness of the running instance
+        cf_U_i: CycleFoldCommittedInstance<C2>, // running instance
+        cf_u_i_x: Vec<C2::ScalarField>,
+        cf_circuit: ProtoGalaxyCycleFoldCircuit<C1, GC1>,
+        rng: &mut impl RngCore,
+    ) -> Result<
+        (
+            CycleFoldWitness<C2>,
+            CycleFoldCommittedInstance<C2>, // u_i
+            CycleFoldWitness<C2>,           // W_i1
+            CycleFoldCommittedInstance<C2>, // U_i1
+            C2,                             // cmT
+            C2::ScalarField,                // r_Fq
+        ),
+        Error,
+    > {
+        fold_cyclefold_circuit::<ProtoGalaxyCycleFoldConfig<C1>, C1, GC1, C2, GC2, CS2, false>(
+            transcript,
+            self.cf_r1cs.clone(),
+            self.cf_cs_params.clone(),
+            self.pp_hash,
+            cf_W_i,
+            cf_U_i,
+            cf_u_i_x,
+            cf_circuit,
+            rng,
+        )
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use ark_bn254::{constraints::GVar, Bn254, Fr, G1Projective as Projective};
+    use ark_grumpkin::{constraints::GVar as GVar2, Projective as Projective2};
+    use ark_std::test_rng;
+    use rayon::prelude::*;
+
+    use crate::{
+        commitment::{kzg::KZG, pedersen::Pedersen},
+        frontend::utils::CubicFCircuit,
+        transcript::poseidon::poseidon_canonical_config,
+    };
+
+    /// This test tests the ProtoGalaxy+CycleFold IVC, and by consequence it is
+    /// also testing the AugmentedFCircuit
+    #[test]
+    fn test_ivc() {
+        let poseidon_config = poseidon_canonical_config::<Fr>();
+
+        let F_circuit = CubicFCircuit::<Fr>::new(()).unwrap();
+
+        // run the test using Pedersen commitments on both sides of the curve cycle
+        test_ivc_opt::<Pedersen<Projective>, Pedersen<Projective2>>(
+            poseidon_config.clone(),
+            F_circuit,
+        );
+        // run the test using KZG for the commitments on the main curve, and Pedersen for the
+        // commitments on the secondary curve
+        test_ivc_opt::<KZG<Bn254>, Pedersen<Projective2>>(poseidon_config, F_circuit);
+    }
+
+    // test_ivc allowing to choose the CommitmentSchemes
+    fn test_ivc_opt<CS1: CommitmentScheme<Projective>, CS2: CommitmentScheme<Projective2>>(
+        poseidon_config: PoseidonConfig<Fr>,
+        F_circuit: CubicFCircuit<Fr>,
+    ) {
+        type PG<CS1, CS2> =
+            ProtoGalaxy<Projective, GVar, Projective2, GVar2, CubicFCircuit<Fr>, CS1, CS2>;
+
+        let params =
+            PG::<CS1, CS2>::preprocess(&mut test_rng(), &(poseidon_config, F_circuit)).unwrap();
+
+        let z_0 = vec![Fr::from(3_u32)];
+        let mut protogalaxy = PG::init(&params, F_circuit, z_0.clone()).unwrap();
+
+        let num_steps: usize = 3;
+        for _ in 0..num_steps {
+            protogalaxy
+                .prove_step(&mut test_rng(), vec![], None)
+                .unwrap();
+        }
+        assert_eq!(Fr::from(num_steps as u32), protogalaxy.i);
+
+        let (running_instance, incoming_instance, cyclefold_instance) = protogalaxy.instances();
+        PG::<CS1, CS2>::verify(
+            params.1,
+            z_0,
+            protogalaxy.z_i,
+            protogalaxy.i,
+            running_instance,
+            incoming_instance,
+            cyclefold_instance,
+        )
+        .unwrap();
+    }
+
+    #[ignore]
+    #[test]
+    fn test_t_bounds() {
+        let d = 2;
+        let k = 1;
+
+        let poseidon_config = poseidon_canonical_config::<Fr>();
+        for state_len in [1, 10, 100] {
+            for external_inputs_len in [1, 10, 100] {
+                let dummy_circuit: DummyCircuit =
+                    FCircuit::<Fr>::new((state_len, external_inputs_len)).unwrap();
+
+                let costs = (1..32)
+                    .into_par_iter()
+                    .map(|t| {
+                        let cs = ConstraintSystem::<Fr>::new_ref();
+                        AugmentedFCircuit::<Projective, Projective2, GVar2, DummyCircuit>::empty(
+                            &poseidon_config,
+                            dummy_circuit.clone(),
+                            t,
+                            d,
+                            k,
+                        )
+                        .generate_constraints(cs.clone())
+                        .unwrap();
+                        cs.num_constraints()
+                    })
+                    .collect::<Vec<_>>();
+
+                for t_lower_bound in log2(costs[0]) as usize..32 {
+                    let num_constraints =
+                        (1 << t_lower_bound) - costs[0] + costs[t_lower_bound - 1];
+                    let t = log2(num_constraints) as usize;
+                    assert!(t == t_lower_bound || t == t_lower_bound + 1);
+                }
+            }
+        }
+    }
+}

@@ -3,9 +3,15 @@ use ark_ec::CurveGroup;
 use ark_ff::PrimeField;
 use ark_r1cs_std::{fields::fp::FpVar, uint8::UInt8, ToConstraintFieldGadget};
 use ark_relations::r1cs::SynthesisError;
+use ark_std::{cfg_iter, log2, rand::RngCore, One, Zero};
+use rayon::prelude::*;
 
-use super::{CommittedInstance, CommittedInstanceVar};
-use crate::transcript::AbsorbNonNative;
+use super::{utils::pow_i, CommittedInstance, CommittedInstanceVar, Witness};
+use crate::{
+    arith::r1cs::{RelaxedR1CS, R1CS},
+    transcript::AbsorbNonNative,
+    Error,
+};
 
 // Implements the trait for absorbing ProtoGalaxy's CommittedInstance.
 impl<C: CurveGroup> Absorb for CommittedInstance<C>
@@ -22,7 +28,6 @@ .to_native_sponge_field_elements_as_vec()
             .to_sponge_field_elements(dest);
         self.betas.to_sponge_field_elements(dest);
         self.e.to_sponge_field_elements(dest);
-        self.u.to_sponge_field_elements(dest);
         self.x.to_sponge_field_elements(dest);
     }
 }
@@ -38,9 +43,74 @@ Ok([
             self.phi.to_constraint_field()?,
             self.betas.to_sponge_field_elements()?,
             self.e.to_sponge_field_elements()?,
-            self.u.to_sponge_field_elements()?,
             self.x.to_sponge_field_elements()?,
         ]
         .concat())
     }
 }
+
+impl<C: CurveGroup> RelaxedR1CS<C, Witness<C::ScalarField>, CommittedInstance<C>>
+    for R1CS<C::ScalarField>
+{
+    fn dummy_running_instance(&self) -> (Witness<C::ScalarField>, CommittedInstance<C>) {
+        let w_len = self.A.n_cols - 1 - self.l;
+        let w_dummy = Witness::new(vec![C::ScalarField::zero(); w_len]);
+        let u_dummy = CommittedInstance::<C>::dummy_running(self.l, log2(self.A.n_rows) as usize);
+        (w_dummy, u_dummy)
+    }
+
+    fn dummy_incoming_instance(&self) -> (Witness<C::ScalarField>, CommittedInstance<C>) {
+        let w_len = self.A.n_cols - 1 - self.l;
+        let w_dummy = Witness::new(vec![C::ScalarField::zero(); w_len]);
+        let u_dummy = CommittedInstance::<C>::dummy_incoming(self.l);
+        (w_dummy, u_dummy)
+    }
+
+    fn is_relaxed(_w: &Witness<C::ScalarField>, u: &CommittedInstance<C>) -> bool {
+        u.e != C::ScalarField::zero() || !u.betas.is_empty()
+    }
+
+    fn extract_z(w: &Witness<C::ScalarField>, u: &CommittedInstance<C>) -> Vec<C::ScalarField> {
+        [&[C::ScalarField::one()][..], &u.x, &w.w].concat()
+    }
+
+    fn check_error_terms(
+        _w: &Witness<C::ScalarField>,
+        u: &CommittedInstance<C>,
+        e: Vec<C::ScalarField>,
+    ) -> Result<(), Error> {
+        if u.betas.len() != log2(e.len()) as usize {
+            return Err(Error::NotSameLength(
+                "instance.betas.len()".to_string(),
+                u.betas.len(),
+                "log2(e.len())".to_string(),
+                log2(e.len()) as usize,
+            ));
+        }
+
+        let r = cfg_iter!(e)
+            .enumerate()
+            .map(|(i, e_i)| pow_i(i, &u.betas) * e_i)
+            .sum();
+        if u.e == r {
+            Ok(())
+        } else {
+            Err(Error::NotSatisfied)
+        }
+    }
+
+    fn sample<CS>(
+        &self,
+        _params: &CS::ProverParams,
+        _rng: impl RngCore,
+    ) -> Result<(Witness<C::ScalarField>, CommittedInstance<C>), Error>
+    where
+        CS: crate::commitment::CommitmentScheme<C, true>,
+    {
+        // Sampling a random pair of witness and committed instance is required
+        // for the zero-knowledge layer for ProtoGalaxy, which is not supported
+        // yet.
+        // Tracking issue: https://github.com/privacy-scaling-explorations/sonobe/issues/82
+        unimplemented!()
+    }
+}

@@ -1,5 +1,6 @@
 use ark_ff::PrimeField;
 use ark_r1cs_std::fields::{fp::FpVar, FieldVar};
+use num_integer::Integer;
 
 /// Returns (b, b^2, b^4, ..., b^{2^{t-1}})
 pub fn exponential_powers<F: PrimeField>(b: F, t: usize) -> Vec<F> {
@@ -70,6 +71,40 @@ .map(|(beta_i, delta_i)| beta_i + alpha * delta_i)
         .collect::<Vec<FpVar<F>>>()
 }
 
+/// Returns the product of selected elements in `betas`.
+/// For every index `j`, whether `betas[j]` is selected is determined by the
+/// `j`-th bit in the binary (little endian) representation of `i`.
+///
+/// If `betas = (β, β^2, β^4, ..., β^{2^{t-1}})`, then the result is equal to
+/// `β^i`.
+pub fn pow_i<F: PrimeField>(mut i: usize, betas: &[F]) -> F {
+    let mut j = 0;
+    let mut r = F::one();
+    while i > 0 {
+        if i.is_odd() {
+            r *= betas[j];
+        }
+        i >>= 1;
+        j += 1;
+    }
+    r
+}
+
+/// The in-circuit version of `pow_i`
+#[allow(dead_code)] // Will remove this once we have the decider circuit for Protogalaxy
+pub fn pow_i_var<F: PrimeField>(mut i: usize, betas: &[FpVar<F>]) -> FpVar<F> {
+    let mut j = 0;
+    let mut r = FieldVar::one();
+    while i > 0 {
+        if i.is_odd() {
+            r *= &betas[j];
+        }
+        i >>= 1;
+        j += 1;
+    }
+    r
+}
+
 #[cfg(test)]
 mod tests {
     use std::error::Error;
@@ -78,6 +113,7 @@ use ark_bn254::Fr;
     use ark_r1cs_std::{alloc::AllocVar, R1CSVar};
     use ark_relations::r1cs::ConstraintSystem;
     use ark_std::{test_rng, UniformRand};
+    use rand::Rng;
 
     use super::*;
 
@@ -138,6 +174,27 @@ let alpha_var = FpVar::new_witness(cs.clone(), || Ok(alpha))?;
 
             let r = betas_star(&betas, &deltas, alpha);
             let r_var = betas_star_var(&betas_var, &deltas_var, &alpha_var);
+            assert_eq!(r, r_var.value()?);
+            assert!(cs.is_satisfied()?);
+        }
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_pow_i() -> Result<(), Box<dyn Error>> {
+        let rng = &mut test_rng();
+
+        for t in 1..10 {
+            let cs = ConstraintSystem::<Fr>::new_ref();
+
+            let betas = (0..t).map(|_| Fr::rand(rng)).collect::<Vec<_>>();
+            let i = rng.gen_range(0..(1 << t));
+
+            let betas_var = Vec::new_witness(cs.clone(), || Ok(betas.clone()))?;
+
+            let r = pow_i(i, &betas);
+            let r_var = pow_i_var(i, &betas_var);
             assert_eq!(r, r_var.value()?);
             assert!(cs.is_satisfied()?);
         }

@@ -258,7 +258,7 @@ CircomFCircuit::<Fr>::new((r1cs_path.into(), wasm_path.into(), 1, 0)).unwrap(); // state_len:1, external_inputs_len:0
 
         // Allocates z_i1 by using step_native function.
         let z_i = vec![Fr::from(3_u32)];
-        let wrapper_circuit = crate::frontend::tests::WrapperCircuit {
+        let wrapper_circuit = crate::frontend::utils::WrapperCircuit {
             FC: circom_fcircuit.clone(),
             z_i: Some(z_i.clone()),
             z_i1: Some(circom_fcircuit.step_native(0, z_i.clone(), vec![]).unwrap()),
@@ -367,7 +367,7 @@ }));
 
         // Allocates z_i1 by using step_native function.
         let z_i = vec![Fr::from(3_u32)];
-        let wrapper_circuit = crate::frontend::tests::WrapperCircuit {
+        let wrapper_circuit = crate::frontend::utils::WrapperCircuit {
             FC: circom_fcircuit.clone(),
             z_i: Some(z_i.clone()),
             z_i1: Some(circom_fcircuit.step_native(0, z_i.clone(), vec![]).unwrap()),

@@ -7,6 +7,7 @@
 pub mod circom;
 pub mod noir;
 pub mod noname;
+pub mod utils;
 
 /// FCircuit defines the trait of the circuit of the F function, which is the one being folded (ie.
 /// inside the agmented F' function).
@@ -53,131 +54,9 @@ #[cfg(test)]
 pub mod tests {
     use super::*;
     use ark_bn254::Fr;
-    use ark_r1cs_std::{alloc::AllocVar, eq::EqGadget};
     use ark_relations::r1cs::{ConstraintSynthesizer, ConstraintSystem};
-    use core::marker::PhantomData;
 
-    /// CubicFCircuit is a struct that implements the FCircuit trait, for the R1CS example circuit
-    /// from https://www.vitalik.ca/general/2016/12/10/qap.html, which checks `x^3 + x + 5 = y`.
-    /// `z_i` is used as `x`, and `z_{i+1}` is used as `y`, and at the next step, `z_{i+1}` will be
-    /// assigned to `z_i`, and a new `z+{i+1}` will be computted.
-    #[derive(Clone, Copy, Debug)]
-    pub struct CubicFCircuit<F: PrimeField> {
-        _f: PhantomData<F>,
-    }
-    impl<F: PrimeField> FCircuit<F> for CubicFCircuit<F> {
-        type Params = ();
-        fn new(_params: Self::Params) -> Result<Self, Error> {
-            Ok(Self { _f: PhantomData })
-        }
-        fn state_len(&self) -> usize {
-            1
-        }
-        fn external_inputs_len(&self) -> usize {
-            0
-        }
-        fn step_native(
-            &self,
-            _i: usize,
-            z_i: Vec<F>,
-            _external_inputs: Vec<F>,
-        ) -> Result<Vec<F>, Error> {
-            Ok(vec![z_i[0] * z_i[0] * z_i[0] + z_i[0] + F::from(5_u32)])
-        }
-        fn generate_step_constraints(
-            &self,
-            cs: ConstraintSystemRef<F>,
-            _i: usize,
-            z_i: Vec<FpVar<F>>,
-            _external_inputs: Vec<FpVar<F>>,
-        ) -> Result<Vec<FpVar<F>>, SynthesisError> {
-            let five = FpVar::<F>::new_constant(cs.clone(), F::from(5u32))?;
-            let z_i = z_i[0].clone();
-
-            Ok(vec![&z_i * &z_i * &z_i + &z_i + &five])
-        }
-    }
-
-    /// CustomFCircuit is a circuit that has the number of constraints specified in the
-    /// `n_constraints` parameter. Note that the generated circuit will have very sparse matrices.
-    #[derive(Clone, Copy, Debug)]
-    pub struct CustomFCircuit<F: PrimeField> {
-        _f: PhantomData<F>,
-        pub n_constraints: usize,
-    }
-    impl<F: PrimeField> FCircuit<F> for CustomFCircuit<F> {
-        type Params = usize;
-
-        fn new(params: Self::Params) -> Result<Self, Error> {
-            Ok(Self {
-                _f: PhantomData,
-                n_constraints: params,
-            })
-        }
-        fn state_len(&self) -> usize {
-            1
-        }
-        fn external_inputs_len(&self) -> usize {
-            0
-        }
-        fn step_native(
-            &self,
-            _i: usize,
-            z_i: Vec<F>,
-            _external_inputs: Vec<F>,
-        ) -> Result<Vec<F>, Error> {
-            let mut z_i1 = F::one();
-            for _ in 0..self.n_constraints - 1 {
-                z_i1 *= z_i[0];
-            }
-            Ok(vec![z_i1])
-        }
-        fn generate_step_constraints(
-            &self,
-            cs: ConstraintSystemRef<F>,
-            _i: usize,
-            z_i: Vec<FpVar<F>>,
-            _external_inputs: Vec<FpVar<F>>,
-        ) -> Result<Vec<FpVar<F>>, SynthesisError> {
-            let mut z_i1 = FpVar::<F>::new_witness(cs.clone(), || Ok(F::one()))?;
-            for _ in 0..self.n_constraints - 1 {
-                z_i1 *= z_i[0].clone();
-            }
-
-            Ok(vec![z_i1])
-        }
-    }
-
-    /// WrapperCircuit is a circuit that wraps any circuit that implements the FCircuit trait. This
-    /// is used to test the `FCircuit.generate_step_constraints` method. This is a similar wrapping
-    /// than the one done in the `AugmentedFCircuit`, but without adding all the extra constraints
-    /// of the AugmentedF circuit logic, in order to run lighter tests when we're not interested in
-    /// the the AugmentedF logic but in the wrapping of the circuits.
-    pub struct WrapperCircuit<F: PrimeField, FC: FCircuit<F>> {
-        pub FC: FC, // F circuit
-        pub z_i: Option<Vec<F>>,
-        pub z_i1: Option<Vec<F>>,
-    }
-    impl<F, FC> ConstraintSynthesizer<F> for WrapperCircuit<F, FC>
-    where
-        F: PrimeField,
-        FC: FCircuit<F>,
-    {
-        fn generate_constraints(self, cs: ConstraintSystemRef<F>) -> Result<(), SynthesisError> {
-            let z_i = Vec::<FpVar<F>>::new_witness(cs.clone(), || {
-                Ok(self.z_i.unwrap_or(vec![F::zero()]))
-            })?;
-            let z_i1 = Vec::<FpVar<F>>::new_input(cs.clone(), || {
-                Ok(self.z_i1.unwrap_or(vec![F::zero()]))
-            })?;
-            let computed_z_i1 =
-                self.FC
-                    .generate_step_constraints(cs.clone(), 0, z_i.clone(), vec![])?;
-
-            computed_z_i1.enforce_equal(&z_i1)?;
-            Ok(())
-        }
-    }
+    use utils::{CubicFCircuit, CustomFCircuit, WrapperCircuit};
 
     #[test]
     fn test_testfcircuit() {

@@ -0,0 +1,181 @@
+use ark_ff::PrimeField;
+use ark_r1cs_std::{alloc::AllocVar, fields::fp::FpVar};
+use ark_relations::r1cs::{ConstraintSystemRef, SynthesisError};
+#[cfg(test)]
+use ark_std::marker::PhantomData;
+use ark_std::{fmt::Debug, Zero};
+
+use super::FCircuit;
+use crate::Error;
+
+/// DummyCircuit is a circuit that has dummy state and external inputs whose
+/// lengths are specified in the `state_len` and `external_inputs_len`
+/// parameters, without any constraints.
+#[derive(Clone, Debug)]
+pub struct DummyCircuit {
+    state_len: usize,
+    external_inputs_len: usize,
+}
+impl<F: PrimeField> FCircuit<F> for DummyCircuit {
+    type Params = (usize, usize);
+
+    fn new((state_len, external_inputs_len): Self::Params) -> Result<Self, Error> {
+        Ok(Self {
+            state_len,
+            external_inputs_len,
+        })
+    }
+    fn state_len(&self) -> usize {
+        self.state_len
+    }
+    fn external_inputs_len(&self) -> usize {
+        self.external_inputs_len
+    }
+    fn step_native(
+        &self,
+        _i: usize,
+        _z_i: Vec<F>,
+        _external_inputs: Vec<F>,
+    ) -> Result<Vec<F>, Error> {
+        Ok(vec![F::zero(); self.state_len])
+    }
+    fn generate_step_constraints(
+        &self,
+        cs: ConstraintSystemRef<F>,
+        _i: usize,
+        _z_i: Vec<FpVar<F>>,
+        _external_inputs: Vec<FpVar<F>>,
+    ) -> Result<Vec<FpVar<F>>, SynthesisError> {
+        Vec::new_witness(cs.clone(), || Ok(vec![Zero::zero(); self.state_len]))
+    }
+}
+
+/// CubicFCircuit is a struct that implements the FCircuit trait, for the R1CS example circuit
+/// from https://www.vitalik.ca/general/2016/12/10/qap.html, which checks `x^3 + x + 5 = y`.
+/// `z_i` is used as `x`, and `z_{i+1}` is used as `y`, and at the next step, `z_{i+1}` will be
+/// assigned to `z_i`, and a new `z+{i+1}` will be computted.
+#[cfg(test)]
+#[derive(Clone, Copy, Debug)]
+pub struct CubicFCircuit<F: PrimeField> {
+    _f: PhantomData<F>,
+}
+
+#[cfg(test)]
+impl<F: PrimeField> FCircuit<F> for CubicFCircuit<F> {
+    type Params = ();
+    fn new(_params: Self::Params) -> Result<Self, Error> {
+        Ok(Self { _f: PhantomData })
+    }
+    fn state_len(&self) -> usize {
+        1
+    }
+    fn external_inputs_len(&self) -> usize {
+        0
+    }
+    fn step_native(
+        &self,
+        _i: usize,
+        z_i: Vec<F>,
+        _external_inputs: Vec<F>,
+    ) -> Result<Vec<F>, Error> {
+        Ok(vec![z_i[0] * z_i[0] * z_i[0] + z_i[0] + F::from(5_u32)])
+    }
+    fn generate_step_constraints(
+        &self,
+        cs: ConstraintSystemRef<F>,
+        _i: usize,
+        z_i: Vec<FpVar<F>>,
+        _external_inputs: Vec<FpVar<F>>,
+    ) -> Result<Vec<FpVar<F>>, SynthesisError> {
+        let five = FpVar::<F>::new_constant(cs.clone(), F::from(5u32))?;
+        let z_i = z_i[0].clone();
+
+        Ok(vec![&z_i * &z_i * &z_i + &z_i + &five])
+    }
+}
+
+/// CustomFCircuit is a circuit that has the number of constraints specified in the
+/// `n_constraints` parameter. Note that the generated circuit will have very sparse matrices.
+#[cfg(test)]
+#[derive(Clone, Copy, Debug)]
+pub struct CustomFCircuit<F: PrimeField> {
+    _f: PhantomData<F>,
+    pub n_constraints: usize,
+}
+
+#[cfg(test)]
+impl<F: PrimeField> FCircuit<F> for CustomFCircuit<F> {
+    type Params = usize;
+
+    fn new(params: Self::Params) -> Result<Self, Error> {
+        Ok(Self {
+            _f: PhantomData,
+            n_constraints: params,
+        })
+    }
+    fn state_len(&self) -> usize {
+        1
+    }
+    fn external_inputs_len(&self) -> usize {
+        0
+    }
+    fn step_native(
+        &self,
+        _i: usize,
+        z_i: Vec<F>,
+        _external_inputs: Vec<F>,
+    ) -> Result<Vec<F>, Error> {
+        let mut z_i1 = F::one();
+        for _ in 0..self.n_constraints - 1 {
+            z_i1 *= z_i[0];
+        }
+        Ok(vec![z_i1])
+    }
+    fn generate_step_constraints(
+        &self,
+        cs: ConstraintSystemRef<F>,
+        _i: usize,
+        z_i: Vec<FpVar<F>>,
+        _external_inputs: Vec<FpVar<F>>,
+    ) -> Result<Vec<FpVar<F>>, SynthesisError> {
+        let mut z_i1 = FpVar::<F>::new_witness(cs.clone(), || Ok(F::one()))?;
+        for _ in 0..self.n_constraints - 1 {
+            z_i1 *= z_i[0].clone();
+        }
+
+        Ok(vec![z_i1])
+    }
+}
+
+/// WrapperCircuit is a circuit that wraps any circuit that implements the FCircuit trait. This
+/// is used to test the `FCircuit.generate_step_constraints` method. This is a similar wrapping
+/// than the one done in the `AugmentedFCircuit`, but without adding all the extra constraints
+/// of the AugmentedF circuit logic, in order to run lighter tests when we're not interested in
+/// the the AugmentedF logic but in the wrapping of the circuits.
+#[cfg(test)]
+pub struct WrapperCircuit<F: PrimeField, FC: FCircuit<F>> {
+    pub FC: FC, // F circuit
+    pub z_i: Option<Vec<F>>,
+    pub z_i1: Option<Vec<F>>,
+}
+
+#[cfg(test)]
+impl<F, FC> ark_relations::r1cs::ConstraintSynthesizer<F> for WrapperCircuit<F, FC>
+where
+    F: PrimeField,
+    FC: FCircuit<F>,
+{
+    fn generate_constraints(self, cs: ConstraintSystemRef<F>) -> Result<(), SynthesisError> {
+        let z_i =
+            Vec::<FpVar<F>>::new_witness(cs.clone(), || Ok(self.z_i.unwrap_or(vec![F::zero()])))?;
+        let z_i1 =
+            Vec::<FpVar<F>>::new_input(cs.clone(), || Ok(self.z_i1.unwrap_or(vec![F::zero()])))?;
+        let computed_z_i1 =
+            self.FC
+                .generate_step_constraints(cs.clone(), 0, z_i.clone(), vec![])?;
+
+        use ark_r1cs_std::eq::EqGadget;
+        computed_z_i1.enforce_equal(&z_i1)?;
+        Ok(())
+    }
+}

@@ -82,7 +82,7 @@ "b.len()".to_string(),
             b.len(),
         ));
     }
-    Ok(a.iter().zip(b.iter()).map(|(x, y)| *x + y).collect())
+    Ok(cfg_iter!(a).zip(b).map(|(x, y)| *x + y).collect())
 }
 
 pub fn vec_sub<F: PrimeField>(a: &[F], b: &[F]) -> Result<Vec<F>, Error> {
@@ -94,15 +94,15 @@ "b.len()".to_string(),
             b.len(),
         ));
     }
-    Ok(a.iter().zip(b.iter()).map(|(x, y)| *x - y).collect())
+    Ok(cfg_iter!(a).zip(b).map(|(x, y)| *x - y).collect())
 }
 
 pub fn vec_scalar_mul<F: PrimeField>(vec: &[F], c: &F) -> Vec<F> {
-    vec.iter().map(|a| *a * c).collect()
+    cfg_iter!(vec).map(|a| *a * c).collect()
 }
 
 pub fn is_zero_vec<F: PrimeField>(vec: &[F]) -> bool {
-    vec.iter().all(|a| a.is_zero())
+    cfg_iter!(vec).all(|a| a.is_zero())
 }
 
 pub fn mat_vec_mul_dense<F: PrimeField>(M: &[Vec<F>], z: &[F]) -> Result<Vec<F>, Error> {
@@ -118,13 +118,9 @@ z.len(),
         ));
     }
 
-    let mut r: Vec<F> = vec![F::zero(); M.len()];
-    for (i, M_i) in M.iter().enumerate() {
-        for (j, M_ij) in M_i.iter().enumerate() {
-            r[i] += *M_ij * z[j];
-        }
-    }
-    Ok(r)
+    Ok(cfg_iter!(M)
+        .map(|row| row.iter().zip(z).map(|(a, b)| *a * b).sum())
+        .collect())
 }
 
 pub fn mat_vec_mul<F: PrimeField>(M: &SparseMatrix<F>, z: &[F]) -> Result<Vec<F>, Error> {
@@ -136,13 +132,9 @@ "z.len()".to_string(),
             z.len(),
         ));
     }
-    let mut res = vec![F::zero(); M.n_rows];
-    for (row_i, row) in M.coeffs.iter().enumerate() {
-        for &(value, col_i) in row.iter() {
-            res[row_i] += value * z[col_i];
-        }
-    }
-    Ok(res)
+    Ok(cfg_iter!(M.coeffs)
+        .map(|row| row.iter().map(|(value, col_i)| *value * z[*col_i]).sum())
+        .collect())
 }
 
 pub fn mat_from_str_mat<F: PrimeField>(str_mat: Vec<Vec<&str>>) -> Result<Vec<Vec<F>>, Error> {