diff --git a/crates/ekore/refs.bib b/crates/ekore/refs.bib
index a8a049ac1..0035b7a39 100644
--- a/crates/ekore/refs.bib
+++ b/crates/ekore/refs.bib
@@ -94,3 +94,16 @@ @article{Buza1996wv
     pages = "301--320",
     year = "1998"
 }
+@article{Bierenbaum2009zt,
+    author = "Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian",
+    title = "{The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production}",
+    eprint = "0901.0669",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY-08-187, SFB-CPP-08-107, IFIC-08-68",
+    doi = "10.1016/j.physletb.2009.01.057",
+    journal = "Phys. Lett. B",
+    volume = "672",
+    pages = "401--406",
+    year = "2009"
+}
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
index 5b5fc6ee6..e35a87472 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as1.rs
@@ -54,10 +54,11 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
 
 #[cfg(test)]
 mod tests {
-    use crate::cmplx;
-    use crate::{anomalous_dimensions::unpolarized::spacelike::as1::*, harmonics::cache::Cache};
-    use float_cmp::assert_approx_eq;
+    use super::*;
+    use crate::harmonics::cache::Cache;
+    use crate::{assert_approx_eq_cmplx, cmplx};
     use num::complex::Complex;
+    use num::Zero;
     const NF: u8 = 5;
 
     #[test]
@@ -65,8 +66,7 @@ mod tests {
         const N: Complex<f64> = cmplx![1., 0.];
         let mut c = Cache::new(N);
         let me = gamma_ns(&mut c, NF);
-        assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
-        assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
+        assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
     }
 
     #[test]
@@ -74,8 +74,7 @@ mod tests {
         const N: Complex<f64> = cmplx![2., 0.];
         let mut c = Cache::new(N);
         let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF);
-        assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
-        assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
+        assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
     }
 
     #[test]
@@ -83,8 +82,7 @@ mod tests {
         const N: Complex<f64> = cmplx![2., 0.];
         let mut c = Cache::new(N);
         let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF);
-        assert_approx_eq!(f64, me.re, 0., epsilon = 1e-12);
-        assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
+        assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12);
     }
 
     #[test]
@@ -92,8 +90,7 @@ mod tests {
         const N: Complex<f64> = cmplx![1., 0.];
         let mut c = Cache::new(N);
         let me = gamma_qg(&mut c, NF);
-        assert_approx_eq!(f64, me.re, -20. / 3.0, ulps = 32);
-        assert_approx_eq!(f64, me.im, 0., epsilon = 1e-12);
+        assert_approx_eq_cmplx!(f64, me, cmplx!(-20. / 3., 0.), ulps = 32, epsilon = 1e-12);
     }
 
     #[test]
@@ -101,8 +98,7 @@ mod tests {
         const N: Complex<f64> = cmplx![0., 1.];
         let mut c = Cache::new(N);
         let me = gamma_gq(&mut c, NF);
-        assert_approx_eq!(f64, me.re, 4. / 3.0, ulps = 32);
-        assert_approx_eq!(f64, me.im, -4. / 3.0, ulps = 32);
+        assert_approx_eq_cmplx!(f64, me, cmplx!(4. / 3.0, -4. / 3.0), ulps = 32);
     }
 
     #[test]
@@ -110,7 +106,12 @@ mod tests {
         const N: Complex<f64> = cmplx![0., 1.];
         let mut c = Cache::new(N);
         let me = gamma_gg(&mut c, NF);
-        assert_approx_eq!(f64, me.re, 5.195725159621, ulps = 32, epsilon = 1e-11);
-        assert_approx_eq!(f64, me.im, 10.52008856962, ulps = 32, epsilon = 1e-11);
+        assert_approx_eq_cmplx!(
+            f64,
+            me,
+            cmplx!(5.195725159621, 10.52008856962),
+            ulps = 32,
+            epsilon = 1e-11
+        );
     }
 }
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
index 857034ab8..f1d0266d5 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as2.rs
@@ -1,6 +1,6 @@
 //! |NLO| |QCD|.
 
-use ::num::complex::Complex;
+use num::complex::Complex;
 use std::f64::consts::LN_2;
 
 use crate::constants::{CA, CF, TR, ZETA2, ZETA3};
@@ -216,11 +216,12 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
 
 #[cfg(test)]
 mod tests {
-    use crate::cmplx;
-    use crate::{anomalous_dimensions::unpolarized::spacelike::as2::*, harmonics::cache::Cache};
-    use float_cmp::assert_approx_eq;
+    use super::*;
+    use crate::harmonics::cache::Cache;
+    use crate::{assert_approx_eq_cmplx, cmplx};
     use num::complex::Complex;
     use num::traits::Pow;
+    use num::Zero;
     use std::f64::consts::PI;
 
     const NF: u8 = 5;
@@ -228,29 +229,42 @@ mod tests {
     #[test]
     fn physical_constraints() {
         // number conservation
-        let mut c = Cache::new(cmplx![1., 0.]);
-        assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, 0.0, epsilon = 2e-6);
+        let mut c = Cache::new(cmplx!(1., 0.));
+        assert_approx_eq_cmplx!(f64, gamma_nsm(&mut c, NF), Complex::zero(), epsilon = 2e-6);
 
         // momentum conservation
-        let mut c = Cache::new(cmplx![2., 0.]);
+        let mut c = Cache::new(cmplx!(2., 0.));
         let gS1 = gamma_singlet(&mut c, NF);
 
         // gluon momentum conservation
-        assert_approx_eq!(f64, (gS1[0][1] + gS1[1][1]).re, 0.0, epsilon = 4e-5);
+        assert_approx_eq_cmplx!(
+            f64,
+            (gS1[0][1] + gS1[1][1]),
+            Complex::zero(),
+            epsilon = 4e-5
+        );
         // quark momentum conservation
-        assert_approx_eq!(f64, (gS1[0][0] + gS1[1][0]).re, 0.0, epsilon = 2e-6);
+        assert_approx_eq_cmplx!(
+            f64,
+            (gS1[0][0] + gS1[1][0]),
+            Complex::zero(),
+            epsilon = 2e-6
+        );
     }
 
     #[test]
     fn N2() {
         // reference values are obtained from MMa
-        let mut c = Cache::new(cmplx![2., 0.]);
+        let mut c = Cache::new(cmplx!(2., 0.));
 
         // ns+
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gamma_nsp(&mut c, NF).re,
-            (-112.0 * CF + 376.0 * CA - 64.0 * (NF as f64)) * CF / 27.0,
+            gamma_nsp(&mut c, NF),
+            cmplx!(
+                (-112.0 * CF + 376.0 * CA - 64.0 * (NF as f64)) * CF / 27.0,
+                0.
+            ),
             epsilon = 2e-6
         );
 
@@ -259,80 +273,102 @@ mod tests {
             + (373.0 / 9.0 - 34.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA
             - 64.0 * (NF as f64) / 27.0)
             * CF;
-        assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, check, epsilon = 2e-6);
+        assert_approx_eq_cmplx!(
+            f64,
+            gamma_nsm(&mut c, NF),
+            cmplx!(check, 0.),
+            epsilon = 2e-6
+        );
 
         // singlet sector
         let gS1 = gamma_singlet(&mut c, NF);
         // ps
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gamma_ps(&mut c, NF).re,
-            -40.0 * CF * (NF as f64) / 27.0
+            gamma_ps(&mut c, NF),
+            cmplx!(-40.0 * CF * (NF as f64) / 27.0, 0.)
         );
         // qg
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gS1[0][1].re,
-            (-74.0 * CF - 35.0 * CA) * (NF as f64) / 27.0
+            gS1[0][1],
+            cmplx!((-74.0 * CF - 35.0 * CA) * (NF as f64) / 27.0, 0.)
         );
         // gq
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gS1[1][0].re,
-            (112.0 * CF - 376.0 * CA + 104.0 * (NF as f64)) * CF / 27.0,
+            gS1[1][0],
+            cmplx!(
+                (112.0 * CF - 376.0 * CA + 104.0 * (NF as f64)) * CF / 27.0,
+                0.
+            ),
             epsilon = 1e-13
         );
     }
 
     #[test]
     fn N3() {
-        let mut c = Cache::new(cmplx![3., 0.]);
+        let mut c = Cache::new(cmplx!(3., 0.));
         // ns+
         let check = ((-34487.0 / 432.0 + 86.0 * PI.pow(2) / 9.0 - 16.0 * ZETA3) * CF
             + (459.0 / 8.0 - 43.0 * PI.pow(2) / 9.0 + 8.0 * ZETA3) * CA
             - 415.0 * (NF as f64) / 108.0)
             * CF;
-        assert_approx_eq!(f64, gamma_nsp(&mut c, NF).re, check, epsilon = 2e-6);
+        assert_approx_eq_cmplx!(
+            f64,
+            gamma_nsp(&mut c, NF),
+            cmplx!(check, 0.),
+            epsilon = 2e-6
+        );
 
         // singlet sector
         let gS1 = gamma_singlet(&mut c, NF);
         // ps
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gamma_ps(&mut c, NF).re,
-            -1391.0 * CF * (NF as f64) / 5400.0
+            gamma_ps(&mut c, NF),
+            cmplx!(-1391.0 * CF * (NF as f64) / 5400.0, 0.)
         );
         // gq
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gS1[1][0].re,
-            (973.0 / 432.0 * CF
-                + (2801.0 / 5400.0 - 7.0 * PI.pow(2) / 9.0) * CA
-                + 61.0 / 54.0 * (NF as f64))
-                * CF
+            gS1[1][0],
+            cmplx!(
+                (973.0 / 432.0 * CF
+                    + (2801.0 / 5400.0 - 7.0 * PI.pow(2) / 9.0) * CA
+                    + 61.0 / 54.0 * (NF as f64))
+                    * CF,
+                0.
+            )
         );
         // gg
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gS1[1][1].re,
-            (-79909.0 / 3375.0 + 194.0 * PI.pow(2) / 45.0 - 8.0 * ZETA3) * CA.pow(2)
-                - 967.0 / 270.0 * CA * (NF as f64)
-                + 541.0 / 216.0 * CF * (NF as f64),
+            gS1[1][1],
+            cmplx!(
+                (-79909.0 / 3375.0 + 194.0 * PI.pow(2) / 45.0 - 8.0 * ZETA3) * CA.pow(2)
+                    - 967.0 / 270.0 * CA * (NF as f64)
+                    + 541.0 / 216.0 * CF * (NF as f64),
+                0.
+            ),
             epsilon = 3e-5
         );
     }
 
     #[test]
     fn N4() {
-        let mut c = Cache::new(cmplx![4., 0.]);
+        let mut c = Cache::new(cmplx!(4., 0.));
         // singlet sector
         let gS1 = gamma_singlet(&mut c, NF);
         // qg
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            gS1[0][1].re,
-            (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * (NF as f64),
+            gS1[0][1],
+            cmplx!(
+                (-56317.0 / 18000.0 * CF + 16387.0 / 9000.0 * CA) * (NF as f64),
+                0.
+            ),
             epsilon = 1e-14
-        )
+        );
     }
 }
diff --git a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
index d491db0e9..223d9406a 100644
--- a/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
+++ b/crates/ekore/src/anomalous_dimensions/unpolarized/spacelike/as3.rs
@@ -425,9 +425,9 @@ pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex<f64>; 2]; 2] {
 
 #[cfg(test)]
 mod tests {
-    use crate::cmplx;
-    use crate::{anomalous_dimensions::unpolarized::spacelike::as3::*, harmonics::cache::Cache};
-    use float_cmp::assert_approx_eq;
+    use super::*;
+    use crate::harmonics::cache::Cache;
+    use crate::{assert_approx_eq_cmplx, cmplx};
     use num::complex::Complex;
 
     const NF: u8 = 5;
@@ -436,20 +436,45 @@ mod tests {
     fn physical_constraints() {
         // number conservation
         let mut c = Cache::new(cmplx![1., 0.]);
-        assert_approx_eq!(f64, gamma_nsv(&mut c, NF).re, -0.000960586, epsilon = 3e-7);
-        assert_approx_eq!(f64, gamma_nsm(&mut c, NF).re, 0.000594225, epsilon = 6e-7);
+        assert_approx_eq_cmplx!(
+            f64,
+            gamma_nsv(&mut c, NF),
+            cmplx!(-0.000960586, 0.),
+            epsilon = 3e-7
+        );
+        assert_approx_eq_cmplx!(
+            f64,
+            gamma_nsm(&mut c, NF),
+            cmplx!(0.000594225, 0.),
+            epsilon = 6e-7
+        );
 
         let mut c = Cache::new(cmplx![2., 0.]);
         let gS2 = gamma_singlet(&mut c, NF);
         // gluon momentum conservation
-        assert_approx_eq!(f64, (gS2[0][1] + gS2[1][1]).re, -0.00388726, epsilon = 2e-6);
+        assert_approx_eq_cmplx!(
+            f64,
+            (gS2[0][1] + gS2[1][1]),
+            cmplx!(-0.00388726, 0.),
+            epsilon = 2e-6
+        );
         // quark momentum conservation
-        assert_approx_eq!(f64, (gS2[0][0] + gS2[1][0]).re, 0.00169375, epsilon = 2e-6);
+        assert_approx_eq_cmplx!(
+            f64,
+            (gS2[0][0] + gS2[1][0]),
+            cmplx!(0.00169375, 0.),
+            epsilon = 2e-6
+        );
     }
 
     #[test]
     fn N2() {
         let mut c = Cache::new(cmplx![2., 0.]);
-        assert_approx_eq!(f64, gamma_nsv(&mut c, NF).re, 188.325593, epsilon = 3e-7);
+        assert_approx_eq_cmplx!(
+            f64,
+            gamma_nsv(&mut c, NF),
+            cmplx!(188.325593, 0.),
+            epsilon = 3e-7
+        );
     }
 }
diff --git a/crates/ekore/src/bib.rs b/crates/ekore/src/bib.rs
index 01b7025f8..3e8a1eb8f 100644
--- a/crates/ekore/src/bib.rs
+++ b/crates/ekore/src/bib.rs
@@ -1,4 +1,4 @@
-//! List of References (autogenerated on 2024-07-18T17:42:47.650964).
+//! List of References (autogenerated on 2024-09-10T11:22:16.645053).
 
 #[allow(non_snake_case)]
 /// The Three loop splitting functions in QCD: The Nonsinglet case
@@ -95,3 +95,15 @@ pub fn Ball2015tna() {}
 ///
 /// DOI: [10.1007/BF01245820](https:dx.doi.org/10.1007/BF01245820)
 pub fn Buza1996wv() {}
+
+#[allow(non_snake_case)]
+/// The Gluonic Operator Matrix Elements at O(alpha(s)**2) for DIS Heavy Flavor Production
+///
+/// Bierenbaum, Isabella and Blumlein, Johannes and Klein, Sebastian
+///
+/// Published in: Phys. Lett. B 672 (2009), 401--406
+///
+/// e-Print: [0901.0669](https://arxiv.org/abs/0901.0669)
+///
+/// DOI: [10.1016/j.physletb.2009.01.057](https:dx.doi.org/10.1016/j.physletb.2009.01.057)
+pub fn Bierenbaum2009zt() {}
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs
index fc17da638..3c57d33ef 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs
@@ -35,6 +35,10 @@ pub fn A_singlet(
     if matching_order_qcd >= 1 {
         A_s[0] = as1::A_singlet(c, nf, L);
     }
+    if matching_order_qcd >= 2 {
+        // TODO recover MSbar mass
+        A_s[1] = as2::A_singlet(c, nf, L, false);
+    }
     A_s
 }
 
@@ -55,5 +59,8 @@ pub fn A_non_singlet(
     if matching_order_qcd >= 1 {
         A_ns[0] = as1::A_ns(c, nf, L);
     }
+    if matching_order_qcd >= 2 {
+        A_ns[1] = as2::A_ns(c, nf, L);
+    }
     A_ns
 }
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
index c2409d2d9..898d3d0c4 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as1.rs
@@ -1,4 +1,4 @@
-//! |NLO| |QCD|
+//! |NLO| |QCD|.
 
 use num::complex::Complex;
 use num::Zero;
@@ -55,63 +55,54 @@ pub fn A_gg(_c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
 }
 
 /// Compute the |NLO| singlet |OME|.
-pub fn A_singlet(c: &mut Cache, _nf: u8, L: f64) -> [[Complex<f64>; 3]; 3] {
+pub fn A_singlet(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 3]; 3] {
     [
-        [A_gg(c, _nf, L), Complex::<f64>::zero(), A_gh(c, _nf, L)],
+        [A_gg(c, nf, L), Complex::<f64>::zero(), A_gh(c, nf, L)],
         [
             Complex::<f64>::zero(),
             Complex::<f64>::zero(),
             Complex::<f64>::zero(),
         ],
-        [A_hg(c, _nf, L), Complex::<f64>::zero(), A_hh(c, _nf, L)],
+        [A_hg(c, nf, L), Complex::<f64>::zero(), A_hh(c, nf, L)],
     ]
 }
 
 /// Compute the |NLO| non-singlet |OME|.
-pub fn A_ns(c: &mut Cache, _nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
+pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
     [
         [Complex::<f64>::zero(), Complex::<f64>::zero()],
-        [Complex::<f64>::zero(), A_hh(c, _nf, L)],
+        [Complex::<f64>::zero(), A_hh(c, nf, L)],
     ]
 }
 
 #[cfg(test)]
 mod tests {
-    use crate::cmplx;
-    use crate::{
-        harmonics::cache::Cache, operator_matrix_elements::unpolarized::spacelike::as1::*,
-    };
-    use float_cmp::assert_approx_eq;
+    use super::*;
+    use crate::harmonics::cache::Cache;
+    use crate::{assert_approx_eq_cmplx, cmplx};
     use num::complex::Complex;
     const NF: u8 = 5;
 
     #[test]
     fn test_momentum_conservation() {
-        const N: Complex<f64> = cmplx![2., 0.];
+        const N: Complex<f64> = cmplx!(2., 0.);
         const L: f64 = 100.;
         let mut c = Cache::new(N);
         let aS1 = A_singlet(&mut c, NF, L);
         // heavy quark momentum conservation
-        assert_approx_eq!(
+        assert_approx_eq_cmplx!(
             f64,
-            (aS1[0][2] + aS1[1][2] + aS1[2][2]).re,
-            0.,
-            epsilon = 1e-10
-        );
-        assert_approx_eq!(
-            f64,
-            (aS1[0][2] + aS1[1][2] + aS1[2][2]).im,
-            0.,
+            (aS1[0][2] + aS1[1][2] + aS1[2][2]),
+            Complex::zero(),
             epsilon = 1e-10
         );
         // gluon momentum conservation
-        assert_approx_eq!(f64, (aS1[0][0] + aS1[1][0] + aS1[2][0]).re, 0.);
-        assert_approx_eq!(f64, (aS1[0][0] + aS1[1][0] + aS1[2][0]).im, 0.);
+        assert_approx_eq_cmplx!(f64, (aS1[0][0] + aS1[1][0] + aS1[2][0]), Complex::zero());
     }
 
     #[test]
     fn test_A1_intrinsic() {
-        const N: Complex<f64> = cmplx![2., 0.];
+        const N: Complex<f64> = cmplx!(2., 0.);
         const L: f64 = 3.0;
         let mut c = Cache::new(N);
         let aNS1i = A_ns(&mut c, NF, L);
@@ -125,16 +116,27 @@ mod tests {
         // Only even moments are available in that code.
         // Note there is a minus sign in the definition of L.
         const L: f64 = 10.;
-        let ref_val_gg = [-6.66667, -6.66667, -6.66667, -6.66667, -6.66667];
-        let ref_val_Hg = [6.66667, 3.66667, 2.61905, 2.05556, 1.69697];
+        let ref_val_gg = [
+            -6.666666667,
+            -6.666666667,
+            -6.666666667,
+            -6.666666667,
+            -6.666666667,
+        ];
+        let ref_val_Hg = [
+            6.666666667,
+            3.666666667,
+            2.61904761905,
+            2.0555555556,
+            1.696969697,
+        ];
 
         for n in 0..4 {
             let N = cmplx![2. * (n as f64) + 2., 0.];
             let mut c = Cache::new(N);
             let aS1 = A_singlet(&mut c, NF, L);
-            // lower numerical accuracy than python?
-            assert_approx_eq!(f64, aS1[0][0].re, ref_val_gg[n], epsilon = 4e-6);
-            assert_approx_eq!(f64, aS1[2][0].re, ref_val_Hg[n], epsilon = 5e-6);
+            assert_approx_eq_cmplx!(f64, aS1[0][0], cmplx!(ref_val_gg[n], 0.), epsilon = 1e-6);
+            assert_approx_eq_cmplx!(f64, aS1[2][0], cmplx!(ref_val_Hg[n], 0.), epsilon = 1e-6);
         }
     }
 }
diff --git a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
index f2cfcf13d..bf74badc8 100644
--- a/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
+++ b/crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as2.rs
@@ -1,4 +1,4 @@
-//! |NNLO| |QCD|
+//! |NNLO| |QCD|.
 
 use num::complex::Complex;
 use num::traits::Pow;
@@ -10,7 +10,8 @@ use crate::harmonics::cache::{Cache, K};
 use crate::operator_matrix_elements::unpolarized::spacelike::as1;
 
 /// |NNLO| light-light non-singlet |OME|.
-/// It is given in Eq.() of
+///
+/// Implements Eq. (B.4) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -40,8 +41,9 @@ pub fn A_qq_ns(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     CF * TR * (a_qq_l2 * L.pow(2) + a_qq_l1 * L + a_qq_l0)
 }
 
-/// |NNLO| heavy-light pure-singlet |OME|
-/// It is given in Eq.() of
+/// |NNLO| heavy-light pure-singlet |OME|.
+///
+/// Implements Eq. (B.1) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S2 = c.get(K::S1);
@@ -79,8 +81,10 @@ pub fn A_hq_ps(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     CF * TR * (a_hq_l2 * L.pow(2) + a_hq_l1 * L + a_hq_l0)
 }
 
-/// |NNLO| heavy-gluon |OME|
-/// It is given in Eq.() of
+/// |NNLO| heavy-gluon |OME|.
+///
+/// Implements Eq. (B.3) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
+/// The expession for ``A_Hg_l0`` comes form [\[Bierenbaum:2009zt\]](crate::bib::Bierenbaum2009zt).
 pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -178,8 +182,9 @@ pub fn A_hg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     a_hg_l2 * L.pow(2) + a_hg_l1 * L + a_hg_l0
 }
 
-/// |NNLO| gluon-quark |OME|
-/// It is given in Eq.() of
+/// |NNLO| gluon-quark |OME|.
+///
+/// Implements Eq. (B.5) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -204,8 +209,9 @@ pub fn A_gq(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     CF * TR * (a_gq_l2 * L.pow(2) + a_gq_l1 * L + a_gq_l0)
 }
 
-/// |NNLO| gluon-gluon |OME|
-/// It is given in Eq.() of
+/// |NNLO| gluon-gluon |OME|.
+///
+/// Implements Eq. (B.7) of [\[Buza:1996wv\]](crate::bib::Buza1996wv).
 pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     let N = c.n();
     let S1 = c.get(K::S1);
@@ -253,17 +259,17 @@ pub fn A_gg(c: &mut Cache, _nf: u8, L: f64) -> Complex<f64> {
     a_gg_l2 * L.pow(2) + a_gg_l1 * L + a_gg_l0
 }
 
-/// |NNLO| singlet |OME|
-pub fn A_singlet(c: &mut Cache, _nf: u8, L: f64, is_msbar: bool) -> [[Complex<f64>; 3]; 3] {
-    let A_hq_2 = A_hq_ps(c, _nf, L);
-    let A_qq_2 = A_qq_ns(c, _nf, L);
-    let mut A_hg_2 = A_hg(c, _nf, L);
-    let A_gq_2 = A_gq(c, _nf, L);
-    let mut A_gg_2 = A_gg(c, _nf, L);
-
-    if is_msbar {
-        A_hg_2 -= 2.0 * 4.0 * CF * as1::A_hg(c, _nf, 1.0);
-        A_gg_2 -= 2.0 * 4.0 * CF * as1::A_gg(c, _nf, 1.0);
+/// |NNLO| singlet |OME|.
+pub fn A_singlet(c: &mut Cache, nf: u8, L: f64, is_msbar_mass: bool) -> [[Complex<f64>; 3]; 3] {
+    let A_hq_2 = A_hq_ps(c, nf, L);
+    let A_qq_2 = A_qq_ns(c, nf, L);
+    let mut A_hg_2 = A_hg(c, nf, L);
+    let A_gq_2 = A_gq(c, nf, L);
+    let mut A_gg_2 = A_gg(c, nf, L);
+
+    if is_msbar_mass {
+        A_hg_2 -= 2.0 * 4.0 * CF * as1::A_hg(c, nf, 1.0);
+        A_gg_2 -= 2.0 * 4.0 * CF * as1::A_gg(c, nf, 1.0);
     }
 
     [
@@ -273,10 +279,10 @@ pub fn A_singlet(c: &mut Cache, _nf: u8, L: f64, is_msbar: bool) -> [[Complex<f6
     ]
 }
 
-/// |NNLO| non-singlet |OME|
-pub fn A_ns(c: &mut Cache, _nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
+/// |NNLO| non-singlet |OME|.
+pub fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
     [
-        [A_qq_ns(c, _nf, L), Complex::<f64>::zero()],
+        [A_qq_ns(c, nf, L), Complex::<f64>::zero()],
         [Complex::<f64>::zero(), Complex::<f64>::zero()],
     ]
 }
diff --git a/crates/ekore/src/util.rs b/crates/ekore/src/util.rs
index 0d815bc98..2f70faea1 100644
--- a/crates/ekore/src/util.rs
+++ b/crates/ekore/src/util.rs
@@ -12,16 +12,8 @@ macro_rules! cmplx {
 #[cfg(test)]
 #[macro_export]
 macro_rules! assert_approx_eq_cmplx {
-    ($size:ty, $ref:expr, $target:expr) => {
-        float_cmp::assert_approx_eq!($size, $ref.re, $target.re);
-        float_cmp::assert_approx_eq!($size, $ref.im, $target.im);
-    };
-    ($size:ty, $ref:expr, $target:expr, ulps=$ulps:expr) => {
-        float_cmp::assert_approx_eq!($size, $ref.re, $target.re, ulps = $ulps);
-        float_cmp::assert_approx_eq!($size, $ref.im, $target.im, ulps = $ulps);
-    };
-    ($size:ty, $ref:expr, $target:expr, epsilon=$epsilon:expr) => {
-        float_cmp::assert_approx_eq!($size, $ref.re, $target.re, epsilon = $epsilon);
-        float_cmp::assert_approx_eq!($size, $ref.im, $target.im, epsilon = $epsilon);
+    ($size:ty, $ref:expr, $target:expr $(, $set:ident = $val:expr)*) => {
+        float_cmp::assert_approx_eq!($size, $ref.re, $target.re $(, $set = $val)*);
+        float_cmp::assert_approx_eq!($size, $ref.im, $target.im $(, $set = $val)*);
     };
 }