diff --git a/README.md b/README.md
index cc38015..0e592dd 100644
--- a/README.md
+++ b/README.md
@@ -21,10 +21,10 @@ Folding schemes implemented:
 
 - [Nova: Recursive Zero-Knowledge Arguments from Folding Schemes](https://eprint.iacr.org/2021/370.pdf), Abhiram Kothapalli, Srinath Setty, Ioanna Tzialla. 2021
 - [CycleFold: Folding-scheme-based recursive arguments over a cycle of elliptic curves](https://eprint.iacr.org/2023/1192.pdf), Abhiram Kothapalli, Srinath Setty. 2023
+- [HyperNova: Recursive arguments for customizable constraint systems](https://eprint.iacr.org/2023/573.pdf), Abhiram Kothapalli, Srinath Setty. 2023
 
 Work in progress:
 
-- [HyperNova: Recursive arguments for customizable constraint systems](https://eprint.iacr.org/2023/573.pdf), Abhiram Kothapalli, Srinath Setty. 2023
 - [ProtoGalaxy: Efficient ProtoStar-style folding of multiple instances](https://eprint.iacr.org/2023/1106.pdf), Liam Eagen, Ariel Gabizon. 2023
 
 ## Available frontends
diff --git a/examples/circom_full_flow.rs b/examples/circom_full_flow.rs
index 69d902f..35a3f49 100644
--- a/examples/circom_full_flow.rs
+++ b/examples/circom_full_flow.rs
@@ -85,7 +85,7 @@ fn main() {
     let nova_params = N::preprocess(&mut rng, &nova_preprocess_params).unwrap();
 
     // initialize the folding scheme engine, in our case we use Nova
-    let mut nova = N::init(nova_params.clone(), f_circuit.clone(), z_0).unwrap();
+    let mut nova = N::init(&nova_params, f_circuit.clone(), z_0).unwrap();
 
     // prepare the Decider prover & verifier params
     let (decider_pp, decider_vp) = D::preprocess(&mut rng, &nova_params, nova.clone()).unwrap();
@@ -93,7 +93,7 @@ fn main() {
     // run n steps of the folding iteration
     for (i, external_inputs_at_step) in external_inputs.iter().enumerate() {
         let start = Instant::now();
-        nova.prove_step(rng, external_inputs_at_step.clone())
+        nova.prove_step(rng, external_inputs_at_step.clone(), None)
             .unwrap();
         println!("Nova::prove_step {}: {:?}", i, start.elapsed());
     }
diff --git a/examples/external_inputs.rs b/examples/external_inputs.rs
index 726e0e9..4939b32 100644
--- a/examples/external_inputs.rs
+++ b/examples/external_inputs.rs
@@ -190,14 +190,13 @@ fn main() {
     let nova_params = N::preprocess(&mut rng, &nova_preprocess_params).unwrap();
 
     println!("Initialize FoldingScheme");
-    let mut folding_scheme =
-        N::init(nova_params.clone(), F_circuit, initial_state.clone()).unwrap();
+    let mut folding_scheme = N::init(&nova_params, F_circuit, initial_state.clone()).unwrap();
 
     // compute a step of the IVC
     for (i, external_inputs_at_step) in external_inputs.iter().enumerate() {
         let start = Instant::now();
         folding_scheme
-            .prove_step(rng, external_inputs_at_step.clone())
+            .prove_step(rng, external_inputs_at_step.clone(), None)
             .unwrap();
         println!("Nova::prove_step {}: {:?}", i, start.elapsed());
     }
diff --git a/examples/full_flow.rs b/examples/full_flow.rs
index 2304b59..6db6733 100644
--- a/examples/full_flow.rs
+++ b/examples/full_flow.rs
@@ -102,7 +102,7 @@ fn main() {
     let nova_params = N::preprocess(&mut rng, &nova_preprocess_params).unwrap();
 
     // initialize the folding scheme engine, in our case we use Nova
-    let mut nova = N::init(nova_params.clone(), f_circuit, z_0).unwrap();
+    let mut nova = N::init(&nova_params, f_circuit, z_0).unwrap();
 
     // prepare the Decider prover & verifier params
     let (decider_pp, decider_vp) = D::preprocess(&mut rng, &nova_params, nova.clone()).unwrap();
@@ -110,7 +110,7 @@ fn main() {
     // run n steps of the folding iteration
     for i in 0..n_steps {
         let start = Instant::now();
-        nova.prove_step(rng, vec![]).unwrap();
+        nova.prove_step(rng, vec![], None).unwrap();
         println!("Nova::prove_step {}: {:?}", i, start.elapsed());
     }
 
diff --git a/examples/multi_inputs.rs b/examples/multi_inputs.rs
index a78b98b..5560e74 100644
--- a/examples/multi_inputs.rs
+++ b/examples/multi_inputs.rs
@@ -144,13 +144,12 @@ fn main() {
     let nova_params = N::preprocess(&mut rng, &nova_preprocess_params).unwrap();
 
     println!("Initialize FoldingScheme");
-    let mut folding_scheme =
-        N::init(nova_params.clone(), F_circuit, initial_state.clone()).unwrap();
+    let mut folding_scheme = N::init(&nova_params, F_circuit, initial_state.clone()).unwrap();
 
     // compute a step of the IVC
     for i in 0..num_steps {
         let start = Instant::now();
-        folding_scheme.prove_step(rng, vec![]).unwrap();
+        folding_scheme.prove_step(rng, vec![], None).unwrap();
         println!("Nova::prove_step {}: {:?}", i, start.elapsed());
     }
 
diff --git a/examples/noname_full_flow.rs b/examples/noname_full_flow.rs
index 6d84aa0..a095475 100644
--- a/examples/noname_full_flow.rs
+++ b/examples/noname_full_flow.rs
@@ -86,7 +86,7 @@ fn main() {
     let nova_params = N::preprocess(&mut rng, &nova_preprocess_params).unwrap();
 
     // initialize the folding scheme engine, in our case we use Nova
-    let mut nova = N::init(nova_params.clone(), f_circuit.clone(), z_0).unwrap();
+    let mut nova = N::init(&nova_params, f_circuit.clone(), z_0).unwrap();
 
     // prepare the Decider prover & verifier params
     let (decider_pp, decider_vp) = D::preprocess(&mut rng, &nova_params, nova.clone()).unwrap();
@@ -94,7 +94,7 @@ fn main() {
     // run n steps of the folding iteration
     for (i, external_inputs_at_step) in external_inputs.iter().enumerate() {
         let start = Instant::now();
-        nova.prove_step(rng, external_inputs_at_step.clone())
+        nova.prove_step(rng, external_inputs_at_step.clone(), None)
             .unwrap();
         println!("Nova::prove_step {}: {:?}", i, start.elapsed());
     }
diff --git a/examples/sha256.rs b/examples/sha256.rs
index 045967b..eea1887 100644
--- a/examples/sha256.rs
+++ b/examples/sha256.rs
@@ -129,12 +129,11 @@ fn main() {
     let nova_params = N::preprocess(&mut rng, &nova_preprocess_params).unwrap();
 
     println!("Initialize FoldingScheme");
-    let mut folding_scheme =
-        N::init(nova_params.clone(), F_circuit, initial_state.clone()).unwrap();
+    let mut folding_scheme = N::init(&nova_params, F_circuit, initial_state.clone()).unwrap();
     // compute a step of the IVC
     for i in 0..num_steps {
         let start = Instant::now();
-        folding_scheme.prove_step(rng, vec![]).unwrap();
+        folding_scheme.prove_step(rng, vec![], None).unwrap();
         println!("Nova::prove_step {}: {:?}", i, start.elapsed());
     }
 
diff --git a/folding-schemes/src/arith/ccs.rs b/folding-schemes/src/arith/ccs.rs
index 80b812f..3dd6b87 100644
--- a/folding-schemes/src/arith/ccs.rs
+++ b/folding-schemes/src/arith/ccs.rs
@@ -67,7 +67,7 @@ impl<F: PrimeField> Arith<F> for CCS<F> {
         Ok(())
     }
 
-    fn params_to_bytes(&self) -> Vec<u8> {
+    fn params_to_le_bytes(&self) -> Vec<u8> {
         [
             self.l.to_le_bytes(),
             self.m.to_le_bytes(),
diff --git a/folding-schemes/src/arith/mod.rs b/folding-schemes/src/arith/mod.rs
index f53b565..e09746d 100644
--- a/folding-schemes/src/arith/mod.rs
+++ b/folding-schemes/src/arith/mod.rs
@@ -11,5 +11,5 @@ pub trait Arith<F: PrimeField> {
 
     /// Returns the bytes that represent the parameters, that is, the matrices sizes, the amount of
     /// public inputs, etc, without the matrices/polynomials values.
-    fn params_to_bytes(&self) -> Vec<u8>;
+    fn params_to_le_bytes(&self) -> Vec<u8>;
 }
diff --git a/folding-schemes/src/arith/r1cs.rs b/folding-schemes/src/arith/r1cs.rs
index ef0c5ba..a09e957 100644
--- a/folding-schemes/src/arith/r1cs.rs
+++ b/folding-schemes/src/arith/r1cs.rs
@@ -28,7 +28,7 @@ impl<F: PrimeField> Arith<F> for R1CS<F> {
         Ok(())
     }
 
-    fn params_to_bytes(&self) -> Vec<u8> {
+    fn params_to_le_bytes(&self) -> Vec<u8> {
         [
             self.l.to_le_bytes(),
             self.A.n_rows.to_le_bytes(),
diff --git a/folding-schemes/src/folding/circuits/cyclefold.rs b/folding-schemes/src/folding/circuits/cyclefold.rs
index 522971e..4a22eb2 100644
--- a/folding-schemes/src/folding/circuits/cyclefold.rs
+++ b/folding-schemes/src/folding/circuits/cyclefold.rs
@@ -28,8 +28,13 @@ use crate::frontend::FCircuit;
 use crate::transcript::{AbsorbNonNativeGadget, Transcript, TranscriptVar};
 use crate::Error;
 
-/// Public inputs length for the CycleFoldCircuit: |[r, p1.x,y, p2.x,y, p3.x,y]|
-pub(crate) const CF_IO_LEN: usize = 7;
+/// Public inputs length for the CycleFoldCircuit:
+/// For Nova this is: |[r, p1.x,y, p2.x,y, p3.x,y]|
+/// In general, |[r * (n_points-1), (p_i.x,y)*n_points, p_folded.x,y]|, thus, io len is:
+/// (n_points-1) + 2*n_points + 2
+pub fn cf_io_len(n_points: usize) -> usize {
+    (n_points - 1) + 2 * n_points + 2
+}
 
 /// CycleFoldCommittedInstanceVar is the CycleFold CommittedInstance representation in the Nova
 /// circuit.
@@ -280,18 +285,24 @@ where
 #[derive(Debug, Clone)]
 pub struct CycleFoldCircuit<C: CurveGroup, GC: CurveVar<C, CF2<C>>> {
     pub _gc: PhantomData<GC>,
-    pub r_bits: Option<Vec<bool>>,
-    pub p1: Option<C>,
-    pub p2: Option<C>,
+    /// number of points being folded
+    pub n_points: usize,
+    /// r_bits is a vector containing the r_bits, one for each point except for the first one. They
+    /// are used for the scalar multiplication of the points. The r_bits are the bit
+    /// representation of each power of r (in Fr, while the CycleFoldCircuit is in Fq).
+    pub r_bits: Option<Vec<Vec<bool>>>,
+    /// points to be folded in the CycleFoldCircuit
+    pub points: Option<Vec<C>>,
     pub x: Option<Vec<CF2<C>>>, // public inputs (cf_u_{i+1}.x)
 }
 impl<C: CurveGroup, GC: CurveVar<C, CF2<C>>> CycleFoldCircuit<C, GC> {
-    pub fn empty() -> Self {
+    /// n_points indicates the number of points being folded in the CycleFoldCircuit
+    pub fn empty(n_points: usize) -> Self {
         Self {
             _gc: PhantomData,
+            n_points,
             r_bits: None,
-            p1: None,
-            p2: None,
+            points: None,
             x: None,
         }
     }
@@ -304,33 +315,64 @@ where
     for<'a> &'a GC: GroupOpsBounds<'a, C, GC>,
 {
     fn generate_constraints(self, cs: ConstraintSystemRef<CF2<C>>) -> Result<(), SynthesisError> {
-        let r_bits: Vec<Boolean<CF2<C>>> = Vec::new_witness(cs.clone(), || {
-            Ok(self.r_bits.unwrap_or(vec![false; N_BITS_RO]))
+        let r_bits: Vec<Vec<Boolean<CF2<C>>>> = self
+            .r_bits
+            // n_points-1, bcs is one for each point except for the first one
+            .unwrap_or(vec![vec![false; N_BITS_RO]; self.n_points - 1])
+            .iter()
+            .map(|r_bits_i| {
+                Vec::<Boolean<CF2<C>>>::new_witness(cs.clone(), || Ok(r_bits_i.clone()))
+            })
+            .collect::<Result<Vec<Vec<Boolean<CF2<C>>>>, SynthesisError>>()?;
+        let points = Vec::<GC>::new_witness(cs.clone(), || {
+            Ok(self.points.unwrap_or(vec![C::zero(); self.n_points]))
         })?;
-        let p1 = GC::new_witness(cs.clone(), || Ok(self.p1.unwrap_or(C::zero())))?;
-        let p2 = GC::new_witness(cs.clone(), || Ok(self.p2.unwrap_or(C::zero())))?;
-        // Fold the original Nova instances natively in CycleFold
-        // For the cmW we're computing: U_i1.cmW = U_i.cmW + r * u_i.cmW
-        // For the cmE we're computing: U_i1.cmE = U_i.cmE + r * cmT + r^2 * u_i.cmE, where u_i.cmE
+
+        #[cfg(test)]
+        {
+            assert_eq!(self.n_points, points.len());
+            assert_eq!(self.n_points - 1, r_bits.len());
+        }
+
+        // Fold the original points of the instances natively in CycleFold.
+        // In Nova,
+        // - for the cmW we're computing: U_i1.cmW = U_i.cmW + r * u_i.cmW
+        // - for the cmE we're computing: U_i1.cmE = U_i.cmE + r * cmT + r^2 * u_i.cmE, where u_i.cmE
         // is assumed to be 0, so, U_i1.cmE = U_i.cmE + r * cmT
-        let p3 = &p1 + p2.scalar_mul_le(r_bits.iter())?;
+        let mut p_folded: GC = points[0].clone();
+        // iter over n_points-1 because the first point is not multiplied by r^i (it is multiplied
+        // by r^0=1)
+        for i in 0..self.n_points - 1 {
+            p_folded += points[i + 1].scalar_mul_le(r_bits[i].iter())?;
+        }
 
         let x = Vec::<FpVar<CF2<C>>>::new_input(cs.clone(), || {
-            Ok(self.x.unwrap_or(vec![CF2::<C>::zero(); CF_IO_LEN]))
+            Ok(self
+                .x
+                .unwrap_or(vec![CF2::<C>::zero(); cf_io_len(self.n_points)]))
         })?;
         #[cfg(test)]
-        assert_eq!(x.len(), CF_IO_LEN); // non-constrained sanity check
-
-        // check that the points coordinates are placed as the public input x: x == [r, p1, p2, p3]
-        let r: FpVar<CF2<C>> = Boolean::le_bits_to_fp_var(&r_bits)?;
-        let points_coords: Vec<FpVar<CF2<C>>> = [
-            vec![r],
-            p1.to_constraint_field()?[..2].to_vec(),
-            p2.to_constraint_field()?[..2].to_vec(),
-            p3.to_constraint_field()?[..2].to_vec(),
+        assert_eq!(x.len(), cf_io_len(self.n_points)); // non-constrained sanity check
+
+        // Check that the points coordinates are placed as the public input x:
+        // 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_0, r_1, r_2, ...,  r_n, p_0, p_1, p_2, ..., p_n, p_folded],
+        // where each p_i is in fact p_i.to_constraint_field()
+        let computed_x: Vec<FpVar<CF2<C>>> = [
+            r_bits
+                .iter()
+                .map(|r_bits_i| Boolean::le_bits_to_fp_var(r_bits_i))
+                .collect::<Result<Vec<FpVar<CF2<C>>>, SynthesisError>>()?,
+            points
+                .iter()
+                .map(|p_i| Ok(p_i.to_constraint_field()?[..2].to_vec()))
+                .collect::<Result<Vec<Vec<FpVar<CF2<C>>>>, SynthesisError>>()?
+                .concat(),
+            p_folded.to_constraint_field()?[..2].to_vec(),
         ]
         .concat();
-        points_coords.enforce_equal(&x)?;
+        computed_x.enforce_equal(&x)?;
 
         Ok(())
     }
@@ -341,6 +383,7 @@ where
 #[allow(clippy::type_complexity)]
 #[allow(clippy::too_many_arguments)]
 pub fn fold_cyclefold_circuit<C1, GC1, C2, GC2, FC, CS1, CS2>(
+    _n_points: usize,
     transcript: &mut impl Transcript<C1::ScalarField>,
     cf_r1cs: R1CS<C2::ScalarField>,
     cf_cs_params: CS2::ProverParams,
@@ -386,7 +429,7 @@ where
     }
 
     #[cfg(test)]
-    assert_eq!(cf_x_i.len(), CF_IO_LEN);
+    assert_eq!(cf_x_i.len(), cf_io_len(_n_points));
 
     // fold cyclefold instances
     let cf_w_i = Witness::<C2>::new(cf_w_i.clone(), cf_r1cs.A.n_rows);
@@ -473,9 +516,9 @@ pub mod tests {
         .concat();
         let cfW_circuit = CycleFoldCircuit::<Projective, GVar> {
             _gc: PhantomData,
-            r_bits: Some(r_bits.clone()),
-            p1: Some(ci1.clone().cmW),
-            p2: Some(ci2.clone().cmW),
+            n_points: 2,
+            r_bits: Some(vec![r_bits.clone()]),
+            points: Some(vec![ci1.clone().cmW, ci2.clone().cmW]),
             x: Some(cfW_u_i_x.clone()),
         };
         cfW_circuit.generate_constraints(cs.clone()).unwrap();
@@ -492,9 +535,9 @@ pub mod tests {
         .concat();
         let cfE_circuit = CycleFoldCircuit::<Projective, GVar> {
             _gc: PhantomData,
-            r_bits: Some(r_bits.clone()),
-            p1: Some(ci1.clone().cmE),
-            p2: Some(cmT),
+            n_points: 2,
+            r_bits: Some(vec![r_bits.clone()]),
+            points: Some(vec![ci1.clone().cmE, cmT]),
             x: Some(cfE_u_i_x.clone()),
         };
         cfE_circuit.generate_constraints(cs.clone()).unwrap();
@@ -550,7 +593,7 @@ pub mod tests {
             u: Fr::zero(),
             cmW: Projective::rand(&mut rng),
             x: std::iter::repeat_with(|| Fr::rand(&mut rng))
-                .take(CF_IO_LEN)
+                .take(7) // 7 = cf_io_len
                 .collect(),
         };
         let U_i = CommittedInstance::<Projective> {
@@ -558,7 +601,7 @@ pub mod tests {
             u: Fr::rand(&mut rng),
             cmW: Projective::rand(&mut rng),
             x: std::iter::repeat_with(|| Fr::rand(&mut rng))
-                .take(CF_IO_LEN)
+                .take(7) // 7 = cf_io_len
                 .collect(),
         };
         let cmT = Projective::rand(&mut rng);
@@ -617,7 +660,7 @@ pub mod tests {
             u: Fr::rand(&mut rng),
             cmW: Projective::rand(&mut rng),
             x: std::iter::repeat_with(|| Fr::rand(&mut rng))
-                .take(CF_IO_LEN)
+                .take(7) // 7 = cf_io_len in Nova
                 .collect(),
         };
         let pp_hash = Fq::from(42u32); // only for test
diff --git a/folding-schemes/src/folding/hypernova/circuits.rs b/folding-schemes/src/folding/hypernova/circuits.rs
index 17b0442..78a28d2 100644
--- a/folding-schemes/src/folding/hypernova/circuits.rs
+++ b/folding-schemes/src/folding/hypernova/circuits.rs
@@ -14,7 +14,7 @@ use ark_r1cs_std::{
     fields::{fp::FpVar, FieldVar},
     groups::GroupOpsBounds,
     prelude::CurveVar,
-    R1CSVar, ToConstraintFieldGadget,
+    R1CSVar, ToBitsGadget, ToConstraintFieldGadget,
 };
 use ark_relations::r1cs::{
     ConstraintSynthesizer, ConstraintSystem, ConstraintSystemRef, Namespace, SynthesisError,
@@ -31,7 +31,7 @@ use super::{
 use crate::constants::N_BITS_RO;
 use crate::folding::{
     circuits::cyclefold::{
-        CycleFoldChallengeGadget, CycleFoldCommittedInstanceVar, NIFSFullGadget, CF_IO_LEN,
+        cf_io_len, CycleFoldChallengeGadget, CycleFoldCommittedInstanceVar, NIFSFullGadget,
     },
     circuits::{
         nonnative::{affine::NonNativeAffineVar, uint::NonNativeUintVar},
@@ -221,18 +221,19 @@ where
     <C as CurveGroup>::BaseField: PrimeField,
 {
     /// Runs (in-circuit) the NIMFS.V, which outputs the new folded LCCCS instance together with
-    /// the rho_bits, which will be used in other parts of the AugmentedFCircuit
+    /// the rho_powers, which will be used in other parts of the AugmentedFCircuit
     #[allow(clippy::type_complexity)]
     pub fn verify<S: CryptographicSponge, T: TranscriptVar<C::ScalarField, S>>(
         cs: ConstraintSystemRef<CF1<C>>,
         // only used the CCS params, not the matrices
         ccs: &CCS<C::ScalarField>,
         transcript: &mut T,
-        running_instances: &[LCCCSVar<C>],
-        new_instances: &[CCCSVar<C>],
+
+        running_instances: &[LCCCSVar<C>], // U
+        new_instances: &[CCCSVar<C>],      // u
         proof: ProofVar<C>,
         enabled: Boolean<C::ScalarField>,
-    ) -> Result<(LCCCSVar<C>, Vec<Boolean<CF1<C>>>), SynthesisError> {
+    ) -> Result<(LCCCSVar<C>, Vec<Vec<Boolean<CF1<C>>>>), SynthesisError> {
         // absorb instances to transcript
         for U_i in running_instances {
             let v = [
@@ -315,18 +316,15 @@ where
         let rho_bits: Vec<Boolean<CF1<C>>> = transcript.get_challenge_nbits(N_BITS_RO)?;
         let rho = Boolean::le_bits_to_fp_var(&rho_bits)?;
 
-        // return the folded instance, together with the rho_bits so they can be used in other
-        // parts of the AugmentedFCircuit
-        Ok((
-            Self::fold(
-                running_instances,
-                new_instances,
-                proof.sigmas_thetas,
-                r_x_prime,
-                rho,
-            )?,
-            rho_bits,
-        ))
+        // Self::fold will return the folded instance, together with the rho's powers vector so
+        // they can be used in other parts of the AugmentedFCircuit
+        Self::fold(
+            running_instances,
+            new_instances,
+            proof.sigmas_thetas,
+            r_x_prime,
+            rho,
+        )
     }
 
     /// Runs (in-circuit) the verifier side of the fold, computing the new folded LCCCS instance
@@ -337,12 +335,14 @@ where
         sigmas_thetas: (Vec<Vec<FpVar<CF1<C>>>>, Vec<Vec<FpVar<CF1<C>>>>),
         r_x_prime: Vec<FpVar<CF1<C>>>,
         rho: FpVar<CF1<C>>,
-    ) -> Result<LCCCSVar<C>, SynthesisError> {
+    ) -> Result<(LCCCSVar<C>, Vec<Vec<Boolean<CF1<C>>>>), SynthesisError> {
         let (sigmas, thetas) = (sigmas_thetas.0.clone(), sigmas_thetas.1.clone());
         let mut u_folded: FpVar<CF1<C>> = FpVar::zero();
         let mut x_folded: Vec<FpVar<CF1<C>>> = vec![FpVar::zero(); lcccs[0].x.len()];
         let mut v_folded: Vec<FpVar<CF1<C>>> = vec![FpVar::zero(); sigmas[0].len()];
 
+        let mut rho_vec: Vec<Vec<Boolean<CF1<C>>>> =
+            vec![vec![Boolean::FALSE; N_BITS_RO]; lcccs.len() + cccs.len() - 1];
         let mut rho_i = FpVar::one();
         for i in 0..(lcccs.len() + cccs.len()) {
             let u: FpVar<CF1<C>>;
@@ -379,16 +379,30 @@ where
                 .map(|(a_i, b_i)| a_i + b_i)
                 .collect();
 
+            // compute the next power of rho
             rho_i *= rho.clone();
+            // crop the size of rho_i to N_BITS_RO
+            let rho_i_bits = rho_i.to_bits_le()?;
+            rho_i = Boolean::le_bits_to_fp_var(&rho_i_bits[..N_BITS_RO])?;
+            if i < lcccs.len() + cccs.len() - 1 {
+                // store the cropped rho_i into the rho_vec
+                rho_vec[i] = rho_i_bits[..N_BITS_RO].to_vec();
+            }
         }
 
-        Ok(LCCCSVar::<C> {
-            C: lcccs[0].C.clone(), // WIP this will come from the cyclefold circuit
-            u: u_folded,
-            x: x_folded,
-            r_x: r_x_prime,
-            v: v_folded,
-        })
+        // return the folded instance, together with the rho's powers vector so they can be used in
+        // other parts of the AugmentedFCircuit
+        Ok((
+            LCCCSVar::<C> {
+                // C this is later overwritten by the U_{i+1}.C value checked by the cyclefold circuit
+                C: lcccs[0].C.clone(),
+                u: u_folded,
+                x: x_folded,
+                r_x: r_x_prime,
+                v: v_folded,
+            },
+            rho_vec,
+        ))
     }
 }
 
@@ -456,16 +470,20 @@ pub struct AugmentedFCircuit<
     pub poseidon_config: PoseidonConfig<CF1<C1>>,
     pub ccs: CCS<C1::ScalarField>, // CCS of the AugmentedFCircuit
     pub pp_hash: Option<CF1<C1>>,
+    pub mu: usize, // max number of LCCCS instances to be folded
+    pub nu: usize, // max number of CCCS instances to be folded
     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_C: Option<C1>, // u_i.C
     pub U_i: Option<LCCCS<C1>>,
-    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 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>>,
 
     // cyclefold verifier on C1
@@ -492,20 +510,29 @@ where
         poseidon_config: &PoseidonConfig<CF1<C1>>,
         F_circuit: FC,
         ccs: CCS<C1::ScalarField>,
-    ) -> Self {
-        Self {
+        mu: usize,
+        nu: usize,
+    ) -> Result<Self, Error> {
+        if mu < 1 || nu < 1 {
+            return Err(Error::CantBeZero("mu,nu".to_string()));
+        }
+        Ok(Self {
             _c2: PhantomData,
             _gc2: PhantomData,
             poseidon_config: poseidon_config.clone(),
             ccs,
             pp_hash: None,
+            mu,
+            nu,
             i: None,
             i_usize: None,
             z_0: None,
             z_i: None,
             external_inputs: None,
-            u_i_C: None,
             U_i: None,
+            Us: None,
+            u_i_C: None,
+            us: None,
             U_i1_C: None,
             F: F_circuit,
             x: None,
@@ -514,13 +541,15 @@ where
             cf_U_i: None,
             cf_x: None,
             cf_cmT: None,
-        }
+        })
     }
 
     pub fn empty(
         poseidon_config: &PoseidonConfig<CF1<C1>>,
-        F_circuit: FC,
+        F: FC, // FCircuit
         ccs: Option<CCS<C1::ScalarField>>,
+        mu: usize,
+        nu: usize,
     ) -> Result<Self, Error> {
         let initial_ccs = CCS {
             // m, n, s, s_prime and M will be overwritten by the `upper_bound_ccs' method
@@ -536,7 +565,7 @@ where
             c: vec![C1::ScalarField::one(), C1::ScalarField::one().neg()],
             M: vec![],
         };
-        let mut augmented_f_circuit = Self::default(poseidon_config, F_circuit, initial_ccs);
+        let mut augmented_f_circuit = Self::default(poseidon_config, F, initial_ccs, mu, nu)?;
         if ccs.is_some() {
             augmented_f_circuit.ccs = ccs.unwrap();
         } else {
@@ -551,51 +580,85 @@ where
     /// feed in as parameter for the AugmentedFCircuit::empty method to avoid computing them there.
     pub fn upper_bound_ccs(&self) -> Result<CCS<C1::ScalarField>, Error> {
         let r1cs = get_r1cs_from_cs::<CF1<C1>>(self.clone()).unwrap();
-        let ccs = CCS::from_r1cs(r1cs.clone());
+        let mut ccs = CCS::from_r1cs(r1cs.clone());
 
         let z_0 = vec![C1::ScalarField::zero(); self.F.state_len()];
-        let W_i = Witness::<C1::ScalarField>::dummy(&ccs);
-        let U_i = LCCCS::<C1>::dummy(ccs.l, ccs.t, ccs.s);
-        let w_i = W_i.clone();
-        let u_i = CCCS::<C1>::dummy(ccs.l);
-
-        let mut transcript_p: PoseidonSponge<C1::ScalarField> =
-            PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config.clone());
-        // since this is only for the number of constraints, no need to absorb the pp_hash here
-        let (nimfs_proof, U_i1, _, _) = NIMFS::<C1, PoseidonSponge<C1::ScalarField>>::prove(
-            &mut transcript_p,
-            &ccs,
-            &[U_i.clone()],
-            &[u_i.clone()],
-            &[W_i.clone()],
-            &[w_i.clone()],
-        )?;
+        let mut W_i = Witness::<C1::ScalarField>::dummy(&ccs);
+        let mut U_i = LCCCS::<C1>::dummy(ccs.l, ccs.t, ccs.s);
+        let mut w_i = W_i.clone();
+        let mut u_i = CCCS::<C1>::dummy(ccs.l);
+
+        let n_iters = 2;
+        for _ in 0..n_iters {
+            let Us = vec![U_i.clone(); self.mu - 1];
+            let Ws = vec![W_i.clone(); self.mu - 1];
+            let us = vec![u_i.clone(); self.nu - 1];
+            let ws = vec![w_i.clone(); self.nu - 1];
+
+            let all_Us = [vec![U_i.clone()], Us.clone()].concat();
+            let all_us = [vec![u_i.clone()], us.clone()].concat();
+            let all_Ws = [vec![W_i.clone()], Ws].concat();
+            let all_ws = [vec![w_i.clone()], ws].concat();
+
+            let mut transcript_p: PoseidonSponge<C1::ScalarField> =
+                PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config.clone());
+            // since this is only for the number of constraints, no need to absorb the pp_hash here
+            let (nimfs_proof, U_i1, _, _) = NIMFS::<C1, PoseidonSponge<C1::ScalarField>>::prove(
+                &mut transcript_p,
+                &ccs,
+                &all_Us,
+                &all_us,
+                &all_Ws,
+                &all_ws,
+            )?;
 
-        let augmented_f_circuit = Self {
-            _c2: PhantomData,
-            _gc2: PhantomData,
-            poseidon_config: self.poseidon_config.clone(),
-            ccs: ccs.clone(),
-            pp_hash: Some(C1::ScalarField::zero()),
-            i: Some(C1::ScalarField::zero()),
-            i_usize: Some(0),
-            z_0: Some(z_0.clone()),
-            z_i: Some(z_0.clone()),
-            external_inputs: Some(vec![]),
-            u_i_C: Some(u_i.C),
-            U_i: Some(U_i.clone()),
-            U_i1_C: Some(U_i1.C),
-            F: self.F.clone(),
-            x: Some(C1::ScalarField::zero()),
-            nimfs_proof: Some(nimfs_proof),
-            // cyclefold values
-            cf_u_i_cmW: None,
-            cf_U_i: None,
-            cf_x: None,
-            cf_cmT: None,
-        };
+            let augmented_f_circuit = Self {
+                _c2: PhantomData,
+                _gc2: PhantomData,
+                poseidon_config: self.poseidon_config.clone(),
+                ccs: ccs.clone(),
+                pp_hash: Some(C1::ScalarField::zero()),
+                mu: self.mu,
+                nu: self.nu,
+                i: Some(C1::ScalarField::zero()),
+                i_usize: Some(0),
+                z_0: Some(z_0.clone()),
+                z_i: Some(z_0.clone()),
+                external_inputs: Some(vec![C1::ScalarField::zero(); self.F.external_inputs_len()]),
+                U_i: Some(U_i.clone()),
+                Us: Some(Us),
+                u_i_C: Some(u_i.C),
+                us: Some(us),
+                U_i1_C: Some(U_i1.C),
+                F: self.F.clone(),
+                x: Some(C1::ScalarField::zero()),
+                nimfs_proof: Some(nimfs_proof),
+                // cyclefold values
+                cf_u_i_cmW: None,
+                cf_U_i: None,
+                cf_x: None,
+                cf_cmT: None,
+            };
+
+            let cs: ConstraintSystem<C1::ScalarField>;
+            (cs, ccs) = augmented_f_circuit.compute_cs_ccs()?;
+            // prepare instances for next loop iteration
+            use crate::arith::r1cs::extract_w_x;
+            let (r1cs_w_i1, r1cs_x_i1) = extract_w_x::<C1::ScalarField>(&cs);
+            u_i = CCCS::<C1> {
+                C: u_i.C,
+                x: r1cs_x_i1,
+            };
+            w_i = Witness::<C1::ScalarField> {
+                w: r1cs_w_i1.clone(),
+                r_w: C1::ScalarField::one(),
+            };
+            W_i = Witness::<C1::ScalarField>::dummy(&ccs);
+            U_i = LCCCS::<C1>::dummy(ccs.l, ccs.t, ccs.s);
+        }
+        Ok(ccs)
 
-        Ok(augmented_f_circuit.compute_cs_ccs()?.1)
+        // Ok(augmented_f_circuit.compute_cs_ccs()?.1)
     }
 
     /// Returns the cs (ConstraintSystem) and the CCS out of the AugmentedFCircuit
@@ -651,20 +714,25 @@ where
         })?;
 
         let U_dummy = LCCCS::<C1>::dummy(self.ccs.l, self.ccs.t, self.ccs.s);
+        let u_dummy = CCCS::<C1>::dummy(self.ccs.l);
 
         let U_i =
             LCCCSVar::<C1>::new_witness(cs.clone(), || Ok(self.U_i.unwrap_or(U_dummy.clone())))?;
+        let Us = Vec::<LCCCSVar<C1>>::new_witness(cs.clone(), || {
+            Ok(self.Us.unwrap_or(vec![U_dummy.clone(); self.mu - 1]))
+        })?;
+        let us = Vec::<CCCSVar<C1>>::new_witness(cs.clone(), || {
+            Ok(self.us.unwrap_or(vec![u_dummy.clone(); self.mu - 1]))
+        })?;
         let U_i1_C = NonNativeAffineVar::new_witness(cs.clone(), || {
             Ok(self.U_i1_C.unwrap_or_else(C1::zero))
         })?;
-        let mu = 1; // Note: at this commit, only 2-to-1 instance fold is supported
-        let nu = 1;
-        let nimfs_proof_dummy = NIMFSProof::<C1>::dummy(&self.ccs, mu, nu);
+        let nimfs_proof_dummy = NIMFSProof::<C1>::dummy(&self.ccs, self.mu, self.nu);
         let nimfs_proof = ProofVar::<C1>::new_witness(cs.clone(), || {
             Ok(self.nimfs_proof.unwrap_or(nimfs_proof_dummy))
         })?;
 
-        let cf_u_dummy = CommittedInstance::dummy(CF_IO_LEN);
+        let cf_u_dummy = CommittedInstance::dummy(cf_io_len(self.mu + self.nu));
         let cf_U_i = CycleFoldCommittedInstanceVar::<C2, GC2>::new_witness(cs.clone(), || {
             Ok(self.cf_U_i.unwrap_or(cf_u_dummy.clone()))
         })?;
@@ -704,18 +772,21 @@ where
             x: vec![u_i_x, cf_u_i_x],
         };
 
+        let all_Us = [vec![U_i.clone()], Us].concat();
+        let all_us = [vec![u_i.clone()], us].concat();
+
         // P.3. NIMFS.verify, obtains U_{i+1} by folding [U_i] & [u_i].
         // Notice that NIMFSGadget::fold_committed_instance does not fold C. We set `U_i1.C` to
         // unconstrained witnesses `U_i1_C` respectively. Its correctness will be checked on the
         // other curve.
         let mut transcript = PoseidonSpongeVar::new(cs.clone(), &self.poseidon_config);
         transcript.absorb(&pp_hash)?;
-        let (mut U_i1, rho_bits) = NIMFSGadget::<C1>::verify(
+        let (mut U_i1, rho_vec) = NIMFSGadget::<C1>::verify(
             cs.clone(),
             &self.ccs.clone(),
             &mut transcript,
-            &[U_i.clone()],
-            &[u_i.clone()],
+            &all_Us,
+            &all_us,
             nimfs_proof,
             is_not_basecase.clone(),
         )?;
@@ -739,23 +810,36 @@ where
         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)?)?;
 
-        // convert rho_bits to a `NonNativeFieldVar`
-        let rho_nonnat = {
-            let mut bits = rho_bits;
-            bits.resize(C1::BaseField::MODULUS_BIT_SIZE as usize, Boolean::FALSE);
-            NonNativeUintVar::from(&bits)
-        };
+        // convert rho_bits of the rho_vec to a `NonNativeFieldVar`
+        let rho_vec_nonnat = rho_vec
+            .iter()
+            .map(|rho_i| {
+                let mut bits = rho_i.clone();
+                bits.resize(C1::BaseField::MODULUS_BIT_SIZE as usize, Boolean::FALSE);
+                NonNativeUintVar::from(&bits)
+            })
+            .collect();
 
         // CycleFold part
         // C.1. Compute cf1_u_i.x and cf2_u_i.x
-        let cf_x = vec![
-            rho_nonnat, U_i.C.x, U_i.C.y, u_i.C.x, u_i.C.y, U_i1.C.x, U_i1.C.y,
-        ];
+        let cf_x: Vec<NonNativeUintVar<CF2<C2>>> = [
+            rho_vec_nonnat,
+            all_Us
+                .iter()
+                .flat_map(|U| vec![U.C.x.clone(), U.C.y.clone()])
+                .collect(),
+            all_us
+                .iter()
+                .flat_map(|u| vec![u.C.x.clone(), u.C.y.clone()])
+                .collect(),
+            vec![U_i1.C.x, U_i1.C.y],
+        ]
+        .concat();
 
         // ensure that cf_u has as public inputs the C from main instances U_i, u_i, U_i+1
         // coordinates of the commitments.
         // C.2. Construct `cf_u_i`
-        let cf_u_i = CycleFoldCommittedInstanceVar {
+        let cf_u_i = CycleFoldCommittedInstanceVar::<C2, GC2> {
             // cf1_u_i.cmE = 0. Notice that we enforce cmE to be equal to 0 since it is allocated
             // as 0.
             cmE: GC2::zero(),
@@ -1091,6 +1175,9 @@ mod tests {
         let poseidon_config = poseidon_canonical_config::<Fr>();
         let sponge = PoseidonSponge::<Fr>::new(&poseidon_config);
 
+        let mu = 3;
+        let nu = 3;
+
         let start = Instant::now();
         let F_circuit = CubicFCircuit::<Fr>::new(()).unwrap();
         let mut augmented_f_circuit = AugmentedFCircuit::<
@@ -1098,7 +1185,7 @@ mod tests {
             Projective2,
             GVar2,
             CubicFCircuit<Fr>,
-        >::empty(&poseidon_config, F_circuit, None)
+        >::empty(&poseidon_config, F_circuit, None, mu, nu)
         .unwrap();
         let ccs = augmented_f_circuit.ccs.clone();
         println!("AugmentedFCircuit & CCS generation: {:?}", start.elapsed());
@@ -1106,7 +1193,7 @@ mod tests {
 
         // CycleFold circuit
         let cs2 = ConstraintSystem::<Fq>::new_ref();
-        let cf_circuit = CycleFoldCircuit::<Projective, GVar>::empty();
+        let cf_circuit = CycleFoldCircuit::<Projective, GVar>::empty(mu + nu);
         cf_circuit.generate_constraints(cs2.clone()).unwrap();
         cs2.finalize();
         let cs2 = cs2
@@ -1114,6 +1201,7 @@ mod tests {
             .ok_or(Error::NoInnerConstraintSystem)
             .unwrap();
         let cf_r1cs = extract_r1cs::<Fq>(&cs2);
+        println!("CF m x n: {} x {}", cf_r1cs.A.n_rows, cf_r1cs.A.n_cols);
 
         let (pedersen_params, _) =
             Pedersen::<Projective>::setup(&mut rng, ccs.n - ccs.l - 1).unwrap();
@@ -1152,6 +1240,16 @@ mod tests {
         for i in 0..n_steps {
             let start = Instant::now();
 
+            // for this test, let Us & us be just an array of copies of the U_i & u_i respectively
+            let Us = vec![U_i.clone(); mu - 1];
+            let Ws = vec![W_i.clone(); mu - 1];
+            let us = vec![u_i.clone(); nu - 1];
+            let ws = vec![w_i.clone(); nu - 1];
+            let all_Us = [vec![U_i.clone()], Us.clone()].concat();
+            let all_us = [vec![u_i.clone()], us.clone()].concat();
+            let all_Ws = [vec![W_i.clone()], Ws].concat();
+            let all_ws = [vec![w_i.clone()], ws].concat();
+
             let z_i1 = F_circuit.step_native(i, z_i.clone(), vec![]).unwrap();
 
             let (U_i1, W_i1);
@@ -1173,13 +1271,17 @@ mod tests {
                         poseidon_config: poseidon_config.clone(),
                         ccs: ccs.clone(),
                         pp_hash: Some(pp_hash),
+                        mu,
+                        nu,
                         i: Some(Fr::zero()),
                         i_usize: Some(0),
                         z_0: Some(z_0.clone()),
                         z_i: Some(z_i.clone()),
                         external_inputs: Some(vec![]),
-                        u_i_C: Some(u_i.C),
                         U_i: Some(U_i.clone()),
+                        Us: Some(Us.clone()),
+                        u_i_C: Some(u_i.C),
+                        us: Some(us.clone()),
                         U_i1_C: Some(U_i1.C),
                         F: F_circuit,
                         x: Some(u_i1_x),
@@ -1195,15 +1297,15 @@ mod tests {
                 let mut transcript_p: PoseidonSponge<Fr> =
                     PoseidonSponge::<Fr>::new(&poseidon_config.clone());
                 transcript_p.absorb(&pp_hash);
-                let (rho_bits, nimfs_proof);
-                (nimfs_proof, U_i1, W_i1, rho_bits) =
+                let (rho_powers, nimfs_proof);
+                (nimfs_proof, U_i1, W_i1, rho_powers) =
                     NIMFS::<Projective, PoseidonSponge<Fr>>::prove(
                         &mut transcript_p,
                         &ccs,
-                        &[U_i.clone()],
-                        &[u_i.clone()],
-                        &[W_i.clone()],
-                        &[w_i.clone()],
+                        &all_Us,
+                        &all_us,
+                        &all_Ws,
+                        &all_ws,
                     )
                     .unwrap();
 
@@ -1213,26 +1315,61 @@ mod tests {
                 let u_i1_x =
                     U_i1.hash(&sponge, pp_hash, iFr + Fr::one(), z_0.clone(), z_i1.clone());
 
-                let rho_Fq = Fq::from_bigint(BigInteger::from_bits_le(&rho_bits)).unwrap();
+                let rho_powers_Fq: Vec<Fq> = rho_powers
+                    .iter()
+                    .map(|rho_i| {
+                        Fq::from_bigint(BigInteger::from_bits_le(&rho_i.into_bigint().to_bits_le()))
+                            .unwrap()
+                    })
+                    .collect();
+                let rho_powers_bits: Vec<Vec<bool>> = rho_powers
+                    .iter()
+                    .map(|rho_i| rho_i.into_bigint().to_bits_le()[..N_BITS_RO].to_vec())
+                    .collect();
+
                 // CycleFold part:
                 // get the vector used as public inputs 'x' in the CycleFold circuit
-                // cyclefold circuit for cmW
                 let cf_u_i_x = [
-                    vec![rho_Fq],
+                    // all values for multiple instances
+                    rho_powers_Fq,
                     get_cm_coordinates(&U_i.C),
+                    Us.iter()
+                        .flat_map(|Us_i| get_cm_coordinates(&Us_i.C))
+                        .collect(),
                     get_cm_coordinates(&u_i.C),
+                    us.iter()
+                        .flat_map(|us_i| get_cm_coordinates(&us_i.C))
+                        .collect(),
                     get_cm_coordinates(&U_i1.C),
                 ]
                 .concat();
 
                 let cf_circuit = CycleFoldCircuit::<Projective, GVar> {
                     _gc: PhantomData,
-                    r_bits: Some(rho_bits.clone()),
-                    p1: Some(U_i.clone().C),
-                    p2: Some(u_i.clone().C),
+                    n_points: mu + nu,
+                    r_bits: Some(rho_powers_bits.clone()),
+                    points: Some(
+                        [
+                            vec![U_i.clone().C],
+                            Us.iter().map(|Us_i| Us_i.C).collect(),
+                            vec![u_i.clone().C],
+                            us.iter().map(|us_i| us_i.C).collect(),
+                        ]
+                        .concat(),
+                    ),
                     x: Some(cf_u_i_x.clone()),
                 };
 
+                // ensure that the CycleFoldCircuit is well defined
+                assert_eq!(
+                    cf_circuit.r_bits.clone().unwrap().len(),
+                    cf_circuit.n_points - 1
+                );
+                assert_eq!(
+                    cf_circuit.points.clone().unwrap().len(),
+                    cf_circuit.n_points
+                );
+
                 let (_cf_w_i, cf_u_i, cf_W_i1, cf_U_i1, cf_cmT, _) = fold_cyclefold_circuit::<
                     Projective,
                     GVar,
@@ -1242,6 +1379,7 @@ mod tests {
                     Pedersen<Projective>,
                     Pedersen<Projective2>,
                 >(
+                    mu + nu,
                     &mut transcript_p,
                     cf_r1cs.clone(),
                     cf_pedersen_params.clone(),
@@ -1264,13 +1402,17 @@ mod tests {
                         poseidon_config: poseidon_config.clone(),
                         ccs: ccs.clone(),
                         pp_hash: Some(pp_hash),
+                        mu,
+                        nu,
                         i: Some(iFr),
                         i_usize: Some(i),
                         z_0: Some(z_0.clone()),
                         z_i: Some(z_i.clone()),
                         external_inputs: Some(vec![]),
-                        u_i_C: Some(u_i.C),
                         U_i: Some(U_i.clone()),
+                        Us: Some(Us.clone()),
+                        u_i_C: Some(u_i.C),
+                        us: Some(us.clone()),
                         U_i1_C: Some(U_i1.C),
                         F: F_circuit,
                         x: Some(u_i1_x),
diff --git a/folding-schemes/src/folding/hypernova/mod.rs b/folding-schemes/src/folding/hypernova/mod.rs
index 0a3db8c..5a15655 100644
--- a/folding-schemes/src/folding/hypernova/mod.rs
+++ b/folding-schemes/src/folding/hypernova/mod.rs
@@ -21,11 +21,8 @@ use cccs::CCCS;
 use lcccs::LCCCS;
 use nimfs::NIMFS;
 
-use crate::arith::{
-    ccs::CCS,
-    r1cs::{extract_w_x, R1CS},
-};
 use crate::commitment::CommitmentScheme;
+use crate::constants::N_BITS_RO;
 use crate::folding::circuits::{
     cyclefold::{fold_cyclefold_circuit, CycleFoldCircuit},
     CF2,
@@ -37,7 +34,13 @@ use crate::folding::nova::{
 use crate::frontend::FCircuit;
 use crate::utils::{get_cm_coordinates, pp_hash};
 use crate::Error;
-use crate::FoldingScheme;
+use crate::{
+    arith::{
+        ccs::CCS,
+        r1cs::{extract_w_x, R1CS},
+    },
+    FoldingScheme, MultiFolding,
+};
 
 /// Witness for the LCCCS & CCCS, containing the w vector, and the r_w used as randomness in the Pedersen commitment.
 #[derive(Debug, Clone, Eq, PartialEq)]
@@ -70,6 +73,8 @@ where
     pub cf_cs_params: CS2::ProverParams,
     // if ccs is set, it will be used, if not, it will be computed at runtime
     pub ccs: Option<CCS<C1::ScalarField>>,
+    pub mu: usize,
+    pub nu: usize,
 }
 
 #[derive(Debug, Clone)]
@@ -136,6 +141,8 @@ where
     pub F: FC,
     /// public params hash
     pub pp_hash: C1::ScalarField,
+    pub mu: usize, // number of LCCCS instances to be folded
+    pub nu: usize, // number of CCCS instances to be folded
     pub i: C1::ScalarField,
     /// initial state
     pub z_0: Vec<C1::ScalarField>,
@@ -152,6 +159,180 @@ where
     pub cf_U_i: CommittedInstance<C2>,
 }
 
+impl<C1, GC1, C2, GC2, FC, CS1, CS2> MultiFolding<C1, C2, FC>
+    for HyperNova<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>,
+{
+    type RunningInstance = (LCCCS<C1>, Witness<C1::ScalarField>);
+    type IncomingInstance = (CCCS<C1>, Witness<C1::ScalarField>);
+    type MultiInstance = (Vec<Self::RunningInstance>, Vec<Self::IncomingInstance>);
+
+    /// Creates a new LCCS instance for the given state, which satisfies the HyperNova.CCS. This
+    /// method can be used to generate the 'other' LCCS instances to be folded in the multi-folding
+    /// step.
+    fn new_running_instance(
+        &self,
+        mut rng: impl RngCore,
+        state: Vec<C1::ScalarField>,
+        external_inputs: Vec<C1::ScalarField>,
+    ) -> Result<Self::RunningInstance, Error> {
+        let r1cs_z = self.new_instance_generic(state, external_inputs)?;
+        // compute committed instances, w_{i+1}, u_{i+1}, which will be used as w_i, u_i, so we
+        // assign them directly to w_i, u_i.
+        let (U_i, W_i) = self
+            .ccs
+            .to_lcccs::<_, _, CS1>(&mut rng, &self.cs_params, &r1cs_z)?;
+
+        #[cfg(test)]
+        U_i.check_relation(&self.ccs, &W_i)?;
+
+        Ok((U_i, W_i))
+    }
+
+    /// Creates a new CCCS instance for the given state, which satisfies the HyperNova.CCS. This
+    /// method can be used to generate the 'other' CCCS instances to be folded in the multi-folding
+    /// step.
+    fn new_incoming_instance(
+        &self,
+        mut rng: impl RngCore,
+        state: Vec<C1::ScalarField>,
+        external_inputs: Vec<C1::ScalarField>,
+    ) -> Result<Self::IncomingInstance, Error> {
+        let r1cs_z = self.new_instance_generic(state, external_inputs)?;
+        // compute committed instances, w_{i+1}, u_{i+1}, which will be used as w_i, u_i, so we
+        // assign them directly to w_i, u_i.
+        let (u_i, w_i) = self
+            .ccs
+            .to_cccs::<_, _, CS1>(&mut rng, &self.cs_params, &r1cs_z)?;
+
+        #[cfg(test)]
+        u_i.check_relation(&self.ccs, &w_i)?;
+
+        Ok((u_i, w_i))
+    }
+}
+
+impl<C1, GC1, C2, GC2, FC, CS1, CS2> HyperNova<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>,
+{
+    /// internal helper for new_running_instance & new_incoming_instance methods, returns the R1CS
+    /// z=[u,x,w] vector to be used to create the LCCCS & CCCS fresh instances.
+    fn new_instance_generic(
+        &self,
+        state: Vec<C1::ScalarField>,
+        external_inputs: Vec<C1::ScalarField>,
+    ) -> Result<Vec<C1::ScalarField>, Error> {
+        // 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): (NovaWitness<C2>, CommittedInstance<C2>) = self.cf_r1cs.dummy_instance();
+
+        let sponge = PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config);
+
+        u_i.x = vec![
+            U_i.hash(
+                &sponge,
+                self.pp_hash,
+                C1::ScalarField::zero(), // i
+                self.z_0.clone(),
+                state.clone(),
+            ),
+            cf_U_i.hash_cyclefold(&sponge, self.pp_hash),
+        ];
+        let us = vec![u_i.clone(); self.nu - 1];
+
+        let z_i1 = self
+            .F
+            .step_native(0, state.clone(), external_inputs.clone())?;
+
+        // compute u_{i+1}.x
+        let U_i1 = LCCCS::dummy(self.ccs.l, self.ccs.t, self.ccs.s);
+        let u_i1_x = U_i1.hash(
+            &sponge,
+            self.pp_hash,
+            C1::ScalarField::one(), // i+1, where i=0
+            self.z_0.clone(),
+            z_i1.clone(),
+        );
+
+        let cf_u_i1_x = cf_U_i.hash_cyclefold(&sponge, self.pp_hash);
+        let augmented_f_circuit = AugmentedFCircuit::<C1, C2, GC2, FC> {
+            _c2: PhantomData,
+            _gc2: PhantomData,
+            poseidon_config: self.poseidon_config.clone(),
+            ccs: self.ccs.clone(),
+            pp_hash: Some(self.pp_hash),
+            mu: self.mu,
+            nu: self.nu,
+            i: Some(C1::ScalarField::zero()),
+            i_usize: Some(0),
+            z_0: Some(self.z_0.clone()),
+            z_i: Some(state.clone()),
+            external_inputs: Some(external_inputs),
+            U_i: Some(U_i.clone()),
+            Us: None,
+            u_i_C: Some(u_i.C),
+            us: Some(us),
+            U_i1_C: Some(U_i1.C),
+            F: self.F.clone(),
+            x: Some(u_i1_x),
+            nimfs_proof: None,
+
+            // cyclefold values
+            cf_u_i_cmW: None,
+            cf_U_i: None,
+            cf_x: Some(cf_u_i1_x),
+            cf_cmT: None,
+        };
+
+        let (cs, _) = augmented_f_circuit.compute_cs_ccs()?;
+
+        #[cfg(test)]
+        assert!(cs.is_satisfied()?);
+
+        let (r1cs_w_i1, r1cs_x_i1) = extract_w_x::<C1::ScalarField>(&cs); // includes 1 and public inputs
+
+        #[cfg(test)]
+        assert_eq!(r1cs_x_i1[0], augmented_f_circuit.x.unwrap());
+
+        let r1cs_z = [
+            vec![C1::ScalarField::one()],
+            r1cs_x_i1.clone(),
+            r1cs_w_i1.clone(),
+        ]
+        .concat();
+        Ok(r1cs_z)
+    }
+}
+
 impl<C1, GC1, C2, GC2, FC, CS1, CS2> FoldingScheme<C1, C2, FC>
     for HyperNova<C1, GC1, C2, GC2, FC, CS1, CS2>
 where
@@ -170,25 +351,37 @@ where
     for<'a> &'a GC1: GroupOpsBounds<'a, C1, GC1>,
     for<'a> &'a GC2: GroupOpsBounds<'a, C2, GC2>,
 {
-    type PreprocessorParam = PreprocessorParam<C1, C2, FC, CS1, CS2>;
+    /// Reuse Nova's PreprocessorParam, together with two usize values, which are mu & nu
+    /// respectively, which indicate the amount of LCCCS & CCCS instances to be folded at each
+    /// folding step.
+    type PreprocessorParam = (PreprocessorParam<C1, C2, FC, CS1, CS2>, usize, usize);
     type ProverParam = ProverParams<C1, C2, CS1, CS2>;
     type VerifierParam = VerifierParams<C1, C2, CS1, CS2>;
     type RunningInstance = (LCCCS<C1>, Witness<C1::ScalarField>);
     type IncomingInstance = (CCCS<C1>, Witness<C1::ScalarField>);
+    type MultiCommittedInstanceWithWitness =
+        (Vec<Self::RunningInstance>, Vec<Self::IncomingInstance>);
     type CFInstance = (CommittedInstance<C2>, NovaWitness<C2>);
 
     fn preprocess(
         mut rng: impl RngCore,
         prep_param: &Self::PreprocessorParam,
     ) -> Result<(Self::ProverParam, Self::VerifierParam), Error> {
+        let (prep_param, mu, nu) = prep_param;
+        if *mu < 1 || *nu < 1 {
+            return Err(Error::CantBeZero("mu,nu".to_string()));
+        }
+
         let augmented_f_circuit = AugmentedFCircuit::<C1, C2, GC2, FC>::empty(
             &prep_param.poseidon_config,
             prep_param.F.clone(),
             None,
+            *mu,
+            *nu,
         )?;
         let ccs = augmented_f_circuit.ccs.clone();
 
-        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty();
+        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty(mu + nu);
         let cf_r1cs = get_r1cs_from_cs::<C2::ScalarField>(cf_circuit)?;
 
         // if cs params exist, use them, if not, generate new ones
@@ -215,6 +408,8 @@ where
             cs_params: cs_pp.clone(),
             cf_cs_params: cf_cs_pp.clone(),
             ccs: Some(ccs.clone()),
+            mu: *mu,
+            nu: *nu,
         };
         let vp = VerifierParams::<C1, C2, CS1, CS2> {
             poseidon_config: prep_param.poseidon_config.clone(),
@@ -228,11 +423,15 @@ where
 
     /// Initializes the HyperNova+CycleFold's IVC for the given parameters and initial state `z_0`.
     fn init(
-        params: (Self::ProverParam, Self::VerifierParam),
+        params: &(Self::ProverParam, Self::VerifierParam),
         F: FC,
         z_0: Vec<C1::ScalarField>,
     ) -> Result<Self, Error> {
         let (pp, vp) = params;
+        if pp.mu < 1 || pp.nu < 1 {
+            return Err(Error::CantBeZero("mu,nu".to_string()));
+        }
+
         // `sponge` is for digest computation.
         let sponge = PoseidonSponge::<C1::ScalarField>::new(&pp.poseidon_config);
 
@@ -242,10 +441,12 @@ where
             &pp.poseidon_config,
             F.clone(),
             pp.ccs.clone(),
+            pp.mu,
+            pp.nu,
         )?;
         let ccs = augmented_f_circuit.ccs.clone();
 
-        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty();
+        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty(pp.mu + pp.nu);
         let cf_r1cs = get_r1cs_from_cs::<C2::ScalarField>(cf_circuit)?;
 
         // compute the public params hash
@@ -282,6 +483,8 @@ where
             cf_cs_params: pp.cf_cs_params.clone(),
             F,
             pp_hash,
+            mu: pp.mu,
+            nu: pp.nu,
             i: C1::ScalarField::zero(),
             z_0: z_0.clone(),
             z_i: z_0,
@@ -300,10 +503,42 @@ where
         &mut self,
         mut rng: impl RngCore,
         external_inputs: Vec<C1::ScalarField>,
+        other_instances: Option<Self::MultiCommittedInstanceWithWitness>,
     ) -> Result<(), Error> {
         // `sponge` is for digest computation.
         let sponge = PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config);
 
+        let other_instances = other_instances.ok_or(Error::MissingOtherInstances)?;
+
+        #[allow(clippy::type_complexity)]
+        let (lcccs, cccs): (
+            Vec<(LCCCS<C1>, Witness<C1::ScalarField>)>,
+            Vec<(CCCS<C1>, Witness<C1::ScalarField>)>,
+        ) = other_instances;
+
+        // recall, mu & nu is the number of all the LCCCS & CCCS respectively, including the
+        // running and incoming instances that are not part of the 'other_instances', hence the +1
+        // in the couple of following checks.
+        if lcccs.len() + 1 != self.mu {
+            return Err(Error::NotSameLength(
+                "other_instances.lcccs.len()".to_string(),
+                lcccs.len(),
+                "hypernova.mu".to_string(),
+                self.mu,
+            ));
+        }
+        if cccs.len() + 1 != self.nu {
+            return Err(Error::NotSameLength(
+                "other_instances.cccs.len()".to_string(),
+                cccs.len(),
+                "hypernova.nu".to_string(),
+                self.nu,
+            ));
+        }
+
+        let (Us, Ws): (Vec<LCCCS<C1>>, Vec<Witness<C1::ScalarField>>) = lcccs.into_iter().unzip();
+        let (us, ws): (Vec<CCCS<C1>>, Vec<Witness<C1::ScalarField>>) = cccs.into_iter().unzip();
+
         let augmented_f_circuit: AugmentedFCircuit<C1, C2, GC2, FC>;
 
         if self.z_i.len() != self.F.state_len() {
@@ -361,13 +596,17 @@ where
                 poseidon_config: self.poseidon_config.clone(),
                 ccs: self.ccs.clone(),
                 pp_hash: Some(self.pp_hash),
+                mu: self.mu,
+                nu: self.nu,
                 i: Some(C1::ScalarField::zero()),
                 i_usize: Some(0),
                 z_0: Some(self.z_0.clone()),
                 z_i: Some(self.z_i.clone()),
                 external_inputs: Some(external_inputs.clone()),
-                u_i_C: Some(self.u_i.C),
                 U_i: Some(self.U_i.clone()),
+                Us: Some(Us.clone()),
+                u_i_C: Some(self.u_i.C),
+                us: Some(us.clone()),
                 U_i1_C: Some(U_i1.C),
                 F: self.F.clone(),
                 x: Some(u_i1_x),
@@ -383,15 +622,15 @@ where
             let mut transcript_p: PoseidonSponge<C1::ScalarField> =
                 PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config);
             transcript_p.absorb(&self.pp_hash);
-            let (rho_bits, nimfs_proof);
-            (nimfs_proof, U_i1, W_i1, rho_bits) =
+            let (rho_powers, nimfs_proof);
+            (nimfs_proof, U_i1, W_i1, rho_powers) =
                 NIMFS::<C1, PoseidonSponge<C1::ScalarField>>::prove(
                     &mut transcript_p,
                     &self.ccs,
-                    &[self.U_i.clone()],
-                    &[self.u_i.clone()],
-                    &[self.W_i.clone()],
-                    &[self.w_i.clone()],
+                    &[vec![self.U_i.clone()], Us.clone()].concat(),
+                    &[vec![self.u_i.clone()], us.clone()].concat(),
+                    &[vec![self.W_i.clone()], Ws].concat(),
+                    &[vec![self.w_i.clone()], ws].concat(),
                 )?;
 
             // sanity check: check the folded instance relation
@@ -406,29 +645,60 @@ where
                 z_i1.clone(),
             );
 
-            let rho_Fq = C2::ScalarField::from_bigint(BigInteger::from_bits_le(&rho_bits))
-                .ok_or(Error::OutOfBounds)?;
+            let rho_powers_Fq: Vec<C1::BaseField> = rho_powers
+                .iter()
+                .map(|rho_i| {
+                    C1::BaseField::from_bigint(BigInteger::from_bits_le(
+                        &rho_i.into_bigint().to_bits_le(),
+                    ))
+                    .unwrap()
+                })
+                .collect();
+            let rho_powers_bits: Vec<Vec<bool>> = rho_powers
+                .iter()
+                .map(|rho_i| rho_i.into_bigint().to_bits_le()[..N_BITS_RO].to_vec())
+                .collect();
+
             // CycleFold part:
-            // get the vector used as public inputs 'x' in the CycleFold circuit
-            // cyclefold circuit for cmW
+            // get the vector used as public inputs 'x' in the CycleFold circuit.
+            // Place the random values and the points coordinates as the public input x:
+            // In Nova, this is: x == [r, p1, p2, p3].
+            // In multifolding schemes such as HyperNova, this is:
+            // computed_x = [r_0, r_1, r_2, ...,  r_n, p_0, p_1, p_2, ..., p_n],
+            // where each p_i is in fact p_i.to_constraint_field()
             let cf_u_i_x = [
-                vec![rho_Fq],
+                rho_powers_Fq,
                 get_cm_coordinates(&self.U_i.C),
+                Us.iter()
+                    .flat_map(|Us_i| get_cm_coordinates(&Us_i.C))
+                    .collect(),
                 get_cm_coordinates(&self.u_i.C),
+                us.iter()
+                    .flat_map(|us_i| get_cm_coordinates(&us_i.C))
+                    .collect(),
                 get_cm_coordinates(&U_i1.C),
             ]
             .concat();
 
             let cf_circuit = CycleFoldCircuit::<C1, GC1> {
                 _gc: PhantomData,
-                r_bits: Some(rho_bits.clone()),
-                p1: Some(self.U_i.clone().C),
-                p2: Some(self.u_i.clone().C),
+                n_points: self.mu + self.nu,
+                r_bits: Some(rho_powers_bits.clone()),
+                points: Some(
+                    [
+                        vec![self.U_i.clone().C],
+                        Us.iter().map(|Us_i| Us_i.C).collect(),
+                        vec![self.u_i.clone().C],
+                        us.iter().map(|us_i| us_i.C).collect(),
+                    ]
+                    .concat(),
+                ),
                 x: Some(cf_u_i_x.clone()),
             };
 
             let (_cf_w_i, cf_u_i, cf_W_i1, cf_U_i1, cf_cmT, _) =
                 fold_cyclefold_circuit::<C1, GC1, C2, GC2, FC, CS1, CS2>(
+                    self.mu + self.nu,
                     &mut transcript_p,
                     self.cf_r1cs.clone(),
                     self.cf_cs_params.clone(),
@@ -447,13 +717,17 @@ where
                 poseidon_config: self.poseidon_config.clone(),
                 ccs: self.ccs.clone(),
                 pp_hash: Some(self.pp_hash),
+                mu: self.mu,
+                nu: self.nu,
                 i: Some(self.i),
                 i_usize: Some(i_usize),
                 z_0: Some(self.z_0.clone()),
                 z_i: Some(self.z_i.clone()),
                 external_inputs: Some(external_inputs),
-                u_i_C: Some(self.u_i.C),
                 U_i: Some(self.U_i.clone()),
+                Us: Some(Us.clone()),
+                u_i_C: Some(self.u_i.C),
+                us: Some(us.clone()),
                 U_i1_C: Some(U_i1.C),
                 F: self.F.clone(),
                 x: Some(u_i1_x),
@@ -618,31 +892,48 @@ mod tests {
 
         type HN<CS1, CS2> =
             HyperNova<Projective, GVar, Projective2, GVar2, CubicFCircuit<Fr>, CS1, CS2>;
+        let (mu, nu) = (2, 3);
 
         let prep_param =
             PreprocessorParam::<Projective, Projective2, CubicFCircuit<Fr>, CS1, CS2>::new(
                 poseidon_config.clone(),
                 F_circuit,
             );
-        let (prover_params, verifier_params) = HN::preprocess(&mut rng, &prep_param).unwrap();
+        let hypernova_params = HN::preprocess(&mut rng, &(prep_param, mu, nu)).unwrap();
 
         let z_0 = vec![Fr::from(3_u32)];
-        let mut hypernova = HN::init(
-            (prover_params, verifier_params.clone()),
-            F_circuit,
-            z_0.clone(),
-        )
-        .unwrap();
+        let mut hypernova = HN::init(&hypernova_params, F_circuit, z_0.clone()).unwrap();
 
         let num_steps: usize = 3;
         for _ in 0..num_steps {
-            hypernova.prove_step(&mut rng, vec![]).unwrap();
+            // prepare some new instances to fold in the multifolding step
+            let mut lcccs = vec![];
+            for j in 0..mu - 1 {
+                let instance_state = vec![Fr::from(j as u32 + 85_u32)];
+                let (U, W) = hypernova
+                    .new_running_instance(&mut rng, instance_state, vec![])
+                    .unwrap();
+                lcccs.push((U, W));
+            }
+            let mut cccs = vec![];
+            for j in 0..nu - 1 {
+                let instance_state = vec![Fr::from(j as u32 + 15_u32)];
+                let (u, w) = hypernova
+                    .new_incoming_instance(&mut rng, instance_state, vec![])
+                    .unwrap();
+                cccs.push((u, w));
+            }
+
+            dbg!(&hypernova.i);
+            hypernova
+                .prove_step(&mut rng, vec![], Some((lcccs, cccs)))
+                .unwrap();
         }
         assert_eq!(Fr::from(num_steps as u32), hypernova.i);
 
         let (running_instance, incoming_instance, cyclefold_instance) = hypernova.instances();
         HN::verify(
-            verifier_params,
+            hypernova_params.1, // verifier_params
             z_0,
             hypernova.z_i,
             hypernova.i,
diff --git a/folding-schemes/src/folding/hypernova/nimfs.rs b/folding-schemes/src/folding/hypernova/nimfs.rs
index ab3bb52..2a010a6 100644
--- a/folding-schemes/src/folding/hypernova/nimfs.rs
+++ b/folding-schemes/src/folding/hypernova/nimfs.rs
@@ -73,7 +73,7 @@ where
         sigmas_thetas: &SigmasThetas<C::ScalarField>,
         r_x_prime: Vec<C::ScalarField>,
         rho: C::ScalarField,
-    ) -> LCCCS<C> {
+    ) -> (LCCCS<C>, Vec<C::ScalarField>) {
         let (sigmas, thetas) = (sigmas_thetas.0.clone(), sigmas_thetas.1.clone());
         let mut C_folded = C::zero();
         let mut u_folded = C::ScalarField::zero();
@@ -81,6 +81,7 @@ where
         let mut v_folded: Vec<C::ScalarField> = vec![C::ScalarField::zero(); sigmas[0].len()];
 
         let mut rho_i = C::ScalarField::one();
+        let mut rho_powers = vec![C::ScalarField::zero(); lcccs.len() + cccs.len() - 1];
         for i in 0..(lcccs.len() + cccs.len()) {
             let c: C;
             let u: C::ScalarField;
@@ -120,16 +121,28 @@ where
                 .map(|(a_i, b_i)| *a_i + b_i)
                 .collect();
 
+            // compute the next power of rho
             rho_i *= rho;
+            // crop the size of rho_i to N_BITS_RO
+            let rho_i_bits = rho_i.into_bigint().to_bits_le();
+            rho_i = C::ScalarField::from_bigint(BigInteger::from_bits_le(&rho_i_bits[..N_BITS_RO]))
+                .unwrap();
+            if i < lcccs.len() + cccs.len() - 1 {
+                // store the cropped rho_i into the rho_powers vector
+                rho_powers[i] = rho_i;
+            }
         }
 
-        LCCCS::<C> {
-            C: C_folded,
-            u: u_folded,
-            x: x_folded,
-            r_x: r_x_prime,
-            v: v_folded,
-        }
+        (
+            LCCCS::<C> {
+                C: C_folded,
+                u: u_folded,
+                x: x_folded,
+                r_x: r_x_prime,
+                v: v_folded,
+            },
+            rho_powers,
+        )
     }
 
     pub fn fold_witness(
@@ -140,8 +153,9 @@ where
         let mut w_folded: Vec<C::ScalarField> = vec![C::ScalarField::zero(); w_lcccs[0].w.len()];
         let mut r_w_folded = C::ScalarField::zero();
 
+        let mut rho_i = C::ScalarField::one();
         for i in 0..(w_lcccs.len() + w_cccs.len()) {
-            let rho_i = rho.pow([i as u64]);
+            // let rho_i = rho.pow([i as u64]);
             let w: Vec<C::ScalarField>;
             let r_w: C::ScalarField;
 
@@ -164,6 +178,13 @@ where
                 .collect();
 
             r_w_folded += rho_i * r_w;
+
+            // compute the next power of rho
+            rho_i *= rho;
+            // crop the size of rho_i to N_BITS_RO
+            let rho_i_bits = rho_i.into_bigint().to_bits_le();
+            rho_i = C::ScalarField::from_bigint(BigInteger::from_bits_le(&rho_i_bits[..N_BITS_RO]))
+                .unwrap();
         }
         Witness {
             w: w_folded,
@@ -183,7 +204,16 @@ where
         new_instances: &[CCCS<C>],
         w_lcccs: &[Witness<C::ScalarField>],
         w_cccs: &[Witness<C::ScalarField>],
-    ) -> Result<(NIMFSProof<C>, LCCCS<C>, Witness<C::ScalarField>, Vec<bool>), Error> {
+    ) -> Result<
+        (
+            NIMFSProof<C>,
+            LCCCS<C>,
+            Witness<C::ScalarField>,
+            // Vec<bool>,
+            Vec<C::ScalarField>,
+        ),
+        Error,
+    > {
         // absorb instances to transcript
         transcript.absorb(&running_instances);
         transcript.absorb(&new_instances);
@@ -247,7 +277,7 @@ where
             C::ScalarField::from_bigint(BigInteger::from_bits_le(&rho_bits)).unwrap();
 
         // Step 7: Create the folded instance
-        let folded_lcccs = Self::fold(
+        let (folded_lcccs, rho_powers) = Self::fold(
             running_instances,
             new_instances,
             &sigmas_thetas,
@@ -265,7 +295,7 @@ where
             },
             folded_lcccs,
             folded_witness,
-            rho_bits,
+            rho_powers,
         ))
     }
 
@@ -372,7 +402,8 @@ where
             &proof.sigmas_thetas,
             r_x_prime,
             rho,
-        ))
+        )
+        .0)
     }
 }
 
@@ -422,7 +453,7 @@ pub mod tests {
         let mut rng = test_rng();
         let rho = Fr::rand(&mut rng);
 
-        let folded = NIMFS::<Projective, PoseidonSponge<Fr>>::fold(
+        let (folded, _) = NIMFS::<Projective, PoseidonSponge<Fr>>::fold(
             &[lcccs],
             &[cccs],
             &sigmas_thetas,
diff --git a/folding-schemes/src/folding/nova/circuits.rs b/folding-schemes/src/folding/nova/circuits.rs
index 3ad2269..49532f3 100644
--- a/folding-schemes/src/folding/nova/circuits.rs
+++ b/folding-schemes/src/folding/nova/circuits.rs
@@ -20,11 +20,11 @@ use ark_relations::r1cs::{ConstraintSynthesizer, ConstraintSystemRef, Namespace,
 use ark_std::{fmt::Debug, One, Zero};
 use core::{borrow::Borrow, marker::PhantomData};
 
-use super::CommittedInstance;
+use super::{CommittedInstance, NOVA_CF_N_POINTS};
 use crate::constants::N_BITS_RO;
 use crate::folding::circuits::{
     cyclefold::{
-        CycleFoldChallengeGadget, CycleFoldCommittedInstanceVar, NIFSFullGadget, CF_IO_LEN,
+        cf_io_len, CycleFoldChallengeGadget, CycleFoldCommittedInstanceVar, NIFSFullGadget,
     },
     nonnative::{affine::NonNativeAffineVar, uint::NonNativeUintVar},
     CF1, CF2,
@@ -337,7 +337,7 @@ where
         let cmT =
             NonNativeAffineVar::new_witness(cs.clone(), || Ok(self.cmT.unwrap_or_else(C1::zero)))?;
 
-        let cf_u_dummy = CommittedInstance::dummy(CF_IO_LEN);
+        let cf_u_dummy = CommittedInstance::dummy(cf_io_len(NOVA_CF_N_POINTS));
         let cf_U_i = CycleFoldCommittedInstanceVar::<C2, GC2>::new_witness(cs.clone(), || {
             Ok(self.cf_U_i.unwrap_or(cf_u_dummy.clone()))
         })?;
diff --git a/folding-schemes/src/folding/nova/decider_eth.rs b/folding-schemes/src/folding/nova/decider_eth.rs
index 4b011e8..a8a0136 100644
--- a/folding-schemes/src/folding/nova/decider_eth.rs
+++ b/folding-schemes/src/folding/nova/decider_eth.rs
@@ -361,12 +361,12 @@ pub mod tests {
         let nova_params = N::preprocess(&mut rng, &prep_param).unwrap();
 
         let start = Instant::now();
-        let mut nova = N::init(nova_params.clone(), F_circuit, z_0.clone()).unwrap();
+        let mut nova = N::init(&nova_params, F_circuit, z_0.clone()).unwrap();
         println!("Nova initialized, {:?}", start.elapsed());
         let start = Instant::now();
-        nova.prove_step(&mut rng, vec![]).unwrap();
+        nova.prove_step(&mut rng, vec![], None).unwrap();
         println!("prove_step, {:?}", start.elapsed());
-        nova.prove_step(&mut rng, vec![]).unwrap(); // do a 2nd step
+        nova.prove_step(&mut rng, vec![], None).unwrap(); // do a 2nd step
 
         let mut rng = rand::rngs::OsRng;
 
diff --git a/folding-schemes/src/folding/nova/decider_eth_circuit.rs b/folding-schemes/src/folding/nova/decider_eth_circuit.rs
index 2a0ab0b..4cd6c65 100644
--- a/folding-schemes/src/folding/nova/decider_eth_circuit.rs
+++ b/folding-schemes/src/folding/nova/decider_eth_circuit.rs
@@ -447,11 +447,12 @@ where
         {
             // imports here instead of at the top of the file, so we avoid having multiple
             // `#[cfg(not(test))]`
+            use super::NOVA_CF_N_POINTS;
             use crate::commitment::pedersen::PedersenGadget;
-            use crate::folding::circuits::cyclefold::{CycleFoldCommittedInstanceVar, CF_IO_LEN};
+            use crate::folding::circuits::cyclefold::{cf_io_len, CycleFoldCommittedInstanceVar};
             use ark_r1cs_std::ToBitsGadget;
 
-            let cf_u_dummy_native = CommittedInstance::<C2>::dummy(CF_IO_LEN);
+            let cf_u_dummy_native = CommittedInstance::<C2>::dummy(cf_io_len(NOVA_CF_N_POINTS));
             let w_dummy_native = Witness::<C2>::new(
                 vec![C2::ScalarField::zero(); self.cf_r1cs.A.n_cols - 1 - self.cf_r1cs.l],
                 self.cf_E_len,
@@ -796,20 +797,15 @@ pub mod tests {
             Pedersen<Projective>,
             Pedersen<Projective2>,
         >::new(poseidon_config, F_circuit);
-        let (prover_params, verifier_params) = N::preprocess(&mut rng, &prep_param).unwrap();
+        let nova_params = N::preprocess(&mut rng, &prep_param).unwrap();
 
         // generate a Nova instance and do a step of it
-        let mut nova = N::init(
-            (prover_params, verifier_params.clone()),
-            F_circuit,
-            z_0.clone(),
-        )
-        .unwrap();
-        nova.prove_step(&mut rng, vec![]).unwrap();
+        let mut nova = N::init(&nova_params, F_circuit, z_0.clone()).unwrap();
+        nova.prove_step(&mut rng, vec![], None).unwrap();
         let ivc_v = nova.clone();
         let (running_instance, incoming_instance, cyclefold_instance) = ivc_v.instances();
         N::verify(
-            verifier_params,
+            nova_params.1, // verifier_params
             z_0,
             ivc_v.z_i,
             Fr::one(),
diff --git a/folding-schemes/src/folding/nova/mod.rs b/folding-schemes/src/folding/nova/mod.rs
index 4847314..99d559d 100644
--- a/folding-schemes/src/folding/nova/mod.rs
+++ b/folding-schemes/src/folding/nova/mod.rs
@@ -37,6 +37,10 @@ use circuits::{AugmentedFCircuit, ChallengeGadget};
 use nifs::NIFS;
 use traits::NovaR1CS;
 
+/// Number of points to be folded in the CycleFold circuit, in Nova's case, this is a fixed amount:
+/// 2 points to be folded.
+const NOVA_CF_N_POINTS: usize = 2_usize;
+
 #[derive(Debug, Clone, Eq, PartialEq, CanonicalSerialize, CanonicalDeserialize)]
 pub struct CommittedInstance<C: CurveGroup> {
     pub cmE: C,
@@ -344,6 +348,7 @@ where
     type VerifierParam = VerifierParams<C1, C2, CS1, CS2>;
     type RunningInstance = (CommittedInstance<C1>, Witness<C1>);
     type IncomingInstance = (CommittedInstance<C1>, Witness<C1>);
+    type MultiCommittedInstanceWithWitness = ();
     type CFInstance = (CommittedInstance<C2>, Witness<C2>);
 
     fn preprocess(
@@ -390,7 +395,7 @@ where
 
     /// Initializes the Nova+CycleFold's IVC for the given parameters and initial state `z_0`.
     fn init(
-        params: (Self::ProverParam, Self::VerifierParam),
+        params: &(Self::ProverParam, Self::VerifierParam),
         F: FC,
         z_0: Vec<C1::ScalarField>,
     ) -> Result<Self, Error> {
@@ -402,7 +407,7 @@ where
 
         let augmented_F_circuit =
             AugmentedFCircuit::<C1, C2, GC2, FC>::empty(&pp.poseidon_config, F.clone());
-        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty();
+        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty(NOVA_CF_N_POINTS);
 
         augmented_F_circuit.generate_constraints(cs.clone())?;
         cs.finalize();
@@ -452,6 +457,8 @@ where
         &mut self,
         _rng: impl RngCore,
         external_inputs: Vec<C1::ScalarField>,
+        // Nova does not support multi-instances folding
+        _other_instances: Option<Self::MultiCommittedInstanceWithWitness>,
     ) -> Result<(), Error> {
         // `sponge` is for digest computation.
         let sponge = PoseidonSponge::<C1::ScalarField>::new(&self.poseidon_config);
@@ -460,6 +467,11 @@ where
 
         let augmented_F_circuit: AugmentedFCircuit<C1, C2, GC2, FC>;
 
+        // Nova does not support (by design) multi-instances folding
+        if _other_instances.is_some() {
+            return Err(Error::NoMultiInstances);
+        }
+
         if self.z_i.len() != self.F.state_len() {
             return Err(Error::NotSameLength(
                 "z_i.len()".to_string(),
@@ -572,16 +584,16 @@ where
 
             let cfW_circuit = CycleFoldCircuit::<C1, GC1> {
                 _gc: PhantomData,
-                r_bits: Some(r_bits.clone()),
-                p1: Some(self.U_i.clone().cmW),
-                p2: Some(self.u_i.clone().cmW),
+                n_points: NOVA_CF_N_POINTS,
+                r_bits: Some(vec![r_bits.clone()]),
+                points: Some(vec![self.U_i.clone().cmW, self.u_i.clone().cmW]),
                 x: Some(cfW_u_i_x.clone()),
             };
             let cfE_circuit = CycleFoldCircuit::<C1, GC1> {
                 _gc: PhantomData,
-                r_bits: Some(r_bits.clone()),
-                p1: Some(self.U_i.clone().cmE),
-                p2: Some(cmT),
+                n_points: NOVA_CF_N_POINTS,
+                r_bits: Some(vec![r_bits.clone()]),
+                points: Some(vec![self.U_i.clone().cmE, cmT]),
                 x: Some(cfE_u_i_x.clone()),
             };
 
@@ -820,6 +832,7 @@ where
         Error,
     > {
         fold_cyclefold_circuit::<C1, GC1, C2, GC2, FC, CS1, CS2>(
+            NOVA_CF_N_POINTS,
             transcript,
             self.cf_r1cs.clone(),
             self.cf_cs_pp.clone(),
@@ -866,7 +879,7 @@ where
 {
     let augmented_F_circuit =
         AugmentedFCircuit::<C1, C2, GC2, FC>::empty(poseidon_config, F_circuit);
-    let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty();
+    let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty(NOVA_CF_N_POINTS);
     let r1cs = get_r1cs_from_cs::<C1::ScalarField>(augmented_F_circuit)?;
     let cf_r1cs = get_r1cs_from_cs::<C2::ScalarField>(cf_circuit)?;
     Ok((r1cs, cf_r1cs))
@@ -944,11 +957,11 @@ pub mod tests {
         let nova_params = N::preprocess(&mut rng, &prep_param).unwrap();
 
         let z_0 = vec![Fr::from(3_u32)];
-        let mut nova = N::init(nova_params.clone(), F_circuit, z_0.clone()).unwrap();
+        let mut nova = N::init(&nova_params, F_circuit, z_0.clone()).unwrap();
 
         let num_steps: usize = 3;
         for _ in 0..num_steps {
-            nova.prove_step(&mut rng, vec![]).unwrap();
+            nova.prove_step(&mut rng, vec![], None).unwrap();
         }
         assert_eq!(Fr::from(num_steps as u32), nova.i);
 
diff --git a/folding-schemes/src/folding/nova/serialize.rs b/folding-schemes/src/folding/nova/serialize.rs
index 6506adb..62dea3f 100644
--- a/folding-schemes/src/folding/nova/serialize.rs
+++ b/folding-schemes/src/folding/nova/serialize.rs
@@ -12,7 +12,10 @@ use std::marker::PhantomData;
 
 use super::{circuits::AugmentedFCircuit, Nova, ProverParams};
 use super::{CommittedInstance, Witness};
-use crate::folding::circuits::{cyclefold::CycleFoldCircuit, CF2};
+use crate::folding::{
+    circuits::{cyclefold::CycleFoldCircuit, CF2},
+    nova::NOVA_CF_N_POINTS,
+};
 use crate::{
     arith::r1cs::extract_r1cs, commitment::CommitmentScheme, folding::circuits::CF1,
     frontend::FCircuit,
@@ -134,7 +137,7 @@ where
         let cs2 = ConstraintSystem::<C1::BaseField>::new_ref();
         let augmented_F_circuit =
             AugmentedFCircuit::<C1, C2, GC2, FC>::empty(&poseidon_config, f_circuit.clone());
-        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty();
+        let cf_circuit = CycleFoldCircuit::<C1, GC1>::empty(NOVA_CF_N_POINTS);
 
         augmented_F_circuit
             .generate_constraints(cs.clone())
@@ -209,11 +212,11 @@ pub mod tests {
         let nova_params = N::preprocess(&mut rng, &prep_param).unwrap();
 
         let z_0 = vec![Fr::from(3_u32)];
-        let mut nova = N::init(nova_params.clone(), F_circuit, z_0.clone()).unwrap();
+        let mut nova = N::init(&nova_params, F_circuit, z_0.clone()).unwrap();
 
         let num_steps: usize = 3;
         for _ in 0..num_steps {
-            nova.prove_step(&mut rng, vec![]).unwrap();
+            nova.prove_step(&mut rng, vec![], None).unwrap();
         }
 
         let mut writer = vec![];
@@ -252,8 +255,10 @@ pub mod tests {
 
         let num_steps: usize = 3;
         for _ in 0..num_steps {
-            deserialized_nova.prove_step(&mut rng, vec![]).unwrap();
-            nova.prove_step(&mut rng, vec![]).unwrap();
+            deserialized_nova
+                .prove_step(&mut rng, vec![], None)
+                .unwrap();
+            nova.prove_step(&mut rng, vec![], None).unwrap();
         }
 
         assert_eq!(deserialized_nova.w_i, nova.w_i);
diff --git a/folding-schemes/src/lib.rs b/folding-schemes/src/lib.rs
index 4ed0176..73b30e8 100644
--- a/folding-schemes/src/lib.rs
+++ b/folding-schemes/src/lib.rs
@@ -69,6 +69,8 @@ pub enum Error {
     NotEnoughSteps,
     #[error("Evaluation failed")]
     EvaluationFail,
+    #[error("{0} can not be zero")]
+    CantBeZero(String),
 
     // Commitment errors
     #[error("Pedersen parameters length is not sufficient (generators.len={0} < vector.len={1} unsatisfied)")]
@@ -97,6 +99,10 @@ pub enum Error {
     BigIntConversionError(String),
     #[error("Failed to serde: {0}")]
     JSONSerdeError(String),
+    #[error("Multi instances folding not supported in this scheme")]
+    NoMultiInstances,
+    #[error("Missing 'other' instances, since this is a multi-instances folding scheme")]
+    MissingOtherInstances,
 }
 
 /// FoldingScheme defines trait that is implemented by the diverse folding schemes. It is defined
@@ -116,6 +122,7 @@ where
     type VerifierParam: Debug + Clone;
     type RunningInstance: Debug; // contains the CommittedInstance + Witness
     type IncomingInstance: Debug; // contains the CommittedInstance + Witness
+    type MultiCommittedInstanceWithWitness: Debug; // type used for the extra instances in the multi-instance folding setting
     type CFInstance: Debug; // CycleFold CommittedInstance & Witness
 
     fn preprocess(
@@ -124,7 +131,7 @@ where
     ) -> Result<(Self::ProverParam, Self::VerifierParam), Error>;
 
     fn init(
-        params: (Self::ProverParam, Self::VerifierParam),
+        params: &(Self::ProverParam, Self::VerifierParam),
         step_circuit: FC,
         z_0: Vec<C1::ScalarField>, // initial state
     ) -> Result<Self, Error>;
@@ -133,6 +140,7 @@ where
         &mut self,
         rng: impl RngCore,
         external_inputs: Vec<C1::ScalarField>,
+        other_instances: Option<Self::MultiCommittedInstanceWithWitness>,
     ) -> Result<(), Error>;
 
     // returns the state at the current step
@@ -160,6 +168,35 @@ where
     ) -> Result<(), Error>;
 }
 
+/// Trait with auxiliary methods for multi-folding schemes (ie. HyperNova, ProtoGalaxy, etc),
+/// allowing to create new instances for the multifold.
+pub trait MultiFolding<C1: CurveGroup, C2: CurveGroup, FC>: Clone + Debug
+where
+    C1: CurveGroup<BaseField = C2::ScalarField, ScalarField = C2::BaseField>,
+    C2::BaseField: PrimeField,
+    FC: FCircuit<C1::ScalarField>,
+{
+    type RunningInstance: Debug;
+    type IncomingInstance: Debug;
+    type MultiInstance: Debug;
+
+    /// Creates a new RunningInstance for the given state, to be folded in the multi-folding step.
+    fn new_running_instance(
+        &self,
+        rng: impl RngCore,
+        state: Vec<C1::ScalarField>,
+        external_inputs: Vec<C1::ScalarField>,
+    ) -> Result<Self::RunningInstance, Error>;
+
+    /// Creates a new IncomingInstance for the given state, to be folded in the multi-folding step.
+    fn new_incoming_instance(
+        &self,
+        rng: impl RngCore,
+        state: Vec<C1::ScalarField>,
+        external_inputs: Vec<C1::ScalarField>,
+    ) -> Result<Self::IncomingInstance, Error>;
+}
+
 pub trait Decider<
     C1: CurveGroup,
     C2: CurveGroup,
diff --git a/folding-schemes/src/utils/mod.rs b/folding-schemes/src/utils/mod.rs
index 764dc18..ca017d5 100644
--- a/folding-schemes/src/utils/mod.rs
+++ b/folding-schemes/src/utils/mod.rs
@@ -60,9 +60,9 @@ where
     hasher.update(C1::ScalarField::MODULUS_BIT_SIZE.to_le_bytes());
     hasher.update(C2::ScalarField::MODULUS_BIT_SIZE.to_le_bytes());
     // AugmentedFCircuit Arith params
-    hasher.update(arith.params_to_bytes());
+    hasher.update(arith.params_to_le_bytes());
     // CycleFold Circuit Arith params
-    hasher.update(cf_arith.params_to_bytes());
+    hasher.update(cf_arith.params_to_le_bytes());
     // cs_vp & cf_cs_vp (commitments setup)
     let mut cs_vp_bytes = Vec::new();
     cs_vp.serialize_uncompressed(&mut cs_vp_bytes)?;
diff --git a/solidity-verifiers/src/verifiers/nova_cyclefold.rs b/solidity-verifiers/src/verifiers/nova_cyclefold.rs
index 9328460..afe7686 100644
--- a/solidity-verifiers/src/verifiers/nova_cyclefold.rs
+++ b/solidity-verifiers/src/verifiers/nova_cyclefold.rs
@@ -324,7 +324,7 @@ mod tests {
         );
         let nova_params = NOVA::preprocess(&mut rng, &prep_param).unwrap();
         let nova = NOVA::init(
-            nova_params.clone(),
+            &nova_params,
             f_circuit.clone(),
             vec![Fr::zero(); f_circuit.state_len()].clone(),
         )
@@ -358,9 +358,9 @@ mod tests {
 
         let mut rng = rand::rngs::OsRng;
 
-        let mut nova = NOVA::<FC>::init(fs_params, f_circuit, z_0).unwrap();
+        let mut nova = NOVA::<FC>::init(&fs_params, f_circuit, z_0).unwrap();
         for _ in 0..n_steps {
-            nova.prove_step(&mut rng, vec![]).unwrap();
+            nova.prove_step(&mut rng, vec![], None).unwrap();
         }
 
         let start = Instant::now();