*
* ..File: xpij1e.f  
*
*
* ..The 1- and 2-loop MS(bar) singlet splitting functions  P_ij^(2) 
*    for the evolution of unpolarized 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.
*
*
* =====================================================================
*
*
* ..The 2-loop pure-singlet splitting functions P_ps^(1)
*
       FUNCTION X1PSA (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 = 4 )
       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)
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       DX = 1.D0/X
*
* ...The harmonic polylogs up to weight 4 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
*
       X1PSA =
     &  + nf*cf * (  - 8.D0 + 24.D0*x - 224.D0/9.D0*x**2 + 80.D0/9.D0*
     &    dx + 4.D0*Hr1(0) + 20.D0*Hr1(0)*x + 32.D0/3.D0*Hr1(0)*x**2 -
     &    8.D0*Hr2(0,0) - 8.D0*Hr2(0,0)*x )
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop gluon->quark splitting functions P_qg^(1)
*
       FUNCTION X1QGA (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 = 4 )
       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 up to weight 
*
       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
*
       gqg1 =
     &  + nf*ca * (  - 8.D0 + 100.D0*x - 872.D0/9.D0*x**2 + 80.D0/9.D0*
     &    dx - 16.D0*z2*x + 4.D0*Hr1(0) + 32.D0*Hr1(0)*x + 176.D0/3.D0*
     &    Hr1(0)*x**2 + 16.D0*Hr1(1)*x - 16.D0*Hr1(1)*x**2 - 8.D0*Hr2(
     &    -1,0) - 16.D0*Hr2(-1,0)*x - 16.D0*Hr2(-1,0)*x**2 - 8.D0*Hr2(0
     &    ,0) - 16.D0*Hr2(0,0)*x - 8.D0*Hr2(1,1) + 16.D0*Hr2(1,1)*x -
     &    16.D0*Hr2(1,1)*x**2 )
       gqg1 = gqg1 + nf*cf * ( 28.D0 - 58.D0*x + 40.D0*x**2 - 8.D0*z2 +
     &    16.D0*z2*x - 16.D0*z2*x**2 + 6.D0*Hr1(0) - 8.D0*Hr1(0)*x + 16.
     &    D0*Hr1(0)*x**2 - 16.D0*Hr1(1)*x + 16.D0*Hr1(1)*x**2 + 4.D0*
     &    Hr2(0,0) - 8.D0*Hr2(0,0)*x + 16.D0*Hr2(0,0)*x**2 + 8.D0*Hr2(0
     &    ,1) - 16.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*x**2 + 8.D0*Hr2(1,0)
     &     - 16.D0*Hr2(1,0)*x + 16.D0*Hr2(1,0)*x**2 + 8.D0*Hr2(1,1) -
     &    16.D0*Hr2(1,1)*x + 16.D0*Hr2(1,1)*x**2 )
*
       X1QGA = gqg1
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 2-loop quark->gluon splitting functions P_gq^(1)
*
       FUNCTION X1GQA (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 = 4 ) 
       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 up to weight 
*
       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
*
      ggq1 =
     &  + cf*ca * ( 76.D0/9.D0 + 148.D0/9.D0*x + 176.D0/9.D0*x**2 + 4.D0
     &    *dx + 16.D0*z2 - 48.D0*Hr1(0) - 20.D0*Hr1(0)*x - 32.D0/3.D0*
     &    Hr1(0)*x**2 + 88.D0/3.D0*Hr1(1) - 68.D0/3.D0*Hr1(1)*x - 88.D0/
     &    3.D0*Hr1(1)*dx + 16.D0*Hr2(-1,0) + 8.D0*Hr2(-1,0)*x + 16.D0*
     &    Hr2(-1,0)*dx + 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(0,1)*dx - 16.D0*Hr2(1,0) +
     &    8.D0*Hr2(1,0)*x + 16.D0*Hr2(1,0)*dx - 16.D0*Hr2(1,1) + 8.D0*
     &    Hr2(1,1)*x + 16.D0*Hr2(1,1)*dx )
      ggq1 = ggq1 + cf**2 * (  - 10.D0 - 14.D0*x + 8.D0*Hr1(0) + 14.D0*
     &    Hr1(0)*x - 24.D0*Hr1(1) + 20.D0*Hr1(1)*x + 24.D0*Hr1(1)*dx -
     &    8.D0*Hr2(0,0) + 4.D0*Hr2(0,0)*x + 16.D0*Hr2(1,1) - 8.D0*Hr2(1
     &    ,1)*x - 16.D0*Hr2(1,1)*dx )
      ggq1 = ggq1 + nf*cf * ( 80.D0/9.D0 - 64.D0/9.D0*x - 80.D0/9.D0*dx
     &     - 16.D0/3.D0*Hr1(1) + 8.D0/3.D0*Hr1(1)*x + 16.D0/3.D0*Hr1(1)
     &    *dx )
*   
       X1GQA = ggq1       
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The regular piece of the 2-loop gg splitting function P_gg^(1) 
*
       FUNCTION X1GGA (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 = 4 ) 
       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 X1GGB and X1GGC
*
       COMMON / P1GSOFT / A2G
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
*
* ...Some abbreviations
*
       DX = 1.D0/X
       DM = 1.D0/(1.D0-X)
       DP = 1.D0/(1.D0+X)
*
* ...The harmonic polylogs up to weight
*
       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)
*
      ggg1 =
     &  + ca**2 * (  - 50.D0/9.D0 - 218.D0/9.D0*x + 268.D0/9.D0*dm + 32.
     &    D0*z2 + 16.D0*z2*x**2 - 8.D0*z2*dp - 8.D0*z2*dm - 100.D0/3.D0
     &    *Hr1(0) + 44.D0/3.D0*Hr1(0)*x - 176.D0/3.D0*Hr1(0)*x**2 + 32.D
     &    0*Hr2(-1,0) + 16.D0*Hr2(-1,0)*x + 16.D0*Hr2(-1,0)*x**2 + 16.D0
     &    *Hr2(-1,0)*dx - 16.D0*Hr2(-1,0)*dp + 32.D0*Hr2(0,0)*x - 16.D0
     &    *Hr2(0,0)*x**2 + 8.D0*Hr2(0,0)*dp + 8.D0*Hr2(0,0)*dm - 32.D0*
     &    Hr2(0,1) + 16.D0*Hr2(0,1)*x - 16.D0*Hr2(0,1)*x**2 + 16.D0*
     &    Hr2(0,1)*dx + 16.D0*Hr2(0,1)*dm - 32.D0*Hr2(1,0) + 16.D0*Hr2(
     &    1,0)*x - 16.D0*Hr2(1,0)*x**2 + 16.D0*Hr2(1,0)*dx + 16.D0*Hr2(
     &    1,0)*dm )
      ggg1 = ggg1 + nf*ca * ( 116.D0/9.D0 - 76.D0/9.D0*x + 92.D0/9.D0*
     &    x**2 - 92.D0/9.D0*dx - 40.D0/9.D0*dm - 8.D0/3.D0*Hr1(0) - 8.D0
     &    /3.D0*Hr1(0)*x )
      ggg1 = ggg1 + nf*cf * (  - 32.D0 + 16.D0*x + 40.D0/3.D0*x**2 + 8.D
     &    0/3.D0*dx - 12.D0*Hr1(0) - 20.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
*
       GGG1L = DM * A2G
*
* ...The regular piece of the coefficient function
*
       X1GGA = GGG1 - GGG1L
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the singular (soft) piece.
*
       FUNCTION X1GGB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P1GSOFT / A2G
*
       X1GGB  = A2G/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the 'local' piece.
*
       FUNCTION X1GGC (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)
*
       P1DELT = 
     ,    - 2.D0*cf*nf
     ,    - 8.D0/3.D0*ca*nf
     ,    + 32.D0/3.D0*ca**2
     ,    + 12.D0*z3*ca**2
*
       X1GGC = DLOG (1.D0-Y) * A2G + P1DELT
*
       RETURN
       END
*
* =====================================================================
*
*
* ..The 1-loop gluon->quark splitting functions P_qg^(0)
*
       FUNCTION X0QGA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       X0QGA = 2.* NF * ( 1. - 2. * X + 2. * X**2 )
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 1-loop quark->gluon splitting functions P_gq^(0)
*
       FUNCTION X0GQA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       CF = 4./3.D0
       X0GQA = 4.*CF * ( - 1. + 0.5 * X + 1./X )
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The regular piece of the 1-loop gg splitting function P_gg^(0)
*
       FUNCTION X0GGA (X,NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       CA = 3.D0
       X0GGA = 4.*CA * ( - 2. + X - X**2 + 1./X )
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the singular (soft) piece.
*
       FUNCTION X0GGB (X,NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       CA = 3.D0
       X0GGB = 4.*CA / (1.D0-X)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the 'local' piece.
*
       FUNCTION X0GGC (X,NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       CA = 3.D0
       X0GGC = 4.*CA * LOG (1.D0-X) - 2./3.D0 * NF + 11./3.D0 * CA
*
       RETURN
       END
*
* =================================================================av==
