diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 3aa029801..0d9082599 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -2,7 +2,7 @@
 # See https://pre-commit.com/hooks.html for more hooks
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.2.0
+    rev: v4.3.0
     hooks:
       - id: trailing-whitespace
       - id: end-of-file-fixer
@@ -12,7 +12,7 @@ repos:
       - id: debug-statements
       - id: fix-encoding-pragma
   - repo: https://github.com/psf/black
-    rev: 22.3.0
+    rev: 22.6.0
     hooks:
       - id: black
   - repo: https://github.com/asottile/blacken-docs
@@ -25,6 +25,13 @@ repos:
       - id: isort
         args: ["--profile", "black"]
   - repo: https://github.com/asottile/pyupgrade
-    rev: v2.32.1
+    rev: v2.37.1
     hooks:
       - id: pyupgrade
+  - repo: https://github.com/pycqa/pydocstyle
+    rev: 6.1.1
+    hooks:
+    -   id: pydocstyle
+        files: ^src/
+        additional_dependencies:
+          - toml
diff --git a/doc/source/refs.bib b/doc/source/refs.bib
index 903de2abf..2864bb671 100644
--- a/doc/source/refs.bib
+++ b/doc/source/refs.bib
@@ -722,3 +722,57 @@ @misc{Rottoli
   author        = "Rottoli, Luca",
   howpublished  = "private communication"
 }
+
+@article{Bonvini:2018xvt,
+    author = "Bonvini, Marco and Marzani, Simone",
+    title = "{Four-loop splitting functions at small $x$}",
+    eprint = "1805.06460",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    doi = "10.1007/JHEP06(2018)145",
+    journal = "JHEP",
+    volume = "06",
+    pages = "145",
+    year = "2018"
+}
+
+@article{Moch:2021qrk,
+    author = "Moch, S. and Ruijl, B. and Ueda, T. and Vermaseren, J. A. M. and Vogt, A.",
+    title = "{Low moments of the four-loop splitting functions in QCD}",
+    eprint = "2111.15561",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY 21-203, NIKHEF 21-030, LTH 1282",
+    doi = "10.1016/j.physletb.2021.136853",
+    journal = "Phys. Lett. B",
+    volume = "825",
+    pages = "136853",
+    year = "2022"
+}
+
+@article{Soar:2009yh,
+    author = "Soar, G. and Moch, S. and Vermaseren, J. A. M. and Vogt, A.",
+    title = "{On Higgs-exchange DIS, physical evolution kernels and fourth-order splitting functions at large x}",
+    eprint = "0912.0369",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "LTH-857, DESY-09-211, SFB-CPP-09-119, NIKHEF-09-031",
+    doi = "10.1016/j.nuclphysb.2010.02.003",
+    journal = "Nucl. Phys. B",
+    volume = "832",
+    pages = "152--227",
+    year = "2010"
+}
+
+@article{Albino:2000cp,
+    author = "Albino, Simon and Ball, Richard D.",
+    title = "{Soft resummation of quark anomalous dimensions and coefficient functions in MS-bar factorization}",
+    eprint = "hep-ph/0011133",
+    archivePrefix = "arXiv",
+    reportNumber = "CERN-TH-2000-332, EDINBURGH-2000-23",
+    doi = "10.1016/S0370-2693(01)00742-0",
+    journal = "Phys. Lett. B",
+    volume = "513",
+    pages = "93--102",
+    year = "2001"
+}
diff --git a/doc/source/shared/abbreviations.rst b/doc/source/shared/abbreviations.rst
index fa2a013e0..5ae7c01a8 100644
--- a/doc/source/shared/abbreviations.rst
+++ b/doc/source/shared/abbreviations.rst
@@ -30,6 +30,15 @@
 .. |N3LO| replace::
    :abbr:`N3LO (Next-to-Next-to-Next-to-Leading Order)`
 
+.. |LL| replace::
+   :abbr:`LL (Leading Log)`
+
+.. |NLL| replace::
+   :abbr:`NLL (Next-to-Leading Log)`
+
+.. |NNLL| replace::
+   :abbr:`NNLL (Next-to-Next-to-Leading Log)`
+
 .. Names
 
 .. |DGLAP| replace::
@@ -60,6 +69,11 @@
 .. |QED| replace::
    :abbr:`QED (Quantum Electrodynamics)`
 
+.. |DIS| replace::
+   :abbr:`DIS (Deep Inelastic Scattering)`
+
+.. |BFKL| replace::
+   :abbr:`BFKL (Balitsky-Fadin-Kuraev-Lipatov)`
 .. external
 
 .. |yadism| replace::
diff --git a/doc/source/theory/N3LO_ad.rst b/doc/source/theory/N3LO_ad.rst
index f734d9aa5..d90863b05 100644
--- a/doc/source/theory/N3LO_ad.rst
+++ b/doc/source/theory/N3LO_ad.rst
@@ -2,23 +2,23 @@ N3LO Anomalous Dimensions
 =========================
 
 The |N3LO| |QCD| anomalous dimensions :math:`\gamma^{(3)}` are not yet fully known,
-since they rely on the calculation of 4-loop DIS integrals.
-Moreover the analytical structure of these function is already known to be complicated
+since they rely on the calculation of 4-loop |DIS| integrals.
+Moreover, the analytical structure of these function is already known to be complicated
 since in Mellin space it will included harmonics sum up to weight 7, for which an
-analytical contribution is not available.
+analytical expression is not available.
 
-Here we describe the various assumptions and limits used in order to reconstruct a parametrization
+Here, we describe the various assumptions and limits used in order to reconstruct a parametrization
 that can approximate their contribution.
 In particular we will take advantage of some known physical constrain,
 such as large-x limit, small-x limit, and sum rules, in order to make our reconstruction reasonable.
 
-Generally we remark that the large-x limit correspond to large-N in Mellin space
+Generally, we remark that the large-x limit correspond to large-N in Mellin space
 where the leading contribution comes from the harmonics :math:`S_1(N)`,
 while the small-x region corresponds to poles at :math:`N=0,1` depending on the type of
 divergence.
 
 In any case |N3LO| |DGLAP| evolution at small-x, especially for singlet-like PDFs, will not be reliable
-until the splitting function resummation will not be available up to NNLL.
+until the splitting function resummation will not be available up to |NNLL|.
 
 Non-singlet sector
 ------------------
@@ -33,6 +33,7 @@ In particular, making explicitly the dependence on :math:`n_f`, the non-singlet
 the following terms:
 
     .. list-table:: non-singlet 4-loop Anomalous Dimensions
+        :align: center
         :header-rows: 1
 
         *   -
@@ -80,7 +81,7 @@ In |EKO| they are implemented as follows:
             This part contains the so called double logarithms:
 
             .. math ::
-                \ln(x)^k \quad k=1,..,6, \quad \mathcal{M}[\ln^k(x)] = \frac{1}{N^{(k+1)}}
+                \ln^k(x), \quad \mathcal{M}[\ln^k(x)] = \frac{1}{N^{k+1}}, \quad k=1,..,6
 
             Note the expressions are evaluated with the exact values of the |QCD|
             Casimir invariants, to better agree with the :cite:`Moch:2017uml` parametrization.
@@ -96,7 +97,7 @@ In |EKO| they are implemented as follows:
             :cite:`Henn:2019swt`, while the :math:`B_4` is fixed by the integral of the 4-loop splitting function.
             :math:`C_4,D_4` instead can be computed directly from lower order splitting functions.
             From large-x resummation :cite:`Davies:2016jie`, it is possible to infer further constrains
-            on sub-leading terms :math:`\frac{\ln(N)^k}{N^2}`, since the non-singlet splitting
+            on sub-leading terms :math:`\frac{\ln^k(N)}{N^2}`, since the non-singlet splitting
             functions contain only terms :math:`(1-x)^a\ln^k(1-x)` with :math:`a \ge 1`.
 
         -   The 8 lowest odd or even N moments provided in :cite:`Moch:2017uml`, where
@@ -106,7 +107,8 @@ In |EKO| they are implemented as follows:
         -   The difference between the known moments and the known limits is parametrized
             in Mellin space. The basis includes:
 
-            .. list-table::
+            .. list-table:: :math:`\gamma_{ns,\pm}^{(3)}` parametrization basis
+                :align: center
                 :header-rows: 1
 
                 *   - x-space
@@ -114,30 +116,27 @@ In |EKO| they are implemented as follows:
                 *   - :math:`\delta(1-x)`
                     - 1
                 *   - :math:`(1-x)\ln(1-x)`
-                    - :math:`\mathcal{M}[(1-x)\ln(1-x)] \approx \frac{S_1(N)}{N^2}`
+                    - :math:`\mathcal{M}[(1-x)\ln(1-x)]`
                 *   - :math:`(1-x)\ln^2(1-x)`
-                    - :math:`\mathcal{M}[(1-x)\ln^2(1-x)] \approx \frac{S_1^2(N)}{N^2}`
+                    - :math:`\mathcal{M}[(1-x)\ln^2(1-x)]`
                 *   - :math:`(1-x)\ln^3(1-x)`
-                    - :math:`\mathcal{M}[(1-x)\ln^3(1-x)] \approx \frac{S_1^3(N)}{N^2}`
-                *   - :math:`- Li_2(x) + \zeta_2`
+                    - :math:`\mathcal{M}[(1-x)\ln^3(1-x)]`
+                *   - :math:`- \rm{Li_2}(x) + \zeta_2`
                     - :math:`\frac{S_1(N)}{N^2}`
-
-            which model the sub-leading differences in the :math:`N\to \infty` limit, and:
-
-            .. list-table::
-                :header-rows: 1
-
-                *   - x-space
-                    - N-space
                 *   - :math:`x\ln(x)`
                     - :math:`\frac{1}{(N+1)^2}`
                 *   - :math:`\frac{x}{2}\ln^2(x)`
                     - :math:`\frac{1}{(N+1)^3}`
+                *   - :math:`x^{2}, x^{3}`
+                    - :math:`\frac{1}{(N-2)},\frac{1}{(N-3)}`
 
+            The first five functions model the sub-leading differences in the :math:`N\to \infty` limit,
+            while the last three help the convergence in the small-N region. Finally, we add a polynomial part
+            :math:`x^{2}` or :math:`x^{3}` respectively for :math:`\gamma_{ns,+},\gamma_{ns,-}`.
+            For large-N we have the limit:
 
-            to help the convergence in the small-N region. Finally we add a polynomial part
-            :math:`x^{2(3)}` which corresponds to simple poles at :math:`N=-2,-3`
-            respectively for :math:`\gamma_{ns,+},\gamma_{ns,-}`.
+                .. math ::
+                    \mathcal{M}[(1-x)\ln^k(1-x)] \approx \frac{S_1^k(N)}{N^2}
 
             Note that the constant coefficient is included in the fit, following the procedure done
             in :cite:`Moch:2017uml` (section 4), to achieve a better accuracy.
@@ -147,7 +146,7 @@ Singlet sector
 --------------
 
 In the singlet sector we construct a parametrization for
-:math:`\gamma_{gg}^{(3)},\gamma_{gq}^{(3)},\gamma_{qq}^{(3)},\gamma_{qg}^{(3)}` where:
+:math:`\gamma_{gg}^{(3)},\gamma_{gq}^{(3)},\gamma_{qg}^{(3)},\gamma_{qq}^{(3)}` where:
 
     .. math ::
         \gamma_{qq}^{(3)} = \gamma_{ns,+}^{(3)} + \gamma_{qq,ps}^{(3)}
@@ -156,6 +155,7 @@ In particular, making explicitly the dependence on :math:`n_f`, the singlet anom
 the following terms:
 
     .. list-table:: singlet 4-loop Anomalous Dimensions
+        :align: center
         :header-rows: 1
 
         *   -
@@ -172,23 +172,123 @@ the following terms:
             - |T|
 
         *   - :math:`\gamma_{gq}^{(3)}`
-            -
+            - |T|
             - |T|
             - |T|
             - |T|
 
-        *   - :math:`\gamma_{qq,ps}^{(3)}`
+        *   - :math:`\gamma_{qg}^{(3)}`
             -
             - |T|
             - |T|
             - |T|
 
-        *   - :math:`\gamma_{qg}^{(3)}`
-            - |T|
+        *   - :math:`\gamma_{qq,ps}^{(3)}`
+            -
             - |T|
             - |T|
             - |T|
 
 Only the parts proportional to :math:`n_f^3` are known analytically
 :cite:`Davies:2016jie` and have been included so far.
-The rest will be approximated using some known limits.
+The other parts are approximated using some known limits:
+
+* The remaining contributions include the following constrains.
+
+    *   The small-x limit, given in the large :math:`N_c` approximation by
+        :cite:`Davies:2022ofz` (see Eq. 5.9, 5.10, 5.11, 5.12) and coming
+        from small-x resummation of double-logarithms which fix the leading terms
+        for the pole at :math:`N=0`:
+
+            .. math ::
+                \ln^k(x), \quad \mathcal{M}[\ln^k(x)] = \frac{1}{N^{k+1}}, \quad k=4,5,6
+
+    *   The small-x limit, coming from |BFKL| resummation
+        :cite:`Bonvini:2018xvt` (see Eq. 2.32, 2.20b, 2.21a, 2.21b)
+        which fix the leading terms (|LL|, |NLL|) for the pole at :math:`N=1`:
+
+            .. math ::
+                \frac{\ln^k(x)}{x}, \quad \mathcal{M}[\frac{\ln^k(x)}{x}] = \frac{1}{(N-1)^{k+1}}, \quad k=4,5
+
+        Note that in principle also the term :math:`\frac{\ln^6(x)}{x}` could be present at |N3LO|,
+        but they are vanishing.
+        These terms are way larger than the previous ones in the small-x limit and
+        are effectively determining the raise of the splitting functions at small-x.
+        In particular only the expansion for :math:`\gamma_{gg}^{(3)}` is known at |NLL|.
+        |LL| terms respect the representation symmetry :
+
+            .. math ::
+                \gamma_{gq} & \approx \frac{C_F}{C_A} \gamma_{gg}  \\
+                \gamma_{qq,ps} & \approx \frac{C_F}{C_A} \gamma_{qg} \\
+
+
+    *   The large-x limit of the singlet splitting function is different for the diagonal part
+        and the off-diagonal.
+        It is known that :cite:`Albino:2000cp,Moch:2021qrk` the diagonal terms diverge in N-space as:
+
+            .. math ::
+                \gamma_{kk} \approx A_4 S_1(N)  + \mathcal{O}(1)
+
+        Where again the coefficient :math:`A_4` is the |QCD| cusp anomalous dimension. However, :math:`\gamma_{qq,ps}^{(3)}`
+        do not constrain any divergence at large-x or constant term so its expansion will start as
+        :math:`\mathcal{O}(\frac{1}{N^2})`.
+        The off-diagonal do not contain any +-distributions or delta distributions but can include divergent logarithms
+        of the type :cite:`Soar:2009yh`:
+
+            .. math ::
+                \ln^k(1-x) \quad k=1,..,6
+
+        where also in this case the term :math:`k=6` vanish. The values of the coefficient for :math:`k=4,5`
+        can be guessed from the lower order splitting functions. These logarithms are not present in the diagonal
+        splitting function, which can include at most term :math:`(1-x)\ln^4(1-x)`.
+
+
+    *   The 4 lowest even N moments provided in :cite:`Moch:2021qrk`, where we can use momentum conservation
+        to fix:
+
+            .. math ::
+                & \gamma_{qg}(2) + \gamma_{gg}(2) = 0 \\
+                & \gamma_{qq}(2) + \gamma_{gq}(2) = 0 \\
+
+    *   Finally difference between the known moments and the known limits is parametrized
+        in Mellin space. The basis used in this approximation is different for each splitting
+        function as listed in the following tables.
+
+
+        .. list-table::  :math:`\gamma_{gg}^{(3)}` parametrization basis
+            :align: center
+
+            *   - :math:`\frac{1}{(N-1)^2}`
+                - :math:`\frac{1}{(N-1)}`
+                - :math:`1`
+                - :math:`\mathcal{M}[\ln(1-x)](N)`
+
+        .. list-table::  :math:`\gamma_{gq}^{(3)}` parametrization basis
+            :align: center
+
+            *   - :math:`\frac{1}{(N-1)^3}`
+                - :math:`\frac{1}{(N-1)^2}`
+                - :math:`\mathcal{M}[\ln^3(1-x)](N)`
+                - :math:`\mathcal{M}[(1-x)\ln^3(1-x)](N)`
+
+        .. list-table::  :math:`\gamma_{qg}^{(3)}` parametrization basis
+            :align: center
+
+            *   - :math:`\frac{1}{(N-1)^2}`
+                - :math:`\frac{1}{(N-1)}`
+                - :math:`\mathcal{M}[\ln^3(1-x)](N)`
+                - :math:`\mathcal{M}[(1-x)\ln^3(1-x)](N)`
+
+        .. list-table::  :math:`\gamma_{qq,ps}^{(3)}` parametrization basis
+            :align: center
+
+            *   - :math:`\frac{1}{(N-1)^2} - \frac{1}{N^2}`
+                - :math:`\frac{1}{(N-1)} - \frac{1}{N}`
+                - :math:`\frac{S_1^3(N)}{N^2}`
+                - :math:`\frac{S_1^2(N)}{N^2}`
+
+        Note that for :math:`\gamma_{qq,ps},\gamma_{qg}` the parts proportional
+        to :math:`n_f^0` are not present.
+        Furthermore for the part :math:`\propto n_f^2` in :math:`\gamma_{gq}^{(3)}`
+        we adopt a slightly different basis to account fot the fact that the leading
+        contribution for the pole at :math:`N=1` is :math:`\frac{1}{(N-1)^2}`.
diff --git a/extras/n3lo_ad.tar.gz b/extras/n3lo_ad.tar.gz
index 7b0cbe595..0cc138fb7 100644
Binary files a/extras/n3lo_ad.tar.gz and b/extras/n3lo_ad.tar.gz differ
diff --git a/src/eko/anomalous_dimensions/as4/ggg.py b/src/eko/anomalous_dimensions/as4/ggg.py
index 939cda60d..c5b2d661a 100644
--- a/src/eko/anomalous_dimensions/as4/ggg.py
+++ b/src/eko/anomalous_dimensions/as4/ggg.py
@@ -1,15 +1,16 @@
 # -*- coding: utf-8 -*-
-"""
-This module contains the anomalous dimension :math:`\\gamma_{gg}^{(3)}`
+r"""This module contains the anomalous dimension :math:`\gamma_{gg}^{(3)}`
 """
 import numba as nb
 import numpy as np
 
+from ...harmonics.log_functions import lm11
+
 
 @nb.njit(cache=True)
 def gamma_gg_nf3(n, sx):
-    """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`.
+    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`.
 
     Parameters
     ----------
@@ -21,7 +22,7 @@ def gamma_gg_nf3(n, sx):
     Returns
     -------
     complex
-            |N3LO| non-singlet anomalous dimension :math:`\\gamma_{gg}^{(3)}|_{nf^3}`
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^3}`
 
     """
     S1 = sx[0][0]
@@ -137,23 +138,101 @@ def gamma_gg_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf2(_n, _sx):
-    return 0
+def gamma_gg_nf1(n, sx):
+    r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{gg}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^1}`
+
+    """
+    S1 = sx[0][0]
+    return (
+        18371.82290926215
+        + 1992.766087237516 / np.power(-1.0 + n, 3)
+        - 10452.576219113906 / np.power(-1.0 + n, 2)
+        + 19915.749052258507 / (-1.0 + n)
+        + 20005.925925925927 / np.power(n, 7)
+        - 19449.679012345678 / np.power(n, 6)
+        + 80274.123066115 / np.power(n, 5)
+        - 11714.245609287387 * S1
+        + 18617.18939338288 * lm11(n, S1)
+    )
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf1(_n, _sx):
-    return 0
+def gamma_gg_nf2(n, sx):
+    r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{gg}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^2}`
+
+    """
+    S1 = sx[0][0]
+    return (
+        -436.18166675733164
+        + 18.346203400819753 / np.power(-1.0 + n, 2)
+        - 728.2679478666736 / (-1.0 + n)
+        - 568.8888888888889 / np.power(n, 7)
+        + 1725.6296296296296 / np.power(n, 6)
+        - 2196.543209876543 / np.power(n, 5)
+        + 440.0487580115612 * S1
+        - 382.0574816096663 * lm11(n, S1)
+    )
 
 
 @nb.njit(cache=True)
-def gamma_gg_nf0(_n, _sx):
-    return 0
+def gamma_gg_nf0(n, sx):
+    r"""Implements the part proportional to :math:`nf^0` of :math:`\gamma_{gg}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{gg}^{(3)}|_{nf^0}`
+
+    """
+    S1 = sx[0][0]
+    return (
+        -61782.868048046716
+        - 49851.703887834694 / np.power(-1.0 + n, 4)
+        + 213823.9810748423 / np.power(-1.0 + n, 3)
+        - 241732.26024803068 / np.power(-1.0 + n, 2)
+        + 88694.20340182834 / (-1.0 + n)
+        - 103680.0 / np.power(n, 7)
+        - 17280.0 / np.power(n, 6)
+        - 627978.8224813186 / np.power(n, 5)
+        + 40880.33011934297 * S1
+        - 13643.320974357539 * lm11(n, S1)
+    )
 
 
 @nb.njit(cache=True)
 def gamma_gg(n, nf, sx):
-    """Computes the |N3LO| gluon-gluon singlet anomalous dimension.
+    r"""Computes the |N3LO| gluon-gluon singlet anomalous dimension.
 
     Parameters
     ----------
@@ -168,14 +247,14 @@ def gamma_gg(n, nf, sx):
     -------
     complex
         |N3LO| gluon-gluon singlet anomalous dimension
-        :math:`\\gamma_{gg}^{(3)}(N)`
+        :math:`\gamma_{gg}^{(3)}(N)`
 
     See Also
     --------
-    gamma_gg_nf0: :math:`\\gamma_{gg}^{(3)}|_{nf^0}`
-    gamma_gg_nf1: :math:`\\gamma_{gg}^{(3)}|_{nf^1}`
-    gamma_gg_nf2: :math:`\\gamma_{gg}^{(3)}|_{nf^2}`
-    gamma_gg_nf3: :math:`\\gamma_{gg}^{(3)}|_{nf^3}`
+    gamma_gg_nf0: :math:`\gamma_{gg}^{(3)}|_{nf^0}`
+    gamma_gg_nf1: :math:`\gamma_{gg}^{(3)}|_{nf^1}`
+    gamma_gg_nf2: :math:`\gamma_{gg}^{(3)}|_{nf^2}`
+    gamma_gg_nf3: :math:`\gamma_{gg}^{(3)}|_{nf^3}`
 
     """
     return (
diff --git a/src/eko/anomalous_dimensions/as4/ggq.py b/src/eko/anomalous_dimensions/as4/ggq.py
index 696f93647..fddbe3fd5 100644
--- a/src/eko/anomalous_dimensions/as4/ggq.py
+++ b/src/eko/anomalous_dimensions/as4/ggq.py
@@ -1,15 +1,16 @@
 # -*- coding: utf-8 -*-
-"""
-This module contains the anomalous dimension :math:`\\gamma_{gq}^{(3)}`
+r"""This module contains the anomalous dimension :math:`\gamma_{gq}^{(3)}`
 """
 import numba as nb
 import numpy as np
 
+from ...harmonics.log_functions import lm13, lm13m1, lm14, lm15
+
 
 @nb.njit(cache=True)
 def gamma_gq_nf3(n, sx):
-    """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`.
+    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`.
 
     Parameters
     ----------
@@ -21,7 +22,7 @@ def gamma_gq_nf3(n, sx):
     Returns
     -------
     complex
-        |N3LO| non-singlet anomalous dimension :math:`\\gamma_{gq}^{(3)}|_{nf^3}`
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{gq}^{(3)}|_{nf^3}`
 
     """
     S1 = sx[0][0]
@@ -52,18 +53,110 @@ def gamma_gq_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_gq_nf2(n, sx):
-    return 0
+def gamma_gq_nf0(n, sx):
+    r"""Implements the part proportional to :math:`nf^0` of :math:`\gamma_{gq}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |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]
+    return (
+        -22156.31283903764 / np.power(-1.0 + n, 4)
+        + 63019.91215580799 / np.power(-1.0 + n, 3)
+        - 52669.92830510009 / np.power(-1.0 + n, 2)
+        - 37609.87654320987 / np.power(n, 7)
+        - 35065.67901234568 / np.power(n, 6)
+        - 175454.58483973087 / np.power(n, 5)
+        - 1600.0895275985124 * lm13(n, S1, S2, S3)
+        + 21359.28988011691 * lm13m1(n, S1, S2, S3)
+        - 375.3983146907502 * lm14(n, S1, S2, S3, S4)
+        - 13.443072702331962 * lm15(n, S1, S2, S3, S4, S5)
+    )
 
 
 @nb.njit(cache=True)
 def gamma_gq_nf1(n, sx):
-    return 0
+    r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{gq}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |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]
+    return (
+        -4989.9438192798825 / np.power(-1.0 + n, 3)
+        + 9496.873384515262 / np.power(-1.0 + n, 2)
+        + 5309.62962962963 / np.power(n, 7)
+        + 221.23456790123456 / np.power(n, 6)
+        + 9092.91243376357 / np.power(n, 5)
+        + 215.7864293930184 * lm13(n, S1, S2, S3)
+        - 4337.106332800272 * lm13m1(n, S1, S2, S3)
+        + 34.49474165523548 * lm14(n, S1, S2, S3, S4)
+        + 0.5486968449931413 * lm15(n, S1, S2, S3, S4, S5)
+    )
+
+
+@nb.njit(cache=True)
+def gamma_gq_nf2(n, sx):
+    r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{gq}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |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]
+    return (
+        -215.9801828033175 / np.power(-1.0 + n, 2)
+        - 18.066114610010438 / (-1.0 + n)
+        + 778.5349794238683 / np.power(n, 5)
+        - 4.756294267863632 * lm13(n, S1, S2, S3)
+        - 44.54796646244799 * lm13m1(n, S1, S2, S3)
+        - 0.877914951989026 * lm14(n, S1, S2, S3, S4)
+    )
 
 
 @nb.njit(cache=True)
 def gamma_gq(n, nf, sx):
-    """Computes the |N3LO| gluon-quark singlet anomalous dimension.
+    r"""Computes the |N3LO| gluon-quark singlet anomalous dimension.
 
     Parameters
     ----------
@@ -78,17 +171,19 @@ def gamma_gq(n, nf, sx):
     -------
     complex
         |N3LO| gluon-quark singlet anomalous dimension
-        :math:`\\gamma_{gq}^{(3)}(N)`
+        :math:`\gamma_{gq}^{(3)}(N)`
 
     See Also
     --------
-    gamma_gq_nf1: :math:`\\gamma_{gq}^{(3)}|_{nf^1}`
-    gamma_gq_nf2: :math:`\\gamma_{gq}^{(3)}|_{nf^2}`
-    gamma_gq_nf3: :math:`\\gamma_{gq}^{(3)}|_{nf^3}`
+    gamma_gq_nf0: :math:`\gamma_{gq}^{(3)}|_{nf^0}`
+    gamma_gq_nf1: :math:`\gamma_{gq}^{(3)}|_{nf^1}`
+    gamma_gq_nf2: :math:`\gamma_{gq}^{(3)}|_{nf^2}`
+    gamma_gq_nf3: :math:`\gamma_{gq}^{(3)}|_{nf^3}`
 
     """
     return (
-        +nf * gamma_gq_nf1(n, sx)
+        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)
     )
diff --git a/src/eko/anomalous_dimensions/as4/gps.py b/src/eko/anomalous_dimensions/as4/gps.py
index 20878f063..9224ef5b9 100644
--- a/src/eko/anomalous_dimensions/as4/gps.py
+++ b/src/eko/anomalous_dimensions/as4/gps.py
@@ -1,6 +1,5 @@
 # -*- coding: utf-8 -*-
-"""
-This module contains the anomalous dimension :math:`\\gamma_{ps}^{(3)}`
+r"""This module contains the anomalous dimension :math:`\gamma_{ps}^{(3)}`
 """
 import numba as nb
 import numpy as np
@@ -8,8 +7,8 @@
 
 @nb.njit(cache=True)
 def gamma_ps_nf3(n, sx):
-    """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`.
+    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`.
 
     Parameters
     ----------
@@ -21,7 +20,7 @@ def gamma_ps_nf3(n, sx):
     Returns
     -------
     complex
-        |N3LO| non-singlet anomalous dimension :math:`\\gamma_{ps}^{(3)}|_{nf^3}`
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{ps}^{(3)}|_{nf^3}`
 
     """
     S1 = sx[0][0]
@@ -79,18 +78,67 @@ def gamma_ps_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_ps_nf2(_n, _sx):
-    return 0
+def gamma_ps_nf1(n, sx):
+    r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{ps}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{ps}^{(3)}|_{nf^1}`
+
+    """
+    S1 = sx[0][0]
+    return (
+        -10185.956311680911 * (1 / (-1.0 + n) - 1.0 / n)
+        - 3498.454512979491 / np.power(-1.0 + n, 3)
+        + 6990.974328540751 / np.power(-1.0 + n, 2)
+        + 5404.444444444444 / np.power(n, 7)
+        + 3425.9753086419755 / np.power(n, 6)
+        + 20515.223982421852 / np.power(n, 5)
+        + (418.1883933687789 * np.power(S1, 2)) / np.power(n, 2)
+        - (74.92761620078642 * np.power(S1, 3)) / np.power(n, 2)
+    )
 
 
 @nb.njit(cache=True)
-def gamma_ps_nf1(_n, _sx):
-    return 0
+def gamma_ps_nf2(n, sx):
+    r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{ps}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{ps}^{(3)}|_{nf^2}`
+
+    """
+    S1 = sx[0][0]
+    return (
+        -547.9815124741111 * (1 / (-1.0 + n) - 1.0 / n)
+        + 383.72024118761027 / np.power(-1.0 + n, 2)
+        - 568.8888888888889 / np.power(n, 7)
+        + 455.1111111111111 / np.power(n, 6)
+        - 1856.79012345679 / np.power(n, 5)
+        + (64.65057318186398 * np.power(S1, 2)) / np.power(n, 2)
+        - (7.067764355654242 * np.power(S1, 3)) / np.power(n, 2)
+    )
 
 
 @nb.njit(cache=True)
 def gamma_ps(n, nf, sx):
-    """Computes the |N3LO| pure singlet quark-quark anomalous dimension.
+    r"""Computes the |N3LO| pure singlet quark-quark anomalous dimension.
 
     Parameters
     ----------
@@ -105,13 +153,13 @@ def gamma_ps(n, nf, sx):
     -------
     complex
         |N3LO| pure singlet quark-quark anomalous dimension
-        :math:`\\gamma_{ps}^{(3)}(N)`
+        :math:`\gamma_{ps}^{(3)}(N)`
 
     See Also
     --------
-    gamma_ps_nf1: :math:`\\gamma_{ps}^{(3)}|_{nf^1}`
-    gamma_ps_nf2: :math:`\\gamma_{ps}^{(3)}|_{nf^2}`
-    gamma_ps_nf3: :math:`\\gamma_{ps}^{(3)}|_{nf^3}`
+    gamma_ps_nf1: :math:`\gamma_{ps}^{(3)}|_{nf^1}`
+    gamma_ps_nf2: :math:`\gamma_{ps}^{(3)}|_{nf^2}`
+    gamma_ps_nf3: :math:`\gamma_{ps}^{(3)}|_{nf^3}`
 
     """
     return (
diff --git a/src/eko/anomalous_dimensions/as4/gqg.py b/src/eko/anomalous_dimensions/as4/gqg.py
index 76047a249..35bc92852 100644
--- a/src/eko/anomalous_dimensions/as4/gqg.py
+++ b/src/eko/anomalous_dimensions/as4/gqg.py
@@ -1,15 +1,16 @@
 # -*- coding: utf-8 -*-
-"""
-This module contains the anomalous dimension :math:`\\gamma_{qg}^{(3)}`
+r"""This module contains the anomalous dimension :math:`\gamma_{qg}^{(3)}`
 """
 import numba as nb
 import numpy as np
 
+from ...harmonics.log_functions import lm13, lm13m1, lm14, lm15
+
 
 @nb.njit(cache=True)
 def gamma_qg_nf3(n, sx):
-    """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`
+    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`.
 
     Parameters
     ----------
@@ -21,7 +22,7 @@ def gamma_qg_nf3(n, sx):
     Returns
     -------
     complex
-        |N3LO| non-singlet anomalous dimension :math:`\\gamma_{qg}^{(3)}|_{nf^3}`
+        |N3LO| non-singlet anomalous dimension :math:`\gamma_{qg}^{(3)}|_{nf^3}`
 
     """
     S1 = sx[0][0]
@@ -329,23 +330,79 @@ def gamma_qg_nf3(n, sx):
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf2(_n, _sx):
-    return 0
+def gamma_qg_nf1(n, sx):
+    r"""Implements the part proportional to :math:`nf^1` of :math:`\gamma_{qg}^{(3)}`.
 
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
 
-@nb.njit(cache=True)
-def gamma_qg_nf1(_n, _sx):
-    return 0
+    Returns
+    -------
+    complex
+        |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]
+    return (
+        -7871.5226542038545 / np.power(-1.0 + n, 3)
+        + 13490.08729769531 / np.power(-1.0 + n, 2)
+        - 8680.177926532133 / (-1.0 + n)
+        + 14103.703703703704 / np.power(n, 7)
+        + 2588.8395061728397 / np.power(n, 6)
+        + 68802.34242841466 / np.power(n, 5)
+        - 92.30369362429062 * lm13(n, S1, S2, S3)
+        - 10782.474427912908 * lm13m1(n, S1, S2, S3)
+        - 35.68779444531073 * lm14(n, S1, S2, S3, S4)
+        - 1.8518518518518519 * lm15(n, S1, S2, S3, S4, S5)
+    )
 
 
 @nb.njit(cache=True)
-def gamma_qg_nf0(_n, _sx):
-    return 0
+def gamma_qg_nf2(n, sx):
+    r"""Implements the part proportional to :math:`nf^2` of :math:`\gamma_{qg}^{(3)}`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    sx : list
+        harmonic sums cache
+
+    Returns
+    -------
+    complex
+        |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]
+    return (
+        237.02932710506118 / np.power(-1.0 + n, 2)
+        + 447.1862550338544 / (-1.0 + n)
+        - 1991.111111111111 / np.power(n, 7)
+        + 2069.3333333333335 / np.power(n, 6)
+        - 7229.376633440217 / np.power(n, 5)
+        + 45.9259903713295 * lm13(n, S1, S2, S3)
+        + 226.61095385854318 * lm13m1(n, S1, S2, S3)
+        + 3.511659807956104 * lm14(n, S1, S2, S3, S4)
+        + 0.411522633744856 * lm15(n, S1, S2, S3, S4, S5)
+    )
 
 
 @nb.njit(cache=True)
 def gamma_qg(n, nf, sx):
-    """Computes the |N3LO| quark-gluon singlet anomalous dimension.
+    r"""Computes the |N3LO| quark-gluon singlet anomalous dimension.
 
     Parameters
     ----------
@@ -360,19 +417,17 @@ def gamma_qg(n, nf, sx):
     -------
     complex
         |N3LO| quark-gluon singlet anomalous dimension
-        :math:`\\gamma_{qg}^{(3)}(N)`
+        :math:`\gamma_{qg}^{(3)}(N)`
 
     See Also
     --------
-    gamma_qg_nf0: :math:`\\gamma_{qg}^{(3)}|_{nf^0}`
-    gamma_qg_nf1: :math:`\\gamma_{qg}^{(3)}|_{nf^1}`
-    gamma_qg_nf2: :math:`\\gamma_{qg}^{(3)}|_{nf^2}`
-    gamma_qg_nf3: :math:`\\gamma_{qg}^{(3)}|_{nf^3}`
+    gamma_qg_nf1: :math:`\gamma_{qg}^{(3)}|_{nf^1}`
+    gamma_qg_nf2: :math:`\gamma_{qg}^{(3)}|_{nf^2}`
+    gamma_qg_nf3: :math:`\gamma_{qg}^{(3)}|_{nf^3}`
 
     """
     return (
-        gamma_qg_nf0(n, sx)
-        + nf * gamma_qg_nf1(n, sx)
+        +nf * gamma_qg_nf1(n, sx)
         + nf**2 * gamma_qg_nf2(n, sx)
         + nf**3 * gamma_qg_nf3(n, sx)
     )
diff --git a/src/eko/harmonics/log_functions.py b/src/eko/harmonics/log_functions.py
index f62dde3ea..a1b1bbe89 100644
--- a/src/eko/harmonics/log_functions.py
+++ b/src/eko/harmonics/log_functions.py
@@ -1,7 +1,10 @@
 # -*- coding: utf-8 -*-
-"""Implementation of some Mellin transformation of functions as:
+r"""Implementation of Mellin transformation of logarithms.
 
-    :math:`(1-x)\\ln^k(1-x), \\quad k = 1,2,3`
+We provide transforms of:
+
+- :math:`(1-x)\ln^k(1-x), \quad k = 1,2,3`
+- :math:`\ln^k(1-x), \quad k = 1,3,4,5`
 
 """
 import numba as nb
@@ -11,39 +14,40 @@
 
 @nb.njit(cache=True)
 def lm11m1(n, S1):
-    """Mellin transform of :math:`(1-x)\\ln(1-x)`
+    r"""Mellin transform of :math:`(1-x)\ln(1-x)`.
 
     Parameters
     ----------
     n : complex
         Mellin moment
-    S1:  complex
+    S1 : complex
         Harmonic sum :math:`S_{1}(N)`
 
     Returns
     -------
     complex
-
+        :math:`\mathcal{M}[(1-x)\ln(1-x)](N)`
     """
     return 1 / (1 + n) ** 2 - S1 / (1 + n) ** 2 - S1 / (n * (1 + n) ** 2)
 
 
 @nb.njit(cache=True)
 def lm12m1(n, S1, S2):
-    """Mellin transform of :math:`(1-x)\\ln^2(1-x)`
+    r"""Mellin transform of :math:`(1-x)\ln^2(1-x)`.
 
     Parameters
     ----------
     n : complex
         Mellin moment
-    S1:  complex
+    S1 : complex
         Harmonic sum :math:`S_{1}(N)`
-    S2:  complex
+    S2 : complex
         Harmonic sum :math:`S_{2}(N)`
 
     Returns
     -------
     complex
+        :math:`\mathcal{M}[(1-x)\ln^2(1-x)](N)`
 
     """
     return (
@@ -58,22 +62,23 @@ def lm12m1(n, S1, S2):
 
 @nb.njit(cache=True)
 def lm13m1(n, S1, S2, S3):
-    """Mellin transform of :math:`(1-x)\\ln^3(1-x)`
+    r"""Mellin transform of :math:`(1-x)\ln^3(1-x)`.
 
     Parameters
     ----------
     n : complex
         Mellin moment
-    S1:  complex
+    S1 : complex
         Harmonic sum :math:`S_{1}(N)`
-    S2:  complex
+    S2 : complex
         Harmonic sum :math:`S_{2}(N)`
-    S3:  complex
+    S3 : complex
         Harmonic sum :math:`S_{3}(N)`
 
     Returns
     -------
     complex
+        :math:`\mathcal{M}[(1-x)\ln^3(1-x)](N)`
 
     """
     return (
@@ -88,3 +93,110 @@ def lm13m1(n, S1, S2, S3):
         - (2 * (6 * (2 * S3 - 2 * zeta3) + zeta3)) / n
         + (2 * (6 * (2 * S3 - 2 * zeta3) + zeta3)) / (1 + n)
     )
+
+
+@nb.njit(cache=True)
+def lm11(n, S1):
+    r"""Mellin transform of :math:`\ln(1-x)`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    S1 : complex
+        Harmonic sum :math:`S_{1}(N)`
+
+    Returns
+    -------
+    complex
+        :math:`\mathcal{M}[\ln(1-x)](N)`
+
+    """
+    return -S1 / n
+
+
+@nb.njit(cache=True)
+def lm13(n, S1, S2, S3):
+    r"""Mellin transform of :math:`\ln^3(1-x)`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    S1 : complex
+        Harmonic sum :math:`S_{1}(N)`
+    S2 : complex
+        Harmonic sum :math:`S_{2}(N)`
+    S3 : complex
+        Harmonic sum :math:`S_{3}(N)`
+
+    Returns
+    -------
+    complex
+        :math:`\mathcal{M}[\ln^3(1-x)](N)`
+
+    """
+    return -((S1**3 + 3 * S1 * S2 + 2 * S3) / n)
+
+
+@nb.njit(cache=True)
+def lm14(n, S1, S2, S3, S4):
+    r"""Mellin transform of :math:`\ln^4(1-x)`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    S1 : complex
+        Harmonic sum :math:`S_{1}(N)`
+    S2 : complex
+        Harmonic sum :math:`S_{2}(N)`
+    S3 : complex
+        Harmonic sum :math:`S_{3}(N)`
+    S4 : complex
+        Harmonic sum :math:`S_{4}(N)`
+
+    Returns
+    -------
+    complex
+        :math:`\mathcal{M}[\ln^4(1-x)](N)`
+
+    """
+    return (S1**4 + 6 * S1**2 * S2 + 3 * S2**2 + 8 * S1 * S3 + 6 * S4) / n
+
+
+@nb.njit(cache=True)
+def lm15(n, S1, S2, S3, S4, S5):
+    r"""Mellin transform of :math:`\ln^5(1-x)`.
+
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    S1 : complex
+        Harmonic sum :math:`S_{1}(N)`
+    S2 : complex
+        Harmonic sum :math:`S_{2}(N)`
+    S3 : complex
+        Harmonic sum :math:`S_{3}(N)`
+    S4 : complex
+        Harmonic sum :math:`S_{4}(N)`
+    S5 : complex
+        Harmonic sum :math:`S_{5}(N)`
+
+    Returns
+    -------
+    complex
+        :math:`\mathcal{M}[\ln^5(1-x)](N)`
+
+    """
+    return (
+        -(
+            S1**5
+            + 10 * S1**3 * S2
+            + 20 * S1**2 * S3
+            + 15 * S1 * (S2**2 + 2 * S4)
+            + 4 * (5 * S2 * S3 + 6 * S5)
+        )
+        / n
+    )
diff --git a/src/eko/interpolation.py b/src/eko/interpolation.py
index 63e679b98..118d8ffa0 100644
--- a/src/eko/interpolation.py
+++ b/src/eko/interpolation.py
@@ -681,7 +681,7 @@ def make_lambert_grid(
     .. math::
         y(x) = 5(1-x)-\log(x)
 
-    This method is implemted in `PineAPPL`, :cite:`Carrazza_2020` eq 2.11 and relative
+    This method is implemented in `PineAPPL`, :cite:`Carrazza_2020` eq 2.11 and relative
     paragraph.
 
     Parameters
diff --git a/src/eko/mellin.py b/src/eko/mellin.py
index e0b9c637a..7f68d9281 100644
--- a/src/eko/mellin.py
+++ b/src/eko/mellin.py
@@ -48,7 +48,7 @@ def Talbot_path(t, r, o):
     """
     theta = np.pi * (2.0 * t - 1.0)
     re = 0.0
-    if t == 0.5:  # treat singular point seperately
+    if t == 0.5:  # treat singular point separately
         re = 1.0
     else:
         re = theta / np.tan(theta)
@@ -81,7 +81,7 @@ def Talbot_jac(t, r, o):  # pylint: disable=unused-argument
     """
     theta = np.pi * (2.0 * t - 1.0)
     re = 0.0
-    if t == 0.5:  # treat singular point seperately
+    if t == 0.5:  # treat singular point separately
         re = 0.0
     else:
         re = 1.0 / np.tan(theta)
diff --git a/tests/eko/test_ad_as4.py b/tests/eko/test_ad_as4.py
index fbddfcbdb..56e94d027 100644
--- a/tests/eko/test_ad_as4.py
+++ b/tests/eko/test_ad_as4.py
@@ -61,13 +61,64 @@ def test_momentum_conservation():
         ggg.gamma_gg_nf3(N, sx_cache) + gqg.gamma_qg_nf3(N, sx_cache), 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),
+        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
+    )
+
+    # 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),
+        0,
+        atol=3e-13,
+    )
+    np.testing.assert_allclose(
+        ggg.gamma_gg_nf2(N, sx_cache) + gqg.gamma_qg_nf2(N, sx_cache),
+        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),
+        0,
+    )
+    np.testing.assert_allclose(
+        ggg.gamma_gg_nf1(N, sx_cache) + gqg.gamma_qg_nf1(N, sx_cache),
+        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),
+        0,
+        atol=3e-11,
+    )
+    np.testing.assert_allclose(
+        ggg.gamma_gg_nf0(N, sx_cache),
+        0,
+        atol=4e-11,
+    )
+
     # total
     g_singlet = gamma_singlet(N, NF, sx_cache)
-    # TODO: can't test for the time being since ns,+ is complete.
-    # np.testing.assert_allclose(
-    #     g_singlet[0, 0] + g_singlet[1, 0],
-    #     0,
-    # )
+    np.testing.assert_allclose(
+        g_singlet[0, 0] + g_singlet[1, 0],
+        0,
+        atol=3e-11,
+    )
     np.testing.assert_allclose(
         g_singlet[0, 1] + g_singlet[1, 1],
         0,
diff --git a/tests/eko/test_harmonics_log_functions.py b/tests/eko/test_harmonics_log_functions.py
new file mode 100644
index 000000000..0f0090d72
--- /dev/null
+++ b/tests/eko/test_harmonics_log_functions.py
@@ -0,0 +1,46 @@
+# -*- coding: utf-8 -*-
+import numpy as np
+from scipy.integrate import quad
+
+from eko import harmonics as h
+
+
+def test_lm1pm1():
+    # test mellin transformation with some random N values
+    def mellin_lm1pm1(x, k):
+        return x ** (N - 1) * (1 - x) * np.log(1 - x) ** k
+
+    Ns = [12.345, 81.113, 27.787]
+    for N in Ns:
+        sx = h.sx(N, 3)
+
+        ref_values = {
+            1: h.log_functions.lm11m1(N, sx[0]),
+            2: h.log_functions.lm12m1(N, sx[0], sx[1]),
+            3: h.log_functions.lm13m1(N, sx[0], sx[1], sx[2]),
+        }
+
+        for k in [1, 2, 3]:
+            test_value = quad(mellin_lm1pm1, 0, 1, args=(k))[0]
+            np.testing.assert_allclose(test_value, ref_values[k], atol=5e-4)
+
+
+def test_lm1p():
+    # test mellin transformation with some random N values
+    def mellin_lm1p(x, k):
+        return x ** (N - 1) * np.log(1 - x) ** k
+
+    Ns = [65.780, 56.185, 94.872]
+    for N in Ns:
+        sx = h.sx(N, 5)
+
+        ref_values = {
+            1: h.log_functions.lm11(N, sx[0]),
+            3: h.log_functions.lm13(N, sx[0], sx[1], sx[2]),
+            4: h.log_functions.lm14(N, sx[0], sx[1], sx[2], sx[3]),
+            5: h.log_functions.lm15(N, sx[0], sx[1], sx[2], sx[3], sx[4]),
+        }
+
+        for k in [1, 3, 4, 5]:
+            test_value = quad(mellin_lm1p, 0, 1, args=(k))[0]
+            np.testing.assert_allclose(test_value, ref_values[k])