diff --git a/benchmarks/FF_LHAPDF_bench.py b/benchmarks/FF_LHAPDF_bench.py new file mode 100644 index 000000000..a807b2c05 --- /dev/null +++ b/benchmarks/FF_LHAPDF_bench.py @@ -0,0 +1,146 @@ +r"""Benchmark FFs from LHAPDF""" + +from banana import register + +from eko import interpolation +from ekomark.benchmark.runner import Runner + +register(__file__) + +base_operator = {"ev_op_iterations": 1, "backward_inversion": "exact"} + +base_theory = { + "Qref": 91.2, + "mc": 1.51, + "mb": 4.92, + "mt": 172.5, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "alphas": 0.118000, + "alphaqed": 0.007496, + "FNS": "ZM-VFNS", + "ModEv": "TRN", +} + +FF_sets_lo = [ + "NNFF10_PIm_lo", + "NNFF10_PIp_lo", + "NNFF10_PIsum_lo", + "NNFF10_KAm_lo", + "NNFF10_KAp_lo", + "NNFF10_KAsum_lo", + "NNFF10_PRm_lo", + "NNFF10_PRp_lo", + "NNFF10_PRsum_lo", +] +FF_sets_nlo = [ + "NNFF10_PIm_nlo", + "NNFF10_PIp_nlo", + "NNFF10_PIsum_nlo", + "NNFF10_KAm_nlo", + "NNFF10_KAp_nlo", + "NNFF10_KAsum_nlo", + "NNFF10_PRm_nlo", + "NNFF10_PRp_nlo", + "NNFF10_PRsum_nlo", + "MAPFF10NLOPIm", + "MAPFF10NLOPIp", + "MAPFF10NLOPIsum", + "MAPFF10NLOKAm", + "MAPFF10NLOKAp", + "MAPFF10NLOKAsum", +] +FF_sets_nnlo = [ + "NNFF10_PIm_nnlo", + "NNFF10_PIp_nnlo", + "NNFF10_PIsum_nnlo", + "NNFF10_KAm_nnlo", + "NNFF10_KAp_nnlo", + "NNFF10_KAsum_nnlo", + "NNFF10_PRm_nnlo", + "NNFF10_PRp_nnlo", + "NNFF10_PRsum_nnlo", + "MAPFF10NNLOPIm", + "MAPFF10NNLOPIp", + "MAPFF10NNLOPIsum", + "MAPFF10NNLOKAm", + "MAPFF10NNLOKAp", + "MAPFF10NNLOKAsum", +] + + +class BenchmarkFF(Runner): + external = "LHAPDF" + rotate_to_evolution_basis = True + + def skip_pdfs(self, _theory): + return [ + -6, + 6, + 22, + "ph", + "T35", + "V35", + ] + + def benchmark_lo(self, Q0=10, Q2grid=(10000, 1000000)): + theory_card = { + **base_theory, + "PTO": 0, + "QED": 0, + "Q0": Q0, + "MaxNfPdf": 5, + "MaxNfAs": 5, + } + + operator_card = { + **base_operator, + "Q2grid": list(Q2grid), + "time_like": [True], + "interpolation_xgrid": interpolation.lambertgrid(100, 0.01), + } + + self.run([theory_card], [operator_card], [FF_sets_lo[7]]) + + def benchmark_nlo(self, Q0=1.65, Q2grid=(10000,)): + theory_card = { + **base_theory, + "PTO": 1, + "QED": 0, + "Q0": Q0, + "MaxNfPdf": 5, + "MaxNfAs": 5, + } + + operator_card = { + **base_operator, + "Q2grid": list(Q2grid), + "time_like": [True], + "interpolation_xgrid": interpolation.lambertgrid(5, 0.01), + } + + self.run([theory_card], [operator_card], [FF_sets_nlo[11]]) + + def benchmark_nnlo(self, Q0=10, Q2grid=(10000,)): + theory_card = { + **base_theory, + "PTO": 2, + "QED": 0, + "Q0": Q0, + "MaxNfPdf": 5, + "MaxNfAs": 5, + } + + operator_card = { + **base_operator, + "Q2grid": list(Q2grid), + "time_like": [True], + "interpolation_xgrid": interpolation.lambertgrid(100, 0.01), + } + + self.run([theory_card], [operator_card], [FF_sets_nnlo[10]]) + + +if __name__ == "__main__": + BenchmarkFF().benchmark_nlo() diff --git a/benchmarks/NNPDF_bench.py b/benchmarks/NNPDF_bench.py index 053c6960d..3cd66f309 100644 --- a/benchmarks/NNPDF_bench.py +++ b/benchmarks/NNPDF_bench.py @@ -4,6 +4,7 @@ import numpy as np from banana import register +from eko import interpolation from ekomark.benchmark.runner import Runner register(__file__) @@ -30,6 +31,21 @@ def skip_pdfs(self, _theory): ] +class BenchmarkNNFF(Runner): + external = "LHAPDF" + rotate_to_evolution_basis = True + + def skip_pdfs(self, _theory): + return [ + -6, + 6, + 22, + "ph", + "T35", + "V35", + ] + + base_operator = {"ev_op_iterations": 1, "backward_inversion": "exact"} base_theory = { @@ -90,6 +106,44 @@ def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)): self.run([theory_card], [operator_card], ["NNPDF40_nnlo_as_01180"]) +class BenchmarkNNFF10(BenchmarkNNFF): + def benchmark_lo(self, Q0=1.65, Q2grid=(100,)): + theory_card = {**base_theory, "PTO": 0, "QED": 0, "Q0": Q0, "MaxNfPdf": 5} + + operator_card = { + **base_operator, + "Q2grid": list(Q2grid), + "time_like": [True], + "interpolation_xgrid": interpolation.lambertgrid(50, 0.01), + } + + self.run([theory_card], [operator_card], ["NNFF10_PIm_lo"]) + + def benchmark_nlo(self, Q0=10, Q2grid=(10000,)): + theory_card = {**base_theory, "PTO": 1, "QED": 0, "Q0": Q0, "MaxNfPdf": 5} + + operator_card = { + **base_operator, + "Q2grid": list(Q2grid), + "time_like": [True], + "interpolation_xgrid": interpolation.lambertgrid(50, 0.01), + } + + self.run([theory_card], [operator_card], ["MAPFF10NLOPIp"]) + + def benchmark_nnlo(self, Q0=100, Q2grid=(1000000,)): + theory_card = {**base_theory, "PTO": 2, "QED": 0, "Q0": Q0, "MaxNfPdf": 5} + + operator_card = { + **base_operator, + "Q2grid": list(Q2grid), + "time_like": [True], + "interpolation_xgrid": interpolation.lambertgrid(100, 0.01), + } + + self.run([theory_card], [operator_card], ["MAPFF10NNLOPIp"]) + + if __name__ == "__main__": low2 = 4**2 high2 = 30**2 @@ -98,6 +152,9 @@ def benchmark_nnlo(self, Q0=1.65, Q2grid=(100,)): # nn31.benchmark_nlo(Q0=np.sqrt(low2), Q2grid=[10]) # # test backward # #nn31.benchmark_nlo(Q0=np.sqrt(high2), Q2grid=[low2]) - nn40 = BenchmarkNNPDF40() + # nn40 = BenchmarkNNPDF40() # nn40.benchmark_nnlo(Q2grid=[100]) - nn40.benchmark_nnlo(Q0=np.sqrt(high2), Q2grid=[low2]) + # nn40.benchmark_nnlo(Q0=np.sqrt(high2), Q2grid=[low2]) + nnff10 = BenchmarkNNFF10() + # nnff10.benchmark_lo() + nnff10.benchmark_nnlo() diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index ca1ca6dc8..8c286b7e0 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -158,6 +158,16 @@ def benchmark_plain(self, pto): cartesian_product(th), operators.build(operators.apfel_config), ["ToyLH"] ) + def benchmark_plain_tl(self, pto): + """Plain configuration""" + + th = self.ffns_theory.copy() + th.update({"PTO": [pto]}) + th["ModEv"] = ["EXA"] # TODO for the time one is sufficient + op = operators.apfel_config.copy() + op["time_like"] = [True] + self.run(cartesian_product(th), operators.build(op), ["ToyLH"]) + def benchmark_sv(self, pto, svmode): """Scale Variation""" @@ -188,11 +198,11 @@ def benchmark_sv(self, pto, svmode): if __name__ == "__main__": - - obj = BenchmarkVFNS() - # obj = BenchmarkFFNS() + # obj = BenchmarkVFNS() + obj = BenchmarkFFNS() # obj.benchmark_plain(2) - obj.benchmark_sv(2, "exponentiated") + # obj.benchmark_sv(2, "exponentiated") # obj.benchmark_kthr(2) # obj.benchmark_msbar(2) + obj.benchmark_plain_tl(1) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index 2864bb671..93756514a 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -776,3 +776,26 @@ @article{Albino:2000cp pages = "93--102", year = "2001" } +@article{Mitov:2006wy, + author = "Mitov, Alexander and Moch, Sven-Olaf", + title = "{QCD Corrections to Semi-Inclusive Hadron Production in Electron-Positron Annihilation at Two Loops}", + eprint = "hep-ph/0604160", + archivePrefix = "arXiv", + reportNumber = "DESY-06-043", + doi = "10.1016/j.nuclphysb.2006.05.018", + journal = "Nucl. Phys. B", + volume = "751", + pages = "18--52", + year = "2006" +} +@article{Gluck:1992zx, + author = "Gluck, M. and Reya, E. and Vogt, A.", + title = "{Parton fragmentation into photons beyond the leading order}", + reportNumber = "DO-TH-92-23", + doi = "10.1103/PhysRevD.51.1427", + journal = "Phys. Rev. D", + volume = "48", + pages = "116", + year = "1993", + note = "[Erratum: Phys.Rev.D 51, 1427 (1995)]" +} \ No newline at end of file diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py index 270d434a8..a2172d3cd 100644 --- a/src/eko/evolution_operator/grid.py +++ b/src/eko/evolution_operator/grid.py @@ -106,6 +106,8 @@ def __init__( if order == (1, 0) and method != "iterate-exact": logger.warning("Evolution: In LO we use the exact solution always!") + logger.info(("Time" if configs.time_like else "Space") + "-like evolution") + self.config = config self.q2_grid = mu2grid self.managers = dict( diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py index 3f241519a..df05dfe32 100644 --- a/src/eko/evolution_operator/operator_matrix_element.py +++ b/src/eko/evolution_operator/operator_matrix_element.py @@ -36,7 +36,6 @@ def build_ome(A, matching_order, a_s, backward_method): ------- ome : numpy.ndarray matching operator matrix - """ # to get the inverse one can use this FORM snippet # Symbol a; @@ -123,7 +122,6 @@ def quad_ker( ------- ker : float evaluated integration kernel - """ ker_base = QuadKerBase(u, is_log, logx, mode0) integrand = ker_base.integrand(areas) diff --git a/src/ekomark/benchmark/external/apfel_utils.py b/src/ekomark/benchmark/external/apfel_utils.py index 0fb2303c4..531a32c2c 100644 --- a/src/ekomark/benchmark/external/apfel_utils.py +++ b/src/ekomark/benchmark/external/apfel_utils.py @@ -66,19 +66,18 @@ def compute_apfel_data( # apfel.SetGridParameters(3, 50, 3, 8e-1) # init evolution + apfel.SetTimeLikeEvolution(operators["time_like"]) apfel.InitializeAPFEL() print(f"Loading APFEL took {(time.perf_counter() - apf_start)} s") # Run apf_tabs = {} for q2 in operators["Q2grid"]: - apfel.EvolveAPFEL(theory["Q0"], np.sqrt(q2)) print(f"Executing APFEL took {(time.perf_counter() - apf_start)} s") tab = {} for pid in br.flavor_basis_pids: - if pid in skip_pdfs: continue diff --git a/src/ekomark/benchmark/external/lhapdf_utils.py b/src/ekomark/benchmark/external/lhapdf_utils.py index fa3095d51..4272c24b1 100644 --- a/src/ekomark/benchmark/external/lhapdf_utils.py +++ b/src/ekomark/benchmark/external/lhapdf_utils.py @@ -26,15 +26,14 @@ def compute_LHAPDF_data(operators, pdf, skip_pdfs, rotate_to_evolution_basis=Fal ref : dict output containing: target_xgrid, values """ + # import pdb; pdb.set_trace() - target_xgrid = operators["xgrid"] + target_xgrid = operators["interpolation_xgrid"] out_tabs = {} for q2 in operators["Q2grid"]: - tab = {} for pid in br.flavor_basis_pids: - if pid in skip_pdfs: continue 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/time_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py index 7d353c93e..3c21d3c3a 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py +++ b/src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py @@ -1,21 +1,92 @@ -r"""The unpolarized, time-like Altarelli-Parisi splitting kernels. +"""The unpolarized time-like Altarelli-Parisi splitting kernels.""" -Normalization is given by +import numba as nb +import numpy as np -.. math:: - \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x) +from . import as1, as2, as3 -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): + r"""Compute the tower of the non-singlet anomalous dimensions. + Parameters + ---------- + order : tuple(int,int) + perturbative orders + mode : 10201 | 10101 | 10200 + sector identifier + n : complex + Mellin variable + nf : int + Number of active flavors -@nb.njit(cache=True) -def gamma_ns(_order, _mode, _n, _nf): - raise NotImplementedError("Polarised is not yet implemented") + Returns + ------- + numpy.ndarray + non-singlet anomalous dimensions + + See Also + -------- + ekore.anomalous_dimensions.unpolarized.time_like.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)` + + """ + gamma_ns = np.zeros(order[0], np.complex_) + gamma_ns[0] = as1.gamma_ns(n, None, False) + if order[0] >= 2: + if mode == 10101: + gamma_ns_1 = as2.gamma_nsp(n, nf, None, 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) + else: + raise NotImplementedError("Non-singlet sector is not implemented") + gamma_ns[1] = gamma_ns_1 + if order[0] >= 3: + if mode == 10101: + gamma_ns_2 = as3.gamma_nsp(n, nf, None, False) + elif mode == 10201: + gamma_ns_2 = as3.gamma_nsm(n, nf, None, False) + elif mode == 10200: + gamma_ns_2 = as3.gamma_nsv(n, nf, None, False) + gamma_ns[2] = gamma_ns_2 + return gamma_ns @nb.njit(cache=True) -def gamma_singlet(_order, _n, _nf): - raise NotImplementedError("Polarised is not yet implemented") +def gamma_singlet(order, n, nf): + r"""Compute the tower of the singlet anomalous dimensions' matrices. + + Parameters + ---------- + order : tuple(int,int) + perturbative orders + n : complex + Mellin variable + nf : int + Number of active flavors + + Returns + ------- + numpy.ndarray + singlet anomalous dimensions matrices + + See Also + -------- + ekore.anomalous_dimensions.unpolarized.time_like.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)` + ekore.anomalous_dimensions.unpolarized.time_like.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)` + + """ + gamma_s = np.zeros((order[0], 2, 2), np.complex_) + gamma_s[0] = as1.gamma_singlet(n, nf, None, True) + if order[0] >= 2: + gamma_s[1] = as2.gamma_singlet(n, nf, None) + if order[0] >= 3: + gamma_s[2] = as3.gamma_singlet(n, nf, None, True) + return gamma_s diff --git a/src/ekore/anomalous_dimensions/unpolarized/time_like/as1.py b/src/ekore/anomalous_dimensions/unpolarized/time_like/as1.py new file mode 100644 index 000000000..57db5e424 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/time_like/as1.py @@ -0,0 +1,175 @@ +"""The unpolarized LO Altarelli-Parisi splitting kernels.""" + +import numba as nb +import numpy as np + +from eko import constants + +from ....harmonics import w1 + +# from ....harmonics import cache as c + + +@nb.njit(cache=True) +def gamma_qq(N, cache, is_singlet=None): + r"""Compute the LO quark-quark anomalous dimension. + + Implements Eqn. (B.3) from :cite:`Mitov:2006wy`. + + Parameters + ---------- + N : complex + Mellin moment + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_qq : complex + LO quark-quark anomalous dimension + :math:`\gamma_{qq}^{(0)}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + result = constants.CF * (-3.0 + (4.0 * s1) - 2.0 / (N * (N + 1.0))) + return result + + +@nb.njit(cache=True) +def gamma_qg(N): + r"""Compute the LO quark-gluon anomalous dimension. + + Implements Eqn. (B.4) from :cite:`Mitov:2006wy` + and Eqn. (A1) from :cite:`Gluck:1992zx`. + + Parameters + ---------- + N : complex + Mellin moment + + Returns + ------- + gamma_qg : complex + LO quark-gluon anomalous dimension + :math:`\gamma_{qg}^{(0)}(N)` + + """ + result = -(N**2 + N + 2.0) / (N * (N + 1.0) * (N + 2.0)) + return result + + +@nb.njit(cache=True) +def gamma_gq(N, nf): + r"""Compute the LO gluon-quark anomalous dimension. + + Implements Eqn. (B.5) from :cite:`Mitov:2006wy` + and Eqn. (A1) from :cite:`Gluck:1992zx`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + + Returns + ------- + gamma_qg : complex + LO quark-gluon anomalous dimension + :math:`\gamma_{gq}^{(0)}(N)` + + """ + result = -4.0 * nf * constants.CF * (N**2 + N + 2.0) / (N * (N - 1.0) * (N + 1.0)) + return result + + +@nb.njit(cache=True) +def gamma_gg(N, nf, cache, is_singlet=None): + r"""Compute the LO gluon-gluon anomalous dimension. + + Implements Eqn. (B.6) from :cite:`Mitov:2006wy`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_qq : complex + LO quark-quark anomalous dimension + :math:`\gamma_{gg}^{(0)}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + result = (2.0 * nf - 11.0 * constants.CA) / 3.0 + 4.0 * constants.CA * ( + s1 - 1.0 / (N * (N - 1.0)) - 1.0 / ((N + 1.0) * (N + 2.0)) + ) + return result + + +@nb.njit(cache=True) +def gamma_ns(N, cache, is_singlet=False): + r"""Compute the LO non-singlet anomalous dimension. + + At LO, :math:`\gamma_{ns}^{(0)} = \gamma_{qq}^{(0)}`. + + Parameters + ---------- + N : complex + Mellin moment + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_ns : complex + LO quark-quark anomalous dimension + :math:`\gamma_{ns}^{(0)}(N)` + + """ + return gamma_qq(N, cache, is_singlet) + + +@nb.njit(cache=True) +def gamma_singlet(N, nf, cache, is_singlet=True): + r"""Compute the LO singlet anomalous dimension matrix. + + Implements Eqn. (2.13) from :cite:`Gluck:1992zx`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_singlet : numpy.ndarray + LO singlet anomalous dimension matrix + :math:`\gamma_{s}^{(0)}` + + """ + result = np.array( + [ + [gamma_qq(N, cache, is_singlet), gamma_gq(N, nf)], + [gamma_qg(N), gamma_gg(N, nf, cache, is_singlet)], + ], + np.complex_, + ) + return result diff --git a/src/ekore/anomalous_dimensions/unpolarized/time_like/as2.py b/src/ekore/anomalous_dimensions/unpolarized/time_like/as2.py new file mode 100644 index 000000000..0579b9f62 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/time_like/as2.py @@ -0,0 +1,741 @@ +"""The unpolarized time-like NLO Altarelli-Parisi splitting kernels.""" + +import numba as nb +import numpy as np + +# from eko.constants import zeta2 +from numpy import power as npp + +from eko import constants + +from ....harmonics import w1, w2, w3 +from ....harmonics.constants import zeta2, zeta3 +from ....harmonics.polygamma import cern_polygamma as polygamma + +# from ....harmonics import cache as c + + +@nb.njit(cache=True) +def gamma_nsp(N, nf, cache=None, is_singlet=None): + r"""Compute the NLO non-singlet positive anomalous dimension. + + Implements Eqn. (B.7) from :cite:`Mitov:2006wy`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_nsp : complex + NLO non-singlet positive anomalous dimension + :math:`\gamma_{ns}^{(1)+}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + s2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + s3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + sm1 = w1.Sm1(N, s1, is_singlet) # to be removed + sm2 = w2.Sm2(N, s2, is_singlet) # c.get(c.Sm2, cache, N, is_singlet) + sm3 = w3.Sm3(N, s3, is_singlet) # c.get(c.Sm3, cache, N, is_singlet) + sm21 = w3.Sm21(N, s1, sm1, is_singlet) # c.get(c.Sm21, cache, N, is_singlet) + + if is_singlet == True: + m1tpN = 1 + elif is_singlet == False: + m1tpN = -1 + else: + m1tpn = npp(-1, N) + + nsp1 = ( + constants.CF + * constants.CF + * ( + ( + 8 + + 24 * N + + 32 * npp(N, 2) + - 11 * npp(N, 3) + - 9 * npp(N, 4) + - 9 * npp(N, 5) + - 3 * npp(N, 6) + + 32 * zeta2 * npp(N, 2) + + 112 * zeta2 * npp(N, 3) + + 176 * zeta2 * npp(N, 4) + + 144 * zeta2 * npp(N, 5) + + 48 * zeta2 * npp(N, 6) + + 32 * npp(N, 3) * m1tpN + ) + / (2 * npp(N, 3) * npp(1 + N, 3)) + - 16 * sm3 + + sm2 * ((16) / (N * (1 + N)) - 32 * s1) + - (4 * (2 + 3 * N + 3 * npp(N, 2)) * s2) / (N * (1 + N)) + + s1 + * ( + 16 * s2 + - ( + 8 + * ( + 1 + + 2 * N + + 4 * zeta2 * npp(N, 2) + + 8 * zeta2 * npp(N, 3) + + 4 * zeta2 * npp(N, 4) + ) + ) + / (npp(N, 2) * npp(1 + N, 2)) + ) + - 16 * s3 + + 32 * sm21 + ) + ) + nsp2 = ( + constants.CF + * constants.CA + * ( + ( + 132 + - 208 * N + - 851 * npp(N, 2) + - 757 * npp(N, 3) + - 153 * npp(N, 4) + - 51 * npp(N, 5) + - 144 * npp(N, 2) * m1tpN + ) + / (18 * npp(N, 2) * npp(1 + N, 3)) + + 8 * sm3 + + s1 * (268 / 9) + + sm2 * (16 * s1 - 8 / (N * (1 + N))) + - s2 * (44 / 3) + + 8 * s3 + - 16 * sm21 + ) + ) + nsp3 = ( + nf + * constants.CF + * ( + (-12 + 20 * N + 47 * npp(N, 2) + 6 * npp(N, 3) + 3 * npp(N, 4)) + / (9 * npp(N, 2) * npp(1 + N, 2)) + - s1 * (40 / 9) + + s2 * (8 / 3) + ) + ) + result = nsp1 + nsp2 + nsp3 + return result + + +@nb.njit(cache=True) +def gamma_nsm(N, nf): + r"""Compute the NLO non-singlet negative anomalous dimension. + + Implements Eqn. (B.8) from :cite:`Mitov:2006wy`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_nsm : complex + NLO non-singlet negative anomalous dimension + :math:`\gamma_{ns}^{(1)-}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + s2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + s3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + sm1 = w1.Sm1(N, s1, is_singlet=None) # to be removed + sm2 = w2.Sm2(N, s2, is_singlet=None) # c.get(c.Sm2, cache, N, is_singlet) + sm3 = w3.Sm3(N, s3, is_singlet=None) # c.get(c.Sm3, cache, N, is_singlet) + sm21 = w3.Sm21(N, s1, sm1, is_singlet=None) # c.get(c.Sm21, cache, N, is_singlet) + + m1tpN = -1 + + nsm1 = ( + constants.CF + * constants.CF + * ( + ( + 40 + + 88 * N + + 96 * npp(N, 2) + + 53 * npp(N, 3) + - 9 * npp(N, 4) + - 9 * npp(N, 5) + - 3 * npp(N, 6) + + 32 * zeta2 * npp(N, 2) + + 112 * zeta2 * npp(N, 3) + + 176 * zeta2 * npp(N, 4) + + 144 * zeta2 * npp(N, 5) + + 48 * zeta2 * npp(N, 6) + + 32 * npp(N, 3) * m1tpN + ) + / (2 * npp(N, 3) * npp(1 + N, 3)) + - 16 * sm3 + + sm2 * (16 / (N * (1 + N)) - 32 * s1) + - (4 * (2 + 3 * N + 3 * npp(N, 2)) * s2) / (N * (1 + N)) + + s1 + * ( + -( + ( + 8 + * ( + 1 + + 2 * N + + 4 * zeta2 * npp(N, 2) + + 8 * zeta2 * npp(N, 3) + + 4 * zeta2 * npp(N, 4) + ) + ) + / (npp(N, 2) * npp(1 + N, 2)) + ) + + 16 * s2 + ) + - 16 * s3 + + 32 * sm21 + ) + ) + nsm2 = ( + constants.CF + * constants.CA + * ( + ( + -144 + - 156 * N + - 496 * npp(N, 2) + - 1139 * npp(N, 3) + - 757 * npp(N, 4) + - 153 * npp(N, 5) + - 51 * npp(N, 6) + - 144 * npp(N, 3) * m1tpN + ) + / (18 * npp(N, 3) * npp(1 + N, 3)) + + 8 * sm3 + + s1 * (268 / 9) + + sm2 * (16 * s1 - 8 / (N * (1 + N))) + - s2 * (44 / 3) + + 8 * s3 + - 16 * sm21 + ) + ) + nsm3 = ( + nf + * constants.CF + * ( + (-12 + 20 * N + 47 * npp(N, 2) + 6 * npp(N, 3) + 3 * npp(N, 4)) + / (9 * npp(N, 2) * npp(1 + N, 2)) + - s1 * (40 / 9) + + s2 * (8 / 3) + ) + ) + result = nsm1 + nsm2 + nsm3 + return result + + +@nb.njit(cache=True) +def gamma_qqs(N, nf): + r"""Compute the NLO quark-quark singlet anomalous dimension. + + Implements Eqn. (B.9) from :cite:`Mitov:2006wy`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + + Returns + ------- + gamma_qqs : complex + NLO quark-quark singlet anomalous dimension + :math:`\gamma_{qq}^{(1)s}(N)` + + """ + qqs1 = ( + constants.CF + * nf + * ( + ( + 4 + * ( + 8 + + 44 * N + + 46 * npp(N, 2) + + 21 * npp(N, 3) + + 14 * npp(N, 4) + + 15 * npp(N, 5) + + 10 * npp(N, 6) + + 2 * npp(N, 7) + ) + ) + / ((-1 + N) * npp(N, 3) * npp(1 + N, 3) * npp(2 + N, 2)) + ) + ) + return qqs1 + + +@nb.njit(cache=True) +def gamma_qg(N, nf): + r"""Compute the NLO quark-gluon anomalous dimension. + + Implements Eqn. (B.10) from :cite:`Mitov:2006wy` + and Eqn. (A1) from :cite:`Gluck:1992zx`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_qg : complex + NLO quark-gluon anomalous dimension + :math:`\gamma_{qg}^{(1)}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + s2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + sm2 = w2.Sm2(N, s2, is_singlet=None) # c.get(c.Sm2, cache, N, is_singlet) + + m1tpN = 1 + + qg1 = ( + nf + * constants.CF + * ( + ( + 2 + * ( + 8 + + 12 * N + + 18 * npp(N, 2) + + 77 * npp(N, 3) + + 127 * npp(N, 4) + + 104 * npp(N, 5) + + 45 * npp(N, 6) + + 9 * npp(N, 7) + ) + ) + / (npp(N, 3) * npp(1 + N, 3) * npp(2 + N, 2)) + - ( + 4 + * s1 + * ( + -8 + - 12 * N + + 22 * npp(N, 2) + + 25 * npp(N, 3) + + 10 * npp(N, 4) + + 3 * npp(N, 5) + ) + ) + / (npp(N, 2) * npp(1 + N, 2) * npp(2 + N, 2)) + + (4 * npp(s1, 2) * (2 + N + npp(N, 2))) / (N * (1 + N) * (2 + N)) + - (20 * s2 * (2 + N + npp(N, 2))) / (N * (1 + N) * (2 + N)) + ) + ) + qg2 = ( + nf + * constants.CA + * ( + ( + 4 + * ( + 144 + + 600 * N + + 980 * npp(N, 2) + + 2366 * npp(N, 3) + + 2564 * npp(N, 4) + + 379 * npp(N, 5) + - 1177 * npp(N, 6) + - 1037 * npp(N, 7) + - 423 * npp(N, 8) + - 76 * npp(N, 9) + - 288 * zeta2 * npp(N, 2) + - 720 * zeta2 * npp(N, 3) + - 504 * zeta2 * npp(N, 4) + + 180 * zeta2 * npp(N, 5) + + 576 * zeta2 * npp(N, 6) + + 504 * zeta2 * npp(N, 7) + + 216 * zeta2 * npp(N, 8) + + 36 * zeta2 * npp(N, 9) + - 180 * m1tpN * npp(N, 3) + - 72 * m1tpN * npp(N, 4) + + 108 * m1tpN * npp(N, 5) + + 108 * m1tpN * npp(N, 6) + + 36 * m1tpN * npp(N, 7) + ) + ) + / (9 * (-1 + N) * npp(N, 3) * npp(1 + N, 3) * npp(2 + N, 3)) + + (8 * sm2 * (2 + N + npp(N, 2))) / (N * (1 + N) * (2 + N)) + + ( + 4 + * s1 + * ( + -48 + - 100 * N + + 40 * npp(N, 2) + + 77 * npp(N, 3) + + 32 * npp(N, 4) + + 11 * npp(N, 5) + ) + ) + / (3 * npp(N, 2) * npp(1 + N, 2) * npp(2 + N, 2)) + - (4 * npp(s1, 2) * (2 + N + npp(N, 2))) / (N * (1 + N) * (2 + N)) + + (12 * s2 * (2 + N + npp(N, 2))) / (N * (1 + N) * (2 + N)) + ) + ) + qg3 = ( + nf + * nf + * ( + ( + 8 + * ( + -12 + - 16 * N + + 37 * npp(N, 2) + + 41 * npp(N, 3) + + 17 * npp(N, 4) + + 5 * npp(N, 5) + ) + ) + / (9 * npp(N, 2) * npp(1 + N, 2) * npp(2 + N, 2)) + - (8 * s1 * (2 + N + npp(N, 2))) / (3 * N * (1 + N) * (2 + N)) + ) + ) + result = (1 / (2 * nf)) * (qg1 + qg2 + qg3) + return result + + +@nb.njit(cache=True) +def gamma_gq(N, nf): + r"""Compute the NLO gluon-quark anomalous dimension. + + Implements Eqn. (B.11) from :cite:`Mitov:2006wy` + and Eqn. (A1) from :cite:`Gluck:1992zx`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_gq : complex + NLO gluon-quark anomalous dimension + :math:`\gamma_{gq}^{(1)}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + s2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + sm2 = w2.Sm2(N, s2, True) # c.get(c.Sm2, cache, N, is_singlet) + + m1tpN = 1 + + gq1 = ( + constants.CF + * constants.CF + * ( + ( + 2 + * ( + -4 + - 4 * N + + 41 * npp(N, 2) + + 83 * npp(N, 3) + + 41 * npp(N, 4) + - 11 * npp(N, 5) + - 10 * npp(N, 6) + - 8 * npp(N, 7) + - 16 * zeta2 * npp(N, 2) + - 24 * zeta2 * npp(N, 3) + + 16 * zeta2 * npp(N, 5) + + 16 * zeta2 * npp(N, 6) + + 8 * zeta2 * npp(N, 7) + ) + ) + / (npp(-1 + N, 2) * npp(N, 3) * npp(1 + N, 3)) + + ( + 8 + * (4 - 2 * N - 16 * npp(N, 2) - npp(N, 3) - 2 * npp(N, 4) + npp(N, 5)) + * s1 + ) + / (npp(-1 + N, 2) * npp(N, 2) * npp(1 + N, 2)) + - (4 * (2 + N + npp(N, 2)) * npp(s1, 2)) / ((-1 + N) * N * (1 + N)) + + (12 * (2 + N + npp(N, 2)) * s2) / ((-1 + N) * N * (1 + N)) + ) + ) + gq2 = ( + constants.CF + * constants.CA + * ( + -( + 4 + * ( + 16 + - 144 * npp(N, 2) + - 156 * npp(N, 3) + - 101 * npp(N, 4) + - 77 * npp(N, 5) + - 75 * npp(N, 6) + - 44 * npp(N, 7) + - npp(N, 8) + + 5 * npp(N, 9) + + npp(N, 10) + + 16 * N * m1tpN + + 32 * npp(N, 2) * m1tpN + - 20 * npp(N, 3) * m1tpN + - 44 * npp(N, 4) * m1tpN + - 26 * npp(N, 5) * m1tpN + + 14 * npp(N, 6) * m1tpN + + 22 * npp(N, 7) * m1tpN + + 6 * npp(N, 8) * m1tpN + ) + ) + / (npp(-1 + N, 3) * npp(N, 3) * npp(1 + N, 3) * npp(2 + N, 2)) + + (8 * (2 + N + npp(N, 2)) * sm2) / ((-1 + N) * N * (1 + N)) + - (8 * (2 - 2 * N - 9 * npp(N, 2) + npp(N, 3) - npp(N, 4) + npp(N, 5)) * s1) + / (npp(-1 + N, 2) * npp(N, 2) * npp(1 + N, 2)) + + (4 * (2 + N + npp(N, 2)) * npp(s1, 2)) / ((-1 + N) * N * (1 + N)) + - (20 * (2 + N + npp(N, 2)) * s2) / ((-1 + N) * N * (1 + N)) + ) + ) + result = (2 * nf) * (gq1 + gq2) + return result + + +@nb.njit(cache=True) +def gamma_gg(N, nf): + r"""Compute the NLO gluon-gluon anomalous dimension. + + Implements Eqn. (B.12) from :cite:`Mitov:2006wy`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_gg : complex + NLO gluon-gluon anomalous dimension + :math:`\gamma_{gg}^{(1)}(N)` + + """ + s1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + s2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + s3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + sm1 = w1.Sm1(N, s1, is_singlet=None) # to be removed + sm2 = w2.Sm2(N, s2, is_singlet=None) # c.get(c.Sm2, cache, N, is_singlet) + sm3 = w3.Sm3(N, s3, is_singlet=None) # c.get(c.Sm3, cache, N, is_singlet) + sm21 = w3.Sm21(N, s1, sm1, is_singlet=None) # c.get(c.Sm21, cache, N, is_singlet) + + m1tpN = 1 + + gg1 = ( + nf + * constants.CF + * ( + ( + 2 + * ( + -16 + + 8 * N + + 108 * npp(N, 2) + + 162 * npp(N, 3) + + 106 * npp(N, 4) + + 11 * npp(N, 5) + - 5 * npp(N, 6) + - 2 * npp(N, 7) + + 6 * npp(N, 8) + + 5 * npp(N, 9) + + npp(N, 10) + ) + ) + / (npp(-1 + N, 2) * npp(N, 3) * npp(1 + N, 3) * npp(2 + N, 2)) + ) + ) + gg2 = ( + constants.CA + * constants.CA + * ( + ( + 2 + * ( + -576 + + 240 * N + + 3824 * npp(N, 2) + + 1240 * npp(N, 3) + + 1928 * npp(N, 4) + + 8303 * npp(N, 5) + + 10651 * npp(N, 6) + + 6614 * npp(N, 7) + + 1238 * npp(N, 8) + - 1133 * npp(N, 9) + - 889 * npp(N, 10) + - 288 * npp(N, 11) + - 48 * npp(N, 12) + + 1152 * zeta2 * npp(N, 2) + + 1248 * zeta2 * npp(N, 3) + - 1296 * zeta2 * npp(N, 4) + - 792 * zeta2 * npp(N, 5) + + 876 * zeta2 * npp(N, 6) + - 1368 * zeta2 * npp(N, 7) + - 2340 * zeta2 * npp(N, 8) + + 120 * zeta2 * npp(N, 9) + + 1476 * zeta2 * npp(N, 10) + + 792 * zeta2 * npp(N, 11) + + 132 * zeta2 * npp(N, 12) + - 576 * m1tpN * N + - 1440 * m1tpN * npp(N, 2) + + 216 * m1tpN * npp(N, 3) + + 1800 * m1tpN * npp(N, 4) + + 1800 * m1tpN * npp(N, 5) + - 72 * m1tpN * npp(N, 6) + - 1008 * m1tpN * npp(N, 7) + - 576 * m1tpN * npp(N, 8) + - 144 * m1tpN * npp(N, 9) + ) + ) + / (9 * npp(-1 + N, 3) * npp(N, 3) * npp(1 + N, 3) * npp(2 + N, 3)) + - (8 * sm3) + + (sm2) + * ( + (32 * (1 + N + npp(N, 2))) / ((-1 + N) * N * (1 + N) * (2 + N)) + - 16 * s1 + ) + - (8 * s2 * (12 - 10 * N + npp(N, 2) + 22 * npp(N, 3) + 11 * npp(N, 4))) + / (3 * (-1 + N) * N * (1 + N) * (2 + N)) + + (s1) + * ( + -( + 4 + * ( + -144 + - 144 * N + + 236 * npp(N, 2) + + 308 * npp(N, 3) + + 829 * npp(N, 4) + + 680 * npp(N, 5) + - 134 * npp(N, 6) + - 268 * npp(N, 7) + - 67 * npp(N, 8) + + 288 * zeta2 * npp(N, 2) + + 288 * zeta2 * npp(N, 3) + - 504 * zeta2 * npp(N, 4) + - 576 * zeta2 * npp(N, 5) + + 144 * zeta2 * npp(N, 6) + + 288 * zeta2 * npp(N, 7) + + 72 * zeta2 * npp(N, 8) + ) + ) + / (9 * npp(-1 + N, 2) * npp(N, 2) * npp(1 + N, 2) * npp(2 + N, 2)) + + 16 * s2 + ) + - (8 * s3) + + (16 * sm21) + ) + ) + gg3 = ( + nf + * constants.CA + * ( + -( + 8 + * ( + -12 + + 26 * N + + 132 * npp(N, 2) + + 85 * npp(N, 3) + + 34 * npp(N, 4) + - 9 * npp(N, 5) + - 25 * npp(N, 6) + - 12 * npp(N, 7) + - 3 * npp(N, 8) + + 24 * zeta2 * npp(N, 2) + + 24 * zeta2 * npp(N, 3) + - 42 * zeta2 * npp(N, 4) + - 48 * zeta2 * npp(N, 5) + + 12 * zeta2 * npp(N, 6) + + 24 * zeta2 * npp(N, 7) + + 6 * zeta2 * npp(N, 8) + ) + ) + / (9 * npp(-1 + N, 2) * npp(N, 2) * npp(1 + N, 2) * npp(2 + N, 2)) + - (s1 * (40 / 9)) + + (s2 * (16 / 3)) + ) + ) + result = gg1 + gg2 + gg3 + return result + + +@nb.njit(cache=True) +def gamma_singlet(N, nf, cache): + r"""Compute the NLO singlet anomalous dimension matrix. + + Implements Eqn. (2.13) from :cite:`Gluck:1992zx`. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_singlet : numpy.ndarray + NLO singlet anomalous dimension matrix + :math:`\gamma_{s}^{(1)}` + + """ + gamma_qq = gamma_nsp(N, nf, cache, True) + gamma_qqs(N, nf) + + result = np.array( + [ + [gamma_qq, gamma_gq(N, nf)], + [gamma_qg(N, nf), gamma_gg(N, nf)], + ], + np.complex_, + ) + return result diff --git a/src/ekore/anomalous_dimensions/unpolarized/time_like/as2mela.py b/src/ekore/anomalous_dimensions/unpolarized/time_like/as2mela.py new file mode 100644 index 000000000..daa61b1e7 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/time_like/as2mela.py @@ -0,0 +1,708 @@ +import numba as nb +import numpy as np + +# from eko.constants import zeta2 +from numpy import power as npp + +from eko import constants + +from ....harmonics import w1, w2, w3 +from ....harmonics.constants import zeta2 as ZETA2 +from ....harmonics.constants import zeta3 as ZETA3 +from ....harmonics.polygamma import cern_polygamma as polygamma + +# NS = N * N +# NT = NS * N +# NFO = NT * N +# NFI = NFO * N +# NSI = NFI * N +# NSE = NSI * N +# NE = NSE * N +# NN = NE * N + +# NM2 = N - 2 +# NM = N - 1 +# N1 = N + 1 +# N2 = N + 2 +# NMS = NM * NM +# NMT = NMS * NM +# N1S = N1 * N1 +# N1T = N1S * N1 +# N2S = N2 * N2 +# N2T = N2S * N2 + +# S1 = w1.S1(N) +# S2 = w2.S2(N) + +# N3 = N + 3 +# N4 = N + 4 +# N5 = N + 5 +# N6 = N + 6 + +# S11 = S1 + 1/N1 +# S12 = S11 + 1/N2 +# S13 = S12 + 1/N3 +# S14 = S13 + 1/N4 +# S15 = S14 + 1/N5 +# S16 = S15 + 1/N6 +# S1M = S1 - 1 / N +# S21 = S2 + 1/N1S +# S22 = S21 + 1/N2S +# S2M = S2 - 1 / NS + +# SPMOM = (1.0000 * (ZETA2 - S1 / N ) / N - +# 0.9992 * (ZETA2 - S11/ N1) / N1 + +# 0.9851 * (ZETA2 - S12/ N2) / N2 - +# 0.9005 * (ZETA2 - S13/ N3) / N3 + +# 0.6621 * (ZETA2 - S14/ N4) / N4 - +# 0.3174 * (ZETA2 - S15/ N5) / N5 + +# 0.0699 * (ZETA2 - S16/ N6) / N6 ) + +# SLC = - 5/8 * ZETA3 +# SLV = (- ZETA2/2* (polygamma(N1/2) - polygamma(N/2)) +# + S1/NS + SPMOM) +# SSCHLM = SLC - SLV +# SSTR2M = ZETA2 - polygamma(N1/2,1) +# SSTR3M = 0.5 * polygamma(N1/2,2) + ZETA3 + +# SSCHLP = SLC + SLV +# SSTR2P = ZETA2 - polygamma(N2/2,1) +# SSTR3P = 0.5 * polygamma(N2/2,2) + ZETA3 + +# DS2NM = - polygamma(NM2/2+1,1) + polygamma(NM/2+1,1) +# DS2N = - polygamma(NM/2+1,1) + polygamma(N/2+1,1) +# DS2N1 = - polygamma(N/2+1,1) + polygamma(N1/2+1,1) +# DS2N2 = - polygamma(N1/2+1,1) + polygamma(N2/2+1,1) + +# PNMA = ( 16* S1 * (2* N + 1) / (NS * N1S) + +# 16* (2* S1 - 1/(N * N1)) * ( S2 - SSTR2M ) + +# 64* SSCHLM + 24* S2 - 3 - 8* SSTR3M - +# 8* (3* NT + NS -1) / (NT * N1T) + +# 16* (2* NS + 2* N +1) / (NT * N1T) ) * (-0.5) +# PNPA = ( 16* S1 * (2* N + 1) / (NS * N1S) + +# 16* (2* S1 - 1/(N * N1)) * ( S2 - SSTR2P ) + +# 64* SSCHLP + 24* S2 - 3 - 8* SSTR3P - +# 8* (3* NT + NS -1) / (NT * N1T) - +# 16* (2* NS + 2* N +1)/(NT * N1T) ) * (-0.5) + +# PNSB = ( S1 * (536/9 + 8* (2* N + 1) / (NS * N1S)) - +# (16* S1 + 52/3- 8/(N * N1)) * S2 - 43/6 - +# (151* NFO + 263* NT + 97* NS + 3* N + 9) * +# 4/ (9* NT * N1T) ) * (-0.5) +# PNSC = ( -160/9* S1 + 32/3.* S2 + 4/3 + +# 16*(11*NS+5*N-3)/(9* NS * N1S))*(-0.5) + +# PPSA = ((5* NFI + 32* NFO + 49* NT+38* NS + 28* N + 8) +# / (NM * NT * N1T * N2S) * 2 ) +# PGGA = (- (2* NFI + 5* NFO + 8* NT + 7* NS- 2* N - 2) +# * 8* S1 / (NMS * NS * N1S * N2S) - 67/9* S1 + 8/3 +# - 4* SSTR2P * (NS + N + 1) / (NM * N * N1 * N2) +# + 2* S1 * SSTR2P - 4* SSCHLP + 0.5 * SSTR3P +# + (457* NN + 2742* NE + 6040* NSE + 6098* NSI +# + 1567* NFI - 2344* NFO - 1632* NT + 560* NS +# + 1488* N + 576) / (18* NMS * NT * N1T * N2T)) +# PGGB = ((38* NFO + 76* NT + 94* NS + 56* N + 12) *(-2) +# / (9* NM * NS * N1S * N2) + 20/9* S1 - 4/3) +# PGGC = ((2* NSI + 4* NFI + NFO - 10* NT - 5* NS - 4* N +# - 4) * (-2) / (NM * NT * N1T * N2) - 1) + +# PPSTL = (-40/9 * 1/NM + 4/NT + 10/NS - 16/N +# + 8/N1 + 112/9 * 1/N2 + 18/N1S +# + 4/N1T + 16/3 * 1/N2S) + +# PQQATL = (( -4 * S1 + 3 + 2/(N*N1) ) +# * ( 2*S2 - 2 * ZETA2 - (2*N + 1)/(NS*N1S) )) +# PQQBTL = (-80/9 * 1/NM + 8/NT + 12/NS - 12/N +# + 8/N1T + 28/N1S - 4/N1 + 32/3 * 1/N2S +# + 224/9 * 1/N2) + +# PQGA = (S11 * (NS + N + 2)/(N * N1 * N2) + 1/NS - 5/3 * 1/N +# - 1/(N * N1) - 2/N1S + 4/3 * 1/N1 + 4/N2S +# - 4/3 * 1/N2) +# PQGB = (( - 2 * S11**2 + 2 * S11 + 10 * S21 ) +# * ( NS + N + 2 ) / ( N * N1 * N2 ) +# + 4 * S11 +# * ( -1/NS + 1/N + 1/(N*N1) + 2/N1S - 4/N2S ) +# - 2/NT + 5/NS - 12/N + 4/(NS*N1) - 12/(N*N1S) +# - 6/(N*N1) + 4/N1T - 4/N1S + 23/N1 - 20/N2) +# PQGC = (( 2 * S11**2 - 10/3 * S11 - 6 * S21 +# + 1 * ( DPSI(N2/2,1) - DPSI(N1/2,1) ) - 6 * ZETA2 ) +# * ( NS + N + 2 ) / ( N * N1 * N2 ) +# - 4 * S11 * +# ( -2/NS + 1/N + 1/(N*N1) + 4/N1S - 6/N2S ) +# - 40/9 * 1/NM + 4/NT + 8/3 * 1/NS +# + 26/9 * 1/N - 8/(NS*N1S) + 22/3 * 1/(N*N1) +# + 16/N1T + 68/3 * 1/N1S - 190/9 * 1/N1 +# + 8/(N1S*N2) - 4/N2S + 356/9 * 1/N2) + +# PGQA = (( S1**2 - 3*S2 - 4 * ZETA2 ) +# * ( NS + N + 2 ) / ( NM * N * N1 ) +# + 2 * S1 +# * ( 4/NMS - 2/(NM*N) - 4/NS + 3/N1S - 1/N1 ) +# - 8/(NMS*N) + 8/(NM*NS) + 2/NT + 8/NS - 1/(2*N) +# + 1/N1T - 5/2 * 1/N1S + 9/2 * 1/N1) + +# PGQB = (( -1 * S1**2 + 5 * S2 +# - 0.5 * ( DPSI(N1/2,1) - DPSI(N/2,1) ) + ZETA2 ) +# * ( NS + N + 2 ) / ( NM * N * N1 ) +# + 2 * S1 +# * ( -2/NMS + 2/(NM*N) + 2/NS - 2/N1S + 1/N1 ) +# - 8/NMT + 6/NMS + 17/9 * 1/NM + 4/(NMS*N) +# - 12/(NM*NS) - 8/NS + 5/N - 2/(NS*N1) - 2/N1T +# - 7/N1S - 1/N1 - 8/3 * 1/N2S - 44/9 * 1/N2) + +# PGGATL = (- 16/3 * 1/NMS + 80/9 * 1/NM + 8/NT +# - 16/NS + 12/N + 8/N1T - 24/N1S + 4/N1 +# - 16/3 * 1/N2S - 224/9 * 1/N2) +# PGGBTL = S2 - 1/NMS + 1/NS - 1/N1S + 1/N2S - ZETA2 + +# PGGCTL = (- 8 * S1 * S2 + 8 * S1 +# * ( 1 / NMS - 1 / NS + 1 / N1S - 1 / N2S + ZETA2 ) +# + ( 8 * S2 - 8 * ZETA2 ) +# * ( 1 / NM - 1 / N + 1 / N1 - 1 / N2 + 11 / 12 ) +# - 8 / NMT + 22 / 3 * 1 / NMS - 8 / ( NMS * N ) +# - 8 / ( NM * NS ) - 8 / NT - 14 / 3 * 1 / NS +# - 8 / N1T + 14 / 3 * 1 / N1S - 8 / ( N1S * N2 ) +# - 8 / ( N1 * N2S ) - 8 / N2T - 22 / 3 * 1 / N2S) + +# PNSTL = (( - 4*S1 + 3 + 2/(N*N1) ) +# * ( 2*S2 - 2*ZETA2 - (2*N + 1)/(NS*N1S) )) + + +@nb.njit(cache=True) +def gamma_nsp(N, nf): + NS = N * N + NT = NS * N + NFO = NT * N + NFI = NFO * N + NSI = NFI * N + NSE = NSI * N + NE = NSE * N + NN = NE * N + + NM2 = N - 2 + NM = N - 1 + N1 = N + 1 + N2 = N + 2 + NMS = NM * NM + NMT = NMS * NM + N1S = N1 * N1 + N1T = N1S * N1 + N2S = N2 * N2 + N2T = N2S * N2 + + S1 = w1.S1(N) + S2 = w2.S2(N) + + N3 = N + 3 + N4 = N + 4 + N5 = N + 5 + N6 = N + 6 + + S11 = S1 + 1 / N1 + S12 = S11 + 1 / N2 + S13 = S12 + 1 / N3 + S14 = S13 + 1 / N4 + S15 = S14 + 1 / N5 + S16 = S15 + 1 / N6 + S1M = S1 - 1 / N + S21 = S2 + 1 / N1S + S22 = S21 + 1 / N2S + S2M = S2 - 1 / NS + + SPMOM = ( + 1.0000 * (ZETA2 - S1 / N) / N + - 0.9992 * (ZETA2 - S11 / N1) / N1 + + 0.9851 * (ZETA2 - S12 / N2) / N2 + - 0.9005 * (ZETA2 - S13 / N3) / N3 + + 0.6621 * (ZETA2 - S14 / N4) / N4 + - 0.3174 * (ZETA2 - S15 / N5) / N5 + + 0.0699 * (ZETA2 - S16 / N6) / N6 + ) + + SLC = -5 / 8 * ZETA3 + SLV = -ZETA2 / 2 * (polygamma(N1 / 2, 0) - polygamma(N / 2, 0)) + S1 / NS + SPMOM + SSCHLM = SLC - SLV + SSTR2M = ZETA2 - polygamma(N1 / 2, 1) + SSTR3M = 0.5 * polygamma(N1 / 2, 2) + ZETA3 + + SSCHLP = SLC + SLV + SSTR2P = ZETA2 - polygamma(N2 / 2, 1) + SSTR3P = 0.5 * polygamma(N2 / 2, 2) + ZETA3 + + DS2NM = -polygamma(NM2 / 2 + 1, 1) + polygamma(NM / 2 + 1, 1) + DS2N = -polygamma(NM / 2 + 1, 1) + polygamma(N / 2 + 1, 1) + DS2N1 = -polygamma(N / 2 + 1, 1) + polygamma(N1 / 2 + 1, 1) + DS2N2 = -polygamma(N1 / 2 + 1, 1) + polygamma(N2 / 2 + 1, 1) + + PNPA = ( + 16 * S1 * (2 * N + 1) / (NS * N1S) + + 16 * (2 * S1 - 1 / (N * N1)) * (S2 - SSTR2P) + + 64 * SSCHLP + + 24 * S2 + - 3 + - 8 * SSTR3P + - 8 * (3 * NT + NS - 1) / (NT * N1T) + - 16 * (2 * NS + 2 * N + 1) / (NT * N1T) + ) * (-0.5) + PNSB = ( + S1 * (536 / 9 + 8 * (2 * N + 1) / (NS * N1S)) + - (16 * S1 + 52 / 3 - 8 / (N * N1)) * S2 + - 43 / 6 + - (151 * NFO + 263 * NT + 97 * NS + 3 * N + 9) * 4 / (9 * NT * N1T) + ) * (-0.5) + PNSC = ( + -160 / 9 * S1 + + 32 / 3.0 * S2 + + 4 / 3 + + 16 * (11 * NS + 5 * N - 3) / (9 * NS * N1S) + ) * (-0.5) + PNSTL = (-4 * S1 + 3 + 2 / (N * N1)) * ( + 2 * S2 - 2 * ZETA2 - (2 * N + 1) / (NS * N1S) + ) + + result = ( + constants.CF + * ( + (constants.CF - constants.CA / 2) * PNPA + + constants.CA * PNSB + + (1 / 2) * nf * PNSC + ) + + constants.CF**2 * PNSTL * 4 + ) + return -result + + +@nb.njit(cache=True) +def gamma_nsm(N, nf): + NS = N * N + NT = NS * N + NFO = NT * N + NFI = NFO * N + NSI = NFI * N + NSE = NSI * N + NE = NSE * N + NN = NE * N + + NM2 = N - 2 + NM = N - 1 + N1 = N + 1 + N2 = N + 2 + NMS = NM * NM + NMT = NMS * NM + N1S = N1 * N1 + N1T = N1S * N1 + N2S = N2 * N2 + N2T = N2S * N2 + + S1 = w1.S1(N) + S2 = w2.S2(N) + + N3 = N + 3 + N4 = N + 4 + N5 = N + 5 + N6 = N + 6 + + S11 = S1 + 1 / N1 + S12 = S11 + 1 / N2 + S13 = S12 + 1 / N3 + S14 = S13 + 1 / N4 + S15 = S14 + 1 / N5 + S16 = S15 + 1 / N6 + S1M = S1 - 1 / N + S21 = S2 + 1 / N1S + S22 = S21 + 1 / N2S + S2M = S2 - 1 / NS + + SPMOM = ( + 1.0000 * (ZETA2 - S1 / N) / N + - 0.9992 * (ZETA2 - S11 / N1) / N1 + + 0.9851 * (ZETA2 - S12 / N2) / N2 + - 0.9005 * (ZETA2 - S13 / N3) / N3 + + 0.6621 * (ZETA2 - S14 / N4) / N4 + - 0.3174 * (ZETA2 - S15 / N5) / N5 + + 0.0699 * (ZETA2 - S16 / N6) / N6 + ) + + SLC = -5 / 8 * ZETA3 + SLV = -ZETA2 / 2 * (polygamma(N1 / 2, 0) - polygamma(N / 2, 0)) + S1 / NS + SPMOM + SSCHLM = SLC - SLV + SSTR2M = ZETA2 - polygamma(N1 / 2, 1) + SSTR3M = 0.5 * polygamma(N1 / 2, 2) + ZETA3 + + SSCHLP = SLC + SLV + SSTR2P = ZETA2 - polygamma(N2 / 2, 1) + SSTR3P = 0.5 * polygamma(N2 / 2, 2) + ZETA3 + + DS2NM = -polygamma(NM2 / 2 + 1, 1) + polygamma(NM / 2 + 1, 1) + DS2N = -polygamma(NM / 2 + 1, 1) + polygamma(N / 2 + 1, 1) + DS2N1 = -polygamma(N / 2 + 1, 1) + polygamma(N1 / 2 + 1, 1) + DS2N2 = -polygamma(N1 / 2 + 1, 1) + polygamma(N2 / 2 + 1, 1) + + PNMA = ( + 16 * S1 * (2 * N + 1) / (NS * N1S) + + 16 * (2 * S1 - 1 / (N * N1)) * (S2 - SSTR2M) + + 64 * SSCHLM + + 24 * S2 + - 3 + - 8 * SSTR3M + - 8 * (3 * NT + NS - 1) / (NT * N1T) + + 16 * (2 * NS + 2 * N + 1) / (NT * N1T) + ) * (-0.5) + PNSB = ( + S1 * (536 / 9 + 8 * (2 * N + 1) / (NS * N1S)) + - (16 * S1 + 52 / 3 - 8 / (N * N1)) * S2 + - 43 / 6 + - (151 * NFO + 263 * NT + 97 * NS + 3 * N + 9) * 4 / (9 * NT * N1T) + ) * (-0.5) + PNSC = ( + -160 / 9 * S1 + + 32 / 3.0 * S2 + + 4 / 3 + + 16 * (11 * NS + 5 * N - 3) / (9 * NS * N1S) + ) * (-0.5) + PNSTL = (-4 * S1 + 3 + 2 / (N * N1)) * ( + 2 * S2 - 2 * ZETA2 - (2 * N + 1) / (NS * N1S) + ) + + result = ( + constants.CF + * ( + (constants.CF - constants.CA / 2) * PNMA + + constants.CA * PNSB + + (1 / 2) * nf * PNSC + ) + + constants.CF**2 * PNSTL * 4 + ) + return -result + + +@nb.njit(cache=True) +def gamma_singlet(N, nf): + NS = N * N + NT = NS * N + NFO = NT * N + NFI = NFO * N + NSI = NFI * N + NSE = NSI * N + NE = NSE * N + NN = NE * N + + NM2 = N - 2 + NM = N - 1 + N1 = N + 1 + N2 = N + 2 + NMS = NM * NM + NMT = NMS * NM + N1S = N1 * N1 + N1T = N1S * N1 + N2S = N2 * N2 + N2T = N2S * N2 + + S1 = w1.S1(N) + S2 = w2.S2(N) + + N3 = N + 3 + N4 = N + 4 + N5 = N + 5 + N6 = N + 6 + + S11 = S1 + 1 / N1 + S12 = S11 + 1 / N2 + S13 = S12 + 1 / N3 + S14 = S13 + 1 / N4 + S15 = S14 + 1 / N5 + S16 = S15 + 1 / N6 + S1M = S1 - 1 / N + S21 = S2 + 1 / N1S + S22 = S21 + 1 / N2S + S2M = S2 - 1 / NS + + SPMOM = ( + 1.0000 * (ZETA2 - S1 / N) / N + - 0.9992 * (ZETA2 - S11 / N1) / N1 + + 0.9851 * (ZETA2 - S12 / N2) / N2 + - 0.9005 * (ZETA2 - S13 / N3) / N3 + + 0.6621 * (ZETA2 - S14 / N4) / N4 + - 0.3174 * (ZETA2 - S15 / N5) / N5 + + 0.0699 * (ZETA2 - S16 / N6) / N6 + ) + + SLC = -5 / 8 * ZETA3 + SLV = -ZETA2 / 2 * (polygamma(N1 / 2, 0) - polygamma(N / 2, 0)) + S1 / NS + SPMOM + SSCHLM = SLC - SLV + SSTR2M = ZETA2 - polygamma(N1 / 2, 1) + SSTR3M = 0.5 * polygamma(N1 / 2, 2) + ZETA3 + + SSCHLP = SLC + SLV + SSTR2P = ZETA2 - polygamma(N2 / 2, 1) + SSTR3P = 0.5 * polygamma(N2 / 2, 2) + ZETA3 + + DS2NM = -polygamma(NM2 / 2 + 1, 1) + polygamma(NM / 2 + 1, 1) + DS2N = -polygamma(NM / 2 + 1, 1) + polygamma(N / 2 + 1, 1) + DS2N1 = -polygamma(N / 2 + 1, 1) + polygamma(N1 / 2 + 1, 1) + DS2N2 = -polygamma(N1 / 2 + 1, 1) + polygamma(N2 / 2 + 1, 1) + + PNMA = ( + 16 * S1 * (2 * N + 1) / (NS * N1S) + + 16 * (2 * S1 - 1 / (N * N1)) * (S2 - SSTR2M) + + 64 * SSCHLM + + 24 * S2 + - 3 + - 8 * SSTR3M + - 8 * (3 * NT + NS - 1) / (NT * N1T) + + 16 * (2 * NS + 2 * N + 1) / (NT * N1T) + ) * (-0.5) + PNPA = ( + 16 * S1 * (2 * N + 1) / (NS * N1S) + + 16 * (2 * S1 - 1 / (N * N1)) * (S2 - SSTR2P) + + 64 * SSCHLP + + 24 * S2 + - 3 + - 8 * SSTR3P + - 8 * (3 * NT + NS - 1) / (NT * N1T) + - 16 * (2 * NS + 2 * N + 1) / (NT * N1T) + ) * (-0.5) + + PNSB = ( + S1 * (536 / 9 + 8 * (2 * N + 1) / (NS * N1S)) + - (16 * S1 + 52 / 3 - 8 / (N * N1)) * S2 + - 43 / 6 + - (151 * NFO + 263 * NT + 97 * NS + 3 * N + 9) * 4 / (9 * NT * N1T) + ) * (-0.5) + PNSC = ( + -160 / 9 * S1 + + 32 / 3.0 * S2 + + 4 / 3 + + 16 * (11 * NS + 5 * N - 3) / (9 * NS * N1S) + ) * (-0.5) + + PPSA = ( + (5 * NFI + 32 * NFO + 49 * NT + 38 * NS + 28 * N + 8) + / (NM * NT * N1T * N2S) + * 2 + ) + PGGA = ( + -(2 * NFI + 5 * NFO + 8 * NT + 7 * NS - 2 * N - 2) + * 8 + * S1 + / (NMS * NS * N1S * N2S) + - 67 / 9 * S1 + + 8 / 3 + - 4 * SSTR2P * (NS + N + 1) / (NM * N * N1 * N2) + + 2 * S1 * SSTR2P + - 4 * SSCHLP + + 0.5 * SSTR3P + + ( + 457 * NN + + 2742 * NE + + 6040 * NSE + + 6098 * NSI + + 1567 * NFI + - 2344 * NFO + - 1632 * NT + + 560 * NS + + 1488 * N + + 576 + ) + / (18 * NMS * NT * N1T * N2T) + ) + PGGB = ( + (38 * NFO + 76 * NT + 94 * NS + 56 * N + 12) * (-2) / (9 * NM * NS * N1S * N2) + + 20 / 9 * S1 + - 4 / 3 + ) + PGGC = (2 * NSI + 4 * NFI + NFO - 10 * NT - 5 * NS - 4 * N - 4) * (-2) / ( + NM * NT * N1T * N2 + ) - 1 + + PPSTL = ( + -40 / 9 * 1 / NM + + 4 / NT + + 10 / NS + - 16 / N + + 8 / N1 + + 112 / 9 * 1 / N2 + + 18 / N1S + + 4 / N1T + + 16 / 3 * 1 / N2S + ) + + PQQATL = (-4 * S1 + 3 + 2 / (N * N1)) * ( + 2 * S2 - 2 * ZETA2 - (2 * N + 1) / (NS * N1S) + ) + PQQBTL = ( + -80 / 9 * 1 / NM + + 8 / NT + + 12 / NS + - 12 / N + + 8 / N1T + + 28 / N1S + - 4 / N1 + + 32 / 3 * 1 / N2S + + 224 / 9 * 1 / N2 + ) + + PQGA = ( + S11 * (NS + N + 2) / (N * N1 * N2) + + 1 / NS + - 5 / 3 * 1 / N + - 1 / (N * N1) + - 2 / N1S + + 4 / 3 * 1 / N1 + + 4 / N2S + - 4 / 3 * 1 / N2 + ) + PQGB = ( + (-2 * S11**2 + 2 * S11 + 10 * S21) * (NS + N + 2) / (N * N1 * N2) + + 4 * S11 * (-1 / NS + 1 / N + 1 / (N * N1) + 2 / N1S - 4 / N2S) + - 2 / NT + + 5 / NS + - 12 / N + + 4 / (NS * N1) + - 12 / (N * N1S) + - 6 / (N * N1) + + 4 / N1T + - 4 / N1S + + 23 / N1 + - 20 / N2 + ) + PQGC = ( + ( + 2 * S11**2 + - 10 / 3 * S11 + - 6 * S21 + + 1 * (polygamma(N2 / 2, 1) - polygamma(N1 / 2, 1)) + - 6 * ZETA2 + ) + * (NS + N + 2) + / (N * N1 * N2) + - 4 * S11 * (-2 / NS + 1 / N + 1 / (N * N1) + 4 / N1S - 6 / N2S) + - 40 / 9 * 1 / NM + + 4 / NT + + 8 / 3 * 1 / NS + + 26 / 9 * 1 / N + - 8 / (NS * N1S) + + 22 / 3 * 1 / (N * N1) + + 16 / N1T + + 68 / 3 * 1 / N1S + - 190 / 9 * 1 / N1 + + 8 / (N1S * N2) + - 4 / N2S + + 356 / 9 * 1 / N2 + ) + + PGQA = ( + (S1**2 - 3 * S2 - 4 * ZETA2) * (NS + N + 2) / (NM * N * N1) + + 2 * S1 * (4 / NMS - 2 / (NM * N) - 4 / NS + 3 / N1S - 1 / N1) + - 8 / (NMS * N) + + 8 / (NM * NS) + + 2 / NT + + 8 / NS + - 1 / (2 * N) + + 1 / N1T + - 5 / 2 * 1 / N1S + + 9 / 2 * 1 / N1 + ) + + PGQB = ( + ( + -1 * S1**2 + + 5 * S2 + - 0.5 * (polygamma(N1 / 2, 1) - polygamma(N / 2, 1)) + + ZETA2 + ) + * (NS + N + 2) + / (NM * N * N1) + + 2 * S1 * (-2 / NMS + 2 / (NM * N) + 2 / NS - 2 / N1S + 1 / N1) + - 8 / NMT + + 6 / NMS + + 17 / 9 * 1 / NM + + 4 / (NMS * N) + - 12 / (NM * NS) + - 8 / NS + + 5 / N + - 2 / (NS * N1) + - 2 / N1T + - 7 / N1S + - 1 / N1 + - 8 / 3 * 1 / N2S + - 44 / 9 * 1 / N2 + ) + + PGGATL = ( + -16 / 3 * 1 / NMS + + 80 / 9 * 1 / NM + + 8 / NT + - 16 / NS + + 12 / N + + 8 / N1T + - 24 / N1S + + 4 / N1 + - 16 / 3 * 1 / N2S + - 224 / 9 * 1 / N2 + ) + PGGBTL = S2 - 1 / NMS + 1 / NS - 1 / N1S + 1 / N2S - ZETA2 + + PGGCTL = ( + -8 * S1 * S2 + + 8 * S1 * (1 / NMS - 1 / NS + 1 / N1S - 1 / N2S + ZETA2) + + (8 * S2 - 8 * ZETA2) * (1 / NM - 1 / N + 1 / N1 - 1 / N2 + 11 / 12) + - 8 / NMT + + 22 / 3 * 1 / NMS + - 8 / (NMS * N) + - 8 / (NM * NS) + - 8 / NT + - 14 / 3 * 1 / NS + - 8 / N1T + + 14 / 3 * 1 / N1S + - 8 / (N1S * N2) + - 8 / (N1 * N2S) + - 8 / N2T + - 22 / 3 * 1 / N2S + ) + + PNSTL = (-4 * S1 + 3 + 2 / (N * N1)) * ( + 2 * S2 - 2 * ZETA2 - (2 * N + 1) / (NS * N1S) + ) + + el1 = -( + ( + constants.CF + * ( + (constants.CF - constants.CA / 2) * PNPA + + constants.CA * PNSB + + (1 / 2) * nf * PNSC + ) + + constants.CF**2 * PNSTL * 4 + ) + + (1 / 2) * nf * constants.CF * PPSTL * 4 + ) + el2 = -(constants.CF**2 * PGQA + constants.CF * constants.CA * PGQB) * 4 * 2 * nf + el3 = -( + ( + 8 / 3 * ((1 / 2) * nf) ** 2 * PQGA + + constants.CF * (1 / 2) * nf * PQGB + + constants.CA * (1 / 2) * nf * PQGC + ) + * 4 + / (2 * nf) + ) + el4 = -( + ( + constants.CA * constants.CA * PGGA + + (1 / 2) * nf * (constants.CA * PGGB + constants.CF * PGGC) + ) + * 4 + + constants.CF * (1 / 2) * nf * 4 * PGGATL + + constants.CA * (1 / 2) * nf * 4 * (-8 / 3) * PGGBTL + + constants.CA * constants.CA * 4 * PGGCTL + ) + result = np.array( + [ + [el1, el2], + [el3, el4], + ], + np.complex_, + ) + return result diff --git a/src/ekore/anomalous_dimensions/unpolarized/time_like/as3.py b/src/ekore/anomalous_dimensions/unpolarized/time_like/as3.py new file mode 100644 index 000000000..a36815938 --- /dev/null +++ b/src/ekore/anomalous_dimensions/unpolarized/time_like/as3.py @@ -0,0 +1,832 @@ +"""The unpolarized time-like NNLO Altarelli-Parisi splitting kernels.""" + +import numba as nb +import numpy as np +from numpy import power as npp + +from ....harmonics import w1, w2, w3, w4 +from ....harmonics.constants import zeta2, zeta3 + +# from eko.constants import zeta2, zeta3 + + +# from ....harmonics import cache as c + + +@nb.njit(cache=True) +def gamma_nsp(N, nf, cache, is_singlet=None): + r"""Compute the NNLO non-singlet positive anomalous dimension. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_nsp : complex + NNLO non-singlet positive anomalous dimension + :math:`\gamma_{ns}^{(2)+}(N)` + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S21 = S2 + N1I2 + S31 = S3 + N1I3 + A0 = -S1M + B1 = -S1 * NI + C0 = NI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + D1 = -NI2 + D11 = -N1I2 + D2 = 2 * NI3 + D3 = -6 * NI2 * NI2 + D4 = 24 * NI2 * NI3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + E11 = S11 * N1I2 + (S21 - zeta2) * N1I + E2 = 2 * (-S1 * NI3 + (zeta2 - S2) * NI2 - (S3 - zeta3) * NI) + E21 = 2 * (-S11 * N1I3 + (zeta2 - S21) * N1I2 - (S31 - zeta3) * N1I) + + PP2 = ( + 1174.898 * A0 + + 1295.625 + - 707.67 * B1 + + 593.9 * C3 + - 1075.3 * C2 + - 4249.4 * C1 + + 1658.7 * C0 + + 1327.5 * D1 + - 189.37 * D2 + - 352 / 9 * D3 + + 128 / 81 * D4 + - 56.907 * E1 + - 559.1 * E11 + - 519.37 * E2 + + nf + * ( + -183.187 * A0 + - 173.935 + + 5120 / 81 * B1 + - 31.84 * C3 + + 181.18 * C2 + + 466.29 * C1 + - 198.10 * C0 + - 168.89 * D1 + - 176 / 81 * D2 + + 64 / 27 * D3 + - 50.758 * E1 + + 85.72 * E11 + + 28.551 * E2 + - 23.102 * E21 + - 39.113 * D11 + ) + ) + PF2 = ( + -( + 17 / 72 + - 2 / 27 * S1 + - 10 / 27 * S2 + + 2 / 9 * S3 + - (12 * npp(N, 4) + 2 * npp(N, 3) - 12 * npp(N, 2) - 2 * N + 3) + / (27 * npp(N, 3) * npp(N1, 3)) + ) + * 32 + / 3 + ) + + result = PP2 + npp(nf, 2) * PF2 + return -result + + +@nb.njit(cache=True) +def gamma_nsm(N, nf, cache, is_singlet=None): + r"""Compute the NNLO non-singlet negative anomalous dimension. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_nsm : complex + NNLO non-singlet negative anomalous dimension + :math:`\gamma_{ns}^{(2)-}(N)` + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + NM = N - 1 + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S12 = S11 + N2I + S21 = S2 + N1I2 + S31 = S3 + N1I3 + A0 = -S1M + B1 = -S1 * NI + C0 = NI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + D1 = -NI2 + D11 = -N1I2 + D2 = 2 * NI3 + D3 = -6 * NI2 * NI2 + D31 = -6 * N1I2 * N1I2 + D4 = 24 * NI2 * NI3 + D41 = 24 * N1I2 * N1I3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + E11 = S11 * N1I2 + (S21 - zeta2) * N1I + E2 = 2 * (-S1 * NI3 + (zeta2 - S2) * NI2 - (S3 - zeta3) * NI) + E21 = 2 * (-S11 * N1I3 + (zeta2 - S21) * N1I2 - (S31 - zeta3) * N1I) + + PM2 = ( + 1174.898 * A0 + + 1295.622 + - 707.94 * B1 + + 407.89 * C3 + - 577.42 * C2 + - 4885.7 * C1 + + 1981.3 * C0 + + 1625.5 * D1 + - 38.298 * D2 + - 3072 / 81 * D3 + - 140 / 81 * D4 + + 4563.2 * E1 + - 5140.6 * E11 + + 1905.4 * E2 + + 1969.5 * E21 + - 437.03 * D31 + - 34.683 * D41 + + nf + * ( + -183.187 * A0 + - 173.9376 + + 5120 / 81 * B1 + - 85.786 * C3 + + 209.19 * C2 + + 511.92 * C1 + - 217.84 * C0 + - 188.99 * D1 + - 784 / 81 * D2 + + 128 / 81 * D3 + + 71.428 * E1 + - 23.722 * E11 + + 30.554 * E2 + - 18.975 * E21 + + 92.453 * D11 + ) + ) + PF2 = ( + -( + 17 / 72 + - 2 / 27 * S1 + - 10 / 27 * S2 + + 2 / 9 * S3 + - (12 * npp(N, 4) + 2 * npp(N, 3) - 12 * npp(N, 2) - 2 * N + 3) + / (27 * npp(N, 3) * npp(N1, 3)) + ) + * 32 + / 3 + ) + + result = PM2 + npp(nf, 2) * PF2 + return -result + + +@nb.njit(cache=True) +def gamma_nsv(N, nf, cache, is_singlet=None): + r"""Compute the NNLO non-singlet valence anomalous dimension. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_nsm : complex + NNLO non-singlet valence anomalous dimension + :math:`\gamma_{ns}^{(2)-}(N)` + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + NM = N - 1 + NMI = 1 / NM + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S12 = S11 + N2I + S21 = S2 + N1I2 + S31 = S3 + N1I3 + B1 = -S1 * NI + if abs(N.imag) < 0.00001 and abs(N.real) < 0.00001: + B1M = -zeta2 + else: + B1M = -S1M * NMI + B11 = -S11 * N1I + B12 = -S12 * N2I + C0 = NI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + C4 = 1 / (N + 4) + D1 = -NI2 + D2 = 2 * NI3 + D3 = -6 * NI2 * NI2 + D4 = 24 * NI2 * NI3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + E2 = 2 * (-S1 * NI3 + (zeta2 - S2) * NI2 - (S3 - zeta3) * NI) + + PS2 = ( + -163.9 * (B1M - B1) + - 7.208 * (B11 - B12) + + 4.82 * (C3 - C4) + - 43.12 * (C2 - C3) + + 44.51 * (C1 - C2) + + 151.49 * (C0 - C1) + + 178.04 * D1 + + 6.892 * D2 + - 40 / 27 * (2 * D3 - D4) + - 173.1 * E1 + + 46.18 * E2 + ) + + result = gamma_nsm(N, nf, cache, is_singlet) + nf * PS2 + return -result + + +@nb.njit(cache=True) +def gamma_qq(N, nf, cache, is_singlet=None): + """Compute the NNLO quark-quark anomalous dimension. + + To be added post benchmark. + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + NM = N - 1 + NMI = 1 / NM + NMI2 = NMI * NMI + NMI3 = NMI * NMI2 + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + N2I2 = N2I * N2I + N2I3 = N2I * N2I2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S4 = w4.S4(N) # c.get(c.S4, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S12 = S11 + N2I + S2M = S2 - NI2 + S21 = S2 + N1I2 + S22 = S21 + N2I2 + S31 = S3 + N1I3 + B1 = -S1 * NI + B11 = -S11 * N1I + B2 = (npp(S1, 2) + S2) * NI + B21 = (npp(S11, 2) + S21) * N1I + B3 = -(npp(S1, 3) + 3 * S1 * S2 + 2 * S3) * NI + B31 = -(npp(S11, 3) + 3 * S11 * S21 + 2 * S31) * N1I + C0 = NI + CM = NMI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + C4 = 1 / (N + 4) + C5 = 1 / (N + 5) + D1 = -NI2 + D1M = -NMI2 + D11 = -N1I2 + D12 = -N2I2 + D2 = 2 * NI3 + D2M = 2 * NMI3 + D21 = 2 * N1I3 + D3 = -6 * NI2 * NI2 + D3M = -6 * NMI2 * NMI2 + D31 = -6 * N1I2 * N1I2 + D32 = -6 * N2I2 * N2I2 + D4 = 24 * NI2 * NI3 + D41 = 24 * N1I2 * N1I3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + E11 = S11 * N1I2 + (S21 - zeta2) * N1I + + PS1 = ( + -256 / 9 * (D3M - D3) + - 128 / 9 * (D2M - D2) + + 324.07 * (D1M - D1) + + 479.87 * (CM - C0) + + 9.072 * (D4 - D41) + + 47.322 * (D3 - D31) + + 425.14 * (D2 - D21) + + 656.49 * (D1 - D11) + - 5.926 * (B3 - B31) + - 9.751 * (B2 - B21) + - 8.650 * (B1 - B11) + - 106.65 * (C0 - C1) + - 848.97 * (C1 - C2) + + 368.79 * (C2 - C3) + - 61.284 * (C3 - C4) + + 96.171 * (E1 - E11) + ) + PS2 = ( + -128 / 81 * (CM - C0) + + 0.019122 * (D4 - D41) + - 1.900 * (D3 - D31) + + 9.1682 * (D2 - D21) + + 57.713 * (D1 - D11) + + 1.778 * (B2 - B21) + + 16.611 * (B1 - B11) + + 87.795 * (C0 - C1) + - 57.688 * (C1 - C2) + - 41.827 * (C2 - C3) + + 25.628 * (C3 - C4) + - 7.9934 * (C4 - C5) + - 2.1031 * (E1 - E11) + + 26.294 * (D11 - D12) + - 7.8645 * (D31 - D32) + ) + + result = nf * (PS1 + nf * PS2) + return -result + + +@nb.njit(cache=True) +def gamma_qg(N, nf, cache, is_singlet=None): + r"""Compute the NNLO quark-gluon anomalous dimension. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_qg : complex + NNLO quark-gluon anomalous dimension + :math:`\gamma_{qg}^{(2)}(N)` + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + NM = N - 1 + NMI = 1 / NM + NMI2 = NMI * NMI + NMI3 = NMI * NMI2 + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + N2I2 = N2I * N2I + N2I3 = N2I * N2I2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S4 = w4.S4(N) # c.get(c.S4, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S12 = S11 + N2I + S2M = S2 - NI2 + S21 = S2 + N1I2 + S22 = S21 + N2I2 + S31 = S3 + N1I3 + B1 = -S1 * NI + B11 = -S11 * N1I + B12 = -S12 * N2I + B2 = (npp(S1, 2) + S2) * NI + B21 = (npp(S11, 2) + S21) * N1I + B22 = (npp(S12, 2) + S22) * N2I + B3 = -(npp(S1, 3) + 3 * S1 * S2 + 2 * S3) * NI + B4 = (npp(S1, 4) + 6 * npp(S1, 2) * S2 + 8 * S1 * S3 + 3 * npp(S2, 2) + 6 * S4) * NI + C0 = NI + CM = NMI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + C4 = 1 / (N + 4) + D1 = -NI2 + D1M = -NMI2 + D11 = -N1I2 + D12 = -N2I2 + D2 = 2 * NI3 + D2M = 2 * NMI3 + D21 = 2 * N1I3 + D22 = 2 * N2I3 + D3 = -6 * NI2 * NI2 + D3M = -6 * NMI2 * NMI2 + D31 = -6 * N1I2 * N1I2 + D4 = 24 * NI2 * NI3 + D41 = 24 * N1I2 * N1I3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + E11 = S11 * N1I2 + (S21 - zeta2) * N1I + E12 = S12 * N2I2 + (S22 - zeta2) * N2I + F1 = 2 * NI * (zeta3 + zeta2 * S1 - 0.5 * NI * (S1 * S1 + S2) - S1 * S2 - S3) + + QG1 = ( + -64 * (D3M + D2M) + + 675.83 * D1M + + 1141.7 * CM + + 42.328 * D4 + + 361.28 * D3 + + 1512 * D2 + + 1864 * D1 + + 100 / 27 * B4 + + 350 / 9 * B3 + + 263.07 * B2 + + 693.84 * B1 + + 603.71 * C0 + - 882.48 * C1 + + 4723.2 * C2 + - 4745.8 * C3 + - 175.28 * C4 + - 1809.4 * E1 + - 107.59 * E11 + - 885.5 * D41 + ) + QG2 = ( + -32 / 27 * D2M + - 3.1752 * D1M + - 2.8986 * CM + + 21.569 * D3 + + 255.62 * D2 + + 619.75 * D1 + - 100 / 27 * B3 + - 35.446 * B2 + - 103.609 * B1 + - 113.81 * C0 + + 341.26 * C1 + - 853.35 * C2 + + 492.1 * C3 + + 14.803 * C4 + + 966.96 * E1 + - 709.1 * E11 + - 1.593 * F1 + - 333.8 * D31 + ) + QG3 = ( + ( + 4 * C0 + + 6 * (D1 + B1) + + 3.8696 * (C0 - 2 * C1 + 2 * C2) + + 4 * (D1 - 2 * D11 + 2 * D12 + B1 - 2 * B11 + 2 * B12) + + 3 * (D2 - 2 * D21 + 2 * D22 + B2 - 2 * B21 + 2 * B22) + + 6 * (E1 - 2 * E11 + 2 * E12) + ) + * 4 + / 9 + ) + + result = (QG1 + nf * (QG2 + nf * QG3)) / 2 + return -result + + +@nb.njit(cache=True) +def gamma_gq(N, nf, cache, is_singlet=None): + r"""Compute the NNLO gluon-quark anomalous dimension. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_gq : complex + NNLO gluon-quark anomalous dimension + :math:`\gamma_{gq}^{(2)}(N)` + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + NM = N - 1 + NMI = 1 / NM + NMI2 = NMI * NMI + NMI3 = NMI * NMI2 + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + N2I2 = N2I * N2I + N2I3 = N2I * N2I2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S4 = w4.S4(N) # c.get(c.S4, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S12 = S11 + N2I + S2M = S2 - NI2 + S21 = S2 + N1I2 + S22 = S21 + N2I2 + S31 = S3 + N1I3 + B1 = -S1 * NI + B2 = (npp(S1, 2) + S2) * NI + B3 = -(npp(S1, 3) + 3 * S1 * S2 + 2 * S3) * NI + B4 = (npp(S1, 4) + 6 * npp(S1, 2) * S2 + 8 * S1 * S3 + 3 * npp(S2, 2) + 6 * S4) * NI + C0 = NI + CM = NMI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + C4 = 1 / (N + 4) + D1 = -NI2 + D1M = -NMI2 + D2 = 2 * NI3 + D2M = 2 * NMI3 + D3 = -6 * NI2 * NI2 + D3M = -6 * NMI2 * NMI2 + D31 = -6 * N1I2 * N1I2 + D4 = 24 * NI2 * NI3 + D4M = 24 * NMI2 * NMI3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + F1 = 2 * NI * (zeta3 + zeta2 * S1 - 0.5 * NI * (S1 * S1 + S2) - S1 * S2 - S3) + + GQ0 = ( + 256 * D4M + + 3712 / 3 * D3M + + 1001.89 * D2M + + 4776.5 * D1M + + 5803.7 * CM + - 30.062 * D4 + - 126.38 * D3 + - 0.71252 * D2 + + 4.4136 * D1 + + 400 / 81 * B4 + + 520 / 27 * B3 + - 220.13 * B2 + - 152.6 * B1 + + 272.85 * C0 + - 7188.7 * C1 + + 5693.2 * C2 + + 146.98 * C3 + + 128.19 * C4 + - 1300.6 * E1 + - 71.23 * F1 + + 543.8 * D31 + ) + GQ1 = ( + 1280 / 81 * D3M + + 2912 / 27 * D2M + + 141.93 * D1M + + 6.0041 * CM + - 48.60 * D3 + - 343.1 * D2 + - 492.0 * D1 + + 80 / 81 * B3 + + 1040 / 81 * B2 + - 16.914 * B1 + - 871.3 * C0 + + 790.13 * C1 + - 241.23 * C2 + + 43.252 * C3 + - 4.3465 * D31 + + 55.048 * E1 + ) + + result = 2 * nf * (GQ0 + nf * GQ1) + return -result + + +@nb.njit(cache=True) +def gamma_gg(N, nf, cache, is_singlet=None): + r"""Compute the NNLO gluon-gluon anomalous dimension. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_gg : complex + NNLO gluon-gluon anomalous dimension + :math:`\gamma_{gg}^{(2)}(N)` + + """ + NI = 1 / N + NI2 = NI * NI + NI3 = NI * NI2 + NM = N - 1 + NMI = 1 / NM + NMI2 = NMI * NMI + NMI3 = NMI * NMI2 + N1 = N + 1 + N1I = 1 / N1 + N1I2 = N1I * N1I + N1I3 = N1I * N1I2 + N2 = N + 2 + N2I = 1 / N2 + N2I2 = N2I * N2I + N2I3 = N2I * N2I2 + S1 = w1.S1(N) # c.get(c.S1, cache, N, is_singlet) + S2 = w2.S2(N) # c.get(c.S2, cache, N, is_singlet) + S3 = w3.S3(N) # c.get(c.S3, cache, N, is_singlet) + S4 = w4.S4(N) # c.get(c.S4, cache, N, is_singlet) + S1M = S1 - NI + S11 = S1 + N1I + S12 = S11 + N2I + S2M = S2 - NI2 + S21 = S2 + N1I2 + S22 = S21 + N2I2 + S31 = S3 + N1I3 + A0 = -S1M + B1 = -S1 * NI + C0 = NI + CM = NMI + C1 = N1I + C2 = N2I + C3 = 1 / (N + 3) + C4 = 1 / (N + 4) + D1 = -NI2 + D1M = -NMI2 + D2 = 2 * NI3 + D2M = 2 * NMI3 + D3 = -6 * NI2 * NI2 + D3M = -6 * NMI2 * NMI2 + D31 = -6 * N1I2 * N1I2 + D4 = 24 * NI2 * NI3 + D4M = 24 * NMI2 * NMI3 + E1 = S1 * NI2 + (S2 - zeta2) * NI + E11 = S11 * N1I2 + (S21 - zeta2) * N1I + E2 = 2 * (-S1 * NI3 + (zeta2 - S2) * NI2 - (S3 - zeta3) * NI) + F1 = 2 * NI * (zeta3 + zeta2 * S1 - 0.5 * NI * (S1 * S1 + S2) - S1 * S2 - S3) + + GG0 = ( + 576 * D4M + + 3168 * D3M + + 3651.1 * D2M + + 10233 * D1M + + 14214.4 * CM + + 191.99 * D4 + + 3281.7 * D3 + + 13528 * D2 + + 12258 * D1 + - 28489 * C0 + + 7469 * C1 + + 30421 * C2 + - 53017 * C3 + + 19556 * C4 + - 186.4 * E1 + - 21328 * E2 + + 5685.8 * D31 + - 3590.1 * B1 + + 4425.451 + + 2643.521 * A0 + ) + GG1 = ( + 448 / 9 * D3M + + 2368 / 9 * D2M + - 5.470 * D1M + - 804.13 * CM + + 18.085 * D4 + + 155.10 * D3 + + 482.94 * D2 + + 4.9934 * D1 + + 248.95 * C0 + + 260.6 * C1 + + 272.79 * C2 + + 2133.2 * C3 + - 926.87 * C4 + + 1266.5 * E1 + - 29.709 * E2 + + 87.771 * F1 + + 485.18 * D31 + + 319.97 * B1 + - 528.719 + - 412.172 * A0 + ) + GG2 = ( + 32 / 27 * D2M + + 368 / 81 * D1M + + 472 / 243 * CM + - 5.0372 * D3 + - 44.80 * D2 + - 69.712 * D1 + - 77.190 * C0 + + 153.27 * C1 + - 106.03 * C2 + + 11.995 * C3 + - 115.01 * E1 + + 96.522 * E11 + - 62.908 * E2 + + 6.4628 + - 16 / 9 * A0 + ) + + result = GG0 + nf * (GG1 + nf * GG2) + return -result + + +@nb.njit(cache=True) +def gamma_singlet(N, nf, cache, is_singlet=True): + r"""Compute the NNLO singlet anomalous dimension matrix. + + Parameters + ---------- + N : complex + Mellin moment + nf : int + No. of active flavors + cache : numpy.ndarray + Harmonic sum cache + is_singlet : boolean + True for singlet, False for non-singlet, None otherwise + + Returns + ------- + gamma_singlet : numpy.ndarray + NNLO singlet anomalous dimension matrix + :math:`\gamma_{s}^{(2)}` + + """ + result = np.array( + [ + [gamma_qq(N, nf, cache, is_singlet), gamma_gq(N, nf, cache, is_singlet)], + [gamma_qg(N, nf, cache, is_singlet), gamma_gg(N, nf, cache, is_singlet)], + ], + np.complex_, + ) + return result 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..81b9f3de8 --- /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, _sx, _nf, _L): + raise NotImplementedError("Polarised, time-like is not yet implemented") + + +@nb.njit(cache=True) +def A_singlet(_matching_order, _n, _sx, _nf, _L, _is_msbar, _sx_ns=None): + raise NotImplementedError("Polarised, time-like is not yet implemented") diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index 2a95edcb3..c0f6d4c2b 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -1,10 +1,10 @@ import numpy as np 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 +from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet def test_modes(): diff --git a/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_as2.py b/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_as2.py new file mode 100644 index 000000000..46559588b --- /dev/null +++ b/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_as2.py @@ -0,0 +1,161 @@ +import numpy as np + +import ekore.anomalous_dimensions.unpolarized.time_like.as2 as as2 +import ekore.anomalous_dimensions.unpolarized.time_like.as2mela as as2m +import ekore.harmonics as h +from eko import constants + +ca = 3 +cf = 4 / 3 +nf = 5 + +# def test_nsp(): +# N = complex(1.0, 0.0) +# # val = [3.11111, 18.0753, 51.3206, 64.5425, 69.9012] +# # for i, j in enumerate[1, 2, 5, 7, 8]: +# # np.testing.assert_almost_equal(as2.gamma_nsp(j, nf), val[i]) +# np.testing.assert_almost_equal(as2.gamma_nsp(N, nf), 3.11111) + +############################# mathematica based unit test ############################ +# def test_nsp(): +# val = [ +# 14.914763950248524, +# 27.170573816296294, +# 36.20843001367968, +# 43.6337052707661, +# 49.77278324218862, +# 55.10505809829998, +# 59.755406196991686, +# 63.922130588519025, +# ] +# for i, j in enumerate(val): +# np.testing.assert_almost_equal(as2.gamma_nsp(i + 2, nf, None), j, decimal=6) + + +# def test_nsm(): +# val = [ +# 14.437397695104497, +# 27.008022375967073, +# 36.13331890256858, +# 43.592750538255814, +# 49.7479719067197, +# 55.08888138952125, +# 59.74427031831462, +# 63.914136685150666, +# ] +# for i, j in enumerate(val): +# np.testing.assert_almost_equal(as2.gamma_nsm(i + 2, nf), j, decimal=6) + + +# def test_qqs(): +# val = [ +# 1280 / 81, +# 4391 / 810, +# 11867 / 4050, +# 186064 / 99225, +# 36490 / 27783, +# 3898835 / 4000752, +# 922333 / 1224720, +# 3969344 / 6615675, +# ] +# for i, j in enumerate(val): +# np.testing.assert_almost_equal(as2.gamma_qqs(i + 2, nf), j, decimal=6) + + +# def test_qg(): +# val = [ +# 4.394040436034569, +# 1.2283097867056794, +# 1.224271622535062, +# 0.5682223219230416, +# 0.5146142100999999, +# 0.20949896137930538, +# 0.15630406923958376, +# -0.015229173167235112, +# ] +# for i, j in enumerate(val): +# np.testing.assert_almost_equal(as2.gamma_qg(i + 2, nf), j, decimal=6) + + +# def test_gq(): +# val = [ +# -307.1723308605102, +# -382.1730799366584, +# -319.66481897017945, +# -226.2889349719677, +# -188.00245135408588, +# -148.32874785657552, +# -127.67168089632983, +# -106.46381918972175, +# ] +# for i, j in enumerate(val): +# np.testing.assert_almost_equal(as2.gamma_gq(i + 2, nf, None), j, decimal=6) + + +# def test_gg(): +# val = [-43.94040436034575, +# -63.52992705212047, +# -49.340293336279544, +# -21.931061116073508, +# -6.6250317132958685, +# 10.14906155587218, +# 21.959513930070727, +# 34.02767982616993] +# for i, j in enumerate(val): +# np.testing.assert_almost_equal(as2.gamma_gg(i + 2, nf), j, decimal=4) + + +def test_nsp(): + for N in 3, 5, 7, 9: + for nf in 3, 4, 5: + np.testing.assert_almost_equal( + as2.gamma_nsp(N, nf, None, None), as2m.gamma_nsp(N, nf), decimal=4 + ) + + +def test_nsm(): + for N in 3, 5, 7, 9: + for nf in 3, 4, 5: + np.testing.assert_almost_equal( + as2.gamma_nsm(N, nf), as2m.gamma_nsm(N, nf), decimal=5 + ) + + +def notest_qq(): + for N in 2, 4, 6, 8: + for nf in 3, 4, 5: + np.testing.assert_almost_equal( + as2.gamma_singlet(N, nf, None)[0, 0], + as2m.gamma_singlet(N, nf)[0, 0], + decimal=5, + ) + + +def test_gq(): + for N in 2, 3, 4, 5, 6, 7, 8, 9: + for nf in 3, 4, 5: + np.testing.assert_almost_equal( + as2.gamma_singlet(N, nf, None)[0, 1], + as2m.gamma_singlet(N, nf)[0, 1], + decimal=6, + ) + + +def test_qg(): + for N in 2, 4, 6, 8: + for nf in 3, 4, 5: + np.testing.assert_almost_equal( + as2.gamma_singlet(N, nf, None)[1, 0], + as2m.gamma_singlet(N, nf)[1, 0], + decimal=5, + ) + + +def test_gg(): + for N in 2, 4, 6, 8: + for nf in 3, 4, 5: + np.testing.assert_almost_equal( + as2.gamma_singlet(N, nf, None)[1, 1], + as2m.gamma_singlet(N, nf)[1, 1], + decimal=4, + ) diff --git a/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_init.py index 9cece76e5..028879bf0 100644 --- a/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_init.py +++ b/tests/ekore/anomalous_dimensions/unpolarized/time_like/test_init.py @@ -2,9 +2,10 @@ import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut +# tests to be added later -def test_init(): - with pytest.raises(NotImplementedError): - ad_ut.gamma_ns((1, 0), 0, 1.0, 4) - with pytest.raises(NotImplementedError): - ad_ut.gamma_singlet((1, 0), 1.0, 4) +# def test_init(): +# with pytest.raises(NotImplementedError): +# ad_ut.gamma_ns((1, 0), 0, 1.0, 4) +# with pytest.raises(NotImplementedError): +# ad_ut.gamma_singlet((1, 0), 1.0, 4)