diff --git a/synir/src/data_structures.rs b/synir/src/data_structures.rs
index 98aade8b..3d19249b 100644
--- a/synir/src/data_structures.rs
+++ b/synir/src/data_structures.rs
@@ -1,11 +1,15 @@
 use crate::IndexType;
+use bitvec::vec::BitVec;
 
+pub mod angle;
 mod clifford_tableau;
+mod pauli_exponential;
 mod pauli_polynomial;
 mod pauli_string;
 
-use bitvec::vec::BitVec;
+pub use angle::Angle;
 pub use clifford_tableau::CliffordTableau;
+pub use pauli_exponential::PauliExponential;
 pub use pauli_polynomial::PauliPolynomial;
 pub use pauli_string::PauliString;
 
diff --git a/synir/src/data_structures/angle.rs b/synir/src/data_structures/angle.rs
new file mode 100644
index 00000000..703e1875
--- /dev/null
+++ b/synir/src/data_structures/angle.rs
@@ -0,0 +1,292 @@
+use std::{
+    f64::consts::PI,
+    ops::{Add, AddAssign, Sub, SubAssign},
+};
+
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub enum Angle {
+    Arbitrary(f64),
+    Pi4Rotations(u8),
+}
+
+impl Angle {
+    pub fn from_angle(rad: f64) -> Self {
+        Angle::Arbitrary(rad)
+    }
+
+    pub fn from_angles(angles: &[f64]) -> Vec<Self> {
+        angles
+            .into_iter()
+            .map(|rad| Angle::from_angle(*rad))
+            .collect()
+    }
+
+    pub fn from_pi4_rotation(n: u8) -> Self {
+        Angle::Pi4Rotations(n % 8)
+    }
+
+    pub fn from_pi4_rotations(ns: &[u8]) -> Vec<Self> {
+        ns.into_iter()
+            .map(|n| Angle::from_pi4_rotation(*n))
+            .collect()
+    }
+
+    pub fn to_radians(&self) -> f64 {
+        match self {
+            Angle::Arbitrary(rad) => *rad,
+            Angle::Pi4Rotations(n) => (*n as f64) * (std::f64::consts::FRAC_PI_4),
+        }
+    }
+
+    pub fn flip(&mut self) {
+        match self {
+            Angle::Arbitrary(rad) => *rad = -*rad,
+            Angle::Pi4Rotations(n) => *n = (8 - *n) % 8,
+        }
+    }
+}
+
+impl AddAssign for Angle {
+    fn add_assign(&mut self, other: Self) {
+        match (self, other) {
+            (Angle::Arbitrary(rad1), Angle::Arbitrary(rad2)) => {
+                *rad1 += rad2;
+            }
+            (Angle::Pi4Rotations(n1), Angle::Pi4Rotations(n2)) => {
+                *n1 = (*n1 + n2) % 8;
+            }
+            (Angle::Arbitrary(rad1), Angle::Pi4Rotations(n2)) => {
+                *rad1 += n2 as f64 * PI / 4.0;
+            }
+            _ => panic!("Cannot add Arbitrary Angle to Pi4 rotation"),
+        }
+    }
+}
+
+impl SubAssign for Angle {
+    fn sub_assign(&mut self, other: Self) {
+        match (self, other) {
+            (Angle::Arbitrary(rad1), Angle::Arbitrary(rad2)) => {
+                *rad1 -= rad2;
+            }
+            (Angle::Pi4Rotations(n1), Angle::Pi4Rotations(n2)) => {
+                *n1 = (*n1 + (8 - n2)) % 8;
+            }
+            (Angle::Arbitrary(rad1), Angle::Pi4Rotations(n2)) => {
+                *rad1 -= n2 as f64 * PI / 4.0;
+            }
+            _ => panic!("Cannot subtract different types of Angles"),
+        }
+    }
+}
+
+impl Add for Angle {
+    type Output = Angle;
+
+    fn add(self, other: Angle) -> Angle {
+        match (self, other) {
+            (Angle::Arbitrary(rad1), Angle::Arbitrary(rad2)) => Angle::Arbitrary(rad1 + rad2),
+            (Angle::Pi4Rotations(n1), Angle::Pi4Rotations(n2)) => {
+                Angle::Pi4Rotations((n1 + n2) % 8)
+            }
+            (Angle::Arbitrary(rad1), Angle::Pi4Rotations(n2)) => {
+                Angle::Arbitrary(rad1 + n2 as f64 * PI / 4.0)
+            }
+            (Angle::Pi4Rotations(n1), Angle::Arbitrary(rad2)) => {
+                Angle::Arbitrary(rad2 + n1 as f64 * PI / 4.0)
+            }
+        }
+    }
+}
+
+impl Sub for Angle {
+    type Output = Angle;
+
+    fn sub(self, other: Angle) -> Angle {
+        match (self, other) {
+            (Angle::Arbitrary(rad1), Angle::Arbitrary(rad2)) => Angle::Arbitrary(rad1 - rad2),
+            (Angle::Pi4Rotations(n1), Angle::Pi4Rotations(n2)) => {
+                Angle::Pi4Rotations((n1 + 8 - n2) % 8)
+            }
+            (Angle::Arbitrary(rad1), Angle::Pi4Rotations(n2)) => {
+                Angle::Arbitrary(rad1 - n2 as f64 * PI / 4.0)
+            }
+            (Angle::Pi4Rotations(n1), Angle::Arbitrary(rad2)) => {
+                Angle::Arbitrary(n1 as f64 * PI / 4.0 - rad2)
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn check_angle_approx(angle1: Angle, angle2: Angle) -> bool {
+        match (angle1, angle2) {
+            (Angle::Arbitrary(a1), Angle::Arbitrary(a2)) => (a1 - a2).abs() < 1e-9,
+            (Angle::Pi4Rotations(a1), Angle::Pi4Rotations(a2)) => a1 == a2,
+            _ => panic!("Not defined for Arbitrary Angles and Pi4 rotations"),
+        }
+    }
+
+    #[test]
+    fn test_angle_simple_add() {
+        let n1 = 1;
+        let n2 = 2;
+
+        let a1 = Angle::from_pi4_rotation(n1);
+        let a2 = Angle::from_pi4_rotation(n2);
+
+        assert_eq!(a1 + a2, Angle::from_pi4_rotation(3));
+
+        let mut a3 = Angle::from_pi4_rotation(n1);
+        a3 += a2;
+
+        assert_eq!(a3, Angle::from_pi4_rotation(3));
+    }
+
+    #[test]
+    fn test_angle_overflow_add() {
+        let n1 = 5;
+        let n2 = 6;
+
+        let a1 = Angle::from_pi4_rotation(n1);
+        let a2 = Angle::from_pi4_rotation(n2);
+
+        assert_eq!(a1 + a2, Angle::from_pi4_rotation(3));
+
+        let mut a3 = Angle::from_pi4_rotation(n1);
+        a3 += a2;
+
+        assert_eq!(a3, Angle::from_pi4_rotation(3));
+    }
+
+    #[test]
+    fn test_angle_simple_sub() {
+        let n1 = 4;
+        let n2 = 2;
+
+        let a1 = Angle::from_pi4_rotation(n1);
+        let a2 = Angle::from_pi4_rotation(n2);
+
+        assert_eq!(a1 - a2, Angle::from_pi4_rotation(2));
+
+        let mut a3 = Angle::from_pi4_rotation(n1);
+        a3 -= a2;
+
+        assert_eq!(a3, Angle::from_pi4_rotation(2));
+    }
+
+    #[test]
+    fn test_angle_overflow_sub() {
+        let n1 = 2;
+        let n2 = 6;
+
+        let a1 = Angle::from_pi4_rotation(n1);
+        let a2 = Angle::from_pi4_rotation(n2);
+
+        assert_eq!(a1 - a2, Angle::from_pi4_rotation(4));
+
+        let mut a3 = Angle::from_pi4_rotation(n1);
+        a3 -= a2;
+
+        assert_eq!(a3, Angle::from_pi4_rotation(4));
+    }
+
+    #[test]
+    fn test_angle_float_simple_add() {
+        let n1 = 0.32;
+        let n2 = 0.64;
+
+        let a1 = Angle::from_angle(n1);
+        let a2 = Angle::from_angle(n2);
+
+        let ref_a = Angle::from_angle(0.96);
+
+        assert!(check_angle_approx(a1 + a2, ref_a));
+
+        let mut a3 = Angle::from_angle(n1);
+        a3 += a2;
+
+        assert!(check_angle_approx(a3, ref_a));
+    }
+
+    #[test]
+    fn test_angle_float_simple_sub() {
+        let n1 = 0.32;
+        let n2 = 0.64;
+
+        let a1 = Angle::from_angle(n1);
+        let a2 = Angle::from_angle(n2);
+
+        let ref_a = Angle::from_angle(-0.32);
+
+        assert!(check_angle_approx(a1 - a2, ref_a));
+
+        let mut a3 = Angle::from_angle(n1);
+        a3 -= a2;
+
+        assert!(check_angle_approx(a3, ref_a));
+    }
+
+    #[test]
+    fn test_angle_mixed_simple_add() {
+        let n1 = 0.32;
+        let n2 = 2;
+
+        let a1 = Angle::from_angle(n1);
+        let a2 = Angle::from_pi4_rotation(n2);
+
+        let ref_a = Angle::from_angle(1.8907963268);
+        assert!(check_angle_approx(a1 + a2, ref_a));
+        assert!(check_angle_approx(a2 + a1, ref_a));
+
+        let mut a3 = Angle::from_angle(n1);
+        a3 += a2;
+
+        assert!(check_angle_approx(a3, ref_a));
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_angle_bad_mixed_simple_add() {
+        let n1 = 0.32;
+        let n2 = 2;
+
+        let mut a2 = Angle::from_pi4_rotation(n2);
+        let a3 = Angle::from_angle(n1);
+        a2 += a3
+    }
+
+    #[test]
+    fn test_angle_mixed_simple_sub() {
+        let n1 = 0.32;
+        let n2 = 2;
+
+        let a1 = Angle::from_angle(n1);
+        let a2 = Angle::from_pi4_rotation(n2);
+
+        let ref_a1 = Angle::from_angle(-1.2507963268);
+        let ref_a2 = Angle::from_angle(1.2507963268);
+
+        assert!(check_angle_approx(a1 - a2, ref_a1));
+        assert!(check_angle_approx(a2 - a1, ref_a2));
+
+        let mut a3 = Angle::from_angle(n1);
+        a3 -= a2;
+
+        assert!(check_angle_approx(a3, ref_a1));
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_angle_bad_mixed_simple_sub() {
+        let n1 = 0.32;
+        let n2 = 2;
+
+        let mut a2 = Angle::from_pi4_rotation(n2);
+        let a3 = Angle::from_angle(n1);
+        a2 -= a3
+    }
+}
diff --git a/synir/src/data_structures/clifford_tableau.rs b/synir/src/data_structures/clifford_tableau.rs
index a5d34cc5..1c031a37 100644
--- a/synir/src/data_structures/clifford_tableau.rs
+++ b/synir/src/data_structures/clifford_tableau.rs
@@ -10,7 +10,7 @@ use super::{
     pauli_string::{cx, PauliString},
     IndexType, PropagateClifford,
 };
-use crate::data_structures::PauliLetter;
+use crate::data_structures::{Angle, PauliLetter};
 
 #[derive(PartialEq, Eq, Debug, Clone, Default)]
 pub struct CliffordTableau {
@@ -191,6 +191,90 @@ impl CliffordTableau {
         }
     }
 
+    /// Composes a gadget onto a Clifford tableau if the angle is Clifford
+    /// Decomposes the Pauli gadget by performing naive decomposition into mapping to Z legs, CNOT walls and Z-rotations
+    pub fn compose_gadget(&mut self, rhs: (PauliString, Angle)) -> Result<(), String> {
+        let (pauli_string, angle) = rhs;
+        let size = self.size();
+        assert_eq!(
+            size,
+            pauli_string.len(),
+            "Cannot compose Clifford tableau with PauliPolynomial of different size"
+        );
+        let pi2rotations = match angle {
+            Angle::Arbitrary(angle) => panic!(
+                "Cannot compose Clifford tableau with non-Clifford angle: {}",
+                angle
+            ),
+            Angle::Pi4Rotations(rotations) => {
+                if rotations % 2 == 1 {
+                    panic!("Cannot compose Clifford tableau with non-Clifford angle: {} pi/4 rotations", rotations);
+                }
+                (rotations >> 1) % 4
+            }
+        };
+        let mut leg_numbers = Vec::with_capacity(size);
+        for i in 0..size {
+            match pauli_string.pauli(i) {
+                PauliLetter::I => {}
+                PauliLetter::X => {
+                    self.h(i);
+                    leg_numbers.push(i);
+                }
+                PauliLetter::Y => {
+                    self.v(i);
+                    leg_numbers.push(i);
+                }
+                PauliLetter::Z => {
+                    leg_numbers.push(i);
+                }
+            }
+        }
+
+        for (control, target) in leg_numbers.iter().tuple_windows() {
+            self.cx(*control, *target);
+        }
+        match pi2rotations {
+            0 => {}
+            1 => {
+                let target = *leg_numbers.last().unwrap();
+                self.s(target);
+            }
+            2 => {
+                let target = *leg_numbers.last().unwrap();
+                self.z(target);
+            }
+            3 => {
+                let target = *leg_numbers.last().unwrap();
+                self.s_dgr(target);
+            }
+            _ => unreachable!(),
+        }
+
+        for (control, target) in leg_numbers
+            .iter()
+            .tuple_windows()
+            .collect_vec()
+            .iter()
+            .rev()
+        {
+            self.cx(**control, **target);
+        }
+        for i in 0..size {
+            match pauli_string.pauli(i) {
+                PauliLetter::I => {}
+                PauliLetter::X => {
+                    self.h(i);
+                }
+                PauliLetter::Y => {
+                    self.v_dgr(i);
+                }
+                PauliLetter::Z => {}
+            }
+        }
+        Ok(())
+    }
+
     pub fn permute(&mut self, permutation_vector: Vec<usize>) {
         assert_eq!(
             permutation_vector
diff --git a/synir/src/data_structures/pauli_exponential.rs b/synir/src/data_structures/pauli_exponential.rs
new file mode 100644
index 00000000..84924dc1
--- /dev/null
+++ b/synir/src/data_structures/pauli_exponential.rs
@@ -0,0 +1,49 @@
+use std::collections::VecDeque;
+
+use crate::data_structures::{CliffordTableau, PauliPolynomial, PropagateClifford};
+
+#[derive(Default)]
+pub struct PauliExponential {
+    pub(crate) pauli_polynomials: VecDeque<PauliPolynomial>,
+    pub(crate) clifford_tableau: CliffordTableau,
+}
+
+impl PauliExponential {
+    pub fn new(
+        pauli_polynomials: VecDeque<PauliPolynomial>,
+        clifford_tableau: CliffordTableau,
+    ) -> Self {
+        PauliExponential {
+            pauli_polynomials,
+            clifford_tableau,
+        }
+    }
+
+    pub fn chains(&mut self) -> &mut VecDeque<PauliPolynomial> {
+        &mut self.pauli_polynomials
+    }
+
+    pub fn size(&self) -> usize {
+        self.clifford_tableau.size()
+    }
+}
+
+impl PropagateClifford for PauliExponential {
+    fn cx(&mut self, control: crate::IndexType, target: crate::IndexType) -> &mut Self {
+        self.pauli_polynomials.cx(control, target);
+        self.clifford_tableau.cx(control, target);
+        self
+    }
+
+    fn s(&mut self, target: crate::IndexType) -> &mut Self {
+        self.pauli_polynomials.s(target);
+        self.clifford_tableau.s(target);
+        self
+    }
+
+    fn v(&mut self, target: crate::IndexType) -> &mut Self {
+        self.pauli_polynomials.v(target);
+        self.clifford_tableau.v(target);
+        self
+    }
+}
diff --git a/synir/src/data_structures/pauli_polynomial.rs b/synir/src/data_structures/pauli_polynomial.rs
index ff38897b..4a5804e3 100644
--- a/synir/src/data_structures/pauli_polynomial.rs
+++ b/synir/src/data_structures/pauli_polynomial.rs
@@ -1,12 +1,11 @@
 use std::iter::zip;
 
+use super::{pauli_string::PauliString, IndexType, MaskedPropagateClifford, PropagateClifford};
+use crate::data_structures::{Angle, PauliLetter};
 use bitvec::vec::BitVec;
 use itertools::zip_eq;
 
-use super::{pauli_string::PauliString, IndexType, MaskedPropagateClifford, PropagateClifford};
-
-// todo: Make this into a union / type Angle
-type Angle = f64;
+mod simplify;
 
 #[derive(Debug, Clone, Default)]
 pub struct PauliPolynomial {
@@ -62,6 +61,47 @@ impl PauliPolynomial {
     pub fn angle(&self, i: usize) -> Angle {
         self.angles[i]
     }
+
+    pub fn extend_z(&mut self, target: usize, angle: f64) {
+        for (i, chain) in self.chains.iter_mut().enumerate() {
+            if i == target {
+                chain.z.push(true);
+            } else {
+                chain.z.push(false);
+            }
+            chain.x.push(false);
+        }
+        self.angles.push(Angle::Arbitrary(angle));
+    }
+
+    pub fn commutes_with(&self, other: &PauliPolynomial) -> bool {
+        let size = self.size();
+        assert_eq!(size, other.size());
+
+        let self_length = self.length();
+        let other_length = other.length();
+
+        for index_1 in 0..self_length {
+            let mut pauli_string = Vec::with_capacity(size);
+            for q1 in 0..size {
+                pauli_string.push(self.chain(q1).pauli(index_1));
+            }
+            for index_2 in 0..other_length {
+                let other_pauli_string = (0..size).map(|q2| other.chain(q2).pauli(index_2));
+                let mut commutes = true;
+                for (p1, p2) in zip(&pauli_string, other_pauli_string) {
+                    if *p1 == PauliLetter::I || p2 == PauliLetter::I || p1 == &p2 {
+                        continue;
+                    }
+                    commutes = !commutes;
+                }
+                if !commutes {
+                    return false;
+                }
+            }
+        }
+        true
+    }
 }
 
 impl PropagateClifford for PauliPolynomial {
@@ -78,7 +118,7 @@ impl PropagateClifford for PauliPolynomial {
         super::pauli_string::cx(control, target);
         for (angle, flip) in zip(self.angles.iter_mut(), bit_mask.iter()) {
             if *flip {
-                *angle *= -1.0;
+                angle.flip();
             }
         }
 
@@ -91,7 +131,7 @@ impl PropagateClifford for PauliPolynomial {
         let y_vec = chains_target.y_bitmask();
         for (angle, flip) in zip(self.angles.iter_mut(), y_vec.iter()) {
             if *flip {
-                *angle *= -1.0;
+                angle.flip();
             }
         }
         chains_target.s();
@@ -105,7 +145,7 @@ impl PropagateClifford for PauliPolynomial {
         let y_vec = chains_target.y_bitmask();
         for (angle, flip) in zip(self.angles.iter_mut(), y_vec.iter()) {
             if *flip {
-                *angle *= -1.0;
+                angle.flip();
             }
         }
         self
@@ -126,7 +166,7 @@ impl MaskedPropagateClifford for PauliPolynomial {
         super::pauli_string::masked_cx(control, target, mask);
         for (angle, flip) in zip(self.angles.iter_mut(), bit_mask.iter()) {
             if *flip {
-                *angle *= -1.0;
+                angle.flip();
             }
         }
 
@@ -140,7 +180,7 @@ impl MaskedPropagateClifford for PauliPolynomial {
         let y_vec = chains_target.masked_y_bitmask(mask);
         for (angle, flip) in zip(self.angles.iter_mut(), y_vec.iter()) {
             if *flip {
-                *angle *= -1.0;
+                angle.flip();
             }
         }
         chains_target.masked_s(mask);
@@ -154,15 +194,17 @@ impl MaskedPropagateClifford for PauliPolynomial {
         let y_vec = chains_target.masked_y_bitmask(mask);
         for (angle, flip) in zip(self.angles.iter_mut(), y_vec.iter()) {
             if *flip {
-                *angle *= -1.0;
+                angle.flip();
             }
         }
         self
     }
 }
+
 #[cfg(test)]
 mod tests {
     use super::*;
+    use itertools::Itertools;
 
     impl PartialEq for PauliPolynomial {
         fn eq(&self, other: &Self) -> bool {
@@ -173,7 +215,11 @@ mod tests {
     #[test]
     fn test_pauli_polynomial_constructor() {
         let size = 3;
-        let ham = vec![("IXYZ", 0.3), ("XXII", 0.7), ("YYII", 0.12)];
+        let ham = vec![
+            ("IXYZ", Angle::from_angle(0.3)),
+            ("XXII", Angle::from_angle(0.7)),
+            ("YYII", Angle::from_angle(0.12)),
+        ];
         let pp = PauliPolynomial::from_hamiltonian(ham);
 
         let pg1_ref = PauliString::from_text("IXY");
@@ -181,7 +227,7 @@ mod tests {
         let pg3_ref = PauliString::from_text("YII");
         let pg4_ref = PauliString::from_text("ZII");
 
-        let angles_ref = vec![0.3, 0.7, 0.12];
+        let angles_ref = Angle::from_angles(&[0.3, 0.7, 0.12]);
 
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref, pg4_ref],
@@ -201,7 +247,11 @@ mod tests {
     #[test]
     #[should_panic]
     fn test_pauli_polynomial_constructor_unequal_strings() {
-        let ham = vec![("IXYZ", 0.3), ("XXI", 0.7), ("YYII", 0.12)];
+        let ham = vec![
+            ("IXYZ", Angle::from_angle(0.3)),
+            ("XXI", Angle::from_angle(0.7)),
+            ("YYII", Angle::from_angle(0.12)),
+        ];
         let _ = PauliPolynomial::from_hamiltonian(ham);
     }
 
@@ -210,7 +260,11 @@ mod tests {
         let pg1_ref = PauliString::from_text("IXY");
         let pg2_ref = PauliString::from_text("ZYX");
         let pg3_ref = PauliString::from_text("YIX");
-        let angles_ref = vec![0.3, 0.7, 0.12];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.12),
+        ];
         PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -236,7 +290,11 @@ mod tests {
         let pg2_ref = PauliString::from_text("ZXY");
         // YIX
         let pg3_ref = PauliString::from_text("YIX");
-        let angles_ref = vec![0.3, -0.7, -0.12];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(6),
+            Angle::from_angle(-0.12),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -261,7 +319,11 @@ mod tests {
         let pg2_ref = PauliString::from_text("YZX");
         // YIX -> ZIX
         let pg3_ref = PauliString::from_text("ZIX");
-        let angles_ref = vec![-0.3, 0.7, 0.12];
+        let angles_ref = vec![
+            Angle::from_angle(-0.3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.12),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -286,7 +348,11 @@ mod tests {
         let pg2_ref = PauliString::from_text("ZXY");
         // YIX -> XI(-Y)
         let pg3_ref = PauliString::from_text("XIY");
-        let angles_ref = vec![0.3, 0.7, 0.12];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.12),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -311,7 +377,11 @@ mod tests {
         let pg2_ref = PauliString::from_text("YZX");
         // YIX -> (-Z)IX
         let pg3_ref = PauliString::from_text("ZIX");
-        let angles_ref = vec![-0.3, -0.7, 0.12];
+        let angles_ref = vec![
+            Angle::from_angle(-0.3),
+            Angle::from_pi4_rotation(6),
+            Angle::from_angle(0.12),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -336,7 +406,11 @@ mod tests {
         let pg2_ref = PauliString::from_text("XYZ");
         // YIX -
         let pg3_ref = PauliString::from_text("YIX");
-        let angles_ref = vec![0.3, -0.7, -0.12];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(6),
+            Angle::from_angle(-0.12),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -358,7 +432,12 @@ mod tests {
         let pg2_ref = PauliString::from_text("IXYZ");
         let pg3_ref = PauliString::from_text("YIXZ");
 
-        let angles_ref = vec![0.3, 0.7, 0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.15),
+        ];
 
         PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
@@ -385,7 +464,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, 1, 1, 1]
-        let angles_ref = vec![0.3, 0.7, 0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -411,7 +495,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, 1, 1, -1]
-        let angles_ref = vec![0.3, 0.7, 0.12, -0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(-0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -437,7 +526,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, 1, -1, 1]
-        let angles_ref = vec![0.3, 0.7, -0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(6),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -463,7 +557,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, 1, 1, 1]
-        let angles_ref = vec![0.3, 0.7, 0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -490,7 +589,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, 1, 1, 1]
-        let angles_ref = vec![0.3, 0.7, 0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -516,7 +620,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, 1, -1, 1]
-        let angles_ref = vec![0.3, 0.7, -0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(6),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -542,7 +651,12 @@ mod tests {
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
         // [1, -1, 1, 1]
-        let angles_ref = vec![0.3, -0.7, 0.12, 0.15];
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(5),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -567,8 +681,13 @@ mod tests {
         let pg2_ref = PauliString::from_text("IXYZ");
         // YIXZ
         let pg3_ref = PauliString::from_text("YIXZ");
-        // [1, 1, -1, 1]
-        let angles_ref = vec![0.3, 0.7, 0.12, 0.15];
+        // [1, 1, 1, 1]
+        let angles_ref = vec![
+            Angle::from_angle(0.3),
+            Angle::from_pi4_rotation(3),
+            Angle::from_pi4_rotation(2),
+            Angle::from_angle(0.15),
+        ];
         let pp_ref = PauliPolynomial {
             chains: vec![pg1_ref, pg2_ref, pg3_ref],
             angles: angles_ref,
@@ -576,4 +695,57 @@ mod tests {
         };
         assert_eq!(pp, pp_ref);
     }
+
+    #[test]
+    fn test_commutes_with_simple() {
+        let pp1s = vec![
+            vec![("I", Angle::from_angle(0.3))],
+            vec![("X", Angle::from_angle(0.5))],
+            vec![("Y", Angle::from_angle(0.7))],
+            vec![("Z", Angle::from_angle(0.9))],
+        ]
+        .into_iter()
+        .map(|ham| PauliPolynomial::from_hamiltonian(ham))
+        .collect::<Vec<_>>();
+
+        let pp2s = pp1s.clone();
+
+        for (i, (pp1, pp2)) in pp1s.iter().cartesian_product(pp2s.iter()).enumerate() {
+            if i <= 5 || i == 8 || i == 10 || i == 12 || i == 15 {
+                assert!(pp1.commutes_with(pp2));
+            } else {
+                assert!(!pp1.commutes_with(pp2));
+            }
+        }
+    }
+
+    #[test]
+    fn test_commutes_with() {
+        let pp1 = PauliPolynomial::from_hamiltonian(vec![
+            ("IYYX", Angle::from_angle(0.3)),
+            ("XXXI", Angle::from_angle(0.5)),
+        ]);
+
+        let pp2 = PauliPolynomial::from_hamiltonian(vec![
+            ("IYZZ", Angle::from_angle(0.7)),
+            ("ZZXI", Angle::from_angle(0.9)),
+        ]);
+
+        assert!(pp1.commutes_with(&pp2));
+    }
+
+    #[test]
+    fn test_not_commutes_with() {
+        let pp1 = PauliPolynomial::from_hamiltonian(vec![
+            ("IYYX", Angle::from_angle(0.3)),
+            ("XXXI", Angle::from_angle(0.5)),
+        ]);
+
+        let pp2 = PauliPolynomial::from_hamiltonian(vec![
+            ("IYZZ", Angle::from_angle(0.7)),
+            ("ZZXY", Angle::from_angle(0.9)),
+        ]);
+
+        assert!(!pp1.commutes_with(&pp2));
+    }
 }
diff --git a/synir/src/data_structures/pauli_polynomial/simplify.rs b/synir/src/data_structures/pauli_polynomial/simplify.rs
new file mode 100644
index 00000000..65e009e5
--- /dev/null
+++ b/synir/src/data_structures/pauli_polynomial/simplify.rs
@@ -0,0 +1,173 @@
+use crate::data_structures::PauliPolynomial;
+use itertools::Itertools;
+use std::collections::HashMap;
+
+/// Check for items in PauliPolynomial that has the same Pauli string
+/// Returns a vector of a tuple of Pauli string identifier and the locations of this repeated string
+pub(crate) fn check_repeats(pp: &PauliPolynomial) -> Vec<(usize, Vec<usize>)> {
+    let size = pp.size();
+    let length = pp.length();
+    let mut repeats = HashMap::<usize, Vec<usize>>::new();
+    for index in 0..length {
+        let mut num = 0;
+        for letter in 0..size {
+            num += (pp.chain(letter).x(index) as usize) << 2 * letter;
+            num += (pp.chain(letter).z(index) as usize) << 2 * letter + 1;
+        }
+        repeats
+            .entry(num)
+            .and_modify(|e: &mut Vec<usize>| e.push(index))
+            .or_insert(vec![index]);
+    }
+    repeats
+        .into_iter()
+        .filter(|(_, v)| v.len() > 1)
+        .sorted()
+        .collect_vec()
+}
+
+/// Takes output of check_repeats and merges items with the same Pauli string
+/// Assumes that all angles are of the same type
+pub(crate) fn _merge_repeats(
+    mut pp: PauliPolynomial,
+    merge_list: Vec<(usize, Vec<usize>)>,
+) -> PauliPolynomial {
+    let mut pp_merge_list = Vec::<usize>::new();
+    // merge all the angles first
+    for (_, angle_merge_list) in merge_list {
+        let merge_index = angle_merge_list[0];
+        let mut angle = pp.angle(merge_index);
+        for angle_index in angle_merge_list.iter().skip(1) {
+            angle += pp.angle(*angle_index);
+        }
+        pp.angles[merge_index] = angle;
+        pp_merge_list.extend_from_slice(&angle_merge_list[1..]);
+    }
+    // remove duplicate entries
+    pp_merge_list.sort_by(|a, b| b.cmp(a));
+    for remove_index in pp_merge_list {
+        pp.angles.remove(remove_index);
+        for chain_index in 0..pp.chains().len() {
+            pp.chains[chain_index].x.remove(remove_index);
+            pp.chains[chain_index].z.remove(remove_index);
+        }
+    }
+    pp
+}
+
+/// Merges all items in PauliPolynomial that share the same Pauli string
+pub fn merge_repeats(mut pp: PauliPolynomial) -> PauliPolynomial {
+    let repeats = check_repeats(&pp);
+    _merge_repeats(pp, repeats)
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::data_structures::Angle;
+    use crate::data_structures::PauliString;
+
+    use super::*;
+
+    #[test]
+    fn test_simple_check_repeats() {
+        // Combined reading from back -> 01 = 1
+        let pp = PauliPolynomial::from_hamiltonian(vec![
+            ("I", Angle::from_angle(1.0)),
+            ("X", Angle::from_angle(2.0)),
+            ("Z", Angle::from_angle(3.0)),
+            ("X", Angle::from_angle(4.0)),
+            ("X", Angle::from_angle(5.0)),
+        ]);
+        let repeats = check_repeats(&pp);
+        assert!(repeats.len() == 1);
+        assert_eq!(repeats, vec![(1, vec![1, 3, 4])]);
+    }
+
+    #[test]
+    fn test_check_repeats() {
+        // XIZY appears twice
+        // Z string -> 0011
+        // X string -> 1001
+        // Combined reading from back -> 11 10 00 01 = 225
+        let pp = PauliPolynomial::from_hamiltonian(vec![
+            ("XIZY", Angle::from_angle(1.0)),
+            ("XIZY", Angle::from_angle(2.0)),
+            ("YZZI", Angle::from_angle(3.0)),
+        ]);
+        let repeats = check_repeats(&pp);
+        assert!(repeats.len() == 1);
+        assert_eq!(repeats, vec![(225, vec![0, 1])]);
+    }
+
+    #[test]
+    fn test_multiple_repeats() {
+        // Combined reading from back -> 01 = 1
+        let pp = PauliPolynomial::from_hamiltonian(vec![
+            ("II", Angle::from_angle(1.0)),
+            ("IX", Angle::from_angle(2.0)),
+            ("ZZ", Angle::from_angle(3.0)),
+            ("IX", Angle::from_angle(4.0)),
+            ("ZZ", Angle::from_angle(5.0)),
+        ]);
+        let repeats = check_repeats(&pp);
+        assert!(repeats.len() == 2);
+        assert_eq!(repeats, vec![(4, vec![1, 3]), (10, vec![2, 4])]);
+    }
+
+    #[test]
+    fn test_simple_merge_repeats() {
+        // Combined reading from back -> 01 = 1
+        let pp = PauliPolynomial::from_hamiltonian(vec![
+            ("I", Angle::from_angle(1.0)),
+            ("X", Angle::from_angle(2.0)),
+            ("Z", Angle::from_angle(3.0)),
+            ("X", Angle::from_angle(4.0)),
+            ("X", Angle::from_angle(5.0)),
+        ]);
+        let repeats = check_repeats(&pp);
+
+        let pp = _merge_repeats(pp, repeats);
+
+        assert!(pp.chain(0).len() == 3);
+        assert_eq!(pp.chain(0), &PauliString::from_text("IXZ"));
+        assert_eq!(pp.angles, Angle::from_angles(&[1.0, 11.0, 3.0]));
+    }
+
+    #[test]
+    fn test_merge_repeats() {
+        let pp = PauliPolynomial::from_hamiltonian(vec![
+            ("XIZY", Angle::from_angle(1.0)),
+            ("XIZY", Angle::from_angle(2.0)),
+            ("YZZI", Angle::from_angle(3.0)),
+        ]);
+        let repeats = check_repeats(&pp);
+        let pp = _merge_repeats(pp, repeats);
+
+        assert!(pp.chain(0).len() == 2);
+        assert_eq!(pp.chain(0), &PauliString::from_text("XY"));
+        assert_eq!(pp.chain(1), &PauliString::from_text("IZ"));
+        assert_eq!(pp.chain(2), &PauliString::from_text("ZZ"));
+        assert_eq!(pp.chain(3), &PauliString::from_text("YI"));
+        assert_eq!(pp.angles, Angle::from_angles(&[3.0, 3.0]));
+    }
+
+    #[test]
+    fn test_multiple_merge_repeats() {
+        // Combined reading from back -> 01 = 1
+        let pp = PauliPolynomial::from_hamiltonian(vec![
+            ("II", Angle::from_angle(1.0)),
+            ("IX", Angle::from_angle(2.0)),
+            ("ZZ", Angle::from_angle(3.0)),
+            ("IX", Angle::from_angle(4.0)),
+            ("ZZ", Angle::from_angle(5.0)),
+        ]);
+
+        let repeats = check_repeats(&pp);
+        let pp = _merge_repeats(pp, repeats);
+
+        assert!(pp.chain(0).len() == 3);
+        assert_eq!(pp.chain(0), &PauliString::from_text("IIZ"));
+        assert_eq!(pp.chain(1), &PauliString::from_text("IXZ"));
+        assert_eq!(pp.angles, Angle::from_angles(&[1.0, 6.0, 8.0]));
+    }
+}
diff --git a/synir/src/data_structures/pauli_string.rs b/synir/src/data_structures/pauli_string.rs
index d66a3f5a..590b4df4 100644
--- a/synir/src/data_structures/pauli_string.rs
+++ b/synir/src/data_structures/pauli_string.rs
@@ -123,6 +123,21 @@ impl PauliString {
         mask &= &self.z;
         mask
     }
+
+    pub(crate) fn commutes_with(&self, other: &PauliString) -> bool {
+        assert_eq!(self.len(), other.len());
+        let length = self.len();
+        let mut commutes = true;
+        for index in 0..length {
+            let p1 = self.pauli(index);
+            let p2 = other.pauli(index);
+            if p1 == PauliLetter::I || p2 == PauliLetter::I || p1 == p2 {
+                continue;
+            }
+            commutes = !commutes;
+        }
+        commutes
+    }
 }
 
 pub(crate) fn cx(control: &mut PauliString, target: &mut PauliString) {
@@ -164,6 +179,7 @@ mod tests {
     use super::*;
     use bitvec::prelude::Lsb0;
     use bitvec::{bits, bitvec};
+    use itertools::Itertools;
 
     #[test]
     fn test_from_basis_int() {
@@ -300,4 +316,49 @@ mod tests {
         let pauli_string = PauliString::from_text("IXYZI");
         assert_eq!(pauli_string.to_string(), String::from("I X Y Z I"));
     }
+
+    #[test]
+    fn test_check_commute_simple() {
+        let ps1 = vec![
+            PauliString::from_text("I"),
+            PauliString::from_text("X"),
+            PauliString::from_text("Y"),
+            PauliString::from_text("Z"),
+        ];
+
+        let ps2 = vec![
+            PauliString::from_text("I"),
+            PauliString::from_text("X"),
+            PauliString::from_text("Y"),
+            PauliString::from_text("Z"),
+        ];
+
+        for (p1, p2) in ps1.iter().cartesian_product(ps2.iter()) {
+            if p1 == &PauliString::from_text("I") || p2 == &PauliString::from_text("I") || p1 == p2
+            {
+                assert!(p1.commutes_with(p2));
+                continue;
+            }
+            assert!(!p1.commutes_with(p2));
+        }
+    }
+
+    #[test]
+    #[should_panic]
+    fn test_bad_commute() {
+        let ps1 = PauliString::from_text("IXXYZ");
+        let ps2 = PauliString::from_text("IYZX");
+        ps1.commutes_with(&ps2);
+    }
+
+    #[test]
+    fn test_check_commute() {
+        let ps1 = PauliString::from_text("IXYZ");
+        let ps2 = PauliString::from_text("IYZX");
+        assert!(!ps1.commutes_with(&ps2));
+
+        let ps3 = PauliString::from_text("IXYZ");
+        let ps4 = PauliString::from_text("IZXI");
+        assert!(ps3.commutes_with(&ps4));
+    }
 }
diff --git a/synir/src/ir/pauli_exponential.rs b/synir/src/ir/pauli_exponential.rs
index 79be5518..9faae0c9 100644
--- a/synir/src/ir/pauli_exponential.rs
+++ b/synir/src/ir/pauli_exponential.rs
@@ -1,7 +1,7 @@
-use std::collections::VecDeque;
-
-use crate::architecture::connectivity::Connectivity;
-use crate::data_structures::{CliffordTableau, HasAdjoint, PauliPolynomial};
+use crate::{
+    architecture::connectivity::Connectivity,
+    data_structures::{HasAdjoint, PauliExponential},
+};
 
 use crate::ir::{CliffordGates, Gates, Synthesizer};
 
@@ -12,24 +12,6 @@ use crate::ir::{
     pauli_polynomial::{naive::NaivePauliPolynomialSynthesizer, PauliPolynomialSynthStrategy},
 };
 
-#[derive(Default)]
-pub struct PauliExponential {
-    pauli_polynomials: VecDeque<PauliPolynomial>,
-    clifford_tableau: CliffordTableau,
-}
-
-impl PauliExponential {
-    pub fn new(
-        pauli_polynomials: VecDeque<PauliPolynomial>,
-        clifford_tableau: CliffordTableau,
-    ) -> Self {
-        PauliExponential {
-            pauli_polynomials,
-            clifford_tableau,
-        }
-    }
-}
-
 #[derive(Default)]
 pub struct PauliExponentialSynthesizer {
     pauli_strategy: PauliPolynomialSynthStrategy,
diff --git a/synir/src/ir/pauli_polynomial/helper.rs b/synir/src/ir/pauli_polynomial/helper.rs
index 127f5ad1..536a747b 100644
--- a/synir/src/ir/pauli_polynomial/helper.rs
+++ b/synir/src/ir/pauli_polynomial/helper.rs
@@ -108,7 +108,7 @@ pub(super) fn push_down_pauli_polynomial_update<G>(
             }
         }
         let last_qubit = *affected_qubits.last().unwrap();
-        repr.rz(last_qubit, pauli_polynomial.angle(col));
+        repr.rz(last_qubit, pauli_polynomial.angle(col).to_radians());
         mask.replace(col, false);
     }
 }
diff --git a/synir/tests/ct_synthesis/ct_compose.rs b/synir/tests/ct_synthesis/ct_compose.rs
new file mode 100644
index 00000000..e7fbae63
--- /dev/null
+++ b/synir/tests/ct_synthesis/ct_compose.rs
@@ -0,0 +1,184 @@
+extern crate rand;
+
+use crate::common::mock_circuit::{check_mock_equals_clifford_tableau, MockCircuit, MockCommand};
+use itertools::Itertools;
+use synir::data_structures::{Angle, CliffordTableau, PropagateClifford};
+use synir::ir::clifford_tableau::CallbackCliffordSynthesizer;
+use synir::ir::Synthesizer;
+
+fn run_synthesizer(clifford_tableau: &CliffordTableau) -> (MockCircuit, CliffordTableau) {
+    let mut mock = MockCircuit::new();
+    // let mut rng = rand::rng(); //TODO make this from seed
+    let custom_columns = (0..clifford_tableau.size()).collect_vec();
+    let custom_rows = (0..clifford_tableau.size()).collect_vec();
+    // custom_rows.shuffle(&mut rng);
+
+    let mut synthesizer = CallbackCliffordSynthesizer::custom_pivot(custom_columns, custom_rows);
+    let new_ct = synthesizer.synthesize(clifford_tableau.clone(), &mut mock);
+    (mock, new_ct)
+}
+
+macro_rules! test_clifford {
+    ($fun:ident, $expected:expr) => {
+        paste::item! {
+                #[test]
+                fn [< synthesize_ $fun>]() {
+                    let clifford_tableau = $fun();
+                    let (mock, new_ct) = run_synthesizer(&clifford_tableau);
+                    if $expected.is_some() {
+                        assert_eq!(mock.commands(), $expected.unwrap());
+                    }
+                    check_mock_equals_clifford_tableau(&clifford_tableau, &mock, new_ct.get_permutation());
+                }
+            }
+    };
+}
+
+fn compose_x_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(2);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("XI"),
+        Angle::from_pi4_rotation(2),
+    ));
+    ct
+}
+
+fn compose_z_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(2);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("ZI"),
+        Angle::from_pi4_rotation(4),
+    ));
+    ct
+}
+
+fn compose_y_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(2);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("YI"),
+        Angle::from_pi4_rotation(6),
+    ));
+    ct
+}
+
+test_clifford!(compose_x_gadget, Some(&vec![MockCommand::V(0),]));
+test_clifford!(compose_z_gadget, Some(&vec![MockCommand::Z(0),]));
+
+test_clifford!(
+    compose_y_gadget,
+    Some(&vec![MockCommand::H(0), MockCommand::Z(0)])
+);
+
+fn compose_xx_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(2);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("XX"),
+        Angle::from_pi4_rotation(2),
+    ));
+    ct
+}
+
+fn compose_zz_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(2);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("ZZ"),
+        Angle::from_pi4_rotation(2),
+    ));
+    ct
+}
+
+fn compose_yy_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(2);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("YY"),
+        Angle::from_pi4_rotation(2),
+    ));
+    ct
+}
+
+test_clifford!(
+    compose_xx_gadget,
+    Some(&vec![
+        MockCommand::V(0),
+        MockCommand::H(1),
+        MockCommand::CX(1, 0),
+        MockCommand::H(1),
+        MockCommand::V(1)
+    ])
+);
+
+test_clifford!(
+    compose_zz_gadget,
+    Some(&vec![
+        MockCommand::S(0),
+        MockCommand::H(1),
+        MockCommand::CX(0, 1),
+        MockCommand::S(1),
+        MockCommand::V(1),
+        MockCommand::Z(1)
+    ])
+);
+
+test_clifford!(
+    compose_yy_gadget,
+    Some(&vec![
+        MockCommand::H(0),
+        MockCommand::S(1),
+        MockCommand::CX(0, 1),
+        MockCommand::H(1),
+        MockCommand::CX(1, 0),
+        MockCommand::S(1),
+        MockCommand::V(1),
+        MockCommand::X(1),
+        MockCommand::Z(0),
+        MockCommand::Z(1)
+    ])
+);
+
+fn compose_complex_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(3);
+    let _ = ct.compose_gadget((
+        synir::data_structures::PauliString::from_text("XYZ"),
+        Angle::from_pi4_rotation(2),
+    ));
+    ct
+}
+
+fn manual_compose_complex_gadget() -> CliffordTableau {
+    let mut ct = CliffordTableau::new(3);
+    ct.h(0);
+    ct.v(1);
+    ct.cx(0, 1);
+    ct.cx(1, 2);
+    ct.s(2);
+    ct.cx(1, 2);
+    ct.cx(0, 1);
+    ct.h(0);
+    ct.v_dgr(1);
+    ct
+}
+
+test_clifford!(
+    compose_complex_gadget,
+    Some(&vec![
+        MockCommand::V(0),
+        MockCommand::V(1),
+        MockCommand::CX(1, 0),
+        MockCommand::CX(2, 0),
+        MockCommand::S(1),
+        MockCommand::H(2),
+        MockCommand::CX(1, 2),
+        MockCommand::V(1),
+        MockCommand::S(2),
+        MockCommand::V(2),
+        MockCommand::X(1),
+        MockCommand::Z(2)
+    ])
+);
+
+#[test]
+fn test_correctness() {
+    let ct1 = compose_complex_gadget();
+    let ct2 = manual_compose_complex_gadget();
+    assert_eq!(ct1, ct2);
+}
diff --git a/synir/tests/ct_synthesis/mod.rs b/synir/tests/ct_synthesis/mod.rs
index 2fff57d6..3d2bebc6 100644
--- a/synir/tests/ct_synthesis/mod.rs
+++ b/synir/tests/ct_synthesis/mod.rs
@@ -1,3 +1,4 @@
+pub mod ct_compose;
 pub mod custom_callback;
 pub mod naive;
 pub mod naive_adjoint;
diff --git a/synir/tests/pauli_exponential.rs b/synir/tests/pauli_exponential.rs
index 324aa30d..cc9c89c6 100644
--- a/synir/tests/pauli_exponential.rs
+++ b/synir/tests/pauli_exponential.rs
@@ -3,15 +3,16 @@ mod common;
 use std::collections::VecDeque;
 
 use common::mock_circuit::{parse_clifford_commands, MockCircuit, MockCommand};
-use synir::data_structures::{CliffordTableau, HasAdjoint, PauliPolynomial};
+use synir::data_structures::{
+    Angle, CliffordTableau, HasAdjoint, PauliExponential, PauliPolynomial,
+};
 use synir::ir::clifford_tableau::{CliffordTableauSynthStrategy, NaiveCliffordSynthesizer};
-use synir::ir::pauli_exponential::PauliExponential;
 use synir::ir::pauli_exponential::PauliExponentialSynthesizer;
 use synir::ir::pauli_polynomial::PauliPolynomialSynthStrategy;
 use synir::ir::Synthesizer;
 
 fn setup_simple_pe() -> PauliExponential {
-    let ham = vec![("IZZZ", 0.3)];
+    let ham = vec![("IZZZ", Angle::from_angle(0.3))];
 
     let pauli_polynomial = PauliPolynomial::from_hamiltonian(ham);
     let clifford_tableau = CliffordTableau::new(4);
@@ -19,7 +20,11 @@ fn setup_simple_pe() -> PauliExponential {
 }
 
 fn setup_complex_pe() -> PauliExponential {
-    let ham = vec![("IXYZ", 0.3), ("XXII", 0.7), ("YYII", 0.12)];
+    let ham = vec![
+        ("IXYZ", Angle::from_angle(0.3)),
+        ("XXII", Angle::from_angle(0.7)),
+        ("YYII", Angle::from_angle(0.12)),
+    ];
 
     let pauli_polynomial = PauliPolynomial::from_hamiltonian(ham);
     let clifford_tableau = CliffordTableau::new(4);
diff --git a/synir/tests/pauli_polynomial.rs b/synir/tests/pauli_polynomial.rs
index f4560e9b..7d1d8103 100644
--- a/synir/tests/pauli_polynomial.rs
+++ b/synir/tests/pauli_polynomial.rs
@@ -3,13 +3,14 @@ mod common;
 use std::collections::VecDeque;
 
 use common::mock_circuit::{parse_clifford_commands, MockCircuit, MockCommand};
+use synir::data_structures::Angle;
 use synir::data_structures::{CliffordTableau, PauliPolynomial};
 use synir::ir::pauli_polynomial::NaivePauliPolynomialSynthesizer;
 use synir::ir::Synthesizer;
 
 fn setup_complex_pp() -> VecDeque<PauliPolynomial> {
-    let ham_1 = vec![("IZZZ", 0.3)];
-    let ham_2 = vec![("XXII", 0.7)];
+    let ham_1 = vec![("IZZZ", Angle::from_angle(0.3))];
+    let ham_2 = vec![("XXII", Angle::from_angle(0.7))];
 
     let pp_1 = PauliPolynomial::from_hamiltonian(ham_1);
     let pp_2 = PauliPolynomial::from_hamiltonian(ham_2);
@@ -17,7 +18,7 @@ fn setup_complex_pp() -> VecDeque<PauliPolynomial> {
 }
 
 fn setup_simple_pp() -> VecDeque<PauliPolynomial> {
-    let ham = vec![("IXYZ", 0.3)];
+    let ham = vec![("IXYZ", Angle::from_angle(0.3))];
 
     let pauli_polynomial = PauliPolynomial::from_hamiltonian(ham);
 
diff --git a/synpy/integration_tests/test_qiskit.py b/synpy/integration_tests/test_qiskit.py
index 95dc3f37..488bf8bf 100644
--- a/synpy/integration_tests/test_qiskit.py
+++ b/synpy/integration_tests/test_qiskit.py
@@ -1,7 +1,8 @@
-from qiskit.quantum_info import Clifford
+from qiskit.quantum_info import Clifford, Operator
 from qiskit import QuantumCircuit
 
-from synpy.qiskit.plugin import SynPyCliffordPlugin
+
+from synpy.qiskit.plugin import SynPyCliffordPlugin, qiskit_to_synir
 from synpy.synpy_rust import QiskitSynIR
 
 
@@ -37,10 +38,28 @@ def test_qiskit_bell() -> None:
     qc = QuantumCircuit(2)
     qc.h(0)
     qc.cx(0, 1)
-
     cliff = Clifford(qc)
 
     plugin = SynPyCliffordPlugin()
     circ = plugin.run(cliff, None, None, [])
 
     assert circ == qc
+
+
+def test_qiskit_loop() -> None:
+    circuit = QuantumCircuit(3)
+    circuit.h(0)
+    circuit.cx(0, 1)
+    circuit.rz(1.5, 1)
+    sample_circuit = circuit.copy()
+
+    pe_wrap = qiskit_to_synir(circuit)
+
+    synir_result = QiskitSynIR(circuit.copy_empty_like())
+    pe_wrap.synthesize_to_qiskit(synir_result)
+    circuit = synir_result.get_circuit()
+
+    op1 = Operator.from_circuit(circuit)
+    op2 = Operator.from_circuit(sample_circuit)
+
+    assert op1.equiv(op2)
diff --git a/synpy/python/synpy/qiskit/__init__.py b/synpy/python/synpy/qiskit/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/synpy/python/synpy/qiskit/plugin.py b/synpy/python/synpy/qiskit/plugin.py
index d41ac09a..9bd6855a 100644
--- a/synpy/python/synpy/qiskit/plugin.py
+++ b/synpy/python/synpy/qiskit/plugin.py
@@ -1,9 +1,9 @@
 from qiskit.transpiler import CouplingMap, Target
 from qiskit.transpiler.passes.synthesis.plugin import HighLevelSynthesisPlugin
 from qiskit.quantum_info import Clifford
-from qiskit import QuantumCircuit
+from qiskit import QuantumCircuit, transpile
 
-from synpy.synpy_rust import PyCliffordTableau
+from synpy.synpy_rust import PyCliffordTableau, PyPauliExponential
 from synpy.utils import pycommand_to_qasm
 
 
@@ -26,3 +26,19 @@ def run(self, clifford: Clifford, coupling_map: CouplingMap, target: Target, qub
         commands = synpy_tableau.synthesize()
         qasm = pycommand_to_qasm(n, commands)
         return QuantumCircuit.from_qasm_str(qasm)
+
+
+def qiskit_to_synir(circuit: QuantumCircuit) -> PyPauliExponential:
+    new_circuit = transpile(circuit, basis_gates=["cx", "h", "rz"])
+    pe = PyPauliExponential(new_circuit.num_qubits)
+
+    for gate in reversed(new_circuit.data):
+        if gate.name == "cx":
+            pe.add_cx(gate.qubits[0]._index, gate.qubits[1]._index)
+        elif gate.name == "h":
+            pe.add_h(gate.qubits[0]._index)
+        elif gate.name == "rz":
+            pe.add_rz(gate.qubits[0]._index, gate.params[0])
+        else:
+            raise Exception("Gate is not supported")
+    return pe
diff --git a/synpy/src/synthesis.rs b/synpy/src/synthesis.rs
index 6e905add..8dc4898e 100644
--- a/synpy/src/synthesis.rs
+++ b/synpy/src/synthesis.rs
@@ -1,17 +1,18 @@
 use pyo3::exceptions::PyException;
 use pyo3::{pyclass, pymethods, PyErr};
 use std::ops::Deref;
+use synir::data_structures::Angle;
 
 use pyo3::{pyfunction, PyRef, PyResult};
-use synir::data_structures::PropagateClifford;
-use synir::data_structures::{CliffordTableau, PauliPolynomial};
-use synir::ir::clifford_tableau::CliffordTableauSynthStrategy;
-use synir::ir::pauli_exponential::{PauliExponential, PauliExponentialSynthesizer};
-use synir::ir::pauli_polynomial::PauliPolynomialSynthStrategy;
-use synir::ir::CliffordGates;
-use synir::ir::Gates;
-use synir::ir::Synthesizer;
-use synir::IndexType;
+use synir::{
+    data_structures::{CliffordTableau, PauliExponential, PauliPolynomial, PropagateClifford},
+    ir::{
+        clifford_tableau::CliffordTableauSynthStrategy,
+        pauli_exponential::PauliExponentialSynthesizer,
+        pauli_polynomial::PauliPolynomialSynthStrategy, CliffordGates, Gates, Synthesizer,
+    },
+    IndexType,
+};
 
 use crate::validation::validate;
 
@@ -210,8 +211,14 @@ pub fn synthesize_pauli_exponential(
 ) -> PyResult<Vec<PyCommand>> {
     let converted_hamiltonian = hamiltonian
         .iter()
-        .map(|inner| inner.iter().map(PyPauliString::as_tuple).collect())
-        .map(|inner: Vec<(&str, f64)>| PauliPolynomial::from_hamiltonian(inner))
+        .map(|inner| {
+            inner
+                .iter()
+                .map(PyPauliString::as_tuple)
+                .map(|(pauli_str, phase)| (pauli_str, Angle::from_angle(phase)))
+                .collect::<Vec<(&str, Angle)>>()
+        })
+        .map(|inner: Vec<(&str, Angle)>| PauliPolynomial::from_hamiltonian(inner))
         .collect();
     let clifford_gates: Vec<PyCommand> = clifford_gates.iter().map(|cmd| *(cmd.deref())).collect();
 
diff --git a/synpy/src/wrapper.rs b/synpy/src/wrapper.rs
index ab24a0e6..6f991a14 100644
--- a/synpy/src/wrapper.rs
+++ b/synpy/src/wrapper.rs
@@ -7,12 +7,13 @@ use std::collections::VecDeque;
 
 use pyo3::prelude::*;
 use synir::{
-    data_structures::CliffordTableau,
+    data_structures::{
+        Angle, CliffordTableau, PauliExponential, PauliPolynomial, PropagateClifford,
+    },
     ir::{
         clifford_tableau::CliffordTableauSynthStrategy,
-        pauli_exponential::{PauliExponential, PauliExponentialSynthesizer},
-        pauli_polynomial::PauliPolynomialSynthStrategy,
-        CliffordGates, Gates, Synthesizer,
+        pauli_exponential::PauliExponentialSynthesizer,
+        pauli_polynomial::PauliPolynomialSynthStrategy, CliffordGates, Gates, Synthesizer,
     },
 };
 
@@ -47,7 +48,6 @@ impl PyPauliExponential {
             _ => panic!("Unknown Pauli polynomial synthesis strategy: {}", strategy),
         }
     }
-
     pub fn set_tableau_strategy(&mut self, strategy: String) {
         match strategy.as_str() {
             "Naive" => self.tableau_strategy = CliffordTableauSynthStrategy::Naive,
@@ -55,6 +55,62 @@ impl PyPauliExponential {
             _ => panic!("Unknown Clifford tableau synthesis strategy: {}", strategy),
         }
     }
+
+    pub fn add_h(&mut self, target: usize) {
+        self.pe.h(target);
+    }
+
+    pub fn add_s(&mut self, target: usize) {
+        self.pe.s(target);
+    }
+
+    pub fn add_s_dgr(&mut self, target: usize) {
+        self.pe.s_dgr(target);
+    }
+
+    pub fn add_x(&mut self, target: usize) {
+        self.pe.x(target);
+    }
+
+    pub fn add_y(&mut self, target: usize) {
+        self.pe.y(target);
+    }
+
+    pub fn add_z(&mut self, target: usize) {
+        self.pe.z(target);
+    }
+
+    pub fn add_cx(&mut self, control: usize, target: usize) {
+        self.pe.cx(control, target);
+    }
+
+    pub fn add_rz(&mut self, target: usize, angle: f64) {
+        let size = self.pe.size();
+        let ppvec = self.pe.chains();
+        if let Some(pp) = ppvec.front_mut() {
+            if pp.chain(target).x_weight() != 0 {
+                pp.extend_z(target, angle);
+            }
+            return;
+        }
+        let newpp = PauliPolynomial::from_hamiltonian(vec![(
+            &to_pauli_component(size, &target, 'Z'),
+            Angle::Arbitrary(angle),
+        )]);
+        ppvec.push_front(newpp);
+    }
+}
+
+fn to_pauli_component(size: usize, target: &usize, pauli: char) -> String {
+    let mut term = String::new();
+    for i in 0..size {
+        if i == *target {
+            term.push(pauli);
+        } else {
+            term.push('I');
+        }
+    }
+    term
 }
 
 pub fn synthesize<G>(pe: &mut PyPauliExponential, circuit: &mut G)
@@ -66,5 +122,5 @@ where
         pe.tableau_strategy.clone(),
     );
     let pe = std::mem::take(&mut pe.pe);
-    synth.synthesize(pe, circuit)
+    synth.synthesize(pe, circuit);
 }
diff --git a/synpy/src/wrapper/qiskit.rs b/synpy/src/wrapper/qiskit.rs
index 7e7f7606..da4a139c 100644
--- a/synpy/src/wrapper/qiskit.rs
+++ b/synpy/src/wrapper/qiskit.rs
@@ -1,7 +1,7 @@
 extern crate pyo3;
 extern crate pyo3_ffi;
 
-use pyo3::prelude::*;
+use pyo3::{intern, prelude::*};
 use synir::ir::{CliffordGates, Gates};
 
 #[pyclass]
@@ -18,6 +18,10 @@ impl QiskitSynIR {
         }
     }
 
+    fn get_circuit(&self, py: Python) -> Py<PyAny> {
+        self.circuit.clone_ref(py)
+    }
+
     pub fn s(&mut self, target: synir::IndexType) {
         Python::attach(|py| -> PyResult<()> {
             self.circuit.call_method1(py, "s", (target,))?;