Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rustify QED #416

Draft
wants to merge 32 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
32 commits
Select commit Hold shift + click to select a range
1762846
start implementation of aem1
tgiani Oct 16, 2024
61220f8
complete aem1.py
tgiani Oct 21, 2024
a40908b
unit tests
tgiani Oct 21, 2024
f0edfae
add struct for charge combinations
tgiani Oct 21, 2024
a8d72c0
some work on gamma_ns_qed
tgiani Oct 23, 2024
e8a3a8b
start as1aem1.rs
tgiani Oct 23, 2024
136dbe8
gamma_qph as1aem1
tgiani Oct 23, 2024
ce11227
some work on as1aem1
tgiani Oct 25, 2024
659fe31
addingv recursive harmonics to cache
tgiani Oct 25, 2024
4007f0d
S1p2 and G3p2. Check that this is correct
tgiani Oct 25, 2024
7c160f1
gamma_nsp
tgiani Oct 25, 2024
37a65e8
rewrite ChargeCombination struct
tgiani Nov 11, 2024
fc5de81
missing entries in as1aem1
tgiani Nov 11, 2024
2a911b1
adding tests
tgiani Nov 11, 2024
bad124a
some work on spacelike.rs and recasting from list to vec
tgiani Nov 20, 2024
d50e9c7
valence qed
tgiani Nov 20, 2024
b0f6722
some unit tests
tgiani Nov 20, 2024
e4c621b
use Vec<> only when dimension is not known, else use normal list
tgiani Nov 20, 2024
81ece22
fix notation row/columns
tgiani Nov 20, 2024
0def938
add unravel functions for qed case. To be checked
tgiani Nov 20, 2024
9b58455
modifying rust_quad_ker_qcd, PyQuadKerQCDT and QuadQCDargs to include…
tgiani Nov 21, 2024
588496a
add qed valence option in rust_quad_ker_qcd
tgiani Nov 21, 2024
3452243
remove useless flag is_qed
tgiani Nov 21, 2024
9335b97
extend QuadQCDargs to include arguments for c_quad_ker_qed
tgiani Nov 21, 2024
3390aca
setting up input vectors for cb_quad_ker_qed
tgiani Nov 28, 2024
6a3af19
cb_quad_ker_qed
tgiani Dec 3, 2024
f0a727d
fix modes for singlet QED and use == in if statements
tgiani Dec 9, 2024
398c89a
small fix
tgiani Dec 9, 2024
ea44cee
cargo fmt
tgiani Dec 9, 2024
5721326
forgot pre-commit
tgiani Dec 9, 2024
83af836
updating patch file
tgiani Dec 11, 2024
1af5cbf
Update crates/ekore/src/constants.rs
tgiani Dec 11, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
//! The $O(a_s^1a_{em}^1)$ Altarelli-Parisi splitting kernels.
use crate::cmplx;
use num::complex::Complex;

use crate::constants::{CA, CF};
use crate::constants::{ed2, eu2, uplike_flavors, CA, CF, NC, TR, ZETA2, ZETA3};
use crate::harmonics::cache::{Cache, K};

/// Compute the $O(a_s^1a_{em}^1)$ photon-quark anomalous dimension.
Expand Down Expand Up @@ -44,6 +45,7 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {
let S1 = c.get(K::S1);
let S2 = c.get(K::S2);

#[rustfmt::skip]
let tmp_const = -2.0
* (4.0
+ 8.0 * N
Expand All @@ -60,3 +62,81 @@ pub fn gamma_qph(c: &mut Cache, nf: u8) -> Complex<f64> {

2.0 * (nf as f64) * CA * CF * (tmp_const + tmp_S1 * S1 + tmp_S12 * S1.powu(2) + tmp_S2 * S2)
}

/// Compute the $O(a_s^1a_{em}^1)$ gluon-photon anomalous dimension.
///
/// Implements Eq. (27) of
pub fn gamma_gph(c: &mut Cache, _nf: u8) -> Complex<f64> {
let N = c.n();
CF * CA
* (8.0 * (-4.0 + N * (-4.0 + N * (-5.0 + N * (-10.0 + N + 2.0 * N.powu(2) * (2.0 + N))))))
/ (N.powu(3) * (1.0 + N).powu(3) * (-2.0 + N + N.powu(2)))
}

/// Compute the $O(a_s^1a_{em}^1)$ photon-gluon anomalous dimension.
///
/// Implements Eq. (30) of
pub fn gamma_phg(c: &mut Cache, nf: u8) -> Complex<f64> {
TR / CF / CA * (NC as f64) * gamma_gph(c, nf)
}

/// Compute the $O(a_s^1a_{em}^1)$ quark-gluon singlet anomalous dimension.
///
/// Implements Eq. (29) of
pub fn gamma_qg(c: &mut Cache, nf: u8) -> Complex<f64> {
TR / CF / CA * (NC as f64) * gamma_qph(c, nf)
}

/// Compute the $O(a_s^1a_{em}^1)$ gluon-quark singlet anomalous dimension.
///
/// Implements Eq. (35) of
pub fn gamma_gq(c: &mut Cache, nf: u8) -> Complex<f64> {
gamma_phq(c, nf)
}

/// Compute the $O(a_s^1a_{em}^1)$ photon-photon singlet anomalous dimension.
///
/// Implements Eq. (28) of
pub fn gamma_phph(_c: &mut Cache, nf: u8) -> Complex<f64> {
let nu = uplike_flavors(nf);
let nd = nf - nu;
cmplx!(4.0 * CF * CA * ((nu as f64) * eu2 + (nd as f64) * ed2), 0.)
}

/// Compute the $O(a_s^1a_{em}^1)$ gluon-gluon singlet anomalous dimension.
///
/// Implements Eq. (31) of
pub fn gamma_gg(_c: &mut Cache, _nf: u8) -> Complex<f64> {
cmplx!(4.0 * TR * (NC as f64), 0.)
}

/// Compute the $O(a_s^1a_{em}^1)$ singlet-like non singlet anomalous dimension.
///
/// Implements Eqs. (33-34) of
pub fn gamma_nsp(c: &mut Cache, _nf: u8) -> Complex<f64> {
let N = c.n();
let S1 = c.get(K::S1);
let S2 = c.get(K::S2);
let S3 = c.get(K::S3);
let S1h = c.get(K::S1h);
let S2h = c.get(K::S2h);
let S3h = c.get(K::S3h);
let S1p1h = c.get(K::S1ph);
let S2p1h = c.get(K::S2ph);
let S3p1h = c.get(K::S3ph);
let g3N = c.get(K::G3);
let g3Np2 = c.get(K::G3p2);

#[rustfmt::skip]
let result = 32.0 * ZETA2 * S1h - 32.0 * ZETA2 * S1p1h
+ 8.0 / (N + N.powu(2)) * S2h
- 4.0 * S3h + (24.0 + 16.0 / (N + N.powu(2))) * S2
- 32.0 * S3 - 8.0 / (N + N.powu(2)) * S2p1h
+ S1 * (16.0 * (3.0 / N.powu(2) - 3.0 / (1.0 + N).powu(2) + 2.0 * ZETA2) - 16.0 * S2h
- 32.0 * S2 + 16.0 * S2p1h )
+ (-8.0 + N * (-32.0 + N * ( -8.0 - 3.0 * N * (3.0 + N) * (3.0 + N.powu(2)) - 48.0 * (1.0 + N).powu(2) * ZETA2)))
/ (N.powu(3) * (1.0 + N).powu(3))
+ 32.0 * (g3N + g3Np2) + 4.0 * S3p1h - 16.0 * ZETA3;

CF * result
}
11 changes: 11 additions & 0 deletions crates/ekore/src/harmonics/cache.rs
Original file line number Diff line number Diff line change
Expand Up @@ -46,6 +46,12 @@ pub enum K {
Sm21e,
/// $S_{-2,1}(N)$ odd moments
Sm21o,
/// recursive harmonics
S1ph,
S2ph,
S3ph,
S1p2,
G3p2,
Comment on lines +49 to +54
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure we want to keep those - can we inline them?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

what do you mean? Could you please provide an example?

Copy link
Contributor

@felixhekhorn felixhekhorn Oct 28, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

instead of

    let g3N = c.get(K::G3);
    let g3Np2 = c.get(K::G3p2);

we do

    let g3N = c.get(K::G3);
    let g3Np2 = g3N + f;

where $f = g_3(N+2)-g_3(N)$ is computed analytically beforehand

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am afraid I don t understand. How do you compute f analytically beforehand? Do you mean to move g_functions::g3(self.n + 2., self.get(K::S1p2)) out of the cache?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

out of the cache?

yes

How do you compute f analytically beforehand?

Mathematica 🙃 let's try to use a short cut: @giacomomagni can you compute quickly this $f$? (also feel free to give your opinion if we should keep that cache entry or not)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

so I think
$$f = \frac{6 (n+1) (2 n+1) H_n+n \left(-\pi ^2 (n+1)^2-6 n\right)}{6 n^2 (n+1)^3}$$
which I computed half with pen and paper and half Mathematica and checked for $N=1$, $N=2$

to be more precise: for $g_3$ we could use the above equation and for the others we inline the expressions

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ok thank you! I ll first try to get things working as they are and then I will add this

}

/// Hold all cached values.
Expand Down Expand Up @@ -90,6 +96,7 @@ impl Cache {
K::S2mh => w2::S2((self.n - 1.) / 2.),
K::S3mh => w3::S3((self.n - 1.) / 2.),
K::G3 => g_functions::g3(self.n, self.get(K::S1)),
K::G3p2 => g_functions::g3(self.n + 2., self.get(K::S1p2)),
K::Sm1e => w1::Sm1e(self.get(K::S1), self.get(K::S1h)),
K::Sm1o => w1::Sm1o(self.get(K::S1), self.get(K::S1mh)),
K::Sm2e => w2::Sm2e(self.get(K::S2), self.get(K::S2h)),
Expand All @@ -98,6 +105,10 @@ impl Cache {
K::Sm3o => w3::Sm3o(self.get(K::S3), self.get(K::S3mh)),
K::Sm21e => w3::Sm21e(self.n, self.get(K::S1), self.get(K::Sm1e)),
K::Sm21o => w3::Sm21o(self.n, self.get(K::S1), self.get(K::Sm1o)),
K::S1ph => recursive_harmonic_sum(self.get(K::S1mh), (self.n - 1.) / 2., 1, 1),
K::S2ph => recursive_harmonic_sum(self.get(K::S2mh), (self.n - 1.) / 2., 1, 2),
K::S3ph => recursive_harmonic_sum(self.get(K::S3mh), (self.n - 1.) / 2., 1, 3),
K::S1p2 => recursive_harmonic_sum(self.get(K::S1), self.n, 2, 1),
};
// insert
self.m.insert(k, val);
Expand Down
Loading