Merge pull request #369 from NNPDF/rustify-ekore-2

Rustify ekore part 2: ad.u.s.{as2,as3}

crates/eko/src/bib.rs

@@ -0,0 +1 @@
+../../ekore/src/bib.rs

crates/eko/src/lib.rs

@@ -4,7 +4,8 @@ use ekore;
 use ekore::harmonics::cache::Cache;
 use std::ffi::c_void;
 
-mod mellin;
+pub mod bib;
+pub mod mellin;
 
 /// QCD intergration kernel inside quad.
 #[no_mangle]

crates/eko/src/mellin.rs

@@ -9,10 +9,13 @@ use std::f64::consts::PI;
 /// Talbot inversion path.
 ///
 /// Implements the algorithm presented in [\[Abate\]](crate::bib::Abate).
-/// $p_{\text{Talbot}}(t) =  o + r \cdot ( \theta \cot(\theta) + i\theta)$ with $\theta = \pi(2t-1)$
-/// The default values for the parameters $r,o$ are given by $r = 1/2, o = 0$ for
-/// the non-singlet integrals and by $r = \frac{2}{5} \frac{16}{1 - \ln(x)}, o = 1$
-/// for the singlet sector. Note that the non-singlet kernels evolve poles only up to
+///
+/// $$p_{\text{Talbot}}(t) =  o + r \cdot ( \theta \cot(\theta) + i\theta) ~ \text{with}~ \theta = \pi(2t-1)$$
+///
+/// The default values for the parameters $r,o$ are given by
+/// $r = \frac{2}{5} \frac{16}{0.1 - \ln(x)}$ and $o = 0$ for
+/// the non-singlet integrals and $o = 1$ for the singlet sector.
+/// Note that the non-singlet kernels evolve poles only up to
 /// $N=0$ whereas the singlet kernels have poles up to $N=1$.
 pub struct TalbotPath {
     /// integration variable
@@ -26,16 +29,20 @@ pub struct TalbotPath {
 }
 
 impl TalbotPath {
-    /// Auxilary angle.
+    /// Auxiliary angle.
     fn theta(&self) -> f64 {
         PI * (2.0 * self.t - 1.0)
     }
 
     /// Constructor from parameters.
     pub fn new(t: f64, logx: f64, is_singlet: bool) -> Self {
+        // The prescription suggested by Abate for r is 0.4 * M / ( - logx)
+        // with M the number of accurate digits; Maria Ubiali suggested in her thesis M = 16.
+        // However, this seems to yield unstable results for the OME in the large x region
+        // so we add an additional regularization, which makes the path less "edgy" there.
         Self {
             t,
-            r: 0.4 * 16.0 / (1.0 - logx),
+            r: 0.4 * 16.0 / (0.1 - logx),
             o: if is_singlet { 1. } else { 0. },
         }
     }

crates/ekore/refs.bib

@@ -50,3 +50,22 @@ @article{Abate
       journal = {International Journal for Numerical Methods in Engineering},
       doi = {10.1002/nme.995}
 }
+@phdthesis{MuselliPhD,
+    author = "Muselli, Claudio",
+    title = "{Resummations of Transverse Momentum Distributions}",
+    doi = "10.13130/muselli-claudio_phd2017-12-01",
+    school = "Università degli studi di Milano, Dipartimento di Fisica",
+    year = "2017"
+}
+@article{Vogt2004ns,
+    author = "Vogt, A.",
+    archivePrefix = "arXiv",
+    doi = "10.1016/j.cpc.2005.03.103",
+    eprint = "hep-ph/0408244",
+    journal = "Comput. Phys. Commun.",
+    pages = "65--92",
+    reportNumber = "NIKHEF-04-011",
+    title = "{Efficient evolution of unpolarized and polarized parton distributions with QCD-PEGASUS}",
+    volume = "170",
+    year = "2005"
+}

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

@@ -1,14 +1,37 @@
 //! The unpolarized, space-like anomalous dimensions at various couplings power.
 
+use crate::constants::{PID_NSM, PID_NSP, PID_NSV};
 use crate::harmonics::cache::Cache;
 use num::complex::Complex;
 use num::Zero;
 pub mod as1;
+pub mod as2;
+pub mod as3;
 
 /// Compute the tower of the non-singlet anomalous dimensions.
-pub fn gamma_ns_qcd(order_qcd: usize, _mode: u16, c: &mut Cache, nf: u8) -> Vec<Complex<f64>> {
+pub fn gamma_ns_qcd(order_qcd: usize, mode: u16, c: &mut Cache, nf: u8) -> Vec<Complex<f64>> {
     let mut gamma_ns = vec![Complex::<f64>::zero(); order_qcd];
     gamma_ns[0] = as1::gamma_ns(c, nf);
+    // NLO and beyond
+    if order_qcd >= 2 {
+        let gamma_ns_1 = match mode {
+            PID_NSP => as2::gamma_nsp(c, nf),
+            // To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
+            PID_NSM | PID_NSV => as2::gamma_nsm(c, nf),
+            _ => panic!("Unkown non-singlet sector element"),
+        };
+        gamma_ns[1] = gamma_ns_1
+    }
+    // NNLO and beyond
+    if order_qcd >= 3 {
+        let gamma_ns_2 = match mode {
+            PID_NSP => as3::gamma_nsp(c, nf),
+            PID_NSM => as3::gamma_nsm(c, nf),
+            PID_NSV => as3::gamma_nsv(c, nf),
+            _ => panic!("Unkown non-singlet sector element"),
+        };
+        gamma_ns[2] = gamma_ns_2
+    }
     gamma_ns
 }
 
@@ -22,5 +45,13 @@ pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Compl
         order_qcd
     ];
     gamma_S[0] = as1::gamma_singlet(c, nf);
+    // NLO and beyond
+    if order_qcd >= 2 {
+        gamma_S[1] = as2::gamma_singlet(c, nf);
+    }
+    // NNLO and beyond
+    if order_qcd >= 3 {
+        gamma_S[2] = as3::gamma_singlet(c, nf);
+    }
     gamma_S
 }