diff --git a/extras/harmonics_w5/src/harmonics_w5/__init__.py b/extras/harmonics_w5/src/harmonics_w5/__init__.py index af6ccfebe..8a599117d 100644 --- a/extras/harmonics_w5/src/harmonics_w5/__init__.py +++ b/extras/harmonics_w5/src/harmonics_w5/__init__.py @@ -2,9 +2,9 @@ import numba as nb +from eko.constants import log2, zeta2, zeta3, zeta4, zeta5 from ekore.harmonics import S5, Sm5 -from ekore.harmonics.constants import log2, zeta2, zeta3, zeta4, zeta5 -from ekore.harmonics.polygamma import cern_polygamma, symmetry_factor +from ekore.harmonics.polygamma import symmetry_factor from . import f_functions as f diff --git a/extras/harmonics_w5/tests/test_f_functions.py b/extras/harmonics_w5/tests/test_f_functions.py index 0259adda7..fcae8c5fd 100644 --- a/extras/harmonics_w5/tests/test_f_functions.py +++ b/extras/harmonics_w5/tests/test_f_functions.py @@ -3,12 +3,8 @@ import harmonics_w5 as w5 import numpy as np -from ekore import harmonics +from ekore.harmonics import cache as c -zeta2 = harmonics.constants.zeta2 -zeta3 = harmonics.constants.zeta3 -zeta4 = harmonics.constants.zeta4 -zeta5 = harmonics.constants.zeta5 log2 = np.log(2) # reference values coming fom Mathematica: @@ -60,26 +56,29 @@ # F19, F20,F21 are not present in that paper. def test_F9(): for N, vals in zip(testN, refvals["S41"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) S41 = w5.S41(N, S1, S2, S3) np.testing.assert_allclose(S41, vals, atol=1e-05) def test_F11(): for N, vals in zip(testN, refvals["S311"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) S311 = w5.S311(N, S1, S2) np.testing.assert_allclose(S311, vals, atol=1e-05) def test_F13(): for N, vals in zip(testN, refvals["S221"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S21 = harmonics.S21(N, S1, S2) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S21 = c.get(c.S21, cache, N, None) S221 = w5.S221(N, S1, S2, S21) np.testing.assert_allclose(S221, vals, atol=1e-05) @@ -87,45 +86,49 @@ def test_F13(): def test_F12_F14(): # here there is a typo in eq (9.25) of Bl_mlein_2009 for N, vals in zip(testN, refvals["Sm221"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - Sm1 = harmonics.Sm1(N, S1) - S21 = harmonics.S21(N, S1, S2) - Sm21 = harmonics.Sm21(N, S1, Sm1) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + Sm1 = c.get(c.Sm1, cache, N, None) + S21 = c.get(c.S21, cache, N, None) + Sm21 = c.get(c.Sm21, cache, N, None) Sm221 = w5.Sm221(N, S1, Sm1, S21, Sm21) np.testing.assert_allclose(Sm221, vals, atol=1e-05) def test_F16(): for N, vals in zip(testN, refvals["S21m2"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) - Sm1 = harmonics.Sm1(N, S1) - Sm2 = harmonics.Sm2(N, S2) - Sm3 = harmonics.Sm3(N, S3) - S21 = harmonics.S21(N, S1, S2) - S2m1 = harmonics.S2m1(N, S2, Sm1, Sm2) - Sm21 = harmonics.Sm21(N, S1, Sm1) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) + Sm1 = c.get(c.Sm1, cache, N, None) + Sm2 = c.get(c.Sm2, cache, N, None) + Sm3 = c.get(c.Sm3, cache, N, None) + S21 = c.get(c.S21, cache, N, None) + S2m1 = c.get(c.S2m1, cache, N, None) + Sm21 = c.get(c.Sm21, cache, N, None) S21m2 = w5.S21m2(N, S1, S2, Sm1, Sm2, Sm3, S21, Sm21, S2m1) np.testing.assert_allclose(S21m2, vals, atol=1e-04) def test_F17(): for N, vals in zip(testN, refvals["S2111"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) S2111 = w5.S2111(N, S1, S2, S3) np.testing.assert_allclose(S2111, vals, atol=1e-05) def test_F18(): for N, vals in zip(testN, refvals["Sm2111"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) - Sm1 = harmonics.Sm1(N, S1) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) + Sm1 = c.get(c.Sm1, cache, N, None) Sm2111 = w5.Sm2111(N, S1, S2, S3, Sm1) np.testing.assert_allclose(Sm2111, vals, atol=1e-05) @@ -133,9 +136,10 @@ def test_F18(): # different parametrization, less accurate def test_F19(): for N, vals in zip(testN, refvals["S23"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) S23 = w5.S23(N, S1, S2, S3) np.testing.assert_allclose(S23, vals, rtol=2e-03) @@ -143,12 +147,13 @@ def test_F19(): # different parametrization, less accurate def test_F20(): for N, vals in zip(testN, refvals["Sm23"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) - Sm3 = harmonics.Sm3(N, S3) - Sm2 = harmonics.Sm2(N, S2) - Sm1 = harmonics.Sm1(N, S1) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) + Sm3 = c.get(c.Sm3, cache, N, None) + Sm2 = c.get(c.Sm2, cache, N, None) + Sm1 = c.get(c.Sm1, cache, N, None) Sm23 = w5.Sm23(N, Sm1, Sm2, Sm3) np.testing.assert_allclose(Sm23, vals, rtol=1e-03) @@ -156,11 +161,12 @@ def test_F20(): # different parametrization, less accurate def test_F21(): for N, vals in zip(testN, refvals["S2m3"]): - S1 = harmonics.S1(N) - S2 = harmonics.S2(N) - S3 = harmonics.S3(N) - Sm3 = harmonics.Sm3(N, S3) - Sm2 = harmonics.Sm2(N, S2) - Sm1 = harmonics.Sm1(N, S1) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) + Sm3 = c.get(c.Sm3, cache, N, None) + Sm2 = c.get(c.Sm2, cache, N, None) + Sm1 = c.get(c.Sm1, cache, N, None) S2m3 = w5.S2m3(N, S2, Sm1, Sm2, Sm3) np.testing.assert_allclose(S2m3, vals, rtol=1e-03) diff --git a/src/eko/beta.py b/src/eko/beta.py index 92212bf58..72c8cc272 100644 --- a/src/eko/beta.py +++ b/src/eko/beta.py @@ -6,8 +6,9 @@ import numba as nb +from eko.constants import zeta3 + from . import constants -from ekore.harmonics.constants import zeta3 @nb.njit(cache=True) diff --git a/src/eko/constants.py b/src/eko/constants.py index b5cfb54f7..3eb57aa1b 100644 --- a/src/eko/constants.py +++ b/src/eko/constants.py @@ -1,6 +1,8 @@ """This files sets the physical constants.""" import numba as nb +import numpy as np +from scipy.special import zeta NC = 3 """The number of colors.""" @@ -29,6 +31,23 @@ ed2 = 1.0 / 9 """Down quarks charge squared.""" +zeta2 = zeta(2) +r""":math:`\zeta(2)`""" + +zeta3 = zeta(3) +r""":math:`\zeta(3)`""" + +zeta4 = zeta(4) +r""":math:`\zeta(4)`""" + +zeta5 = zeta(5) +r""":math:`\zeta(5)`""" + +log2 = np.log(2) +r""":math:`\ln(2)`""" + +li4half = 0.517479 +""":math:`Li_{4}(1/2)`""" def update_colors(nc): """Updates the number of colors to :math:`NC = nc`. diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py index 3f241519a..ef4df8a62 100644 --- a/src/eko/evolution_operator/operator_matrix_element.py +++ b/src/eko/evolution_operator/operator_matrix_element.py @@ -21,22 +21,25 @@ def build_ome(A, matching_order, a_s, backward_method): r"""Construct the matching expansion in :math:`a_s` with the appropriate method. - Parameters - ---------- - A : numpy.ndarray - list of |OME| - matching_order : tuple(int,int) - perturbation matching order - a_s : float - strong coupling, needed only for the exact inverse - backward_method : ["exact", "expanded" or ""] - empty or method for inverting the matching condition (exact or expanded) - - Returns - ------- - ome : numpy.ndarray - matching operator matrix + Parameters + ---------- + A : numpy.ndarray + list of |OME| + matching_order : tuple(int,int) + perturbation matching order + a_s : float + strong coupling, needed only for the exact inverse + backward_method : ["exact", "expanded" or ""] + empty or method for inverting the matching condition (exact or expanded) + Returns + ------- + ome : numpy.ndarray + matching operator matrix + <<<<<<< HEAD:src/eko/evolution_operator/operator_matrix_element.py + ======= + + >>>>>>> master:src/eko/matching_conditions/operator_matrix_element.py """ # to get the inverse one can use this FORM snippet # Symbol a; @@ -88,42 +91,45 @@ def quad_ker( ): r"""Raw kernel inside quad. - Parameters - ---------- - u : float - quad argument - order : tuple(int,int) - perturbation matching order - mode0 : int - pid for first element in the singlet sector - mode1 : int - pid for second element in the singlet sector - is_log : boolean - logarithmic interpolation - logx : float - Mellin inversion point - areas : tuple - basis function configuration - a_s : float - strong coupling, needed only for the exact inverse - nf: int - number of active flavor below threshold - L : float - :math:``\ln(\mu_F^2 / m_h^2)`` - backward_method : ["exact", "expanded" or ""] - empty or method for inverting the matching condition (exact or expanded) - is_msbar: bool - add the |MSbar| contribution - is_polarized : boolean - is polarized evolution ? - is_time_like : boolean - is time-like evolution ? + Parameters + ---------- + u : float + quad argument + order : tuple(int,int) + perturbation matching order + mode0 : int + pid for first element in the singlet sector + mode1 : int + pid for second element in the singlet sector + is_log : boolean + logarithmic interpolation + logx : float + Mellin inversion point + areas : tuple + basis function configuration + a_s : float + strong coupling, needed only for the exact inverse + nf: int + number of active flavor below threshold + L : float + :math:``\ln(\mu_F^2 / m_h^2)`` + backward_method : ["exact", "expanded" or ""] + empty or method for inverting the matching condition (exact or expanded) + is_msbar: bool + add the |MSbar| contribution + is_polarized : boolean + is polarized evolution ? + is_time_like : boolean + is time-like evolution ? - Returns - ------- - ker : float - evaluated integration kernel + Returns + ------- + ker : float + evaluated integration kernel + <<<<<<< HEAD:src/eko/evolution_operator/operator_matrix_element.py + ======= + >>>>>>> master:src/eko/matching_conditions/operator_matrix_element.py """ ker_base = QuadKerBase(u, is_log, logx, mode0) integrand = ker_base.integrand(areas) @@ -131,27 +137,7 @@ def quad_ker( return 0.0 max_weight_dict = {1: 2, 2: 3, 3: 5} - sx = harmonics.compute_cache( - ker_base.n, max_weight_dict[order[0]], ker_base.is_singlet - ) - sx_ns = sx.copy() - if order[0] == 3 and ( - (backward_method != "" and ker_base.is_singlet) - or (mode0 == 100 and mode1 == 100) - ): - # At N3LO for A_qq singlet or backward you need to compute - # both the singlet and non-singlet like harmonics - # avoiding recomputing all of them ... - smx_ns = harmonics.smx(ker_base.n, np.array([s[0] for s in sx]), False) - for w, sm in enumerate(smx_ns): - sx_ns[w][-1] = sm - sx_ns[2][2] = harmonics.S2m1(ker_base.n, sx[0][1], smx_ns[0], smx_ns[1], False) - sx_ns[2][3] = harmonics.Sm21(ker_base.n, sx[0][0], smx_ns[0], False) - sx_ns[3][5] = harmonics.Sm31(ker_base.n, sx[0][0], smx_ns[0], smx_ns[1], False) - sx_ns[3][4] = harmonics.Sm211(ker_base.n, sx[0][0], sx[0][1], smx_ns[0], False) - sx_ns[3][3] = harmonics.Sm22( - ker_base.n, sx[0][0], sx[0][1], smx_ns[1], sx_ns[3][5], False - ) + # compute the ome if ker_base.is_singlet: @@ -160,24 +146,24 @@ def quad_ker( if is_time_like: raise NotImplementedError("Polarized, time-like is not implemented") else: - A = ome_ps.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns) + A = ome_ps.A_singlet(order, ker_base.n, nf, L, is_msbar) else: if is_time_like: - A = ome_ut.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns) + A = ome_ut.A_singlet(order, ker_base.n, nf, L, is_msbar) else: - A = ome_us.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns) + A = ome_us.A_singlet(order, ker_base.n, nf, L, is_msbar) else: indices = {200: 0, 91: 1} if is_polarized: if is_time_like: raise NotImplementedError("Polarized, time-like is not implemented") else: - A = ome_ps.A_non_singlet(order, ker_base.n, sx, nf, L) + A = ome_ps.A_non_singlet(order, ker_base.n, nf, L) else: if is_time_like: - A = ome_ut.A_non_singlet(order, ker_base.n, sx, nf, L) + A = ome_ut.A_non_singlet(order, ker_base.n, nf, L) else: - A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L) + A = ome_us.A_non_singlet(order, ker_base.n, nf, L) # build the expansion in alpha_s depending on the strategy ker = build_ome(A, order, a_s, backward_method) diff --git a/src/eko/gamma.py b/src/eko/gamma.py index b2c4a857b..6e67d2845 100644 --- a/src/eko/gamma.py +++ b/src/eko/gamma.py @@ -5,7 +5,7 @@ """ import numba as nb -from ekore.harmonics.constants import zeta3, zeta4, zeta5 +from eko.constants import zeta3, zeta4, zeta5 @nb.njit(cache=True) diff --git a/src/ekore/anomalous_dimensions/polarized/time_like/__init__.py b/src/ekore/anomalous_dimensions/polarized/time_like/__init__.py new file mode 100644 index 000000000..74aa85770 --- /dev/null +++ b/src/ekore/anomalous_dimensions/polarized/time_like/__init__.py @@ -0,0 +1,21 @@ +r"""The polarized, time-like Altarelli-Parisi splitting kernels. + +Normalization is given by + +.. math:: + \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x) + +with :math:`a_s = \frac{\alpha_S(\mu^2)}{4\pi}`. +""" + +import numba as nb + + +@nb.njit(cache=True) +def gamma_ns(_order, _mode, _n, _nf): + raise NotImplementedError("Polarised, time-like is not yet implemented") + + +@nb.njit(cache=True) +def gamma_singlet(_order, _n, _nf): + raise NotImplementedError("Polarised, time-like is not yet implemented") diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index dc4ad9dbf..7e8bd0bb3 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py @@ -19,6 +19,7 @@ import numpy as np from .... import harmonics +from ....harmonics import cache as c from . import aem1, aem2, as1, as2, as3, as4 @@ -42,62 +43,39 @@ def gamma_ns(order, mode, n, nf): numpy.ndarray non-singlet anomalous dimensions - See Also - -------- - ekore.anomalous_dimensions.unpolarized.space_like.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_nsp : :math:`\gamma_{ns,+}^{(3)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_nsm : :math:`\gamma_{ns,-}^{(3)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_nsv : :math:`\gamma_{ns,v}^{(3)}(N)` - """ - # cache the s-es - if order[0] >= 4: - full_sx_cache = harmonics.compute_cache(n, 5, is_singlet=False) - sx = np.array( - [ - full_sx_cache[0][0], - full_sx_cache[1][0], - full_sx_cache[2][0], - full_sx_cache[3][0], - ] - ) - else: - sx = harmonics.sx(n, max_weight=order[0] + 1) - # now combine + gamma_ns = np.zeros(order[0], np.complex_) - gamma_ns[0] = as1.gamma_ns(n, sx[0]) + + cache = c.reset() + gamma_ns[0] = as1.gamma_ns(n, cache) # NLO and beyond if order[0] >= 2: if mode == 10101: - gamma_ns_1 = as2.gamma_nsp(n, nf, sx) + gamma_ns_1 = as2.gamma_nsp(n, nf, cache, False) # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here elif mode in [10201, 10200]: - gamma_ns_1 = as2.gamma_nsm(n, nf, sx) + gamma_ns_1 = as2.gamma_nsm(n, nf, cache, False) else: raise NotImplementedError("Non-singlet sector is not implemented") gamma_ns[1] = gamma_ns_1 # NNLO and beyond if order[0] >= 3: if mode == 10101: - gamma_ns_2 = as3.gamma_nsp(n, nf, sx) + gamma_ns_2 = as3.gamma_nsp(n, nf, cache, False) elif mode == 10201: - gamma_ns_2 = as3.gamma_nsm(n, nf, sx) + gamma_ns_2 = as3.gamma_nsm(n, nf, cache, False) elif mode == 10200: - gamma_ns_2 = as3.gamma_nsv(n, nf, sx) + gamma_ns_2 = as3.gamma_nsv(n, nf, cache, False) gamma_ns[2] = gamma_ns_2 # N3LO if order[0] >= 4: if mode == 10101: - gamma_ns_3 = as4.gamma_nsp(n, nf, full_sx_cache) + gamma_ns_3 = as4.gamma_nsp(n, nf, cache, False) elif mode == 10201: - gamma_ns_3 = as4.gamma_nsm(n, nf, full_sx_cache) + gamma_ns_3 = as4.gamma_nsm(n, nf, cache, False) elif mode == 10200: - gamma_ns_3 = as4.gamma_nsv(n, nf, full_sx_cache) + gamma_ns_3 = as4.gamma_nsv(n, nf, cache, False) gamma_ns[3] = gamma_ns_3 return gamma_ns @@ -120,37 +98,15 @@ def gamma_singlet(order, n, nf): numpy.ndarray singlet anomalous dimensions matrices - See Also - -------- - ekore.anomalous_dimensions.unpolarized.space_like.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)` - ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_singlet : :math:`\gamma_{S}^{(3)}(N)` - """ - # cache the s-es - if order[0] >= 4: - full_sx_cache = harmonics.compute_cache(n, 5, is_singlet=False) - sx = np.array( - [ - full_sx_cache[0][0], - full_sx_cache[1][0], - full_sx_cache[2][0], - full_sx_cache[3][0], - ] - ) - elif order[0] >= 3: - # here we need only S1,S2,S3,S4 - sx = harmonics.sx(n, max_weight=order[0] + 1) - else: - sx = harmonics.sx(n, max_weight=order[0]) + cache = c.reset() gamma_s = np.zeros((order[0], 2, 2), np.complex_) - gamma_s[0] = as1.gamma_singlet(n, sx[0], nf) + gamma_s[0] = as1.gamma_singlet(n, nf, cache) if order[0] >= 2: - gamma_s[1] = as2.gamma_singlet(n, nf, sx) + gamma_s[1] = as2.gamma_singlet(n, nf, cache, True) if order[0] >= 3: - gamma_s[2] = as3.gamma_singlet(n, nf, sx) + gamma_s[2] = as3.gamma_singlet(n, nf, cache, True) if order[0] >= 4: - gamma_s[3] = as4.gamma_singlet(n, nf, full_sx_cache) + gamma_s[3] = as4.gamma_singlet(n, nf, cache, True) return gamma_s diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py index d26683969..79c2a0163 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py @@ -5,6 +5,7 @@ import numba as nb from eko import constants + from . import as1 @@ -75,7 +76,7 @@ def gamma_phph(nf): @nb.njit(cache=True) -def gamma_ns(N, s1): +def gamma_ns(N, cache): """ Computes the leading-order non-singlet QED anomalous dimension. @@ -85,12 +86,12 @@ def gamma_ns(N, s1): ---------- N : complex Mellin moment - s1 : complex - S1(N) + cache : numpy.ndarray + Harmonic sum cache Returns ------- gamma_ns : complex Leading-order non-singlet QED anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` """ - return as1.gamma_ns(N, s1) / constants.CF + return as1.gamma_ns(N, cache) / constants.CF diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py index 83f40e456..5e5613022 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py @@ -5,6 +5,8 @@ import numba as nb from eko import constants + +from ....harmonics import cache as c from . import as1aem1 @@ -35,7 +37,7 @@ def gamma_phph(N, nf): @nb.njit(cache=True) -def gamma_uph(N, nf, sx): +def gamma_uph(N, nf, cache, is_singlet): """Computes the O(aem2) quark-photon anomalous dimension for up quarks. Implements Eq. (55) of :cite:`deFlorian:2016gvk` for q=u. @@ -46,8 +48,10 @@ def gamma_uph(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -55,11 +59,11 @@ def gamma_uph(N, nf, sx): O(aem2) quark-photon anomalous dimension :math:`\\gamma_{u \\gamma}^{(0,2)}(N)` """ - return constants.eu2 * as1aem1.gamma_qph(N, nf, sx) / constants.CF + return constants.eu2 * as1aem1.gamma_qph(N, nf, cache, is_singlet) / constants.CF @nb.njit(cache=True) -def gamma_dph(N, nf, sx): +def gamma_dph(N, nf, cache, is_singlet): """Computes the O(aem2) quark-photon anomalous dimension for down quarks. Implements Eq. (55) of :cite:`deFlorian:2016gvk` for q=d. @@ -70,8 +74,10 @@ def gamma_dph(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -79,11 +85,11 @@ def gamma_dph(N, nf, sx): O(aem2) quark-photon anomalous dimension :math:`\\gamma_{d \\gamma}^{(0,2)}(N)` """ - return constants.ed2 * as1aem1.gamma_qph(N, nf, sx) / constants.CF + return constants.ed2 * as1aem1.gamma_qph(N, nf, cache, is_singlet) / constants.CF @nb.njit(cache=True) -def gamma_phu(N, nf, sx): +def gamma_phu(N, nf, cache, is_singlet): """Computes the O(aem2) photon-quark anomalous dimension for up quarks. Implements Eq. (56) of :cite:`deFlorian:2016gvk` for q=u. @@ -94,8 +100,10 @@ def gamma_phu(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -105,18 +113,21 @@ def gamma_phu(N, nf, sx): """ nu = constants.uplike_flavors(nf) nd = nf - nu - S1 = sx[0] + S1 = c.get(c.S1, cache, N, is_singlet) tmp = (-16 * (-16 - 27 * N - 13 * N**2 - 8 * N**3)) / ( 9.0 * (-1 + N) * N * (1 + N) ** 2 ) - 16 * (2 + 3 * N + 2 * N**2 + N**3) / ( 3.0 * (-1 + N) * N * (1 + N) ** 2 ) * S1 eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) - return constants.eu2 * as1aem1.gamma_phq(N, sx) / constants.CF + eSigma2 * tmp + return ( + constants.eu2 * as1aem1.gamma_phq(N, cache, is_singlet) / constants.CF + + eSigma2 * tmp + ) @nb.njit(cache=True) -def gamma_phd(N, nf, sx): +def gamma_phd(N, nf, cache, is_singlet): """Computes the O(aem2) photon-quark anomalous dimension for down quarks. Implements Eq. (56) of :cite:`deFlorian:2016gvk` for q=d. @@ -127,8 +138,10 @@ def gamma_phd(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -138,18 +151,21 @@ def gamma_phd(N, nf, sx): """ nu = constants.uplike_flavors(nf) nd = nf - nu - S1 = sx[0] + S1 = c.get(c.S1, cache, N, is_singlet) tmp = (-16 * (-16 - 27 * N - 13 * N**2 - 8 * N**3)) / ( 9.0 * (-1 + N) * N * (1 + N) ** 2 ) - 16 * (2 + 3 * N + 2 * N**2 + N**3) / ( 3.0 * (-1 + N) * N * (1 + N) ** 2 ) * S1 eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) - return constants.ed2 * as1aem1.gamma_phq(N, sx) / constants.CF + eSigma2 * tmp + return ( + constants.ed2 * as1aem1.gamma_phq(N, cache, is_singlet) / constants.CF + + eSigma2 * tmp + ) @nb.njit(cache=True) -def gamma_nspu(N, nf, sx): +def gamma_nspu(N, nf, cache, is_singlet): """Computes the O(aem2) singlet-like non-singlet anomalous dimension for up quarks. Implements sum of Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=u. @@ -160,9 +176,10 @@ def gamma_nspu(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- gamma_nspu : complex @@ -170,8 +187,8 @@ def gamma_nspu(N, nf, sx): :math:`\\gamma_{ns,+,u}^{(0,2)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) nu = constants.uplike_flavors(nf) nd = nf - nu eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) @@ -182,11 +199,13 @@ def gamma_nspu(N, nf, sx): - 80 / 9 * S1 + 16 / 3 * S2 ) * eSigma2 - return constants.eu2 * as1aem1.gamma_nsp(N, sx) / constants.CF / 2 + tmp + return ( + constants.eu2 * as1aem1.gamma_nsp(N, cache, is_singlet) / constants.CF / 2 + tmp + ) @nb.njit(cache=True) -def gamma_nspd(N, nf, sx): +def gamma_nspd(N, nf, cache, is_singlet): """Computes the O(aem2) singlet-like non-singlet anomalous dimension for down quarks. Implements sum of Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=d. @@ -197,8 +216,10 @@ def gamma_nspd(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -207,8 +228,8 @@ def gamma_nspd(N, nf, sx): :math:`\\gamma_{ns,+,d}^{(0,2)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) nu = constants.uplike_flavors(nf) nd = nf - nu eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) @@ -219,11 +240,13 @@ def gamma_nspd(N, nf, sx): - 80 / 9 * S1 + 16 / 3 * S2 ) * eSigma2 - return constants.ed2 * as1aem1.gamma_nsp(N, sx) / constants.CF / 2 + tmp + return ( + constants.ed2 * as1aem1.gamma_nsp(N, cache, is_singlet) / constants.CF / 2 + tmp + ) @nb.njit(cache=True) -def gamma_nsmu(N, nf, sx): +def gamma_nsmu(N, nf, cache, is_singlet): """Computes the O(aem2) valence-like non-singlet anomalous dimension for up quarks. Implements difference between Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=u. @@ -234,8 +257,10 @@ def gamma_nsmu(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -244,8 +269,8 @@ def gamma_nsmu(N, nf, sx): :math:`\\gamma_{ns,-,u}^{(0,2)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) nu = constants.uplike_flavors(nf) nd = nf - nu eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) @@ -256,11 +281,13 @@ def gamma_nsmu(N, nf, sx): - 80 / 9 * S1 + 16 / 3 * S2 ) * eSigma2 - return constants.eu2 * as1aem1.gamma_nsm(N, sx) / constants.CF / 2 + tmp + return ( + constants.eu2 * as1aem1.gamma_nsm(N, cache, is_singlet) / constants.CF / 2 + tmp + ) @nb.njit(cache=True) -def gamma_nsmd(N, nf, sx): +def gamma_nsmd(N, nf, cache, is_singlet): """Computes the O(aem2) valence-like non-singlet anomalous dimension for down quarks. Implements difference between Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=d. @@ -271,8 +298,10 @@ def gamma_nsmd(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -281,8 +310,8 @@ def gamma_nsmd(N, nf, sx): :math:`\\gamma_{ns,-,d}^{(0,2)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) nu = constants.uplike_flavors(nf) nd = nf - nu eSigma2 = constants.NC * (nu * constants.eu2 + nd * constants.ed2) @@ -293,7 +322,9 @@ def gamma_nsmd(N, nf, sx): - 80 / 9 * S1 + 16 / 3 * S2 ) * eSigma2 - return constants.ed2 * as1aem1.gamma_nsm(N, sx) / constants.CF / 2 + tmp + return ( + constants.ed2 * as1aem1.gamma_nsm(N, cache, is_singlet) / constants.CF / 2 + tmp + ) @nb.njit(cache=True) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py index 28afb3a00..ad83127de 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py @@ -5,9 +5,11 @@ from eko import constants +from ....harmonics import cache as c + @nb.njit(cache=True) -def gamma_ns(N, s1): +def gamma_ns(N, cache): """ Computes the leading-order non-singlet anomalous dimension. @@ -17,15 +19,13 @@ def gamma_ns(N, s1): ---------- N : complex Mellin moment - s1 : complex - harmonic sum :math:`S_{1}` Returns ------- gamma_ns : complex Leading-order non-singlet anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)` """ - gamma = -(3.0 - 4.0 * s1 + 2.0 / N / (N + 1.0)) + gamma = -(3.0 - 4.0 * c.get(c.S1, cache, N, False) + 2.0 / N / (N + 1.0)) result = constants.CF * gamma return result @@ -77,7 +77,7 @@ def gamma_gq(N): @nb.njit(cache=True) -def gamma_gg(N, s1, nf): +def gamma_gg(N, nf, cache): """ Computes the leading-order gluon-gluon anomalous dimension @@ -87,8 +87,6 @@ def gamma_gg(N, s1, nf): ---------- N : complex Mellin moment - s1 : complex - harmonic sum :math:`S_{1}` nf : int Number of active flavors @@ -97,13 +95,15 @@ def gamma_gg(N, s1, nf): gamma_gg : complex Leading-order gluon-gluon anomalous dimension :math:`\\gamma_{gg}^{(0)}(N)` """ - gamma = s1 - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0) + gamma = ( + c.get(c.S1, cache, N, True) - 1.0 / N / (N - 1.0) - 1.0 / (N + 1.0) / (N + 2.0) + ) result = constants.CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * constants.TR * nf return result @nb.njit(cache=True) -def gamma_singlet(N, s1, nf): +def gamma_singlet(N, nf, cache): r""" Computes the leading-order singlet anomalous dimension matrix @@ -134,8 +134,9 @@ def gamma_singlet(N, s1, nf): gamma_gq : :math:`\gamma_{gq}^{(0)}` gamma_gg : :math:`\gamma_{gg}^{(0)}` """ - gamma_qq = gamma_ns(N, s1) + gamma_qq = gamma_ns(N, cache) gamma_S_0 = np.array( - [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, s1, nf)]], np.complex_ + [[gamma_qq, gamma_qg(N, nf)], [gamma_gq(N), gamma_gg(N, nf, cache)]], + np.complex_, ) return gamma_S_0 diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py index e16133a05..4c0db5b10 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py @@ -5,12 +5,13 @@ import numba as nb from eko import constants -from .... import harmonics -from ....harmonics.constants import zeta2, zeta3 +from eko.constants import zeta2, zeta3 + +from ....harmonics import cache as c @nb.njit(cache=True) -def gamma_phq(N, sx): +def gamma_phq(N, cache, is_singlet): """Computes the O(as1aem1) photon-quark anomalous dimension Implements Eq. (36) of :cite:`deFlorian:2015ujt`. @@ -19,8 +20,10 @@ def gamma_phq(N, sx): ---------- N : complex Mellin moment - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -28,8 +31,8 @@ def gamma_phq(N, sx): O(as1aem1) photon-quark anomalous dimension :math:`\\gamma_{\\gamma q}^{(1,1)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) tmp_const = ( 2 * (-4 - 12 * N - N**2 + 28 * N**3 + 43 * N**4 + 30 * N**5 + 12 * N**6) @@ -47,7 +50,7 @@ def gamma_phq(N, sx): @nb.njit(cache=True) -def gamma_qph(N, nf, sx): +def gamma_qph(N, nf, cache, is_singlet): """Computes the O(as1aem1) quark-photon anomalous dimension Implements Eq. (26) of :cite:`deFlorian:2015ujt`. @@ -58,8 +61,10 @@ def gamma_qph(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -67,8 +72,8 @@ def gamma_qph(N, nf, sx): O(as1aem1) quark-photon anomalous dimension :math:`\\gamma_{q \\gamma}^{(1,1)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) tmp_const = ( -2 * ( @@ -140,7 +145,7 @@ def gamma_phg(N): @nb.njit(cache=True) -def gamma_qg(N, nf, sx): +def gamma_qg(N, nf, cache, is_singlet): """Computes the O(as1aem1) quark-gluon singlet anomalous dimension. Implements Eq. (29) of :cite:`deFlorian:2015ujt`. @@ -151,8 +156,10 @@ def gamma_qg(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -161,11 +168,13 @@ def gamma_qg(N, nf, sx): :math:`\\gamma_{qg}^{(1,1)}(N)` """ - return constants.TR / constants.CF / constants.CA * gamma_qph(N, nf, sx) + return ( + constants.TR / constants.CF / constants.CA * gamma_qph(N, nf, cache, is_singlet) + ) @nb.njit(cache=True) -def gamma_gq(N, sx): +def gamma_gq(N, cache, is_singlet): """Computes the O(as1aem1) gluon-quark singlet anomalous dimension. Implements Eq. (35) of :cite:`deFlorian:2015ujt`. @@ -174,8 +183,10 @@ def gamma_gq(N, sx): ---------- N : complex Mellin moment - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -184,7 +195,7 @@ def gamma_gq(N, sx): :math:`\\gamma_{gq}^{(1,1)}(N)` """ - return gamma_phq(N, sx) + return gamma_phq(N, cache, is_singlet) @nb.njit(cache=True) @@ -230,7 +241,7 @@ def gamma_gg(): @nb.njit(cache=True) -def gamma_nsp(N, sx): +def gamma_nsp(N, cache, is_singlet): """Computes the O(as1aem1) singlet-like non-singlet anomalous dimension. Implements sum of Eqs. (33-34) of :cite:`deFlorian:2015ujt`. @@ -239,8 +250,10 @@ def gamma_nsp(N, sx): ---------- N : complex Mellin moment - sx : np array - List of harmonic sums + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -249,18 +262,18 @@ def gamma_nsp(N, sx): :math:`\\gamma_{ns,+}^{(1)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] - S1h = harmonics.S1(N / 2) - S2h = harmonics.S2(N / 2) - S3h = harmonics.S3(N / 2) - S1p1h = harmonics.S1((N + 1.0) / 2) - S2p1h = harmonics.S2((N + 1) / 2) - S3p1h = harmonics.S3((N + 1) / 2) - g3N = harmonics.g_functions.mellin_g3(N, S1) - S1p2 = harmonics.polygamma.recursive_harmonic_sum(S1, N, 2, 1) - g3Np2 = harmonics.g_functions.mellin_g3(N + 2, S1p2) + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) + S3 = c.get(c.S3, cache, N, is_singlet) + S1h = c.get(c.S1h, cache, N, is_singlet) + S2h = c.get(c.S2h, cache, N, is_singlet) + S3h = c.get(c.S3h, cache, N, is_singlet) + S1p1h = c.get(c.S1ph, cache, N, is_singlet) + S2p1h = c.get(c.S2ph, cache, N, is_singlet) + S3p1h = c.get(c.S3ph, cache, N, is_singlet) + g3N = c.get(c.g3, cache, N, is_singlet) + S1p2 = c.get(c.S1p2, cache, N, is_singlet) + g3Np2 = c.get(c.g3p2, cache, N, is_singlet) result = ( +32 * zeta2 * S1h - 32 * zeta2 * S1p1h @@ -293,7 +306,7 @@ def gamma_nsp(N, sx): @nb.njit(cache=True) -def gamma_nsm(N, sx): +def gamma_nsm(N, cache, is_singlet): """Computes the O(as1aem1) valence-like non-singlet anomalous dimension. Implements difference between Eqs. (33-34) of :cite:`deFlorian:2015ujt`. @@ -302,9 +315,10 @@ def gamma_nsm(N, sx): ---------- N : complex Mellin moment - sx : np array - List of harmonic sums - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- gamma_nsm : complex @@ -312,18 +326,18 @@ def gamma_nsm(N, sx): :math:`\\gamma_{ns,-}^{(1,1)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] - S1h = harmonics.S1(N / 2) - S2h = harmonics.S2(N / 2) - S3h = harmonics.S3(N / 2) - S1p1h = harmonics.S1((N + 1.0) / 2) - S2p1h = harmonics.S2((N + 1) / 2) - S3p1h = harmonics.S3((N + 1) / 2) - g3N = harmonics.g_functions.mellin_g3(N, S1) - S1p2 = harmonics.polygamma.recursive_harmonic_sum(S1, N, 2, 1) - g3Np2 = harmonics.g_functions.mellin_g3(N + 2, S1p2) + S1 = c.get(c.S1, cache, N, is_singlet) + S2 = c.get(c.S2, cache, N, is_singlet) + S3 = c.get(c.S3, cache, N, is_singlet) + S1h = c.get(c.S1h, cache, N, is_singlet) + S2h = c.get(c.S2h, cache, N, is_singlet) + S3h = c.get(c.S3h, cache, N, is_singlet) + S1p1h = c.get(c.S1ph, cache, N, is_singlet) + S2p1h = c.get(c.S2ph, cache, N, is_singlet) + S3p1h = c.get(c.S3ph, cache, N, is_singlet) + g3N = c.get(c.g3, cache, N, is_singlet) + S1p2 = c.get(c.S1p2, cache, N, is_singlet) + g3Np2 = c.get(c.g3p2, cache, N, is_singlet) result = ( -32.0 * zeta2 * S1h diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py index 8691a54b2..88adc7eb4 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py @@ -10,12 +10,13 @@ import numpy as np from eko import constants -from .... import harmonics -from ....harmonics.constants import log2, zeta2, zeta3 +from eko.constants import log2, zeta2, zeta3 + +from ....harmonics import cache as c @nb.njit(cache=True) -def gamma_nsm(n, nf, sx): +def gamma_nsm(n, nf, cache, is_singlet): """ Computes the |NLO| valence-like non-singlet anomalous dimension. @@ -27,8 +28,10 @@ def gamma_nsm(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : numpy.ndarray - List of harmonic sums: :math:`S_{1},S_{2}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -36,14 +39,14 @@ def gamma_nsm(n, nf, sx): |NLO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(1)}(N)` """ - S1 = sx[0] - S2 = sx[1] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) # Here, Sp refers to S' ("s-prime") (german: "s-strich" or in Pegasus language: SSTR) # of :cite:`Gluck:1989ze` and NOT to the Spence function a.k.a. dilogarithm - Sp1m = harmonics.S1((n - 1) / 2) - Sp2m = harmonics.S2((n - 1) / 2) - Sp3m = harmonics.S3((n - 1) / 2) - g3n = harmonics.g_functions.mellin_g3(n, S1) + Sp1m = c.get(c.S1mh, cache, n, is_singlet) + Sp2m = c.get(c.S2mh, cache, n, is_singlet) + Sp3m = c.get(c.S3mh, cache, n, is_singlet) + g3n = c.get(c.g3, cache, n, is_singlet) # fmt: off gqq1m_cfca = 16*g3n - (144 + n*(1 + n)*(156 + n*(340 + n*(655 + 51*n*(2 + n)))))/(18.*np.power(n,3)*np.power(1 + n,3)) + (-14.666666666666666 + 8/n - 8/(1 + n))*S2 - (4*Sp2m)/(n + np.power(n,2)) + S1*(29.77777777777778 + 16/np.power(n,2) - 16*S2 + 8*Sp2m) + 2*Sp3m + 10*zeta3 + zeta2*(16*S1 - 16*Sp1m - (16*(1 + n*log2))/n) # pylint: disable=line-too-long gqq1m_cfcf = -32*g3n + (24 - n*(-32 + 3*n*(-8 + n*(3 + n)*(3 + np.power(n,2)))))/(2.*np.power(n,3)*np.power(1 + n,3)) + (12 - 8/n + 8/(1 + n))*S2 + S1*(-24/np.power(n,2) - 8/np.power(1 + n,2) + 16*S2 - 16*Sp2m) + (8*Sp2m)/(n + np.power(n,2)) - 4*Sp3m - 20*zeta3 + zeta2*(-32*S1 + 32*Sp1m + 32*(1/n + log2)) # pylint: disable=line-too-long @@ -58,7 +61,7 @@ def gamma_nsm(n, nf, sx): @nb.njit(cache=True) -def gamma_nsp(n, nf, sx): +def gamma_nsp(n, nf, cache, is_singlet): """ Computes the |NLO| singlet-like non-singlet anomalous dimension. @@ -70,8 +73,10 @@ def gamma_nsp(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : numpy.ndarray - List of harmonic sums: :math:`S_{1},S_{2}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -79,12 +84,12 @@ def gamma_nsp(n, nf, sx): |NLO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(1)}(N)` """ - S1 = sx[0] - S2 = sx[1] - Sp1p = harmonics.S1(n / 2) - Sp2p = harmonics.S2(n / 2) - Sp3p = harmonics.S3(n / 2) - g3n = harmonics.g_functions.mellin_g3(n, S1) + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sp1p = c.get(c.S1h, cache, n, is_singlet) + Sp2p = c.get(c.S2h, cache, n, is_singlet) + Sp3p = c.get(c.S3h, cache, n, is_singlet) + g3n = c.get(c.g3, cache, n, is_singlet) # fmt: off gqq1p_cfca = -16*g3n + (132 - n*(340 + n*(655 + 51*n*(2 + n))))/(18.*np.power(n,2)*np.power(1 + n,2)) + (-14.666666666666666 + 8/n - 8/(1 + n))*S2 - (4*Sp2p)/(n + np.power(n,2)) + S1*(29.77777777777778 - 16/np.power(n,2) - 16*S2 + 8*Sp2p) + 2*Sp3p + 10*zeta3 + zeta2*(16*S1 - 16*Sp1p + 16*(1/n - log2)) # pylint: disable=line-too-long gqq1p_cfcf = 32*g3n - (8 + n*(32 + n*(40 + 3*n*(3 + n)*(3 + np.power(n,2)))))/(2.*np.power(n,3)*np.power(1 + n,3)) + (12 - 8/n + 8/(1 + n))*S2 + S1*(40/np.power(n,2) - 8/np.power(1 + n,2) + 16*S2 - 16*Sp2p) + (8*Sp2p)/(n + np.power(n,2)) - 4*Sp3p - 20*zeta3 + zeta2*(-32*S1 + 32*Sp1p + 32*(-(1/n) + log2)) # pylint: disable=line-too-long @@ -126,7 +131,7 @@ def gamma_ps(n, nf): @nb.njit(cache=True) -def gamma_qg(n, nf, sx): +def gamma_qg(n, nf, cache, is_singlet): """ Computes the |NLO| quark-gluon singlet anomalous dimension. @@ -138,8 +143,10 @@ def gamma_qg(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : numpy.ndarray - List of harmonic sums: :math:`S_{1},S_{2}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -147,9 +154,9 @@ def gamma_qg(n, nf, sx): |NLO| quark-gluon singlet anomalous dimension :math:`\\gamma_{qg}^{(1)}(N)` """ - S1 = sx[0] - S2 = sx[1] - Sp2p = harmonics.S2(n / 2) + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sp2p = c.get(c.S2h, cache, n, is_singlet) # fmt: off gqg1_nfca = (-4*(16 + n*(64 + n*(104 + n*(128 + n*(85 + n*(36 + n*(25 + n*(15 + n*(6 + n))))))))))/((-1 + n)*np.power(n,3)*np.power(1 + n,3)*np.power(2 + n,3)) - (16*(3 + 2*n)*S1)/np.power(2 + 3*n + np.power(n,2),2) + (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(2 + 3*n + np.power(n,2))) - (4*(2 + n + np.power(n,2))*S2)/(n*(2 + 3*n + np.power(n,2))) + (4*(2 + n + np.power(n,2))*Sp2p)/(n*(2 + 3*n + np.power(n,2))) # pylint: disable=line-too-long gqg1_nfcf = (-2*(4 + n*(8 + n*(1 + n)*(25 + n*(26 + 5*n*(2 + n))))))/(np.power(n,3)*np.power(1 + n,3)*(2 + n)) + (8*S1)/np.power(n,2) - (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(2 + 3*n + np.power(n,2))) + (4*(2 + n + np.power(n,2))*S2)/(n*(2 + 3*n + np.power(n,2))) # pylint: disable=line-too-long @@ -161,7 +168,7 @@ def gamma_qg(n, nf, sx): @nb.njit(cache=True) -def gamma_gq(n, nf, sx): +def gamma_gq(n, nf, cache, is_singlet): """ Computes the |NLO| gluon-quark singlet anomalous dimension. @@ -173,8 +180,10 @@ def gamma_gq(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : numpy.ndarray - List of harmonic sums: :math:`S_{1},S_{2}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -182,9 +191,9 @@ def gamma_gq(n, nf, sx): |NLO| gluon-quark singlet anomalous dimension :math:`\\gamma_{gq}^{(1)}(N)` """ - S1 = sx[0] - S2 = sx[1] - Sp2p = harmonics.S2(n / 2) + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sp2p = c.get(c.S2h, cache, n, is_singlet) # fmt: off ggq1_cfcf = (-8 + 2*n*(-12 + n*(-1 + n*(28 + n*(43 + 6*n*(5 + 2*n))))))/((-1 + n)*np.power(n,3)*np.power(1 + n,3)) - (4*(10 + n*(17 + n*(8 + 5*n)))*S1)/((-1 + n)*n*np.power(1 + n,2)) + (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(-1 + np.power(n,2))) + (4*(2 + n + np.power(n,2))*S2)/(n*(-1 + np.power(n,2))) # pylint: disable=line-too-long ggq1_cfca = (-4*(144 + n*(432 + n*(-152 + n*(-1304 + n*(-1031 + n*(695 + n*(1678 + n*(1400 + n*(621 + 109*n))))))))))/(9.*np.power(n,3)*np.power(1 + n,3)*np.power(-2 + n + np.power(n,2),2)) + (4*(-12 + n*(-22 + 41*n + 17*np.power(n,3)))*S1)/(3.*np.power(-1 + n,2)*np.power(n,2)*(1 + n)) + ((8 + 4*n + 4*np.power(n,2))*np.power(S1,2))/(n - np.power(n,3)) + ((8 + 4*n + 4*np.power(n,2))*S2)/(n - np.power(n,3)) + (4*(2 + n + np.power(n,2))*Sp2p)/(n*(-1 + np.power(n,2))) # pylint: disable=line-too-long @@ -199,7 +208,7 @@ def gamma_gq(n, nf, sx): @nb.njit(cache=True) -def gamma_gg(n, nf, sx): +def gamma_gg(n, nf, cache, is_singlet): """ Computes the |NLO| gluon-gluon singlet anomalous dimension. @@ -211,8 +220,10 @@ def gamma_gg(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : numpy.ndarray - List of harmonic sums: :math:`S_{1},S_{2}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -220,11 +231,11 @@ def gamma_gg(n, nf, sx): |NLO| gluon-gluon singlet anomalous dimension :math:`\\gamma_{gg}^{(1)}(N)` """ - S1 = sx[0] - Sp1p = harmonics.S1(n / 2) - Sp2p = harmonics.S2(n / 2) - Sp3p = harmonics.S3(n / 2) - g3n = harmonics.g_functions.mellin_g3(n, S1) + S1 = c.get(c.S1, cache, n, is_singlet) + Sp1p = c.get(c.S1h, cache, n, is_singlet) + Sp2p = c.get(c.S2h, cache, n, is_singlet) + Sp3p = c.get(c.S3h, cache, n, is_singlet) + g3n = c.get(c.g3, cache, n, is_singlet) # fmt: off ggg1_caca = 16*g3n - (2*(576 + n*(1488 + n*(560 + n*(-1248 + n*(-1384 + n*(1663 + n*(4514 + n*(4744 + n*(3030 + n*(1225 + 48*n*(7 + n))))))))))))/(9.*np.power(-1 + n,2)*np.power(n,3)*np.power(1 + n,3)*np.power(2 + n,3)) + S1*(29.77777777777778 + 16/np.power(-1 + n,2) + 16/np.power(1 + n,2) - 16/np.power(2 + n,2) - 8*Sp2p) + (16*(1 + n + np.power(n,2))*Sp2p)/(n*(1 + n)*(-2 + n + np.power(n,2))) - 2*Sp3p - 10*zeta3 + zeta2*(-16*S1 + 16*Sp1p + 16*(-(1/n) + log2)) # pylint: disable=line-too-long ggg1_canf = (8*(6 + n*(1 + n)*(28 + n*(1 + n)*(13 + 3*n*(1 + n)))))/(9.*np.power(n,2)*np.power(1 + n,2)*(-2 + n + np.power(n,2))) - (40*S1)/9. # pylint: disable=line-too-long @@ -237,7 +248,7 @@ def gamma_gg(n, nf, sx): @nb.njit(cache=True) -def gamma_singlet(n, nf, sx): +def gamma_singlet(n, nf, cache, is_singlet): r""" Computes the next-leading-order singlet anomalous dimension matrix @@ -251,10 +262,12 @@ def gamma_singlet(n, nf, sx): ---------- N : complex Mellin moment - sx : numpy.ndarray - List of harmonic sums: :math:`S_{1},S_{2}` nf : int Number of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -269,9 +282,12 @@ def gamma_singlet(n, nf, sx): gamma_gq : :math:`\gamma_{gq}^{(1)}` gamma_gg : :math:`\gamma_{gg}^{(1)}` """ - gamma_qq = gamma_nsp(n, nf, sx) + gamma_ps(n, nf) + gamma_qq = gamma_nsp(n, nf, cache, is_singlet) + gamma_ps(n, nf) gamma_S_0 = np.array( - [[gamma_qq, gamma_qg(n, nf, sx)], [gamma_gq(n, nf, sx), gamma_gg(n, nf, sx)]], + [ + [gamma_qq, gamma_qg(n, nf, cache, is_singlet)], + [gamma_gq(n, nf, cache, is_singlet), gamma_gg(n, nf, cache, is_singlet)], + ], np.complex_, ) return gamma_S_0 diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py index a2387bd08..e0ef10e35 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py @@ -8,11 +8,13 @@ import numba as nb import numpy as np -from ....harmonics.constants import zeta2, zeta3 +from eko.constants import zeta2, zeta3 + +from ....harmonics import cache as c @nb.njit(cache=True) -def gamma_nsm(n, nf, sx): +def gamma_nsm(n, nf, cache, is_singlet): """ Computes the |NNLO| valence-like non-singlet anomalous dimension. @@ -24,8 +26,10 @@ def gamma_nsm(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -33,9 +37,9 @@ def gamma_nsm(n, nf, sx): |NNLO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E2 = 2.0 * (-S1 / n**3 + (zeta2 - S2) / n**2 - (S3 - zeta3) / n) @@ -89,7 +93,7 @@ def gamma_nsm(n, nf, sx): @nb.njit(cache=True) -def gamma_nsp(n, nf, sx): +def gamma_nsp(n, nf, cache, is_singlet): """ Computes the |NNLO| singlet-like non-singlet anomalous dimension. @@ -101,8 +105,10 @@ def gamma_nsp(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -110,9 +116,9 @@ def gamma_nsp(n, nf, sx): |NNLO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E2 = 2.0 * (-S1 / n**3 + (zeta2 - S2) / n**2 - (S3 - zeta3) / n) @@ -166,7 +172,7 @@ def gamma_nsp(n, nf, sx): @nb.njit(cache=True) -def gamma_nsv(n, nf, sx): +def gamma_nsv(n, nf, cache, is_singlet): """ Computes the |NNLO| valence non-singlet anomalous dimension. @@ -178,8 +184,10 @@ def gamma_nsv(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -187,9 +195,9 @@ def gamma_nsv(n, nf, sx): |NNLO| valence non-singlet anomalous dimension :math:`\\gamma_{ns,v}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E2 = 2.0 * (-S1 / n**3 + (zeta2 - S2) / n**2 - (S3 - zeta3) / n) @@ -216,12 +224,12 @@ def gamma_nsv(n, nf, sx): + 46.18 * E2 ) - result = gamma_nsm(n, nf, sx) + nf * ps2 + result = gamma_nsm(n, nf, cache, is_singlet) + nf * ps2 return result @nb.njit(cache=True) -def gamma_ps(n, nf, sx): +def gamma_ps(n, nf, cache, is_singlet): """ Computes the |NNLO| pure-singlet quark-quark anomalous dimension. @@ -233,8 +241,10 @@ def gamma_ps(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -242,9 +252,9 @@ def gamma_ps(n, nf, sx): |NNLO| pure-singlet quark-quark anomalous dimension :math:`\\gamma_{ps}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E11 = (S1 + 1.0 / (n + 1.0)) / (n + 1.0) ** 2 + ( @@ -293,7 +303,7 @@ def gamma_ps(n, nf, sx): @nb.njit(cache=True) -def gamma_qg(n, nf, sx): +def gamma_qg(n, nf, cache, is_singlet): """ Computes the |NNLO| quark-gluon singlet anomalous dimension. @@ -305,8 +315,10 @@ def gamma_qg(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -314,10 +326,10 @@ def gamma_qg(n, nf, sx): |NNLO| quark-gluon singlet anomalous dimension :math:`\\gamma_{qg}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] - S4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E2 = 2.0 * (-S1 / n**3 + (zeta2 - S2) / n**2 - (S3 - zeta3) / n) @@ -367,7 +379,7 @@ def gamma_qg(n, nf, sx): @nb.njit(cache=True) -def gamma_gq(n, nf, sx): +def gamma_gq(n, nf, cache, is_singlet): """ Computes the |NNLO| gluon-quark singlet anomalous dimension. @@ -379,8 +391,10 @@ def gamma_gq(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -388,10 +402,10 @@ def gamma_gq(n, nf, sx): |NNLO| gluon-quark singlet anomalous dimension :math:`\\gamma_{gq}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] - S4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E2 = 2.0 * (-S1 / n**3 + (zeta2 - S2) / n**2 - (S3 - zeta3) / n) @@ -457,7 +471,7 @@ def gamma_gq(n, nf, sx): @nb.njit(cache=True) -def gamma_gg(n, nf, sx): +def gamma_gg(n, nf, cache, is_singlet): """ Computes the |NNLO| gluon-gluon singlet anomalous dimension. @@ -469,8 +483,10 @@ def gamma_gg(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -478,9 +494,9 @@ def gamma_gg(n, nf, sx): |NNLO| gluon-gluon singlet anomalous dimension :math:`\\gamma_{gg}^{(2)}(N)` """ - S1 = sx[0] - S2 = sx[1] - S3 = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) E1 = S1 / n**2 + (S2 - zeta2) / n E11 = (S1 + 1.0 / (n + 1.0)) / (n + 1.0) ** 2 + ( @@ -545,7 +561,7 @@ def gamma_gg(n, nf, sx): @nb.njit(cache=True) -def gamma_singlet(N, nf, sx): +def gamma_singlet(N, nf, cache, is_singlet): r""" Computes the |NNLO| singlet anomalous dimension matrix @@ -561,8 +577,10 @@ def gamma_singlet(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : np.ndarray - List of harmonic sums: :math:`S_{1},S_{2},S_{3}` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns @@ -579,9 +597,12 @@ def gamma_singlet(N, nf, sx): gamma_gq : :math:`\gamma_{gq}^{(2)}` gamma_gg : :math:`\gamma_{gg}^{(2)}` """ - gamma_qq = gamma_nsp(N, nf, sx) + gamma_ps(N, nf, sx) + gamma_qq = gamma_nsp(N, nf, cache, is_singlet) + gamma_ps(N, nf, cache, is_singlet) gamma_S_0 = np.array( - [[gamma_qq, gamma_qg(N, nf, sx)], [gamma_gq(N, nf, sx), gamma_gg(N, nf, sx)]], + [ + [gamma_qq, gamma_qg(N, nf, cache, is_singlet)], + [gamma_gq(N, nf, cache, is_singlet), gamma_gg(N, nf, cache, is_singlet)], + ], np.complex_, ) return gamma_S_0 diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py index e6860191e..fa67a1486 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py @@ -16,7 +16,7 @@ @nb.njit(cache=True) -def gamma_singlet(N, nf, sx): +def gamma_singlet(N, nf, cache, is_singlet): r"""Computes the |N3LO| singlet anomalous dimension matrix .. math:: @@ -31,9 +31,10 @@ def gamma_singlet(N, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- numpy.ndarray @@ -49,9 +50,12 @@ def gamma_singlet(N, nf, sx): gamma_gg : :math:`\gamma_{gg}^{(3)}` """ - gamma_qq = gamma_nsp(N, nf, sx) + gamma_ps(N, nf, sx) + gamma_qq = gamma_nsp(N, nf, cache, is_singlet) + gamma_ps(N, nf, cache, is_singlet) gamma_S_0 = np.array( - [[gamma_qq, gamma_qg(N, nf, sx)], [gamma_gq(N, nf, sx), gamma_gg(N, nf, sx)]], + [ + [gamma_qq, gamma_qg(N, nf, cache, is_singlet)], + [gamma_gq(N, nf, cache, is_singlet), gamma_gg(N, nf, cache, is_singlet)], + ], np.complex_, ) return gamma_S_0 diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py index c627b474f..66d5b6159 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py @@ -3,11 +3,12 @@ import numba as nb import numpy as np +from .....harmonics import cache as c from .....harmonics.log_functions import lm11 @nb.njit(cache=True) -def gamma_gg_nf3(n, sx): +def gamma_gg_nf3(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^3` of :math:`\gamma_{gg}^{(3)}`, the expression is copied exact from Eq. 3.14 of :cite:`Davies:2016jie`. @@ -15,18 +16,20 @@ def gamma_gg_nf3(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- complex |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^3}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3, S21, _, _, _, _ = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) return 3.0 * ( -0.0205761316872428 + 2.599239604033225 / (-1.0 + n) @@ -137,23 +140,24 @@ def gamma_gg_nf3(n, sx): @nb.njit(cache=True) -def gamma_gg_nf1(n, sx): +def gamma_gg_nf1(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{gg}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- complex |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) return ( 18371.82290926215 + 1992.766087237516 / np.power(-1.0 + n, 3) @@ -168,15 +172,17 @@ def gamma_gg_nf1(n, sx): @nb.njit(cache=True) -def gamma_gg_nf2(n, sx): +def gamma_gg_nf2(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{gg}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -184,7 +190,7 @@ def gamma_gg_nf2(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^2}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) return ( -436.18166675733164 + 18.346203400819753 / np.power(-1.0 + n, 2) @@ -198,15 +204,17 @@ def gamma_gg_nf2(n, sx): @nb.njit(cache=True) -def gamma_gg_nf0(n, sx): +def gamma_gg_nf0(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^0` of :math:`\gamma_{gg}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -214,7 +222,7 @@ def gamma_gg_nf0(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^0}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) return ( -61782.868048046716 - 49851.703887834694 / np.power(-1.0 + n, 4) @@ -230,7 +238,7 @@ def gamma_gg_nf0(n, sx): @nb.njit(cache=True) -def gamma_gg(n, nf, sx): +def gamma_gg(n, nf, cache, is_singlet): r"""Computes the |N3LO| gluon-gluon singlet anomalous dimension. Parameters @@ -239,8 +247,10 @@ def gamma_gg(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -257,8 +267,8 @@ def gamma_gg(n, nf, sx): """ return ( - gamma_gg_nf0(n, sx) - + nf * gamma_gg_nf1(n, sx) - + nf**2 * gamma_gg_nf2(n, sx) - + nf**3 * gamma_gg_nf3(n, sx) + gamma_gg_nf0(n, cache, is_singlet) + + nf * gamma_gg_nf1(n, cache, is_singlet) + + nf**2 * gamma_gg_nf2(n, cache, is_singlet) + + nf**3 * gamma_gg_nf3(n, cache, is_singlet) ) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py index 4752328dc..83e976234 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py @@ -3,11 +3,12 @@ import numba as nb import numpy as np +from .....harmonics import cache as c from .....harmonics.log_functions import lm13, lm13m1, lm14, lm15 @nb.njit(cache=True) -def gamma_gq_nf3(n, sx): +def gamma_gq_nf3(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^3` of :math:`\gamma_{gq}^{(3)}`, the expression is copied exact from Eq. 3.13 of :cite:`Davies:2016jie`. @@ -15,18 +16,19 @@ def gamma_gq_nf3(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- complex |N3LO| non-singlet anomalous dimension :math:`\gamma_{gq}^{(3)}|_{nf^3}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) return 1.3333333333333333 * ( -11.39728026699467 / (-1.0 + n) + 11.39728026699467 / n @@ -52,15 +54,17 @@ def gamma_gq_nf3(n, sx): @nb.njit(cache=True) -def gamma_gq_nf0(n, sx): +def gamma_gq_nf0(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^0` of :math:`\gamma_{gq}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -68,11 +72,11 @@ def gamma_gq_nf0(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{gq}^{(3)}|_{nf^0}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] - S5 = sx[4][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S5 = c.get(c.S5, cache, n, is_singlet) return ( -22156.31283903764 / np.power(-1.0 + n, 4) + 63019.91215580799 / np.power(-1.0 + n, 3) @@ -88,15 +92,17 @@ def gamma_gq_nf0(n, sx): @nb.njit(cache=True) -def gamma_gq_nf1(n, sx): +def gamma_gq_nf1(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{gq}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -104,11 +110,11 @@ def gamma_gq_nf1(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{gq}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] - S5 = sx[4][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S5 = c.get(c.S5, cache, n, is_singlet) return ( -4989.9438192798825 / np.power(-1.0 + n, 3) + 9496.873384515262 / np.power(-1.0 + n, 2) @@ -123,15 +129,17 @@ def gamma_gq_nf1(n, sx): @nb.njit(cache=True) -def gamma_gq_nf2(n, sx): +def gamma_gq_nf2(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{gq}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -139,10 +147,10 @@ def gamma_gq_nf2(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{gq}^{(3)}|_{nf^2}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) return ( -215.9801828033175 / np.power(-1.0 + n, 2) - 18.066114610010438 / (-1.0 + n) @@ -154,7 +162,7 @@ def gamma_gq_nf2(n, sx): @nb.njit(cache=True) -def gamma_gq(n, nf, sx): +def gamma_gq(n, nf, cache, is_singlet): r"""Computes the |N3LO| gluon-quark singlet anomalous dimension. Parameters @@ -163,8 +171,10 @@ def gamma_gq(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -181,8 +191,8 @@ def gamma_gq(n, nf, sx): """ return ( - gamma_gq_nf0(n, sx) - + nf * gamma_gq_nf1(n, sx) - + nf**2 * gamma_gq_nf2(n, sx) - + nf**3 * gamma_gq_nf3(n, sx) + gamma_gq_nf0(n, cache, is_singlet) + + nf * gamma_gq_nf1(n, cache, is_singlet) + + nf**2 * gamma_gq_nf2(n, cache, is_singlet) + + nf**3 * gamma_gq_nf3(n, cache, is_singlet) ) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py index 336d5b568..5371a3c11 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py @@ -3,13 +3,14 @@ """ import numba as nb -from eko.constants import CF -from .....harmonics.constants import zeta3 +from eko.constants import CF, zeta3 + +from .....harmonics import cache as c from .....harmonics.log_functions import lm11m1, lm12m1, lm13m1 @nb.njit(cache=True) -def gamma_ns_nf3(n, sx): +def gamma_ns_nf3(n, cache, is_singlet): """Implements the common part proportional to :math:`nf^3`, of :math:`\\gamma_{ns,+}^{(3)},\\gamma_{ns,-}^{(3)},\\gamma_{ns,v}^{(3)}`. The expression is copied exact from Eq. 3.6. of :cite:`Davies:2016jie`. @@ -18,8 +19,10 @@ def gamma_ns_nf3(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -27,10 +30,10 @@ def gamma_ns_nf3(n, sx): |N3LO| non-singlet anomalous dimension :math:`\\gamma_{ns}^{(3)}|_{nf^3}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) eta = 1 / n * 1 / (n + 1) g_ns_nf3 = CF * ( -32 / 27 * zeta3 * eta @@ -52,15 +55,17 @@ def gamma_ns_nf3(n, sx): @nb.njit(cache=True) -def gamma_nsm_nf2(n, sx): +def gamma_nsm_nf2(n, cache, is_singlet): """Implements the parametrized valence-like non-singlet part proportional to :math:`nf^2`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -69,10 +74,12 @@ def gamma_nsm_nf2(n, sx): :math:`\\gamma_{ns,-}^{(3)}|_{nf^2}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) Lm11m1 = lm11m1(n, S1) - Lm12m1 = lm12m1(n, S1, sx[1][0]) - Lm13m1 = lm13m1(n, S1, sx[1][0], sx[2][0]) + Lm12m1 = lm12m1(n, S1, S2) + Lm13m1 = lm13m1(n, S1, S2, S3) return ( -193.85692903712987 - 18.962962962962962 / n**5 @@ -93,15 +100,17 @@ def gamma_nsm_nf2(n, sx): @nb.njit(cache=True) -def gamma_nsm_nf1(n, sx): +def gamma_nsm_nf1(n, cache, is_singlet): """Implements the parametrized valence-like non-singlet part proportional to :math:`nf^1`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -110,10 +119,12 @@ def gamma_nsm_nf1(n, sx): :math:`\\gamma_{ns,-}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) Lm11m1 = lm11m1(n, S1) - Lm12m1 = lm12m1(n, S1, sx[1][0]) - Lm13m1 = lm13m1(n, S1, sx[1][0], sx[2][0]) + Lm12m1 = lm12m1(n, S1, S2) + Lm13m1 = lm13m1(n, S1, S2, S3) return ( 5549.7951398549685 - 126.41975308641975 / n**6 @@ -135,15 +146,17 @@ def gamma_nsm_nf1(n, sx): @nb.njit(cache=True) -def gamma_nsm_nf0(n, sx): +def gamma_nsm_nf0(n, cache, is_singlet): """Implements the parametrized valence-like non-singlet part proportional to :math:`nf^0`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -151,10 +164,12 @@ def gamma_nsm_nf0(n, sx): |N3LO| valence-like non-singlet anomalous dimension :math:`\\gamma_{ns,-}^{(3)}|_{nf^0}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) Lm11m1 = lm11m1(n, S1) - Lm12m1 = lm12m1(n, S1, sx[1][0]) - Lm13m1 = lm13m1(n, S1, sx[1][0], sx[2][0]) + Lm12m1 = lm12m1(n, S1, S2) + Lm13m1 = lm13m1(n, S1, S2, S3) return ( -23383.08164724965 - 252.8395061728395 / n**7 @@ -177,7 +192,7 @@ def gamma_nsm_nf0(n, sx): @nb.njit(cache=True) -def gamma_nsm(n, nf, sx): +def gamma_nsm(n, nf, cache, is_singlet): """Computes the |N3LO| valence-like non-singlet anomalous dimension. Parameters @@ -186,8 +201,10 @@ def gamma_nsm(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -204,8 +221,8 @@ def gamma_nsm(n, nf, sx): """ return ( - gamma_nsm_nf0(n, sx) - + nf * gamma_nsm_nf1(n, sx) - + nf**2 * gamma_nsm_nf2(n, sx) - + nf**3 * gamma_ns_nf3(n, sx) + gamma_nsm_nf0(n, cache, is_singlet) + + nf * gamma_nsm_nf1(n, cache, is_singlet) + + nf**2 * gamma_nsm_nf2(n, cache, is_singlet) + + nf**3 * gamma_ns_nf3(n, cache, is_singlet) ) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py index 9f5be5c4d..84b4e78a1 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py @@ -3,20 +3,23 @@ """ import numba as nb +from .....harmonics import cache as c from .....harmonics.log_functions import lm11m1, lm12m1, lm13m1 from .gnsm import gamma_ns_nf3 @nb.njit(cache=True) -def gamma_nsp_nf2(n, sx): +def gamma_nsp_nf2(n, cache, is_singlet): """Implements the parametrized singlet-like non-singlet part proportional to :math:`nf^2`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -24,10 +27,12 @@ def gamma_nsp_nf2(n, sx): |N3LO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(3)}|_{nf^2}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) Lm11m1 = lm11m1(n, S1) - Lm12m1 = lm12m1(n, S1, sx[1][0]) - Lm13m1 = lm13m1(n, S1, sx[1][0], sx[2][0]) + Lm12m1 = lm12m1(n, S1, S2) + Lm13m1 = lm13m1(n, S1, S2, S3) return ( -193.85479604848626 - 18.962962962962962 / n**5 @@ -48,15 +53,17 @@ def gamma_nsp_nf2(n, sx): @nb.njit(cache=True) -def gamma_nsp_nf1(n, sx): +def gamma_nsp_nf1(n, cache, is_singlet): """Implements the parametrized singlet-like non-singlet part proportional to :math:`nf^1`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -64,10 +71,12 @@ def gamma_nsp_nf1(n, sx): |N3LO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) Lm11m1 = lm11m1(n, S1) - Lm12m1 = lm12m1(n, S1, sx[1][0]) - Lm13m1 = lm13m1(n, S1, sx[1][0], sx[2][0]) + Lm12m1 = lm12m1(n, S1, S2) + Lm13m1 = lm13m1(n, S1, S2, S3) return ( 5550.063827367692 - 126.41975308641975 / n**6 @@ -89,15 +98,17 @@ def gamma_nsp_nf1(n, sx): @nb.njit(cache=True) -def gamma_nsp_nf0(n, sx): +def gamma_nsp_nf0(n, cache, is_singlet): """Implements the parametrized singlet-like non-singlet part proportional to :math:`nf^0`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -105,10 +116,12 @@ def gamma_nsp_nf0(n, sx): |N3LO| singlet-like non-singlet anomalous dimension :math:`\\gamma_{ns,+}^{(3)}|_{nf^0}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) Lm11m1 = lm11m1(n, S1) - Lm12m1 = lm12m1(n, S1, sx[1][0]) - Lm13m1 = lm13m1(n, S1, sx[1][0], sx[2][0]) + Lm12m1 = lm12m1(n, S1, S2) + Lm13m1 = lm13m1(n, S1, S2, S3) return ( -23391.854890259732 - 252.8395061728395 / n**7 @@ -131,7 +144,7 @@ def gamma_nsp_nf0(n, sx): @nb.njit(cache=True) -def gamma_nsp(n, nf, sx): +def gamma_nsp(n, nf, cache, is_singlet): """Computes the |N3LO| singlet-like non-singlet anomalous dimension. Parameters @@ -140,8 +153,10 @@ def gamma_nsp(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -158,8 +173,8 @@ def gamma_nsp(n, nf, sx): """ return ( - gamma_nsp_nf0(n, sx) - + nf * gamma_nsp_nf1(n, sx) - + nf**2 * gamma_nsp_nf2(n, sx) - + nf**3 * gamma_ns_nf3(n, sx) + gamma_nsp_nf0(n, cache, is_singlet) + + nf * gamma_nsp_nf1(n, cache, is_singlet) + + nf**2 * gamma_nsp_nf2(n, cache, is_singlet) + + nf**3 * gamma_ns_nf3(n, cache, is_singlet) ) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsv.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsv.py index 5d831590f..99a57c4a2 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsv.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsv.py @@ -3,11 +3,12 @@ """ import numba as nb +from .....harmonics import cache as c from .gnsm import gamma_nsm @nb.njit(cache=True) -def gamma_nss_nf2(n, sx): +def gamma_nss_nf2(n, cache, is_singlet): """Implements the sea non-singlet part proportional to :math:`nf^2` as in Eq. 3.5 of :cite:`Davies:2016jie`. @@ -15,8 +16,10 @@ def gamma_nss_nf2(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -24,10 +27,18 @@ def gamma_nss_nf2(n, sx): |N3LO| sea non-singlet anomalous dimension :math:`\\gamma_{ns,s}^{(3)}|_{nf^2}` """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, _, _, Sm21, _, Sm3 = sx[2] - S4, S31, _, Sm22, Sm211, Sm31, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm31 = c.get(c.Sm31, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) return ( 160 / 27 @@ -163,7 +174,7 @@ def gamma_nss_nf2(n, sx): @nb.njit(cache=True) -def gamma_nss_nf1(n, sx): +def gamma_nss_nf1(n, cache, is_singlet): """Implements the sea non-singlet part proportional to :math:`nf^1`. The expression is the average of the Mellin transform of Eq. 4.19, 4.20 of :cite:`Moch:2017uml` @@ -172,8 +183,10 @@ def gamma_nss_nf1(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -181,9 +194,9 @@ def gamma_nss_nf1(n, sx): |N3LO| sea non-singlet anomalous dimension :math:`\\gamma_{ns,s}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) return ( 3803.04 / n**6 - 3560.16 / n**5 @@ -218,7 +231,7 @@ def gamma_nss_nf1(n, sx): @nb.njit(cache=True) -def gamma_nsv(n, nf, sx): +def gamma_nsv(n, nf, cache, is_singlet): """Computes the |N3LO| valence non-singlet anomalous dimension. Parameters @@ -227,8 +240,10 @@ def gamma_nsv(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -244,7 +259,7 @@ def gamma_nsv(n, nf, sx): """ return ( - gamma_nsm(n, nf, sx) - + nf * gamma_nss_nf1(n, sx) - + nf**2 * gamma_nss_nf2(n, sx) + gamma_nsm(n, nf, cache, is_singlet) + + nf * gamma_nss_nf1(n, cache, is_singlet) + + nf**2 * gamma_nss_nf2(n, cache, is_singlet) ) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gps.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gps.py index a563a693c..bf777ab21 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gps.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gps.py @@ -3,9 +3,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def gamma_ps_nf3(n, sx): +def gamma_ps_nf3(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^3` of :math:`\gamma_{ps}^{(3)}`, the expression is copied exact from Eq. 3.10 of :cite:`Davies:2016jie`. @@ -13,8 +15,10 @@ def gamma_ps_nf3(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -22,9 +26,9 @@ def gamma_ps_nf3(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{ps}^{(3)}|_{nf^3}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) return 1.3333333333333333 * ( 16.305796943701882 / (-1.0 + n) + 3.5555555555555554 / np.power(n, 5) @@ -77,23 +81,24 @@ def gamma_ps_nf3(n, sx): @nb.njit(cache=True) -def gamma_ps_nf1(n, sx): +def gamma_ps_nf1(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{ps}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache - + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- complex |N3LO| non-singlet anomalous dimension :math:`\gamma_{ps}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) return ( -10185.956311680911 * (1 / (-1.0 + n) - 1.0 / n) - 3498.454512979491 / np.power(-1.0 + n, 3) @@ -107,15 +112,17 @@ def gamma_ps_nf1(n, sx): @nb.njit(cache=True) -def gamma_ps_nf2(n, sx): +def gamma_ps_nf2(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{ps}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -123,7 +130,7 @@ def gamma_ps_nf2(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{ps}^{(3)}|_{nf^2}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) return ( -547.9815124741111 * (1 / (-1.0 + n) - 1.0 / n) + 383.72024118761027 / np.power(-1.0 + n, 2) @@ -136,7 +143,7 @@ def gamma_ps_nf2(n, sx): @nb.njit(cache=True) -def gamma_ps(n, nf, sx): +def gamma_ps(n, nf, cache, is_singlet): r"""Computes the |N3LO| pure singlet quark-quark anomalous dimension. Parameters @@ -145,8 +152,10 @@ def gamma_ps(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -162,7 +171,7 @@ def gamma_ps(n, nf, sx): """ return ( - +nf * gamma_ps_nf1(n, sx) - + nf**2 * gamma_ps_nf2(n, sx) - + nf**3 * gamma_ps_nf3(n, sx) + +nf * gamma_ps_nf1(n, cache, is_singlet) + + nf**2 * gamma_ps_nf2(n, cache, is_singlet) + + nf**3 * gamma_ps_nf3(n, cache, is_singlet) ) diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py index 7f146eebc..8add77e47 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py +++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py @@ -3,11 +3,12 @@ import numba as nb import numpy as np +from .....harmonics import cache as c from .....harmonics.log_functions import lm13, lm13m1, lm14, lm15 @nb.njit(cache=True) -def gamma_qg_nf3(n, sx): +def gamma_qg_nf3(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^3` of :math:`\gamma_{qg}^{(3)}`, the expression is copied exact from Eq. 3.12 of :cite:`Davies:2016jie`. @@ -15,8 +16,10 @@ def gamma_qg_nf3(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -24,10 +27,16 @@ def gamma_qg_nf3(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{qg}^{(3)}|_{nf^3}` """ - S1 = sx[0][0] - S2, Sm2 = sx[1] - S3, S21, _, _, _, Sm3 = sx[2] - S4, S31, S211, _, _, _, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) return 1.3333333333333333 * ( 44.56685134331718 / (-1.0 + n) - 82.37037037037037 / np.power(n, 5) @@ -329,15 +338,17 @@ def gamma_qg_nf3(n, sx): @nb.njit(cache=True) -def gamma_qg_nf1(n, sx): +def gamma_qg_nf1(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{qg}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -345,11 +356,11 @@ def gamma_qg_nf1(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{qg}^{(3)}|_{nf^1}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] - S5 = sx[4][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S5 = c.get(c.S5, cache, n, is_singlet) return ( -7871.5226542038545 / np.power(-1.0 + n, 3) + 13490.08729769531 / np.power(-1.0 + n, 2) @@ -365,15 +376,17 @@ def gamma_qg_nf1(n, sx): @nb.njit(cache=True) -def gamma_qg_nf2(n, sx): +def gamma_qg_nf2(n, cache, is_singlet): r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{qg}^{(3)}`. Parameters ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -381,11 +394,11 @@ def gamma_qg_nf2(n, sx): |N3LO| non-singlet anomalous dimension :math:`\gamma_{qg}^{(3)}|_{nf^2}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] - S5 = sx[4][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S5 = c.get(c.S5, cache, n, is_singlet) return ( 237.02932710506118 / np.power(-1.0 + n, 2) + 447.1862550338544 / (-1.0 + n) @@ -400,7 +413,7 @@ def gamma_qg_nf2(n, sx): @nb.njit(cache=True) -def gamma_qg(n, nf, sx): +def gamma_qg(n, nf, cache, is_singlet): r"""Computes the |N3LO| quark-gluon singlet anomalous dimension. Parameters @@ -409,8 +422,10 @@ def gamma_qg(n, nf, sx): Mellin moment nf : int Number of active flavors - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -426,7 +441,7 @@ def gamma_qg(n, nf, sx): """ return ( - +nf * gamma_qg_nf1(n, sx) - + nf**2 * gamma_qg_nf2(n, sx) - + nf**3 * gamma_qg_nf3(n, sx) + +nf * gamma_qg_nf1(n, cache, is_singlet) + + nf**2 * gamma_qg_nf2(n, cache, is_singlet) + + nf**3 * gamma_qg_nf3(n, cache, is_singlet) ) diff --git a/src/ekore/harmonics/__init__.py b/src/ekore/harmonics/__init__.py index 919ca559b..b99ae4dea 100644 --- a/src/ekore/harmonics/__init__.py +++ b/src/ekore/harmonics/__init__.py @@ -11,16 +11,12 @@ from .w4 import S4, S31, S211, Sm4, Sm22, Sm31, Sm211 from .w5 import S5, Sm5 - @nb.njit(cache=True) def base_harmonics_cache(n, is_singlet, max_weight=5, n_max_sums_weight=7): r"""Get the harmonics sums S basic cache. - Only single index harmonics are computed and stored in the first (:math:`S_{n}`) or in the last column (:math:`S_{-n}`). - Multi indices harmonics sums can be stored in the middle columns. - Parameters ---------- n : complex @@ -32,36 +28,31 @@ def base_harmonics_cache(n, is_singlet, max_weight=5, n_max_sums_weight=7): max harmonics weight, max value 5 (default) n_max_sums_weight : int max number of harmonics sums for a given weight - Returns ------- np.ndarray list of harmonics sums: (weights, n_max_sums_weight) - """ h_cache = np.zeros((max_weight, n_max_sums_weight), dtype=np.complex_) h_cache[:, 0] = sx(n, max_weight) - if n_max_sums_weight > 1: - h_cache[:, -1] = smx(n, h_cache[:, 0], is_singlet) +# if n_max_sums_weight > 1: +# h_cache[:, -1] = smx(n, h_cache[:, 0], is_singlet) return h_cache @nb.njit(cache=True) def sx(n, max_weight=5): """Get the harmonics sums S cache. - Parameters ---------- n : complex Mellin moment max_weight : int max harmonics weight, max value 5 (default) - Returns ------- np.ndarray list of harmonics sums (:math:`S_{1,..,w}`) - """ sx = np.zeros(max_weight, dtype=np.complex_) if max_weight >= 1: @@ -77,45 +68,10 @@ def sx(n, max_weight=5): return sx -@nb.njit(cache=True) -def smx(n, sx, is_singlet): - r"""Get the harmonics S-minus cache. - - Parameters - ---------- - n : complex - Mellin moment - sx : numpy.ndarray - List of harmonics sums: :math:`S_{1},\dots,S_{w}` - is_singlet : bool - symmetry factor: True for singlet like quantities (:math:`\eta=(-1)^N = 1`), - False for non-singlet like quantities (:math:`\eta=(-1)^N=-1`) - - Returns - ------- - np.ndarray - list of harmonics sums (:math:`S_{-1,..,-w}`) - - """ - max_weight = sx.size - smx = np.zeros(max_weight, dtype=np.complex_) - if max_weight >= 1: - smx[0] = Sm1(n, sx[0], is_singlet) - if max_weight >= 2: - smx[1] = Sm2(n, sx[1], is_singlet) - if max_weight >= 3: - smx[2] = Sm3(n, sx[2], is_singlet) - if max_weight >= 4: - smx[3] = Sm4(n, sx[3], is_singlet) - if max_weight >= 5: - smx[4] = Sm5(n, sx[4], is_singlet) - return smx - @nb.njit(cache=True) def s3x(n, sx, smx, is_singlet): r"""Compute the weight 3 multi indices harmonics sums cache. - Parameters ---------- n : complex @@ -127,12 +83,10 @@ def s3x(n, sx, smx, is_singlet): is_singlet : bool symmetry factor: True for singlet like quantities (:math:`\eta=(-1)^N = 1`), False for non-singlet like quantities (:math:`\eta=(-1)^N=-1`) - Returns ------- np.ndarray list containing: :math:`S_{2,1},S_{2,-1},S_{-2,1},S_{-2,-1}` - """ return np.array( [ @@ -147,7 +101,6 @@ def s3x(n, sx, smx, is_singlet): @nb.njit(cache=True) def s4x(n, sx, smx, is_singlet): r"""Compute the weight 4 multi indices harmonics sums cache. - Parameters ---------- n : complex @@ -159,12 +112,10 @@ def s4x(n, sx, smx, is_singlet): is_singlet : bool symmetry factor: True for singlet like quantities (:math:`\eta=(-1)^N = 1`), False for non-singlet like quantities (:math:`\eta=(-1)^N=-1`) - Returns ------- np.ndarray list containing: :math:`S_{3,1},S_{2,1,1},S_{-2,2},S_{-2,1,1},S_{-3,1}` - """ sm31 = Sm31(n, sx[0], smx[0], smx[1], is_singlet) return np.array( @@ -181,7 +132,6 @@ def s4x(n, sx, smx, is_singlet): @nb.njit(cache=True) def compute_cache(n, max_weight, is_singlet): r"""Get the harmonics sums cache. - Parameters ---------- n : complex @@ -191,19 +141,16 @@ def compute_cache(n, max_weight, is_singlet): is_singlet : bool symmetry factor: True for singlet like quantities (:math:`\eta=(-1)^N = 1`), False for non-singlet like quantities (:math:`\eta=(-1)^N=-1`) - Returns ------- list harmonic sums cache. At |N3LO| it contains: - .. math :: [[S_1,S_{-1}], [S_2,S_{-2}], [S_{3}, S_{2,1}, S_{2,-1}, S_{-2,1}, S_{-2,-1}, S_{-3}], [S_{4}, S_{3,1}, S_{2,1,1}, S_{-2,-2}, S_{-3, 1}, S_{-4}],] [S_{5}, S_{-5}] - """ # max number of harmonics sum of a given weight for a given max weight. n_max_sums_weight = {2: 1, 3: 3, 4: 7, 5: 7} @@ -216,4 +163,4 @@ def compute_cache(n, max_weight, is_singlet): sx[2, 1:-2] = s3x(n, sx[:, 0], sx[:, -1], is_singlet) sx[3, 1:-1] = s4x(n, sx[:, 0], sx[:, -1], is_singlet) # return list of list keeping the non zero values - return [[el for el in sx_list if el != 0] for sx_list in sx] + return [[el for el in sx_list if el != 0] for sx_list in sx] \ No newline at end of file diff --git a/src/ekore/harmonics/cache.py b/src/ekore/harmonics/cache.py new file mode 100644 index 000000000..845e14ebc --- /dev/null +++ b/src/ekore/harmonics/cache.py @@ -0,0 +1,173 @@ +"""Caching (complicated) harmonic sums across :mod:`ekore`.""" +import numba as nb +import numpy as np +import numpy.typing as npt + +from . import w1, w2, w3, w4, w5 +from .g_functions import mellin_g3, mellin_g3p1, mellin_g3p2 +from .polygamma import recursive_harmonic_sum + +# here a register of all possible functions +S1 = 0 # = S_1(N) +S2 = 1 +S3 = 2 +S4 = 3 +S5 = 4 +Sm1 = 5 +Sm2 = 6 +Sm3 = 7 +Sm4 = 8 +Sm5 = 9 +S21 = 10 +S2m1 = 11 +Sm21 = 12 +Sm2m1 = 13 +S31 = 14 +Sm31 = 15 +Sm22 = 16 +S211 = 17 +Sm211 = 18 +S1h = 19 +S2h = 20 +S3h = 21 +S4h = 22 +S5h = 23 +S1mh = 24 +S2mh = 25 +S3mh = 26 +S4mh = 27 +S5mh = 28 +S1ph = 29 +S2ph = 30 +S3ph = 31 +S4ph = 32 +S5ph = 33 +g3 = 34 +S1p2 = 35 +g3p1 = 36 +g3p2 = 37 + +# this could also be S1h = S1(N/2) +# the only requirement is that they are insubsequent order +# and reset knows the maximum size (to fit them all) + +# this is a plain list holding the values +@nb.njit(cache=True) +def reset(): + """Return the cache placeholder array.""" + return np.full(38, np.nan, np.complex_) + + +@nb.njit(cache=True) +def get(key: int, cache: npt.ArrayLike, + n: complex, is_singlet: bool) -> complex: + """Retrieve an element of the cache. + + Parameters + ---------- + key : + harmonic sum key + cache : + cache list holding all elements + n : + complex evaluation point + is_singlet : + symmetry factor: True for singlet like quantities + (:math:`\eta=(-1)^N = 1`), + False for non-singlet like quantities + (:math:`\eta=(-1)^N=-1`) + + Returns + ------- + complex : + requested harmonic sum evaluated at n + + """ + # Maybe improve error + if key < 0 or key > len(cache): + raise RuntimeError + # load the thing + s = cache[key] + # found? i.e. not NaN? + if not np.isnan(s): + return s + # compute it now ... + if key == S1: + s = w1.S1(n) + elif key == S2: + s = w2.S2(n) + elif key == S3: + s = w3.S3(n) + elif key == S4: + s = w4.S4(n) + elif key == S5: + s = w5.S5(n) + elif key == Sm1: + s = w1.Sm1(n, get(S1, cache, n, is_singlet), get(S1mh, cache, n, is_singlet), get(S1h, cache, n, is_singlet), is_singlet) + elif key == Sm2: + s = w2.Sm2(n, get(S2, cache, n, is_singlet), get(S2mh, cache, n, is_singlet), get(S2h, cache, n, is_singlet), is_singlet) + elif key == Sm3: + s = w3.Sm3(n, get(S3, cache, n, is_singlet), get(S3mh, cache, n, is_singlet), get(S3h, cache, n, is_singlet), is_singlet) + elif key == Sm4: + s = w4.Sm4(n, get(S4, cache, n, is_singlet), get(S4mh, cache, n, is_singlet), get(S4h, cache, n, is_singlet), is_singlet) + elif key == Sm5: + s = w5.Sm5(n, get(S5, cache, n, is_singlet), get(S5mh, cache, n, is_singlet), get(S5h, cache, n, is_singlet), is_singlet) + elif key == S21: + s = w3.S21(n, get(S1, cache, n, is_singlet), get(S2, cache, n, is_singlet)) + elif key == S2m1: + s = w3.S2m1(n, get(S2, cache, n, is_singlet), get(Sm1, cache, n, is_singlet), get(Sm2, cache, n, is_singlet), is_singlet) + elif key == Sm21: + s = w3.Sm21(n, get(S1, cache, n, is_singlet), get(Sm1, cache, n, is_singlet), get(g3p1, cache, n, is_singlet), is_singlet) + elif key == Sm2m1: + s = w3.Sm2m1(n, get(S1, cache, n, is_singlet), get(S2, cache, n, is_singlet), get(Sm2, cache, n, is_singlet)) + elif key == S31: + s = w4.S31(n, get(S1, cache, n, is_singlet), get(S2, cache, n, is_singlet), get(S3, cache, n, is_singlet), get(S4, cache, n, is_singlet)) + elif key == Sm31: + s = w4.Sm31(n, get(S1, cache, n, is_singlet), get(Sm1, cache, n, is_singlet), get(Sm2, cache, n, is_singlet), is_singlet) + elif key == Sm22: + s = w4.Sm22(n, get(S1, cache, n, is_singlet), get(S2, cache, n, is_singlet), get(Sm2, cache, n, is_singlet), get(Sm31, cache, n, is_singlet), is_singlet) + elif key == S211: + s = w4.S211(n, get(S1, cache, n, is_singlet), get(S2, cache, n, is_singlet), get(S3, cache, n, is_singlet)) + elif key == Sm211: + s = w4.Sm211(n, get(S1, cache, n, is_singlet), get(S2, cache, n, is_singlet), get(Sm1, cache, n, is_singlet), is_singlet) + elif key == S1h: + s = w1.S1(n/2) + elif key == S2h: + s = w2.S2(n/2) + elif key == S3h: + s = w3.S3(n/2) + elif key == S4h: + s = w4.S4(n/2) + elif key == S5h: + s = w5.S5(n/2) + elif key == S1mh: + s = w1.S1((n-1)/2) + elif key == S2mh: + s = w2.S2((n-1)/2) + elif key == S3mh: + s = w3.S3((n-1)/2) + elif key == S4mh: + s = w4.S4((n-1)/2) + elif key == S5mh: + s = w5.S5((n-1)/2) + elif key == S1ph: + s = recursive_harmonic_sum(get(S1mh, cache, n, is_singlet), n, 1, 1) + elif key == S2ph: + s = recursive_harmonic_sum(get(S2mh, cache, n, is_singlet), n, 1, 2) + elif key == S3ph: + s = recursive_harmonic_sum(get(S3mh, cache, n, is_singlet), n, 1, 3) + elif key == S4ph: + s = recursive_harmonic_sum(get(S4mh, cache, n, is_singlet), n, 1, 4) + elif key == S5ph: + s = recursive_harmonic_sum(get(S5mh, cache, n, is_singlet), n, 1, 5) + elif key == g3: + s = mellin_g3(n, get(S1, cache, n, is_singlet)) + elif key == S1p2: + s = recursive_harmonic_sum(get(S1, cache, n, is_singlet), n, 2, 1) + elif key == g3p1: + s = mellin_g3p1(n, get(S1, cache, n, is_singlet), get(g3, cache, n, is_singlet)) + elif key == g3p2: + s = mellin_g3p2(n, get(S1, cache, n, is_singlet), get(g3, cache, n, is_singlet)) + # store and return + cache[key] = s + return s diff --git a/src/ekore/harmonics/constants.py b/src/ekore/harmonics/constants.py deleted file mode 100644 index df64544dd..000000000 --- a/src/ekore/harmonics/constants.py +++ /dev/null @@ -1,22 +0,0 @@ -"""Zeta function values and other constants.""" - -import numpy as np -from scipy.special import zeta - -zeta2 = zeta(2) -r""":math:`\zeta(2)`""" - -zeta3 = zeta(3) -r""":math:`\zeta(3)`""" - -zeta4 = zeta(4) -r""":math:`\zeta(4)`""" - -zeta5 = zeta(5) -r""":math:`\zeta(5)`""" - -log2 = np.log(2) -r""":math:`\ln(2)`""" - -li4half = 0.517479 -""":math:`Li_{4}(1/2)`""" diff --git a/src/ekore/harmonics/g_functions.py b/src/ekore/harmonics/g_functions.py index dc7073175..0561d299b 100644 --- a/src/ekore/harmonics/g_functions.py +++ b/src/ekore/harmonics/g_functions.py @@ -5,8 +5,8 @@ import numba as nb import numpy as np +from eko.constants import log2, zeta2, zeta3 from . import w1 -from .constants import log2, zeta2, zeta3 from .polygamma import recursive_harmonic_sum as s a1 = np.array( @@ -95,6 +95,54 @@ def mellin_g3(N, S1): return g3 +@nb.njit(cache=True) +def mellin_g3p1(N, s1, g3): + """ + Computes :math:`g_{3}(N+1)` + + Parameters + ---------- + N : complex + Mellin moment + s1 : complex + Harmonic sum :math:`S_{1}(N)` + mellin_g3 : complex + :math:`g_{3}(N)` + + Returns + ------- + mellin_g3p1 : complex + :math:`g_{3}(N+1)` + + """ + g3p1 = - s1/np.power(N, 2) - g3 + zeta2/N + return g3p1 + +@nb.njit(cache=True) +def mellin_g3p2(N, s1, g3): + """ + Computes :math:`g_{3}(N+2)` + + Parameters + ---------- + N : complex + Mellin moment + s1 : complex + Harmonic sum :math:`S_{1}(N)` + mellin_g3 : complex + :math:`g_{3}(N)` + + Returns + ------- + mellin_g3p2 : complex + :math:`g_{3}(N+2)` + + """ + g3p2 = (s1/np.power(N, 2) + g3 - zeta2/N + zeta2/(N + 1) + - s1/np.power(N + 1, 2) - 1/np.power(N + 1, 3)) + return g3p2 + + @nb.njit(cache=True) def mellin_g4(N): r""" diff --git a/src/ekore/harmonics/log_functions.py b/src/ekore/harmonics/log_functions.py index 8cf1bf857..62a563ccf 100644 --- a/src/ekore/harmonics/log_functions.py +++ b/src/ekore/harmonics/log_functions.py @@ -8,8 +8,7 @@ """ import numba as nb -from .constants import zeta3 - +from eko.constants import zeta3 @nb.njit(cache=True) def lm11m1(n, S1): diff --git a/src/ekore/harmonics/w1.py b/src/ekore/harmonics/w1.py index 10f7cc3d4..5d5914e51 100644 --- a/src/ekore/harmonics/w1.py +++ b/src/ekore/harmonics/w1.py @@ -34,7 +34,7 @@ def S1(N): @nb.njit(cache=True) -def Sm1(N, hS1, is_singlet=None): +def Sm1(N, hS1, hS1mh, hS1h, is_singlet=None): r"""Analytic continuation of harmonic sum :math:`S_{-1}(N)`. .. math:: @@ -60,10 +60,10 @@ def Sm1(N, hS1, is_singlet=None): """ if is_singlet is None: return ( - (1 - (-1) ** N) / 2 * S1((N - 1) / 2) - + ((-1) ** N + 1) / 2 * S1(N / 2) + (1 - (-1) ** N) / 2 * hS1mh + + ((-1) ** N + 1) / 2 * hS1h - hS1 ) if is_singlet: - return S1(N / 2) - hS1 - return S1((N - 1) / 2) - hS1 + return hS1h - hS1 + return hS1mh - hS1 diff --git a/src/ekore/harmonics/w2.py b/src/ekore/harmonics/w2.py index 69c0029eb..d0561a3c0 100644 --- a/src/ekore/harmonics/w2.py +++ b/src/ekore/harmonics/w2.py @@ -2,7 +2,7 @@ import numba as nb -from .constants import zeta2 +from eko.constants import zeta2 from .polygamma import cern_polygamma @@ -34,7 +34,7 @@ def S2(N): @nb.njit(cache=True) -def Sm2(N, hS2, is_singlet=None): +def Sm2(N, hS2, hS2mh, hS2h, is_singlet=None): r"""Analytic continuation of harmonic sum :math:`S_{-2}(N)`. .. math:: @@ -64,9 +64,10 @@ def Sm2(N, hS2, is_singlet=None): return ( 1 / 2 - * ((1 - (-1) ** N) / 2 * S2((N - 1) / 2) + ((-1) ** N + 1) / 2 * S2(N / 2)) + * ((1 - (-1) ** N) / 2 * hS2mh + + ((-1) ** N + 1) / 2 * hS2h) - hS2 ) if is_singlet: - return 1 / 2 * S2(N / 2) - hS2 - return 1 / 2 * S2((N - 1) / 2) - hS2 + return 1 / 2 * hS2h - hS2 + return 1 / 2 * hS2mh - hS2 diff --git a/src/ekore/harmonics/w3.py b/src/ekore/harmonics/w3.py index dcd52428c..789ff7a9d 100644 --- a/src/ekore/harmonics/w3.py +++ b/src/ekore/harmonics/w3.py @@ -2,8 +2,8 @@ import numba as nb +from eko.constants import log2, zeta2, zeta3 from . import g_functions as gf -from .constants import log2, zeta2, zeta3 from .polygamma import cern_polygamma, symmetry_factor @@ -36,7 +36,7 @@ def S3(N): @nb.njit(cache=True) -def Sm3(N, hS3, is_singlet=None): +def Sm3(N, hS3, hS3mh, hS3h, is_singlet=None): r"""Analytic continuation of harmonic sum :math:`S_{-3}(N)`. .. math:: @@ -66,12 +66,13 @@ def Sm3(N, hS3, is_singlet=None): return ( 1 / 2**2 - * ((1 - (-1) ** N) / 2 * S3((N - 1) / 2) + ((-1) ** N + 1) / 2 * S3(N / 2)) + * ((1 - (-1) ** N) / 2 * hS3mh + + ((-1) ** N + 1) / 2 * hS3h) - hS3 ) if is_singlet: - return 1 / 2**2 * S3(N / 2) - hS3 - return 1 / 2**2 * S3((N - 1) / 2) - hS3 + return 1 / 2**2 * hS3h - hS3 + return 1 / 2**2 * hS3mh - hS3 @nb.njit(cache=True) @@ -104,7 +105,7 @@ def S21(N, S1, S2): @nb.njit(cache=True) -def Sm21(N, S1, Sm1, is_singlet=None): +def Sm21(N, S1, Sm1, g3p1, is_singlet=None): r"""Analytic continuation of harmonic sum :math:`S_{-2,1}(N)`. As implemented in eq B.5.75 of :cite:`MuselliPhD` and eq 22 of @@ -135,7 +136,7 @@ def Sm21(N, S1, Sm1, is_singlet=None): # Note mellin g3 was integrated following x^(N-1) convention. eta = symmetry_factor(N, is_singlet) return ( - -eta * gf.mellin_g3(N + 1, S1 + 1 / (N + 1)) + -eta * g3p1 + zeta2 * Sm1 - 5 / 8 * zeta3 + zeta2 * log2 diff --git a/src/ekore/harmonics/w4.py b/src/ekore/harmonics/w4.py index 855602ec7..1eca1708f 100644 --- a/src/ekore/harmonics/w4.py +++ b/src/ekore/harmonics/w4.py @@ -2,8 +2,9 @@ import numba as nb +from eko.constants import li4half, log2, zeta2, zeta3, zeta4 +from . import cache as c from . import g_functions as gf -from .constants import li4half, log2, zeta2, zeta3, zeta4 from .polygamma import cern_polygamma, symmetry_factor @@ -36,7 +37,7 @@ def S4(N): @nb.njit(cache=True) -def Sm4(N, hS4, is_singlet=None): +def Sm4(N, hS4, hS4mh, hS4h, is_singlet=None): r"""Analytic continuation of harmonic sum :math:`S_{-4}(N)`. .. math:: @@ -66,12 +67,13 @@ def Sm4(N, hS4, is_singlet=None): return ( 1 / 2**3 - * ((1 - (-1) ** N) / 2 * S4((N - 1) / 2) + ((-1) ** N + 1) / 2 * S4(N / 2)) + * ((1 - (-1) ** N) / 2 * hS4mh + + ((-1) ** N + 1) / 2 * hS4h) - hS4 ) if is_singlet: - return 1 / 2**3 * S4(N / 2) - hS4 - return 1 / 2**3 * S4((N - 1) / 2) - hS4 + return 1 / 2**3 * hS4h - hS4 + return 1 / 2**3 * hS4mh - hS4 @nb.njit(cache=True) diff --git a/src/ekore/harmonics/w5.py b/src/ekore/harmonics/w5.py index 0ca6a37b8..7ad815f0f 100644 --- a/src/ekore/harmonics/w5.py +++ b/src/ekore/harmonics/w5.py @@ -2,7 +2,8 @@ import numba as nb -from .constants import zeta5 +from eko.constants import zeta5 + from .polygamma import cern_polygamma @@ -35,7 +36,7 @@ def S5(N): @nb.njit(cache=True) -def Sm5(N, hS5, is_singlet=None): +def Sm5(N, hS5, hS5mh, hS5h, is_singlet=None): r"""Analytic continuation of harmonic sum :math:`S_{-5}(N)`. .. math:: @@ -63,11 +64,9 @@ def Sm5(N, hS5, is_singlet=None): """ if is_singlet is None: return ( - 1 - / 2**4 - * ((1 - (-1) ** N) / 2 * S5((N - 1) / 2) + ((-1) ** N + 1) / 2 * S5(N / 2)) + 1 / 2**4 * ((1 - (-1) ** N) / 2 * hS5mh + ((-1) ** N + 1) / 2 * hS5h) - hS5 ) if is_singlet: - return 1 / 2**4 * S5(N / 2) - hS5 - return 1 / 2**4 * S5((N - 1) / 2) - hS5 + return 1 / 2**4 * hS5h - hS5 + return 1 / 2**4 * hS5mh - hS5 diff --git a/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py b/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py index 5d5f771f7..0d961d2b9 100644 --- a/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py +++ b/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py @@ -4,10 +4,10 @@ @nb.njit(cache=True) -def A_non_singlet(_matching_order, _n, _sx, _nf, _L): +def A_non_singlet(_matching_order, _n, _nf, _L): raise NotImplementedError("Polarised, space-like is not yet implemented") @nb.njit(cache=True) -def A_singlet(_matching_order, _n, _sx, _nf, _L, _is_msbar, _sx_ns=None): +def A_singlet(_matching_order, _n, _nf, _L, _is_msbar): raise NotImplementedError("Polarised, space-like is not yet implemented") diff --git a/src/ekore/operator_matrix_elements/polarized/time_like/__init__.py b/src/ekore/operator_matrix_elements/polarized/time_like/__init__.py new file mode 100644 index 000000000..5573f69b0 --- /dev/null +++ b/src/ekore/operator_matrix_elements/polarized/time_like/__init__.py @@ -0,0 +1,13 @@ +r"""The polarized, time-like |OME|.""" + +import numba as nb + + +@nb.njit(cache=True) +def A_non_singlet(_matching_order, _n, _nf, _L): + raise NotImplementedError("Polarised, time-like is not yet implemented") + + +@nb.njit(cache=True) +def A_singlet(_matching_order, _n, _nf, _L, _is_msbar): + raise NotImplementedError("Polarised, time-like is not yet implemented") diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py index a1d2f8ec6..8e1f519fc 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py @@ -6,11 +6,12 @@ import numba as nb import numpy as np +from ....harmonics import cache as c from . import as1, as2, as3 @nb.njit(cache=True) -def A_singlet(matching_order, n, sx, nf, L, is_msbar, sx_ns=None): +def A_singlet(matching_order, n, nf, L, is_msbar): r""" Computes the tower of the singlet |OME|. @@ -20,41 +21,33 @@ def A_singlet(matching_order, n, sx, nf, L, is_msbar, sx_ns=None): perturbative matching_order n : complex Mellin variable - sx : list - singlet like harmonic sums cache nf: int number of active flavor below threshold L : float :math:``\ln(\mu_F^2 / m_h^2)`` is_msbar: bool add the |MSbar| contribution - sx_ns : list - non-singlet like harmonic sums cache Returns ------- A_singlet : numpy.ndarray singlet |OME| - See Also - -------- - ekore.matching_conditions.nlo.A_singlet_1 : :math:`A^{S,(1)}(N)` - ekore.matching_conditions.nlo.A_hh_1 : :math:`A_{HH}^{(1)}(N)` - ekore.matching_conditions.nlo.A_gh_1 : :math:`A_{gH}^{(1)}(N)` - ekore.matching_conditions.nnlo.A_singlet_2 : :math:`A_{S,(2)}(N)` """ A_s = np.zeros((matching_order[0], 3, 3), np.complex_) + + cache = c.reset() if matching_order[0] >= 1: - A_s[0] = as1.A_singlet(n, sx, L) + A_s[0] = as1.A_singlet(n, L, cache, True) if matching_order[0] >= 2: - A_s[1] = as2.A_singlet(n, sx, L, is_msbar) + A_s[1] = as2.A_singlet(n, L, cache, True, is_msbar) if matching_order[0] >= 3: - A_s[2] = as3.A_singlet(n, sx, sx_ns, nf, L) + A_s[2] = as3.A_singlet(n, nf, L, cache) return A_s @nb.njit(cache=True) -def A_non_singlet(matching_order, n, sx, nf, L): +def A_non_singlet(matching_order, n, nf, L): r""" Computes the tower of the non-singlet |OME| @@ -64,8 +57,6 @@ def A_non_singlet(matching_order, n, sx, nf, L): perturbative matching_order n : complex Mellin variable - sx : list - harmonic sums cache nf: int number of active flavor below threshold L : float @@ -76,16 +67,14 @@ def A_non_singlet(matching_order, n, sx, nf, L): A_non_singlet : numpy.ndarray non-singlet |OME| - See Also - -------- - ekore.matching_conditions.nlo.A_hh_1 : :math:`A_{HH}^{(1)}(N)` - ekore.matching_conditions.nnlo.A_ns_2 : :math:`A_{qq,H}^{NS,(2)}` """ A_ns = np.zeros((matching_order[0], 2, 2), np.complex_) + + cache = c.reset() if matching_order[0] >= 1: - A_ns[0] = as1.A_ns(n, sx, L) + A_ns[0] = as1.A_ns(n, L, cache, False) if matching_order[0] >= 2: - A_ns[1] = as2.A_ns(n, sx, L) + A_ns[1] = as2.A_ns(n, L, cache, False) if matching_order[0] >= 3: - A_ns[2] = as3.A_ns(n, sx, nf, L) + A_ns[2] = as3.A_ns(n, nf, L, cache, False) return A_ns diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py index c373f8091..41fea4c8b 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py @@ -11,9 +11,11 @@ from eko.constants import CF +from ....harmonics import cache as c + @nb.njit(cache=True) -def A_hh(n, sx, L): +def A_hh(n, L, cache, is_singlet): r""" |NLO| heavy-heavy |OME| :math:`A_{HH}^{(1)}` defined as the mellin transform of :math:`K_{hh}` given in Eq. (20a) of :cite:`Ball_2016`. @@ -22,18 +24,20 @@ def A_hh(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_hh : complex |NLO| heavy-heavy |OME| :math:`A_{HH}^{(1)}` """ - S1m = sx[0][0] - 1 / n # harmonics.S1(n - 1) - S2m = sx[1][0] - 1 / n**2 # harmonics.S2(n - 1) + S1m = c.get(c.S1, cache, n, is_singlet) - 1 / n # harmonics.S1(n - 1) + S2m = c.get(c.S2, cache, n, is_singlet) - 1 / n**2 # harmonics.S2(n - 1) ahh_l = (2 + n - 3 * n**2) / (n * (1 + n)) + 4 * S1m ahh = 2 * ( 2 @@ -115,7 +119,7 @@ def A_gg(L): @nb.njit(cache=True) -def A_singlet(n, sx, L): +def A_singlet(n, L, cache, is_singlet): r""" Computes the |NLO| singlet |OME|. @@ -130,10 +134,12 @@ def A_singlet(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache containing: [[:math:`S_1`][:math:`S_2`]] L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -151,7 +157,7 @@ def A_singlet(n, sx, L): [ [A_gg(L), 0.0, A_gh(n, L)], [0 + 0j, 0 + 0j, 0 + 0j], - [A_hg(n, L), 0.0, A_hh(n, sx, L)], + [A_hg(n, L), 0.0, A_hh(n, L, cache, is_singlet)], ], np.complex_, ) @@ -159,7 +165,7 @@ def A_singlet(n, sx, L): @nb.njit(cache=True) -def A_ns(n, sx, L): +def A_ns(n, L, cache, is_singlet): r""" Computes the |NLO| non-singlet |OME| with intrinsic contributions. @@ -173,10 +179,12 @@ def A_ns(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache containing: [[:math:`S_1`][:math:`S_2`]] L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_NS : numpy.ndarray @@ -186,4 +194,6 @@ def A_ns(n, sx, L): -------- A_hh : :math:`A_{HH}^{(1)}` """ - return np.array([[0 + 0j, 0 + 0j], [0 + 0j, A_hh(n, sx, L)]], np.complex_) + return np.array( + [[0 + 0j, 0 + 0j], [0 + 0j, A_hh(n, L, cache, is_singlet)]], np.complex_ + ) diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py index 2197fbfd9..bb3dcdb75 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py @@ -12,14 +12,15 @@ import numpy as np from eko import constants +from eko.constants import zeta2, zeta3 -from ....harmonics.constants import zeta2, zeta3 +from ....harmonics import cache as c from .as1 import A_gg as A_gg_1 from .as1 import A_hg as A_hg_1 @nb.njit(cache=True) -def A_qq_ns(n, sx, L): +def A_qq_ns(n, L, cache, is_singlet): r""" |NNLO| light-light non-singlet |OME| :math:`A_{qq,H}^{NS,(2)}` given in Eq. (B.4) of :cite:`Buza_1998`. @@ -28,19 +29,21 @@ def A_qq_ns(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_qq_ns : complex |NNLO| light-light non-singlet |OME| :math:`A_{qq,H}^{NS,(2)}` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) S1m = S1 - 1 / n # harmonic_S1(n - 1) S2m = S2 - 1 / n**2 # harmonic_S2(n - 1) @@ -74,7 +77,7 @@ def A_qq_ns(n, sx, L): @nb.njit(cache=True) -def A_hq_ps(n, sx, L): +def A_hq_ps(n, L, cache, is_singlet): r""" |NNLO| heavy-light pure-singlet |OME| :math:`A_{Hq}^{PS,(2)}` given in Eq. (B.1) of :cite:`Buza_1998`. @@ -83,17 +86,19 @@ def A_hq_ps(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_hq_ps : complex |NNLO| heavy-light pure-singlet |OME| :math:`A_{Hq}^{PS,(2)}` """ - S2 = sx[1][0] + S2 = c.get(c.S2, cache, n, is_singlet) F1M = 1.0 / (n - 1.0) * (zeta2 - (S2 - 1.0 / n**2)) F11 = 1.0 / (n + 1.0) * (zeta2 - (S2 + 1.0 / (n + 1.0) ** 2)) @@ -137,7 +142,7 @@ def A_hq_ps(n, sx, L): @nb.njit(cache=True) -def A_hg(n, sx, L): +def A_hg(n, L, cache, is_singlet): r""" |NNLO| heavy-gluon |OME| :math:`A_{Hg}^{S,(2)}` given in Eq. (B.3) of :cite:`Buza_1998`. @@ -147,22 +152,24 @@ def A_hg(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_hg : complex |NNLO| heavy-gluon |OME| :math:`A_{Hg}^{S,(2)}` """ - S1 = sx[0][0] - S2, Sm2 = sx[1] - if len(sx[2]) == 3: - S3, Sm21, Sm3 = sx[2] - else: - S3, _, _, Sm21, _, Sm3 = sx[2] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) S1m = S1 - 1 / n S2m = S2 - 1 / n**2 @@ -277,7 +284,7 @@ def A_hg(n, sx, L): @nb.njit(cache=True) -def A_gq(n, sx, L): +def A_gq(n, L, cache, is_singlet): r""" |NNLO| gluon-quark |OME| :math:`A_{gq,H}^{S,(2)}` given in Eq. (B.5) of :cite:`Buza_1998`. @@ -286,18 +293,20 @@ def A_gq(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_gq : complex |NNLO| gluon-quark |OME| :math:`A_{gq,H}^{S,(2)}` """ - S1 = sx[0][0] - S2 = sx[1][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) S1m = S1 - 1 / n # harmonic_S1(n - 1) B2M = ((S1 - 1.0 / n) ** 2 + S2 - 1.0 / n**2) / (n - 1.0) @@ -327,7 +336,7 @@ def A_gq(n, sx, L): @nb.njit(cache=True) -def A_gg(n, sx, L): +def A_gg(n, L, cache, is_singlet): r""" |NNLO| gluon-gluon |OME| :math:`A_{gg,H}^{S,(2)} ` given in Eq. (B.7) of :cite:`Buza_1998`. @@ -336,17 +345,19 @@ def A_gg(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- A_gg : complex |NNLO| gluon-gluon |OME| :math:`A_{gg,H}^{S,(2)}` """ - S1 = sx[0][0] + S1 = c.get(c.S1, cache, n, is_singlet) S1m = S1 - 1 / n # harmonic_S1(n - 1) D1 = -1.0 / n**2 @@ -414,7 +425,7 @@ def A_gg(n, sx, L): @nb.njit(cache=True) -def A_singlet(n, sx, L, is_msbar=False): +def A_singlet(n, L, cache, is_singlet, is_msbar=False): r""" Computes the |NNLO| singlet |OME|. @@ -429,11 +440,12 @@ def A_singlet(n, sx, L, is_msbar=False): ---------- n : complex Mellin moment - sx : list - harmonic sums cache containing: - [[:math:`S_1,S_{-1}`],[:math:`S_2,S_{-2}`],[:math:`S_3,S_{-2,1},S_{-3}`]] L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise is_msbar: bool add the |MSbar| contribution @@ -450,11 +462,11 @@ def A_singlet(n, sx, L, is_msbar=False): A_gq : :math:`A_{gq, H}^{S,(2)}` A_gg : :math:`A_{gg, H}^{S,(2)}` """ - A_hq_2 = A_hq_ps(n, sx, L) - A_qq_2 = A_qq_ns(n, sx, L) - A_hg_2 = A_hg(n, sx, L) - A_gq_2 = A_gq(n, sx, L) - A_gg_2 = A_gg(n, sx, L) + A_hq_2 = A_hq_ps(n, L, cache, is_singlet) + A_qq_2 = A_qq_ns(n, L, cache, is_singlet) + A_hg_2 = A_hg(n, L, cache, is_singlet) + A_gq_2 = A_gq(n, L, cache, is_singlet) + A_gg_2 = A_gg(n, L, cache, is_singlet) if is_msbar: A_hg_2 -= 2.0 * 4.0 * constants.CF * A_hg_1(n, L=1.0) A_gg_2 -= 2.0 * 4.0 * constants.CF * A_gg_1(L=1.0) @@ -465,7 +477,7 @@ def A_singlet(n, sx, L, is_msbar=False): @nb.njit(cache=True) -def A_ns(n, sx, L): +def A_ns(n, L, cache, is_singlet): r""" Computes the |NNLO| non-singlet |OME|. @@ -479,11 +491,12 @@ def A_ns(n, sx, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache containing: - [[:math:`S_1,S_{-1}`],[:math:`S_2,S_{-2}`],[:math:`S_3,S_{-2,1},S_{-3}`]] L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -494,4 +507,6 @@ def A_ns(n, sx, L): -------- A_qq_ns : :math:`A_{qq,H}^{NS,(2)}` """ - return np.array([[A_qq_ns(n, sx, L), 0.0], [0 + 0j, 0 + 0j]], np.complex_) + return np.array( + [[A_qq_ns(n, L, cache, is_singlet), 0.0], [0 + 0j, 0 + 0j]], np.complex_ + ) diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py index 06feccd38..da92e0e34 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py @@ -69,7 +69,7 @@ @nb.njit(cache=True) -def A_singlet(n, sx_singlet, sx_non_singlet, nf, L): +def A_singlet(n, nf, L, cache): r"""Computes the |N3LO| singlet |OME|. .. math:: @@ -86,21 +86,13 @@ def A_singlet(n, sx_singlet, sx_non_singlet, nf, L): ---------- n : complex Mellin moment - sx_singlet : list - singlet like harmonic sums cache containing: - - .. math :: - [[S_1,S_{-1}], - [S_2,S_{-2}], - [S_{3}, S_{2,1}, S_{2,-1}, S_{-2,1}, S_{-2,-1}, S_{-3}], - [S_{4}, S_{3,1}, S_{2,1,1}, S_{-2,-2}, S_{-3, 1}, S_{-4}],] - - sx_non_singlet: list - same as sx_singlet but now for non-singlet like harmonics nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + Returns ------- @@ -108,15 +100,15 @@ def A_singlet(n, sx_singlet, sx_non_singlet, nf, L): |NNLO| singlet |OME| :math:`A^{S,(3)}(N)` """ - A_hq_3 = A_Hq(n, sx_singlet, nf, L) - A_hg_3 = A_Hg(n, sx_singlet, nf, L) + A_hq_3 = A_Hq(n, nf, L, cache, True) + A_hg_3 = A_Hg(n, nf, L, cache, True) - A_gq_3 = A_gq(n, sx_singlet, nf, L) - A_gg_3 = A_gg(n, sx_singlet, nf, L) + A_gq_3 = A_gq(n, nf, L, cache, True) + A_gg_3 = A_gg(n, nf, L, cache, True) - A_qq_ps_3 = A_qqPS(n, sx_singlet, nf, L) - A_qq_ns_3 = A_qqNS(n, sx_non_singlet, nf, L) - A_qg_3 = A_qg(n, sx_singlet, nf, L) + A_qq_ps_3 = A_qqPS(n, nf, L, cache, True) + A_qq_ns_3 = A_qqNS(n, nf, L, cache, False) + A_qg_3 = A_qg(n, nf, L, cache, True) A_S = np.array( [ @@ -130,7 +122,7 @@ def A_singlet(n, sx_singlet, sx_non_singlet, nf, L): @nb.njit(cache=True) -def A_ns(n, sx_all, nf, L): +def A_ns(n, nf, L, cache, is_singlet): r"""Computes the |N3LO| non-singlet |OME|. .. math:: @@ -146,20 +138,14 @@ def A_ns(n, sx_all, nf, L): ---------- n : complex Mellin moment - sx_all : list - harmonic sums cache containing: - - .. math :: - [[S_1,S_{-1}], - [S_2,S_{-2}], - [S_{3}, S_{2,1}, S_{2,-1}, S_{-2,1}, S_{-2,-1}, S_{-3}], - [S_{4}, S_{3,1}, S_{2,1,1}, S_{-2,-2}, S_{-3, 1}, S_{-4}],], - [S_{5}, S_{-5}] - nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -171,4 +157,6 @@ def A_ns(n, sx_all, nf, L): A_qqNS_3 : :math:`A_{qq,H}^{NS,(3))}` """ - return np.array([[A_qqNS(n, sx_all, nf, L), 0.0], [0 + 0j, 0 + 0j]], np.complex_) + return np.array( + [[A_qqNS(n, nf, L, cache, is_singlet), 0.0], [0 + 0j, 0 + 0j]], np.complex_ + ) diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py index c676ce558..be45b3be3 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py @@ -2,11 +2,12 @@ import numba as nb import numpy as np +from .....harmonics import cache as c from .aHgstfac import A_Hgstfac @nb.njit(cache=True) -def A_Hg(n, sx, nf, L): # pylint: disable=too-many-locals +def A_Hg(n, nf, L, cache, is_singlet): # pylint: disable=too-many-locals r"""Computes the |N3LO| singlet |OME| :math:`A_{Hg}^{S,(3)}(N)`. The expression is presented in :cite:`Bierenbaum:2009mv`. @@ -17,12 +18,14 @@ def A_Hg(n, sx, nf, L): # pylint: disable=too-many-locals ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -35,12 +38,22 @@ def A_Hg(n, sx, nf, L): # pylint: disable=too-many-locals Incomplete part of the |OME|. """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, S21, _, Sm21, _, Sm3 = sx[2] - S4, S31, S211, Sm22, Sm211, Sm31, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm31 = c.get(c.Sm31, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) a_Hg_l0 = ( - A_Hgstfac(n, sx, nf) + A_Hgstfac(n, nf, cache, is_singlet) + (1.0684950250307503 * (2.0 + n + np.power(n, 2))) / (n * (1.0 + n) * (2.0 + n)) + 0.3333333333333333 diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py index 4eed1890a..fb135fb25 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py @@ -2,9 +2,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_Hgstfac(n, sx, nf): +def A_Hgstfac(n, nf, cache, is_singlet): r"""Computes the approximate incomplete part of :math:`A_{Hg}^{S,(3)}(N)` proportional to :math:`T_{F}`. The expression is presented in :cite:`Blumlein:2017wxd` (eq 3.1). @@ -23,18 +25,30 @@ def A_Hgstfac(n, sx, nf): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- complex """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, S21, _, Sm21, _, Sm3 = sx[2] - S4, S31, S211, Sm22, Sm211, Sm31, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm31 = c.get(c.Sm31, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) momentum_conservation_shift = ( -136.47358087568801 / n**3 * nf - 375.88217393160176 / n**2 ) diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py index 6ea2968ec..a52fdfaa9 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py @@ -2,9 +2,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_Hq(n, sx, nf, L): # pylint: disable=too-many-locals +def A_Hq(n, nf, L, cache, is_singlet): # pylint: disable=too-many-locals r"""Computes the |N3LO| singlet |OME| :math:`A_{Hq}^{S,(3)}(N)`. The expression is presented in :cite:`Ablinger_2015` (eq 5.1) and :cite:`Blumlein:2017wxd` (eq 3.1). @@ -22,12 +24,14 @@ def A_Hq(n, sx, nf, L): # pylint: disable=too-many-locals ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -35,11 +39,21 @@ def A_Hq(n, sx, nf, L): # pylint: disable=too-many-locals :math:`A_{Hq}^{S,(3)}(N)` """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, S21, _, Sm21, _, Sm3 = sx[2] - S4, S31, S211, Sm22, Sm211, Sm31, Sm4 = sx[3] - S5, _ = sx[4] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm31 = c.get(c.Sm31, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) + S5 = c.get(c.S5, cache, n, is_singlet) # fit of: # 2^-N * ( H1 + 8.41439832211716) diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py index b4783d83e..ee2e43798 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py @@ -2,11 +2,12 @@ import numba as nb import numpy as np +from .....harmonics import cache as c from .aggTF2 import A_ggTF2 @nb.njit(cache=True) -def A_gg(n, sx, nf, L): # pylint: disable=too-many-locals +def A_gg(n, nf, L, cache, is_singlet): # pylint: disable=too-many-locals r"""Computes the |N3LO| singlet |OME| :math:`A_{gg}^{S,(3)}(N)`. The expression is presented in :cite:`Bierenbaum:2009mv`. @@ -17,12 +18,14 @@ def A_gg(n, sx, nf, L): # pylint: disable=too-many-locals ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -35,13 +38,23 @@ def A_gg(n, sx, nf, L): # pylint: disable=too-many-locals Incomplete part proportional to :math:`T_{F}^2`. """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, S21, _, Sm21, _, Sm3 = sx[2] - S4, S31, S211, Sm22, Sm211, Sm31, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm31 = c.get(c.Sm31, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) a_gg_l0 = ( -0.35616500834358344 - + A_ggTF2(n, sx) + + A_ggTF2(n, cache, is_singlet) + 0.75 * ( (-19.945240467240673 * (1.0 + n + np.power(n, 2))) diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py index 299006dc6..a3856dbb6 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py @@ -1,9 +1,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_ggTF2(n, sx): +def A_ggTF2(n, cache, is_singlet): r"""Computes the approximate incomplete part of :math:`A_{gg}^{S,(3)}(N)` proportional to :math:`T_{F}^2`. The expression is presented in :cite:`Ablinger:2014uka` (eq 4.2). @@ -16,8 +18,10 @@ def A_ggTF2(n, sx): ---------- n : complex Mellin moment - sx : list - harmonic sums cache + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -25,12 +29,12 @@ def A_ggTF2(n, sx): :math:`A_{gg,T_{F}^2}^{S,(3)}(N)` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - S4 = sx[3][0] - S5 = sx[4][0] - S21 = sx[2][1] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S5 = c.get(c.S5, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) # Parametrization of: # 4^(1-n) Binomial[2 n,n] ( # -7 Zeta[3]+ Sum[(4^x(x!)^2(S[1,x] x-1))/((2x)! x^3),{x,1,n}] diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py index 906e08e52..791a9dcb8 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py @@ -1,9 +1,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_gq(n, sx, nf, L): # pylint: disable=too-many-locals +def A_gq(n, nf, L, cache, is_singlet): # pylint: disable=too-many-locals r"""Computes the |N3LO| singlet |OME| :math:`A_{gq}^{S,(3)}(N)`. The expression is presented in :cite:`Ablinger_2014` (eq 6.3). @@ -14,12 +16,14 @@ def A_gq(n, sx, nf, L): # pylint: disable=too-many-locals ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -27,10 +31,23 @@ def A_gq(n, sx, nf, L): # pylint: disable=too-many-locals :math:`A_{gq}^{S,(3)}(N)` """ - S1, Sm1 = sx[0] - S2, Sm2 = sx[1] - S3, S21, S2m1, Sm21, Sm2m1, Sm3 = sx[2] - S4, S31, S211, Sm22, Sm211, Sm31, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + Sm1 = c.get(c.Sm1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + S2m1 = c.get(c.S2m1, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm2m1 = c.get(c.Sm2m1, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm31 = c.get(c.Sm31, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) a_gq_l0 = ( 0.3333333333333333 * ( diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py index 3c2cc36e4..b510fae85 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py @@ -1,9 +1,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_qg(n, sx, nf, L): +def A_qg(n, nf, L, cache, is_singlet): r"""Computes the |N3LO| singlet |OME| :math:`A_{qg}^{S,(3)}(N)`. The expression is presented in :cite:`Bierenbaum:2009mv`. @@ -14,12 +16,14 @@ def A_qg(n, sx, nf, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -27,10 +31,16 @@ def A_qg(n, sx, nf, L): :math:`A_{qg}^{S,(3)}(N)` """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, S21, _, _, _, Sm3 = sx[2] - S4, S31, S211, _, _, _, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + S21 = c.get(c.S21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + S211 = c.get(c.S211, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) a_qg_l0 = 0.3333333333333333 * nf * ( ( -8.547960200246003 diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py index 40be9ed69..7246c1b7f 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py @@ -1,9 +1,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_qqNS(n, sx, nf, L): +def A_qqNS(n, nf, L, cache, is_singlet): r"""Computes the |N3LO| singlet |OME| :math:`A_{qq}^{NS,(3)}(N)`. The expression is presented in :cite:`Bierenbaum:2009mv` and :cite:`Ablinger:2014vwa`. It contains some weight 5 harmonics sums. @@ -21,12 +23,14 @@ def A_qqNS(n, sx, nf, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -34,10 +38,17 @@ def A_qqNS(n, sx, nf, L): :math:`A_{qq}^{NS,(3)}(N)` """ - S1, _ = sx[0] - S2, Sm2 = sx[1] - S3, _, _, Sm21, _, Sm3 = sx[2] - S4, S31, _, Sm22, Sm211, _, Sm4 = sx[3] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + Sm2 = c.get(c.Sm2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) + Sm21 = c.get(c.Sm21, cache, n, is_singlet) + Sm3 = c.get(c.Sm3, cache, n, is_singlet) + S4 = c.get(c.S4, cache, n, is_singlet) + S31 = c.get(c.S31, cache, n, is_singlet) + Sm22 = c.get(c.Sm22, cache, n, is_singlet) + Sm211 = c.get(c.Sm211, cache, n, is_singlet) + Sm4 = c.get(c.Sm4, cache, n, is_singlet) a_qqNS_l0_nf1 = ( 0.3333333333333333 * nf diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py index ee8f136cd..2fd442f84 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py +++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py @@ -1,9 +1,11 @@ import numba as nb import numpy as np +from .....harmonics import cache as c + @nb.njit(cache=True) -def A_qqPS(n, sx, nf, L): +def A_qqPS(n, nf, L, cache, is_singlet): r"""Computes the |N3LO| singlet |OME| :math:`A_{qq}^{PS,(3)}(N)`. The expression is presented in :cite:`Bierenbaum:2009mv`. @@ -14,12 +16,14 @@ def A_qqPS(n, sx, nf, L): ---------- n : complex Mellin moment - sx : list - harmonic sums cache nf : int number of active flavor below the threshold L : float :math:`\ln(\mu_F^2 / m_h^2)` + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise Returns ------- @@ -27,9 +31,9 @@ def A_qqPS(n, sx, nf, L): :math:`A_{qq}^{PS,(3)}(N)` """ - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] + S1 = c.get(c.S1, cache, n, is_singlet) + S2 = c.get(c.S2, cache, n, is_singlet) + S3 = c.get(c.S3, cache, n, is_singlet) a_qqPS_l0 = ( 0.3333333333333333 diff --git a/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py index 53c68e8d6..f8d5417bc 100644 --- a/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py +++ b/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py @@ -4,10 +4,10 @@ @nb.njit(cache=True) -def A_non_singlet(_matching_order, _n, _sx, _nf, _L): +def A_non_singlet(_matching_order, _n, _nf, _L): raise NotImplementedError("Time-like is not yet implemented") @nb.njit(cache=True) -def A_singlet(_matching_order, _n, _sx, _nf, _L, _is_msbar, _sx_ns=None): +def A_singlet(_matching_order, _n, _nf, _L, _is_msbar): raise NotImplementedError("Time-like is not yet implemented") diff --git a/tests/eko/evolution_operator/test_ome.py b/tests/eko/evolution_operator/test_ome.py index d6f3267f1..08c96b52b 100644 --- a/tests/eko/evolution_operator/test_ome.py +++ b/tests/eko/evolution_operator/test_ome.py @@ -14,7 +14,6 @@ from eko.io.runcards import OperatorCard, TheoryCard from eko.io.types import InversionMethod from eko.runner import legacy -from ekore.harmonics import compute_cache from ekore.operator_matrix_elements.unpolarized.space_like import ( A_non_singlet, A_singlet, @@ -31,13 +30,9 @@ def test_build_ome_as(): nf = 3 is_msbar = False for o in [1, 2, 3]: - sx_singlet = compute_cache(N, max_weight_dict[o], True) - sx_ns = sx_singlet - if o == 3: - sx_ns = compute_cache(N, max_weight_dict[o], False) - - aNS = A_non_singlet((o, 0), N, sx_ns, nf, L) - aS = A_singlet((o, 0), N, sx_singlet, nf, L, is_msbar, sx_ns) + + aNS = A_non_singlet((o, 0), N, nf, L) + aS = A_singlet((o, 0), N, nf, L, is_msbar) for a in [aNS, aS]: for method in ["", "expanded", "exact"]: @@ -58,8 +53,8 @@ def test_build_ome_nlo(): is_msbar = False sx = [[1], [1], [1]] nf = 4 - aNSi = A_non_singlet((1, 0), N, sx, nf, L) - aSi = A_singlet((1, 0), N, sx, nf, L, is_msbar) + aNSi = A_non_singlet((1, 0), N, nf, L) + aSi = A_singlet((1, 0), N, nf, L, is_msbar) for a in [aNSi, aSi]: for method in ["", "expanded", "exact"]: dim = len(a[0]) diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 2a95edcb3..6a2744ff3 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -2,7 +2,9 @@ from eko import basis_rotation as br from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet + from eko.beta import beta_qcd_as2, beta_qcd_as3 + from eko.kernels import non_singlet, singlet from eko.scale_variations import Modes, expanded, exponentiated diff --git a/tests/eko/test_beta.py b/tests/eko/test_beta.py index fb089bb62..4621336de 100644 --- a/tests/eko/test_beta.py +++ b/tests/eko/test_beta.py @@ -6,7 +6,7 @@ import pytest from eko import beta -from ekore.harmonics.constants import zeta3 +from eko.constants import zeta3 def _flav_test(function): diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py index 38935ec00..54d74b949 100644 --- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py +++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py @@ -12,14 +12,16 @@ gqg, ) from eko.constants import CA, CF -from ekore.harmonics import compute_cache + +from ekore.harmonics import cache as c NF = 5 def test_quark_number_conservation(): N = 1 - sx_cache = compute_cache(N, 5, False) + + cache = c.reset() # (ns,s) # the exact expression (nf^2 part) has an nonphysical pole at N=1, @@ -29,90 +31,90 @@ def test_quark_number_conservation(): # np.testing.assert_allclose(gamma_nsv(N, NF, sx_cache), 0, rtol=3e-7) # nf^1 part - np.testing.assert_allclose(gnsv.gamma_nss_nf1(N, sx_cache), 0.000400625, atol=2e-6) + np.testing.assert_allclose(gnsv.gamma_nss_nf1(N, cache, None), 0.000400625, atol=2e-6) # (ns,-) # nf^3 part - np.testing.assert_allclose(gnsp.gamma_ns_nf3(N, sx_cache), 0, atol=3e-15) + np.testing.assert_allclose(gnsp.gamma_ns_nf3(N, cache, None), 0, atol=3e-15) # nf^2 part - np.testing.assert_allclose(gnsm.gamma_nsm_nf2(N, sx_cache), 0, atol=3e-13) + np.testing.assert_allclose(gnsm.gamma_nsm_nf2(N, cache, None), 0, atol=3e-13) # nf^1 part - np.testing.assert_allclose(gnsm.gamma_nsm_nf1(N, sx_cache), 0, atol=2e-11) + np.testing.assert_allclose(gnsm.gamma_nsm_nf1(N, cache, None), 0, atol=2e-11) # nf^0 part - np.testing.assert_allclose(gnsm.gamma_nsm_nf0(N, sx_cache), 0, atol=2e-10) + np.testing.assert_allclose(gnsm.gamma_nsm_nf0(N, cache, None), 0, atol=2e-10) # total - np.testing.assert_allclose(gnsm.gamma_nsm(N, NF, sx_cache), 0, atol=1e-10) + np.testing.assert_allclose(gnsm.gamma_nsm(N, NF, cache, None), 0, atol=1e-10) def test_momentum_conservation(): N = 2 - sx_cache = compute_cache(N, 5, True) + cache = c.reset() # nf^3 part np.testing.assert_allclose( - gnsp.gamma_ns_nf3(N, sx_cache) - + gps.gamma_ps_nf3(N, sx_cache) - + ggq.gamma_gq_nf3(N, sx_cache), + gnsp.gamma_ns_nf3(N, cache, None) + + gps.gamma_ps_nf3(N, cache, None) + + ggq.gamma_gq_nf3(N, cache, None), 0, atol=3e-15, ) np.testing.assert_allclose( - ggg.gamma_gg_nf3(N, sx_cache) + gqg.gamma_qg_nf3(N, sx_cache), 0, atol=2e-7 + ggg.gamma_gg_nf3(N, cache, None) + gqg.gamma_qg_nf3(N, cache, None), 0, atol=2e-7 ) # nf^3 part np.testing.assert_allclose( - gnsp.gamma_ns_nf3(N, sx_cache) - + gps.gamma_ps_nf3(N, sx_cache) - + ggq.gamma_gq_nf3(N, sx_cache), + gnsp.gamma_ns_nf3(N, cache, None) + + gps.gamma_ps_nf3(N, cache, None) + + ggq.gamma_gq_nf3(N, cache, None), 0, atol=3e-15, ) np.testing.assert_allclose( - ggg.gamma_gg_nf3(N, sx_cache) + gqg.gamma_qg_nf3(N, sx_cache), 0, atol=2e-7 + ggg.gamma_gg_nf3(N, cache, None) + gqg.gamma_qg_nf3(N, cache, None), 0, atol=2e-7 ) # nf^2 part np.testing.assert_allclose( - gnsp.gamma_nsp_nf2(N, sx_cache) - + gps.gamma_ps_nf2(N, sx_cache) - + ggq.gamma_gq_nf2(N, sx_cache), + gnsp.gamma_nsp_nf2(N, cache, None) + + gps.gamma_ps_nf2(N, cache, None) + + ggq.gamma_gq_nf2(N, cache, None), 0, atol=3e-13, ) np.testing.assert_allclose( - ggg.gamma_gg_nf2(N, sx_cache) + gqg.gamma_qg_nf2(N, sx_cache), + ggg.gamma_gg_nf2(N, cache, None) + gqg.gamma_qg_nf2(N, cache, None), 0, atol=4e-13, ) # nf^1 part np.testing.assert_allclose( - gnsp.gamma_nsp_nf1(N, sx_cache) - + gps.gamma_ps_nf1(N, sx_cache) - + ggq.gamma_gq_nf1(N, sx_cache), + gnsp.gamma_nsp_nf1(N, cache, None) + + gps.gamma_ps_nf1(N, cache, None) + + ggq.gamma_gq_nf1(N, cache, None), 0, ) np.testing.assert_allclose( - ggg.gamma_gg_nf1(N, sx_cache) + gqg.gamma_qg_nf1(N, sx_cache), + ggg.gamma_gg_nf1(N, cache, None) + gqg.gamma_qg_nf1(N, cache, None), 0, atol=5e-11, ) # nf^0 part np.testing.assert_allclose( - gnsp.gamma_nsp_nf0(N, sx_cache) + ggq.gamma_gq_nf0(N, sx_cache), + gnsp.gamma_nsp_nf0(N, cache, None) + ggq.gamma_gq_nf0(N, cache, None), 0, atol=3e-11, ) np.testing.assert_allclose( - ggg.gamma_gg_nf0(N, sx_cache), + ggg.gamma_gg_nf0(N, cache, None), 0, atol=4e-11, ) # total - g_singlet = gamma_singlet(N, NF, sx_cache) + g_singlet = gamma_singlet(N, NF, cache, None) np.testing.assert_allclose( g_singlet[0, 0] + g_singlet[1, 0], 0, @@ -148,18 +150,18 @@ def test_non_singlet_reference_moments(): 2.90857799, ] for N in [3.0, 5.0, 7.0, 9.0, 11.0, 13.0, 15.0, 17.0]: - sx_cache = compute_cache(N, 5, False) + cache = c.reset() idx = int((N - 3) / 2) if N != 17: np.testing.assert_allclose( - gnsm.gamma_nsm(N, NF, sx_cache), nsm_nf4_refs[idx] + gnsm.gamma_nsm(N, NF, cache, None), nsm_nf4_refs[idx] ) np.testing.assert_allclose( - gnsv.gamma_nsv(N, NF, sx_cache), nss_nf4_refs[idx] + nsm_nf4_refs[idx] + gnsv.gamma_nsv(N, NF, cache, None), nss_nf4_refs[idx] + nsm_nf4_refs[idx] ) gamma_nss = ( - gnsv.gamma_nss_nf1(N, sx_cache) * NF - + gnsv.gamma_nss_nf2(N, sx_cache) * NF**2 + gnsv.gamma_nss_nf1(N, cache, None) * NF + + gnsv.gamma_nss_nf2(N, cache, None) * NF**2 ) np.testing.assert_allclose(gamma_nss, nss_nf4_refs[idx], atol=4e-4) @@ -177,9 +179,9 @@ def test_singlet_reference_moments(): 8119.044600816003, ] for N in [2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0]: - sx_cache = compute_cache(N, 5, True) + cache = c.reset() np.testing.assert_allclose( - gnsp.gamma_nsp(N, NF, sx_cache), nsp_nf4_refs[int((N - 2) / 2)] + gnsp.gamma_nsp(N, NF, cache, None), nsp_nf4_refs[int((N - 2) / 2)] ) @@ -249,15 +251,15 @@ def deltaB3(n, sx): diff = [] ref_vals = [] for N in range(10, 51): - sx_cache = compute_cache(N, 5, not bool(N % 2)) - diff.append(gnsp.gamma_nsp_nf2(N, sx_cache) - gnsm.gamma_nsm_nf2(N, sx_cache)) - ref_vals.append(deltaB3(N, sx_cache)) + cache = c.reset() + diff.append(gnsp.gamma_nsp_nf2(N, cache, None) - gnsm.gamma_nsm_nf2(N, cache, None)) + ref_vals.append(deltaB3(N, cache, None)) np.testing.assert_allclose(diff, ref_vals, atol=5e-4) diff = [] ref_vals = [] for N in range(4, 10): - sx_cache = compute_cache(N, 5, not bool(N % 2)) - diff.append(gnsp.gamma_nsp_nf2(N, sx_cache) - gnsm.gamma_nsm_nf2(N, sx_cache)) - ref_vals.append(deltaB3(N, sx_cache)) + cache = c.reset() + diff.append(gnsp.gamma_nsp_nf2(N, cache, None) - gnsm.gamma_nsm_nf2(N, cache, None)) + ref_vals.append(deltaB3(N, cache, None)) np.testing.assert_allclose(diff, ref_vals, atol=2e-2) diff --git a/tests/ekore/harmonics/test_g_functions.py b/tests/ekore/harmonics/test_g_functions.py index b6a895e6e..3d1851ab5 100644 --- a/tests/ekore/harmonics/test_g_functions.py +++ b/tests/ekore/harmonics/test_g_functions.py @@ -4,8 +4,6 @@ from ekore import harmonics as h -zeta3 = h.constants.zeta3 -log2 = h.constants.log2 # Reference values comes form Mathematica, and are diff --git a/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py index eed992f3d..2408f6ec5 100644 --- a/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py +++ b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py @@ -1,9 +1,9 @@ # Test N3LO OME import numpy as np -from ekore.harmonics import compute_cache from ekore.operator_matrix_elements.unpolarized.space_like import as3 from ekore.operator_matrix_elements.unpolarized.space_like.as3 import A_ns, A_qqNS, A_singlet +from ekore.harmonics import cache as c def test_A_3(): @@ -12,15 +12,15 @@ def test_A_3(): for L in logs: N = 1.0 - sx_cache = compute_cache(N, 5, False) - aNSqq3 = A_qqNS(N, sx_cache, nf, L) + cache = c.reset() + aNSqq3 = A_qqNS(N, nf, L, cache, None) # quark number conservation # the accuracy of this test depends directly on the precision of the # fitted part of aNSqq3 np.testing.assert_allclose(aNSqq3, 0.0, atol=5e-5) N = 2.0 - sx_cache = compute_cache(N, 5, True) + cache = c.reset() # reference value comes form Mathematica, Hg is not fully complete # thus the reference value is not 0.0 but 145.14813631128334. # Here we have imposed a small shift in A_Hgstfac such that @@ -28,9 +28,9 @@ def test_A_3(): # The accuracy of this test depends on the approximation of AggTF2. atol = 2e-4 if L == 0 else 2e-3 np.testing.assert_allclose( - as3.A_gg(N, sx_cache, nf, L) - + as3.A_qg(N, sx_cache, nf, L) - + as3.A_Hg(N, sx_cache, nf, L), + as3.A_gg(N, nf, L, cache, None) + + as3.A_qg(N, nf, L, cache, None) + + as3.A_Hg(N, nf, L, cache, None), 0, atol=atol, ) @@ -45,27 +45,26 @@ def test_A_3(): eps = 1e-6 atol = 2e-5 N = 2.0 + eps - sx_cache = compute_cache(N, 5, True) + cache = c.reset() np.testing.assert_allclose( - as3.A_gq(N, sx_cache, nf, L) - + as3.A_qqNS(N, sx_cache, nf, L) - + as3.A_qqPS(N, sx_cache, nf, L) - + as3.A_Hq(N, sx_cache, nf, L), + as3.A_gq(N, nf, L, cache, None) + + as3.A_qqNS(N, nf, L, cache, None) + + as3.A_qqPS(N, nf, L, cache, None) + + as3.A_Hq(N, nf, L, cache, None), 0, atol=atol, ) N = 3 + 2j - sx_cache = compute_cache(np.random.rand(), 5, True) - ns_sx_cache = compute_cache(np.random.rand(), 5, False) - aS3 = A_singlet(N, sx_cache, ns_sx_cache, nf, L) - aNS3 = A_ns(N, ns_sx_cache, nf, L) + cache = c.reset() + aS3 = A_singlet(N, nf, L, cache) + aNS3 = A_ns(N, nf, L, cache, None) assert aNS3.shape == (2, 2) assert aS3.shape == (3, 3) np.testing.assert_allclose(aS3[:, 2], np.zeros(3)) np.testing.assert_allclose(aNS3[1, 1], 0) - np.testing.assert_allclose(aNS3[0, 0], as3.A_qqNS(N, ns_sx_cache, nf, L)) + np.testing.assert_allclose(aNS3[0, 0], as3.A_qqNS(N, nf, L, cache, None)) def test_Blumlein_3(): @@ -151,9 +150,8 @@ def test_Blumlein_3(): for i, N in enumerate([4.0, 6.0, 10.0, 100.0]): idx = i + 1 for L in [0, 10]: - sx_cache = compute_cache(N, 5, True) - ns_sx_cache = compute_cache(N, 5, False) - aS3 = A_singlet(N, sx_cache, ns_sx_cache, nf, L) + cache = c.reset() + aS3 = A_singlet(N, nf, L, cache) np.testing.assert_allclose( aS3[0, 0], ref_val_gg[L][idx] + ref_val_ggTF2[L][idx], rtol=3e-6 @@ -179,11 +177,11 @@ def test_Blumlein_3(): # Here we test the critical parts for idx, N in enumerate([2.0, 4.0, 6.0, 10.0, 100.0]): - sx_cache = compute_cache(N, 5, True) - Aggtf2 = as3.aggTF2.A_ggTF2(N, sx_cache) + cache = c.reset() + Aggtf2 = as3.aggTF2.A_ggTF2(N, cache, None) np.testing.assert_allclose(Aggtf2, ref_val_ggTF2[0][idx], rtol=3e-6) np.testing.assert_allclose( - as3.agg.A_gg(N, sx_cache, nf, L=0) - Aggtf2, + as3.agg.A_gg(N, nf, 0, cache, None) - Aggtf2, ref_val_gg[0][idx], rtol=3e-6, ) @@ -191,9 +189,9 @@ def test_Blumlein_3(): # odd numbers of qqNS ref_qqNS_odd = [-40.94998646588999, -21.598793547423504, 6.966325573931755] for N, ref in zip([3.0, 15.0, 101.0], ref_qqNS_odd): - sx_cache = compute_cache(N, 5, False) + cache = c.reset() np.testing.assert_allclose( - as3.aqqNS.A_qqNS(N, sx_cache, nf, L=0), ref, rtol=1e-4 + as3.aqqNS.A_qqNS(N, nf, 0, cache, None), ref, rtol=1e-4 ) @@ -245,7 +243,7 @@ def test_AHq_asymptotic(): # Ns = [31.,32.,33.,34.,35.,36.,37.,38.,39.] nf = 3 for N, r in zip(Ns, refs): - sx_cache = compute_cache(N, 5, True) + cache = c.reset() np.testing.assert_allclose( - as3.aHq.A_Hq(N, sx_cache, nf, L=0), r, rtol=7e-6, atol=1e-5 + as3.aHq.A_Hq(N, nf, 0, cache, None), r, rtol=7e-6, atol=1e-5 ) diff --git a/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py index 953d797c7..66daa4f88 100644 --- a/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py +++ b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py @@ -1,16 +1,16 @@ # Test NLO OME import numpy as np -from ekore.harmonics import compute_cache from ekore.operator_matrix_elements.unpolarized.space_like.as1 import A_ns, A_singlet +from ekore.harmonics import cache as c def test_A_1_intrinsic(): + cache = c.reset() L = 100.0 N = 2 - sx = compute_cache(N, 2, True) - aS1 = A_singlet(N, sx, L) + aS1 = A_singlet(N, L, cache, None) # heavy quark momentum conservation np.testing.assert_allclose(aS1[0, 2] + aS1[1, 2] + aS1[2, 2], 0.0, atol=1e-10) @@ -20,11 +20,11 @@ def test_A_1_intrinsic(): def test_A_1_shape(): + cache = c.reset() N = 2 L = 3.0 - sx = compute_cache(N, 2, (-1) ** N == 1) - aNS1i = A_ns(N, sx, L) - aS1i = A_singlet(N, sx, L) + aNS1i = A_ns(N, L, cache, None) + aS1i = A_singlet(N, L, cache, None) assert aNS1i.shape == (2, 2) assert aS1i.shape == (3, 3) @@ -46,8 +46,8 @@ def test_Blumlein_1(): for n in range(N_vals): N = 2 * n + 2 - sx = compute_cache(N, 2, True) + cache = c.reset() for L, ref_gg in ref_val_gg.items(): - aS1 = A_singlet(N, sx, L) + aS1 = A_singlet(N, L, cache, None) np.testing.assert_allclose(aS1[0, 0], ref_gg[n], rtol=1e-6) np.testing.assert_allclose(aS1[2, 0], ref_val_Hg[L][n], rtol=3e-6) diff --git a/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py index 648568d1c..46ffab46f 100644 --- a/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py +++ b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py @@ -2,8 +2,9 @@ import numpy as np -from ekore.harmonics import compute_cache, constants +from eko import constants from ekore.operator_matrix_elements.unpolarized.space_like.as2 import A_ns, A_qq_ns, A_singlet +from ekore.harmonics import cache as c def test_A_2(): @@ -11,21 +12,20 @@ def test_A_2(): for L in logs: N = 1 - sx = compute_cache(N, 3, False) - aNSqq2 = A_qq_ns(N, sx, L) + cache = c.reset() + aNSqq2 = A_qq_ns(N, L, cache, None) # quark number conservation np.testing.assert_allclose(aNSqq2, 0.0, atol=2e-11) N = 2 - sx = compute_cache(N, 3, True) - aS2 = A_singlet(N, sx, L) + aS2 = A_singlet(N, L, cache, None) # gluon momentum conservation np.testing.assert_allclose(aS2[0, 0] + aS2[1, 0] + aS2[2, 0], 0.0, atol=2e-6) # quark momentum conservation np.testing.assert_allclose(aS2[0, 1] + aS2[1, 1] + aS2[2, 1], 0.0, atol=1e-11) - aNS2 = A_ns(N, sx, L) + aNS2 = A_ns(N, L, cache, False) assert aNS2.shape == (2, 2) assert aS2.shape == (3, 3) @@ -37,9 +37,9 @@ def test_A_2(): def test_A_2_shape(): N = np.random.rand() L = 3 - sx = compute_cache(N, 3, (-1) ** N == 1) - aNS2 = A_ns(N, sx, L) - aS2 = A_singlet(N, sx, L) + cache = c.reset() + aNS2 = A_ns(N, L, cache, False) + aS2 = A_singlet(N, L, cache, None) assert aNS2.shape == (2, 2) assert aS2.shape == (3, 3) @@ -53,9 +53,9 @@ def test_pegasus_sign(): # reference value come from Pegasus code transalted Mathematica ref_val = -21133.9 N = 2 - sx = compute_cache(N, 3, True) L = 100.0 - aS2 = A_singlet(N, sx, L) + cache = c.reset() + aS2 = A_singlet(N, L, cache, None) np.testing.assert_allclose(aS2[0, 0], ref_val, rtol=4e-5) @@ -117,9 +117,9 @@ def test_Blumlein_2(): ], } for N in range(2, 11): + cache = c.reset() for L, ref_Hg in ref_val_Hg.items(): - sx = compute_cache(N, 3, True) - aS2 = A_singlet(N, sx, L) + aS2 = A_singlet(N, L, cache, None) if N % 2 == 0: idx = int(N / 2 - 1) np.testing.assert_allclose(aS2[0, 0], ref_val_gg[L][idx], rtol=2e-6) @@ -136,11 +136,11 @@ def test_Hg2_pegasus(): L = 0 for N in range(3, 20): - sx = compute_cache(N, 3, True) - S1 = sx[0][0] - S2 = sx[1][0] - S3 = sx[2][0] - aS2 = A_singlet(N, sx, L) + cache = c.reset() + S1 = c.get(c.S1, cache, N, None) + S2 = c.get(c.S2, cache, N, None) + S3 = c.get(c.S3, cache, N, None) + aS2 = A_singlet(N, L, cache, None) E2 = ( 2.0 / N * (constants.zeta3 - S3 + 1.0 / N * (constants.zeta2 - S2 - S1 / N)) @@ -168,7 +168,7 @@ def test_msbar_matching(): for L in logs: N = 2 - sx = compute_cache(N, 3, True) - aS2 = A_singlet(N, sx, L, True) + cache = c.reset() + aS2 = A_singlet(N, L, cache, None, True) # gluon momentum conservation np.testing.assert_allclose(aS2[0, 0] + aS2[1, 0] + aS2[2, 0], 0.0, atol=2e-6)