diff --git a/book/src/zksnark/subsec2.md b/book/src/zksnark/subsec2.md
index 04562ec..361f24c 100644
--- a/book/src/zksnark/subsec2.md
+++ b/book/src/zksnark/subsec2.md
@@ -145,6 +145,50 @@ O &= \begin{bmatrix}
 \end{bmatrix}
 \end{align*}
 
+**Implementation:**
+
+```rust
+fn dot<F: Field>(a: &Vec<F>, b: &Vec<F>) -> F {
+    let mut result = F::zero();
+    for (a_i, b_i) in a.iter().zip(b.iter()) {
+        result = result + a_i.clone() * b_i.clone();
+    }
+    result
+}
+
+#[derive(Debug, Clone)]
+pub struct R1CS<F: Field> {
+    pub left: Vec<Vec<F>>,
+    pub right: Vec<Vec<F>>,
+    pub out: Vec<Vec<F>>,
+    pub m: usize,
+    pub d: usize,
+}
+
+impl<F: Field> R1CS<F> {
+    pub fn new(left: Vec<Vec<F>>, right: Vec<Vec<F>>, out: Vec<Vec<F>>) -> Self {
+        let d = left.len();
+        let m = if d == 0 { 0 } else { left[0].len() };
+        R1CS {
+            left,
+            right,
+            out,
+            m,
+            d,
+        }
+    }
+
+    pub fn is_satisfied(&self, a: &Vec<F>) -> bool {
+        let zero = F::zero();
+        self.left
+            .iter()
+            .zip(self.right.iter())
+            .zip(self.out.iter())
+            .all(|((l, r), o)| dot(&l, &a) * dot(&r, &a) - dot(&o, &a) == zero)
+    }
+}
+```
+
 ## Quadratic Arithmetic Program (QAP)
 
 Recall that the prover aims to demonstrate knowledge of a witness \\(w\\) without revealing it. This is equivalent to knowing a vector \\(a\\) that satisfies \\((L \cdot a) \circ (R \cdot a) = O \cdot a\\), where \\(\circ\\) denotes the Hadamard (element-wise) product. However, evaluating this equivalence directly requires \\(\Omega(d)\\) operations, where \\(d\\) is the number of rows. To improve efficiency, we can convert this matrix comparison to a polynomial comparison, leveraging the Schwartz-Zippel Lemma, which allows us to check polynomial equality with \\(\Omega(1)\\) evaluations.
@@ -210,3 +254,57 @@ To address this discrepancy, we introduce a degree \\(d\\) polynomial \\(t(x) =
 \end{equation}
 
 where \\(h(x) = \frac{\ell(x) \cdot r(x) - o(x)}{t(x)}\\). This formulation allows us to maintain the desired polynomial relationships while accounting for the degree differences.
+
+**Implementation:**
+
+```rust
+#[derive(Debug, Clone)]
+pub struct QAP<'a, F: Field> {
+    pub r1cs: &'a R1CS<F>,
+    pub t: Polynomial<F>,
+}
+
+impl<'a, F: Field> QAP<'a, F> {
+    fn new(r1cs: &'a R1CS<F>) -> Self {
+        QAP {
+            r1cs: r1cs,
+            t: Polynomial::<F>::from_monomials(
+                &(1..=r1cs.d).map(|i| F::from_value(i)).collect::<Vec<F>>(),
+            ),
+        }
+    }
+
+    fn generate_polynomials(&self, a: &Vec<F>) -> (Polynomial<F>, Polynomial<F>, Polynomial<F>) {
+        let left_dot_products = self
+            .r1cs
+            .left
+            .iter()
+            .map(|v| dot(&v, &a))
+            .collect::<Vec<F>>();
+        let right_dot_products = self
+            .r1cs
+            .right
+            .iter()
+            .map(|v| dot(&v, &a))
+            .collect::<Vec<F>>();
+        let out_dot_products = self
+            .r1cs
+            .out
+            .iter()
+            .map(|v| dot(&v, &a))
+            .collect::<Vec<F>>();
+
+        let x = (1..=self.r1cs.m)
+            .map(|i| F::from_value(i))
+            .collect::<Vec<F>>();
+        let left_interpolated_polynomial = Polynomial::<F>::interpolate(&x, &left_dot_products);
+        let right_interpolated_polynomial = Polynomial::<F>::interpolate(&x, &right_dot_products);
+        let out_interpolated_polynomial = Polynomial::<F>::interpolate(&x, &out_dot_products);
+        (
+            left_interpolated_polynomial,
+            right_interpolated_polynomial,
+            out_interpolated_polynomial,
+        )
+    }
+}
+```
diff --git a/myzkp/src/modules/mod.rs b/myzkp/src/modules/mod.rs
index c9e6315..0ac1676 100644
--- a/myzkp/src/modules/mod.rs
+++ b/myzkp/src/modules/mod.rs
@@ -4,6 +4,8 @@ pub mod educational_protocols;
 pub mod efield;
 pub mod field;
 pub mod polynomial;
+pub mod qap;
+pub mod r1cs;
 pub mod ring;
 pub mod utils;
 // pub mod snark;
diff --git a/myzkp/src/modules/qap.rs b/myzkp/src/modules/qap.rs
new file mode 100644
index 0000000..01bfb75
--- /dev/null
+++ b/myzkp/src/modules/qap.rs
@@ -0,0 +1,84 @@
+use num_traits::One;
+use num_traits::Zero;
+
+use crate::modules::field::Field;
+use crate::modules::polynomial::Polynomial;
+use crate::modules::r1cs::{dot, R1CS};
+
+#[derive(Debug, Clone)]
+pub struct QAP<'a, F: Field> {
+    pub r1cs: &'a R1CS<F>,
+    pub t: Polynomial<F>,
+}
+
+impl<'a, F: Field> QAP<'a, F> {
+    fn new(r1cs: &'a R1CS<F>) -> Self {
+        QAP {
+            r1cs: r1cs,
+            t: Polynomial::<F>::from_monomials(
+                &(1..=r1cs.d).map(|i| F::from_value(i)).collect::<Vec<F>>(),
+            ),
+        }
+    }
+
+    fn generate_polynomials(&self, a: &Vec<F>) -> (Polynomial<F>, Polynomial<F>, Polynomial<F>) {
+        let left_dot_products = self
+            .r1cs
+            .left
+            .iter()
+            .map(|v| dot(&v, &a))
+            .collect::<Vec<F>>();
+        let right_dot_products = self
+            .r1cs
+            .right
+            .iter()
+            .map(|v| dot(&v, &a))
+            .collect::<Vec<F>>();
+        let out_dot_products = self
+            .r1cs
+            .out
+            .iter()
+            .map(|v| dot(&v, &a))
+            .collect::<Vec<F>>();
+
+        let x = (1..=self.r1cs.m)
+            .map(|i| F::from_value(i))
+            .collect::<Vec<F>>();
+        let left_interpolated_polynomial = Polynomial::<F>::interpolate(&x, &left_dot_products);
+        let right_interpolated_polynomial = Polynomial::<F>::interpolate(&x, &right_dot_products);
+        let out_interpolated_polynomial = Polynomial::<F>::interpolate(&x, &out_dot_products);
+        (
+            left_interpolated_polynomial,
+            right_interpolated_polynomial,
+            out_interpolated_polynomial,
+        )
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use crate::modules::field::{FiniteFieldElement, ModEIP197};
+    use crate::modules::ring::Ring;
+
+    type F = FiniteFieldElement<ModEIP197>;
+
+    #[test]
+    fn test_qap_single_multiplication() {
+        // z = x * y
+        // (1, z, x, y) = (1, 3690, 82, 45)
+        let left = vec![vec![F::zero(), F::zero(), F::one(), F::zero()]];
+        let right = vec![vec![F::zero(), F::zero(), F::zero(), F::one()]];
+        let out = vec![vec![F::zero(), F::one(), F::zero(), F::zero()]];
+        let a = vec![
+            F::one(),
+            F::from_value(3690),
+            F::from_value(82),
+            F::from_value(45),
+        ];
+        let r1cs = R1CS::new(left, right, out);
+        let qap = QAP::new(&r1cs);
+        let (_, _, _) = qap.generate_polynomials(&a);
+    }
+}
diff --git a/myzkp/src/modules/r1cs.rs b/myzkp/src/modules/r1cs.rs
new file mode 100644
index 0000000..1299d44
--- /dev/null
+++ b/myzkp/src/modules/r1cs.rs
@@ -0,0 +1,256 @@
+use num_traits::One;
+use num_traits::Zero;
+
+use crate::modules::field::Field;
+
+pub fn dot<F: Field>(a: &Vec<F>, b: &Vec<F>) -> F {
+    let mut result = F::zero();
+    for (a_i, b_i) in a.iter().zip(b.iter()) {
+        result = result + a_i.clone() * b_i.clone();
+    }
+    result
+}
+
+#[derive(Debug, Clone)]
+pub struct R1CS<F: Field> {
+    pub left: Vec<Vec<F>>,
+    pub right: Vec<Vec<F>>,
+    pub out: Vec<Vec<F>>,
+    pub m: usize,
+    pub d: usize,
+}
+
+impl<F: Field> R1CS<F> {
+    pub fn new(left: Vec<Vec<F>>, right: Vec<Vec<F>>, out: Vec<Vec<F>>) -> Self {
+        let m = left.len();
+        let d = if m == 0 { 0 } else { left[0].len() };
+        R1CS {
+            left,
+            right,
+            out,
+            m,
+            d,
+        }
+    }
+
+    pub fn is_satisfied(&self, a: &Vec<F>) -> bool {
+        let zero = F::zero();
+        self.left
+            .iter()
+            .zip(self.right.iter())
+            .zip(self.out.iter())
+            .all(|((l, r), o)| dot(&l, &a) * dot(&r, &a) - dot(&o, &a) == zero)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use crate::modules::field::{FiniteFieldElement, ModEIP197};
+    use crate::modules::ring::Ring;
+
+    type F = FiniteFieldElement<ModEIP197>;
+
+    #[test]
+    fn test_r1cs_single_multiplication() {
+        // z = x * y
+        // (1, z, x, y) = (1, 3690, 82, 45)
+        let left = vec![vec![F::zero(), F::zero(), F::one(), F::zero()]];
+        let right = vec![vec![F::zero(), F::zero(), F::zero(), F::one()]];
+        let out = vec![vec![F::zero(), F::one(), F::zero(), F::zero()]];
+        let a = vec![
+            F::one(),
+            F::from_value(3690),
+            F::from_value(82),
+            F::from_value(45),
+        ];
+        let a_prime = vec![
+            F::one(),
+            F::from_value(3690),
+            F::from_value(82),
+            F::from_value(46),
+        ];
+        let r1cs = R1CS::new(left, right, out);
+        assert!(r1cs.is_satisfied(&a));
+        assert!(!r1cs.is_satisfied(&a_prime));
+    }
+
+    #[test]
+    fn test_r1cs_addition_with_constant() {
+        // z = x * y + 3
+        // (1, z, x, y) = (1, 9, 2, 3)
+        let left = vec![vec![F::zero(), F::zero(), F::one(), F::zero()]];
+        let right = vec![vec![F::zero(), F::zero(), F::zero(), F::one()]];
+        let out = vec![vec![F::from_value(-3_i32), F::one(), F::zero(), F::zero()]];
+        let a = vec![
+            F::one(),
+            F::from_value(9),
+            F::from_value(2),
+            F::from_value(3),
+        ];
+        let a_prime = vec![
+            F::one(),
+            F::from_value(9),
+            F::from_value(2),
+            F::from_value(4),
+        ];
+        let r1cs = R1CS::new(left, right, out);
+        assert!(r1cs.is_satisfied(&a));
+        assert!(!r1cs.is_satisfied(&a_prime));
+    }
+
+    #[test]
+    fn test_r1cs_multiplication_with_constant() {
+        // z = 3x^2 + y
+        // (1, z, x, y) = (1, 17, 2, 5)
+        let left = vec![vec![F::zero(), F::zero(), F::from_value(3_u32), F::zero()]];
+        let right = vec![vec![F::zero(), F::zero(), F::one(), F::zero()]];
+        let out = vec![vec![F::zero(), F::one(), F::zero(), F::from_value(-1_i32)]];
+        let a = vec![
+            F::one(),
+            F::from_value(17),
+            F::from_value(2),
+            F::from_value(5),
+        ];
+        let a_prime = vec![
+            F::one(),
+            F::from_value(17),
+            F::from_value(2),
+            F::from_value(7),
+        ];
+        let r1cs = R1CS::new(left, right, out);
+        assert!(r1cs.is_satisfied(&a));
+        assert!(!r1cs.is_satisfied(&a_prime));
+    }
+
+    #[test]
+    fn test_r1cs_multiple_constraints() {
+        // v = a * b
+        // u = c * d
+        // r = v * u
+        // (1, r, a, b, c, d, v, u) = (1, 210, 2, 3, 5, 7, 6, 35)
+        let left = vec![
+            vec![
+                F::zero(),
+                F::zero(),
+                F::one(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+            ],
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+            ],
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+                F::zero(),
+            ],
+        ];
+
+        let right = vec![
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+            ],
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+                F::zero(),
+                F::zero(),
+            ],
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+            ],
+        ];
+
+        let out = vec![
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+                F::zero(),
+            ],
+            vec![
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::one(),
+            ],
+            vec![
+                F::zero(),
+                F::one(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+                F::zero(),
+            ],
+        ];
+
+        let a = vec![
+            F::one(),
+            F::from_value(210),
+            F::from_value(2),
+            F::from_value(3),
+            F::from_value(5),
+            F::from_value(7),
+            F::from_value(6),
+            F::from_value(35),
+        ];
+        let a_prime = vec![
+            F::one(),
+            F::from_value(211),
+            F::from_value(2),
+            F::from_value(3),
+            F::from_value(5),
+            F::from_value(7),
+            F::from_value(6),
+            F::from_value(35),
+        ];
+        let r1cs = R1CS::new(left, right, out);
+        assert!(r1cs.is_satisfied(&a));
+        assert!(!r1cs.is_satisfied(&a_prime));
+    }
+}