Mock external libome

Cargo.toml

@@ -39,6 +39,7 @@ yaml-rust2 = "0.8"
 
 # Dev dependencies
 assert_fs = "1.1.2"
+cc = "1.0"
 float-cmp = "0.9.0"
 GSL = "7.0"
 predicates = "3.1.2"

crates/ekore/Cargo.toml

@@ -20,6 +20,9 @@ rustdoc-args = ["--html-in-header", "doc-header.html"]
 [dependencies]
 num.workspace = true
 
+[build-dependencies]
+cc.workspace = true
+
 [dev-dependencies]
 float-cmp.workspace = true
 GSL.workspace = true

crates/ekore/build.rs

@@ -0,0 +1,26 @@
+//! Build script to compile the mock `libome` C++ library from `extras/gsoc/libome`
+//! and link it statically into `ekore`.
+
+use std::path::PathBuf;
+
+fn main() {
+    let manifest_dir = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
+    let libome_dir = manifest_dir.join("../../extras/gsoc/libome");
+
+    println!(
+        "cargo:rerun-if-changed={}",
+        libome_dir.join("ome.cpp").display()
+    );
+    println!(
+        "cargo:rerun-if-changed={}",
+        libome_dir.join("ome.h").display()
+    );
+
+    cc::Build::new()
+        .cpp(true)
+        .flag_if_supported("-std=c++11")
+        .flag_if_supported("/std:c++14")
+        .include(&libome_dir)
+        .file(libome_dir.join("ome.cpp"))
+        .compile("ome");
+}

crates/ekore/src/operator_matrix_elements/unpolarized/spacelike.rs

@@ -5,20 +5,22 @@ use num::Zero;
 use num::complex::Complex;
 mod as1;
 mod as2;
+mod as3;
 
 /// Compute the tower of the unpolarized, space-like singlet |OME|.
 ///
 /// This function computes the first `matching_order_qcd` entries at the Mellin variable set in `cache`
 /// using `nf` light flavors and `L`, logarithm of the squared ratio of factorization scale and mass
 /// $\ln(\mu_F^2 / m^2)$.
 ///
-/// Returns an array of shape `(MAX_ORDER_QCD - 2, 3, 3)`. Only the first `matching_order_qcd`
+/// Returns an array of shape `(MAX_ORDER_QCD - 1, 3, 3)`. Only the first `matching_order_qcd`
 /// entries along the outer axis are filled; remaining slots are zero.
 ///
 /// # Available perturbative orders:
 ///
 /// - $a_s^1$ [\[Ball:2015tna\]](crate::bib::Ball2015tna) [\[Buza:1996wv\]](crate::bib::Buza1996wv)
 /// - $a_s^2$ [\[Buza:1996wv\]](crate::bib::Buza1996wv) [\[Bierenbaum:2009zt\]](crate::bib::Bierenbaum2009zt)
+/// - $a_s^3$ (via external `libome`)
 ///
 /// # Panics
 ///
@@ -28,18 +30,21 @@ pub fn A_singlet(
     cache: &mut Cache,
     nf: u8,
     L: f64,
-) -> [[[Complex<f64>; 3]; 3]; MAX_ORDER_QCD - 2] {
-    if matching_order_qcd >= 3 {
-        panic!("OME beyond NLO is not yet implemented");
+) -> [[[Complex<f64>; 3]; 3]; MAX_ORDER_QCD - 1] {
+    if matching_order_qcd >= 4 {
+        panic!("OME beyond N3LO is not yet implemented");
     }
-    let mut A_s = [[[Complex::<f64>::zero(); 3]; 3]; MAX_ORDER_QCD - 2];
+    let mut A_s = [[[Complex::<f64>::zero(); 3]; 3]; MAX_ORDER_QCD - 1];
     if matching_order_qcd >= 1 {
         A_s[0] = as1::A_singlet(cache, nf, L);
     }
     if matching_order_qcd >= 2 {
         // TODO recover MSbar mass
         A_s[1] = as2::A_singlet(cache, nf, L, false);
     }
+    if matching_order_qcd >= 3 {
+        A_s[2] = as3::A_singlet(cache, nf, L);
+    }
     A_s
 }
 
@@ -49,13 +54,14 @@ pub fn A_singlet(
 /// using `nf` light flavors and `L`, logarithm of the squared ratio of factorization scale and mass
 /// $\ln(\mu_F^2 / m^2)$.
 ///
-/// Returns an array of shape `(MAX_ORDER_QCD - 2, 2, 2)`. Only the first `matching_order_qcd`
+/// Returns an array of shape `(MAX_ORDER_QCD - 1, 2, 2)`. Only the first `matching_order_qcd`
 /// entries along the outer axis are filled; remaining slots are zero.
 ///
 /// # Available perturbative orders:
 ///
 /// - $a_s^1$ [\[Ball:2015tna\]](crate::bib::Ball2015tna)
 /// - $a_s^2$ [\[Buza:1996wv\]](crate::bib::Buza1996wv)
+/// - $a_s^3$ (via external `libome`)
 ///
 /// # Panics
 ///
@@ -65,17 +71,20 @@ pub fn A_non_singlet(
     cache: &mut Cache,
     nf: u8,
     L: f64,
-) -> [[[Complex<f64>; 2]; 2]; MAX_ORDER_QCD - 2] {
-    if matching_order_qcd >= 3 {
-        panic!("OME beyond NLO is not yet implemented");
+) -> [[[Complex<f64>; 2]; 2]; MAX_ORDER_QCD - 1] {
+    if matching_order_qcd >= 4 {
+        panic!("OME beyond N3LO is not yet implemented");
     }
-    let mut A_ns = [[[Complex::<f64>::zero(); 2]; 2]; MAX_ORDER_QCD - 2];
+    let mut A_ns = [[[Complex::<f64>::zero(); 2]; 2]; MAX_ORDER_QCD - 1];
     if matching_order_qcd >= 1 {
         A_ns[0] = as1::A_ns(cache, nf, L);
     }
     if matching_order_qcd >= 2 {
         A_ns[1] = as2::A_ns(cache, nf, L);
     }
+    if matching_order_qcd >= 3 {
+        A_ns[2] = as3::A_ns(cache, nf, L);
+    }
     A_ns
 }
 
@@ -90,18 +99,18 @@ mod tests {
         const NF: u8 = 4;
         const N: Complex<f64> = cmplx!(0., 1.);
         const L: f64 = 0.0;
-        for matching_order_qcd in 1..=2usize {
+        for matching_order_qcd in 1..=3usize {
             let mut cache = Cache::new(N);
             let a_s = A_singlet(matching_order_qcd, &mut cache, NF, L);
-            assert_eq!(a_s.len(), MAX_ORDER_QCD - 2);
+            assert_eq!(a_s.len(), MAX_ORDER_QCD - 1);
             assert_eq!(a_s[0].len(), 3);
             assert_eq!(a_s[0][0].len(), 3);
             // slots beyond order_qcd must be zero
             for item in a_s.iter().skip(matching_order_qcd) {
                 assert_approx_eq_cmplx_2d!(f64, item, [[cmplx!(0., 0.); 3]; 3], 3);
             }
             let a_ns = A_non_singlet(matching_order_qcd, &mut cache, NF, L);
-            assert_eq!(a_ns.len(), MAX_ORDER_QCD - 2);
+            assert_eq!(a_ns.len(), MAX_ORDER_QCD - 1);
             assert_eq!(a_ns[0].len(), 2);
             assert_eq!(a_ns[0][0].len(), 2);
             // slots beyond order_qcd must be zero
@@ -118,11 +127,13 @@ mod tests {
         const N: Complex<f64> = cmplx!(1., 0.);
         const L: f64 = 0.0;
         let mut cache = Cache::new(N);
-        let a_ns = A_non_singlet(2, &mut cache, NF, L);
+        let a_ns = A_non_singlet(3, &mut cache, NF, L);
         // LO
         assert_approx_eq_cmplx_2d!(f64, a_ns[0], [[cmplx!(0., 0.); 2]; 2], 2, epsilon = 1e-14);
         // NNLO
         assert_approx_eq_cmplx_2d!(f64, a_ns[1], [[cmplx!(0., 0.); 2]; 2], 2, epsilon = 1e-14);
+        // N3LO
+        assert_approx_eq_cmplx_2d!(f64, a_ns[2], [[cmplx!(0., 0.); 2]; 2], 2, epsilon = 1e-14);
     }
 
     #[test]
@@ -132,7 +143,7 @@ mod tests {
         const N: Complex<f64> = cmplx!(2., 0.);
         const L: f64 = 100.;
         let mut cache = Cache::new(N);
-        let a_s = A_singlet(2, &mut cache, NF, L);
+        let a_s = A_singlet(3, &mut cache, NF, L);
         // LO
         assert_approx_eq_cmplx!(
             f64,
@@ -158,25 +169,27 @@ mod tests {
             Complex::zero(),
             epsilon = 1e-11
         );
+        // N3LO
+        assert_approx_eq_cmplx_2d!(f64, a_s[2], [[cmplx!(0., 0.); 3]; 3], 3, epsilon = 1e-14);
     }
 
     #[test]
-    #[should_panic(expected = "OME beyond NLO is not yet implemented")]
-    fn test_a_singlet_order3_panics() {
+    #[should_panic(expected = "OME beyond N3LO is not yet implemented")]
+    fn test_a_singlet_order4_panics() {
         const NF: u8 = 4;
         const N: Complex<f64> = cmplx!(1.234, 0.);
         const L: f64 = 0.0;
         let mut cache = Cache::new(N);
-        A_singlet(3, &mut cache, NF, L);
+        A_singlet(4, &mut cache, NF, L);
     }
 
     #[test]
-    #[should_panic(expected = "OME beyond NLO is not yet implemented")]
-    fn test_a_non_singlet_order3_panics() {
+    #[should_panic(expected = "OME beyond N3LO is not yet implemented")]
+    fn test_a_non_singlet_order4_panics() {
         const NF: u8 = 4;
         const N: Complex<f64> = cmplx!(1.234, 0.);
         const L: f64 = 0.0;
         let mut cache = Cache::new(N);
-        A_non_singlet(3, &mut cache, NF, L);
+        A_non_singlet(4, &mut cache, NF, L);
     }
 }

crates/ekore/src/operator_matrix_elements/unpolarized/spacelike/as3.rs

@@ -0,0 +1,94 @@
+//! |N3LO| unpolarized, space-like |OME| via external `libome` C ABI.
+
+use std::ffi::{c_int, c_uint};
+
+use num::Zero;
+use num::complex::Complex;
+
+use crate::harmonics::cache::Cache;
+
+#[repr(C)]
+#[derive(Clone, Copy, Debug)]
+struct OmeComplex {
+    re: f64,
+    im: f64,
+}
+
+impl From<Complex<f64>> for OmeComplex {
+    fn from(c: Complex<f64>) -> Self {
+        Self { re: c.re, im: c.im }
+    }
+}
+
+impl From<OmeComplex> for Complex<f64> {
+    fn from(c: OmeComplex) -> Self {
+        Complex::new(c.re, c.im)
+    }
+}
+
+unsafe extern "C" {
+    fn ome_as3_Agg(n: OmeComplex, nf: c_uint, L: f64) -> OmeComplex;
+    fn ome_as3_Agq(n: OmeComplex, nf: c_uint, L: f64) -> OmeComplex;
+    fn ome_as3_Aqg(n: OmeComplex, nf: c_uint, L: f64) -> OmeComplex;
+    fn ome_as3_AHg(n: OmeComplex, nf: c_uint, L: f64) -> OmeComplex;
+    fn ome_as3_AHq(n: OmeComplex, nf: c_uint, L: f64) -> OmeComplex;
+    fn ome_as3_AqqPS(n: OmeComplex, nf: c_uint, L: f64) -> OmeComplex;
+    fn ome_as3_AqqNS(n: OmeComplex, nf: c_uint, L: f64, eta: c_int) -> OmeComplex;
+}
+
+/// Compute the |N3LO| singlet |OME|.
+pub(super) fn A_singlet(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 3]; 3] {
+    let n: OmeComplex = c.n().into();
+    let nf_c = c_uint::from(nf);
+
+    let a_gg = unsafe { Complex::from(ome_as3_Agg(n, nf_c, L)) };
+    let a_gq = unsafe { Complex::from(ome_as3_Agq(n, nf_c, L)) };
+    let a_qg = unsafe { Complex::from(ome_as3_Aqg(n, nf_c, L)) };
+    let a_hg = unsafe { Complex::from(ome_as3_AHg(n, nf_c, L)) };
+    let a_hq = unsafe { Complex::from(ome_as3_AHq(n, nf_c, L)) };
+    let a_qq_ps = unsafe { Complex::from(ome_as3_AqqPS(n, nf_c, L)) };
+    let a_qq_ns = unsafe { Complex::from(ome_as3_AqqNS(n, nf_c, L, 1 as c_int)) };
+
+    [
+        [a_gg, a_gq, Complex::<f64>::zero()],
+        [a_qg, a_qq_ps + a_qq_ns, Complex::<f64>::zero()],
+        [a_hg, a_hq, Complex::<f64>::zero()],
+    ]
+}
+
+/// Compute the |N3LO| non-singlet |OME|.
+pub(super) fn A_ns(c: &mut Cache, nf: u8, L: f64) -> [[Complex<f64>; 2]; 2] {
+    let n: OmeComplex = c.n().into();
+    let nf_c = c_uint::from(nf);
+    let a_qq_ns = unsafe { Complex::from(ome_as3_AqqNS(n, nf_c, L, -1 as c_int)) };
+
+    [
+        [a_qq_ns, Complex::<f64>::zero()],
+        [Complex::<f64>::zero(), Complex::<f64>::zero()],
+    ]
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::cmplx;
+
+    #[test]
+    fn test_as3_calls() {
+        const NF: u8 = 4;
+        let n = cmplx!(2.0, 1.5);
+        let mut c = Cache::new(n);
+        let a_s = A_singlet(&mut c, NF, 0.0);
+        for row in a_s.iter() {
+            for entry in row.iter() {
+                assert_eq!(*entry, Complex::zero());
+            }
+        }
+        let a_ns = A_ns(&mut c, NF, 0.0);
+        for row in a_ns.iter() {
+            for entry in row.iter() {
+                assert_eq!(*entry, Complex::zero());
+            }
+        }
+    }
+}