diff --git a/.github/workflows/lha_bot_rust.yml b/.github/workflows/lha_bot_rust.yml index 5de94c7e0..04d1293de 100644 --- a/.github/workflows/lha_bot_rust.yml +++ b/.github/workflows/lha_bot_rust.yml @@ -17,11 +17,6 @@ jobs: lhabench: name: LHA paper Benchmarks runs-on: ubuntu-latest - # container: - # image: ghcr.io/nnpdf/bench-evol:v2 - # credentials: - # username: ${{ github.repository_owner }} - # password: ${{ secrets.GITHUB_TOKEN }} strategy: matrix: @@ -30,9 +25,6 @@ jobs: steps: - uses: actions/checkout@v2 - # with: - # # tags needed for dynamic versioning - # fetch-depth: 0 - name: Set up Python ${{ matrix.python-version }} 🐍 id: setup-python uses: actions/setup-python@v5 @@ -53,5 +45,4 @@ jobs: ./rustify.sh poe compile poe lha -m "nnlo and sv" - # TODO wait for polarized to reactivate - # poe lha -m "ffns_pol and sv" + poe lha -m "ffns_pol and sv" diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs index 4e6347494..a6c5fa9ae 100644 --- a/crates/eko/src/lib.rs +++ b/crates/eko/src/lib.rs @@ -70,22 +70,21 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 { ); } } else if is_singlet { + let gamma_singlet_qcd = match args.is_polarized { + true => ekore::anomalous_dimensions::polarized::spacelike::gamma_singlet_qcd, + false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd, + }; raw = unravel( - ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd( - args.order_qcd, - &mut c, - args.nf, - ), + gamma_singlet_qcd(args.order_qcd, &mut c, args.nf), args.order_qcd, ); } else { // we can not do 1D - let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd( - args.order_qcd, - args.mode0, - &mut c, - args.nf, - ); + let gamma_ns_qcd = match args.is_polarized { + true => ekore::anomalous_dimensions::polarized::spacelike::gamma_ns_qcd, + false => ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd, + }; + let res = gamma_ns_qcd(args.order_qcd, args.mode0, &mut c, args.nf); for el in res.iter().take(args.order_qcd) { raw.re.push(el.re); raw.im.push(el.im); diff --git a/crates/ekore/src/anomalous_dimensions.rs b/crates/ekore/src/anomalous_dimensions.rs index 55a4dad43..ecb2e4038 100644 --- a/crates/ekore/src/anomalous_dimensions.rs +++ b/crates/ekore/src/anomalous_dimensions.rs @@ -1,3 +1,4 @@ //! The anomalous dimensions for |DGLAP| evolution. +pub mod polarized; pub mod unpolarized; diff --git a/crates/ekore/src/anomalous_dimensions/polarized.rs b/crates/ekore/src/anomalous_dimensions/polarized.rs new file mode 100644 index 000000000..501138f09 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/polarized.rs @@ -0,0 +1,3 @@ +//! The polarized anomalous dimensions for space-like kinematics. + +pub mod spacelike; diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs new file mode 100644 index 000000000..bb394deda --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike.rs @@ -0,0 +1,57 @@ +//! The polarized, 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> { + let mut gamma_ns = vec![Complex::::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 +} + +/// Compute the tower of the singlet anomalous dimension matrices. +pub fn gamma_singlet_qcd(order_qcd: usize, c: &mut Cache, nf: u8) -> Vec<[[Complex; 2]; 2]> { + let mut gamma_S = vec![ + [ + [Complex::::zero(), Complex::::zero()], + [Complex::::zero(), Complex::::zero()] + ]; + 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 +} diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs new file mode 100644 index 000000000..2eec91867 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as1.rs @@ -0,0 +1,85 @@ +//! |LO| |QCD|. + +use num::complex::Complex; + +use super::super::super::unpolarized::spacelike::as1::gamma_ns as unpol; +use crate::constants::{CA, CF, TR}; +use crate::harmonics::cache::{Cache, K}; + +/// Compute the non-singlet anomalous dimension. +/// +/// Identical to the unpolarized counterpart. +pub fn gamma_ns(c: &mut Cache, nf: u8) -> Complex { + unpol(c, nf) +} + +/// Compute the quark-gluon anomalous dimension. +/// +/// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let gamma = -(N - 1.) / N / (N + 1.); + 2.0 * TR * 2.0 * (nf as f64) * gamma +} + +/// Compute the gluon-quark anomalous dimension. +/// +/// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_gq(c: &mut Cache, _nf: u8) -> Complex { + let N = c.n(); + let gamma = -(N + 2.) / N / (N + 1.); + 2.0 * CF * gamma +} + +/// Compute the gluon-gluon anomalous dimension. +/// +/// Implements Eq. (A.1) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { + let N = c.n(); + let S1 = c.get(K::S1); + let gamma = -S1 + 2. / N / (N + 1.); + CA * (-4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * TR * (nf as f64) +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let gamma_qq = gamma_ns(c, nf); + [ + [gamma_qq, gamma_qg(c, nf)], + [gamma_gq(c, nf), gamma_gg(c, nf)], + ] +} + +#[cfg(test)] +mod tests { + 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] + fn quark_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + let me = gamma_ns(&mut c, NF) + gamma_gq(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, cmplx!((4. * CF) / 3., 0.), epsilon = 1e-12); + } + + #[test] + fn gluon_momentum_conservation() { + const N: Complex = cmplx!(2., 0.); + let mut c = Cache::new(N); + let me = gamma_qg(&mut c, NF) + gamma_gg(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, cmplx!(3. + (NF as f64) / 3., 0.), epsilon = 1e-12); + } + + #[test] + fn qg_helicity_conservation() { + const N: Complex = cmplx!(1., 0.); + let mut c = Cache::new(N); + let me = gamma_qg(&mut c, NF); + assert_approx_eq_cmplx!(f64, me, Complex::zero(), epsilon = 1e-12); + } +} diff --git a/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs new file mode 100644 index 000000000..8dcdf6733 --- /dev/null +++ b/crates/ekore/src/anomalous_dimensions/polarized/spacelike/as2.rs @@ -0,0 +1,195 @@ +//! |NLO| |QCD|. + +use num::complex::Complex; + +use super::super::super::unpolarized::spacelike::as2::gamma_nsm as unpol_nsm; +use super::super::super::unpolarized::spacelike::as2::gamma_nsp as unpol_nsp; +use crate::constants::{CA, CF, TR, ZETA2, ZETA3}; +use crate::harmonics::cache::{Cache, K}; + +pub fn gamma_nsm(c: &mut Cache, nf: u8) -> Complex { + unpol_nsp(c, nf) +} + +pub fn gamma_nsp(c: &mut Cache, nf: u8) -> Complex { + unpol_nsm(c, nf) +} + +/// Compute the pure-singlet quark-quark anomalous dimension. +/// +/// Implements Eq. (A.3) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_ps(c: &mut Cache, nf: u8) -> Complex { + let n = c.n(); + let gqqps1_nfcf = (2. * (n + 2.) * (1. + 2. * n + n.powu(3))) / ((1. + n).powu(3) * n.powu(3)); + 4.0 * TR * (nf as f64) * CF * gqqps1_nfcf +} + +/// Compute the quark-gluon singlet anomalous dimension. +/// +/// Implements Eq. (A.4) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex { + let n = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let Sp2m = c.get(K::S2mh); + + #[rustfmt::skip] + let gqg1_nfca = ( + (S1.powu(2) - S2 + Sp2m) * (n - 1.) / (n * (n + 1.)) + - 4. * S1 / (n * (1. + n).powu(2)) + - (-2. - 7. * n + 3. * n.powu(2) - 4. * n.powu(3) + n.powu(4) + n.powu(5)) / (n.powu(3) * (1. + n).powu(3)) + ) * 2.0; + #[rustfmt::skip] + let gqg1_nfcf = ( + (-(S1.powu(2)) + S2 + 2. * S1 / n) * (n - 1.) / (n * (n + 1.)) + - (n - 1.) + * (1. + 3.5 * n + 4. * n.powu(2) + 5. * n.powu(3) + 2.5 * n.powu(4)) + / (n.powu(3) * (1. + n).powu(3)) + + 4. * (n - 1.) / (n.powu(2) * (1. + n).powu(2)) + ) * 2.; + 4.0 * TR * (nf as f64) * (CA * gqg1_nfca + CF * gqg1_nfcf) +} + +/// Compute the gluon-quark singlet anomalous dimension. +/// +/// Implements Eq. (A.5) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex { + let n = c.n(); + let S1 = c.get(K::S1); + let S2 = c.get(K::S2); + let Sp2m = c.get(K::S2mh); + #[rustfmt::skip] + let ggq1_cfcf = ( + (2. * (S1.powu(2) + S2) * (n + 2.)) / (n * (n + 1.)) + - (2. * S1 * (n + 2.) * (1. + 3. * n)) / (n * (1. + n).powu(2)) + - ((n + 2.) * (2. + 15. * n + 8. * n.powu(2) - 12.0 * n.powu(3) - 9.0 * n.powu(4))) + / (n.powu(3) * (1. + n).powu(3)) + + 8. * (n + 2.) / (n.powu(2) * (1. + n).powu(2)) + ) * 0.5; + #[rustfmt::skip] + let ggq1_cfca = -( + -(-(S1.powu(2)) - S2 + Sp2m) * (n + 2.) / (n * (n + 1.)) + - S1 * (12. + 22. * n + 11. * n.powu(2)) / (3. * n.powu(2) * (n + 1.)) + + (36. + 72. * n + 41. * n.powu(2) + 254. * n.powu(3) + 271. * n.powu(4) + 76. * n.powu(5)) + / (9. * n.powu(3) * (1. + n).powu(3)) + ); + #[rustfmt::skip] + let ggq1_cfnf =(-S1 * (n + 2.)) / (3. * n * (n + 1.)) + ((n + 2.) * (2. + 5. * n)) / ( + 9. * n * (1. + n).powu(2) + ); + 4. * CF * (CA * ggq1_cfca + CF * ggq1_cfcf + 4.0 * TR * (nf as f64) * ggq1_cfnf) +} + +/// Compute the gluon-gluon singlet anomalous dimension. +/// +/// Implements Eq. (A.6) of [\[Gluck:1995yr\]][crate::bib::Gluck1995yr]. +pub fn gamma_gg(c: &mut Cache, nf: u8) -> Complex { + let n = c.n(); + let S1 = c.get(K::S1); + let Sp1m = c.get(K::S1mh); + let Sp2m = c.get(K::S2mh); + let Sp3m = c.get(K::S3mh); + let S1h = c.get(K::S1h); + let g3 = c.get(K::G3); + let SSCHLM = ZETA2 / 2. * (Sp1m - S1h + 2. / n) - S1 / n.powu(2) - g3 - 5. * ZETA3 / 8.; + #[rustfmt::skip] + let ggg1_caca = ( + -4. * S1 * Sp2m + - Sp3m + + 8. * SSCHLM + + 8. * Sp2m / (n * (n + 1.)) + + 2.0 + * S1 + * (72. + 144. * n + 67. * n.powu(2) + 134. * n.powu(3) + 67. * n.powu(4)) + / (9. * n.powu(2) * (n + 1.).powu(2)) + - (144. + 258. * n + 7. * n.powu(2) + 698. * n.powu(3) + 469. * n.powu(4) + 144. * n.powu(5) + 48. * n.powu(6)) + / (9. * n.powu(3) * (1. + n).powu(3)) + ) * 0.5; + #[rustfmt::skip] + let ggg1_canf = ( + -5. * S1 / 9. + + (-3. + 13. * n + 16. * n.powu(2) + 6.* n.powu(3) + 3. * n.powu(4)) / (9. * n.powu(2) * (1. + n).powu(2)) + ) * 4.; + #[rustfmt::skip] + let ggg1_cfnf = (4. + 2. * n - 8. * n.powu(2) + n.powu(3) + 5. * n.powu(4) + 3. * n.powu(5) + n.powu(6)) / ( + n.powu(3) * (1. + n).powu(3) + ); + 4. * (CA * CA * ggg1_caca + TR * (nf as f64) * (CA * ggg1_canf + CF * ggg1_cfnf)) +} + +/// Compute the singlet anomalous dimension matrix. +pub fn gamma_singlet(c: &mut Cache, nf: u8) -> [[Complex; 2]; 2] { + let gamma_qq = gamma_nsp(c, nf) + gamma_ps(c, nf); + [ + [gamma_qq, gamma_qg(c, nf)], + [gamma_gq(c, nf), gamma_gg(c, nf)], + ] +} + +#[cfg(test)] +mod tests { + 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; + + #[test] + fn physical_constraints() { + // qg_helicity_conservation + let mut c = Cache::new(cmplx!(1., 0.)); + assert_approx_eq_cmplx!(f64, gamma_qg(&mut c, NF), Complex::zero(), epsilon = 2e-6); + + // qg momentum + let mut c = Cache::new(cmplx!(1., 0.)); + let gS1 = gamma_singlet(&mut c, NF); + assert_approx_eq_cmplx!( + f64, + gS1[0][0], + cmplx!(12. * TR * (NF as f64) * CF, 0.), + epsilon = 2e-6 + ); + } + + #[test] + fn N2() { + let mut c = Cache::new(cmplx!(2., 0.)); + + // singlet sector + let gS1 = gamma_singlet(&mut c, NF); + // ps + assert_approx_eq_cmplx!( + f64, + -gamma_ps(&mut c, NF), + cmplx!(-4.0 * CF * TR * (NF as f64) * 13. / 27.0, 0.) + ); + // qg + assert_approx_eq_cmplx!( + f64, + -gS1[0][1], + cmplx!( + 4. * (NF as f64) + * (0.574074074 * CF - 2. * CA * (-7. / 18. + 1. / 6. * (5. - PI.pow(2) / 3.))) + * TR, + 0. + ), + epsilon = 1e-9 + ); + // gq + assert_approx_eq_cmplx!( + f64, + -gS1[1][0], + cmplx!( + 4. * (-2.074074074074074 * CF.pow(2) + + CA * CF * (29. / 54. - 2. / 3. * (1. / 2. - PI.pow(2) / 3.)) + + (4. * CF * (NF as f64) * TR) / 27.), + 0. + ), + epsilon = 1e-13 + ); + } +}