MNT4/6 curves and recursive SNARKs (#150)

* Add mnt6_753 curve
Generalize mnt6 curve model

* Add mnt4 curves

* Use resampled generators

* Calculate correct G2 cofactors

* Add fields to r1cs-std

* Add pairings

* Improve reusing of Fq/Fr among MNT curves

* Add instantiations of curves
Fix Fp6_2over3
Rebase code to current master

* Add test for recursive NIZK proof verification

* Address comments in PR

* Improve test case and port to GM17
Also fix a minor bug in to_field_vec

crypto-primitives/Cargo.toml

@@ -46,6 +46,6 @@ std = ["r1cs", "algebra-core/std", "r1cs-core/std", "r1cs-std/std"]
 parallel = ["std", "rayon", "gm17/parallel", "groth16/parallel", "ff-fft/parallel"]
 
 [dev-dependencies]
-algebra = { path = "../algebra", default-features = false, features = [ "jubjub", "bls12_377" ] }
-r1cs-std = { path = "../r1cs-std", default-features = false, features = [ "jubjub", "bls12_377" ] }
+algebra = { path = "../algebra", default-features = false, features = [ "jubjub", "bls12_377", "mnt4_298", "mnt6_298" ] }
+r1cs-std = { path = "../r1cs-std", default-features = false, features = [ "jubjub", "bls12_377", "mnt4_298", "mnt6_298" ] }
 rand_xorshift = { version = "0.2" }

crypto-primitives/src/nizk/gm17/constraints.rs

@@ -532,3 +532,266 @@ mod test {
         }
     }
 }
+
+#[cfg(test)]
+mod test_recursive {
+    use gm17::*;
+    use r1cs_core::{ConstraintSynthesizer, ConstraintSystem, SynthesisError};
+
+    use super::*;
+    use algebra::{
+        fields::FpParameters,
+        mnt4_298::{Fq as MNT4Fq, FqParameters as MNT4FqParameters, Fr as MNT4Fr, MNT4_298},
+        mnt6_298::{Fq as MNT6Fq, FqParameters as MNT6FqParameters, Fr as MNT6Fr, MNT6_298},
+        test_rng, BigInteger, PrimeField,
+    };
+    use r1cs_std::{
+        fields::fp::FpGadget, mnt4_298::PairingGadget as MNT4_298PairingGadget,
+        mnt6_298::PairingGadget as MNT6_298PairingGadget,
+        test_constraint_system::TestConstraintSystem, uint8::UInt8,
+    };
+    use rand::Rng;
+
+    type TestProofSystem1 = Gm17<MNT6_298, Bench<MNT4Fq>, MNT6Fr>;
+    type TestVerifierGadget1 = Gm17VerifierGadget<MNT6_298, MNT6Fq, MNT6_298PairingGadget>;
+    type TestProofGadget1 = ProofGadget<MNT6_298, MNT6Fq, MNT6_298PairingGadget>;
+    type TestVkGadget1 = VerifyingKeyGadget<MNT6_298, MNT6Fq, MNT6_298PairingGadget>;
+
+    type TestProofSystem2 = Gm17<MNT4_298, Wrapper, MNT4Fr>;
+    type TestVerifierGadget2 = Gm17VerifierGadget<MNT4_298, MNT4Fq, MNT4_298PairingGadget>;
+    type TestProofGadget2 = ProofGadget<MNT4_298, MNT4Fq, MNT4_298PairingGadget>;
+    type TestVkGadget2 = VerifyingKeyGadget<MNT4_298, MNT4Fq, MNT4_298PairingGadget>;
+
+    #[derive(Clone)]
+    struct Bench<F: Field> {
+        inputs:          Vec<Option<F>>,
+        num_constraints: usize,
+    }
+
+    impl<F: Field> ConstraintSynthesizer<F> for Bench<F> {
+        fn generate_constraints<CS: ConstraintSystem<F>>(
+            self,
+            cs: &mut CS,
+        ) -> Result<(), SynthesisError> {
+            assert!(self.inputs.len() >= 2);
+            assert!(self.num_constraints >= self.inputs.len());
+
+            let mut variables: Vec<_> = Vec::with_capacity(self.inputs.len());
+            for (i, input) in self.inputs.into_iter().enumerate() {
+                let input_var = cs.alloc_input(
+                    || format!("Input {}", i),
+                    || input.ok_or(SynthesisError::AssignmentMissing),
+                )?;
+                variables.push((input, input_var));
+            }
+
+            for i in 0..self.num_constraints {
+                let new_entry = {
+                    let (input_1_val, input_1_var) = variables[i];
+                    let (input_2_val, input_2_var) = variables[i + 1];
+                    let result_val = input_1_val
+                        .and_then(|input_1| input_2_val.map(|input_2| input_1 * &input_2));
+                    let result_var = cs.alloc(
+                        || format!("Result {}", i),
+                        || result_val.ok_or(SynthesisError::AssignmentMissing),
+                    )?;
+                    cs.enforce(
+                        || format!("Enforce constraint {}", i),
+                        |lc| lc + input_1_var,
+                        |lc| lc + input_2_var,
+                        |lc| lc + result_var,
+                    );
+                    (result_val, result_var)
+                };
+                variables.push(new_entry);
+            }
+            Ok(())
+        }
+    }
+
+    struct Wrapper {
+        inputs: Vec<Option<MNT4Fq>>,
+        params: Parameters<MNT6_298>,
+        proof:  Proof<MNT6_298>,
+    }
+
+    impl ConstraintSynthesizer<MNT6Fq> for Wrapper {
+        fn generate_constraints<CS: ConstraintSystem<MNT6Fq>>(
+            self,
+            cs: &mut CS,
+        ) -> Result<(), SynthesisError> {
+            let params = self.params;
+            let proof = self.proof;
+            let inputs: Vec<_> = self
+                .inputs
+                .into_iter()
+                .map(|input| input.unwrap())
+                .collect();
+            let input_gadgets;
+
+            {
+                let mut cs = cs.ns(|| "Allocate Input");
+                // Chain all input values in one large byte array.
+                let input_bytes = inputs
+                    .clone()
+                    .into_iter()
+                    .flat_map(|input| {
+                        input
+                            .into_repr()
+                            .as_ref()
+                            .iter()
+                            .flat_map(|l| l.to_le_bytes().to_vec())
+                            .collect::<Vec<_>>()
+                    })
+                    .collect::<Vec<_>>();
+
+                // Allocate this byte array as input packed into field elements.
+                let input_bytes = UInt8::alloc_input_vec(cs.ns(|| "Input"), &input_bytes[..])?;
+                // 40 byte
+                let element_size = <MNT4FqParameters as FpParameters>::BigInt::NUM_LIMBS * 8;
+                input_gadgets = input_bytes
+                    .chunks(element_size)
+                    .map(|chunk| {
+                        chunk
+                            .iter()
+                            .flat_map(|byte| byte.into_bits_le())
+                            .collect::<Vec<_>>()
+                    })
+                    .collect::<Vec<_>>();
+            }
+
+            let vk_gadget = TestVkGadget1::alloc(cs.ns(|| "Vk"), || Ok(&params.vk))?;
+            let proof_gadget =
+                TestProofGadget1::alloc(cs.ns(|| "Proof"), || Ok(proof.clone())).unwrap();
+            <TestVerifierGadget1 as NIZKVerifierGadget<TestProofSystem1, MNT6Fq>>::check_verify(
+                cs.ns(|| "Verify"),
+                &vk_gadget,
+                input_gadgets.iter(),
+                &proof_gadget,
+            )?;
+            Ok(())
+        }
+    }
+
+    #[test]
+    fn gm17_recursive_verifier_test() {
+        let num_inputs = 100;
+        let num_constraints = num_inputs;
+        let rng = &mut test_rng();
+        let mut inputs: Vec<Option<MNT4Fq>> = Vec::with_capacity(num_inputs);
+        for _ in 0..num_inputs {
+            inputs.push(Some(rng.gen()));
+        }
+
+        // Generate inner params and proof.
+        let inner_params = {
+            let c = Bench::<MNT4Fq> {
+                inputs: vec![None; num_inputs],
+                num_constraints,
+            };
+
+            generate_random_parameters(c, rng).unwrap()
+        };
+
+        let inner_proof = {
+            // Create an instance of our circuit (with the
+            // witness)
+            let c = Bench {
+                inputs: inputs.clone(),
+                num_constraints,
+            };
+            // Create a groth16 proof with our parameters.
+            create_random_proof(c, &inner_params, rng).unwrap()
+        };
+
+        // Generate outer params and proof.
+        let params = {
+            let c = Wrapper {
+                inputs: inputs.clone(),
+                params: inner_params.clone(),
+                proof:  inner_proof.clone(),
+            };
+
+            generate_random_parameters(c, rng).unwrap()
+        };
+
+        {
+            let proof = {
+                // Create an instance of our circuit (with the
+                // witness)
+                let c = Wrapper {
+                    inputs: inputs.clone(),
+                    params: inner_params.clone(),
+                    proof:  inner_proof.clone(),
+                };
+                // Create a groth16 proof with our parameters.
+                create_random_proof(c, &params, rng).unwrap()
+            };
+
+            let mut cs = TestConstraintSystem::<MNT4Fq>::new();
+
+            let inputs: Vec<_> = inputs.into_iter().map(|input| input.unwrap()).collect();
+            let mut input_gadgets = Vec::new();
+
+            {
+                let bigint_size = <MNT4FqParameters as FpParameters>::BigInt::NUM_LIMBS * 64;
+                let mut input_bits = Vec::new();
+                let mut cs = cs.ns(|| "Allocate Input");
+                for (i, input) in inputs.into_iter().enumerate() {
+                    let input_gadget =
+                        FpGadget::alloc_input(cs.ns(|| format!("Input {}", i)), || Ok(input))
+                            .unwrap();
+                    let mut fp_bits = input_gadget
+                        .to_bits(cs.ns(|| format!("To bits {}", i)))
+                        .unwrap();
+
+                    // FpGadget::to_bits outputs a big-endian binary representation of
+                    // fe_gadget's value, so we have to reverse it to get the little-endian
+                    // form.
+                    fp_bits.reverse();
+
+                    // Use 320 bits per element.
+                    for _ in fp_bits.len()..bigint_size {
+                        fp_bits.push(Boolean::constant(false));
+                    }
+                    input_bits.extend_from_slice(&fp_bits);
+                }
+
+                // Pack input bits into field elements of the underlying circuit.
+                let max_size = 8 * (<MNT6FqParameters as FpParameters>::CAPACITY / 8) as usize;
+                let max_size = max_size as usize;
+                let bigint_size = <MNT6FqParameters as FpParameters>::BigInt::NUM_LIMBS * 64;
+                for chunk in input_bits.chunks(max_size) {
+                    let mut chunk = chunk.to_vec();
+                    let len = chunk.len();
+                    for _ in len..bigint_size {
+                        chunk.push(Boolean::constant(false));
+                    }
+                    input_gadgets.push(chunk);
+                }
+                // assert!(!verify_proof(&pvk, &proof, &[a]).unwrap());
+            }
+
+            let vk_gadget = TestVkGadget2::alloc_input(cs.ns(|| "Vk"), || Ok(&params.vk)).unwrap();
+            let proof_gadget =
+                TestProofGadget2::alloc(cs.ns(|| "Proof"), || Ok(proof.clone())).unwrap();
+            println!("Time to verify!\n\n\n\n");
+            <TestVerifierGadget2 as NIZKVerifierGadget<TestProofSystem2, MNT4Fq>>::check_verify(
+                cs.ns(|| "Verify"),
+                &vk_gadget,
+                input_gadgets.iter(),
+                &proof_gadget,
+            )
+            .unwrap();
+            if !cs.is_satisfied() {
+                println!("=========================================================");
+                println!("Unsatisfied constraints:");
+                println!("{:?}", cs.which_is_unsatisfied().unwrap());
+                println!("=========================================================");
+            }
+
+            // cs.print_named_objects();
+            assert!(cs.is_satisfied());
+        }
+    }
+}

crypto-primitives/src/nizk/groth16/constraints.rs

@@ -479,3 +479,266 @@ mod test {
         }
     }
 }
+
+#[cfg(test)]
+mod test_recursive {
+    use groth16::*;
+    use r1cs_core::{ConstraintSynthesizer, ConstraintSystem, SynthesisError};
+
+    use super::*;
+    use algebra::{
+        fields::FpParameters,
+        mnt4_298::{Fq as MNT4Fq, FqParameters as MNT4FqParameters, Fr as MNT4Fr, MNT4_298},
+        mnt6_298::{Fq as MNT6Fq, FqParameters as MNT6FqParameters, Fr as MNT6Fr, MNT6_298},
+        test_rng, BigInteger, PrimeField,
+    };
+    use r1cs_std::{
+        fields::fp::FpGadget, mnt4_298::PairingGadget as MNT4_298PairingGadget,
+        mnt6_298::PairingGadget as MNT6_298PairingGadget,
+        test_constraint_system::TestConstraintSystem, uint8::UInt8,
+    };
+    use rand::Rng;
+
+    type TestProofSystem1 = Groth16<MNT6_298, Bench<MNT4Fq>, MNT6Fr>;
+    type TestVerifierGadget1 = Groth16VerifierGadget<MNT6_298, MNT6Fq, MNT6_298PairingGadget>;
+    type TestProofGadget1 = ProofGadget<MNT6_298, MNT6Fq, MNT6_298PairingGadget>;
+    type TestVkGadget1 = VerifyingKeyGadget<MNT6_298, MNT6Fq, MNT6_298PairingGadget>;
+
+    type TestProofSystem2 = Groth16<MNT4_298, Wrapper, MNT4Fr>;
+    type TestVerifierGadget2 = Groth16VerifierGadget<MNT4_298, MNT4Fq, MNT4_298PairingGadget>;
+    type TestProofGadget2 = ProofGadget<MNT4_298, MNT4Fq, MNT4_298PairingGadget>;
+    type TestVkGadget2 = VerifyingKeyGadget<MNT4_298, MNT4Fq, MNT4_298PairingGadget>;
+
+    #[derive(Clone)]
+    struct Bench<F: Field> {
+        inputs:          Vec<Option<F>>,
+        num_constraints: usize,
+    }
+
+    impl<F: Field> ConstraintSynthesizer<F> for Bench<F> {
+        fn generate_constraints<CS: ConstraintSystem<F>>(
+            self,
+            cs: &mut CS,
+        ) -> Result<(), SynthesisError> {
+            assert!(self.inputs.len() >= 2);
+            assert!(self.num_constraints >= self.inputs.len());
+
+            let mut variables: Vec<_> = Vec::with_capacity(self.inputs.len());
+            for (i, input) in self.inputs.into_iter().enumerate() {
+                let input_var = cs.alloc_input(
+                    || format!("Input {}", i),
+                    || input.ok_or(SynthesisError::AssignmentMissing),
+                )?;
+                variables.push((input, input_var));
+            }
+
+            for i in 0..self.num_constraints {
+                let new_entry = {
+                    let (input_1_val, input_1_var) = variables[i];
+                    let (input_2_val, input_2_var) = variables[i + 1];
+                    let result_val = input_1_val
+                        .and_then(|input_1| input_2_val.map(|input_2| input_1 * &input_2));
+                    let result_var = cs.alloc(
+                        || format!("Result {}", i),
+                        || result_val.ok_or(SynthesisError::AssignmentMissing),
+                    )?;
+                    cs.enforce(
+                        || format!("Enforce constraint {}", i),
+                        |lc| lc + input_1_var,
+                        |lc| lc + input_2_var,
+                        |lc| lc + result_var,
+                    );
+                    (result_val, result_var)
+                };
+                variables.push(new_entry);
+            }
+            Ok(())
+        }
+    }
+
+    struct Wrapper {
+        inputs: Vec<Option<MNT4Fq>>,
+        params: Parameters<MNT6_298>,
+        proof:  Proof<MNT6_298>,
+    }
+
+    impl ConstraintSynthesizer<MNT6Fq> for Wrapper {
+        fn generate_constraints<CS: ConstraintSystem<MNT6Fq>>(
+            self,
+            cs: &mut CS,
+        ) -> Result<(), SynthesisError> {
+            let params = self.params;
+            let proof = self.proof;
+            let inputs: Vec<_> = self
+                .inputs
+                .into_iter()
+                .map(|input| input.unwrap())
+                .collect();
+            let input_gadgets;
+
+            {
+                let mut cs = cs.ns(|| "Allocate Input");
+                // Chain all input values in one large byte array.
+                let input_bytes = inputs
+                    .clone()
+                    .into_iter()
+                    .flat_map(|input| {
+                        input
+                            .into_repr()
+                            .as_ref()
+                            .iter()
+                            .flat_map(|l| l.to_le_bytes().to_vec())
+                            .collect::<Vec<_>>()
+                    })
+                    .collect::<Vec<_>>();
+
+                // Allocate this byte array as input packed into field elements.
+                let input_bytes = UInt8::alloc_input_vec(cs.ns(|| "Input"), &input_bytes[..])?;
+                // 40 byte
+                let element_size = <MNT4FqParameters as FpParameters>::BigInt::NUM_LIMBS * 8;
+                input_gadgets = input_bytes
+                    .chunks(element_size)
+                    .map(|chunk| {
+                        chunk
+                            .iter()
+                            .flat_map(|byte| byte.into_bits_le())
+                            .collect::<Vec<_>>()
+                    })
+                    .collect::<Vec<_>>();
+            }
+
+            let vk_gadget = TestVkGadget1::alloc(cs.ns(|| "Vk"), || Ok(&params.vk))?;
+            let proof_gadget =
+                TestProofGadget1::alloc(cs.ns(|| "Proof"), || Ok(proof.clone())).unwrap();
+            <TestVerifierGadget1 as NIZKVerifierGadget<TestProofSystem1, MNT6Fq>>::check_verify(
+                cs.ns(|| "Verify"),
+                &vk_gadget,
+                input_gadgets.iter(),
+                &proof_gadget,
+            )?;
+            Ok(())
+        }
+    }
+
+    #[test]
+    fn groth16_recursive_verifier_test() {
+        let num_inputs = 100;
+        let num_constraints = num_inputs;
+        let rng = &mut test_rng();
+        let mut inputs: Vec<Option<MNT4Fq>> = Vec::with_capacity(num_inputs);
+        for _ in 0..num_inputs {
+            inputs.push(Some(rng.gen()));
+        }
+
+        // Generate inner params and proof.
+        let inner_params = {
+            let c = Bench::<MNT4Fq> {
+                inputs: vec![None; num_inputs],
+                num_constraints,
+            };
+
+            generate_random_parameters(c, rng).unwrap()
+        };
+
+        let inner_proof = {
+            // Create an instance of our circuit (with the
+            // witness)
+            let c = Bench {
+                inputs: inputs.clone(),
+                num_constraints,
+            };
+            // Create a groth16 proof with our parameters.
+            create_random_proof(c, &inner_params, rng).unwrap()
+        };
+
+        // Generate outer params and proof.
+        let params = {
+            let c = Wrapper {
+                inputs: inputs.clone(),
+                params: inner_params.clone(),
+                proof:  inner_proof.clone(),
+            };
+
+            generate_random_parameters(c, rng).unwrap()
+        };
+
+        {
+            let proof = {
+                // Create an instance of our circuit (with the
+                // witness)
+                let c = Wrapper {
+                    inputs: inputs.clone(),
+                    params: inner_params.clone(),
+                    proof:  inner_proof.clone(),
+                };
+                // Create a groth16 proof with our parameters.
+                create_random_proof(c, &params, rng).unwrap()
+            };
+
+            let mut cs = TestConstraintSystem::<MNT4Fq>::new();
+
+            let inputs: Vec<_> = inputs.into_iter().map(|input| input.unwrap()).collect();
+            let mut input_gadgets = Vec::new();
+
+            {
+                let bigint_size = <MNT4FqParameters as FpParameters>::BigInt::NUM_LIMBS * 64;
+                let mut input_bits = Vec::new();
+                let mut cs = cs.ns(|| "Allocate Input");
+                for (i, input) in inputs.into_iter().enumerate() {
+                    let input_gadget =
+                        FpGadget::alloc_input(cs.ns(|| format!("Input {}", i)), || Ok(input))
+                            .unwrap();
+                    let mut fp_bits = input_gadget
+                        .to_bits(cs.ns(|| format!("To bits {}", i)))
+                        .unwrap();
+
+                    // FpGadget::to_bits outputs a big-endian binary representation of
+                    // fe_gadget's value, so we have to reverse it to get the little-endian
+                    // form.
+                    fp_bits.reverse();
+
+                    // Use 320 bits per element.
+                    for _ in fp_bits.len()..bigint_size {
+                        fp_bits.push(Boolean::constant(false));
+                    }
+                    input_bits.extend_from_slice(&fp_bits);
+                }
+
+                // Pack input bits into field elements of the underlying circuit.
+                let max_size = 8 * (<MNT6FqParameters as FpParameters>::CAPACITY / 8) as usize;
+                let max_size = max_size as usize;
+                let bigint_size = <MNT6FqParameters as FpParameters>::BigInt::NUM_LIMBS * 64;
+                for chunk in input_bits.chunks(max_size) {
+                    let mut chunk = chunk.to_vec();
+                    let len = chunk.len();
+                    for _ in len..bigint_size {
+                        chunk.push(Boolean::constant(false));
+                    }
+                    input_gadgets.push(chunk);
+                }
+                // assert!(!verify_proof(&pvk, &proof, &[a]).unwrap());
+            }
+
+            let vk_gadget = TestVkGadget2::alloc_input(cs.ns(|| "Vk"), || Ok(&params.vk)).unwrap();
+            let proof_gadget =
+                TestProofGadget2::alloc(cs.ns(|| "Proof"), || Ok(proof.clone())).unwrap();
+            println!("Time to verify!\n\n\n\n");
+            <TestVerifierGadget2 as NIZKVerifierGadget<TestProofSystem2, MNT4Fq>>::check_verify(
+                cs.ns(|| "Verify"),
+                &vk_gadget,
+                input_gadgets.iter(),
+                &proof_gadget,
+            )
+            .unwrap();
+            if !cs.is_satisfied() {
+                println!("=========================================================");
+                println!("Unsatisfied constraints:");
+                println!("{:?}", cs.which_is_unsatisfied().unwrap());
+                println!("=========================================================");
+            }
+
+            // cs.print_named_objects();
+            assert!(cs.is_satisfied());
+        }
+    }
+}

r1cs-std/Cargo.toml

@@ -35,12 +35,16 @@ algebra = { path = "../algebra", default-features = false, features = [ "bls12_3
 
 [features]
 default = ["std"]
-full = [ "bls12_377", "jubjub", "edwards_bls12", "edwards_sw6", ]
+full = [ "bls12_377", "jubjub", "edwards_bls12", "edwards_sw6", "mnt4_298", "mnt4_753", "mnt6_298", "mnt6_753" ]
 
 bls12_377 = [ "algebra/bls12_377" ]
 jubjub = [ "algebra/jubjub" ]
 edwards_bls12 = [ "algebra/edwards_bls12" ]
 edwards_sw6 = [ "algebra/edwards_sw6" ]
+mnt4_298 = [ "algebra/mnt4_298" ]
+mnt4_753 = [ "algebra/mnt4_753" ]
+mnt6_298 = [ "algebra/mnt6_298" ]
+mnt6_753 = [ "algebra/mnt6_753" ]
 
 std = [ "algebra/std" ]
 parallel = [ "std", "algebra/parallel" ]

r1cs-std/src/fields/fp/mod.rs

@@ -270,6 +270,10 @@ impl<F: PrimeField> FieldGadget<F, F> for FpGadget<F> {
         1
     }
 
+    fn cost_of_mul_equals() -> usize {
+        1
+    }
+
     fn cost_of_inv() -> usize {
         1
     }

r1cs-std/src/fields/fp12.rs

@@ -642,6 +642,10 @@ where
         unimplemented!()
     }
 
+    fn cost_of_mul_equals() -> usize {
+        Self::cost_of_mul()
+    }
+
     fn cost_of_inv() -> usize {
         Self::cost_of_mul() + <Self as EqGadget<ConstraintF>>::cost()
     }

r1cs-std/src/fields/fp2.rs

@@ -453,6 +453,10 @@ impl<P: Fp2Parameters<Fp = ConstraintF>, ConstraintF: PrimeField> FieldGadget<Fp
         3
     }
 
+    fn cost_of_mul_equals() -> usize {
+        Self::cost_of_mul()
+    }
+
     fn cost_of_inv() -> usize {
         2
     }

r1cs-std/src/fields/fp3.rs

@@ -0,0 +1,932 @@
+use algebra::{
+    fields::{
+        fp3::{Fp3, Fp3Parameters},
+        Field,
+    },
+    PrimeField, SquareRootField,
+};
+use core::{borrow::Borrow, marker::PhantomData};
+use r1cs_core::{ConstraintSystem, ConstraintVar, SynthesisError};
+
+use crate::{fields::fp::FpGadget, prelude::*, Assignment, Vec};
+
+#[derive(Derivative)]
+#[derivative(Debug(
+    bound = "P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField"
+))]
+#[must_use]
+pub struct Fp3Gadget<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+{
+    pub c0:  FpGadget<ConstraintF>,
+    pub c1:  FpGadget<ConstraintF>,
+    pub c2:  FpGadget<ConstraintF>,
+    #[derivative(Debug = "ignore")]
+    _params: PhantomData<P>,
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    Fp3Gadget<P, ConstraintF>
+{
+    #[inline]
+    pub fn new(
+        c0: FpGadget<ConstraintF>,
+        c1: FpGadget<ConstraintF>,
+        c2: FpGadget<ConstraintF>,
+    ) -> Self {
+        Self {
+            c0,
+            c1,
+            c2,
+            _params: PhantomData,
+        }
+    }
+    /// Multiply a FpGadget by quadratic nonresidue P::NONRESIDUE.
+    #[inline]
+    pub fn mul_fp_gadget_by_nonresidue<CS: ConstraintSystem<ConstraintF>>(
+        cs: CS,
+        fe: &FpGadget<ConstraintF>,
+    ) -> Result<FpGadget<ConstraintF>, SynthesisError> {
+        fe.mul_by_constant(cs, &P::NONRESIDUE)
+    }
+
+    /// Multiply a Fp3Gadget by an element of fp.
+    #[inline]
+    pub fn mul_by_fp_constant_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        fe: &P::Fp,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.mul_by_constant_in_place(cs.ns(|| "c0"), fe)?;
+        self.c1.mul_by_constant_in_place(cs.ns(|| "c1"), fe)?;
+        self.c2.mul_by_constant_in_place(cs.ns(|| "c2"), fe)?;
+        Ok(self)
+    }
+
+    /// Multiply a Fp3Gadget by an element of fp.
+    #[inline]
+    pub fn mul_by_fp_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        fe: &P::Fp,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.mul_by_fp_constant_in_place(cs, fe)?;
+        Ok(result)
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    FieldGadget<Fp3<P>, ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    type Variable = (
+        ConstraintVar<ConstraintF>,
+        ConstraintVar<ConstraintF>,
+        ConstraintVar<ConstraintF>,
+    );
+
+    #[inline]
+    fn get_value(&self) -> Option<Fp3<P>> {
+        match (
+            self.c0.get_value(),
+            self.c1.get_value(),
+            self.c2.get_value(),
+        ) {
+            (Some(c0), Some(c1), Some(c2)) => Some(Fp3::new(c0, c1, c2)),
+            (..) => None,
+        }
+    }
+
+    #[inline]
+    fn get_variable(&self) -> Self::Variable {
+        (
+            self.c0.get_variable(),
+            self.c1.get_variable(),
+            self.c2.get_variable(),
+        )
+    }
+
+    #[inline]
+    fn zero<CS: ConstraintSystem<ConstraintF>>(mut cs: CS) -> Result<Self, SynthesisError> {
+        let c0 = FpGadget::<ConstraintF>::zero(cs.ns(|| "c0"))?;
+        let c1 = FpGadget::<ConstraintF>::zero(cs.ns(|| "c1"))?;
+        let c2 = FpGadget::<ConstraintF>::zero(cs.ns(|| "c2"))?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn one<CS: ConstraintSystem<ConstraintF>>(mut cs: CS) -> Result<Self, SynthesisError> {
+        let c0 = FpGadget::<ConstraintF>::one(cs.ns(|| "c0"))?;
+        let c1 = FpGadget::<ConstraintF>::zero(cs.ns(|| "c1"))?;
+        let c2 = FpGadget::<ConstraintF>::zero(cs.ns(|| "c2"))?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn conditionally_add_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        bit: &Boolean,
+        coeff: Fp3<P>,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self
+            .c0
+            .conditionally_add_constant(cs.ns(|| "c0"), bit, coeff.c0)?;
+        let c1 = self
+            .c1
+            .conditionally_add_constant(cs.ns(|| "c1"), bit, coeff.c1)?;
+        let c2 = self
+            .c2
+            .conditionally_add_constant(cs.ns(|| "c2"), bit, coeff.c2)?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn add<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.add(&mut cs.ns(|| "add c0"), &other.c0)?;
+        let c1 = self.c1.add(&mut cs.ns(|| "add c1"), &other.c1)?;
+        let c2 = self.c2.add(&mut cs.ns(|| "add c2"), &other.c2)?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn sub<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.sub(&mut cs.ns(|| "sub c0"), &other.c0)?;
+        let c1 = self.c1.sub(&mut cs.ns(|| "sub c1"), &other.c1)?;
+        let c2 = self.c2.sub(&mut cs.ns(|| "sub c2"), &other.c2)?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn negate<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.negate(&mut cs.ns(|| "negate c0"))?;
+        let c1 = self.c1.negate(&mut cs.ns(|| "negate c1"))?;
+        let c2 = self.c2.negate(&mut cs.ns(|| "negate c2"))?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn negate_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.negate_in_place(&mut cs.ns(|| "negate c0"))?;
+        self.c1.negate_in_place(&mut cs.ns(|| "negate c1"))?;
+        self.c2.negate_in_place(&mut cs.ns(|| "negate c2"))?;
+        Ok(self)
+    }
+
+    /// Use the Toom-Cook-3x method to compute multiplication.
+    #[inline]
+    fn mul<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        // Uses Toom-Cook-3x multiplication from
+        //
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        //    Devegili, OhEigeartaigh, Scott, Dahab
+
+        // v0 = a(0)b(0)   = a0 * b0
+        let v0 = self.c0.mul(&mut cs.ns(|| "Calc v0"), &other.c0)?;
+
+        // v1 = a(1)b(1)   = (a0 + a1 + a2)(b0 + b1 + b2)
+        let v1 = {
+            let mut v1_cs = cs.ns(|| "compute v1");
+            let a0_plus_a1_plus_a2 = self
+                .c0
+                .add(v1_cs.ns(|| "a0 + a1"), &self.c1)?
+                .add(v1_cs.ns(|| "a0 + a1 + a2"), &self.c2)?;
+            let b0_plus_b1_plus_b2 = other
+                .c0
+                .add(v1_cs.ns(|| "b0 + b1"), &other.c1)?
+                .add(v1_cs.ns(|| "b0 + b1 + b2"), &other.c2)?;
+
+            a0_plus_a1_plus_a2.mul(
+                v1_cs.ns(|| "(a0 + a1 + a2)(b0 + b1 + b2)"),
+                &b0_plus_b1_plus_b2,
+            )?
+        };
+
+        // v2 = a(−1)b(−1) = (a0 − a1 + a2)(b0 − b1 + b2)
+        let v2 = {
+            let mut v2_cs = cs.ns(|| "compute v2");
+
+            let a0_minus_a1_plus_a2 = self
+                .c0
+                .sub(v2_cs.ns(|| "a0 - a1"), &self.c1)?
+                .add(v2_cs.ns(|| "a0 - a1 + a2"), &self.c2)?;
+
+            let b0_minus_b1_plus_b2 = other
+                .c0
+                .sub(v2_cs.ns(|| "b0 - b1"), &other.c1)?
+                .add(v2_cs.ns(|| "b0 - b1 + b2"), &other.c2)?;
+
+            a0_minus_a1_plus_a2.mul(
+                v2_cs.ns(|| "(a0 - a1 + a2)(b0 - b1 + b2)"),
+                &b0_minus_b1_plus_b2,
+            )?
+        };
+
+        // v3 = a(2)b(2)   = (a0 + 2a1 + 4a2)(b0 + 2b1 + 4b2)
+        let v3 = {
+            let v3_cs = &mut cs.ns(|| "compute v3");
+
+            let a1_double = self.c1.double(v3_cs.ns(|| "2 * a1"))?;
+            let a2_quad = self
+                .c2
+                .double(v3_cs.ns(|| "2 * a2"))?
+                .double(v3_cs.ns(|| "4 * a2"))?;
+
+            let a0_plus_2_a1_plus_4_a2 = self
+                .c0
+                .add(v3_cs.ns(|| "a0 + 2a1"), &a1_double)?
+                .add(v3_cs.ns(|| "a0 + 2a1 + 4a2"), &a2_quad)?;
+
+            let b1_double = other.c1.double(v3_cs.ns(|| "2 * b1"))?;
+            let b2_quad = other
+                .c2
+                .double(v3_cs.ns(|| "2 * b2"))?
+                .double(v3_cs.ns(|| "4 * b2"))?;
+            let b0_plus_2_b1_plus_4_b2 = other
+                .c0
+                .add(v3_cs.ns(|| "b0 + 2b1"), &b1_double)?
+                .add(v3_cs.ns(|| "b0 + 2b1 + 4b2"), &b2_quad)?;
+
+            a0_plus_2_a1_plus_4_a2.mul(
+                v3_cs.ns(|| "(a0 + 2a1 + 4a2)(b0 + 2b1 + 4b2)"),
+                &b0_plus_2_b1_plus_4_b2,
+            )?
+        };
+
+        // v4 = a(∞)b(∞)   = a2 * b2
+        let v4 = self.c2.mul(cs.ns(|| "v2: a2 * b2"), &other.c2)?;
+
+        let two = P::Fp::one().double();
+        let six = two.double() + &two;
+        let mut two_and_six = [two, six];
+        algebra::fields::batch_inversion(&mut two_and_six);
+        let (two_inverse, six_inverse) = (two_and_six[0], two_and_six[1]);
+
+        let half_v0 = v0.mul_by_constant(cs.ns(|| "half_v0"), &two_inverse)?;
+        let half_v1 = v1.mul_by_constant(cs.ns(|| "half_v1"), &two_inverse)?;
+        let one_sixth_v2 = v2.mul_by_constant(cs.ns(|| "v2_by_six"), &six_inverse)?;
+        let one_sixth_v3 = v3.mul_by_constant(cs.ns(|| "v3_by_six"), &six_inverse)?;
+        let two_v4 = v4.double(cs.ns(|| "2 * v4"))?;
+
+        // c0 = v0 + β((1/2)v0 − (1/2)v1 − (1/6)v2 + (1/6)v3 − 2v4)
+        let c0 = {
+            let c0_cs = &mut cs.ns(|| "c0");
+
+            // No constraints, only get a linear combination back.
+            let temp = half_v0
+                .sub(c0_cs.ns(|| "sub1"), &half_v1)?
+                .sub(c0_cs.ns(|| "sub2"), &one_sixth_v2)?
+                .add(c0_cs.ns(|| "add3"), &one_sixth_v3)?
+                .sub(c0_cs.ns(|| "sub4"), &two_v4)?;
+            let non_residue_times_inner =
+                temp.mul_by_constant(&mut c0_cs.ns(|| "mul5"), &P::NONRESIDUE)?;
+            v0.add(c0_cs.ns(|| "add6"), &non_residue_times_inner)?
+        };
+
+        // −(1/2)v0 + v1 − (1/3)v2 − (1/6)v3 + 2v4 + βv4
+        let c1 = {
+            let c1_cs = &mut cs.ns(|| "c1");
+            let one_third_v2 = one_sixth_v2.double(&mut c1_cs.ns(|| "v2_by_3"))?;
+            let non_residue_v4 =
+                v4.mul_by_constant(&mut c1_cs.ns(|| "mul_by_beta"), &P::NONRESIDUE)?;
+
+            let result = half_v0
+                .negate(c1_cs.ns(|| "neg1"))?
+                .add(c1_cs.ns(|| "add2"), &v1)?
+                .sub(c1_cs.ns(|| "sub3"), &one_third_v2)?
+                .sub(c1_cs.ns(|| "sub4"), &one_sixth_v3)?
+                .add(c1_cs.ns(|| "sub5"), &two_v4)?
+                .add(c1_cs.ns(|| "sub6"), &non_residue_v4)?;
+            result
+        };
+
+        // -v0 + (1/2)v1 + (1/2)v2 −v4
+        let c2 = {
+            let c2_cs = &mut cs.ns(|| "c2");
+            let half_v2 = v2.mul_by_constant(&mut c2_cs.ns(|| "mul1"), &two_inverse)?;
+            let result = half_v1
+                .add(c2_cs.ns(|| "add1"), &half_v2)?
+                .sub(c2_cs.ns(|| "sub1"), &v4)?
+                .sub(c2_cs.ns(|| "sub2"), &v0)?;
+            result
+        };
+
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    /// Use the Toom-Cook-3x method to compute multiplication.
+    #[inline]
+    fn square<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        // Uses Toom-Cook-3x multiplication from
+        //
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        //    Devegili, OhEigeartaigh, Scott, Dahab
+
+        // v0 = a(0)^2 = a0^2
+        let v0 = self.c0.square(&mut cs.ns(|| "Calc v0"))?;
+
+        // v1 = a(1)^2 = (a0 + a1 + a2)^2
+        let v1 = {
+            let a0_plus_a1_plus_a2 = self
+                .c0
+                .add(cs.ns(|| "a0 + a1"), &self.c1)?
+                .add(cs.ns(|| "a0 + a1 + a2"), &self.c2)?;
+            a0_plus_a1_plus_a2.square(&mut cs.ns(|| "(a0 + a1 + a2)^2"))?
+        };
+
+        // v2 = a(−1)^2 = (a0 − a1 + a2)^2
+        let v2 = {
+            let a0_minus_a1_plus_a2 = self
+                .c0
+                .sub(cs.ns(|| "a0 - a1"), &self.c1)?
+                .add(cs.ns(|| "a0 - a2 + a2"), &self.c2)?;
+            a0_minus_a1_plus_a2.square(&mut cs.ns(|| "(a0 - a1 + a2)^2"))?
+        };
+
+        // v3 = a(2)^2 = (a0 + 2a1 + 4a2)^2
+        let v3 = {
+            let a1_double = self.c1.double(cs.ns(|| "2a1"))?;
+            let a2_quad = self.c2.double(cs.ns(|| "2a2"))?.double(cs.ns(|| "4a2"))?;
+            let a0_plus_2_a1_plus_4_a2 = self
+                .c0
+                .add(cs.ns(|| "a0 + 2a1"), &a1_double)?
+                .add(cs.ns(|| "a0 + 2a1 + 4a2"), &a2_quad)?;
+
+            a0_plus_2_a1_plus_4_a2.square(&mut cs.ns(|| "(a0 + 2a1 + 4a2)^2"))?
+        };
+
+        // v4 = a(∞)^2 = a2^2
+        let v4 = self.c2.square(&mut cs.ns(|| "a2^2"))?;
+
+        let two = P::Fp::one().double();
+        let six = two.double() + &two;
+        let mut two_and_six = [two, six];
+        algebra::fields::batch_inversion(&mut two_and_six);
+        let (two_inverse, six_inverse) = (two_and_six[0], two_and_six[1]);
+
+        let half_v0 = v0.mul_by_constant(cs.ns(|| "half_v0"), &two_inverse)?;
+        let half_v1 = v1.mul_by_constant(cs.ns(|| "half_v1"), &two_inverse)?;
+        let one_sixth_v2 = v2.mul_by_constant(cs.ns(|| "one_sixth_v2"), &six_inverse)?;
+        let one_sixth_v3 = v3.mul_by_constant(cs.ns(|| "one_sixth_v3"), &six_inverse)?;
+        let two_v4 = v4.double(cs.ns(|| "double_v4"))?;
+
+        // c0 = v0 + β((1/2)v0 − (1/2)v1 − (1/6)v2 + (1/6)v3 − 2v4)
+        let c0 = {
+            let mut c0_cs = cs.ns(|| "c0");
+            // No constraints, only get a linear combination back.
+            let inner = half_v0
+                .sub(c0_cs.ns(|| "sub1"), &half_v1)?
+                .sub(c0_cs.ns(|| "sub2"), &one_sixth_v2)?
+                .add(c0_cs.ns(|| "add3"), &one_sixth_v3)?
+                .sub(c0_cs.ns(|| "sub4"), &two_v4)?;
+            let non_residue_times_inner =
+                inner.mul_by_constant(c0_cs.ns(|| "mul_by_res"), &P::NONRESIDUE)?;
+            v0.add(c0_cs.ns(|| "add5"), &non_residue_times_inner)?
+        };
+
+        // −(1/2)v0 + v1 − (1/3)v2 − (1/6)v3 + 2v4 + βv4
+        let c1 = {
+            let mut c1_cs = cs.ns(|| "c1");
+            let one_third_v2 = one_sixth_v2.double(c1_cs.ns(|| "v2_by_3"))?;
+            let non_residue_v4 = v4.mul_by_constant(c1_cs.ns(|| "mul_by_res"), &P::NONRESIDUE)?;
+
+            half_v0
+                .negate(c1_cs.ns(|| "neg1"))?
+                .add(c1_cs.ns(|| "add1"), &v1)?
+                .sub(c1_cs.ns(|| "sub2"), &one_third_v2)?
+                .sub(c1_cs.ns(|| "sub3"), &one_sixth_v3)?
+                .add(c1_cs.ns(|| "add4"), &two_v4)?
+                .add(c1_cs.ns(|| "add5"), &non_residue_v4)?
+        };
+
+        // -v0 + (1/2)v1 + (1/2)v2 −v4
+        let c2 = {
+            let mut c2_cs = cs.ns(|| "c2");
+            let half_v2 = v2.mul_by_constant(c2_cs.ns(|| "half_v2"), &two_inverse)?;
+            half_v1
+                .add(c2_cs.ns(|| "add1"), &half_v2)?
+                .sub(c2_cs.ns(|| "sub1"), &v4)?
+                .sub(c2_cs.ns(|| "sub2"), &v0)?
+        };
+
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    // 18 constraints, we can probably do better but not sure it's worth it.
+    #[inline]
+    fn inverse<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        let inverse = Self::alloc(&mut cs.ns(|| "alloc inverse"), || {
+            self.get_value().and_then(|val| val.inverse()).get()
+        })?;
+        let one = Self::one(cs.ns(|| "one"))?;
+        inverse.mul_equals(cs.ns(|| "check inverse"), &self, &one)?;
+        Ok(inverse)
+    }
+
+    #[inline]
+    fn add_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Fp3<P>,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.add_constant(cs.ns(|| "c0"), &other.c0)?;
+        let c1 = self.c1.add_constant(cs.ns(|| "c1"), &other.c1)?;
+        let c2 = self.c2.add_constant(cs.ns(|| "c2"), &other.c2)?;
+
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn add_constant_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        other: &Fp3<P>,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.add_constant_in_place(cs.ns(|| "c0"), &other.c0)?;
+        self.c1.add_constant_in_place(cs.ns(|| "c1"), &other.c1)?;
+        self.c2.add_constant_in_place(cs.ns(|| "c2"), &other.c2)?;
+        Ok(self)
+    }
+
+    /// Use the Toom-Cook-3x method to compute multiplication.
+    #[inline]
+    fn mul_by_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Fp3<P>,
+    ) -> Result<Self, SynthesisError> {
+        // Uses Toom-Cook-3x multiplication from
+        //
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        //    Devegili, OhEigeartaigh, Scott, Dahab
+
+        // v0 = a(0)b(0)   = a0 * b0
+        let v0 = self.c0.mul_by_constant(cs.ns(|| "v0"), &other.c0)?;
+
+        // v1 = a(1)b(1)   = (a0 + a1 + a2)(b0 + b1 + b2)
+        let v1 = {
+            let mut v1_cs = cs.ns(|| "v1");
+            let mut a0_plus_a1_plus_a2 = self
+                .c0
+                .add(v1_cs.ns(|| "a0 + a1"), &self.c1)?
+                .add(v1_cs.ns(|| "a0 + a1 + a2"), &self.c2)?;
+            let b0_plus_b1_plus_b2 = other.c0 + &other.c1 + &other.c2;
+
+            a0_plus_a1_plus_a2.mul_by_constant_in_place(
+                v1_cs.ns(|| "(a0 + a1 + a2)*(b0 + b1 + b2)"),
+                &b0_plus_b1_plus_b2,
+            )?;
+            a0_plus_a1_plus_a2
+        };
+
+        // v2 = a(−1)b(−1) = (a0 − a1 + a2)(b0 − b1 + b2)
+        let mut v2 = {
+            let mut v2_cs = cs.ns(|| "v2");
+            let mut a0_minus_a1_plus_a2 = self
+                .c0
+                .sub(v2_cs.ns(|| "sub1"), &self.c1)?
+                .add(v2_cs.ns(|| "add2"), &self.c2)?;
+            let b0_minus_b1_plus_b2 = other.c0 - &other.c1 + &other.c2;
+            a0_minus_a1_plus_a2.mul_by_constant_in_place(
+                v2_cs.ns(|| "(a0 - a1 + a2)*(b0 - b1 + b2)"),
+                &b0_minus_b1_plus_b2,
+            )?;
+            a0_minus_a1_plus_a2
+        };
+
+        // v3 = a(2)b(2)   = (a0 + 2a1 + 4a2)(b0 + 2b1 + 4b2)
+        let mut v3 = {
+            let mut v3_cs = cs.ns(|| "v3");
+            let a1_double = self.c1.double(v3_cs.ns(|| "2a1"))?;
+            let a2_quad = self
+                .c2
+                .double(v3_cs.ns(|| "2a2"))?
+                .double(v3_cs.ns(|| "4a2"))?;
+            let mut a0_plus_2_a1_plus_4_a2 = self
+                .c0
+                .add(v3_cs.ns(|| "a0 + 2a1"), &a1_double)?
+                .add(v3_cs.ns(|| "a0 + 2a1 + 4a2"), &a2_quad)?;
+
+            let b1_double = other.c1.double();
+            let b2_quad = other.c2.double().double();
+            let b0_plus_2_b1_plus_4_b2 = other.c0 + &b1_double + &b2_quad;
+
+            a0_plus_2_a1_plus_4_a2.mul_by_constant_in_place(
+                v3_cs.ns(|| "(a0 + 2a1 + 4a2)*(b0 + 2b1 + 4b2)"),
+                &b0_plus_2_b1_plus_4_b2,
+            )?;
+            a0_plus_2_a1_plus_4_a2
+        };
+
+        // v4 = a(∞)b(∞)   = a2 * b2
+        let v4 = self.c2.mul_by_constant(cs.ns(|| "v4"), &other.c2)?;
+
+        let two = P::Fp::one().double();
+        let six = two.double() + &two;
+        let mut two_and_six = [two, six];
+        algebra::fields::batch_inversion(&mut two_and_six);
+        let (two_inverse, six_inverse) = (two_and_six[0], two_and_six[1]);
+
+        let mut half_v0 = v0.mul_by_constant(cs.ns(|| "half_v0"), &two_inverse)?;
+        let half_v1 = v1.mul_by_constant(cs.ns(|| "half_v1"), &two_inverse)?;
+        let mut one_sixth_v2 = v2.mul_by_constant(cs.ns(|| "v2_by_6"), &six_inverse)?;
+        let one_sixth_v3 = v3.mul_by_constant_in_place(cs.ns(|| "v3_by_6"), &six_inverse)?;
+        let two_v4 = v4.double(cs.ns(|| "2v4"))?;
+
+        // c0 = v0 + β((1/2)v0 − (1/2)v1 − (1/6)v2 + (1/6)v3 − 2v4)
+        let c0 = {
+            let mut c0_cs = cs.ns(|| "c0");
+
+            // No constraints, only get a linear combination back.
+            let mut inner = half_v0
+                .sub(c0_cs.ns(|| "sub1"), &half_v1)?
+                .sub(c0_cs.ns(|| "sub2"), &one_sixth_v2)?
+                .add(c0_cs.ns(|| "add3"), &one_sixth_v3)?
+                .sub(c0_cs.ns(|| "sub4"), &two_v4)?;
+            let non_residue_times_inner =
+                inner.mul_by_constant_in_place(&mut c0_cs, &P::NONRESIDUE)?;
+            v0.add(c0_cs.ns(|| "add5"), non_residue_times_inner)?
+        };
+
+        // −(1/2)v0 + v1 − (1/3)v2 − (1/6)v3 + 2v4 + βv4
+        let c1 = {
+            let mut c1_cs = cs.ns(|| "c1");
+            let one_third_v2 = one_sixth_v2.double_in_place(c1_cs.ns(|| "double1"))?;
+            let non_residue_v4 =
+                v4.mul_by_constant(c1_cs.ns(|| "mul_by_const1"), &P::NONRESIDUE)?;
+
+            half_v0
+                .negate_in_place(c1_cs.ns(|| "neg1"))?
+                .add(c1_cs.ns(|| "add1"), &v1)?
+                .sub(c1_cs.ns(|| "sub2"), one_third_v2)?
+                .sub(c1_cs.ns(|| "sub3"), &one_sixth_v3)?
+                .add(c1_cs.ns(|| "add4"), &two_v4)?
+                .add(c1_cs.ns(|| "add5"), &non_residue_v4)?
+        };
+
+        // -v0 + (1/2)v1 + (1/2)v2 −v4
+        let c2 = {
+            let mut c2_cs = cs.ns(|| "c2");
+            let half_v2 = v2.mul_by_constant_in_place(c2_cs.ns(|| "half_v2"), &two_inverse)?;
+            half_v1
+                .add(c2_cs.ns(|| "add1"), half_v2)?
+                .sub(c2_cs.ns(|| "sub2"), &v4)?
+                .sub(c2_cs.ns(|| "sub3"), &v0)?
+        };
+
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    fn frobenius_map<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        power: usize,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.frobenius_map_in_place(cs, power)?;
+        Ok(result)
+    }
+
+    fn frobenius_map_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        power: usize,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c1.mul_by_constant_in_place(
+            cs.ns(|| "c1_power"),
+            &P::FROBENIUS_COEFF_FP3_C1[power % 3],
+        )?;
+        self.c2.mul_by_constant_in_place(
+            cs.ns(|| "c2_power"),
+            &P::FROBENIUS_COEFF_FP3_C2[power % 3],
+        )?;
+
+        Ok(self)
+    }
+
+    fn cost_of_mul() -> usize {
+        5 * FpGadget::<ConstraintF>::cost_of_mul()
+    }
+
+    fn cost_of_inv() -> usize {
+        Self::cost_of_mul() + <Self as EqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField> PartialEq
+    for Fp3Gadget<P, ConstraintF>
+{
+    fn eq(&self, other: &Self) -> bool {
+        self.c0 == other.c0 && self.c1 == other.c1 && self.c2 == other.c2
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField> Eq
+    for Fp3Gadget<P, ConstraintF>
+{
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    EqGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    ConditionalEqGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    #[inline]
+    fn conditional_enforce_equal<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+        condition: &Boolean,
+    ) -> Result<(), SynthesisError> {
+        self.c0
+            .conditional_enforce_equal(&mut cs.ns(|| "c0"), &other.c0, condition)?;
+        self.c1
+            .conditional_enforce_equal(&mut cs.ns(|| "c1"), &other.c1, condition)?;
+        self.c2
+            .conditional_enforce_equal(&mut cs.ns(|| "c2"), &other.c2, condition)?;
+        Ok(())
+    }
+
+    fn cost() -> usize {
+        3 * <FpGadget<ConstraintF> as ConditionalEqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    NEqGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    #[inline]
+    fn enforce_not_equal<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<(), SynthesisError> {
+        self.c0.enforce_not_equal(&mut cs.ns(|| "c0"), &other.c0)?;
+        self.c1.enforce_not_equal(&mut cs.ns(|| "c1"), &other.c1)?;
+        self.c2.enforce_not_equal(&mut cs.ns(|| "c2"), &other.c2)?;
+        Ok(())
+    }
+
+    fn cost() -> usize {
+        3 * <FpGadget<ConstraintF> as NEqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    ToBitsGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    fn to_bits<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<Boolean>, SynthesisError> {
+        let mut c0 = self.c0.to_bits(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_bits(cs.ns(|| "c1"))?;
+        let mut c2 = self.c2.to_bits(cs.ns(|| "c2"))?;
+
+        c0.append(&mut c1);
+        c0.append(&mut c2);
+
+        Ok(c0)
+    }
+
+    fn to_non_unique_bits<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<Boolean>, SynthesisError> {
+        let mut c0 = self.c0.to_non_unique_bits(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_non_unique_bits(cs.ns(|| "c1"))?;
+        let mut c2 = self.c2.to_non_unique_bits(cs.ns(|| "c2"))?;
+
+        c0.append(&mut c1);
+        c0.append(&mut c2);
+
+        Ok(c0)
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    ToBytesGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    fn to_bytes<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c0 = self.c0.to_bytes(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_bytes(cs.ns(|| "c1"))?;
+        let mut c2 = self.c2.to_bytes(cs.ns(|| "c2"))?;
+
+        c0.append(&mut c1);
+        c0.append(&mut c2);
+
+        Ok(c0)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c0 = self.c0.to_non_unique_bytes(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_non_unique_bytes(cs.ns(|| "c1"))?;
+        let mut c2 = self.c2.to_non_unique_bytes(cs.ns(|| "c2"))?;
+
+        c0.append(&mut c1);
+        c0.append(&mut c2);
+
+        Ok(c0)
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField> Clone
+    for Fp3Gadget<P, ConstraintF>
+{
+    fn clone(&self) -> Self {
+        Self::new(self.c0.clone(), self.c1.clone(), self.c2.clone())
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    CondSelectGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    #[inline]
+    fn conditionally_select<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        cond: &Boolean,
+        first: &Self,
+        second: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = FpGadget::<ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c0"),
+            cond,
+            &first.c0,
+            &second.c0,
+        )?;
+        let c1 = FpGadget::<ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c1"),
+            cond,
+            &first.c1,
+            &second.c1,
+        )?;
+        let c2 = FpGadget::<ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c2"),
+            cond,
+            &first.c2,
+            &second.c2,
+        )?;
+
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    fn cost() -> usize {
+        3 * <FpGadget<ConstraintF> as CondSelectGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    TwoBitLookupGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    type TableConstant = Fp3<P>;
+    fn two_bit_lookup<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        b: &[Boolean],
+        c: &[Self::TableConstant],
+    ) -> Result<Self, SynthesisError> {
+        let c0s = c.iter().map(|f| f.c0).collect::<Vec<_>>();
+        let c1s = c.iter().map(|f| f.c1).collect::<Vec<_>>();
+        let c2s = c.iter().map(|f| f.c2).collect::<Vec<_>>();
+        let c0 = FpGadget::<ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c0"), b, &c0s)?;
+        let c1 = FpGadget::<ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c1"), b, &c1s)?;
+        let c2 = FpGadget::<ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c2"), b, &c2s)?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    fn cost() -> usize {
+        3 * <FpGadget<ConstraintF> as TwoBitLookupGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    ThreeBitCondNegLookupGadget<ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    type TableConstant = Fp3<P>;
+
+    fn three_bit_cond_neg_lookup<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        b: &[Boolean],
+        b0b1: &Boolean,
+        c: &[Self::TableConstant],
+    ) -> Result<Self, SynthesisError> {
+        let c0s = c.iter().map(|f| f.c0).collect::<Vec<_>>();
+        let c1s = c.iter().map(|f| f.c1).collect::<Vec<_>>();
+        let c2s = c.iter().map(|f| f.c2).collect::<Vec<_>>();
+        let c0 = FpGadget::<ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c0"),
+            b,
+            b0b1,
+            &c0s,
+        )?;
+        let c1 = FpGadget::<ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c1"),
+            b,
+            b0b1,
+            &c1s,
+        )?;
+        let c2 = FpGadget::<ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c2"),
+            b,
+            b0b1,
+            &c2s,
+        )?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    fn cost() -> usize {
+        3 * <FpGadget<ConstraintF> as ThreeBitCondNegLookupGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P: Fp3Parameters<Fp = ConstraintF>, ConstraintF: PrimeField + SquareRootField>
+    AllocGadget<Fp3<P>, ConstraintF> for Fp3Gadget<P, ConstraintF>
+{
+    #[inline]
+    fn alloc<F, T, CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        value_gen: F,
+    ) -> Result<Self, SynthesisError>
+    where
+        F: FnOnce() -> Result<T, SynthesisError>,
+        T: Borrow<Fp3<P>>,
+    {
+        let (c0, c1, c2) = match value_gen() {
+            Ok(fe) => {
+                let fe = *fe.borrow();
+                (Ok(fe.c0), Ok(fe.c1), Ok(fe.c2))
+            },
+            _ => (
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+            ),
+        };
+
+        let c0 = FpGadget::<ConstraintF>::alloc(&mut cs.ns(|| "c0"), || c0)?;
+        let c1 = FpGadget::<ConstraintF>::alloc(&mut cs.ns(|| "c1"), || c1)?;
+        let c2 = FpGadget::<ConstraintF>::alloc(&mut cs.ns(|| "c2"), || c2)?;
+        Ok(Self::new(c0, c1, c2))
+    }
+
+    #[inline]
+    fn alloc_input<F, T, CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        value_gen: F,
+    ) -> Result<Self, SynthesisError>
+    where
+        F: FnOnce() -> Result<T, SynthesisError>,
+        T: Borrow<Fp3<P>>,
+    {
+        let (c0, c1, c2) = match value_gen() {
+            Ok(fe) => {
+                let fe = *fe.borrow();
+                (Ok(fe.c0), Ok(fe.c1), Ok(fe.c2))
+            },
+            _ => (
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+            ),
+        };
+
+        let c0 = FpGadget::<ConstraintF>::alloc_input(&mut cs.ns(|| "c0"), || c0)?;
+        let c1 = FpGadget::<ConstraintF>::alloc_input(&mut cs.ns(|| "c1"), || c1)?;
+        let c2 = FpGadget::<ConstraintF>::alloc_input(&mut cs.ns(|| "c2"), || c2)?;
+        Ok(Self::new(c0, c1, c2))
+    }
+}

r1cs-std/src/fields/fp4.rs

@@ -0,0 +1,772 @@
+use algebra::{
+    fields::{Fp2, Fp2Parameters, Fp4, Fp4Parameters},
+    BigInteger, Field, One, PrimeField,
+};
+use core::{borrow::Borrow, marker::PhantomData};
+use r1cs_core::{ConstraintSystem, SynthesisError};
+
+use crate::{prelude::*, Assignment, Vec};
+
+type Fp2Gadget<P, ConstraintF> =
+    super::fp2::Fp2Gadget<<P as Fp4Parameters>::Fp2Params, ConstraintF>;
+type Fp2GadgetVariable<P, ConstraintF> = <Fp2Gadget<P, ConstraintF> as FieldGadget<
+    Fp2<<P as Fp4Parameters>::Fp2Params>,
+    ConstraintF,
+>>::Variable;
+
+#[derive(Derivative)]
+#[derivative(Debug(bound = "ConstraintF: PrimeField"))]
+#[must_use]
+pub struct Fp4Gadget<P, ConstraintF: PrimeField>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    pub c0:  Fp2Gadget<P, ConstraintF>,
+    pub c1:  Fp2Gadget<P, ConstraintF>,
+    #[derivative(Debug = "ignore")]
+    _params: PhantomData<P>,
+}
+
+impl<P, ConstraintF: PrimeField> Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    pub fn new(c0: Fp2Gadget<P, ConstraintF>, c1: Fp2Gadget<P, ConstraintF>) -> Self {
+        Self {
+            c0,
+            c1,
+            _params: PhantomData,
+        }
+    }
+
+    /// Multiply a Fp2Gadget by quadratic nonresidue P::NONRESIDUE.
+    #[inline]
+    pub fn mul_fp2_gadget_by_nonresidue<CS: ConstraintSystem<ConstraintF>>(
+        cs: CS,
+        fe: &Fp2Gadget<P, ConstraintF>,
+    ) -> Result<Fp2Gadget<P, ConstraintF>, SynthesisError> {
+        let new_c0 = Fp2Gadget::<P, ConstraintF>::mul_fp_gadget_by_nonresidue(cs, &fe.c1)?;
+        let new_c1 = fe.c0.clone();
+        Ok(Fp2Gadget::<P, ConstraintF>::new(new_c0, new_c1))
+    }
+
+    /// Multiply a Fp4Gadget by an element of fp.
+    #[inline]
+    pub fn mul_by_fp_constant_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        fe: &<<P as Fp4Parameters>::Fp2Params as Fp2Parameters>::Fp,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.mul_by_fp_constant_in_place(cs.ns(|| "c0"), fe)?;
+        self.c1.mul_by_fp_constant_in_place(cs.ns(|| "c1"), fe)?;
+        Ok(self)
+    }
+
+    /// Multiply a Fp4Gadget by an element of fp.
+    #[inline]
+    pub fn mul_by_fp_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        fe: &<<P as Fp4Parameters>::Fp2Params as Fp2Parameters>::Fp,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.mul_by_fp_constant_in_place(cs, fe)?;
+        Ok(result)
+    }
+
+    pub fn unitary_inverse<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        Ok(Self::new(self.c0.clone(), self.c1.negate(cs)?))
+    }
+
+    #[inline]
+    pub fn cyclotomic_exp<CS: ConstraintSystem<ConstraintF>, B: BigInteger>(
+        &self,
+        mut cs: CS,
+        exponent: &B,
+    ) -> Result<Self, SynthesisError> {
+        let mut res = Self::one(cs.ns(|| "one"))?;
+        let self_inverse = self.unitary_inverse(cs.ns(|| "unitary inverse"))?;
+
+        let mut found_nonzero = false;
+        let naf = exponent.find_wnaf();
+
+        for (i, &value) in naf.iter().rev().enumerate() {
+            if found_nonzero {
+                res.square_in_place(cs.ns(|| format!("square {}", i)))?;
+            }
+
+            if value != 0 {
+                found_nonzero = true;
+
+                if value > 0 {
+                    res.mul_in_place(cs.ns(|| format!("res *= self {}", i)), &self)?;
+                } else {
+                    res.mul_in_place(
+                        cs.ns(|| format!("res *= self_inverse {}", i)),
+                        &self_inverse,
+                    )?;
+                }
+            }
+        }
+
+        Ok(res)
+    }
+}
+
+impl<P, ConstraintF: PrimeField> FieldGadget<Fp4<P>, ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    type Variable = (
+        Fp2GadgetVariable<P, ConstraintF>,
+        Fp2GadgetVariable<P, ConstraintF>,
+    );
+
+    #[inline]
+    fn get_value(&self) -> Option<Fp4<P>> {
+        match (self.c0.get_value(), self.c1.get_value()) {
+            (Some(c0), Some(c1)) => Some(Fp4::new(c0, c1)),
+            (..) => None,
+        }
+    }
+
+    #[inline]
+    fn get_variable(&self) -> Self::Variable {
+        (self.c0.get_variable(), self.c1.get_variable())
+    }
+
+    #[inline]
+    fn zero<CS: ConstraintSystem<ConstraintF>>(mut cs: CS) -> Result<Self, SynthesisError> {
+        let c0 = Fp2Gadget::<P, ConstraintF>::zero(cs.ns(|| "c0"))?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::zero(cs.ns(|| "c1"))?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn one<CS: ConstraintSystem<ConstraintF>>(mut cs: CS) -> Result<Self, SynthesisError> {
+        let c0 = Fp2Gadget::<P, ConstraintF>::one(cs.ns(|| "c0"))?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::zero(cs.ns(|| "c1"))?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn conditionally_add_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        bit: &Boolean,
+        coeff: Fp4<P>,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self
+            .c0
+            .conditionally_add_constant(cs.ns(|| "c0"), bit, coeff.c0)?;
+        let c1 = self
+            .c1
+            .conditionally_add_constant(cs.ns(|| "c1"), bit, coeff.c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn add<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.add(&mut cs.ns(|| "add c0"), &other.c0)?;
+        let c1 = self.c1.add(&mut cs.ns(|| "add c1"), &other.c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn sub<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.sub(&mut cs.ns(|| "sub c0"), &other.c0)?;
+        let c1 = self.c1.sub(&mut cs.ns(|| "sub c1"), &other.c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn double<CS: ConstraintSystem<ConstraintF>>(&self, cs: CS) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.double_in_place(cs)?;
+        Ok(result)
+    }
+
+    #[inline]
+    fn double_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.double_in_place(&mut cs.ns(|| "double c0"))?;
+        self.c1.double_in_place(&mut cs.ns(|| "double c1"))?;
+        Ok(self)
+    }
+
+    #[inline]
+    fn negate<CS: ConstraintSystem<ConstraintF>>(&self, cs: CS) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.negate_in_place(cs)?;
+        Ok(result)
+    }
+
+    #[inline]
+    fn negate_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.negate_in_place(&mut cs.ns(|| "negate c0"))?;
+        self.c1.negate_in_place(&mut cs.ns(|| "negate c1"))?;
+        Ok(self)
+    }
+
+    #[inline]
+    fn mul<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        // Karatsuba multiplication for Fp4:
+        //     v0 = A.c0 * B.c0
+        //     v1 = A.c1 * B.c1
+        //     result.c0 = v0 + non_residue * v1
+        //     result.c1 = (A.c0 + A.c1) * (B.c0 + B.c1) - v0 - v1
+        // Enforced with 3 constraints:
+        //     A.c1 * B.c1 = v1
+        //     A.c0 * B.c0 = result.c0 - non_residue * v1
+        //     (A.c0+A.c1)*(B.c0+B.c1) = result.c1 + result.c0 + (1 - non_residue) * v1
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        // Devegili, OhEigeartaigh, Scott, Dahab
+        let mul_cs = &mut cs.ns(|| "mul");
+
+        let v0 = self.c0.mul(mul_cs.ns(|| "v0"), &other.c0)?;
+        let v1 = self.c1.mul(mul_cs.ns(|| "v1"), &other.c1)?;
+        let c0 = {
+            let non_residue_times_v1 =
+                Self::mul_fp2_gadget_by_nonresidue(mul_cs.ns(|| "first mul_by_nr"), &v1)?;
+            v0.add(mul_cs.ns(|| "v0 + beta * v1"), &non_residue_times_v1)?
+        };
+        let c1 = {
+            let a0_plus_a1 = self.c0.add(mul_cs.ns(|| "a0 + a1"), &self.c1)?;
+            let b0_plus_b1 = other.c0.add(mul_cs.ns(|| "b0 + b1"), &other.c1)?;
+            let a0_plus_a1_times_b0_plus_b1 =
+                a0_plus_a1.mul(&mut mul_cs.ns(|| "(a0 + a1) * (b0 + b1)"), &b0_plus_b1)?;
+            a0_plus_a1_times_b0_plus_b1
+                .sub(mul_cs.ns(|| "res - v0"), &v0)?
+                .sub(mul_cs.ns(|| "res - v0 - v1"), &v1)?
+        };
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn square<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        // From Libsnark/fp4_gadget.tcc
+        // Complex multiplication for Fp4:
+        //     v0 = A.c0 * A.c1
+        //     result.c0 = (A.c0 + A.c1) * (A.c0 + non_residue * A.c1) - (1 +
+        // non_residue) * v0     result.c1 = 2 * v0
+        // Enforced with 2 constraints:
+        //     (2*A.c0) * A.c1 = result.c1
+        //     (A.c0 + A.c1) * (A.c0 + non_residue * A.c1) = result.c0 + result.c1 * (1
+        // + non_residue)/2 Reference:
+        //     "Multiplication and Squaring on Pairing-Friendly Fields"
+        //     Devegili, OhEigeartaigh, Scott, Dahab
+
+        let mut v0 = self.c0.mul(cs.ns(|| "v0"), &self.c1)?;
+        let a0_plus_a1 = self.c0.add(cs.ns(|| "a0 + a1"), &self.c1)?;
+
+        let non_residue_c1 =
+            Self::mul_fp2_gadget_by_nonresidue(cs.ns(|| "non_residue * a1"), &self.c1)?;
+        let a0_plus_non_residue_c1 = self
+            .c0
+            .add(cs.ns(|| "a0 + non_residue * a1"), &non_residue_c1)?;
+        let one_plus_non_residue_v0 =
+            Self::mul_fp2_gadget_by_nonresidue(cs.ns(|| "non_residue * v0"), &v0)?
+                .add(cs.ns(|| "plus v0"), &v0)?;
+
+        let c0 = a0_plus_a1
+            .mul(
+                cs.ns(|| "(a0 + a1) * (a0 + non_residue * a1)"),
+                &a0_plus_non_residue_c1,
+            )?
+            .sub(cs.ns(|| "- (1 + non_residue) v0"), &one_plus_non_residue_v0)?;
+
+        v0.double_in_place(cs.ns(|| "2v0"))?;
+        let c1 = v0;
+
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn inverse<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        let inverse = Self::alloc(&mut cs.ns(|| "alloc inverse"), || {
+            self.get_value().and_then(|val| val.inverse()).get()
+        })?;
+
+        // Karatsuba multiplication for Fp4 with the inverse:
+        //     v0 = A.c0 * B.c0
+        //     v1 = A.c1 * B.c1
+        //
+        //      1 = v0 + non_residue * v1
+        //  => v0 = 1 - non_residue * v1
+        //
+        //      0 = result.c1 = (A.c0 + A.c1) * (B.c0 + B.c1) - v0 - v1
+        //  => v0 + v1 = (A.c0 + A.c1) * (B.c0 + B.c1)
+        //  => 1 + (1 - non_residue) * v1 = (A.c0 + A.c1) * (B.c0 + B.c1)
+        // Enforced with 2 constraints:
+        //     A.c1 * B.c1 = v1
+        //  => 1 + (1 - non_residue) * v1 = (A.c0 + A.c1) * (B.c0 + B.c1)
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        // Devegili, OhEigeartaigh, Scott, Dahab
+
+        // Constraint 1
+        let v1 = self.c1.mul(cs.ns(|| "inv_constraint_1"), &inverse.c1)?;
+
+        // Constraint 2
+        let a0_plus_a1 = self.c0.add(cs.ns(|| "a0 + a1"), &self.c1)?;
+        let b0_plus_b1 = inverse.c0.add(cs.ns(|| "b0 + b1"), &inverse.c1)?;
+
+        let one = Fp2::<<P as Fp4Parameters>::Fp2Params>::one();
+        let rhs = Self::mul_fp2_gadget_by_nonresidue(cs.ns(|| "nr * v1"), &v1)?
+            .sub(cs.ns(|| "sub v1"), &v1)?
+            .negate(cs.ns(|| "negate it"))?
+            .add_constant(cs.ns(|| "add one"), &one)?;
+        a0_plus_a1.mul_equals(cs.ns(|| "inv_constraint_2"), &b0_plus_b1, &rhs)?;
+        Ok(inverse)
+    }
+
+    fn mul_equals<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+        result: &Self,
+    ) -> Result<(), SynthesisError> {
+        // Karatsuba multiplication for Fp4:
+        //     v0 = A.c0 * B.c0
+        //     v1 = A.c1 * B.c1
+        //     result.c0 = v0 + non_residue * v1
+        //     result.c1 = (A.c0 + A.c1) * (B.c0 + B.c1) - v0 - v1
+        // Enforced with 3 constraints:
+        //     A.c1 * B.c1 = v1
+        //     A.c0 * B.c0 = result.c0 - non_residue * v1
+        //     (A.c0+A.c1)*(B.c0+B.c1) = result.c1 + result.c0 + (1 - non_residue) * v1
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        // Devegili, OhEigeartaigh, Scott, Dahab
+        let mul_cs = &mut cs.ns(|| "mul");
+
+        // Compute v1
+        let mut v1 = self.c1.mul(mul_cs.ns(|| "v1"), &other.c1)?;
+
+        // Perform second check
+        let non_residue_times_v1 =
+            Self::mul_fp2_gadget_by_nonresidue(mul_cs.ns(|| "nr * v1"), &v1)?;
+        let rhs = result
+            .c0
+            .sub(mul_cs.ns(|| "sub from result.c0"), &non_residue_times_v1)?;
+        self.c0
+            .mul_equals(mul_cs.ns(|| "second check"), &other.c0, &rhs)?;
+
+        // Last check
+        let a0_plus_a1 = self.c0.add(mul_cs.ns(|| "a0 + a1"), &self.c1)?;
+        let b0_plus_b1 = other.c0.add(mul_cs.ns(|| "b0 + b1"), &other.c1)?;
+        let one_minus_non_residue_v1 =
+            v1.sub_in_place(mul_cs.ns(|| "sub from v1"), &non_residue_times_v1)?;
+
+        let result_c1_plus_result_c0_plus_one_minus_non_residue_v1 = result
+            .c1
+            .add(mul_cs.ns(|| "c1 + c0"), &result.c0)?
+            .add(mul_cs.ns(|| "rest of stuff"), one_minus_non_residue_v1)?;
+
+        a0_plus_a1.mul_equals(
+            mul_cs.ns(|| "third check"),
+            &b0_plus_b1,
+            &result_c1_plus_result_c0_plus_one_minus_non_residue_v1,
+        )?;
+
+        Ok(())
+    }
+
+    fn frobenius_map<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        power: usize,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        let _ = result.frobenius_map_in_place(cs, power)?;
+        Ok(result)
+    }
+
+    fn frobenius_map_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        power: usize,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0
+            .frobenius_map_in_place(cs.ns(|| "frob_map1"), power)?;
+        self.c1
+            .frobenius_map_in_place(cs.ns(|| "frob_map2"), power)?;
+        self.c1
+            .mul_by_fp_constant_in_place(cs.ns(|| "mul"), &P::FROBENIUS_COEFF_FP4_C1[power % 4])?;
+        Ok(self)
+    }
+
+    #[inline]
+    fn add_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        other: &Fp4<P>,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        let _ = result.add_constant_in_place(cs, other)?;
+        Ok(result)
+    }
+
+    #[inline]
+    fn add_constant_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        other: &Fp4<P>,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.add_constant_in_place(cs.ns(|| "c0"), &other.c0)?;
+        self.c1.add_constant_in_place(cs.ns(|| "c1"), &other.c1)?;
+        Ok(self)
+    }
+
+    fn mul_by_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        fe: &Fp4<P>,
+    ) -> Result<Self, SynthesisError> {
+        // Karatsuba multiplication (see mul above).
+        // Doesn't need any constraints; returns linear combinations of
+        // `self`'s variables.
+        //
+        // (The operations below are guaranteed to return linear combinations)
+        let (a0, a1) = (&self.c0, &self.c1);
+        let (b0, b1) = (fe.c0, fe.c1);
+        let mut v0 = a0.mul_by_constant(&mut cs.ns(|| "v0"), &b0)?;
+        let mut v1 = Self::mul_fp2_gadget_by_nonresidue(&mut cs.ns(|| "v1"), a1)?;
+        let beta_v1 = v1.mul_by_constant_in_place(&mut cs.ns(|| "beta * v1"), &b1)?;
+
+        v0.add_in_place(&mut cs.ns(|| "c0"), &beta_v1)?;
+        let c0 = v0;
+
+        let mut a0b1 = a0.mul_by_constant(&mut cs.ns(|| "a0b1"), &b1)?;
+        let a1b0 = a1.mul_by_constant(&mut cs.ns(|| "a1b0"), &b0)?;
+        a0b1.add_in_place(&mut cs.ns(|| "c1"), &a1b0)?;
+        let c1 = a0b1;
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost_of_mul() -> usize {
+        3 * <Fp2Gadget<P, ConstraintF> as FieldGadget<Fp2<P::Fp2Params>, ConstraintF>>::cost_of_mul(
+        )
+    }
+
+    fn cost_of_mul_equals() -> usize {
+        Self::cost_of_mul()
+    }
+
+    fn cost_of_inv() -> usize {
+        unimplemented!()
+    }
+}
+
+impl<P, ConstraintF: PrimeField> PartialEq for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    fn eq(&self, other: &Self) -> bool {
+        self.c0 == other.c0 && self.c1 == other.c1
+    }
+}
+
+impl<P, ConstraintF: PrimeField> Eq for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+}
+
+impl<P, ConstraintF: PrimeField> EqGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+}
+
+impl<P, ConstraintF: PrimeField> ConditionalEqGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn conditional_enforce_equal<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+        condition: &Boolean,
+    ) -> Result<(), SynthesisError> {
+        self.c0
+            .conditional_enforce_equal(&mut cs.ns(|| "c0"), &other.c0, condition)?;
+        self.c1
+            .conditional_enforce_equal(&mut cs.ns(|| "c1"), &other.c1, condition)?;
+        Ok(())
+    }
+
+    fn cost() -> usize {
+        2 * <Fp2Gadget<P, ConstraintF> as ConditionalEqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField> NEqGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn enforce_not_equal<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<(), SynthesisError> {
+        self.c0.enforce_not_equal(&mut cs.ns(|| "c0"), &other.c0)?;
+        self.c1.enforce_not_equal(&mut cs.ns(|| "c1"), &other.c1)?;
+        Ok(())
+    }
+
+    fn cost() -> usize {
+        2 * <Fp2Gadget<P, ConstraintF> as NEqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField> ToBitsGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    fn to_bits<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<Boolean>, SynthesisError> {
+        let mut c0 = self.c0.to_bits(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_bits(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+
+    fn to_non_unique_bits<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<Boolean>, SynthesisError> {
+        let mut c0 = self.c0.to_non_unique_bits(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_non_unique_bits(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+}
+
+impl<P, ConstraintF: PrimeField> ToBytesGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    fn to_bytes<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c0 = self.c0.to_bytes(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_bytes(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c0 = self.c0.to_non_unique_bytes(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_non_unique_bytes(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+}
+
+impl<P, ConstraintF: PrimeField> Clone for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    fn clone(&self) -> Self {
+        Self {
+            c0:      self.c0.clone(),
+            c1:      self.c1.clone(),
+            _params: PhantomData,
+        }
+    }
+}
+
+impl<P, ConstraintF: PrimeField> CondSelectGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn conditionally_select<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        cond: &Boolean,
+        first: &Self,
+        second: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = Fp2Gadget::<P, ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c0"),
+            cond,
+            &first.c0,
+            &second.c0,
+        )?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c1"),
+            cond,
+            &first.c1,
+            &second.c1,
+        )?;
+
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost() -> usize {
+        2 * <Fp2Gadget<P, ConstraintF> as CondSelectGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField> TwoBitLookupGadget<ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    type TableConstant = Fp4<P>;
+    fn two_bit_lookup<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        b: &[Boolean],
+        c: &[Self::TableConstant],
+    ) -> Result<Self, SynthesisError> {
+        let c0s = c.iter().map(|f| f.c0).collect::<Vec<_>>();
+        let c1s = c.iter().map(|f| f.c1).collect::<Vec<_>>();
+        let c0 = Fp2Gadget::<P, ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c0"), b, &c0s)?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c1"), b, &c1s)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost() -> usize {
+        2 * <Fp2Gadget<P, ConstraintF> as TwoBitLookupGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField> ThreeBitCondNegLookupGadget<ConstraintF>
+    for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    type TableConstant = Fp4<P>;
+
+    fn three_bit_cond_neg_lookup<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        b: &[Boolean],
+        b0b1: &Boolean,
+        c: &[Self::TableConstant],
+    ) -> Result<Self, SynthesisError> {
+        let c0s = c.iter().map(|f| f.c0).collect::<Vec<_>>();
+        let c1s = c.iter().map(|f| f.c1).collect::<Vec<_>>();
+        let c0 = Fp2Gadget::<P, ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c0"),
+            b,
+            b0b1,
+            &c0s,
+        )?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c1"),
+            b,
+            b0b1,
+            &c1s,
+        )?;
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost() -> usize {
+        2 * <Fp2Gadget<P, ConstraintF> as ThreeBitCondNegLookupGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField> AllocGadget<Fp4<P>, ConstraintF> for Fp4Gadget<P, ConstraintF>
+where
+    P: Fp4Parameters,
+    P::Fp2Params: Fp2Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn alloc<F, T, CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        value_gen: F,
+    ) -> Result<Self, SynthesisError>
+    where
+        F: FnOnce() -> Result<T, SynthesisError>,
+        T: Borrow<Fp4<P>>,
+    {
+        let (c0, c1) = match value_gen() {
+            Ok(fe) => {
+                let fe = *fe.borrow();
+                (Ok(fe.c0), Ok(fe.c1))
+            },
+            Err(_) => (
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+            ),
+        };
+
+        let c0 = Fp2Gadget::<P, ConstraintF>::alloc(&mut cs.ns(|| "c0"), || c0)?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::alloc(&mut cs.ns(|| "c1"), || c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn alloc_input<F, T, CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        value_gen: F,
+    ) -> Result<Self, SynthesisError>
+    where
+        F: FnOnce() -> Result<T, SynthesisError>,
+        T: Borrow<Fp4<P>>,
+    {
+        let (c0, c1) = match value_gen() {
+            Ok(fe) => {
+                let fe = *fe.borrow();
+                (Ok(fe.c0), Ok(fe.c1))
+            },
+            Err(_) => (
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+            ),
+        };
+
+        let c0 = Fp2Gadget::<P, ConstraintF>::alloc_input(&mut cs.ns(|| "c0"), || c0)?;
+        let c1 = Fp2Gadget::<P, ConstraintF>::alloc_input(&mut cs.ns(|| "c1"), || c1)?;
+        Ok(Self::new(c0, c1))
+    }
+}

r1cs-std/src/fields/fp6_2over3.rs

@@ -0,0 +1,763 @@
+use algebra::{
+    fields::{
+        fp6_2over3::{Fp6, Fp6Parameters},
+        Fp3, Fp3Parameters,
+    },
+    BigInteger, Field, One, PrimeField, SquareRootField,
+};
+use core::{borrow::Borrow, marker::PhantomData};
+use r1cs_core::{ConstraintSystem, SynthesisError};
+
+use crate::{prelude::*, Assignment, Vec};
+
+type Fp3Gadget<P, ConstraintF> =
+    super::fp3::Fp3Gadget<<P as Fp6Parameters>::Fp3Params, ConstraintF>;
+type Fp3GadgetVariable<P, ConstraintF> = <Fp3Gadget<P, ConstraintF> as FieldGadget<
+    Fp3<<P as Fp6Parameters>::Fp3Params>,
+    ConstraintF,
+>>::Variable;
+
+#[derive(Derivative)]
+#[derivative(Debug(bound = "ConstraintF: PrimeField + SquareRootField"))]
+#[must_use]
+pub struct Fp6Gadget<P, ConstraintF: PrimeField + SquareRootField>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    pub c0:  Fp3Gadget<P, ConstraintF>,
+    pub c1:  Fp3Gadget<P, ConstraintF>,
+    #[derivative(Debug = "ignore")]
+    _params: PhantomData<P>,
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    pub fn new(c0: Fp3Gadget<P, ConstraintF>, c1: Fp3Gadget<P, ConstraintF>) -> Self {
+        Self {
+            c0,
+            c1,
+            _params: PhantomData,
+        }
+    }
+
+    /// Multiply a Fp3Gadget by quadratic nonresidue P::NONRESIDUE.
+    #[inline]
+    pub fn mul_fp3_gadget_by_nonresidue<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        fe: &Fp3Gadget<P, ConstraintF>,
+    ) -> Result<Fp3Gadget<P, ConstraintF>, SynthesisError> {
+        let mut res = Fp3Gadget::<P, ConstraintF>::new(fe.c2.clone(), fe.c0.clone(), fe.c1.clone());
+        res.c0.mul_by_constant_in_place(
+            cs.ns(|| "res * non_residue"),
+            &<P::Fp3Params as Fp3Parameters>::NONRESIDUE,
+        )?;
+        Ok(res)
+    }
+
+    pub fn unitary_inverse<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        Ok(Self::new(self.c0.clone(), self.c1.negate(cs)?))
+    }
+
+    #[inline]
+    pub fn cyclotomic_exp<CS: ConstraintSystem<ConstraintF>, B: BigInteger>(
+        &self,
+        mut cs: CS,
+        exponent: &B,
+    ) -> Result<Self, SynthesisError> {
+        let mut res = Self::one(cs.ns(|| "one"))?;
+        let self_inverse = self.unitary_inverse(cs.ns(|| "unitary inverse"))?;
+
+        let mut found_nonzero = false;
+        let naf = exponent.find_wnaf();
+
+        for (i, &value) in naf.iter().rev().enumerate() {
+            if found_nonzero {
+                res.square_in_place(cs.ns(|| format!("square {}", i)))?;
+            }
+
+            if value != 0 {
+                found_nonzero = true;
+
+                if value > 0 {
+                    res.mul_in_place(cs.ns(|| format!("res *= self {}", i)), &self)?;
+                } else {
+                    res.mul_in_place(
+                        cs.ns(|| format!("res *= self_inverse {}", i)),
+                        &self_inverse,
+                    )?;
+                }
+            }
+        }
+
+        Ok(res)
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> FieldGadget<Fp6<P>, ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    type Variable = (
+        Fp3GadgetVariable<P, ConstraintF>,
+        Fp3GadgetVariable<P, ConstraintF>,
+    );
+
+    #[inline]
+    fn get_value(&self) -> Option<Fp6<P>> {
+        match (self.c0.get_value(), self.c1.get_value()) {
+            (Some(c0), Some(c1)) => Some(Fp6::new(c0, c1)),
+            (..) => None,
+        }
+    }
+
+    #[inline]
+    fn get_variable(&self) -> Self::Variable {
+        (self.c0.get_variable(), self.c1.get_variable())
+    }
+
+    #[inline]
+    fn zero<CS: ConstraintSystem<ConstraintF>>(mut cs: CS) -> Result<Self, SynthesisError> {
+        let c0 = Fp3Gadget::<P, ConstraintF>::zero(cs.ns(|| "c0"))?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::zero(cs.ns(|| "c1"))?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn one<CS: ConstraintSystem<ConstraintF>>(mut cs: CS) -> Result<Self, SynthesisError> {
+        let c0 = Fp3Gadget::<P, ConstraintF>::one(cs.ns(|| "c0"))?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::zero(cs.ns(|| "c1"))?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn conditionally_add_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        bit: &Boolean,
+        coeff: Fp6<P>,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self
+            .c0
+            .conditionally_add_constant(cs.ns(|| "c0"), bit, coeff.c0)?;
+        let c1 = self
+            .c1
+            .conditionally_add_constant(cs.ns(|| "c1"), bit, coeff.c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn add<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.add(&mut cs.ns(|| "add c0"), &other.c0)?;
+        let c1 = self.c1.add(&mut cs.ns(|| "add c1"), &other.c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn sub<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = self.c0.sub(&mut cs.ns(|| "sub c0"), &other.c0)?;
+        let c1 = self.c1.sub(&mut cs.ns(|| "sub c1"), &other.c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn double<CS: ConstraintSystem<ConstraintF>>(&self, cs: CS) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.double_in_place(cs)?;
+        Ok(result)
+    }
+
+    #[inline]
+    fn double_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.double_in_place(&mut cs.ns(|| "double c0"))?;
+        self.c1.double_in_place(&mut cs.ns(|| "double c1"))?;
+        Ok(self)
+    }
+
+    #[inline]
+    fn negate<CS: ConstraintSystem<ConstraintF>>(&self, cs: CS) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        result.negate_in_place(cs)?;
+        Ok(result)
+    }
+
+    #[inline]
+    fn negate_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.negate_in_place(&mut cs.ns(|| "negate c0"))?;
+        self.c1.negate_in_place(&mut cs.ns(|| "negate c1"))?;
+        Ok(self)
+    }
+
+    #[inline]
+    fn mul<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<Self, SynthesisError> {
+        // Karatsuba multiplication for Fp6:
+        //     v0 = A.c0 * B.c0
+        //     v1 = A.c1 * B.c1
+        //     result.c0 = v0 + non_residue * v1
+        //     result.c1 = (A.c0 + A.c1) * (B.c0 + B.c1) - v0 - v1
+        // Enforced with 3 constraints:
+        //     A.c1 * B.c1 = v1
+        //     A.c0 * B.c0 = result.c0 - non_residue * v1
+        //     (A.c0+A.c1)*(B.c0+B.c1) = result.c1 + result.c0 + (1 - non_residue) * v1
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        // Devegili, OhEigeartaigh, Scott, Dahab
+        let mul_cs = &mut cs.ns(|| "mul");
+
+        let v0 = self.c0.mul(mul_cs.ns(|| "v0"), &other.c0)?;
+        let v1 = self.c1.mul(mul_cs.ns(|| "v1"), &other.c1)?;
+        let c0 = {
+            let non_residue_times_v1 =
+                Self::mul_fp3_gadget_by_nonresidue(mul_cs.ns(|| "first mul_by_nr"), &v1)?;
+            v0.add(mul_cs.ns(|| "v0 + beta * v1"), &non_residue_times_v1)?
+        };
+        let c1 = {
+            let a0_plus_a1 = self.c0.add(mul_cs.ns(|| "a0 + a1"), &self.c1)?;
+            let b0_plus_b1 = other.c0.add(mul_cs.ns(|| "b0 + b1"), &other.c1)?;
+            let a0_plus_a1_times_b0_plus_b1 =
+                a0_plus_a1.mul(&mut mul_cs.ns(|| "(a0 + a1) * (b0 + b1)"), &b0_plus_b1)?;
+            a0_plus_a1_times_b0_plus_b1
+                .sub(mul_cs.ns(|| "res - v0"), &v0)?
+                .sub(mul_cs.ns(|| "res - v0 - v1"), &v1)?
+        };
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn square<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        // From Libsnark/fp4_gadget.tcc
+        // Complex multiplication for Fp6:
+        //     v0 = A.c0 * A.c1
+        //     result.c0 = (A.c0 + A.c1) * (A.c0 + non_residue * A.c1) - (1 +
+        // non_residue) * v0     result.c1 = 2 * v0
+        // Enforced with 2 constraints:
+        //     (2*A.c0) * A.c1 = result.c1
+        //     (A.c0 + A.c1) * (A.c0 + non_residue * A.c1) = result.c0 + result.c1 * (1
+        // + non_residue)/2 Reference:
+        //     "Multiplication and Squaring on Pairing-Friendly Fields"
+        //     Devegili, OhEigeartaigh, Scott, Dahab
+
+        let mut v0 = self.c0.mul(cs.ns(|| "v0"), &self.c1)?;
+        let a0_plus_a1 = self.c0.add(cs.ns(|| "a0 + a1"), &self.c1)?;
+
+        let non_residue_c1 =
+            Self::mul_fp3_gadget_by_nonresidue(cs.ns(|| "non_residue * a1"), &self.c1)?;
+        let a0_plus_non_residue_c1 = self
+            .c0
+            .add(cs.ns(|| "a0 + non_residue * a1"), &non_residue_c1)?;
+        let one_plus_non_residue_v0 =
+            Self::mul_fp3_gadget_by_nonresidue(cs.ns(|| "non_residue * v0"), &v0)?
+                .add(cs.ns(|| "plus v0"), &v0)?;
+
+        let c0 = a0_plus_a1
+            .mul(
+                cs.ns(|| "(a0 + a1) * (a0 + non_residue * a1)"),
+                &a0_plus_non_residue_c1,
+            )?
+            .sub(cs.ns(|| "- (1 + non_residue) v0"), &one_plus_non_residue_v0)?;
+
+        v0.double_in_place(cs.ns(|| "2v0"))?;
+        let c1 = v0;
+
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn inverse<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Self, SynthesisError> {
+        let inverse = Self::alloc(&mut cs.ns(|| "alloc inverse"), || {
+            self.get_value().and_then(|val| val.inverse()).get()
+        })?;
+
+        // Karatsuba multiplication for Fp6 with the inverse:
+        //     v0 = A.c0 * B.c0
+        //     v1 = A.c1 * B.c1
+        //
+        //      1 = v0 + non_residue * v1
+        //  => v0 = 1 - non_residue * v1
+        //
+        //      0 = result.c1 = (A.c0 + A.c1) * (B.c0 + B.c1) - v0 - v1
+        //  => v0 + v1 = (A.c0 + A.c1) * (B.c0 + B.c1)
+        //  => 1 + (1 - non_residue) * v1 = (A.c0 + A.c1) * (B.c0 + B.c1)
+        // Enforced with 2 constraints:
+        //     A.c1 * B.c1 = v1
+        //  => 1 + (1 - non_residue) * v1 = (A.c0 + A.c1) * (B.c0 + B.c1)
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        // Devegili, OhEigeartaigh, Scott, Dahab
+
+        // Constraint 1
+        let v1 = self.c1.mul(cs.ns(|| "inv_constraint_1"), &inverse.c1)?;
+
+        // Constraint 2
+        let a0_plus_a1 = self.c0.add(cs.ns(|| "a0 + a1"), &self.c1)?;
+        let b0_plus_b1 = inverse.c0.add(cs.ns(|| "b0 + b1"), &inverse.c1)?;
+
+        let one = Fp3::<<P as Fp6Parameters>::Fp3Params>::one();
+        let rhs = Self::mul_fp3_gadget_by_nonresidue(cs.ns(|| "nr * v1"), &v1)?
+            .sub(cs.ns(|| "sub v1"), &v1)?
+            .negate(cs.ns(|| "negate it"))?
+            .add_constant(cs.ns(|| "add one"), &one)?;
+        a0_plus_a1.mul_equals(cs.ns(|| "inv_constraint_2"), &b0_plus_b1, &rhs)?;
+        Ok(inverse)
+    }
+
+    fn mul_equals<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+        result: &Self,
+    ) -> Result<(), SynthesisError> {
+        // Karatsuba multiplication for Fp6:
+        //     v0 = A.c0 * B.c0
+        //     v1 = A.c1 * B.c1
+        //     result.c0 = v0 + non_residue * v1
+        //     result.c1 = (A.c0 + A.c1) * (B.c0 + B.c1) - v0 - v1
+        // Enforced with 3 constraints:
+        //     A.c1 * B.c1 = v1
+        //     A.c0 * B.c0 = result.c0 - non_residue * v1
+        //     (A.c0+A.c1)*(B.c0+B.c1) = result.c1 + result.c0 + (1 - non_residue) * v1
+        // Reference:
+        // "Multiplication and Squaring on Pairing-Friendly Fields"
+        // Devegili, OhEigeartaigh, Scott, Dahab
+        let mul_cs = &mut cs.ns(|| "mul");
+
+        // Compute v1
+        let mut v1 = self.c1.mul(mul_cs.ns(|| "v1"), &other.c1)?;
+
+        // Perform second check
+        let non_residue_times_v1 =
+            Self::mul_fp3_gadget_by_nonresidue(mul_cs.ns(|| "nr * v1"), &v1)?;
+        let rhs = result
+            .c0
+            .sub(mul_cs.ns(|| "sub from result.c0"), &non_residue_times_v1)?;
+        self.c0
+            .mul_equals(mul_cs.ns(|| "second check"), &other.c0, &rhs)?;
+
+        // Last check
+        let a0_plus_a1 = self.c0.add(mul_cs.ns(|| "a0 + a1"), &self.c1)?;
+        let b0_plus_b1 = other.c0.add(mul_cs.ns(|| "b0 + b1"), &other.c1)?;
+        let one_minus_non_residue_v1 =
+            v1.sub_in_place(mul_cs.ns(|| "sub from v1"), &non_residue_times_v1)?;
+
+        let result_c1_plus_result_c0_plus_one_minus_non_residue_v1 = result
+            .c1
+            .add(mul_cs.ns(|| "c1 + c0"), &result.c0)?
+            .add(mul_cs.ns(|| "rest of stuff"), one_minus_non_residue_v1)?;
+
+        a0_plus_a1.mul_equals(
+            mul_cs.ns(|| "third check"),
+            &b0_plus_b1,
+            &result_c1_plus_result_c0_plus_one_minus_non_residue_v1,
+        )?;
+
+        Ok(())
+    }
+
+    fn frobenius_map<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        power: usize,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        let _ = result.frobenius_map_in_place(cs, power)?;
+        Ok(result)
+    }
+
+    fn frobenius_map_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        power: usize,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0
+            .frobenius_map_in_place(cs.ns(|| "frob_map1"), power)?;
+        self.c1
+            .frobenius_map_in_place(cs.ns(|| "frob_map2"), power)?;
+        self.c1
+            .mul_by_fp_constant_in_place(cs.ns(|| "mul"), &P::FROBENIUS_COEFF_FP6_C1[power % 6])?;
+        Ok(self)
+    }
+
+    #[inline]
+    fn add_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        cs: CS,
+        other: &Fp6<P>,
+    ) -> Result<Self, SynthesisError> {
+        let mut result = self.clone();
+        let _ = result.add_constant_in_place(cs, other)?;
+        Ok(result)
+    }
+
+    #[inline]
+    fn add_constant_in_place<CS: ConstraintSystem<ConstraintF>>(
+        &mut self,
+        mut cs: CS,
+        other: &Fp6<P>,
+    ) -> Result<&mut Self, SynthesisError> {
+        self.c0.add_constant_in_place(cs.ns(|| "c0"), &other.c0)?;
+        self.c1.add_constant_in_place(cs.ns(|| "c1"), &other.c1)?;
+        Ok(self)
+    }
+
+    fn mul_by_constant<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        fe: &Fp6<P>,
+    ) -> Result<Self, SynthesisError> {
+        // Karatsuba multiplication (see mul above).
+        // Doesn't need any constraints; returns linear combinations of
+        // `self`'s variables.
+        //
+        // (The operations below are guaranteed to return linear combinations)
+        let (a0, a1) = (&self.c0, &self.c1);
+        let (b0, b1) = (fe.c0, fe.c1);
+        let mut v0 = a0.mul_by_constant(&mut cs.ns(|| "v0"), &b0)?;
+        let mut v1 = Self::mul_fp3_gadget_by_nonresidue(&mut cs.ns(|| "v1"), a1)?;
+        let beta_v1 = v1.mul_by_constant_in_place(&mut cs.ns(|| "beta * v1"), &b1)?;
+
+        v0.add_in_place(&mut cs.ns(|| "c0"), &beta_v1)?;
+        let c0 = v0;
+
+        let mut a0b1 = a0.mul_by_constant(&mut cs.ns(|| "a0b1"), &b1)?;
+        let a1b0 = a1.mul_by_constant(&mut cs.ns(|| "a1b0"), &b0)?;
+        a0b1.add_in_place(&mut cs.ns(|| "c1"), &a1b0)?;
+        let c1 = a0b1;
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost_of_mul() -> usize {
+        3 * <Fp3Gadget<P, ConstraintF> as FieldGadget<Fp3<P::Fp3Params>, ConstraintF>>::cost_of_mul(
+        )
+    }
+
+    fn cost_of_mul_equals() -> usize {
+        Self::cost_of_mul()
+    }
+
+    fn cost_of_inv() -> usize {
+        unimplemented!()
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> PartialEq for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    fn eq(&self, other: &Self) -> bool {
+        self.c0 == other.c0 && self.c1 == other.c1
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> Eq for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> EqGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> ConditionalEqGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn conditional_enforce_equal<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+        condition: &Boolean,
+    ) -> Result<(), SynthesisError> {
+        self.c0
+            .conditional_enforce_equal(&mut cs.ns(|| "c0"), &other.c0, condition)?;
+        self.c1
+            .conditional_enforce_equal(&mut cs.ns(|| "c1"), &other.c1, condition)?;
+        Ok(())
+    }
+
+    fn cost() -> usize {
+        2 * <Fp3Gadget<P, ConstraintF> as ConditionalEqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> NEqGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn enforce_not_equal<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+        other: &Self,
+    ) -> Result<(), SynthesisError> {
+        self.c0.enforce_not_equal(&mut cs.ns(|| "c0"), &other.c0)?;
+        self.c1.enforce_not_equal(&mut cs.ns(|| "c1"), &other.c1)?;
+        Ok(())
+    }
+
+    fn cost() -> usize {
+        2 * <Fp3Gadget<P, ConstraintF> as NEqGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> ToBitsGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    fn to_bits<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<Boolean>, SynthesisError> {
+        let mut c0 = self.c0.to_bits(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_bits(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+
+    fn to_non_unique_bits<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<Boolean>, SynthesisError> {
+        let mut c0 = self.c0.to_non_unique_bits(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_non_unique_bits(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> ToBytesGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    fn to_bytes<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c0 = self.c0.to_bytes(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_bytes(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<ConstraintF>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c0 = self.c0.to_non_unique_bytes(cs.ns(|| "c0"))?;
+        let mut c1 = self.c1.to_non_unique_bytes(cs.ns(|| "c1"))?;
+        c0.append(&mut c1);
+        Ok(c0)
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> Clone for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    fn clone(&self) -> Self {
+        Self {
+            c0:      self.c0.clone(),
+            c1:      self.c1.clone(),
+            _params: PhantomData,
+        }
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> CondSelectGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn conditionally_select<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        cond: &Boolean,
+        first: &Self,
+        second: &Self,
+    ) -> Result<Self, SynthesisError> {
+        let c0 = Fp3Gadget::<P, ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c0"),
+            cond,
+            &first.c0,
+            &second.c0,
+        )?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::conditionally_select(
+            &mut cs.ns(|| "c1"),
+            cond,
+            &first.c1,
+            &second.c1,
+        )?;
+
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost() -> usize {
+        2 * <Fp3Gadget<P, ConstraintF> as CondSelectGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> TwoBitLookupGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    type TableConstant = Fp6<P>;
+    fn two_bit_lookup<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        b: &[Boolean],
+        c: &[Self::TableConstant],
+    ) -> Result<Self, SynthesisError> {
+        let c0s = c.iter().map(|f| f.c0).collect::<Vec<_>>();
+        let c1s = c.iter().map(|f| f.c1).collect::<Vec<_>>();
+        let c0 = Fp3Gadget::<P, ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c0"), b, &c0s)?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::two_bit_lookup(cs.ns(|| "Lookup c1"), b, &c1s)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost() -> usize {
+        2 * <Fp3Gadget<P, ConstraintF> as TwoBitLookupGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> ThreeBitCondNegLookupGadget<ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    type TableConstant = Fp6<P>;
+
+    fn three_bit_cond_neg_lookup<CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        b: &[Boolean],
+        b0b1: &Boolean,
+        c: &[Self::TableConstant],
+    ) -> Result<Self, SynthesisError> {
+        let c0s = c.iter().map(|f| f.c0).collect::<Vec<_>>();
+        let c1s = c.iter().map(|f| f.c1).collect::<Vec<_>>();
+        let c0 = Fp3Gadget::<P, ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c0"),
+            b,
+            b0b1,
+            &c0s,
+        )?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::three_bit_cond_neg_lookup(
+            cs.ns(|| "Lookup c1"),
+            b,
+            b0b1,
+            &c1s,
+        )?;
+        Ok(Self::new(c0, c1))
+    }
+
+    fn cost() -> usize {
+        2 * <Fp3Gadget<P, ConstraintF> as ThreeBitCondNegLookupGadget<ConstraintF>>::cost()
+    }
+}
+
+impl<P, ConstraintF: PrimeField + SquareRootField> AllocGadget<Fp6<P>, ConstraintF>
+    for Fp6Gadget<P, ConstraintF>
+where
+    P: Fp6Parameters,
+    P::Fp3Params: Fp3Parameters<Fp = ConstraintF>,
+{
+    #[inline]
+    fn alloc<F, T, CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        value_gen: F,
+    ) -> Result<Self, SynthesisError>
+    where
+        F: FnOnce() -> Result<T, SynthesisError>,
+        T: Borrow<Fp6<P>>,
+    {
+        let (c0, c1) = match value_gen() {
+            Ok(fe) => {
+                let fe = *fe.borrow();
+                (Ok(fe.c0), Ok(fe.c1))
+            },
+            Err(_) => (
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+            ),
+        };
+
+        let c0 = Fp3Gadget::<P, ConstraintF>::alloc(&mut cs.ns(|| "c0"), || c0)?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::alloc(&mut cs.ns(|| "c1"), || c1)?;
+        Ok(Self::new(c0, c1))
+    }
+
+    #[inline]
+    fn alloc_input<F, T, CS: ConstraintSystem<ConstraintF>>(
+        mut cs: CS,
+        value_gen: F,
+    ) -> Result<Self, SynthesisError>
+    where
+        F: FnOnce() -> Result<T, SynthesisError>,
+        T: Borrow<Fp6<P>>,
+    {
+        let (c0, c1) = match value_gen() {
+            Ok(fe) => {
+                let fe = *fe.borrow();
+                (Ok(fe.c0), Ok(fe.c1))
+            },
+            Err(_) => (
+                Err(SynthesisError::AssignmentMissing),
+                Err(SynthesisError::AssignmentMissing),
+            ),
+        };
+
+        let c0 = Fp3Gadget::<P, ConstraintF>::alloc_input(&mut cs.ns(|| "c0"), || c0)?;
+        let c1 = Fp3Gadget::<P, ConstraintF>::alloc_input(&mut cs.ns(|| "c1"), || c1)?;
+        Ok(Self::new(c0, c1))
+    }
+}

r1cs-std/src/fields/fp6_3over2.rs

@@ -259,7 +259,7 @@ where
         mut cs: CS,
         other: &Self,
     ) -> Result<Self, SynthesisError> {
-        // Uses Toom-Cool-3x multiplication from
+        // Uses Toom-Cook-3x multiplication from
         //
         // Reference:
         // "Multiplication and Squaring on Pairing-Friendly Fields"
@@ -404,7 +404,7 @@ where
         &self,
         mut cs: CS,
     ) -> Result<Self, SynthesisError> {
-        // Uses Toom-Cool-3x multiplication from
+        // Uses Toom-Cook-3x multiplication from
         //
         // Reference:
         // "Multiplication and Squaring on Pairing-Friendly Fields"
@@ -500,7 +500,7 @@ where
         Ok(Self::new(c0, c1, c2))
     }
 
-    // 18 constaints, we can probably do better but not sure it's worth it.
+    // 18 constraints, we can probably do better but not sure it's worth it.
     #[inline]
     fn inverse<CS: ConstraintSystem<ConstraintF>>(
         &self,

r1cs-std/src/fields/mod.rs

@@ -7,6 +7,9 @@ use crate::prelude::*;
 pub mod fp;
 pub mod fp12;
 pub mod fp2;
+pub mod fp3;
+pub mod fp4;
+pub mod fp6_2over3;
 pub mod fp6_3over2;
 
 pub trait FieldGadget<F: Field, ConstraintF: Field>:
@@ -248,6 +251,10 @@ pub trait FieldGadget<F: Field, ConstraintF: Field>:
 
     fn cost_of_mul() -> usize;
 
+    fn cost_of_mul_equals() -> usize {
+        Self::cost_of_mul() + <Self as EqGadget<ConstraintF>>::cost()
+    }
+
     fn cost_of_inv() -> usize;
 }
 

r1cs-std/src/groups/curves/short_weierstrass/mnt4/mod.rs

@@ -0,0 +1,421 @@
+use algebra::{
+    curves::mnt4::{
+        g2::{AteAdditionCoefficients, AteDoubleCoefficients},
+        G1Prepared, G2Prepared, MNT4Parameters,
+    },
+    Field,
+};
+use r1cs_core::{ConstraintSystem, SynthesisError};
+
+use crate::{
+    fields::{fp::FpGadget, fp2::Fp2Gadget, FieldGadget},
+    groups::curves::short_weierstrass::AffineGadget,
+    pairing::mnt4::PairingGadget,
+    prelude::*,
+    Vec,
+};
+
+pub type G1Gadget<P> = AffineGadget<
+    <P as MNT4Parameters>::G1Parameters,
+    <P as MNT4Parameters>::Fp,
+    FpGadget<<P as MNT4Parameters>::Fp>,
+>;
+
+pub type G2Gadget<P> =
+    AffineGadget<<P as MNT4Parameters>::G2Parameters, <P as MNT4Parameters>::Fp, Fp2G<P>>;
+
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT4Parameters"), Debug(bound = "P: MNT4Parameters"))]
+pub struct G1PreparedGadget<P: MNT4Parameters> {
+    pub x:       FpGadget<P::Fp>,
+    pub y:       FpGadget<P::Fp>,
+    pub x_twist: Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub y_twist: Fp2Gadget<P::Fp2Params, P::Fp>,
+}
+
+impl<P: MNT4Parameters> G1PreparedGadget<P> {
+    pub fn get_value(&self) -> Option<G1Prepared<P>> {
+        match (
+            self.x.get_value(),
+            self.y.get_value(),
+            self.x_twist.get_value(),
+            self.y_twist.get_value(),
+        ) {
+            (Some(x), Some(y), Some(x_twist), Some(y_twist)) => Some(G1Prepared {
+                x,
+                y,
+                x_twist,
+                y_twist,
+            }),
+            _ => None,
+        }
+    }
+
+    pub fn from_affine<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        q: &G1Gadget<P>,
+    ) -> Result<Self, SynthesisError> {
+        let x_twist = Fp2Gadget::new(
+            q.x.mul_by_constant(cs.ns(|| "g1.x * twist.c0"), &P::TWIST.c0)?,
+            q.x.mul_by_constant(cs.ns(|| "g1.x * twist.c1"), &P::TWIST.c1)?,
+        );
+        let y_twist = Fp2Gadget::new(
+            q.y.mul_by_constant(cs.ns(|| "g1.y * twist.c0"), &P::TWIST.c0)?,
+            q.y.mul_by_constant(cs.ns(|| "g1.y * twist.c1"), &P::TWIST.c1)?,
+        );
+        Ok(G1PreparedGadget {
+            x: q.x.clone(),
+            y: q.y.clone(),
+            x_twist,
+            y_twist,
+        })
+    }
+}
+
+impl<P: MNT4Parameters> ToBytesGadget<P::Fp> for G1PreparedGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_twist = self.x_twist.to_bytes(&mut cs.ns(|| "x_twist to bytes"))?;
+        let mut y_twist = self.y_twist.to_bytes(&mut cs.ns(|| "y_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_twist);
+        x.append(&mut y_twist);
+        Ok(x)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_non_unique_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_non_unique_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_twist = self
+            .x_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "x_twist to bytes"))?;
+        let mut y_twist = self
+            .y_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "y_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_twist);
+        x.append(&mut y_twist);
+        Ok(x)
+    }
+}
+
+type Fp2G<P> = Fp2Gadget<<P as MNT4Parameters>::Fp2Params, <P as MNT4Parameters>::Fp>;
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT4Parameters"), Debug(bound = "P: MNT4Parameters"))]
+pub struct G2PreparedGadget<P: MNT4Parameters> {
+    pub x:                     Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub y:                     Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub x_over_twist:          Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub y_over_twist:          Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub double_coefficients:   Vec<AteDoubleCoefficientsGadget<P>>,
+    pub addition_coefficients: Vec<AteAdditionCoefficientsGadget<P>>,
+}
+
+impl<P: MNT4Parameters> ToBytesGadget<P::Fp> for G2PreparedGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_over_twist = self
+            .x_over_twist
+            .to_bytes(&mut cs.ns(|| "x_over_twist to bytes"))?;
+        let mut y_over_twist = self
+            .y_over_twist
+            .to_bytes(&mut cs.ns(|| "y_over_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_over_twist);
+        x.append(&mut y_over_twist);
+
+        for (i, coeff) in self.double_coefficients.iter().enumerate() {
+            x.extend_from_slice(&coeff.to_bytes(cs.ns(|| format!("double_coefficients {}", i)))?);
+        }
+        for (i, coeff) in self.addition_coefficients.iter().enumerate() {
+            x.extend_from_slice(&coeff.to_bytes(cs.ns(|| format!("addition_coefficients {}", i)))?);
+        }
+        Ok(x)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_non_unique_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_non_unique_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_over_twist = self
+            .x_over_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "x_over_twist to bytes"))?;
+        let mut y_over_twist = self
+            .y_over_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "y_over_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_over_twist);
+        x.append(&mut y_over_twist);
+
+        for (i, coeff) in self.double_coefficients.iter().enumerate() {
+            x.extend_from_slice(
+                &coeff.to_non_unique_bytes(cs.ns(|| format!("double_coefficients {}", i)))?,
+            );
+        }
+        for (i, coeff) in self.addition_coefficients.iter().enumerate() {
+            x.extend_from_slice(
+                &coeff.to_non_unique_bytes(cs.ns(|| format!("addition_coefficients {}", i)))?,
+            );
+        }
+        Ok(x)
+    }
+}
+
+impl<P: MNT4Parameters> G2PreparedGadget<P> {
+    pub fn get_value(&self) -> Option<G2Prepared<P>> {
+        match (
+            self.x.get_value(),
+            self.y.get_value(),
+            self.x_over_twist.get_value(),
+            self.y_over_twist.get_value(),
+            self.double_coefficients
+                .iter()
+                .map(|coeff| coeff.get_value())
+                .collect::<Option<Vec<AteDoubleCoefficients<P>>>>(),
+            self.addition_coefficients
+                .iter()
+                .map(|coeff| coeff.get_value())
+                .collect::<Option<Vec<AteAdditionCoefficients<P>>>>(),
+        ) {
+            (
+                Some(x),
+                Some(y),
+                Some(x_over_twist),
+                Some(y_over_twist),
+                Some(double_coefficients),
+                Some(addition_coefficients),
+            ) => Some(G2Prepared {
+                x,
+                y,
+                x_over_twist,
+                y_over_twist,
+                double_coefficients,
+                addition_coefficients,
+            }),
+            _ => None,
+        }
+    }
+
+    pub fn from_affine<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        q: &G2Gadget<P>,
+    ) -> Result<Self, SynthesisError> {
+        let twist_inv = P::TWIST.inverse().unwrap();
+
+        let mut g2p = G2PreparedGadget {
+            x:                     q.x.clone(),
+            y:                     q.y.clone(),
+            x_over_twist:          q.x.mul_by_constant(cs.ns(|| "x over twist"), &twist_inv)?,
+            y_over_twist:          q.y.mul_by_constant(cs.ns(|| "y over twist"), &twist_inv)?,
+            double_coefficients:   vec![],
+            addition_coefficients: vec![],
+        };
+
+        let fp2_one = Fp2G::<P>::one(cs.ns(|| "one"))?;
+        let mut r = G2ProjectiveExtendedGadget {
+            x: q.x.clone(),
+            y: q.y.clone(),
+            z: fp2_one.clone(),
+            t: fp2_one,
+        };
+
+        for (idx, value) in P::ATE_LOOP_COUNT.iter().rev().enumerate() {
+            let mut tmp = *value;
+            let skip_extraneous_bits = 64 - value.leading_zeros();
+            let mut v = Vec::with_capacity(16);
+            for i in 0..64 {
+                if idx == 0 && (i == 0 || i >= skip_extraneous_bits) {
+                    continue;
+                }
+                v.push(tmp & 1 == 1);
+                tmp >>= 1;
+            }
+
+            let mut cs = cs.ns(|| format!("ate loop iteration {}", idx));
+
+            for (j, bit) in v.iter().rev().enumerate() {
+                let (r2, coeff) = PairingGadget::<P>::doubling_step_for_flipped_miller_loop(
+                    cs.ns(|| format!("doubling step {}", j)),
+                    &r,
+                )?;
+                g2p.double_coefficients.push(coeff);
+                r = r2;
+
+                if *bit {
+                    let (r2, coeff) =
+                        PairingGadget::<P>::mixed_addition_step_for_flipped_miller_loop(
+                            cs.ns(|| format!("mixed addition step {}", j)),
+                            &q.x,
+                            &q.y,
+                            &r,
+                        )?;
+                    g2p.addition_coefficients.push(coeff);
+                    r = r2;
+                }
+
+                tmp >>= 1;
+            }
+        }
+
+        if P::ATE_IS_LOOP_COUNT_NEG {
+            let rz_inv = r.z.inverse(cs.ns(|| "inverse r.z"))?;
+            let rz2_inv = rz_inv.square(cs.ns(|| "rz_inv^2"))?;
+            let rz3_inv = rz_inv.mul(cs.ns(|| "rz_inv * rz_inv^2"), &rz2_inv)?;
+
+            let minus_r_affine_x = r.x.mul(cs.ns(|| "r.x * rz2_inv"), &rz2_inv)?;
+            let minus_r_affine_y =
+                r.y.negate(cs.ns(|| "-r.y"))?
+                    .mul(cs.ns(|| "-r.y * rz3_inv"), &rz3_inv)?;
+
+            let add_result = PairingGadget::<P>::mixed_addition_step_for_flipped_miller_loop(
+                cs.ns(|| "mixed_addition step"),
+                &minus_r_affine_x,
+                &minus_r_affine_y,
+                &r,
+            )?;
+            g2p.addition_coefficients.push(add_result.1);
+        }
+
+        Ok(g2p)
+    }
+}
+
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT4Parameters"), Debug(bound = "P: MNT4Parameters"))]
+pub struct AteDoubleCoefficientsGadget<P: MNT4Parameters> {
+    pub c_h:  Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub c_4c: Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub c_j:  Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub c_l:  Fp2Gadget<P::Fp2Params, P::Fp>,
+}
+
+impl<P: MNT4Parameters> ToBytesGadget<P::Fp> for AteDoubleCoefficientsGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_h = self.c_h.to_bytes(&mut cs.ns(|| "c_h to bytes"))?;
+        let mut c_4c = self.c_4c.to_bytes(&mut cs.ns(|| "c_4c to bytes"))?;
+        let mut c_j = self.c_j.to_bytes(&mut cs.ns(|| "c_j to bytes"))?;
+        let mut c_l = self.c_l.to_bytes(&mut cs.ns(|| "c_l to bytes"))?;
+
+        c_h.append(&mut c_4c);
+        c_h.append(&mut c_j);
+        c_h.append(&mut c_l);
+        Ok(c_h)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_h = self
+            .c_h
+            .to_non_unique_bytes(&mut cs.ns(|| "c_h to bytes"))?;
+        let mut c_4c = self
+            .c_4c
+            .to_non_unique_bytes(&mut cs.ns(|| "c_4c to bytes"))?;
+        let mut c_j = self
+            .c_j
+            .to_non_unique_bytes(&mut cs.ns(|| "c_j to bytes"))?;
+        let mut c_l = self
+            .c_l
+            .to_non_unique_bytes(&mut cs.ns(|| "c_l to bytes"))?;
+
+        c_h.append(&mut c_4c);
+        c_h.append(&mut c_j);
+        c_h.append(&mut c_l);
+        Ok(c_h)
+    }
+}
+
+impl<P: MNT4Parameters> AteDoubleCoefficientsGadget<P> {
+    pub fn get_value(&self) -> Option<AteDoubleCoefficients<P>> {
+        match (
+            self.c_h.get_value(),
+            self.c_4c.get_value(),
+            self.c_j.get_value(),
+            self.c_l.get_value(),
+        ) {
+            (Some(c_h), Some(c_4c), Some(c_j), Some(c_l)) => Some(AteDoubleCoefficients {
+                c_h,
+                c_4c,
+                c_j,
+                c_l,
+            }),
+            _ => None,
+        }
+    }
+}
+
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT4Parameters"), Debug(bound = "P: MNT4Parameters"))]
+pub struct AteAdditionCoefficientsGadget<P: MNT4Parameters> {
+    pub c_l1: Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub c_rz: Fp2Gadget<P::Fp2Params, P::Fp>,
+}
+
+impl<P: MNT4Parameters> ToBytesGadget<P::Fp> for AteAdditionCoefficientsGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_l1 = self.c_l1.to_bytes(&mut cs.ns(|| "c_l1 to bytes"))?;
+        let mut c_rz = self.c_rz.to_bytes(&mut cs.ns(|| "c_rz to bytes"))?;
+
+        c_l1.append(&mut c_rz);
+        Ok(c_l1)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_l1 = self
+            .c_l1
+            .to_non_unique_bytes(&mut cs.ns(|| "c_l1 to bytes"))?;
+        let mut c_rz = self
+            .c_rz
+            .to_non_unique_bytes(&mut cs.ns(|| "c_rz to bytes"))?;
+
+        c_l1.append(&mut c_rz);
+        Ok(c_l1)
+    }
+}
+
+impl<P: MNT4Parameters> AteAdditionCoefficientsGadget<P> {
+    pub fn get_value(&self) -> Option<AteAdditionCoefficients<P>> {
+        match (self.c_l1.get_value(), self.c_rz.get_value()) {
+            (Some(c_l1), Some(c_rz)) => Some(AteAdditionCoefficients { c_l1, c_rz }),
+            _ => None,
+        }
+    }
+}
+
+pub struct G2ProjectiveExtendedGadget<P: MNT4Parameters> {
+    pub x: Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub y: Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub z: Fp2Gadget<P::Fp2Params, P::Fp>,
+    pub t: Fp2Gadget<P::Fp2Params, P::Fp>,
+}

r1cs-std/src/groups/curves/short_weierstrass/mnt6/mod.rs

@@ -0,0 +1,423 @@
+use algebra::{
+    curves::mnt6::{
+        g2::{AteAdditionCoefficients, AteDoubleCoefficients},
+        G1Prepared, G2Prepared, MNT6Parameters,
+    },
+    Field,
+};
+use r1cs_core::{ConstraintSystem, SynthesisError};
+
+use crate::{
+    fields::{fp::FpGadget, fp3::Fp3Gadget, FieldGadget},
+    groups::curves::short_weierstrass::AffineGadget,
+    pairing::mnt6::PairingGadget,
+    prelude::*,
+    Vec,
+};
+
+pub type G1Gadget<P> = AffineGadget<
+    <P as MNT6Parameters>::G1Parameters,
+    <P as MNT6Parameters>::Fp,
+    FpGadget<<P as MNT6Parameters>::Fp>,
+>;
+
+pub type G2Gadget<P> =
+    AffineGadget<<P as MNT6Parameters>::G2Parameters, <P as MNT6Parameters>::Fp, Fp3G<P>>;
+
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT6Parameters"), Debug(bound = "P: MNT6Parameters"))]
+pub struct G1PreparedGadget<P: MNT6Parameters> {
+    pub x:       FpGadget<P::Fp>,
+    pub y:       FpGadget<P::Fp>,
+    pub x_twist: Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub y_twist: Fp3Gadget<P::Fp3Params, P::Fp>,
+}
+
+impl<P: MNT6Parameters> G1PreparedGadget<P> {
+    pub fn get_value(&self) -> Option<G1Prepared<P>> {
+        match (
+            self.x.get_value(),
+            self.y.get_value(),
+            self.x_twist.get_value(),
+            self.y_twist.get_value(),
+        ) {
+            (Some(x), Some(y), Some(x_twist), Some(y_twist)) => Some(G1Prepared {
+                x,
+                y,
+                x_twist,
+                y_twist,
+            }),
+            _ => None,
+        }
+    }
+
+    pub fn from_affine<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        q: &G1Gadget<P>,
+    ) -> Result<Self, SynthesisError> {
+        let x_twist = Fp3Gadget::new(
+            q.x.mul_by_constant(cs.ns(|| "g1.x * twist.c0"), &P::TWIST.c0)?,
+            q.x.mul_by_constant(cs.ns(|| "g1.x * twist.c1"), &P::TWIST.c1)?,
+            q.x.mul_by_constant(cs.ns(|| "g1.x * twist.c2"), &P::TWIST.c2)?,
+        );
+        let y_twist = Fp3Gadget::new(
+            q.y.mul_by_constant(cs.ns(|| "g1.y * twist.c0"), &P::TWIST.c0)?,
+            q.y.mul_by_constant(cs.ns(|| "g1.y * twist.c1"), &P::TWIST.c1)?,
+            q.y.mul_by_constant(cs.ns(|| "g1.y * twist.c2"), &P::TWIST.c2)?,
+        );
+        Ok(G1PreparedGadget {
+            x: q.x.clone(),
+            y: q.y.clone(),
+            x_twist,
+            y_twist,
+        })
+    }
+}
+
+impl<P: MNT6Parameters> ToBytesGadget<P::Fp> for G1PreparedGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_twist = self.x_twist.to_bytes(&mut cs.ns(|| "x_twist to bytes"))?;
+        let mut y_twist = self.y_twist.to_bytes(&mut cs.ns(|| "y_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_twist);
+        x.append(&mut y_twist);
+        Ok(x)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_non_unique_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_non_unique_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_twist = self
+            .x_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "x_twist to bytes"))?;
+        let mut y_twist = self
+            .y_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "y_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_twist);
+        x.append(&mut y_twist);
+        Ok(x)
+    }
+}
+
+type Fp3G<P> = Fp3Gadget<<P as MNT6Parameters>::Fp3Params, <P as MNT6Parameters>::Fp>;
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT6Parameters"), Debug(bound = "P: MNT6Parameters"))]
+pub struct G2PreparedGadget<P: MNT6Parameters> {
+    pub x:                     Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub y:                     Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub x_over_twist:          Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub y_over_twist:          Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub double_coefficients:   Vec<AteDoubleCoefficientsGadget<P>>,
+    pub addition_coefficients: Vec<AteAdditionCoefficientsGadget<P>>,
+}
+
+impl<P: MNT6Parameters> ToBytesGadget<P::Fp> for G2PreparedGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_over_twist = self
+            .x_over_twist
+            .to_bytes(&mut cs.ns(|| "x_over_twist to bytes"))?;
+        let mut y_over_twist = self
+            .y_over_twist
+            .to_bytes(&mut cs.ns(|| "y_over_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_over_twist);
+        x.append(&mut y_over_twist);
+
+        for (i, coeff) in self.double_coefficients.iter().enumerate() {
+            x.extend_from_slice(&coeff.to_bytes(cs.ns(|| format!("double_coefficients {}", i)))?);
+        }
+        for (i, coeff) in self.addition_coefficients.iter().enumerate() {
+            x.extend_from_slice(&coeff.to_bytes(cs.ns(|| format!("addition_coefficients {}", i)))?);
+        }
+        Ok(x)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut x = self.x.to_non_unique_bytes(&mut cs.ns(|| "x to bytes"))?;
+        let mut y = self.y.to_non_unique_bytes(&mut cs.ns(|| "y to bytes"))?;
+        let mut x_over_twist = self
+            .x_over_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "x_over_twist to bytes"))?;
+        let mut y_over_twist = self
+            .y_over_twist
+            .to_non_unique_bytes(&mut cs.ns(|| "y_over_twist to bytes"))?;
+
+        x.append(&mut y);
+        x.append(&mut x_over_twist);
+        x.append(&mut y_over_twist);
+
+        for (i, coeff) in self.double_coefficients.iter().enumerate() {
+            x.extend_from_slice(
+                &coeff.to_non_unique_bytes(cs.ns(|| format!("double_coefficients {}", i)))?,
+            );
+        }
+        for (i, coeff) in self.addition_coefficients.iter().enumerate() {
+            x.extend_from_slice(
+                &coeff.to_non_unique_bytes(cs.ns(|| format!("addition_coefficients {}", i)))?,
+            );
+        }
+        Ok(x)
+    }
+}
+
+impl<P: MNT6Parameters> G2PreparedGadget<P> {
+    pub fn get_value(&self) -> Option<G2Prepared<P>> {
+        match (
+            self.x.get_value(),
+            self.y.get_value(),
+            self.x_over_twist.get_value(),
+            self.y_over_twist.get_value(),
+            self.double_coefficients
+                .iter()
+                .map(|coeff| coeff.get_value())
+                .collect::<Option<Vec<AteDoubleCoefficients<P>>>>(),
+            self.addition_coefficients
+                .iter()
+                .map(|coeff| coeff.get_value())
+                .collect::<Option<Vec<AteAdditionCoefficients<P>>>>(),
+        ) {
+            (
+                Some(x),
+                Some(y),
+                Some(x_over_twist),
+                Some(y_over_twist),
+                Some(double_coefficients),
+                Some(addition_coefficients),
+            ) => Some(G2Prepared {
+                x,
+                y,
+                x_over_twist,
+                y_over_twist,
+                double_coefficients,
+                addition_coefficients,
+            }),
+            _ => None,
+        }
+    }
+
+    pub fn from_affine<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        q: &G2Gadget<P>,
+    ) -> Result<Self, SynthesisError> {
+        let twist_inv = P::TWIST.inverse().unwrap();
+
+        let mut g2p = G2PreparedGadget {
+            x:                     q.x.clone(),
+            y:                     q.y.clone(),
+            x_over_twist:          q.x.mul_by_constant(cs.ns(|| "x over twist"), &twist_inv)?,
+            y_over_twist:          q.y.mul_by_constant(cs.ns(|| "y over twist"), &twist_inv)?,
+            double_coefficients:   vec![],
+            addition_coefficients: vec![],
+        };
+
+        let fp2_one = Fp3G::<P>::one(cs.ns(|| "one"))?;
+        let mut r = G2ProjectiveExtendedGadget {
+            x: q.x.clone(),
+            y: q.y.clone(),
+            z: fp2_one.clone(),
+            t: fp2_one,
+        };
+
+        for (idx, value) in P::ATE_LOOP_COUNT.iter().rev().enumerate() {
+            let mut tmp = *value;
+            let skip_extraneous_bits = 64 - value.leading_zeros();
+            let mut v = Vec::with_capacity(16);
+            for i in 0..64 {
+                if idx == 0 && (i == 0 || i >= skip_extraneous_bits) {
+                    continue;
+                }
+                v.push(tmp & 1 == 1);
+                tmp >>= 1;
+            }
+
+            let mut cs = cs.ns(|| format!("ate loop iteration {}", idx));
+
+            for (j, bit) in v.iter().rev().enumerate() {
+                let (r2, coeff) = PairingGadget::<P>::doubling_step_for_flipped_miller_loop(
+                    cs.ns(|| format!("doubling step {}", j)),
+                    &r,
+                )?;
+                g2p.double_coefficients.push(coeff);
+                r = r2;
+
+                if *bit {
+                    let (r2, coeff) =
+                        PairingGadget::<P>::mixed_addition_step_for_flipped_miller_loop(
+                            cs.ns(|| format!("mixed addition step {}", j)),
+                            &q.x,
+                            &q.y,
+                            &r,
+                        )?;
+                    g2p.addition_coefficients.push(coeff);
+                    r = r2;
+                }
+
+                tmp >>= 1;
+            }
+        }
+
+        if P::ATE_IS_LOOP_COUNT_NEG {
+            let rz_inv = r.z.inverse(cs.ns(|| "inverse r.z"))?;
+            let rz2_inv = rz_inv.square(cs.ns(|| "rz_inv^2"))?;
+            let rz3_inv = rz_inv.mul(cs.ns(|| "rz_inv * rz_inv^2"), &rz2_inv)?;
+
+            let minus_r_affine_x = r.x.mul(cs.ns(|| "r.x * rz2_inv"), &rz2_inv)?;
+            let minus_r_affine_y =
+                r.y.negate(cs.ns(|| "-r.y"))?
+                    .mul(cs.ns(|| "-r.y * rz3_inv"), &rz3_inv)?;
+
+            let add_result = PairingGadget::<P>::mixed_addition_step_for_flipped_miller_loop(
+                cs.ns(|| "mixed_addition step"),
+                &minus_r_affine_x,
+                &minus_r_affine_y,
+                &r,
+            )?;
+            g2p.addition_coefficients.push(add_result.1);
+        }
+
+        Ok(g2p)
+    }
+}
+
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT6Parameters"), Debug(bound = "P: MNT6Parameters"))]
+pub struct AteDoubleCoefficientsGadget<P: MNT6Parameters> {
+    pub c_h:  Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub c_4c: Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub c_j:  Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub c_l:  Fp3Gadget<P::Fp3Params, P::Fp>,
+}
+
+impl<P: MNT6Parameters> ToBytesGadget<P::Fp> for AteDoubleCoefficientsGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_h = self.c_h.to_bytes(&mut cs.ns(|| "c_h to bytes"))?;
+        let mut c_4c = self.c_4c.to_bytes(&mut cs.ns(|| "c_4c to bytes"))?;
+        let mut c_j = self.c_j.to_bytes(&mut cs.ns(|| "c_j to bytes"))?;
+        let mut c_l = self.c_l.to_bytes(&mut cs.ns(|| "c_l to bytes"))?;
+
+        c_h.append(&mut c_4c);
+        c_h.append(&mut c_j);
+        c_h.append(&mut c_l);
+        Ok(c_h)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_h = self
+            .c_h
+            .to_non_unique_bytes(&mut cs.ns(|| "c_h to bytes"))?;
+        let mut c_4c = self
+            .c_4c
+            .to_non_unique_bytes(&mut cs.ns(|| "c_4c to bytes"))?;
+        let mut c_j = self
+            .c_j
+            .to_non_unique_bytes(&mut cs.ns(|| "c_j to bytes"))?;
+        let mut c_l = self
+            .c_l
+            .to_non_unique_bytes(&mut cs.ns(|| "c_l to bytes"))?;
+
+        c_h.append(&mut c_4c);
+        c_h.append(&mut c_j);
+        c_h.append(&mut c_l);
+        Ok(c_h)
+    }
+}
+
+impl<P: MNT6Parameters> AteDoubleCoefficientsGadget<P> {
+    pub fn get_value(&self) -> Option<AteDoubleCoefficients<P>> {
+        match (
+            self.c_h.get_value(),
+            self.c_4c.get_value(),
+            self.c_j.get_value(),
+            self.c_l.get_value(),
+        ) {
+            (Some(c_h), Some(c_4c), Some(c_j), Some(c_l)) => Some(AteDoubleCoefficients {
+                c_h,
+                c_4c,
+                c_j,
+                c_l,
+            }),
+            _ => None,
+        }
+    }
+}
+
+#[derive(Derivative)]
+#[derivative(Clone(bound = "P: MNT6Parameters"), Debug(bound = "P: MNT6Parameters"))]
+pub struct AteAdditionCoefficientsGadget<P: MNT6Parameters> {
+    pub c_l1: Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub c_rz: Fp3Gadget<P::Fp3Params, P::Fp>,
+}
+
+impl<P: MNT6Parameters> ToBytesGadget<P::Fp> for AteAdditionCoefficientsGadget<P> {
+    #[inline]
+    fn to_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_l1 = self.c_l1.to_bytes(&mut cs.ns(|| "c_l1 to bytes"))?;
+        let mut c_rz = self.c_rz.to_bytes(&mut cs.ns(|| "c_rz to bytes"))?;
+
+        c_l1.append(&mut c_rz);
+        Ok(c_l1)
+    }
+
+    fn to_non_unique_bytes<CS: ConstraintSystem<P::Fp>>(
+        &self,
+        mut cs: CS,
+    ) -> Result<Vec<UInt8>, SynthesisError> {
+        let mut c_l1 = self
+            .c_l1
+            .to_non_unique_bytes(&mut cs.ns(|| "c_l1 to bytes"))?;
+        let mut c_rz = self
+            .c_rz
+            .to_non_unique_bytes(&mut cs.ns(|| "c_rz to bytes"))?;
+
+        c_l1.append(&mut c_rz);
+        Ok(c_l1)
+    }
+}
+
+impl<P: MNT6Parameters> AteAdditionCoefficientsGadget<P> {
+    pub fn get_value(&self) -> Option<AteAdditionCoefficients<P>> {
+        match (self.c_l1.get_value(), self.c_rz.get_value()) {
+            (Some(c_l1), Some(c_rz)) => Some(AteAdditionCoefficients { c_l1, c_rz }),
+            _ => None,
+        }
+    }
+}
+
+pub struct G2ProjectiveExtendedGadget<P: MNT6Parameters> {
+    pub x: Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub y: Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub z: Fp3Gadget<P::Fp3Params, P::Fp>,
+    pub t: Fp3Gadget<P::Fp3Params, P::Fp>,
+}

r1cs-std/src/groups/curves/short_weierstrass/mod.rs

@@ -11,6 +11,8 @@ use r1cs_core::{ConstraintSystem, SynthesisError};
 use crate::{prelude::*, Assignment, Vec};
 
 pub mod bls12;
+pub mod mnt4;
+pub mod mnt6;
 
 #[derive(Derivative)]
 #[derivative(Debug, Clone)]
@@ -321,7 +323,7 @@ where
     }
 
     fn cost_of_add() -> usize {
-        3 * F::cost_of_mul() + F::cost_of_inv()
+        3 * F::cost_of_mul_equals() + F::cost_of_inv()
     }
 
     fn cost_of_double() -> usize {

r1cs-std/src/groups/mod.rs

@@ -6,7 +6,7 @@ use core::{borrow::Borrow, fmt::Debug};
 
 pub mod curves;
 
-pub use self::curves::short_weierstrass::bls12;
+pub use self::curves::short_weierstrass::{bls12, mnt4, mnt6};
 
 pub trait GroupGadget<G: Group, ConstraintF: Field>:
     Sized

r1cs-std/src/instantiated/mnt4_298/curves.rs

@@ -0,0 +1,198 @@
+use crate::groups::mnt4;
+use algebra::mnt4_298::Parameters;
+
+pub type G1Gadget = mnt4::G1Gadget<Parameters>;
+pub type G2Gadget = mnt4::G2Gadget<Parameters>;
+
+pub type G1PreparedGadget = mnt4::G1PreparedGadget<Parameters>;
+pub type G2PreparedGadget = mnt4::G2PreparedGadget<Parameters>;
+
+#[cfg(test)]
+mod test {
+    use rand::Rng;
+
+    use super::{G1Gadget, G2Gadget};
+    use crate::{prelude::*, test_constraint_system::TestConstraintSystem, Vec};
+    use algebra::{mnt4_298::*, test_rng, AffineCurve, BitIterator, PrimeField, ProjectiveCurve};
+    use r1cs_core::ConstraintSystem;
+
+    #[test]
+    fn mnt4_298_g1_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G1Projective = rng.gen();
+        let b: G1Projective = rng.gen();
+        let gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G1Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G1Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G1Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt4_298_g2_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G2Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G2Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G2Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt4_298_g1_gadget_test() {
+        use algebra::UniformRand;
+        use rand::SeedableRng;
+        use rand_xorshift::XorShiftRng;
+        let mut rng = XorShiftRng::seed_from_u64(1231275789u64);
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let a = G1Projective::rand(&mut rng);
+        let b = G1Projective::rand(&mut rng);
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+        let mut gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.value.unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.value.unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.value.unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.value.unwrap(), b_affine.y);
+
+        // Check addition
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+
+        let ab_val = gadget_ab
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(ab_val, ab_affine, "Result of addition is unequal");
+
+        // Check doubling
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+        let aa_val = gadget_a
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(
+            aa_val, aa_affine,
+            "Gadget and native values are unequal after double."
+        );
+
+        // Check mul_bits
+        let scalar = Fr::rand(&mut rng);
+        let native_result = aa.into_affine().mul(scalar) + &b;
+        let native_result = native_result.into_affine();
+
+        let mut scalar: Vec<bool> = BitIterator::new(scalar.into_repr()).collect();
+        // Get the scalar bits into little-endian form.
+        scalar.reverse();
+        let input = Vec::<Boolean>::alloc(cs.ns(|| "Input"), || Ok(scalar)).unwrap();
+        let result = gadget_a
+            .mul_bits(cs.ns(|| "mul_bits"), &gadget_b, input.iter())
+            .unwrap();
+        let result_val = result.get_value().unwrap().into_affine();
+        assert_eq!(
+            result_val, native_result,
+            "gadget & native values are diff. after scalar mul"
+        );
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+
+    #[test]
+    fn mnt4_298_g2_gadget_test() {
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+
+        let mut gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.get_value().unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.get_value().unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.get_value().unwrap(), b_affine.y);
+
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+        assert_eq!(gadget_ab.x.get_value().unwrap(), ab_affine.x);
+        assert_eq!(gadget_ab.y.get_value().unwrap(), ab_affine.y);
+
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+
+        assert_eq!(gadget_a.x.get_value().unwrap(), aa_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), aa_affine.y);
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+}

r1cs-std/src/instantiated/mnt4_298/fields.rs

@@ -0,0 +1,23 @@
+use algebra::mnt4_298::{Fq, Fq2Parameters, Fq4Parameters};
+
+use crate::fields::{fp::FpGadget, fp2::Fp2Gadget, fp4::Fp4Gadget};
+
+pub type FqGadget = FpGadget<Fq>;
+pub type Fq2Gadget = Fp2Gadget<Fq2Parameters, Fq>;
+pub type Fq4Gadget = Fp4Gadget<Fq4Parameters, Fq>;
+
+#[test]
+fn mnt4_298_field_gadgets_test() {
+    use super::*;
+    use crate::fields::tests::*;
+    use algebra::mnt4_298::{Fq, Fq2, Fq4};
+
+    field_test::<_, Fq, FqGadget>();
+    frobenius_tests::<Fq, Fq, FqGadget>(13);
+
+    field_test::<_, Fq, Fq2Gadget>();
+    frobenius_tests::<Fq2, Fq, Fq2Gadget>(13);
+
+    field_test::<_, Fq, Fq4Gadget>();
+    frobenius_tests::<Fq4, Fq, Fq4Gadget>(13);
+}

r1cs-std/src/instantiated/mnt4_298/mod.rs

@@ -0,0 +1,7 @@
+mod curves;
+mod fields;
+mod pairing;
+
+pub use curves::*;
+pub use fields::*;
+pub use pairing::*;

r1cs-std/src/instantiated/mnt4_298/pairing.rs

@@ -0,0 +1,8 @@
+use algebra::mnt4_298::Parameters;
+
+pub type PairingGadget = crate::pairing::mnt4::PairingGadget<Parameters>;
+
+#[test]
+fn test() {
+    crate::pairing::tests::bilinearity_test::<algebra::MNT4_298, _, PairingGadget>()
+}

r1cs-std/src/instantiated/mnt4_753/curves.rs

@@ -0,0 +1,198 @@
+use crate::groups::mnt4;
+use algebra::mnt4_753::Parameters;
+
+pub type G1Gadget = mnt4::G1Gadget<Parameters>;
+pub type G2Gadget = mnt4::G2Gadget<Parameters>;
+
+pub type G1PreparedGadget = mnt4::G1PreparedGadget<Parameters>;
+pub type G2PreparedGadget = mnt4::G2PreparedGadget<Parameters>;
+
+#[cfg(test)]
+mod test {
+    use rand::Rng;
+
+    use super::{G1Gadget, G2Gadget};
+    use crate::{prelude::*, test_constraint_system::TestConstraintSystem, Vec};
+    use algebra::{mnt4_753::*, test_rng, AffineCurve, BitIterator, PrimeField, ProjectiveCurve};
+    use r1cs_core::ConstraintSystem;
+
+    #[test]
+    fn mnt4_753_g1_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G1Projective = rng.gen();
+        let b: G1Projective = rng.gen();
+        let gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G1Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G1Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G1Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt4_753_g2_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G2Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G2Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G2Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt4_753_g1_gadget_test() {
+        use algebra::UniformRand;
+        use rand::SeedableRng;
+        use rand_xorshift::XorShiftRng;
+        let mut rng = XorShiftRng::seed_from_u64(1231275789u64);
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let a = G1Projective::rand(&mut rng);
+        let b = G1Projective::rand(&mut rng);
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+        let mut gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.value.unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.value.unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.value.unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.value.unwrap(), b_affine.y);
+
+        // Check addition
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+
+        let ab_val = gadget_ab
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(ab_val, ab_affine, "Result of addition is unequal");
+
+        // Check doubling
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+        let aa_val = gadget_a
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(
+            aa_val, aa_affine,
+            "Gadget and native values are unequal after double."
+        );
+
+        // Check mul_bits
+        let scalar = Fr::rand(&mut rng);
+        let native_result = aa.into_affine().mul(scalar) + &b;
+        let native_result = native_result.into_affine();
+
+        let mut scalar: Vec<bool> = BitIterator::new(scalar.into_repr()).collect();
+        // Get the scalar bits into little-endian form.
+        scalar.reverse();
+        let input = Vec::<Boolean>::alloc(cs.ns(|| "Input"), || Ok(scalar)).unwrap();
+        let result = gadget_a
+            .mul_bits(cs.ns(|| "mul_bits"), &gadget_b, input.iter())
+            .unwrap();
+        let result_val = result.get_value().unwrap().into_affine();
+        assert_eq!(
+            result_val, native_result,
+            "gadget & native values are diff. after scalar mul"
+        );
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+
+    #[test]
+    fn mnt4_753_g2_gadget_test() {
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+
+        let mut gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.get_value().unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.get_value().unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.get_value().unwrap(), b_affine.y);
+
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+        assert_eq!(gadget_ab.x.get_value().unwrap(), ab_affine.x);
+        assert_eq!(gadget_ab.y.get_value().unwrap(), ab_affine.y);
+
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+
+        assert_eq!(gadget_a.x.get_value().unwrap(), aa_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), aa_affine.y);
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+}

r1cs-std/src/instantiated/mnt4_753/fields.rs

@@ -0,0 +1,23 @@
+use algebra::mnt4_753::{Fq, Fq2Parameters, Fq4Parameters};
+
+use crate::fields::{fp::FpGadget, fp2::Fp2Gadget, fp4::Fp4Gadget};
+
+pub type FqGadget = FpGadget<Fq>;
+pub type Fq2Gadget = Fp2Gadget<Fq2Parameters, Fq>;
+pub type Fq4Gadget = Fp4Gadget<Fq4Parameters, Fq>;
+
+#[test]
+fn mnt4_753_field_gadgets_test() {
+    use super::*;
+    use crate::fields::tests::*;
+    use algebra::mnt4_753::{Fq, Fq2, Fq4};
+
+    field_test::<_, Fq, FqGadget>();
+    frobenius_tests::<Fq, Fq, FqGadget>(13);
+
+    field_test::<_, Fq, Fq2Gadget>();
+    frobenius_tests::<Fq2, Fq, Fq2Gadget>(13);
+
+    field_test::<_, Fq, Fq4Gadget>();
+    frobenius_tests::<Fq4, Fq, Fq4Gadget>(13);
+}

r1cs-std/src/instantiated/mnt4_753/mod.rs

@@ -0,0 +1,7 @@
+mod curves;
+mod fields;
+mod pairing;
+
+pub use curves::*;
+pub use fields::*;
+pub use pairing::*;

r1cs-std/src/instantiated/mnt4_753/pairing.rs

@@ -0,0 +1,8 @@
+use algebra::mnt4_753::Parameters;
+
+pub type PairingGadget = crate::pairing::mnt4::PairingGadget<Parameters>;
+
+#[test]
+fn test() {
+    crate::pairing::tests::bilinearity_test::<algebra::MNT4_753, _, PairingGadget>()
+}

r1cs-std/src/instantiated/mnt6_298/curves.rs

@@ -0,0 +1,198 @@
+use crate::groups::mnt6;
+use algebra::mnt6_298::Parameters;
+
+pub type G1Gadget = mnt6::G1Gadget<Parameters>;
+pub type G2Gadget = mnt6::G2Gadget<Parameters>;
+
+pub type G1PreparedGadget = mnt6::G1PreparedGadget<Parameters>;
+pub type G2PreparedGadget = mnt6::G2PreparedGadget<Parameters>;
+
+#[cfg(test)]
+mod test {
+    use rand::Rng;
+
+    use super::{G1Gadget, G2Gadget};
+    use crate::{prelude::*, test_constraint_system::TestConstraintSystem, Vec};
+    use algebra::{mnt6_298::*, test_rng, AffineCurve, BitIterator, PrimeField, ProjectiveCurve};
+    use r1cs_core::ConstraintSystem;
+
+    #[test]
+    fn mnt6_298_g1_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G1Projective = rng.gen();
+        let b: G1Projective = rng.gen();
+        let gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G1Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G1Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G1Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt6_298_g2_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G2Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G2Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G2Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt6_298_g1_gadget_test() {
+        use algebra::UniformRand;
+        use rand::SeedableRng;
+        use rand_xorshift::XorShiftRng;
+        let mut rng = XorShiftRng::seed_from_u64(1231275789u64);
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let a = G1Projective::rand(&mut rng);
+        let b = G1Projective::rand(&mut rng);
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+        let mut gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.value.unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.value.unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.value.unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.value.unwrap(), b_affine.y);
+
+        // Check addition
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+
+        let ab_val = gadget_ab
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(ab_val, ab_affine, "Result of addition is unequal");
+
+        // Check doubling
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+        let aa_val = gadget_a
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(
+            aa_val, aa_affine,
+            "Gadget and native values are unequal after double."
+        );
+
+        // Check mul_bits
+        let scalar = Fr::rand(&mut rng);
+        let native_result = aa.into_affine().mul(scalar) + &b;
+        let native_result = native_result.into_affine();
+
+        let mut scalar: Vec<bool> = BitIterator::new(scalar.into_repr()).collect();
+        // Get the scalar bits into little-endian form.
+        scalar.reverse();
+        let input = Vec::<Boolean>::alloc(cs.ns(|| "Input"), || Ok(scalar)).unwrap();
+        let result = gadget_a
+            .mul_bits(cs.ns(|| "mul_bits"), &gadget_b, input.iter())
+            .unwrap();
+        let result_val = result.get_value().unwrap().into_affine();
+        assert_eq!(
+            result_val, native_result,
+            "gadget & native values are diff. after scalar mul"
+        );
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+
+    #[test]
+    fn mnt6_298_g2_gadget_test() {
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+
+        let mut gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.get_value().unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.get_value().unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.get_value().unwrap(), b_affine.y);
+
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+        assert_eq!(gadget_ab.x.get_value().unwrap(), ab_affine.x);
+        assert_eq!(gadget_ab.y.get_value().unwrap(), ab_affine.y);
+
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+
+        assert_eq!(gadget_a.x.get_value().unwrap(), aa_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), aa_affine.y);
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+}

r1cs-std/src/instantiated/mnt6_298/fields.rs

@@ -0,0 +1,23 @@
+use algebra::mnt6_298::{Fq, Fq3Parameters, Fq6Parameters};
+
+use crate::fields::{fp::FpGadget, fp3::Fp3Gadget, fp6_2over3::Fp6Gadget};
+
+pub type FqGadget = FpGadget<Fq>;
+pub type Fq3Gadget = Fp3Gadget<Fq3Parameters, Fq>;
+pub type Fq6Gadget = Fp6Gadget<Fq6Parameters, Fq>;
+
+#[test]
+fn mnt6_298_field_gadgets_test() {
+    use super::*;
+    use crate::fields::tests::*;
+    use algebra::mnt6_298::{Fq, Fq3, Fq6};
+
+    field_test::<_, Fq, FqGadget>();
+    frobenius_tests::<Fq, Fq, FqGadget>(13);
+
+    field_test::<_, Fq, Fq3Gadget>();
+    frobenius_tests::<Fq3, Fq, Fq3Gadget>(13);
+
+    field_test::<_, Fq, Fq6Gadget>();
+    frobenius_tests::<Fq6, Fq, Fq6Gadget>(13);
+}

r1cs-std/src/instantiated/mnt6_298/mod.rs

@@ -0,0 +1,7 @@
+mod curves;
+mod fields;
+mod pairing;
+
+pub use curves::*;
+pub use fields::*;
+pub use pairing::*;

r1cs-std/src/instantiated/mnt6_298/pairing.rs

@@ -0,0 +1,8 @@
+use algebra::mnt6_298::Parameters;
+
+pub type PairingGadget = crate::pairing::mnt6::PairingGadget<Parameters>;
+
+#[test]
+fn test() {
+    crate::pairing::tests::bilinearity_test::<algebra::MNT6_298, _, PairingGadget>()
+}

r1cs-std/src/instantiated/mnt6_753/curves.rs

@@ -0,0 +1,198 @@
+use crate::groups::mnt6;
+use algebra::mnt6_753::Parameters;
+
+pub type G1Gadget = mnt6::G1Gadget<Parameters>;
+pub type G2Gadget = mnt6::G2Gadget<Parameters>;
+
+pub type G1PreparedGadget = mnt6::G1PreparedGadget<Parameters>;
+pub type G2PreparedGadget = mnt6::G2PreparedGadget<Parameters>;
+
+#[cfg(test)]
+mod test {
+    use rand::Rng;
+
+    use super::{G1Gadget, G2Gadget};
+    use crate::{prelude::*, test_constraint_system::TestConstraintSystem, Vec};
+    use algebra::{mnt6_753::*, test_rng, AffineCurve, BitIterator, PrimeField, ProjectiveCurve};
+    use r1cs_core::ConstraintSystem;
+
+    #[test]
+    fn mnt6_753_g1_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G1Projective = rng.gen();
+        let b: G1Projective = rng.gen();
+        let gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G1Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G1Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G1Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt6_753_g2_constraint_costs() {
+        use crate::boolean::AllocatedBit;
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let bit = AllocatedBit::alloc(&mut cs.ns(|| "bool"), || Ok(true))
+            .unwrap()
+            .into();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        let alloc_cost = cs.num_constraints();
+        let _ = G2Gadget::conditionally_select(
+            &mut cs.ns(|| "cond_select"),
+            &bit,
+            &gadget_a,
+            &gadget_b,
+        )
+        .unwrap();
+        let cond_select_cost = cs.num_constraints() - alloc_cost;
+
+        let _ = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let add_cost = cs.num_constraints() - cond_select_cost - alloc_cost;
+
+        assert!(cs.is_satisfied());
+        assert_eq!(cond_select_cost, <G2Gadget as CondSelectGadget<Fq>>::cost());
+        assert_eq!(add_cost, G2Gadget::cost_of_add());
+    }
+
+    #[test]
+    fn mnt6_753_g1_gadget_test() {
+        use algebra::UniformRand;
+        use rand::SeedableRng;
+        use rand_xorshift::XorShiftRng;
+        let mut rng = XorShiftRng::seed_from_u64(1231275789u64);
+
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let a = G1Projective::rand(&mut rng);
+        let b = G1Projective::rand(&mut rng);
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+        let mut gadget_a = G1Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G1Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.value.unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.value.unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.value.unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.value.unwrap(), b_affine.y);
+
+        // Check addition
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+
+        let ab_val = gadget_ab
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(ab_val, ab_affine, "Result of addition is unequal");
+
+        // Check doubling
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+        let aa_val = gadget_a
+            .get_value()
+            .expect("Doubling should be successful")
+            .into_affine();
+        assert_eq!(
+            aa_val, aa_affine,
+            "Gadget and native values are unequal after double."
+        );
+
+        // Check mul_bits
+        let scalar = Fr::rand(&mut rng);
+        let native_result = aa.into_affine().mul(scalar) + &b;
+        let native_result = native_result.into_affine();
+
+        let mut scalar: Vec<bool> = BitIterator::new(scalar.into_repr()).collect();
+        // Get the scalar bits into little-endian form.
+        scalar.reverse();
+        let input = Vec::<Boolean>::alloc(cs.ns(|| "Input"), || Ok(scalar)).unwrap();
+        let result = gadget_a
+            .mul_bits(cs.ns(|| "mul_bits"), &gadget_b, input.iter())
+            .unwrap();
+        let result_val = result.get_value().unwrap().into_affine();
+        assert_eq!(
+            result_val, native_result,
+            "gadget & native values are diff. after scalar mul"
+        );
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+
+    #[test]
+    fn mnt6_753_g2_gadget_test() {
+        let mut cs = TestConstraintSystem::<Fq>::new();
+
+        let mut rng = test_rng();
+        let a: G2Projective = rng.gen();
+        let b: G2Projective = rng.gen();
+        let a_affine = a.into_affine();
+        let b_affine = b.into_affine();
+
+        let mut gadget_a = G2Gadget::alloc(&mut cs.ns(|| "a"), || Ok(a)).unwrap();
+        let gadget_b = G2Gadget::alloc(&mut cs.ns(|| "b"), || Ok(b)).unwrap();
+        assert_eq!(gadget_a.x.get_value().unwrap(), a_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), a_affine.y);
+        assert_eq!(gadget_b.x.get_value().unwrap(), b_affine.x);
+        assert_eq!(gadget_b.y.get_value().unwrap(), b_affine.y);
+
+        let ab = a + &b;
+        let ab_affine = ab.into_affine();
+        let gadget_ab = gadget_a.add(&mut cs.ns(|| "ab"), &gadget_b).unwrap();
+        let gadget_ba = gadget_b.add(&mut cs.ns(|| "ba"), &gadget_a).unwrap();
+        gadget_ba
+            .enforce_equal(&mut cs.ns(|| "b + a == a + b?"), &gadget_ab)
+            .unwrap();
+        assert_eq!(gadget_ab.x.get_value().unwrap(), ab_affine.x);
+        assert_eq!(gadget_ab.y.get_value().unwrap(), ab_affine.y);
+
+        let aa = a.double();
+        let aa_affine = aa.into_affine();
+        gadget_a.double_in_place(&mut cs.ns(|| "2a")).unwrap();
+
+        assert_eq!(gadget_a.x.get_value().unwrap(), aa_affine.x);
+        assert_eq!(gadget_a.y.get_value().unwrap(), aa_affine.y);
+
+        if !cs.is_satisfied() {
+            println!("{:?}", cs.which_is_unsatisfied().unwrap());
+        }
+
+        assert!(cs.is_satisfied());
+    }
+}

r1cs-std/src/instantiated/mnt6_753/fields.rs

@@ -0,0 +1,23 @@
+use algebra::mnt6_753::{Fq, Fq3Parameters, Fq6Parameters};
+
+use crate::fields::{fp::FpGadget, fp3::Fp3Gadget, fp6_2over3::Fp6Gadget};
+
+pub type FqGadget = FpGadget<Fq>;
+pub type Fq3Gadget = Fp3Gadget<Fq3Parameters, Fq>;
+pub type Fq6Gadget = Fp6Gadget<Fq6Parameters, Fq>;
+
+#[test]
+fn mnt6_753_field_gadgets_test() {
+    use super::*;
+    use crate::fields::tests::*;
+    use algebra::mnt6_753::{Fq, Fq3, Fq6};
+
+    field_test::<_, Fq, FqGadget>();
+    frobenius_tests::<Fq, Fq, FqGadget>(13);
+
+    field_test::<_, Fq, Fq3Gadget>();
+    frobenius_tests::<Fq3, Fq, Fq3Gadget>(13);
+
+    field_test::<_, Fq, Fq6Gadget>();
+    frobenius_tests::<Fq6, Fq, Fq6Gadget>(13);
+}

r1cs-std/src/instantiated/mnt6_753/mod.rs

@@ -0,0 +1,7 @@
+mod curves;
+mod fields;
+mod pairing;
+
+pub use curves::*;
+pub use fields::*;
+pub use pairing::*;

r1cs-std/src/instantiated/mnt6_753/pairing.rs

@@ -0,0 +1,8 @@
+use algebra::mnt6_753::Parameters;
+
+pub type PairingGadget = crate::pairing::mnt6::PairingGadget<Parameters>;
+
+#[test]
+fn test() {
+    crate::pairing::tests::bilinearity_test::<algebra::MNT6_753, _, PairingGadget>()
+}

r1cs-std/src/instantiated/mod.rs

@@ -9,3 +9,15 @@ pub mod edwards_sw6;
 
 #[cfg(feature = "jubjub")]
 pub mod jubjub;
+
+#[cfg(feature = "mnt4_298")]
+pub mod mnt4_298;
+
+#[cfg(feature = "mnt4_753")]
+pub mod mnt4_753;
+
+#[cfg(feature = "mnt6_298")]
+pub mod mnt6_298;
+
+#[cfg(feature = "mnt6_753")]
+pub mod mnt6_753;

r1cs-std/src/lib.rs

@@ -57,6 +57,18 @@ pub use instantiated::edwards_sw6;
 #[cfg(feature = "jubjub")]
 pub use instantiated::jubjub;
 
+#[cfg(feature = "mnt4_298")]
+pub use instantiated::mnt4_298;
+
+#[cfg(feature = "mnt4_753")]
+pub use instantiated::mnt4_753;
+
+#[cfg(feature = "mnt6_298")]
+pub use instantiated::mnt6_298;
+
+#[cfg(feature = "mnt6_753")]
+pub use instantiated::mnt6_753;
+
 pub mod pairing;
 
 pub mod alloc;

r1cs-std/src/pairing/mnt4/mod.rs

@@ -0,0 +1,338 @@
+use r1cs_core::{ConstraintSystem, SynthesisError};
+
+use super::PairingGadget as PG;
+
+use crate::{
+    fields::{fp::FpGadget, fp2::Fp2Gadget, fp4::Fp4Gadget, FieldGadget},
+    groups::mnt4::{
+        AteAdditionCoefficientsGadget, AteDoubleCoefficientsGadget, G1Gadget, G1PreparedGadget,
+        G2Gadget, G2PreparedGadget, G2ProjectiveExtendedGadget,
+    },
+};
+use algebra::{
+    curves::mnt4::{MNT4Parameters, MNT4},
+    fields::BitIterator,
+};
+use core::marker::PhantomData;
+
+pub struct PairingGadget<P: MNT4Parameters>(PhantomData<P>);
+
+type Fp2G<P> = Fp2Gadget<<P as MNT4Parameters>::Fp2Params, <P as MNT4Parameters>::Fp>;
+type Fp4G<P> = Fp4Gadget<<P as MNT4Parameters>::Fp4Params, <P as MNT4Parameters>::Fp>;
+pub type GTGadget<P> = Fp4G<P>;
+
+impl<P: MNT4Parameters> PairingGadget<P> {
+    pub(crate) fn doubling_step_for_flipped_miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        r: &G2ProjectiveExtendedGadget<P>,
+    ) -> Result<
+        (
+            G2ProjectiveExtendedGadget<P>,
+            AteDoubleCoefficientsGadget<P>,
+        ),
+        SynthesisError,
+    > {
+        let a = r.t.square(cs.ns(|| "r.t^2"))?;
+        let b = r.x.square(cs.ns(|| "r.x^2"))?;
+        let c = r.y.square(cs.ns(|| "r.y^2"))?;
+        let d = c.square(cs.ns(|| "c^2"))?;
+        let mut e = r.x.add(cs.ns(|| "r.x + c"), &c)?;
+        e.square_in_place(cs.ns(|| "(r.x + c)^2"))?;
+        e.sub_in_place(cs.ns(|| "(r.x + c)^2 - b"), &b)?;
+        e.sub_in_place(cs.ns(|| "(r.x + c)^2 - b - d"), &d)?;
+
+        let mut f = b.double(cs.ns(|| "b + b"))?;
+        f.add_in_place(cs.ns(|| "b + b + b"), &b)?;
+        let twist_a = a.mul_by_constant(cs.ns(|| "TWIST_COEFF_A * a"), &P::TWIST_COEFF_A)?;
+        f.add_in_place(cs.ns(|| "(b + b + b) + (TWIST_COEFF_A * a)"), &twist_a)?;
+        let g = f.square(cs.ns(|| "f^2"))?;
+
+        let d_eight = d
+            .double(cs.ns(|| "2 * d"))?
+            .double(cs.ns(|| "4 * d"))?
+            .double(cs.ns(|| "8 * d"))?;
+
+        let e2 = e.double(cs.ns(|| "2 * e"))?;
+        let e4 = e2.double(cs.ns(|| "4 * e"))?;
+        let x = g.sub(cs.ns(|| "- (e + e + e + e) + g"), &e4)?;
+
+        let mut y = e2.sub(cs.ns(|| "e + e - x"), &x)?;
+        y.mul_in_place(cs.ns(|| "f * (e + e - x)"), &f)?;
+        y.sub_in_place(cs.ns(|| "- d_eight + f * (e + e - x)"), &d_eight)?;
+        let mut z = r.y.add(cs.ns(|| "r.y + r.z"), &r.z)?;
+        z.square_in_place(cs.ns(|| "(r.y + r.z)^2"))?;
+        z.sub_in_place(cs.ns(|| "(r.y + r.z)^2 - c"), &c)?;
+        let z2 = r.z.square(cs.ns(|| "r.z^2"))?;
+        z.sub_in_place(cs.ns(|| "(r.y + r.z)^2 - c - r.z^2"), &z2)?;
+        let t = z.square(cs.ns(|| "z^2"))?;
+
+        let r2 = G2ProjectiveExtendedGadget { x, y, z, t };
+
+        let c_h =
+            r2.z.add(cs.ns(|| "r2.z + r.t"), &r.t)?
+                .square(cs.ns(|| "(r2.z + r.t)^2"))?
+                .sub(cs.ns(|| "(r2.z + r.t)^2 - r2.t"), &r2.t)?
+                .sub(cs.ns(|| "(r2.z + r.t)^2 - r2.t - a"), &a)?;
+        let c_4c = c.double(cs.ns(|| "2 * c"))?.double(cs.ns(|| "4 * c"))?;
+        let mut c_j = f.add(cs.ns(|| "f + r.t"), &r.t)?;
+        c_j.square_in_place(cs.ns(|| "(f + r.t)^2"))?;
+        c_j.sub_in_place(cs.ns(|| "(f + r.t)^2 - g"), &g)?;
+        c_j.sub_in_place(cs.ns(|| "(f + r.t)^2 - g - a"), &a)?;
+        let mut c_l = f.add(cs.ns(|| "f + r.x"), &r.x)?;
+        c_l.square_in_place(cs.ns(|| "(f + r.x)^2"))?;
+        c_l.sub_in_place(cs.ns(|| "(f + r.x)^2 - g"), &g)?;
+        c_l.sub_in_place(cs.ns(|| "(f + r.x)^2 - g - b"), &b)?;
+        let coeff = AteDoubleCoefficientsGadget {
+            c_h,
+            c_4c,
+            c_j,
+            c_l,
+        };
+
+        Ok((r2, coeff))
+    }
+
+    pub(crate) fn mixed_addition_step_for_flipped_miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        x: &Fp2G<P>,
+        y: &Fp2G<P>,
+        r: &G2ProjectiveExtendedGadget<P>,
+    ) -> Result<
+        (
+            G2ProjectiveExtendedGadget<P>,
+            AteAdditionCoefficientsGadget<P>,
+        ),
+        SynthesisError,
+    > {
+        let a = y.square(cs.ns(|| "y^2"))?;
+        let b = r.t.mul(cs.ns(|| "r.t * x"), &x)?;
+        let mut d = r.z.add(cs.ns(|| "r.z + y"), &y)?;
+        d.square_in_place(cs.ns(|| "(r.z + y)^2"))?;
+        d.sub_in_place(cs.ns(|| "(r.z + y)^2 - a"), &a)?;
+        d.sub_in_place(cs.ns(|| "(r.z + y)^2 - a - r.t"), &r.t)?;
+        d.mul_in_place(cs.ns(|| "((r.z + y)^2 - a - r.t) * r.t"), &r.t)?;
+        let h = b.sub(cs.ns(|| "b - r.x"), &r.x)?;
+        let i = h.square(cs.ns(|| "h^2"))?;
+        let e = i.double(cs.ns(|| "2 * i"))?.double(cs.ns(|| "4 * i"))?;
+        let j = h.mul(cs.ns(|| "h * e"), &e)?;
+        let v = r.x.mul(cs.ns(|| "r.x * e"), &e)?;
+        let ry2 = r.y.double(cs.ns(|| "r.y + r.y"))?;
+        let l1 = d.sub(cs.ns(|| "d - (r.y + r.y)"), &ry2)?;
+
+        let v2 = v.double(cs.ns(|| "v + v"))?;
+        let x = l1
+            .square(cs.ns(|| "l1^2"))?
+            .sub(cs.ns(|| "l1^2 - j"), &j)?
+            .sub(cs.ns(|| "l1^2 - j - (v + v)"), &v2)?;
+        let v_minus_x = v.sub(cs.ns(|| "v - x"), &x)?;
+        let j_ry2 = j.mul(cs.ns(|| "j * (r.y + r.y)"), &ry2)?;
+        let y = l1
+            .mul(cs.ns(|| "l1 * (v - x)"), &v_minus_x)?
+            .sub(cs.ns(|| "l1 * (v - x) - (j * (r.y + r.y)"), &j_ry2)?;
+        let mut z = r.z.add(cs.ns(|| "r.z + h"), &h)?;
+        z.square_in_place(cs.ns(|| "(r.z + h)^2"))?;
+        z.sub_in_place(cs.ns(|| "(r.z + h)^2 - r.t"), &r.t)?;
+        z.sub_in_place(cs.ns(|| "(r.z + h)^2 - r.t - i"), &i)?;
+        let t = z.square(cs.ns(|| "z^2"))?;
+
+        let r2 = G2ProjectiveExtendedGadget {
+            x,
+            y,
+            z: z.clone(),
+            t,
+        };
+        let coeff = AteAdditionCoefficientsGadget { c_l1: l1, c_rz: z };
+
+        Ok((r2, coeff))
+    }
+
+    pub fn ate_miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        p: &G1PreparedGadget<P>,
+        q: &G2PreparedGadget<P>,
+    ) -> Result<Fp4G<P>, SynthesisError> {
+        let mut l1_coeff = Fp2G::<P>::new(p.x.clone(), FpGadget::<P::Fp>::zero(cs.ns(|| "zero"))?);
+        l1_coeff.sub_in_place(cs.ns(|| "l1_coeff"), &q.x_over_twist)?;
+
+        let mut f = Fp4G::<P>::one(cs.ns(|| "one"))?;
+
+        let mut dbl_idx: usize = 0;
+        let mut add_idx: usize = 0;
+
+        let mut found_one = false;
+
+        for (j, bit) in BitIterator::new(P::ATE_LOOP_COUNT).enumerate() {
+            // code below gets executed for all bits (EXCEPT the MSB itself) of
+            // mnt6_param_p (skipping leading zeros) in MSB to LSB order
+            if !found_one && bit {
+                found_one = true;
+                continue;
+            } else if !found_one {
+                continue;
+            }
+
+            let mut cs = cs.ns(|| format!("bit {}", j));
+
+            let dc = &q.double_coefficients[dbl_idx];
+            dbl_idx += 1;
+
+            let c_j_x_twist = dc.c_j.mul(cs.ns(|| "dc.c_j * p.x_twist"), &p.x_twist)?;
+            let c0 = dc.c_l.sub(cs.ns(|| "-dc.c_4c + dc.c_l"), &dc.c_4c)?.sub(
+                cs.ns(|| "-dc.c_4c - (dc.c_j * p.x_twist) + dc.c_l"),
+                &c_j_x_twist,
+            )?;
+            let c1 = dc.c_h.mul(cs.ns(|| "dc.c_h * p.y_twist"), &p.y_twist)?;
+            let g_rr_at_p = Fp4G::<P>::new(c0, c1);
+
+            f = f
+                .square(cs.ns(|| "f^2"))?
+                .mul(cs.ns(|| "f^2 * g_rr_at_p"), &g_rr_at_p)?;
+
+            if bit {
+                let ac = &q.addition_coefficients[add_idx];
+                add_idx += 1;
+
+                let l1_coeff_c_l1 = l1_coeff.mul(cs.ns(|| "l1_coeff * ac.c_l1"), &ac.c_l1)?;
+                let g_rq_at_p = Fp4G::<P>::new(
+                    ac.c_rz.mul(cs.ns(|| "ac.c_rz * p.y_twist"), &p.y_twist)?,
+                    q.y_over_twist
+                        .mul(cs.ns(|| "q.y_over_twist * ac.c_rz"), &ac.c_rz)?
+                        .add(
+                            cs.ns(|| "q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1)"),
+                            &l1_coeff_c_l1,
+                        )?
+                        .negate(cs.ns(|| "-(q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1))"))?,
+                );
+                f.mul_in_place(cs.ns(|| "f *= g_rq_at_p"), &g_rq_at_p)?;
+            }
+        }
+
+        if P::ATE_IS_LOOP_COUNT_NEG {
+            let ac = &q.addition_coefficients[add_idx];
+
+            let l1_coeff_c_l1 = l1_coeff.mul(cs.ns(|| "l1_coeff * ac.c_l1"), &ac.c_l1)?;
+            let g_rnegr_at_p = Fp4G::<P>::new(
+                ac.c_rz.mul(cs.ns(|| "ac.c_rz * p.y_twist"), &p.y_twist)?,
+                q.y_over_twist
+                    .mul(cs.ns(|| "q.y_over_twist * ac.c_rz"), &ac.c_rz)?
+                    .add(
+                        cs.ns(|| "q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1)"),
+                        &l1_coeff_c_l1,
+                    )?
+                    .negate(cs.ns(|| "-(q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1))"))?,
+            );
+            f = f
+                .mul(cs.ns(|| "f * g_rnegr_at_p"), &g_rnegr_at_p)?
+                .inverse(cs.ns(|| "inverse f"))?;
+        }
+
+        Ok(f)
+    }
+
+    pub fn final_exponentiation<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        value: &Fp4G<P>,
+    ) -> Result<GTGadget<P>, SynthesisError> {
+        let value_inv = value.inverse(cs.ns(|| "value inverse"))?;
+        let value_to_first_chunk = Self::final_exponentiation_first_chunk(
+            cs.ns(|| "value_to_first_chunk"),
+            value,
+            &value_inv,
+        )?;
+        let value_inv_to_first_chunk = Self::final_exponentiation_first_chunk(
+            cs.ns(|| "value_inv_to_first_chunk"),
+            &value_inv,
+            value,
+        )?;
+        Self::final_exponentiation_last_chunk(
+            cs.ns(|| "final_exp_last_chunk"),
+            &value_to_first_chunk,
+            &value_inv_to_first_chunk,
+        )
+    }
+
+    fn final_exponentiation_first_chunk<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        elt: &Fp4G<P>,
+        elt_inv: &Fp4G<P>,
+    ) -> Result<Fp4G<P>, SynthesisError> {
+        // (q^2-1)
+
+        // elt_q2 = elt^(q^2)
+        let mut elt_q2 = elt.clone();
+        elt_q2.frobenius_map_in_place(cs.ns(|| "frobenius 2"), 2)?;
+        // elt_q2_over_elt = elt^(q^2-1)
+        elt_q2.mul(cs.ns(|| "elt_q2 * elt_inv"), elt_inv)
+    }
+
+    fn final_exponentiation_last_chunk<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        elt: &Fp4G<P>,
+        elt_inv: &Fp4G<P>,
+    ) -> Result<Fp4G<P>, SynthesisError> {
+        let elt_clone = elt.clone();
+        let elt_inv_clone = elt_inv.clone();
+
+        let mut elt_q = elt.clone();
+        elt_q.frobenius_map_in_place(cs.ns(|| "frobenius 1"), 1)?;
+
+        let w1_part = elt_q.cyclotomic_exp(cs.ns(|| "w1_part"), &P::FINAL_EXPONENT_LAST_CHUNK_1)?;
+        let w0_part;
+        if P::FINAL_EXPONENT_LAST_CHUNK_W0_IS_NEG {
+            w0_part = elt_inv_clone
+                .cyclotomic_exp(cs.ns(|| "w0_part"), &P::FINAL_EXPONENT_LAST_CHUNK_ABS_OF_W0)?;
+        } else {
+            w0_part = elt_clone
+                .cyclotomic_exp(cs.ns(|| "w0_part"), &P::FINAL_EXPONENT_LAST_CHUNK_ABS_OF_W0)?;
+        }
+
+        w1_part.mul(cs.ns(|| "w1_part * w0_part"), &w0_part)
+    }
+}
+
+impl<P: MNT4Parameters> PG<MNT4<P>, P::Fp> for PairingGadget<P> {
+    type G1Gadget = G1Gadget<P>;
+    type G2Gadget = G2Gadget<P>;
+    type G1PreparedGadget = G1PreparedGadget<P>;
+    type G2PreparedGadget = G2PreparedGadget<P>;
+    type GTGadget = GTGadget<P>;
+
+    fn miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        ps: &[Self::G1PreparedGadget],
+        qs: &[Self::G2PreparedGadget],
+    ) -> Result<Self::GTGadget, SynthesisError> {
+        let mut result = Fp4G::<P>::one(cs.ns(|| "one"))?;
+        for (i, (p, q)) in ps.iter().zip(qs.iter()).enumerate() {
+            let miller =
+                Self::ate_miller_loop(cs.ns(|| format!("ate miller loop iteration {}", i)), p, q)?;
+            result.mul_in_place(
+                cs.ns(|| format!("mul ate miller loop iteration {}", i)),
+                &miller,
+            )?;
+        }
+
+        Ok(result)
+    }
+
+    fn final_exponentiation<CS: ConstraintSystem<P::Fp>>(
+        cs: CS,
+        r: &Self::GTGadget,
+    ) -> Result<Self::GTGadget, SynthesisError> {
+        Self::final_exponentiation(cs, r)
+    }
+
+    fn prepare_g1<CS: ConstraintSystem<P::Fp>>(
+        cs: CS,
+        p: &Self::G1Gadget,
+    ) -> Result<Self::G1PreparedGadget, SynthesisError> {
+        Self::G1PreparedGadget::from_affine(cs, p)
+    }
+
+    fn prepare_g2<CS: ConstraintSystem<P::Fp>>(
+        cs: CS,
+        q: &Self::G2Gadget,
+    ) -> Result<Self::G2PreparedGadget, SynthesisError> {
+        Self::G2PreparedGadget::from_affine(cs, q)
+    }
+}

r1cs-std/src/pairing/mnt6/mod.rs

@@ -0,0 +1,344 @@
+use r1cs_core::{ConstraintSystem, SynthesisError};
+
+use super::PairingGadget as PG;
+
+use crate::{
+    fields::{fp::FpGadget, fp3::Fp3Gadget, fp6_2over3::Fp6Gadget, FieldGadget},
+    groups::mnt6::{
+        AteAdditionCoefficientsGadget, AteDoubleCoefficientsGadget, G1Gadget, G1PreparedGadget,
+        G2Gadget, G2PreparedGadget, G2ProjectiveExtendedGadget,
+    },
+};
+use algebra::{
+    curves::mnt6::{MNT6Parameters, MNT6},
+    fields::BitIterator,
+};
+use core::marker::PhantomData;
+
+pub struct PairingGadget<P: MNT6Parameters>(PhantomData<P>);
+
+type Fp3G<P> = Fp3Gadget<<P as MNT6Parameters>::Fp3Params, <P as MNT6Parameters>::Fp>;
+type Fp6G<P> = Fp6Gadget<<P as MNT6Parameters>::Fp6Params, <P as MNT6Parameters>::Fp>;
+pub type GTGadget<P> = Fp6G<P>;
+
+impl<P: MNT6Parameters> PairingGadget<P> {
+    pub(crate) fn doubling_step_for_flipped_miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        r: &G2ProjectiveExtendedGadget<P>,
+    ) -> Result<
+        (
+            G2ProjectiveExtendedGadget<P>,
+            AteDoubleCoefficientsGadget<P>,
+        ),
+        SynthesisError,
+    > {
+        let a = r.t.square(cs.ns(|| "r.t^2"))?;
+        let b = r.x.square(cs.ns(|| "r.x^2"))?;
+        let c = r.y.square(cs.ns(|| "r.y^2"))?;
+        let d = c.square(cs.ns(|| "c^2"))?;
+        let mut e = r.x.add(cs.ns(|| "r.x + c"), &c)?;
+        e.square_in_place(cs.ns(|| "(r.x + c)^2"))?;
+        e.sub_in_place(cs.ns(|| "(r.x + c)^2 - b"), &b)?;
+        e.sub_in_place(cs.ns(|| "(r.x + c)^2 - b - d"), &d)?;
+
+        let mut f = b.double(cs.ns(|| "b + b"))?;
+        f.add_in_place(cs.ns(|| "b + b + b"), &b)?;
+        let twist_a = a.mul_by_constant(cs.ns(|| "TWIST_COEFF_A * a"), &P::TWIST_COEFF_A)?;
+        f.add_in_place(cs.ns(|| "(b + b + b) + (TWIST_COEFF_A * a)"), &twist_a)?;
+        let g = f.square(cs.ns(|| "f^2"))?;
+
+        let d_eight = d
+            .double(cs.ns(|| "2 * d"))?
+            .double(cs.ns(|| "4 * d"))?
+            .double(cs.ns(|| "8 * d"))?;
+
+        let e2 = e.double(cs.ns(|| "2 * e"))?;
+        let e4 = e2.double(cs.ns(|| "4 * e"))?;
+        let x = g.sub(cs.ns(|| "- (e + e + e + e) + g"), &e4)?;
+
+        let mut y = e2.sub(cs.ns(|| "e + e - x"), &x)?;
+        y.mul_in_place(cs.ns(|| "f * (e + e - x)"), &f)?;
+        y.sub_in_place(cs.ns(|| "- d_eight + f * (e + e - x)"), &d_eight)?;
+        let mut z = r.y.add(cs.ns(|| "r.y + r.z"), &r.z)?;
+        z.square_in_place(cs.ns(|| "(r.y + r.z)^2"))?;
+        z.sub_in_place(cs.ns(|| "(r.y + r.z)^2 - c"), &c)?;
+        let z2 = r.z.square(cs.ns(|| "r.z^2"))?;
+        z.sub_in_place(cs.ns(|| "(r.y + r.z)^2 - c - r.z^2"), &z2)?;
+        let t = z.square(cs.ns(|| "z^2"))?;
+
+        let r2 = G2ProjectiveExtendedGadget { x, y, z, t };
+
+        let c_h =
+            r2.z.add(cs.ns(|| "r2.z + r.t"), &r.t)?
+                .square(cs.ns(|| "(r2.z + r.t)^2"))?
+                .sub(cs.ns(|| "(r2.z + r.t)^2 - r2.t"), &r2.t)?
+                .sub(cs.ns(|| "(r2.z + r.t)^2 - r2.t - a"), &a)?;
+        let c_4c = c.double(cs.ns(|| "2 * c"))?.double(cs.ns(|| "4 * c"))?;
+        let mut c_j = f.add(cs.ns(|| "f + r.t"), &r.t)?;
+        c_j.square_in_place(cs.ns(|| "(f + r.t)^2"))?;
+        c_j.sub_in_place(cs.ns(|| "(f + r.t)^2 - g"), &g)?;
+        c_j.sub_in_place(cs.ns(|| "(f + r.t)^2 - g - a"), &a)?;
+        let mut c_l = f.add(cs.ns(|| "f + r.x"), &r.x)?;
+        c_l.square_in_place(cs.ns(|| "(f + r.x)^2"))?;
+        c_l.sub_in_place(cs.ns(|| "(f + r.x)^2 - g"), &g)?;
+        c_l.sub_in_place(cs.ns(|| "(f + r.x)^2 - g - b"), &b)?;
+        let coeff = AteDoubleCoefficientsGadget {
+            c_h,
+            c_4c,
+            c_j,
+            c_l,
+        };
+
+        Ok((r2, coeff))
+    }
+
+    pub(crate) fn mixed_addition_step_for_flipped_miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        x: &Fp3G<P>,
+        y: &Fp3G<P>,
+        r: &G2ProjectiveExtendedGadget<P>,
+    ) -> Result<
+        (
+            G2ProjectiveExtendedGadget<P>,
+            AteAdditionCoefficientsGadget<P>,
+        ),
+        SynthesisError,
+    > {
+        let a = y.square(cs.ns(|| "y^2"))?;
+        let b = r.t.mul(cs.ns(|| "r.t * x"), &x)?;
+        let mut d = r.z.add(cs.ns(|| "r.z + y"), &y)?;
+        d.square_in_place(cs.ns(|| "(r.z + y)^2"))?;
+        d.sub_in_place(cs.ns(|| "(r.z + y)^2 - a"), &a)?;
+        d.sub_in_place(cs.ns(|| "(r.z + y)^2 - a - r.t"), &r.t)?;
+        d.mul_in_place(cs.ns(|| "((r.z + y)^2 - a - r.t) * r.t"), &r.t)?;
+        let h = b.sub(cs.ns(|| "b - r.x"), &r.x)?;
+        let i = h.square(cs.ns(|| "h^2"))?;
+        let e = i.double(cs.ns(|| "2 * i"))?.double(cs.ns(|| "4 * i"))?;
+        let j = h.mul(cs.ns(|| "h * e"), &e)?;
+        let v = r.x.mul(cs.ns(|| "r.x * e"), &e)?;
+        let ry2 = r.y.double(cs.ns(|| "r.y + r.y"))?;
+        let l1 = d.sub(cs.ns(|| "d - (r.y + r.y)"), &ry2)?;
+
+        let v2 = v.double(cs.ns(|| "v + v"))?;
+        let x = l1
+            .square(cs.ns(|| "l1^2"))?
+            .sub(cs.ns(|| "l1^2 - j"), &j)?
+            .sub(cs.ns(|| "l1^2 - j - (v + v)"), &v2)?;
+        let v_minus_x = v.sub(cs.ns(|| "v - x"), &x)?;
+        let j_ry2 = j.mul(cs.ns(|| "j * (r.y + r.y)"), &ry2)?;
+        let y = l1
+            .mul(cs.ns(|| "l1 * (v - x)"), &v_minus_x)?
+            .sub(cs.ns(|| "l1 * (v - x) - (j * (r.y + r.y)"), &j_ry2)?;
+        let mut z = r.z.add(cs.ns(|| "r.z + h"), &h)?;
+        z.square_in_place(cs.ns(|| "(r.z + h)^2"))?;
+        z.sub_in_place(cs.ns(|| "(r.z + h)^2 - r.t"), &r.t)?;
+        z.sub_in_place(cs.ns(|| "(r.z + h)^2 - r.t - i"), &i)?;
+        let t = z.square(cs.ns(|| "z^2"))?;
+
+        let r2 = G2ProjectiveExtendedGadget {
+            x,
+            y,
+            z: z.clone(),
+            t,
+        };
+        let coeff = AteAdditionCoefficientsGadget { c_l1: l1, c_rz: z };
+
+        Ok((r2, coeff))
+    }
+
+    pub fn ate_miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        p: &G1PreparedGadget<P>,
+        q: &G2PreparedGadget<P>,
+    ) -> Result<Fp6G<P>, SynthesisError> {
+        let zero = FpGadget::<P::Fp>::zero(cs.ns(|| "zero"))?;
+        let mut l1_coeff = Fp3G::<P>::new(p.x.clone(), zero.clone(), zero);
+        l1_coeff.sub_in_place(cs.ns(|| "l1_coeff"), &q.x_over_twist)?;
+
+        let mut f = Fp6G::<P>::one(cs.ns(|| "one"))?;
+
+        let mut dbl_idx: usize = 0;
+        let mut add_idx: usize = 0;
+
+        let mut found_one = false;
+
+        for (j, bit) in BitIterator::new(P::ATE_LOOP_COUNT).enumerate() {
+            // code below gets executed for all bits (EXCEPT the MSB itself) of
+            // mnt6_param_p (skipping leading zeros) in MSB to LSB order
+            if !found_one && bit {
+                found_one = true;
+                continue;
+            } else if !found_one {
+                continue;
+            }
+
+            let mut cs = cs.ns(|| format!("bit {}", j));
+
+            let dc = &q.double_coefficients[dbl_idx];
+            dbl_idx += 1;
+
+            let c_j_x_twist = dc.c_j.mul(cs.ns(|| "dc.c_j * p.x_twist"), &p.x_twist)?;
+            let c0 = dc.c_l.sub(cs.ns(|| "-dc.c_4c + dc.c_l"), &dc.c_4c)?.sub(
+                cs.ns(|| "-dc.c_4c - (dc.c_j * p.x_twist) + dc.c_l"),
+                &c_j_x_twist,
+            )?;
+            let c1 = dc.c_h.mul(cs.ns(|| "dc.c_h * p.y_twist"), &p.y_twist)?;
+            let g_rr_at_p = Fp6G::<P>::new(c0, c1);
+
+            f = f
+                .square(cs.ns(|| "f^2"))?
+                .mul(cs.ns(|| "f^2 * g_rr_at_p"), &g_rr_at_p)?;
+
+            if bit {
+                let ac = &q.addition_coefficients[add_idx];
+                add_idx += 1;
+
+                let l1_coeff_c_l1 = l1_coeff.mul(cs.ns(|| "l1_coeff * ac.c_l1"), &ac.c_l1)?;
+                let g_rq_at_p = Fp6G::<P>::new(
+                    ac.c_rz.mul(cs.ns(|| "ac.c_rz * p.y_twist"), &p.y_twist)?,
+                    q.y_over_twist
+                        .mul(cs.ns(|| "q.y_over_twist * ac.c_rz"), &ac.c_rz)?
+                        .add(
+                            cs.ns(|| "q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1)"),
+                            &l1_coeff_c_l1,
+                        )?
+                        .negate(cs.ns(|| "-(q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1))"))?,
+                );
+                f.mul_in_place(cs.ns(|| "f *= g_rq_at_p"), &g_rq_at_p)?;
+            }
+        }
+
+        if P::ATE_IS_LOOP_COUNT_NEG {
+            let ac = &q.addition_coefficients[add_idx];
+
+            let l1_coeff_c_l1 = l1_coeff.mul(cs.ns(|| "l1_coeff * ac.c_l1"), &ac.c_l1)?;
+            let g_rnegr_at_p = Fp6G::<P>::new(
+                ac.c_rz.mul(cs.ns(|| "ac.c_rz * p.y_twist"), &p.y_twist)?,
+                q.y_over_twist
+                    .mul(cs.ns(|| "q.y_over_twist * ac.c_rz"), &ac.c_rz)?
+                    .add(
+                        cs.ns(|| "q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1)"),
+                        &l1_coeff_c_l1,
+                    )?
+                    .negate(cs.ns(|| "-(q.y_over_twist * ac.c_rz + (l1_coeff * ac.c_l1))"))?,
+            );
+            f = f
+                .mul(cs.ns(|| "f * g_rnegr_at_p"), &g_rnegr_at_p)?
+                .inverse(cs.ns(|| "inverse f"))?;
+        }
+
+        Ok(f)
+    }
+
+    pub fn final_exponentiation<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        value: &Fp6G<P>,
+    ) -> Result<GTGadget<P>, SynthesisError> {
+        let value_inv = value.inverse(cs.ns(|| "value inverse"))?;
+        let value_to_first_chunk = Self::final_exponentiation_first_chunk(
+            cs.ns(|| "value_to_first_chunk"),
+            value,
+            &value_inv,
+        )?;
+        let value_inv_to_first_chunk = Self::final_exponentiation_first_chunk(
+            cs.ns(|| "value_inv_to_first_chunk"),
+            &value_inv,
+            value,
+        )?;
+        Self::final_exponentiation_last_chunk(
+            cs.ns(|| "final_exp_last_chunk"),
+            &value_to_first_chunk,
+            &value_inv_to_first_chunk,
+        )
+    }
+
+    fn final_exponentiation_first_chunk<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        elt: &Fp6G<P>,
+        elt_inv: &Fp6G<P>,
+    ) -> Result<Fp6G<P>, SynthesisError> {
+        // (q^3-1)*(q+1)
+
+        // elt_q3 = elt^(q^3)
+        let mut elt_q3 = elt.clone();
+        elt_q3.frobenius_map_in_place(cs.ns(|| "frobenius 3"), 3)?;
+        // elt_q3_over_elt = elt^(q^3-1)
+        let elt_q3_over_elt = elt_q3.mul(cs.ns(|| "elt_q3 * elt_inv"), elt_inv)?;
+        // alpha = elt^((q^3-1) * q)
+        let mut alpha = elt_q3_over_elt.clone();
+        alpha.frobenius_map_in_place(cs.ns(|| "frobenius 1"), 1)?;
+        // beta = elt^((q^3-1)*(q+1)
+        alpha.mul(cs.ns(|| "alpha * elt_q3_over_elt"), &elt_q3_over_elt)
+    }
+
+    fn final_exponentiation_last_chunk<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        elt: &Fp6G<P>,
+        elt_inv: &Fp6G<P>,
+    ) -> Result<Fp6G<P>, SynthesisError> {
+        let elt_clone = elt.clone();
+        let elt_inv_clone = elt_inv.clone();
+
+        let mut elt_q = elt.clone();
+        elt_q.frobenius_map_in_place(cs.ns(|| "frobenius 1"), 1)?;
+
+        let w1_part = elt_q.cyclotomic_exp(cs.ns(|| "w1_part"), &P::FINAL_EXPONENT_LAST_CHUNK_1)?;
+        let w0_part;
+        if P::FINAL_EXPONENT_LAST_CHUNK_W0_IS_NEG {
+            w0_part = elt_inv_clone
+                .cyclotomic_exp(cs.ns(|| "w0_part"), &P::FINAL_EXPONENT_LAST_CHUNK_ABS_OF_W0)?;
+        } else {
+            w0_part = elt_clone
+                .cyclotomic_exp(cs.ns(|| "w0_part"), &P::FINAL_EXPONENT_LAST_CHUNK_ABS_OF_W0)?;
+        }
+
+        w1_part.mul(cs.ns(|| "w1_part * w0_part"), &w0_part)
+    }
+}
+
+impl<P: MNT6Parameters> PG<MNT6<P>, P::Fp> for PairingGadget<P> {
+    type G1Gadget = G1Gadget<P>;
+    type G2Gadget = G2Gadget<P>;
+    type G1PreparedGadget = G1PreparedGadget<P>;
+    type G2PreparedGadget = G2PreparedGadget<P>;
+    type GTGadget = GTGadget<P>;
+
+    fn miller_loop<CS: ConstraintSystem<P::Fp>>(
+        mut cs: CS,
+        ps: &[Self::G1PreparedGadget],
+        qs: &[Self::G2PreparedGadget],
+    ) -> Result<Self::GTGadget, SynthesisError> {
+        let mut result = Fp6G::<P>::one(cs.ns(|| "one"))?;
+        for (i, (p, q)) in ps.iter().zip(qs.iter()).enumerate() {
+            let miller =
+                Self::ate_miller_loop(cs.ns(|| format!("ate miller loop iteration {}", i)), p, q)?;
+            result.mul_in_place(
+                cs.ns(|| format!("mul ate miller loop iteration {}", i)),
+                &miller,
+            )?;
+        }
+
+        Ok(result)
+    }
+
+    fn final_exponentiation<CS: ConstraintSystem<P::Fp>>(
+        cs: CS,
+        r: &Self::GTGadget,
+    ) -> Result<Self::GTGadget, SynthesisError> {
+        Self::final_exponentiation(cs, r)
+    }
+
+    fn prepare_g1<CS: ConstraintSystem<P::Fp>>(
+        cs: CS,
+        p: &Self::G1Gadget,
+    ) -> Result<Self::G1PreparedGadget, SynthesisError> {
+        Self::G1PreparedGadget::from_affine(cs, p)
+    }
+
+    fn prepare_g2<CS: ConstraintSystem<P::Fp>>(
+        cs: CS,
+        q: &Self::G2Gadget,
+    ) -> Result<Self::G2PreparedGadget, SynthesisError> {
+        Self::G2PreparedGadget::from_affine(cs, q)
+    }
+}

r1cs-std/src/pairing/mod.rs

@@ -4,6 +4,8 @@ use core::fmt::Debug;
 use r1cs_core::{ConstraintSystem, SynthesisError};
 
 pub mod bls12;
+pub mod mnt4;
+pub mod mnt6;
 
 pub trait PairingGadget<PairingE: PairingEngine, ConstraintF: Field> {
     type G1Gadget: GroupGadget<PairingE::G1Projective, ConstraintF>;