*
* ..File: xdPij1M.f  
*
*
* ..The NLO MS(bar) [MvN version] splitting functions  DeltaP_ij^(1)
*    for the evolution of polarized partons densities, mu_r = mu_f.
*    The expansion parameter is alpha_s/(4 pi).
* 
* ..The distributions (in the mathematical sense) are given as in eq.
*    (B.26) of Floratos, Kounnas, Lacaze: Nucl. Phys. B192 (1981) 417.
*    The name-endings A, B, and C of the functions below correspond to 
*    the kernel superscripts [2], [3], and [1] in that equation.
*
* ..The code uses the package of Gehrmann and Remiddi for the harmonic
*    polylogarithms published in hep-ph/0107173 = CPC 141 (2001) 296.
*
* ..These are well-known functions, first obtained 1995 by Mertig and 
*    van Neerven and by Vogelsang, which were re-calulated during the
*    determinations of their NNLO counterparts presented in
*
*    S. Moch, J. Vermaseren and A. Vogt, arXiv:1409.5131
*
* =====================================================================
*
*
* ..2-loop quark-quark splitting function  DeltaP_qq^(1), regular part.
*    DeltaP_ps^(1) is obtained by subtracting the unpolarized P_NS^(1)-
*
       FUNCTION DP1QQMA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 )
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2),
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2),
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2),
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2)
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0 )
*
* ..The soft coefficient for use in DP1QQB and DP1QQC
*
       COMMON / P1SOFT / A2
*
* ...Colour factors and abbreviations
*
       CF  = 4./3.D0
       CA  = 3.D0
*
       DX = 1.D0/X
       DM = 1.D0/(1.D0-X)
       DP = 1.D0/(1.D0+X)
*
* ...The harmonic polylogs by Gehrmann and Remiddi
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2)
*
* ...The splitting function in terms of the harmonic polylogs
*    (without the delta(1-x) part, but with the soft contribution)
*
      dPqq1M =
     &  + cf*ca * ( 178.D0/9.D0 - 446.D0/9.D0*x + 268.D0/9.D0*dm + 8.D0
     &    *z2 - 8.D0*z2*dp - 8.D0*z2*dm + 2.D0/3.D0*Hr1(0) + 2.D0/3.D0*
     &    Hr1(0)*x + 44.D0/3.D0*Hr1(0)*dm + 8.D0*Hr2(-1,0) - 8.D0*Hr2(
     &    -1,0)*x - 16.D0*Hr2(-1,0)*dp - 8.D0*Hr2(0,0) + 8.D0*Hr2(0,0)*
     &    dp + 8.D0*Hr2(0,0)*dm )
      dPqq1M = dPqq1M + cf**2 * (  - 36.D0 + 36.D0*x - 8.D0*z2 + 8.D0
     &    *z2*x + 16.D0*z2*dp - 8.D0*Hr1(0) - 16.D0*Hr1(0)*x - 12.D0*
     &    Hr1(0)*dm - 16.D0*Hr2(-1,0) + 16.D0*Hr2(-1,0)*x + 32.D0*Hr2(
     &    -1,0)*dp + 4.D0*Hr2(0,0) - 12.D0*Hr2(0,0)*x - 16.D0*Hr2(0,0)*
     &    dp - 8.D0*Hr2(0,1) - 8.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*dm - 8.D
     &    0*Hr2(1,0) - 8.D0*Hr2(1,0)*x + 16.D0*Hr2(1,0)*dm )
      dPqq1M = dPqq1M + nf*cf * ( 32.D0/9.D0 + 8.D0/9.D0*x - 40.D0/9.D
     &    0*dm - 8.D0/3.D0*Hr1(0) + 40.D0/3.D0*Hr1(0)*x - 8.D0/3.D0*
     &    Hr1(0)*dm - 8.D0*Hr2(0,0) - 8.D0*Hr2(0,0)*x )
*
* ...The soft (`+'-distribution) part of the splitting function
*
       A2 = - 40.D0/9.D0*cf*nf + 268.D0/9.D0*ca*cf - 8.D0*z2*ca*cf
       DPQQ1L = DM * A2
*
* ...The regular piece of the splitting function
*
       DP1QQMA = dPqq1M - DPQQ1L
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The singular (soft) piece of DeltaP_qq^(1)
*
       FUNCTION DP1QQB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P1SOFT / A2
*
       DP1QQB  = A2/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The 'local' piece of DeltaP_qq^(1)
*
       FUNCTION DP1QQC (Y, NF)
*
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF, NF2
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0 )
*
       COMMON / P1SOFT / A2
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
*
* ...The coefficient of delta(1-x)
*
       P1DELT = 
     &     - 1.D0/3.D0*cf*nf
     &     + 3.D0/2.D0*cf**2
     &     + 17.D0/6.D0*ca*cf
     &     + 24.D0*z3*cf**2
     &     - 12.D0*z3*ca*cf
     &     - 8.D0/3.D0*z2*cf*nf
     &     - 12.D0*z2*cf**2
     &     + 44.D0/3.D0*z2*ca*cf
*
       DP1QQC = LOG (1.D0-Y) * A2 + P1DELT
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop gluon->quark splitting function DeltaP_qg^(1)
*
       FUNCTION DP1QGMA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 )
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2),
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2),
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2),
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2)
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0 )
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       DX = 1.D0/X
*
* ...The harmonic polylogs by Gehrmann and Remiddi
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2)
*
* ...The splitting function in terms of the harmonic polylogs
*
      dPqg1M =
     &  + nf*ca * ( 48.D0 - 44.D0*x - 8.D0*z2 + 4.D0*Hr1(0) + 32.D0*
     &    Hr1(0)*x + 16.D0*Hr1(1) - 16.D0*Hr1(1)*x - 8.D0*Hr2(-1,0) -
     &    16.D0*Hr2(-1,0)*x - 8.D0*Hr2(0,0) - 16.D0*Hr2(0,0)*x + 8.D0*
     &    Hr2(1,1) - 16.D0*Hr2(1,1)*x )
      dPqg1M = dPqg1M + nf*cf * (  - 44.D0 + 54.D0*x + 8.D0*z2 - 16.D0
     &    *z2*x - 18.D0*Hr1(0) - 16.D0*Hr1(1) + 16.D0*Hr1(1)*x - 4.D0*
     &    Hr2(0,0) + 8.D0*Hr2(0,0)*x - 8.D0*Hr2(0,1) + 16.D0*Hr2(0,1)*x
     &     - 8.D0*Hr2(1,0) + 16.D0*Hr2(1,0)*x - 8.D0*Hr2(1,1) + 16.D0*
     &    Hr2(1,1)*x )
*
       DP1QGMA = dPqg1M
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop quark->gluon splitting function DeltaP_gq^(1)
*
       FUNCTION DP1GQMA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 )
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2),
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2),
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2)
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2),
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2)
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0 )
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       DX = 1.D0/X
*
* ...The harmonic polylogs by Gehrmann and Remiddi
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2)
*
* ...The splitting function in terms of the harmonic polylogs
*
      dPgq1M =
     &  + cf*ca * ( 164.D0/9.D0 + 140.D0/9.D0*x + 8.D0*z2*x + 16.D0*
     &    Hr1(0) - 52.D0*Hr1(0)*x - 40.D0/3.D0*Hr1(1) - 4.D0/3.D0*Hr1(1
     &    )*x + 16.D0*Hr2(-1,0) + 8.D0*Hr2(-1,0)*x + 16.D0*Hr2(0,0) + 8.
     &    D0*Hr2(0,0)*x + 16.D0*Hr2(0,1) - 8.D0*Hr2(0,1)*x + 16.D0*Hr2(
     &    1,0) - 8.D0*Hr2(1,0)*x + 16.D0*Hr2(1,1) - 8.D0*Hr2(1,1)*x )
      dPgq1M = dPgq1M + cf**2 * (  - 34.D0 + 16.D0*x - 8.D0*Hr1(0) +
     &    2.D0*Hr1(0)*x + 8.D0*Hr1(1) + 4.D0*Hr1(1)*x + 8.D0*Hr2(0,0)
     &     - 4.D0*Hr2(0,0)*x - 16.D0*Hr2(1,1) + 8.D0*Hr2(1,1)*x )
      dPgq1M = dPgq1M + nf*cf * (  - 32.D0/9.D0 - 8.D0/9.D0*x + 16.D0/
     &    3.D0*Hr1(1) - 8.D0/3.D0*Hr1(1)*x )
*
       DP1GQMA = dPgq1M
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..2-loop gluon-gluon splitting function DeltaP_gg^(1), regular part
*
       FUNCTION DP1GGMA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4 
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N
       PARAMETER ( N1 = -1, N2 = 1, NW = 2 ) 
       DIMENSION HC1(N1:N2),HC2(N1:N2,N1:N2),HC3(N1:N2,N1:N2,N1:N2), 
     ,           HC4(N1:N2,N1:N2,N1:N2,N1:N2) 
       DIMENSION HR1(N1:N2),HR2(N1:N2,N1:N2),HR3(N1:N2,N1:N2,N1:N2), 
     ,           HR4(N1:N2,N1:N2,N1:N2,N1:N2) 
       DIMENSION HI1(N1:N2),HI2(N1:N2,N1:N2),HI3(N1:N2,N1:N2,N1:N2), 
     ,           HI4(N1:N2,N1:N2,N1:N2,N1:N2) 
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0 )
*
* ...The soft coefficient for use in DP1GGB and DP1GGC
*
       COMMON / P1GSOFT / A2G
*
* ...Colour factors and abbreviations
*
       CF  = 4./3.D0
       CA  = 3.D0
*
       DX = 1.D0/X
       DM = 1.D0/(1.D0-X)
       DP = 1.D0/(1.D0+X)
*
* ...The harmonic polylogs by Gehrmann and Remiddi
*
       CALL HPLOG (X, NW, HC1,HC2,HC3,HC4, HR1,HR2,HR3,HR4,
     ,            HI1,HI2,HI3,HI4, N1, N2) 
*
* ...The splitting function in terms of the harmonic polylogs
*    (without the delta(1-x) part, but with the soft contribution)
*
      dPgg1M =
     &  + ca**2 * (  - 74.D0/9.D0 - 194.D0/9.D0*x + 268.D0/9.D0*dm + 32.
     &    D0*z2*x + 8.D0*z2*dp - 8.D0*z2*dm + 116.D0/3.D0*Hr1(0) - 268.D
     &    0/3.D0*Hr1(0)*x + 16.D0*Hr2(-1,0) + 32.D0*Hr2(-1,0)*x + 16.D0
     &    *Hr2(-1,0)*dp + 32.D0*Hr2(0,0) - 8.D0*Hr2(0,0)*dp + 8.D0*Hr2(
     &    0,0)*dm + 16.D0*Hr2(0,1) - 32.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*
     &    dm + 16.D0*Hr2(1,0) - 32.D0*Hr2(1,0)*x + 16.D0*Hr2(1,0)*dm )
      dPgg1M = dPgg1M + nf*ca * (  - 112.D0/9.D0 + 152.D0/9.D0*x - 40.
     &    D0/9.D0*dm - 8.D0/3.D0*Hr1(0) - 8.D0/3.D0*Hr1(0)*x )
      dPgg1M = dPgg1M + nf*cf * (  - 20.D0 + 20.D0*x - 20.D0*Hr1(0)
     &     + 4.D0*Hr1(0)*x - 8.D0*Hr2(0,0) - 8.D0*Hr2(0,0)*x )
*
* ...The soft (`+'-distribution) part of the splitting function
*
       A2G = - 40.D0/9.D0*ca*nf + 268.D0/9.D0*ca**2 - 8.D0*z2*ca**2
       DPGG1L = DM * A2G
*
* ...The regular piece of the coefficient function
*
       DP1GGMA = dPgg1M - DPGG1L
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The singular (soft) piece of DeltaP_gg^(1)
*
       FUNCTION DP1GGB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P1GSOFT / A2G
*
       DP1GGB  = A2G/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The 'local' piece of DeltaP_gg^(1)
*
       FUNCTION DP1GGC (Y, NF)
*
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0 )
*
       COMMON / P1GSOFT / A2G
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
*
* ...The coefficient of delta(1-x)
*
       P1GDELT = 
     ,    - 2.D0*cf*nf
     ,    - 8.D0/3.D0*ca*nf
     ,    + 32.D0/3.D0*ca**2
     ,    + 12.D0*z3*ca**2
*
       DP1GGC = DLOG (1.D0-Y) * A2G + P1GDELT
*
       RETURN
       END
*
* =================================================================av==
