*
* ..File: xdPij2M.f  
*
*                  __
* ..The exact NNLO MS (`M' version, as defined in hep-ph/9803439)
*    flavour-singlet splitting functions  DeltaP^(2)  for the evolution
*    of helicity-difference (polarized) parton densities at 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.
*
* ..Reference:   S. Moch, J. Vermaseren and A. Vogt,
*                The Three-Loop Splitting Functions in QCD: 
*                The Helicity-Dependent Case 
*                arXiv:1409.5131 
*
* =====================================================================
*
*
* ..3-loop quark-quark splitting function  DeltaP_qq^(2), regular part.
*
       FUNCTION DP2QQMA (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,
     ,             Z4 = 1.0823 23233 71113 81915 D0,
     ,             Z5 = 1.0369 27755 14336 99263 D0 )
*
* ..The soft coefficient for use in DP2QQB and DP2QQC
*
       COMMON / P2SOFT / A3
*
* ...Colour factors and abbreviations
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
       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)
*
      dPqq2M =
     &  + cf*ca**2 * ( 923.D0/27.D0 - 5333.D0/27.D0*x - 121.D0*z4 - 31.D
     &    0*z4*x + 80.D0*dp*z4 + 490.D0/3.D0*dm + 192.D0*dm*z4 + 112.D0
     &    *z3 + 48.D0*z3*x + 88.D0*z3*dp + 16.D0*z3*dm + 1504.D0/9.D0*
     &    z2 + 224.D0/3.D0*z2*x - 1072.D0/9.D0*z2*dp - 1072.D0/9.D0*z2*
     &    dm + 192.D0*Hr1(-1)*z3 - 192.D0*Hr1(-1)*z3*x - 384.D0*Hr1(-1)
     &    *z3*dp - 208.D0/3.D0*Hr1(-1)*z2 - 560.D0/3.D0*Hr1(-1)*z2*x - 
     &    352.D0/3.D0*Hr1(-1)*z2*dp - 1378.D0/27.D0*Hr1(0) - 2266.D0/27.
     &    D0*Hr1(0)*x + 24.D0*Hr1(0)*dp + 4172.D0/27.D0*Hr1(0)*dm - 144.
     &    D0*Hr1(0)*z3 + 128.D0*Hr1(0)*z3*dp + 128.D0*Hr1(0)*z3*dm + 
     &    332.D0/3.D0*Hr1(0)*z2 + 188.D0*Hr1(0)*z2*x + 16.D0/3.D0*Hr1(0
     &    )*z2*dp - 248.D0/3.D0*Hr1(0)*z2*dm - 592.D0/3.D0*Hr1(1) + 592.
     &    D0/3.D0*Hr1(1)*x - 144.D0*Hr1(1)*z3 - 144.D0*Hr1(1)*z3*x + 
     &    288.D0*Hr1(1)*z3*dm - 32.D0*Hr1(1)*z2 + 32.D0*Hr1(1)*z2*x - 
     &    256.D0*Hr2(-1,-1)*z2 + 256.D0*Hr2(-1,-1)*z2*x + 512.D0*Hr2(-1
     &    ,-1)*z2*dp + 832.D0/9.D0*Hr2(-1,0) - 1312.D0/9.D0*Hr2(-1,0)*x
     &     - 2144.D0/9.D0*Hr2(-1,0)*dp + 176.D0*Hr2(-1,0)*z2 - 176.D0*
     &    Hr2(-1,0)*z2*x - 352.D0*Hr2(-1,0)*z2*dp + 136.D0*Hr2(0,-1)*z2
     &     - 136.D0*Hr2(0,-1)*z2*x - 256.D0*Hr2(0,-1)*z2*dp - 130.D0*
     &    Hr2(0,0) + 238.D0/9.D0*Hr2(0,0)*x + 1072.D0/9.D0*Hr2(0,0)*dp
     &     + 1556.D0/9.D0*Hr2(0,0)*dm - 68.D0*Hr2(0,0)*z2 + 36.D0*Hr2(0
     &    ,0)*z2*x + 96.D0*Hr2(0,0)*z2*dp - 224.D0/3.D0*Hr2(0,1) - 224.D
     &    0/3.D0*Hr2(0,1)*x + 24.D0*Hr2(0,1)*z2 - 40.D0*Hr2(0,1)*z2*x
     &     - 64.D0*Hr2(0,1)*z2*dp + 16.D0*Hr2(1,0)*z2 + 16.D0*Hr2(1,0)*
     &    z2*x - 32.D0*Hr2(1,0)*z2*dm + 64.D0*Hr3(-1,-1,0) + 64.D0*Hr3(
     &    -1,-1,0)*x + 232.D0/3.D0*Hr3(-1,0,0) + 56.D0/3.D0*Hr3(-1,0,0)
     &    *x - 176.D0/3.D0*Hr3(-1,0,0)*dp + 304.D0/3.D0*Hr3(-1,0,1) + 
     &    656.D0/3.D0*Hr3(-1,0,1)*x + 352.D0/3.D0*Hr3(-1,0,1)*dp + 328.D
     &    0/3.D0*Hr3(0,-1,0) - 376.D0/3.D0*Hr3(0,-1,0)*x - 176.D0/3.D0*
     &    Hr3(0,-1,0)*dp - 48.D0*Hr3(0,-1,0)*dm - 416.D0/3.D0*Hr3(0,0,0
     &    ) + 248.D0/3.D0*Hr3(0,0,0)*dp + 248.D0/3.D0*Hr3(0,0,0)*dm - 4.
     &    D0/3.D0*Hr3(0,0,1) - 188.D0*Hr3(0,0,1)*x - 176.D0/3.D0*Hr3(0,
     &    0,1)*dp + 88.D0/3.D0*Hr3(0,0,1)*dm + 12.D0*Hr3(1,0,0) + 76.D0
     &    *Hr3(1,0,0)*x - 88.D0*Hr3(1,0,0)*dm + 128.D0*Hr4(-1,-1,0,0)
     &     - 128.D0*Hr4(-1,-1,0,0)*x - 256.D0*Hr4(-1,-1,0,0)*dp + 256.D0
     &    *Hr4(-1,-1,0,1) - 256.D0*Hr4(-1,-1,0,1)*x - 512.D0*Hr4(-1,-1,
     &    0,1)*dp - 48.D0*Hr4(-1,0,0,0) + 48.D0*Hr4(-1,0,0,0)*x + 96.D0
     &    *Hr4(-1,0,0,0)*dp - 128.D0*Hr4(-1,0,0,1) + 128.D0*Hr4(-1,0,0,
     &    1)*x + 256.D0*Hr4(-1,0,0,1)*dp - 48.D0*Hr4(0,-1,-1,0) - 80.D0
     &    *Hr4(0,-1,-1,0)*x + 128.D0*Hr4(0,-1,-1,0)*dm - 56.D0*Hr4(0,-1
     &    ,0,0) + 88.D0*Hr4(0,-1,0,0)*x + 128.D0*Hr4(0,-1,0,0)*dp - 32.D
     &    0*Hr4(0,-1,0,0)*dm - 160.D0*Hr4(0,-1,0,1) + 96.D0*Hr4(0,-1,0,
     &    1)*x + 256.D0*Hr4(0,-1,0,1)*dp + 64.D0*Hr4(0,-1,0,1)*dm - 40.D
     &    0*Hr4(0,0,-1,0) + 8.D0*Hr4(0,0,-1,0)*x + 32.D0*Hr4(0,0,-1,0)*
     &    dp + 32.D0*Hr4(0,0,-1,0)*dm + 40.D0*Hr4(0,0,0,0) - 32.D0*Hr4(
     &    0,0,0,0)*dp - 32.D0*Hr4(0,0,0,0)*dm + 28.D0*Hr4(0,0,0,1) - 36.
     &    D0*Hr4(0,0,0,1)*x - 64.D0*Hr4(0,0,0,1)*dp + 32.D0*Hr4(0,0,0,1
     &    )*dm + 28.D0*Hr4(0,1,0,0) + 28.D0*Hr4(0,1,0,0)*x - 64.D0*Hr4(
     &    0,1,0,0)*dm - 96.D0*Hr4(1,0,-1,0) - 96.D0*Hr4(1,0,-1,0)*x + 
     &    192.D0*Hr4(1,0,-1,0)*dm + 48.D0*Hr4(1,0,0,0) + 48.D0*Hr4(1,0,
     &    0,0)*x - 96.D0*Hr4(1,0,0,0)*dm - 64.D0*Hr4(1,0,0,1) - 64.D0*
     &    Hr4(1,0,0,1)*x + 128.D0*Hr4(1,0,0,1)*dm + 64.D0*Hr4(1,1,0,0)
     &     + 64.D0*Hr4(1,1,0,0)*x - 128.D0*Hr4(1,1,0,0)*dm )
      dPqq2M = dPqq2M + cf**2*ca * ( 1516.D0/3.D0 - 1516.D0/3.D0*x + 
     &    28.D0*z4 - 12.D0*z4*x - 20.D0*dp*z4 - 276.D0*dm*z4 - 832.D0/3.
     &    D0*z3 - 880.D0/3.D0*z3*x - 248.D0*z3*dp + 80.D0/3.D0*z3*dm - 
     &    3880.D0/9.D0*z2 - 152.D0/9.D0*z2*x + 2144.D0/9.D0*z2*dp - 672.
     &    D0*Hr1(-1)*z3 + 672.D0*Hr1(-1)*z3*x + 1344.D0*Hr1(-1)*z3*dp
     &     + 704.D0/3.D0*Hr1(-1)*z2 + 1696.D0/3.D0*Hr1(-1)*z2*x + 992.D0
     &    /3.D0*Hr1(-1)*z2*dp + 3404.D0/9.D0*Hr1(0) - 3740.D0/9.D0*Hr1(
     &    0)*x - 72.D0*Hr1(0)*dp - 302.D0/3.D0*Hr1(0)*dm + 400.D0*Hr1(0
     &    )*z3 - 144.D0*Hr1(0)*z3*x - 464.D0*Hr1(0)*z3*dp - 272.D0*Hr1(
     &    0)*z3*dm - 904.D0/3.D0*Hr1(0)*z2 - 568.D0*Hr1(0)*z2*x - 104.D0
     &    /3.D0*Hr1(0)*z2*dp + 328.D0/3.D0*Hr1(0)*z2*dm + 2008.D0/3.D0*
     &    Hr1(1) - 2008.D0/3.D0*Hr1(1)*x + 384.D0*Hr1(1)*z3 + 384.D0*
     &    Hr1(1)*z3*x - 768.D0*Hr1(1)*z3*dm + 112.D0*Hr1(1)*z2 - 112.D0
     &    *Hr1(1)*z2*x + 896.D0*Hr2(-1,-1)*z2 - 896.D0*Hr2(-1,-1)*z2*x
     &     - 1792.D0*Hr2(-1,-1)*z2*dp - 1232.D0/9.D0*Hr2(-1,0) + 3056.D0
     &    /9.D0*Hr2(-1,0)*x + 4288.D0/9.D0*Hr2(-1,0)*dp - 672.D0*Hr2(-1
     &    ,0)*z2 + 672.D0*Hr2(-1,0)*z2*x + 1344.D0*Hr2(-1,0)*z2*dp - 
     &    480.D0*Hr2(0,-1)*z2 + 576.D0*Hr2(0,-1)*z2*x + 1024.D0*Hr2(0,
     &    -1)*z2*dp - 96.D0*Hr2(0,-1)*z2*dm + 2168.D0/9.D0*Hr2(0,0) - 
     &    944.D0/3.D0*Hr2(0,0)*x - 2144.D0/9.D0*Hr2(0,0)*dp - 104.D0*
     &    Hr2(0,0)*dm + 352.D0*Hr2(0,0)*z2 - 96.D0*Hr2(0,0)*z2*x - 416.D
     &    0*Hr2(0,0)*z2*dp - 128.D0*Hr2(0,0)*z2*dm + 2648.D0/9.D0*Hr2(0
     &    ,1) + 152.D0/9.D0*Hr2(0,1)*x + 2144.D0/9.D0*Hr2(0,1)*dm - 64.D
     &    0*Hr2(0,1)*z2 + 160.D0*Hr2(0,1)*z2*x + 224.D0*Hr2(0,1)*z2*dp
     &     - 64.D0*Hr2(0,1)*z2*dm + 176.D0/9.D0*Hr2(1,0) - 2320.D0/9.D0
     &    *Hr2(1,0)*x + 2144.D0/9.D0*Hr2(1,0)*dm + 32.D0*Hr2(1,0)*z2 + 
     &    32.D0*Hr2(1,0)*z2*x - 64.D0*Hr2(1,0)*z2*dm - 224.D0*Hr3(-1,-1
     &    ,0) - 224.D0*Hr3(-1,-1,0)*x - 872.D0/3.D0*Hr3(-1,0,0) - 376.D0
     &    /3.D0*Hr3(-1,0,0)*x + 496.D0/3.D0*Hr3(-1,0,0)*dp - 1040.D0/3.D
     &    0*Hr3(-1,0,1) - 2032.D0/3.D0*Hr3(-1,0,1)*x - 992.D0/3.D0*Hr3(
     &    -1,0,1)*dp - 752.D0/3.D0*Hr3(0,-1,0) + 944.D0/3.D0*Hr3(0,-1,0
     &    )*x + 496.D0/3.D0*Hr3(0,-1,0)*dp + 96.D0*Hr3(0,-1,0)*dm + 868.
     &    D0/3.D0*Hr3(0,0,0) - 36.D0*Hr3(0,0,0)*x - 712.D0/3.D0*Hr3(0,0
     &    ,0)*dp - 184.D0/3.D0*Hr3(0,0,0)*dm + 152.D0/3.D0*Hr3(0,0,1)
     &     + 568.D0*Hr3(0,0,1)*x + 496.D0/3.D0*Hr3(0,0,1)*dp + 64.D0/3.D
     &    0*Hr3(0,0,1)*dm - 80.D0/3.D0*Hr3(0,1,0) - 80.D0/3.D0*Hr3(0,1,
     &    0)*x + 352.D0/3.D0*Hr3(0,1,0)*dm - 304.D0/3.D0*Hr3(1,0,0) - 
     &    688.D0/3.D0*Hr3(1,0,0)*x + 992.D0/3.D0*Hr3(1,0,0)*dm - 576.D0
     &    *Hr4(-1,-1,0,0) + 576.D0*Hr4(-1,-1,0,0)*x + 1152.D0*Hr4(-1,-1
     &    ,0,0)*dp - 896.D0*Hr4(-1,-1,0,1) + 896.D0*Hr4(-1,-1,0,1)*x + 
     &    1792.D0*Hr4(-1,-1,0,1)*dp - 64.D0*Hr4(-1,0,-1,0) + 64.D0*Hr4(
     &    -1,0,-1,0)*x + 128.D0*Hr4(-1,0,-1,0)*dp + 272.D0*Hr4(-1,0,0,0
     &    ) - 272.D0*Hr4(-1,0,0,0)*x - 544.D0*Hr4(-1,0,0,0)*dp + 512.D0
     &    *Hr4(-1,0,0,1) - 512.D0*Hr4(-1,0,0,1)*x - 1024.D0*Hr4(-1,0,0,
     &    1)*dp + 32.D0*Hr4(-1,0,1,0) - 32.D0*Hr4(-1,0,1,0)*x - 64.D0*
     &    Hr4(-1,0,1,0)*dp + 128.D0*Hr4(0,-1,-1,0) + 320.D0*Hr4(0,-1,-1
     &    ,0)*x + 128.D0*Hr4(0,-1,-1,0)*dp - 448.D0*Hr4(0,-1,-1,0)*dm
     &     + 304.D0*Hr4(0,-1,0,0) - 464.D0*Hr4(0,-1,0,0)*x - 672.D0*
     &    Hr4(0,-1,0,0)*dp + 160.D0*Hr4(0,-1,0,0)*dm + 544.D0*Hr4(0,-1,
     &    0,1) - 416.D0*Hr4(0,-1,0,1)*x - 960.D0*Hr4(0,-1,0,1)*dp - 128.
     &    D0*Hr4(0,-1,0,1)*dm + 128.D0*Hr4(0,0,-1,0) - 128.D0*Hr4(0,0,
     &    -1,0)*x - 160.D0*Hr4(0,0,-1,0)*dp - 32.D0*Hr4(0,0,-1,0)*dm - 
     &    224.D0*Hr4(0,0,0,0) + 160.D0*Hr4(0,0,0,0)*dp + 160.D0*Hr4(0,0
     &    ,0,0)*dm - 224.D0*Hr4(0,0,0,1) + 96.D0*Hr4(0,0,0,1)*x + 320.D0
     &    *Hr4(0,0,0,1)*dp + 32.D0*Hr4(0,0,0,1)*dm - 32.D0*Hr4(0,0,1,0)
     &     + 32.D0*Hr4(0,0,1,0)*dp + 32.D0*Hr4(0,0,1,0)*dm - 48.D0*Hr4(
     &    0,1,0,0) - 48.D0*Hr4(0,1,0,0)*x + 160.D0*Hr4(0,1,0,0)*dm + 
     &    256.D0*Hr4(1,0,-1,0) + 256.D0*Hr4(1,0,-1,0)*x - 512.D0*Hr4(1,
     &    0,-1,0)*dm - 176.D0*Hr4(1,0,0,0) - 176.D0*Hr4(1,0,0,0)*x + 
     &    352.D0*Hr4(1,0,0,0)*dm + 128.D0*Hr4(1,0,0,1) + 128.D0*Hr4(1,0
     &    ,0,1)*x - 256.D0*Hr4(1,0,0,1)*dm - 128.D0*Hr4(1,1,0,0) - 128.D
     &    0*Hr4(1,1,0,0)*x + 256.D0*Hr4(1,1,0,0)*dm )
      dPqq2M = dPqq2M + cf**3 * (  - 302.D0 + 302.D0*x + 252.D0*z4 - 
     &    28.D0*z4*x - 280.D0*dp*z4 + 72.D0*dm*z4 + 48.D0*z3 + 336.D0*
     &    z3*x + 144.D0*z3*dp + 204.D0*z2 + 12.D0*z2*x + 576.D0*Hr1(-1)
     &    *z3 - 576.D0*Hr1(-1)*z3*x - 1152.D0*Hr1(-1)*z3*dp - 192.D0*
     &    Hr1(-1)*z2 - 384.D0*Hr1(-1)*z2*x - 192.D0*Hr1(-1)*z2*dp - 328.
     &    D0*Hr1(0) + 464.D0*Hr1(0)*x + 48.D0*Hr1(0)*dp - 6.D0*Hr1(0)*
     &    dm - 224.D0*Hr1(0)*z3 + 288.D0*Hr1(0)*z3*x + 416.D0*Hr1(0)*z3
     &    *dp + 32.D0*Hr1(0)*z3*dm + 192.D0*Hr1(0)*z2 + 544.D0*Hr1(0)*
     &    z2*x + 48.D0*Hr1(0)*z2*dp - 48.D0*Hr1(0)*z2*dm - 400.D0*Hr1(1
     &    ) + 400.D0*Hr1(1)*x - 192.D0*Hr1(1)*z3 - 192.D0*Hr1(1)*z3*x
     &     + 384.D0*Hr1(1)*z3*dm - 96.D0*Hr1(1)*z2 + 96.D0*Hr1(1)*z2*x
     &     - 768.D0*Hr2(-1,-1)*z2 + 768.D0*Hr2(-1,-1)*z2*x + 1536.D0*
     &    Hr2(-1,-1)*z2*dp - 96.D0*Hr2(-1,0) - 96.D0*Hr2(-1,0)*x + 640.D
     &    0*Hr2(-1,0)*z2 - 640.D0*Hr2(-1,0)*z2*x - 1280.D0*Hr2(-1,0)*z2
     &    *dp + 416.D0*Hr2(0,-1)*z2 - 608.D0*Hr2(0,-1)*z2*x - 1024.D0*
     &    Hr2(0,-1)*z2*dp + 192.D0*Hr2(0,-1)*z2*dm - 172.D0*Hr2(0,0) + 
     &    4.D0*Hr2(0,0)*x + 52.D0*Hr2(0,0)*dm - 240.D0*Hr2(0,0)*z2 + 
     &    240.D0*Hr2(0,0)*z2*x + 448.D0*Hr2(0,0)*z2*dp - 300.D0*Hr2(0,1
     &    ) - 12.D0*Hr2(0,1)*x + 96.D0*Hr2(0,1)*z2 - 96.D0*Hr2(0,1)*z2*
     &    x - 192.D0*Hr2(0,1)*z2*dp - 144.D0*Hr2(1,0) + 144.D0*Hr2(1,0)
     &    *x + 192.D0*Hr3(-1,-1,0) + 192.D0*Hr3(-1,-1,0)*x + 272.D0*
     &    Hr3(-1,0,0) + 176.D0*Hr3(-1,0,0)*x - 96.D0*Hr3(-1,0,0)*dp + 
     &    288.D0*Hr3(-1,0,1) + 480.D0*Hr3(-1,0,1)*x + 192.D0*Hr3(-1,0,1
     &    )*dp + 64.D0*Hr3(0,-1,0) - 128.D0*Hr3(0,-1,0)*x - 96.D0*Hr3(0
     &    ,-1,0)*dp - 168.D0*Hr3(0,0,0) - 168.D0*Hr3(0,0,0)*x + 144.D0*
     &    Hr3(0,0,0)*dp - 128.D0*Hr3(0,0,1) - 544.D0*Hr3(0,0,1)*x - 96.D
     &    0*Hr3(0,0,1)*dp + 8.D0*Hr3(0,1,0) - 56.D0*Hr3(0,1,0)*x - 96.D0
     &    *Hr3(0,1,0)*dm + 96.D0*Hr3(1,0,0) + 96.D0*Hr3(1,0,0)*x - 192.D
     &    0*Hr3(1,0,0)*dm + 640.D0*Hr4(-1,-1,0,0) - 640.D0*Hr4(-1,-1,0,
     &    0)*x - 1280.D0*Hr4(-1,-1,0,0)*dp + 768.D0*Hr4(-1,-1,0,1) - 
     &    768.D0*Hr4(-1,-1,0,1)*x - 1536.D0*Hr4(-1,-1,0,1)*dp + 128.D0*
     &    Hr4(-1,0,-1,0) - 128.D0*Hr4(-1,0,-1,0)*x - 256.D0*Hr4(-1,0,-1
     &    ,0)*dp - 352.D0*Hr4(-1,0,0,0) + 352.D0*Hr4(-1,0,0,0)*x + 704.D
     &    0*Hr4(-1,0,0,0)*dp - 512.D0*Hr4(-1,0,0,1) + 512.D0*Hr4(-1,0,0
     &    ,1)*x + 1024.D0*Hr4(-1,0,0,1)*dp - 64.D0*Hr4(-1,0,1,0) + 64.D0
     &    *Hr4(-1,0,1,0)*x + 128.D0*Hr4(-1,0,1,0)*dp - 64.D0*Hr4(0,-1,
     &    -1,0) - 320.D0*Hr4(0,-1,-1,0)*x - 256.D0*Hr4(0,-1,-1,0)*dp + 
     &    384.D0*Hr4(0,-1,-1,0)*dm - 384.D0*Hr4(0,-1,0,0) + 576.D0*Hr4(
     &    0,-1,0,0)*x + 832.D0*Hr4(0,-1,0,0)*dp - 192.D0*Hr4(0,-1,0,0)*
     &    dm - 448.D0*Hr4(0,-1,0,1) + 448.D0*Hr4(0,-1,0,1)*x + 896.D0*
     &    Hr4(0,-1,0,1)*dp - 96.D0*Hr4(0,0,-1,0) + 224.D0*Hr4(0,0,-1,0)
     &    *x + 192.D0*Hr4(0,0,-1,0)*dp - 64.D0*Hr4(0,0,-1,0)*dm + 176.D0
     &    *Hr4(0,0,0,0) - 112.D0*Hr4(0,0,0,0)*x - 192.D0*Hr4(0,0,0,0)*
     &    dp - 64.D0*Hr4(0,0,0,0)*dm + 144.D0*Hr4(0,0,0,1) - 240.D0*
     &    Hr4(0,0,0,1)*x - 384.D0*Hr4(0,0,0,1)*dp + 64.D0*Hr4(0,0,0,1)*
     &    dm - 64.D0*Hr4(0,0,1,0) - 128.D0*Hr4(0,0,1,0)*x - 64.D0*Hr4(0
     &    ,0,1,0)*dp + 128.D0*Hr4(0,0,1,0)*dm - 64.D0*Hr4(0,0,1,1) - 64.
     &    D0*Hr4(0,0,1,1)*x + 128.D0*Hr4(0,0,1,1)*dm - 80.D0*Hr4(0,1,0,
     &    0) - 80.D0*Hr4(0,1,0,0)*x + 64.D0*Hr4(0,1,0,0)*dm - 64.D0*
     &    Hr4(0,1,0,1) - 64.D0*Hr4(0,1,0,1)*x + 128.D0*Hr4(0,1,0,1)*dm
     &     - 64.D0*Hr4(0,1,1,0) - 64.D0*Hr4(0,1,1,0)*x + 128.D0*Hr4(0,1
     &    ,1,0)*dm - 128.D0*Hr4(1,0,-1,0) - 128.D0*Hr4(1,0,-1,0)*x + 
     &    256.D0*Hr4(1,0,-1,0)*dm + 64.D0*Hr4(1,0,0,0) + 64.D0*Hr4(1,0,
     &    0,0)*x - 128.D0*Hr4(1,0,0,0)*dm - 128.D0*Hr4(1,0,0,1) - 128.D0
     &    *Hr4(1,0,0,1)*x + 256.D0*Hr4(1,0,0,1)*dm - 64.D0*Hr4(1,0,1,0)
     &     - 64.D0*Hr4(1,0,1,0)*x + 128.D0*Hr4(1,0,1,0)*dm )
      dPqq2M = dPqq2M + nf*cf*ca * (  - 38614.D0/27.D0 + 13150.D0/9.D0
     &    *x + 234.D0*z4 + 214.D0*z4*x - 836.D0/27.D0*dm + 1496.D0/3.D0
     &    *z3 - 1408.D0/3.D0*z3*x - 16.D0*z3*dp - 48.D0*z3*dm + 92.D0/3.
     &    D0*z2 - 560.D0/9.D0*z2*x + 160.D0/9.D0*z2*dp + 160.D0/9.D0*z2
     &    *dm - 56.D0/3.D0*Hr1(-1)*z2 + 8.D0/3.D0*Hr1(-1)*z2*x + 64.D0/
     &    3.D0*Hr1(-1)*z2*dp - 9944.D0/27.D0*Hr1(0) - 34880.D0/27.D0*
     &    Hr1(0)*x - 1336.D0/27.D0*Hr1(0)*dm + 144.D0*Hr1(0)*z3 + 208.D0
     &    *Hr1(0)*z3*x - 124.D0/3.D0*Hr1(0)*z2 - 212.D0/3.D0*Hr1(0)*z2*
     &    x - 16.D0/3.D0*Hr1(0)*z2*dp + 32.D0/3.D0*Hr1(0)*z2*dm + 4304.D
     &    0/9.D0*Hr1(1) - 4304.D0/9.D0*Hr1(1)*x + 40.D0*Hr1(1)*z2 - 40.D
     &    0*Hr1(1)*z2*x + 2600.D0/9.D0*Hr2(-1,0) + 2920.D0/9.D0*Hr2(-1,
     &    0)*x + 320.D0/9.D0*Hr2(-1,0)*dp - 16.D0*Hr2(0,-1)*z2 + 16.D0*
     &    Hr2(0,-1)*z2*x - 1604.D0/9.D0*Hr2(0,0) + 508.D0*Hr2(0,0)*x - 
     &    160.D0/9.D0*Hr2(0,0)*dp - 112.D0/3.D0*Hr2(0,0)*dm - 8.D0*Hr2(
     &    0,0)*z2 + 8.D0*Hr2(0,0)*z2*x + 2324.D0/9.D0*Hr2(0,1) + 560.D0/
     &    9.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*z2 + 16.D0*Hr2(0,1)*z2*x + 
     &    136.D0*Hr2(1,0) - 136.D0*Hr2(1,0)*x + 260.D0/3.D0*Hr2(1,1) - 
     &    260.D0/3.D0*Hr2(1,1)*x - 80.D0*Hr3(-1,-1,0) - 80.D0*Hr3(-1,-1
     &    ,0)*x - 424.D0/3.D0*Hr3(-1,0,0) - 392.D0/3.D0*Hr3(-1,0,0)*x
     &     + 32.D0/3.D0*Hr3(-1,0,0)*dp - 64.D0/3.D0*Hr3(-1,0,1) - 128.D0
     &    /3.D0*Hr3(-1,0,1)*x - 64.D0/3.D0*Hr3(-1,0,1)*dp + 200.D0/3.D0
     &    *Hr3(0,-1,0) - 440.D0/3.D0*Hr3(0,-1,0)*x + 32.D0/3.D0*Hr3(0,
     &    -1,0)*dp - 272.D0/3.D0*Hr3(0,0,0) - 592.D0/3.D0*Hr3(0,0,0)*x
     &     - 32.D0/3.D0*Hr3(0,0,0)*dp - 32.D0/3.D0*Hr3(0,0,0)*dm + 108.D
     &    0*Hr3(0,0,1) + 212.D0/3.D0*Hr3(0,0,1)*x + 32.D0/3.D0*Hr3(0,0,
     &    1)*dp - 16.D0/3.D0*Hr3(0,0,1)*dm + 48.D0*Hr3(0,1,0) + 8.D0*
     &    Hr3(0,1,0)*x + 152.D0/3.D0*Hr3(0,1,1) - 88.D0/3.D0*Hr3(0,1,1)
     &    *x - 140.D0*Hr3(1,0,0) + 124.D0*Hr3(1,0,0)*x + 16.D0*Hr3(1,0,
     &    0)*dm + 80.D0*Hr3(1,1,0) - 80.D0*Hr3(1,1,0)*x + 80.D0*Hr3(1,1
     &    ,1) - 80.D0*Hr3(1,1,1)*x - 32.D0*Hr4(0,-1,-1,0) + 32.D0*Hr4(0
     &    ,-1,-1,0)*x - 48.D0*Hr4(0,-1,0,0) + 48.D0*Hr4(0,-1,0,0)*x - 
     &    16.D0*Hr4(0,0,-1,0) + 80.D0*Hr4(0,0,-1,0)*x - 48.D0*Hr4(0,0,0
     &    ,0) + 64.D0*Hr4(0,0,0,0)*x - 8.D0*Hr4(0,0,0,1) - 8.D0*Hr4(0,0
     &    ,0,1)*x + 16.D0*Hr4(0,0,1,1) + 16.D0*Hr4(0,0,1,1)*x - 56.D0*
     &    Hr4(0,1,0,0) - 56.D0*Hr4(0,1,0,0)*x + 32.D0*Hr4(0,1,1,0) + 32.
     &    D0*Hr4(0,1,1,0)*x + 32.D0*Hr4(0,1,1,1) + 32.D0*Hr4(0,1,1,1)*x
     &     )
      dPqq2M = dPqq2M + nf*cf**2 * ( 61.D0 - 73.D0/3.D0*x - 148.D0*z4
     &     - 148.D0*z4*x - 110.D0/3.D0*dm - 656.D0/3.D0*z3 + 976.D0/3.D0
     &    *z3*x + 32.D0*z3*dp + 160.D0/3.D0*z3*dm + 2464.D0/9.D0*z2 + 
     &    812.D0/9.D0*z2*x - 320.D0/9.D0*z2*dp + 64.D0/3.D0*Hr1(-1)*z2
     &     - 64.D0/3.D0*Hr1(-1)*z2*x - 128.D0/3.D0*Hr1(-1)*z2*dp - 392.D
     &    0/9.D0*Hr1(0) + 4580.D0/9.D0*Hr1(0)*x + 20.D0/3.D0*Hr1(0)*dm
     &     - 112.D0*Hr1(0)*z3 - 112.D0*Hr1(0)*z3*x + 448.D0/3.D0*Hr1(0)
     &    *z2 - 80.D0*Hr1(0)*z2*x + 32.D0/3.D0*Hr1(0)*z2*dp + 32.D0/3.D0
     &    *Hr1(0)*z2*dm - 1564.D0/3.D0*Hr1(1) + 1564.D0/3.D0*Hr1(1)*x
     &     + 128.D0/9.D0*Hr2(-1,0) - 512.D0/9.D0*Hr2(-1,0)*x - 640.D0/9.
     &    D0*Hr2(-1,0)*dp - 1136.D0/9.D0*Hr2(0,0) - 124.D0/3.D0*Hr2(0,0
     &    )*x + 320.D0/9.D0*Hr2(0,0)*dp + 16.D0*Hr2(0,0)*dm + 96.D0*
     &    Hr2(0,0)*z2 + 96.D0*Hr2(0,0)*z2*x - 2336.D0/9.D0*Hr2(0,1) - 
     &    812.D0/9.D0*Hr2(0,1)*x - 320.D0/9.D0*Hr2(0,1)*dm - 368.D0/9.D0
     &    *Hr2(1,0) + 688.D0/9.D0*Hr2(1,0)*x - 320.D0/9.D0*Hr2(1,0)*dm
     &     - 52.D0*Hr2(1,1) + 52.D0*Hr2(1,1)*x + 32.D0/3.D0*Hr3(-1,0,0)
     &     - 32.D0/3.D0*Hr3(-1,0,0)*x - 64.D0/3.D0*Hr3(-1,0,0)*dp - 64.D
     &    0/3.D0*Hr3(-1,0,1) + 64.D0/3.D0*Hr3(-1,0,1)*x + 128.D0/3.D0*
     &    Hr3(-1,0,1)*dp + 32.D0/3.D0*Hr3(0,-1,0) - 32.D0/3.D0*Hr3(0,-1
     &    ,0)*x - 64.D0/3.D0*Hr3(0,-1,0)*dp - 232.D0/3.D0*Hr3(0,0,0) + 
     &    120.D0*Hr3(0,0,0)*x + 64.D0/3.D0*Hr3(0,0,0)*dp - 32.D0/3.D0*
     &    Hr3(0,0,0)*dm - 416.D0/3.D0*Hr3(0,0,1) + 80.D0*Hr3(0,0,1)*x
     &     - 64.D0/3.D0*Hr3(0,0,1)*dp - 64.D0/3.D0*Hr3(0,0,1)*dm - 160.D
     &    0/3.D0*Hr3(0,1,0) + 368.D0/3.D0*Hr3(0,1,0)*x - 64.D0/3.D0*
     &    Hr3(0,1,0)*dm - 32.D0*Hr3(0,1,1) + 48.D0*Hr3(0,1,1)*x + 184.D0
     &    /3.D0*Hr3(1,0,0) - 56.D0/3.D0*Hr3(1,0,0)*x - 128.D0/3.D0*Hr3(
     &    1,0,0)*dm - 80.D0*Hr3(1,1,0) + 80.D0*Hr3(1,1,0)*x - 80.D0*
     &    Hr3(1,1,1) + 80.D0*Hr3(1,1,1)*x - 64.D0*Hr4(0,0,0,0) - 64.D0*
     &    Hr4(0,0,0,0)*x - 96.D0*Hr4(0,0,0,1) - 96.D0*Hr4(0,0,0,1)*x - 
     &    64.D0*Hr4(0,0,1,0) - 64.D0*Hr4(0,0,1,0)*x - 16.D0*Hr4(0,0,1,1
     &    ) - 16.D0*Hr4(0,0,1,1)*x + 16.D0*Hr4(0,1,0,0) + 16.D0*Hr4(0,1
     &    ,0,0)*x - 32.D0*Hr4(0,1,1,0) - 32.D0*Hr4(0,1,1,0)*x - 32.D0*
     &    Hr4(0,1,1,1) - 32.D0*Hr4(0,1,1,1)*x )
      dPqq2M = dPqq2M + nf2*cf * ( 800.D0/27.D0 - 784.D0/27.D0*x - 
     &    16.D0/27.D0*dm + 16.D0/3.D0*z3 + 16.D0/3.D0*z3*x + 80.D0/9.D0
     &    *z2 - 64.D0/9.D0*z2*x + 176.D0/9.D0*Hr1(0) - 16.D0*Hr1(0)*x
     &     + 80.D0/27.D0*Hr1(0)*dm + 32.D0/3.D0*Hr1(0)*z2 + 32.D0/3.D0*
     &    Hr1(0)*z2*x + 16.D0/9.D0*Hr1(1) - 16.D0/9.D0*Hr1(1)*x + 176.D0
     &    /9.D0*Hr2(0,0) + 128.D0/9.D0*Hr2(0,0)*x + 16.D0/9.D0*Hr2(0,0)
     &    *dm - 80.D0/9.D0*Hr2(0,1) + 64.D0/9.D0*Hr2(0,1)*x + 40.D0/3.D0
     &    *Hr2(1,1) - 40.D0/3.D0*Hr2(1,1)*x + 16.D0/3.D0*Hr3(0,0,0) + 
     &    16.D0/3.D0*Hr3(0,0,0)*x - 32.D0/3.D0*Hr3(0,0,1) - 32.D0/3.D0*
     &    Hr3(0,0,1)*x + 16.D0/3.D0*Hr3(0,1,1) + 16.D0/3.D0*Hr3(0,1,1)*
     &    x )
*
* ...The soft (`+'-distribution) part of the splitting function
*
       A3 =
     ,      ca**2*cf * ( 490.D0/3.D0 - 1072.D0/9.D0*z2 + 88.D0/3.D0*z3
     ,                 + 176.D0/5.D0*z2**2 )
     ,    + ca*cf*nf * ( - 836./27.D0 + 160./9.D0*z2 - 112./3.D0*z3 )
     ,    + cf**2*nf * ( - 110./3.D0 + 32.*z3 ) - cf*nf2 * 16./27.D0
*
       DPQQ2S = DM * A3
*
* ...The regular piece of the splitting function
*
       DP2QQMA = dPqq2M - DPQQ2S
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The singular (soft) piece of DeltaP_qq^(2)
*
       FUNCTION DP2QQB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P2SOFT / A3
*
       DP2QQB  = A3/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The 'local' piece of DeltaP_qq^(2)
*
       FUNCTION DP2QQC (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,
     ,             Z4 = 1.0823 23233 71113 81915 D0,
     ,             Z5 = 1.0369 27755 14336 99263 D0 )
*
       COMMON / P2SOFT / A3
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
* ...The coefficient of delta(1-x)
*
       P2DELT =
     &     + 29.D0/2.D0*cf**3
     &     + 151.D0/4.D0*ca*cf**2
     &     - 1657.D0/36.D0*ca**2*cf
     &     - 240.D0*z5*cf**3
     &     + 120.D0*z5*ca*cf**2
     &     + 40.D0*z5*ca**2*cf
     &     + 68.D0*z3*cf**3
     &     + 844.D0/3.D0*z3*ca*cf**2
     &     - 1552.D0/9.D0*z3*ca**2*cf
     &     + 18.D0*z2*cf**3
     &     - 410.D0/3.D0*z2*ca*cf**2
     &     + 4496.D0/27.D0*z2*ca**2*cf
     &     - 32.D0*z2*z3*cf**3
     &     + 16.D0*z2*z3*ca*cf**2
     &     + 288.D0/5.D0*z2**2*cf**3
     &     - 988.D0/15.D0*z2**2*ca*cf**2
     &     - 2.D0*z2**2*ca**2*cf
*
     &     - 1336.D0/27.D0*z2*ca*cf*nf
     &     + 4.D0/5.D0*z2**2*ca*cf*nf
     &     + 200.D0/9.D0*z3*ca*cf*nf
     &     + 20.D0*ca*cf*nf
     &     + 20.D0/3.D0*z2*cf**2*nf
     &     + 232.D0/15.D0*z2**2*cf**2*nf
     &     - 136.D0/3.D0*z3*cf**2*nf
     &     - 23.D0*cf**2*nf
     &     + 80.D0/27.D0*z2*cf*nf2
     &     - 16.D0/9.D0*z3*cf*nf2
     &     - 17.D0/9.D0*cf*nf2
*
       DP2QQC = LOG (1.D0-Y) * A3 + P2DELT
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..3-loop pure-singlet quark splitting function  DeltaP_ps^(2).
*
       FUNCTION DP2PSMA (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,
     ,             Z4 = 1.0823 23233 71113 81915 D0 )
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       NF2 = NF*NF
*
* ...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
*
      dPps2M =
     &  + nf*cf*ca * (  - 35512.D0/27.D0 + 35512.D0/27.D0*x + 234.D0*z4
     &     + 214.D0*z4*x + 1400.D0/3.D0*z3 - 1456.D0/3.D0*z3*x + 364.D0/
     &    9.D0*z2 - 536.D0/9.D0*z2*x - 8.D0*Hr1(-1)*z2 - 8.D0*Hr1(-1)*
     &    z2*x - 9028.D0/27.D0*Hr1(0) - 35596.D0/27.D0*Hr1(0)*x + 144.D0
     &    *Hr1(0)*z3 + 208.D0*Hr1(0)*z3*x - 116.D0/3.D0*Hr1(0)*z2 - 188.
     &    D0/3.D0*Hr1(0)*z2*x + 4256.D0/9.D0*Hr1(1) - 4256.D0/9.D0*Hr1(
     &    1)*x + 40.D0*Hr1(1)*z2 - 40.D0*Hr1(1)*z2*x + 296.D0*Hr2(-1,0)
     &     + 296.D0*Hr2(-1,0)*x - 16.D0*Hr2(0,-1)*z2 + 16.D0*Hr2(0,-1)*
     &    z2*x - 1708.D0/9.D0*Hr2(0,0) + 4604.D0/9.D0*Hr2(0,0)*x - 8.D0
     &    *Hr2(0,0)*z2 + 8.D0*Hr2(0,0)*z2*x + 2300.D0/9.D0*Hr2(0,1) +
     &    536.D0/9.D0*Hr2(0,1)*x + 16.D0*Hr2(0,1)*z2 + 16.D0*Hr2(0,1)*
     &    z2*x + 136.D0*Hr2(1,0) - 136.D0*Hr2(1,0)*x + 260.D0/3.D0*Hr2(
     &    1,1) - 260.D0/3.D0*Hr2(1,1)*x - 80.D0*Hr3(-1,-1,0) - 80.D0*
     &    Hr3(-1,-1,0)*x - 136.D0*Hr3(-1,0,0) - 136.D0*Hr3(-1,0,0)*x -
     &    32.D0*Hr3(-1,0,1) - 32.D0*Hr3(-1,0,1)*x + 72.D0*Hr3(0,-1,0)
     &     - 152.D0*Hr3(0,-1,0)*x - 304.D0/3.D0*Hr3(0,0,0) - 592.D0/3.D0
     &    *Hr3(0,0,0)*x + 332.D0/3.D0*Hr3(0,0,1) + 188.D0/3.D0*Hr3(0,0,
     &    1)*x + 48.D0*Hr3(0,1,0) + 8.D0*Hr3(0,1,0)*x + 152.D0/3.D0*
     &    Hr3(0,1,1) - 88.D0/3.D0*Hr3(0,1,1)*x - 132.D0*Hr3(1,0,0) +
     &    132.D0*Hr3(1,0,0)*x + 80.D0*Hr3(1,1,0) - 80.D0*Hr3(1,1,0)*x
     &     + 80.D0*Hr3(1,1,1) - 80.D0*Hr3(1,1,1)*x - 32.D0*Hr4(0,-1,-1,
     &    0) + 32.D0*Hr4(0,-1,-1,0)*x - 48.D0*Hr4(0,-1,0,0) + 48.D0*
     &    Hr4(0,-1,0,0)*x - 16.D0*Hr4(0,0,-1,0) + 80.D0*Hr4(0,0,-1,0)*x
     &     - 48.D0*Hr4(0,0,0,0) + 64.D0*Hr4(0,0,0,0)*x - 8.D0*Hr4(0,0,0
     &    ,1) - 8.D0*Hr4(0,0,0,1)*x + 16.D0*Hr4(0,0,1,1) + 16.D0*Hr4(0,
     &    0,1,1)*x - 56.D0*Hr4(0,1,0,0) - 56.D0*Hr4(0,1,0,0)*x + 32.D0*
     &    Hr4(0,1,1,0) + 32.D0*Hr4(0,1,1,0)*x + 32.D0*Hr4(0,1,1,1) + 32.
     &    D0*Hr4(0,1,1,1)*x )
      dPps2M = dPps2M + nf*cf**2 * (  - 48.D0 + 48.D0*x - 148.D0*z4
     &     - 148.D0*z4*x - 176.D0*z3 + 336.D0*z3*x + 256.D0*z2 + 108.D0
     &    *z2*x - 100.D0*Hr1(0) + 456.D0*Hr1(0)*x - 112.D0*Hr1(0)*z3 -
     &    112.D0*Hr1(0)*z3*x + 160.D0*Hr1(0)*z2 - 80.D0*Hr1(0)*z2*x -
     &    500.D0*Hr1(1) + 500.D0*Hr1(1)*x - 128.D0*Hr2(0,0) - 100.D0*
     &    Hr2(0,0)*x + 96.D0*Hr2(0,0)*z2 + 96.D0*Hr2(0,0)*z2*x - 256.D0
     &    *Hr2(0,1) - 108.D0*Hr2(0,1)*x - 48.D0*Hr2(1,0) + 48.D0*Hr2(1,
     &    0)*x - 52.D0*Hr2(1,1) + 52.D0*Hr2(1,1)*x - 80.D0*Hr3(0,0,0)
     &     + 96.D0*Hr3(0,0,0)*x - 160.D0*Hr3(0,0,1) + 80.D0*Hr3(0,0,1)*
     &    x - 64.D0*Hr3(0,1,0) + 112.D0*Hr3(0,1,0)*x - 32.D0*Hr3(0,1,1)
     &     + 48.D0*Hr3(0,1,1)*x + 40.D0*Hr3(1,0,0) - 40.D0*Hr3(1,0,0)*x
     &     - 80.D0*Hr3(1,1,0) + 80.D0*Hr3(1,1,0)*x - 80.D0*Hr3(1,1,1)
     &     + 80.D0*Hr3(1,1,1)*x - 64.D0*Hr4(0,0,0,0) - 64.D0*Hr4(0,0,0,
     &    0)*x - 96.D0*Hr4(0,0,0,1) - 96.D0*Hr4(0,0,0,1)*x - 64.D0*Hr4(
     &    0,0,1,0) - 64.D0*Hr4(0,0,1,0)*x - 16.D0*Hr4(0,0,1,1) - 16.D0*
     &    Hr4(0,0,1,1)*x + 16.D0*Hr4(0,1,0,0) + 16.D0*Hr4(0,1,0,0)*x -
     &    32.D0*Hr4(0,1,1,0) - 32.D0*Hr4(0,1,1,0)*x - 32.D0*Hr4(0,1,1,1
     &    ) - 32.D0*Hr4(0,1,1,1)*x )
      dPps2M = dPps2M + nf2*cf * ( 688.D0/27.D0 - 688.D0/27.D0*x +
     &    16.D0/3.D0*z3 + 16.D0/3.D0*z3*x + 80.D0/9.D0*z2 - 64.D0/9.D0*
     &    z2*x + 520.D0/27.D0*Hr1(0) - 344.D0/27.D0*Hr1(0)*x + 32.D0/3.D
     &    0*Hr1(0)*z2 + 32.D0/3.D0*Hr1(0)*z2*x + 16.D0/9.D0*Hr1(1) - 16.
     &    D0/9.D0*Hr1(1)*x + 184.D0/9.D0*Hr2(0,0) + 136.D0/9.D0*Hr2(0,0
     &    )*x - 80.D0/9.D0*Hr2(0,1) + 64.D0/9.D0*Hr2(0,1)*x + 40.D0/3.D0
     &    *Hr2(1,1) - 40.D0/3.D0*Hr2(1,1)*x + 16.D0/3.D0*Hr3(0,0,0) +
     &    16.D0/3.D0*Hr3(0,0,0)*x - 32.D0/3.D0*Hr3(0,0,1) - 32.D0/3.D0*
     &    Hr3(0,0,1)*x + 16.D0/3.D0*Hr3(0,1,1) + 16.D0/3.D0*Hr3(0,1,1)*
     &    x )
*
       DP2PSMA = dPps2M
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 3-loop gluon->quark splitting function DeltaP_qg^(2)
*
       FUNCTION DP2QGMA (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,
     ,             Z4 = 1.0823 23233 71113 81915 D0 )
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       NF2 = NF*NF
*
       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
*
      dPqg2M =
     &  + nf*ca**2 * (  - 32198.D0/9.D0 + 32692.D0/9.D0*x + 495.D0*z4
     &     + 538.D0*z4*x + 1260.D0*z3 - 1072.D0*z3*x + 2932.D0/9.D0*z2
     &     + 16.D0/3.D0*z2*x + 24.D0*z2*x**2 - 40.D0*Hr1(-1)*z3 - 80.D0
     &    *Hr1(-1)*z3*x - 104.D0*Hr1(-1)*z2 - 104.D0*Hr1(-1)*z2*x - 
     &    22936.D0/27.D0*Hr1(0) - 81286.D0/27.D0*Hr1(0)*x + 288.D0*Hr1(
     &    0)*z3 + 384.D0*Hr1(0)*z3*x + 28.D0*Hr1(0)*z2 - 124.D0*Hr1(0)*
     &    z2*x + 10060.D0/27.D0*Hr1(1) - 9110.D0/27.D0*Hr1(1)*x + 248.D0
     &    *Hr1(1)*z3 - 496.D0*Hr1(1)*z3*x + 364.D0/3.D0*Hr1(1)*z2 - 320.
     &    D0/3.D0*Hr1(1)*z2*x + 16.D0*Hr2(-1,-1)*z2 + 32.D0*Hr2(-1,-1)*
     &    z2*x + 6208.D0/9.D0*Hr2(-1,0) + 5672.D0/9.D0*Hr2(-1,0)*x + 24.
     &    D0*Hr2(-1,0)*x**2 - 8.D0*Hr2(0,-1)*z2 + 112.D0*Hr2(0,-1)*z2*x
     &     - 300.D0*Hr2(0,0) + 9424.D0/9.D0*Hr2(0,0)*x - 24.D0*Hr2(0,0)
     &    *x**2 + 28.D0*Hr2(0,0)*z2 - 24.D0*Hr2(0,0)*z2*x + 364.D0*Hr2(
     &    0,1) - 16.D0/3.D0*Hr2(0,1)*x + 40.D0*Hr2(0,1)*z2 + 48.D0*Hr2(
     &    0,1)*z2*x + 288.D0*Hr2(1,0) - 280.D0*Hr2(1,0)*x + 80.D0*Hr2(1
     &    ,0)*z2 - 160.D0*Hr2(1,0)*z2*x - 16.D0/9.D0*Hr2(1,1) - 520.D0/
     &    9.D0*Hr2(1,1)*x + 16.D0*Hr2(1,1)*z2 - 32.D0*Hr2(1,1)*z2*x - 
     &    336.D0*Hr3(-1,-1,0) - 336.D0*Hr3(-1,-1,0)*x - 832.D0/3.D0*
     &    Hr3(-1,0,0) - 920.D0/3.D0*Hr3(-1,0,0)*x - 64.D0*Hr3(-1,0,1)
     &     - 64.D0*Hr3(-1,0,1)*x + 168.D0*Hr3(0,-1,0) - 216.D0*Hr3(0,-1
     &    ,0)*x - 656.D0/3.D0*Hr3(0,0,0) - 640.D0/3.D0*Hr3(0,0,0)*x + 
     &    140.D0*Hr3(0,0,1) + 124.D0*Hr3(0,0,1)*x + 96.D0*Hr3(0,1,0) - 
     &    32.D0*Hr3(0,1,0)*x + 112.D0*Hr3(0,1,1) - 160.D0*Hr3(0,1,1)*x
     &     - 848.D0/3.D0*Hr3(1,0,0) + 892.D0/3.D0*Hr3(1,0,0)*x + 140.D0/
     &    3.D0*Hr3(1,0,1) - 184.D0/3.D0*Hr3(1,0,1)*x + 628.D0/3.D0*Hr3(
     &    1,1,0) - 584.D0/3.D0*Hr3(1,1,0)*x - 44.D0/3.D0*Hr3(1,1,1) + 
     &    88.D0/3.D0*Hr3(1,1,1)*x + 96.D0*Hr4(-1,-1,-1,0) + 192.D0*Hr4(
     &    -1,-1,-1,0)*x + 32.D0*Hr4(-1,-1,0,0) + 64.D0*Hr4(-1,-1,0,0)*x
     &     + 32.D0*Hr4(-1,-1,0,1) + 64.D0*Hr4(-1,-1,0,1)*x - 16.D0*Hr4(
     &    -1,0,-1,0) - 32.D0*Hr4(-1,0,-1,0)*x - 24.D0*Hr4(-1,0,0,0) - 
     &    48.D0*Hr4(-1,0,0,0)*x - 32.D0*Hr4(-1,0,0,1) - 64.D0*Hr4(-1,0,
     &    0,1)*x + 32.D0*Hr4(-1,0,1,1) + 64.D0*Hr4(-1,0,1,1)*x - 80.D0*
     &    Hr4(0,-1,-1,0) + 96.D0*Hr4(0,-1,-1,0)*x - 104.D0*Hr4(0,-1,0,0
     &    ) + 48.D0*Hr4(0,-1,0,0)*x - 32.D0*Hr4(0,-1,0,1) - 64.D0*Hr4(0
     &    ,-1,0,1)*x - 8.D0*Hr4(0,0,-1,0) + 144.D0*Hr4(0,0,-1,0)*x - 
     &    120.D0*Hr4(0,0,0,0) + 128.D0*Hr4(0,0,0,0)*x - 36.D0*Hr4(0,0,0
     &    ,1) + 24.D0*Hr4(0,0,0,1)*x + 64.D0*Hr4(0,0,1,1) + 96.D0*Hr4(0
     &    ,0,1,1)*x - 100.D0*Hr4(0,1,0,0) - 168.D0*Hr4(0,1,0,0)*x + 64.D
     &    0*Hr4(0,1,1,0) + 64.D0*Hr4(0,1,1,0)*x + 32.D0*Hr4(0,1,1,1) + 
     &    32.D0*Hr4(0,1,1,1)*x + 48.D0*Hr4(1,0,-1,0) - 96.D0*Hr4(1,0,-1
     &    ,0)*x - 40.D0*Hr4(1,0,0,0) + 80.D0*Hr4(1,0,0,0)*x - 48.D0*
     &    Hr4(1,0,0,1) + 96.D0*Hr4(1,0,0,1)*x + 32.D0*Hr4(1,0,1,0) - 64.
     &    D0*Hr4(1,0,1,0)*x + 16.D0*Hr4(1,1,0,0) - 32.D0*Hr4(1,1,0,0)*x
     &     + 32.D0*Hr4(1,1,0,1) - 64.D0*Hr4(1,1,0,1)*x + 64.D0*Hr4(1,1,
     &    1,0) - 128.D0*Hr4(1,1,1,0)*x - 32.D0*Hr4(1,1,1,1) + 64.D0*
     &    Hr4(1,1,1,1)*x )
      dPqg2M = dPqg2M + nf*cf*ca * ( 69407.D0/36.D0 - 34495.D0/18.D0*
     &    x - 202.D0*z4 - 292.D0*z4*x - 2960.D0/3.D0*z3 + 2344.D0/3.D0*
     &    z3*x - 120.D0*z3*x**2 - 808.D0/9.D0*z2 - 340.D0/9.D0*z2*x - 
     &    72.D0*z2*x**2 + 136.D0*Hr1(-1)*z3 + 272.D0*Hr1(-1)*z3*x + 184.
     &    D0*Hr1(-1)*z2 + 256.D0*Hr1(-1)*z2*x + 72.D0*Hr1(-1)*z2*x**2
     &     + 15259.D0/27.D0*Hr1(0) + 25645.D0/27.D0*Hr1(0)*x - 328.D0*
     &    Hr1(0)*z3 - 176.D0*Hr1(0)*z3*x + 72.D0*Hr1(0)*z2 - 280.D0*
     &    Hr1(0)*z2*x - 96.D0*Hr1(0)*z2*x**2 + 9560.D0/27.D0*Hr1(1) - 
     &    11464.D0/27.D0*Hr1(1)*x - 312.D0*Hr1(1)*z3 + 624.D0*Hr1(1)*z3
     &    *x - 236.D0*Hr1(1)*z2 + 248.D0*Hr1(1)*z2*x - 24.D0*Hr1(1)*z2*
     &    x**2 - 144.D0*Hr2(-1,-1)*z2 - 288.D0*Hr2(-1,-1)*z2*x - 64.D0*
     &    Hr2(-1,0) - 136.D0*Hr2(-1,0)*x - 72.D0*Hr2(-1,0)*x**2 + 144.D0
     &    *Hr2(-1,0)*z2 + 288.D0*Hr2(-1,0)*z2*x + 64.D0*Hr2(0,-1)*z2 + 
     &    224.D0*Hr2(0,-1)*z2*x - 701.D0/9.D0*Hr2(0,0) - 2714.D0/9.D0*
     &    Hr2(0,0)*x + 72.D0*Hr2(0,0)*x**2 + 48.D0*Hr2(0,0)*z2 + 192.D0
     &    *Hr2(0,0)*z2*x + 232.D0/9.D0*Hr2(0,1) + 340.D0/9.D0*Hr2(0,1)*
     &    x - 32.D0*Hr2(0,1)*z2 + 32.D0*Hr2(0,1)*z2*x - 136.D0/9.D0*
     &    Hr2(1,0) + 368.D0/9.D0*Hr2(1,0)*x - 144.D0*Hr2(1,0)*z2 + 288.D
     &    0*Hr2(1,0)*z2*x + 2536.D0/9.D0*Hr2(1,1) - 2504.D0/9.D0*Hr2(1,
     &    1)*x - 96.D0*Hr2(1,1)*z2 + 192.D0*Hr2(1,1)*z2*x + 208.D0*Hr3(
     &    -1,-1,0) + 256.D0*Hr3(-1,-1,0)*x + 48.D0*Hr3(-1,-1,0)*x**2 - 
     &    120.D0*Hr3(-1,0,0) - 144.D0*Hr3(-1,0,0)*x - 48.D0*Hr3(-1,0,0)
     &    *x**2 - 80.D0*Hr3(-1,0,1) - 128.D0*Hr3(-1,0,1)*x - 48.D0*Hr3(
     &    -1,0,1)*x**2 - 16.D0*Hr3(0,-1,0) - 256.D0*Hr3(0,-1,0)*x - 48.D
     &    0*Hr3(0,-1,0)*x**2 - 356.D0/3.D0*Hr3(0,0,0) + 520.D0/3.D0*
     &    Hr3(0,0,0)*x + 48.D0*Hr3(0,0,0)*x**2 - 88.D0*Hr3(0,0,1) + 280.
     &    D0*Hr3(0,0,1)*x + 48.D0*Hr3(0,0,1)*x**2 - 80.D0*Hr3(0,1,0) + 
     &    280.D0*Hr3(0,1,0)*x + 164.D0/3.D0*Hr3(0,1,1) + 368.D0/3.D0*
     &    Hr3(0,1,1)*x + 488.D0/3.D0*Hr3(1,0,0) - 472.D0/3.D0*Hr3(1,0,0
     &    )*x + 132.D0*Hr3(1,0,1) - 120.D0*Hr3(1,0,1)*x - 132.D0*Hr3(1,
     &    1,0) + 120.D0*Hr3(1,1,0)*x + 296.D0/3.D0*Hr3(1,1,1) - 376.D0/
     &    3.D0*Hr3(1,1,1)*x - 96.D0*Hr4(-1,-1,-1,0) - 192.D0*Hr4(-1,-1,
     &    -1,0)*x + 96.D0*Hr4(-1,-1,0,0) + 192.D0*Hr4(-1,-1,0,0)*x + 96.
     &    D0*Hr4(-1,-1,0,1) + 192.D0*Hr4(-1,-1,0,1)*x + 48.D0*Hr4(-1,0,
     &    -1,0) + 96.D0*Hr4(-1,0,-1,0)*x - 88.D0*Hr4(-1,0,0,0) - 176.D0
     &    *Hr4(-1,0,0,0)*x - 96.D0*Hr4(-1,0,0,1) - 192.D0*Hr4(-1,0,0,1)
     &    *x - 32.D0*Hr4(-1,0,1,0) - 64.D0*Hr4(-1,0,1,0)*x - 32.D0*Hr4(
     &    -1,0,1,1) - 64.D0*Hr4(-1,0,1,1)*x + 64.D0*Hr4(0,-1,-1,0) + 64.
     &    D0*Hr4(0,-1,-1,0)*x - 48.D0*Hr4(0,-1,0,0) - 128.D0*Hr4(0,-1,0
     &    ,0)*x - 32.D0*Hr4(0,-1,0,1) - 192.D0*Hr4(0,-1,0,1)*x - 16.D0*
     &    Hr4(0,0,-1,0) + 96.D0*Hr4(0,0,-1,0)*x - 32.D0*Hr4(0,0,0,0) - 
     &    64.D0*Hr4(0,0,0,0)*x - 64.D0*Hr4(0,0,0,1) - 192.D0*Hr4(0,0,0,
     &    1)*x - 96.D0*Hr4(0,0,1,0) - 128.D0*Hr4(0,0,1,0)*x - 48.D0*
     &    Hr4(0,0,1,1) - 128.D0*Hr4(0,0,1,1)*x + 8.D0*Hr4(0,1,0,0) + 
     &    176.D0*Hr4(0,1,0,0)*x - 96.D0*Hr4(0,1,1,0) + 8.D0*Hr4(0,1,1,1
     &    ) - 112.D0*Hr4(0,1,1,1)*x - 48.D0*Hr4(1,0,-1,0) + 96.D0*Hr4(1
     &    ,0,-1,0)*x + 24.D0*Hr4(1,0,0,0) - 48.D0*Hr4(1,0,0,0)*x + 96.D0
     &    *Hr4(1,0,0,1) - 192.D0*Hr4(1,0,0,1)*x + 96.D0*Hr4(1,0,1,1) - 
     &    192.D0*Hr4(1,0,1,1)*x - 16.D0*Hr4(1,1,0,0) + 32.D0*Hr4(1,1,0,
     &    0)*x + 48.D0*Hr4(1,1,0,1) - 96.D0*Hr4(1,1,0,1)*x - 48.D0*Hr4(
     &    1,1,1,0) + 96.D0*Hr4(1,1,1,0)*x + 64.D0*Hr4(1,1,1,1) - 128.D0
     &    *Hr4(1,1,1,1)*x )
      dPqg2M = dPqg2M + nf*cf**2 * (  - 1357.D0/4.D0 + 681.D0/2.D0*x
     &     - 48.D0*z4 + 160.D0*z4*x + 200.D0*z3 - 88.D0*z3*x + 236.D0*
     &    z2 - 144.D0*z2*x + 48.D0*z2*x**2 - 32.D0*Hr1(-1)*z2 - 32.D0*
     &    Hr1(-1)*z2*x - 77.D0*Hr1(0) + 397.D0*Hr1(0)*x + 104.D0*Hr1(0)
     &    *z3 - 112.D0*Hr1(0)*z3*x + 92.D0*Hr1(0)*z2 - 16.D0*Hr1(0)*z2*
     &    x - 624.D0*Hr1(1) + 658.D0*Hr1(1)*x + 112.D0*Hr1(1)*z3 - 224.D
     &    0*Hr1(1)*z3*x + 208.D0*Hr1(1)*z2 - 208.D0*Hr1(1)*z2*x - 160.D0
     &    *Hr2(-1,0) - 144.D0*Hr2(-1,0)*x + 48.D0*Hr2(-1,0)*x**2 + 64.D0
     &    *Hr2(0,-1)*z2*x + 59.D0*Hr2(0,0) + 458.D0*Hr2(0,0)*x - 48.D0*
     &    Hr2(0,0)*x**2 + 56.D0*Hr2(0,0)*z2 - 112.D0*Hr2(0,0)*z2*x - 
     &    396.D0*Hr2(0,1) + 144.D0*Hr2(0,1)*x + 64.D0*Hr2(0,1)*z2 - 64.D
     &    0*Hr2(0,1)*z2*x - 272.D0*Hr2(1,0) + 352.D0*Hr2(1,0)*x + 112.D0
     &    *Hr2(1,0)*z2 - 224.D0*Hr2(1,0)*z2*x - 232.D0*Hr2(1,1) + 288.D0
     &    *Hr2(1,1)*x + 80.D0*Hr2(1,1)*z2 - 160.D0*Hr2(1,1)*z2*x - 64.D0
     &    *Hr3(-1,-1,0) - 64.D0*Hr3(-1,-1,0)*x - 64.D0*Hr3(-1,0,0) - 64.
     &    D0*Hr3(-1,0,0)*x - 80.D0*Hr3(0,-1,0) + 160.D0*Hr3(0,-1,0)*x
     &     + 60.D0*Hr3(0,0,0) - 16.D0*Hr3(0,0,0)*x - 172.D0*Hr3(0,0,1)
     &     + 16.D0*Hr3(0,0,1)*x - 160.D0*Hr3(0,1,0) + 16.D0*Hr3(0,1,0)*
     &    x - 140.D0*Hr3(0,1,1) + 56.D0*Hr3(0,1,1)*x - 176.D0*Hr3(1,0,0
     &    ) + 176.D0*Hr3(1,0,0)*x - 176.D0*Hr3(1,0,1) + 176.D0*Hr3(1,0,
     &    1)*x - 80.D0*Hr3(1,1,0) + 80.D0*Hr3(1,1,0)*x - 84.D0*Hr3(1,1,
     &    1) + 96.D0*Hr3(1,1,1)*x + 64.D0*Hr4(-1,-1,0,0) + 128.D0*Hr4(
     &    -1,-1,0,0)*x + 32.D0*Hr4(-1,0,-1,0) + 64.D0*Hr4(-1,0,-1,0)*x
     &     - 16.D0*Hr4(-1,0,0,0) - 32.D0*Hr4(-1,0,0,0)*x + 128.D0*Hr4(0
     &    ,-1,-1,0)*x - 32.D0*Hr4(0,-1,0,0) - 128.D0*Hr4(0,-1,0,0)*x - 
     &    128.D0*Hr4(0,0,-1,0)*x + 32.D0*Hr4(0,0,0,0)*x - 56.D0*Hr4(0,0
     &    ,0,1) + 112.D0*Hr4(0,0,0,1)*x - 64.D0*Hr4(0,0,1,0) + 128.D0*
     &    Hr4(0,0,1,0)*x - 80.D0*Hr4(0,0,1,1) + 160.D0*Hr4(0,0,1,1)*x
     &     - 48.D0*Hr4(0,1,0,0) + 32.D0*Hr4(0,1,0,0)*x - 64.D0*Hr4(0,1,
     &    0,1) + 128.D0*Hr4(0,1,0,1)*x - 32.D0*Hr4(0,1,1,0) + 64.D0*
     &    Hr4(0,1,1,0)*x - 40.D0*Hr4(0,1,1,1) + 80.D0*Hr4(0,1,1,1)*x + 
     &    32.D0*Hr4(1,0,-1,0) - 64.D0*Hr4(1,0,-1,0)*x - 48.D0*Hr4(1,0,0
     &    ,0) + 96.D0*Hr4(1,0,0,0)*x - 112.D0*Hr4(1,0,0,1) + 224.D0*
     &    Hr4(1,0,0,1)*x - 96.D0*Hr4(1,0,1,0) + 192.D0*Hr4(1,0,1,0)*x
     &     - 96.D0*Hr4(1,0,1,1) + 192.D0*Hr4(1,0,1,1)*x - 64.D0*Hr4(1,1
     &    ,0,0) + 128.D0*Hr4(1,1,0,0)*x - 80.D0*Hr4(1,1,0,1) + 160.D0*
     &    Hr4(1,1,0,1)*x - 16.D0*Hr4(1,1,1,0) + 32.D0*Hr4(1,1,1,0)*x - 
     &    32.D0*Hr4(1,1,1,1) + 64.D0*Hr4(1,1,1,1)*x )
      dPqg2M = dPqg2M + nf2*ca * ( 944.D0/9.D0 - 1054.D0/9.D0*x - 8.
     &    D0*z3 + 32.D0/9.D0*z2 + 8.D0/3.D0*z2*x + 2134.D0/27.D0*Hr1(0)
     &     + 400.D0/27.D0*Hr1(0)*x - 184.D0/27.D0*Hr1(1) + 32.D0/27.D0*
     &    Hr1(1)*x + 8.D0/3.D0*Hr1(1)*z2 - 16.D0/3.D0*Hr1(1)*z2*x + 32.D
     &    0/9.D0*Hr2(-1,0) + 112.D0/9.D0*Hr2(-1,0)*x + 44.D0*Hr2(0,0)
     &     + 272.D0/9.D0*Hr2(0,0)*x - 8.D0/3.D0*Hr2(0,1)*x - 32.D0/9.D0
     &    *Hr2(1,1) + 112.D0/9.D0*Hr2(1,1)*x + 16.D0/3.D0*Hr3(-1,0,0)
     &     + 32.D0/3.D0*Hr3(-1,0,0)*x + 32.D0/3.D0*Hr3(0,0,0) - 32.D0/3.
     &    D0*Hr3(0,0,0)*x + 8.D0/3.D0*Hr3(1,0,0) - 16.D0/3.D0*Hr3(1,0,0
     &    )*x - 8.D0/3.D0*Hr3(1,0,1) + 16.D0/3.D0*Hr3(1,0,1)*x + 8.D0/3.
     &    D0*Hr3(1,1,0) - 16.D0/3.D0*Hr3(1,1,0)*x + 8.D0/3.D0*Hr3(1,1,1
     &    ) - 16.D0/3.D0*Hr3(1,1,1)*x )
      dPqg2M = dPqg2M + nf2*cf * ( 4193.D0/18.D0 - 2113.D0/9.D0*x
     &     + 8.D0/3.D0*z3 - 16.D0/3.D0*z3*x - 32.D0/9.D0*z2 + 112.D0/9.D
     &    0*z2*x + 6434.D0/27.D0*Hr1(0) - 118.D0/27.D0*Hr1(0)*x + 112.D0
     &    /27.D0*Hr1(1) + 76.D0/27.D0*Hr1(1)*x + 1802.D0/9.D0*Hr2(0,0)
     &     + 476.D0/9.D0*Hr2(0,0)*x + 32.D0/9.D0*Hr2(0,1) - 112.D0/9.D0
     &    *Hr2(0,1)*x + 88.D0/9.D0*Hr2(1,0) - 128.D0/9.D0*Hr2(1,0)*x + 
     &    32.D0/9.D0*Hr2(1,1) - 112.D0/9.D0*Hr2(1,1)*x + 344.D0/3.D0*
     &    Hr3(0,0,0) + 80.D0/3.D0*Hr3(0,0,0)*x - 8.D0/3.D0*Hr3(0,1,1)
     &     + 16.D0/3.D0*Hr3(0,1,1)*x + 16.D0/3.D0*Hr3(1,0,0) - 32.D0/3.D
     &    0*Hr3(1,0,0)*x - 8.D0/3.D0*Hr3(1,1,1) + 16.D0/3.D0*Hr3(1,1,1)
     &    *x + 32.D0*Hr4(0,0,0,0) - 64.D0*Hr4(0,0,0,0)*x )
*
       DP2QGMA = dPqg2M
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 3-loop quark->gluon splitting function DeltaP_gq^(2)
*
       FUNCTION DP2GQMA (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,
     ,             Z4 = 1.0823 23233 71113 81915 D0 )
*
* ...Colour factors and abbreviation
*
       CF = 4./3.D0
       CA = 3.D0
       NF2 = NF*NF
*
       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
*
      dPgq2M =
     &  + cf*ca**2 * ( 7436.D0/3.D0 - 6698.D0/3.D0*x - 366.D0*z4 - 193.D
     &    0*z4*x - 1288.D0/3.D0*z3 + 1552.D0/3.D0*z3*x - 3328.D0/9.D0*
     &    z2 - 1168.D0/9.D0*z2*x - 368.D0*Hr1(-1)*z3 - 184.D0*Hr1(-1)*
     &    z3*x + 104.D0/3.D0*Hr1(-1)*z2 - 32.D0/3.D0*Hr1(-1)*z2*x - 72.D
     &    0*Hr1(-1)*z2*dx + 23996.D0/27.D0*Hr1(0) + 37298.D0/27.D0*Hr1(
     &    0)*x - 160.D0*Hr1(0)*z3 - 112.D0*Hr1(0)*z3*x - 188.D0/3.D0*
     &    Hr1(0)*z2 + 1676.D0/3.D0*Hr1(0)*z2*x - 5768.D0/27.D0*Hr1(1)
     &     + 1612.D0/27.D0*Hr1(1)*x + 112.D0*Hr1(1)*z3 - 56.D0*Hr1(1)*
     &    z3*x - 368.D0/3.D0*Hr1(1)*z2 + 484.D0/3.D0*Hr1(1)*z2*x - 24.D0
     &    *Hr1(1)*z2*dx + 416.D0*Hr2(-1,-1)*z2 + 208.D0*Hr2(-1,-1)*z2*x
     &     - 848.D0/3.D0*Hr2(-1,0) - 1384.D0/3.D0*Hr2(-1,0)*x - 24.D0*
     &    Hr2(-1,0)*dx - 320.D0*Hr2(-1,0)*z2 - 160.D0*Hr2(-1,0)*z2*x - 
     &    112.D0*Hr2(0,-1)*z2 - 216.D0*Hr2(0,-1)*z2*x + 5776.D0/9.D0*
     &    Hr2(0,0) - 9404.D0/9.D0*Hr2(0,0)*x - 184.D0*Hr2(0,0)*z2 - 52.D
     &    0*Hr2(0,0)*z2*x + 784.D0/9.D0*Hr2(0,1) + 1168.D0/9.D0*Hr2(0,1
     &    )*x - 112.D0*Hr2(0,1)*z2 - 40.D0*Hr2(0,1)*z2*x - 40.D0/9.D0*
     &    Hr2(1,0) + 1184.D0/9.D0*Hr2(1,0)*x - 64.D0*Hr2(1,0)*z2 + 32.D0
     &    *Hr2(1,0)*z2*x + 100.D0/9.D0*Hr2(1,1) + 1162.D0/9.D0*Hr2(1,1)
     &    *x - 32.D0*Hr2(1,1)*z2 + 16.D0*Hr2(1,1)*z2*x + 848.D0/3.D0*
     &    Hr3(-1,-1,0) + 640.D0/3.D0*Hr3(-1,-1,0)*x - 48.D0*Hr3(-1,-1,0
     &    )*dx + 1064.D0/3.D0*Hr3(-1,0,0) + 1096.D0/3.D0*Hr3(-1,0,0)*x
     &     + 48.D0*Hr3(-1,0,0)*dx + 320.D0/3.D0*Hr3(-1,0,1) + 352.D0/3.D
     &    0*Hr3(-1,0,1)*x + 48.D0*Hr3(-1,0,1)*dx - 24.D0*Hr3(0,-1,0) + 
     &    976.D0/3.D0*Hr3(0,-1,0)*x + 376.D0*Hr3(0,0,0) - 40.D0/3.D0*
     &    Hr3(0,0,0)*x + 116.D0/3.D0*Hr3(0,0,1) - 1676.D0/3.D0*Hr3(0,0,
     &    1)*x + 304.D0/3.D0*Hr3(0,1,0) - 1112.D0/3.D0*Hr3(0,1,0)*x - 
     &    56.D0*Hr3(0,1,1) - 140.D0*Hr3(0,1,1)*x + 980.D0/3.D0*Hr3(1,0,
     &    0) - 892.D0/3.D0*Hr3(1,0,0)*x - 56.D0/3.D0*Hr3(1,0,1) - 164.D0
     &    /3.D0*Hr3(1,0,1)*x - 56.D0/3.D0*Hr3(1,1,0) - 164.D0/3.D0*Hr3(
     &    1,1,0)*x - 440.D0/3.D0*Hr3(1,1,1) + 220.D0/3.D0*Hr3(1,1,1)*x
     &     + 192.D0*Hr4(-1,-1,-1,0) + 96.D0*Hr4(-1,-1,-1,0)*x - 384.D0*
     &    Hr4(-1,-1,0,0) - 192.D0*Hr4(-1,-1,0,0)*x - 320.D0*Hr4(-1,-1,0
     &    ,1) - 160.D0*Hr4(-1,-1,0,1)*x - 160.D0*Hr4(-1,0,-1,0) - 80.D0
     &    *Hr4(-1,0,-1,0)*x + 176.D0*Hr4(-1,0,0,0) + 88.D0*Hr4(-1,0,0,0
     &    )*x + 256.D0*Hr4(-1,0,0,1) + 128.D0*Hr4(-1,0,0,1)*x + 64.D0*
     &    Hr4(-1,0,1,0) + 32.D0*Hr4(-1,0,1,0)*x + 64.D0*Hr4(-1,0,1,1)
     &     + 32.D0*Hr4(-1,0,1,1)*x + 32.D0*Hr4(0,-1,-1,0) - 176.D0*Hr4(
     &    0,-1,-1,0)*x + 304.D0*Hr4(0,-1,0,0) + 88.D0*Hr4(0,-1,0,0)*x
     &     + 128.D0*Hr4(0,-1,0,1) + 128.D0*Hr4(0,-1,0,1)*x + 16.D0*Hr4(
     &    0,0,-1,0) - 72.D0*Hr4(0,0,-1,0)*x + 240.D0*Hr4(0,0,0,0) - 64.D
     &    0*Hr4(0,0,0,0)*x + 200.D0*Hr4(0,0,0,1) + 52.D0*Hr4(0,0,0,1)*x
     &     + 256.D0*Hr4(0,0,1,0) + 32.D0*Hr4(0,0,1,0)*x + 192.D0*Hr4(0,
     &    0,1,1) - 32.D0*Hr4(0,0,1,1)*x + 168.D0*Hr4(0,1,0,0) + 52.D0*
     &    Hr4(0,1,0,0)*x + 96.D0*Hr4(0,1,0,1) - 48.D0*Hr4(0,1,0,1)*x + 
     &    96.D0*Hr4(0,1,1,0) - 48.D0*Hr4(0,1,1,0)*x + 32.D0*Hr4(0,1,1,1
     &    ) - 64.D0*Hr4(0,1,1,1)*x - 32.D0*Hr4(1,0,-1,0) + 16.D0*Hr4(1,
     &    0,-1,0)*x + 48.D0*Hr4(1,0,0,0) - 24.D0*Hr4(1,0,0,0)*x + 128.D0
     &    *Hr4(1,0,0,1) - 64.D0*Hr4(1,0,0,1)*x + 192.D0*Hr4(1,0,1,0) - 
     &    96.D0*Hr4(1,0,1,0)*x + 128.D0*Hr4(1,0,1,1) - 64.D0*Hr4(1,0,1,
     &    1)*x + 64.D0*Hr4(1,1,0,0) - 32.D0*Hr4(1,1,0,0)*x + 128.D0*
     &    Hr4(1,1,0,1) - 64.D0*Hr4(1,1,0,1)*x + 128.D0*Hr4(1,1,1,0) - 
     &    64.D0*Hr4(1,1,1,0)*x + 64.D0*Hr4(1,1,1,1) - 32.D0*Hr4(1,1,1,1
     &    )*x )
      dPgq2M = dPgq2M + cf**2*ca * (  - 2911.D0/6.D0 + 11273.D0/36.D0
     &    *x - 184.D0*z4 - 196.D0*z4*x + 896.D0/3.D0*z3 - 700.D0/3.D0*
     &    z3*x - 944.D0/9.D0*z2 + 1168.D0/9.D0*z2*x + 176.D0*Hr1(-1)*z3
     &     + 88.D0*Hr1(-1)*z3*x - 176.D0*Hr1(-1)*z2 - 80.D0*Hr1(-1)*z2*
     &    x - 13037.D0/27.D0*Hr1(0) + 4423.D0/27.D0*Hr1(0)*x + 128.D0*
     &    Hr1(0)*z3 - 32.D0*Hr1(0)*z3*x - 208.D0/3.D0*Hr1(0)*z2 - 64.D0/
     &    3.D0*Hr1(0)*z2*x + 1172.D0/27.D0*Hr1(1) + 1226.D0/27.D0*Hr1(1
     &    )*x + 48.D0*Hr1(1)*z3 - 24.D0*Hr1(1)*z3*x + 296.D0/3.D0*Hr1(1
     &    )*z2 - 244.D0/3.D0*Hr1(1)*z2*x - 160.D0*Hr2(-1,-1)*z2 - 80.D0
     &    *Hr2(-1,-1)*z2*x - 48.D0*Hr2(-1,0) - 8.D0*Hr2(-1,0)*x + 72.D0
     &    *Hr2(-1,0)*dx - 32.D0*Hr2(-1,0)*z2 - 16.D0*Hr2(-1,0)*z2*x + 
     &    48.D0*Hr2(0,-1)*z2*x + 730.D0/9.D0*Hr2(0,0) - 821.D0/9.D0*
     &    Hr2(0,0)*x - 32.D0*Hr2(0,0)*z2 + 512.D0/9.D0*Hr2(0,1) - 1168.D
     &    0/9.D0*Hr2(0,1)*x - 64.D0*Hr2(0,1)*z2 + 80.D0*Hr2(0,1)*z2*x
     &     + 320.D0/3.D0*Hr2(1,0) - 40.D0/3.D0*Hr2(1,0)*x + 32.D0*Hr2(1
     &    ,0)*z2 - 16.D0*Hr2(1,0)*z2*x - 124.D0/9.D0*Hr2(1,1) - 1096.D0/
     &    9.D0*Hr2(1,1)*x - 64.D0*Hr2(1,1)*z2 + 32.D0*Hr2(1,1)*z2*x - 
     &    224.D0*Hr3(-1,-1,0) - 128.D0*Hr3(-1,-1,0)*x + 96.D0*Hr3(-1,0,
     &    0) + 48.D0*Hr3(-1,0,0)*x + 64.D0*Hr3(-1,0,1) + 16.D0*Hr3(-1,0
     &    ,1)*x - 48.D0*Hr3(0,-1,0) + 72.D0*Hr3(0,-1,0)*x + 712.D0/3.D0
     &    *Hr3(0,0,0) - 704.D0/3.D0*Hr3(0,0,0)*x + 64.D0/3.D0*Hr3(0,0,1
     &    ) + 64.D0/3.D0*Hr3(0,0,1)*x - 80.D0/3.D0*Hr3(0,1,0) + 16.D0/3.
     &    D0*Hr3(0,1,0)*x - 64.D0/3.D0*Hr3(0,1,1) + 344.D0/3.D0*Hr3(0,1
     &    ,1)*x - 80.D0*Hr3(1,0,0) + 24.D0*Hr3(1,0,0)*x + 40.D0/3.D0*
     &    Hr3(1,0,1) + 52.D0/3.D0*Hr3(1,0,1)*x - 280.D0/3.D0*Hr3(1,1,0)
     &     + 284.D0/3.D0*Hr3(1,1,0)*x + 584.D0/3.D0*Hr3(1,1,1) - 184.D0/
     &    3.D0*Hr3(1,1,1)*x - 192.D0*Hr4(-1,-1,-1,0) - 96.D0*Hr4(-1,-1,
     &    -1,0)*x + 128.D0*Hr4(-1,-1,0,0) + 64.D0*Hr4(-1,-1,0,0)*x + 64.
     &    D0*Hr4(-1,-1,0,1) + 32.D0*Hr4(-1,-1,0,1)*x + 96.D0*Hr4(-1,0,
     &    -1,0) + 48.D0*Hr4(-1,0,-1,0)*x + 48.D0*Hr4(-1,0,0,0) + 24.D0*
     &    Hr4(-1,0,0,0)*x - 64.D0*Hr4(-1,0,1,1) - 32.D0*Hr4(-1,0,1,1)*x
     &     + 96.D0*Hr4(0,-1,-1,0)*x - 48.D0*Hr4(0,-1,0,0)*x + 32.D0*
     &    Hr4(0,0,-1,0) - 112.D0*Hr4(0,0,-1,0)*x + 128.D0*Hr4(0,0,0,0)
     &     + 32.D0*Hr4(0,0,0,0)*x + 64.D0*Hr4(0,0,0,1) + 64.D0*Hr4(0,0,
     &    1,0) - 32.D0*Hr4(0,0,1,0)*x - 16.D0*Hr4(0,0,1,1) - 56.D0*Hr4(
     &    0,0,1,1)*x + 128.D0*Hr4(0,1,0,0) - 112.D0*Hr4(0,1,0,0)*x + 64.
     &    D0*Hr4(0,1,0,1) - 32.D0*Hr4(0,1,0,1)*x + 32.D0*Hr4(0,1,1,0)
     &     - 16.D0*Hr4(0,1,1,0)*x - 16.D0*Hr4(0,1,1,1) + 56.D0*Hr4(0,1,
     &    1,1)*x + 160.D0*Hr4(1,0,-1,0) - 80.D0*Hr4(1,0,-1,0)*x + 80.D0
     &    *Hr4(1,0,0,0) - 40.D0*Hr4(1,0,0,0)*x - 64.D0*Hr4(1,0,1,1) + 
     &    32.D0*Hr4(1,0,1,1)*x + 160.D0*Hr4(1,1,0,0) - 80.D0*Hr4(1,1,0,
     &    0)*x - 32.D0*Hr4(1,1,0,1) + 16.D0*Hr4(1,1,0,1)*x - 96.D0*Hr4(
     &    1,1,1,0) + 48.D0*Hr4(1,1,1,0)*x - 128.D0*Hr4(1,1,1,1) + 64.D0
     &    *Hr4(1,1,1,1)*x )
      dPgq2M = dPgq2M + cf**3 * ( 2331.D0/2.D0 - 4547.D0/4.D0*x + 60.D
     &    0*z4 + 186.D0*z4*x - 416.D0*z3 + 168.D0*z3*x - 100.D0*z2 - 92.
     &    D0*z2*x + 64.D0*Hr1(-1)*z2 + 64.D0*Hr1(-1)*z2*x + 389.D0*Hr1(
     &    0) + 721.D0*Hr1(0)*x - 96.D0*Hr1(0)*z3 + 96.D0*Hr1(0)*z3*x - 
     &    4.D0*Hr1(0)*z2*x + 64.D0*Hr1(1) - 46.D0*Hr1(1)*x - 256.D0*
     &    Hr1(1)*z3 + 128.D0*Hr1(1)*z3*x - 32.D0*Hr1(1)*z2 + 8.D0*Hr1(1
     &    )*z2*x - 208.D0*Hr2(-1,0) - 192.D0*Hr2(-1,0)*x - 48.D0*Hr2(-1
     &    ,0)*dx + 64.D0*Hr2(-1,0)*z2 + 32.D0*Hr2(-1,0)*z2*x - 32.D0*
     &    Hr2(0,-1)*z2*x - 78.D0*Hr2(0,0) - 93.D0*Hr2(0,0)*x - 48.D0*
     &    Hr2(0,0)*z2 + 24.D0*Hr2(0,0)*z2*x - 108.D0*Hr2(0,1) + 92.D0*
     &    Hr2(0,1)*x + 32.D0*Hr2(0,1)*z2 - 48.D0*Hr2(0,1)*z2*x - 160.D0
     &    *Hr2(1,0) + 64.D0*Hr2(1,0)*x - 64.D0*Hr2(1,0)*z2 + 32.D0*Hr2(
     &    1,0)*z2*x - 16.D0*Hr2(1,1) + 14.D0*Hr2(1,1)*x + 96.D0*Hr2(1,1
     &    )*z2 - 48.D0*Hr2(1,1)*z2*x + 128.D0*Hr3(-1,-1,0) + 128.D0*
     &    Hr3(-1,-1,0)*x + 32.D0*Hr3(-1,0,0) + 32.D0*Hr3(-1,0,0)*x - 32.
     &    D0*Hr3(0,-1,0) + 16.D0*Hr3(0,-1,0)*x - 184.D0*Hr3(0,0,0) - 24.
     &    D0*Hr3(0,0,0)*x - 32.D0*Hr3(0,0,1) + 4.D0*Hr3(0,0,1)*x + 16.D0
     &    *Hr3(0,1,0)*x + 32.D0*Hr3(0,1,1) + 28.D0*Hr3(0,1,1)*x + 48.D0
     &    *Hr3(1,0,0) - 12.D0*Hr3(1,0,0)*x - 32.D0*Hr3(1,0,1) + 56.D0*
     &    Hr3(1,0,1)*x + 64.D0*Hr3(1,1,0) - 16.D0*Hr3(1,1,0)*x - 48.D0*
     &    Hr3(1,1,1) - 12.D0*Hr3(1,1,1)*x - 128.D0*Hr4(-1,-1,0,0) - 64.D
     &    0*Hr4(-1,-1,0,0)*x - 64.D0*Hr4(-1,0,-1,0) - 32.D0*Hr4(-1,0,-1
     &    ,0)*x + 32.D0*Hr4(-1,0,0,0) + 16.D0*Hr4(-1,0,0,0)*x - 64.D0*
     &    Hr4(0,-1,-1,0)*x + 64.D0*Hr4(0,-1,0,0) + 64.D0*Hr4(0,-1,0,0)*
     &    x + 128.D0*Hr4(0,0,-1,0)*x - 64.D0*Hr4(0,0,0,0) - 16.D0*Hr4(0
     &    ,0,0,0)*x + 48.D0*Hr4(0,0,0,1) - 24.D0*Hr4(0,0,0,1)*x - 48.D0
     &    *Hr4(0,0,1,1) + 24.D0*Hr4(0,0,1,1)*x - 16.D0*Hr4(0,1,0,0) + 
     &    40.D0*Hr4(0,1,0,0)*x - 32.D0*Hr4(0,1,0,1) + 16.D0*Hr4(0,1,0,1
     &    )*x - 16.D0*Hr4(0,1,1,1) + 8.D0*Hr4(0,1,1,1)*x - 192.D0*Hr4(1
     &    ,0,-1,0) + 96.D0*Hr4(1,0,-1,0)*x - 64.D0*Hr4(1,0,1,0) + 32.D0
     &    *Hr4(1,0,1,0)*x - 64.D0*Hr4(1,0,1,1) + 32.D0*Hr4(1,0,1,1)*x
     &     - 96.D0*Hr4(1,1,0,0) + 48.D0*Hr4(1,1,0,0)*x - 96.D0*Hr4(1,1,
     &    0,1) + 48.D0*Hr4(1,1,0,1)*x - 32.D0*Hr4(1,1,1,0) + 16.D0*Hr4(
     &    1,1,1,0)*x + 64.D0*Hr4(1,1,1,1) - 32.D0*Hr4(1,1,1,1)*x )
      dPgq2M = dPgq2M + nf*cf*ca * ( 246.D0 - 2572.D0/9.D0*x - 368.D0/
     &    3.D0*z3 + 224.D0/3.D0*z3*x - 464.D0/9.D0*z2 + 364.D0/9.D0*z2*
     &    x - 128.D0/3.D0*Hr1(-1)*z2 - 64.D0/3.D0*Hr1(-1)*z2*x + 1552.D0
     &    /27.D0*Hr1(0) + 6286.D0/27.D0*Hr1(0)*x - 160.D0/3.D0*Hr1(0)*
     &    z2 - 32.D0/3.D0*Hr1(0)*z2*x + 728.D0/27.D0*Hr1(1) + 944.D0/27.
     &    D0*Hr1(1)*x - 16.D0/3.D0*Hr1(1)*z2 + 8.D0/3.D0*Hr1(1)*z2*x - 
     &    160.D0/3.D0*Hr2(-1,0) - 128.D0/3.D0*Hr2(-1,0)*x + 16.D0*Hr2(
     &    -1,0)*dx - 328.D0/9.D0*Hr2(0,0) - 592.D0/9.D0*Hr2(0,0)*x - 16.
     &    D0/9.D0*Hr2(0,1) - 364.D0/9.D0*Hr2(0,1)*x - 80.D0/9.D0*Hr2(1,
     &    0) - 80.D0/9.D0*Hr2(1,0)*x - 640.D0/9.D0*Hr2(1,1) + 296.D0/9.D
     &    0*Hr2(1,1)*x - 128.D0/3.D0*Hr3(-1,-1,0) - 64.D0/3.D0*Hr3(-1,
     &    -1,0)*x + 64.D0/3.D0*Hr3(-1,0,0) + 32.D0/3.D0*Hr3(-1,0,0)*x
     &     + 64.D0/3.D0*Hr3(-1,0,1) + 32.D0/3.D0*Hr3(-1,0,1)*x - 32.D0*
     &    Hr3(0,-1,0) + 80.D0/3.D0*Hr3(0,-1,0)*x + 16.D0/3.D0*Hr3(0,0,0
     &    )*x + 64.D0/3.D0*Hr3(0,0,1) + 32.D0/3.D0*Hr3(0,0,1)*x + 32.D0/
     &    3.D0*Hr3(0,1,0) - 16.D0/3.D0*Hr3(0,1,0)*x + 16.D0*Hr3(0,1,1)
     &     - 24.D0*Hr3(0,1,1)*x - 32.D0/3.D0*Hr3(1,0,0) + 16.D0/3.D0*
     &    Hr3(1,0,0)*x + 80.D0/3.D0*Hr3(1,0,1) - 40.D0/3.D0*Hr3(1,0,1)*
     &    x + 80.D0/3.D0*Hr3(1,1,0) - 40.D0/3.D0*Hr3(1,1,0)*x + 80.D0/3.
     &    D0*Hr3(1,1,1) - 40.D0/3.D0*Hr3(1,1,1)*x )
      dPgq2M = dPgq2M + nf*cf**2 * (  - 1733.D0/3.D0 + 10043.D0/18.D0
     &    *x + 496.D0/3.D0*z3 - 248.D0/3.D0*z3*x + 704.D0/9.D0*z2 + 32.D
     &    0/9.D0*z2*x - 6310.D0/27.D0*Hr1(0) - 7786.D0/27.D0*Hr1(0)*x
     &     + 160.D0/3.D0*Hr1(0)*z2 + 16.D0/3.D0*Hr1(0)*z2*x - 248.D0/27.
     &    D0*Hr1(1) - 248.D0/27.D0*Hr1(1)*x - 32.D0/3.D0*Hr1(1)*z2 + 16.
     &    D0/3.D0*Hr1(1)*z2*x + 64.D0*Hr2(-1,0) + 80.D0*Hr2(-1,0)*x - 
     &    16.D0*Hr2(-1,0)*dx - 1348.D0/9.D0*Hr2(0,0) + 914.D0/9.D0*Hr2(
     &    0,0)*x - 128.D0/9.D0*Hr2(0,1) - 32.D0/9.D0*Hr2(0,1)*x - 32.D0/
     &    3.D0*Hr2(1,0) - 8.D0/3.D0*Hr2(1,0)*x + 328.D0/9.D0*Hr2(1,1)
     &     - 32.D0/9.D0*Hr2(1,1)*x + 64.D0*Hr3(0,-1,0) - 32.D0*Hr3(0,-1
     &    ,0)*x - 400.D0/3.D0*Hr3(0,0,0) - 232.D0/3.D0*Hr3(0,0,0)*x + 
     &    32.D0/3.D0*Hr3(0,0,1) - 16.D0/3.D0*Hr3(0,0,1)*x + 32.D0/3.D0*
     &    Hr3(0,1,0) - 16.D0/3.D0*Hr3(0,1,0)*x + 16.D0/3.D0*Hr3(0,1,1)
     &     - 8.D0/3.D0*Hr3(0,1,1)*x + 16.D0*Hr3(1,0,0) - 8.D0*Hr3(1,0,0
     &    )*x + 32.D0/3.D0*Hr3(1,0,1) - 16.D0/3.D0*Hr3(1,0,1)*x + 64.D0/
     &    3.D0*Hr3(1,1,0) - 32.D0/3.D0*Hr3(1,1,0)*x - 80.D0/3.D0*Hr3(1,
     &    1,1) + 40.D0/3.D0*Hr3(1,1,1)*x - 64.D0*Hr4(0,0,0,0) + 32.D0*
     &    Hr4(0,0,0,0)*x )
      dPgq2M = dPgq2M + nf2*cf * (  - 16.D0/3.D0 + 40.D0/9.D0*x - 
     &    32.D0/9.D0*Hr1(1) - 8.D0/9.D0*Hr1(1)*x + 16.D0/3.D0*Hr2(1,1)
     &     - 8.D0/3.D0*Hr2(1,1)*x )
*
       DP2GQMA = dPgq2M
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..3-loop gluon-gluon splitting function DeltaP_gg^(2), regular part
*
       FUNCTION DP2GGMA (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,
     ,             Z4 = 1.0823 23233 71113 81915 D0,
     ,             Z5 = 1.0369 27755 14336 99263 D0 )
*
* ...The soft coefficient for use in DP2GGB and DP2GGC
*
       COMMON / P2GSOFT / A3G
*
* ...Colour factors and abbreviations
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
       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)
*
      dPgg2M =
     &  + ca**3 * ( 46102.D0/9.D0 - 47572.D0/9.D0*x - 576.D0*z4 - 680.D0
     &    *z4*x - 220.D0*dp*z4 + 490.D0/3.D0*dm - 12.D0*dm*z4 - 2320.D0/
     &    3.D0*z3 + 2792.D0/3.D0*z3*x + 440.D0/3.D0*z3*x**2 + 176.D0/3.D
     &    0*z3*dm - 9136.D0/9.D0*z2 + 304.D0*z2*x + 1072.D0/9.D0*z2*dp
     &     - 1072.D0/9.D0*z2*dm - 192.D0*Hr1(-1)*z3 - 384.D0*Hr1(-1)*z3
     &    *x - 192.D0*Hr1(-1)*z3*dp - 200.D0*Hr1(-1)*z2 - 200.D0*Hr1(-1
     &    )*z2*x - 88.D0*Hr1(-1)*z2*x**2 - 88.D0*Hr1(-1)*z2*dx + 49322.D
     &    0/27.D0*Hr1(0) + 79994.D0/27.D0*Hr1(0)*x + 8.D0/3.D0*Hr1(0)*
     &    dm - 160.D0*Hr1(0)*z3 + 80.D0*Hr1(0)*z3*dp - 112.D0*Hr1(0)*z3
     &    *dm - 1064.D0/3.D0*Hr1(0)*z2 + 3136.D0/3.D0*Hr1(0)*z2*x + 352.
     &    D0/3.D0*Hr1(0)*z2*x**2 + 88.D0/3.D0*Hr1(0)*z2*dp + 8.D0/9.D0*
     &    Hr1(1) - 8.D0/9.D0*Hr1(1)*x - 96.D0*Hr1(1)*z3 + 192.D0*Hr1(1)
     &    *z3*x - 96.D0*Hr1(1)*z3*dm - 168.D0*Hr1(1)*z2 + 168.D0*Hr1(1)
     &    *z2*x + 88.D0/3.D0*Hr1(1)*z2*x**2 - 88.D0/3.D0*Hr1(1)*z2*dx
     &     + 256.D0*Hr2(-1,-1)*z2 + 512.D0*Hr2(-1,-1)*z2*x + 256.D0*
     &    Hr2(-1,-1)*z2*dp - 8176.D0/9.D0*Hr2(-1,0) - 6032.D0/9.D0*Hr2(
     &    -1,0)*x + 2144.D0/9.D0*Hr2(-1,0)*dp - 288.D0*Hr2(-1,0)*z2 - 
     &    576.D0*Hr2(-1,0)*z2*x - 288.D0*Hr2(-1,0)*z2*dp - 96.D0*Hr2(0,
     &    -1)*z2 - 768.D0*Hr2(0,-1)*z2*x - 256.D0*Hr2(0,-1)*z2*dp + 96.D
     &    0*Hr2(0,-1)*z2*dm + 4076.D0/3.D0*Hr2(0,0) - 22156.D0/9.D0*
     &    Hr2(0,0)*x - 1072.D0/9.D0*Hr2(0,0)*dp + 1072.D0/9.D0*Hr2(0,0)
     &    *dm - 352.D0*Hr2(0,0)*z2 + 96.D0*Hr2(0,0)*z2*x + 128.D0*Hr2(0
     &    ,0)*z2*dp - 128.D0*Hr2(0,0)*z2*dm + 320.D0/3.D0*Hr2(0,1) - 
     &    304.D0*Hr2(0,1)*x + 2144.D0/9.D0*Hr2(0,1)*dm - 160.D0*Hr2(0,1
     &    )*z2 - 32.D0*Hr2(0,1)*z2*dp - 64.D0*Hr2(0,1)*z2*dm + 776.D0/9.
     &    D0*Hr2(1,0) - 2920.D0/9.D0*Hr2(1,0)*x + 2144.D0/9.D0*Hr2(1,0)
     &    *dm - 96.D0*Hr2(1,0)*z2 + 192.D0*Hr2(1,0)*z2*x - 96.D0*Hr2(1,
     &    0)*z2*dm + 336.D0*Hr3(-1,-1,0) + 336.D0*Hr3(-1,-1,0)*x - 176.D
     &    0/3.D0*Hr3(-1,-1,0)*x**2 - 176.D0/3.D0*Hr3(-1,-1,0)*dx + 1040.
     &    D0*Hr3(-1,0,0) + 1040.D0*Hr3(-1,0,0)*x + 176.D0/3.D0*Hr3(-1,0
     &    ,0)*x**2 + 176.D0/3.D0*Hr3(-1,0,0)*dx + 368.D0*Hr3(-1,0,1) + 
     &    368.D0*Hr3(-1,0,1)*x + 176.D0/3.D0*Hr3(-1,0,1)*x**2 + 176.D0/
     &    3.D0*Hr3(-1,0,1)*dx - 112.D0*Hr3(0,-1,0) + 2224.D0/3.D0*Hr3(0
     &    ,-1,0)*x + 176.D0/3.D0*Hr3(0,-1,0)*x**2 + 176.D0/3.D0*Hr3(0,
     &    -1,0)*dm + 2528.D0/3.D0*Hr3(0,0,0) - 464.D0*Hr3(0,0,0)*x - 
     &    176.D0/3.D0*Hr3(0,0,0)*x**2 + 728.D0/3.D0*Hr3(0,0,1) - 3136.D0
     &    /3.D0*Hr3(0,0,1)*x - 176.D0/3.D0*Hr3(0,0,1)*x**2 + 88.D0/3.D0
     &    *Hr3(0,0,1)*dm + 1232.D0/3.D0*Hr3(0,1,0) - 1840.D0/3.D0*Hr3(0
     &    ,1,0)*x + 1864.D0/3.D0*Hr3(1,0,0) - 1952.D0/3.D0*Hr3(1,0,0)*x
     &     + 88.D0/3.D0*Hr3(1,0,0)*dm - 384.D0*Hr4(-1,-1,0,0) - 768.D0*
     &    Hr4(-1,-1,0,0)*x - 384.D0*Hr4(-1,-1,0,0)*dp - 256.D0*Hr4(-1,
     &    -1,0,1) - 512.D0*Hr4(-1,-1,0,1)*x - 256.D0*Hr4(-1,-1,0,1)*dp
     &     - 128.D0*Hr4(-1,0,-1,0) - 256.D0*Hr4(-1,0,-1,0)*x - 128.D0*
     &    Hr4(-1,0,-1,0)*dp + 256.D0*Hr4(-1,0,0,0) + 512.D0*Hr4(-1,0,0,
     &    0)*x + 256.D0*Hr4(-1,0,0,0)*dp + 256.D0*Hr4(-1,0,0,1) + 512.D0
     &    *Hr4(-1,0,0,1)*x + 256.D0*Hr4(-1,0,0,1)*dp + 64.D0*Hr4(-1,0,1
     &    ,0) + 128.D0*Hr4(-1,0,1,0)*x + 64.D0*Hr4(-1,0,1,0)*dp + 64.D0
     &    *Hr4(0,-1,-1,0) - 512.D0*Hr4(0,-1,-1,0)*x - 128.D0*Hr4(0,-1,
     &    -1,0)*dp + 64.D0*Hr4(0,-1,-1,0)*dm + 416.D0*Hr4(0,-1,0,0) + 
     &    512.D0*Hr4(0,-1,0,0)*x + 288.D0*Hr4(0,-1,0,0)*dp - 64.D0*Hr4(
     &    0,-1,0,0)*dm + 128.D0*Hr4(0,-1,0,1) + 512.D0*Hr4(0,-1,0,1)*x
     &     + 192.D0*Hr4(0,-1,0,1)*dp - 64.D0*Hr4(0,-1,0,1)*dm + 64.D0*
     &    Hr4(0,0,-1,0) - 64.D0*Hr4(0,0,-1,0)*x + 64.D0*Hr4(0,0,-1,0)*
     &    dp - 64.D0*Hr4(0,0,-1,0)*dm + 448.D0*Hr4(0,0,0,0) - 256.D0*
     &    Hr4(0,0,0,0)*x - 64.D0*Hr4(0,0,0,0)*dp + 64.D0*Hr4(0,0,0,0)*
     &    dm + 416.D0*Hr4(0,0,0,1) - 96.D0*Hr4(0,0,0,1)*x - 128.D0*Hr4(
     &    0,0,0,1)*dp + 128.D0*Hr4(0,0,0,1)*dm + 384.D0*Hr4(0,0,1,0) - 
     &    128.D0*Hr4(0,0,1,0)*x - 32.D0*Hr4(0,0,1,0)*dp + 160.D0*Hr4(0,
     &    0,1,0)*dm + 128.D0*Hr4(0,0,1,1) - 256.D0*Hr4(0,0,1,1)*x + 128.
     &    D0*Hr4(0,0,1,1)*dm + 320.D0*Hr4(0,1,0,0) - 160.D0*Hr4(0,1,0,0
     &    )*x + 160.D0*Hr4(0,1,0,0)*dm + 128.D0*Hr4(0,1,0,1) - 256.D0*
     &    Hr4(0,1,0,1)*x + 128.D0*Hr4(0,1,0,1)*dm + 128.D0*Hr4(0,1,1,0)
     &     - 256.D0*Hr4(0,1,1,0)*x + 128.D0*Hr4(0,1,1,0)*dm - 64.D0*
     &    Hr4(1,0,-1,0) + 128.D0*Hr4(1,0,-1,0)*x - 64.D0*Hr4(1,0,-1,0)*
     &    dm + 128.D0*Hr4(1,0,0,0) - 256.D0*Hr4(1,0,0,0)*x + 128.D0*
     &    Hr4(1,0,0,0)*dm + 128.D0*Hr4(1,0,0,1) - 256.D0*Hr4(1,0,0,1)*x
     &     + 128.D0*Hr4(1,0,0,1)*dm + 128.D0*Hr4(1,0,1,0) - 256.D0*Hr4(
     &    1,0,1,0)*x + 128.D0*Hr4(1,0,1,0)*dm + 128.D0*Hr4(1,1,0,0) - 
     &    256.D0*Hr4(1,1,0,0)*x + 128.D0*Hr4(1,1,0,0)*dm )
      dPgg2M = dPgg2M + nf*ca**2 * ( 23740.D0/27.D0 - 22904.D0/27.D0*
     &    x + 78.D0*z4 + 138.D0*z4*x - 836.D0/27.D0*dm - 752.D0/3.D0*z3
     &     + 496.D0/3.D0*z3*x + 80.D0/3.D0*z3*x**2 - 128.D0/3.D0*z3*dm
     &     - 1996.D0/9.D0*z2 - 1844.D0/9.D0*z2*x + 16.D0*z2*x**2 - 160.D
     &    0/9.D0*z2*dp + 160.D0/9.D0*z2*dm + 56.D0*Hr1(-1)*z2 + 56.D0*
     &    Hr1(-1)*z2*x - 16.D0*Hr1(-1)*z2*x**2 - 16.D0*Hr1(-1)*z2*dx + 
     &    4772.D0/27.D0*Hr1(0) + 15362.D0/27.D0*Hr1(0)*x - 8.D0/3.D0*
     &    Hr1(0)*dm - 96.D0*Hr1(0)*z3 - 220.D0/3.D0*Hr1(0)*z2 + 212.D0/
     &    3.D0*Hr1(0)*z2*x + 64.D0/3.D0*Hr1(0)*z2*x**2 - 16.D0/3.D0*
     &    Hr1(0)*z2*dp + 2174.D0/9.D0*Hr1(1) - 2174.D0/9.D0*Hr1(1)*x - 
     &    88.D0*Hr1(1)*z2 + 88.D0*Hr1(1)*z2*x + 16.D0/3.D0*Hr1(1)*z2*
     &    x**2 - 16.D0/3.D0*Hr1(1)*z2*dx - 1688.D0/9.D0*Hr2(-1,0) - 
     &    2008.D0/9.D0*Hr2(-1,0)*x + 16.D0*Hr2(-1,0)*x**2 + 16.D0*Hr2(
     &    -1,0)*dx - 320.D0/9.D0*Hr2(-1,0)*dp + 48.D0*Hr2(0,-1)*z2 - 48.
     &    D0*Hr2(0,-1)*z2*x - 584.D0/9.D0*Hr2(0,0) - 268.D0/9.D0*Hr2(0,
     &    0)*x - 16.D0*Hr2(0,0)*x**2 + 160.D0/9.D0*Hr2(0,0)*dp - 160.D0/
     &    9.D0*Hr2(0,0)*dm - 40.D0*Hr2(0,0)*z2 - 24.D0*Hr2(0,0)*z2*x + 
     &    308.D0/9.D0*Hr2(0,1) + 1844.D0/9.D0*Hr2(0,1)*x - 320.D0/9.D0*
     &    Hr2(0,1)*dm - 48.D0*Hr2(0,1)*z2 - 48.D0*Hr2(0,1)*z2*x - 608.D0
     &    /9.D0*Hr2(1,0) + 928.D0/9.D0*Hr2(1,0)*x - 320.D0/9.D0*Hr2(1,0
     &    )*dm + 176.D0*Hr3(-1,-1,0) + 176.D0*Hr3(-1,-1,0)*x - 32.D0/3.D
     &    0*Hr3(-1,-1,0)*x**2 - 32.D0/3.D0*Hr3(-1,-1,0)*dx - 72.D0*Hr3(
     &    -1,0,0) - 72.D0*Hr3(-1,0,0)*x + 32.D0/3.D0*Hr3(-1,0,0)*x**2
     &     + 32.D0/3.D0*Hr3(-1,0,0)*dx + 32.D0*Hr3(-1,0,1) + 32.D0*Hr3(
     &    -1,0,1)*x + 32.D0/3.D0*Hr3(-1,0,1)*x**2 + 32.D0/3.D0*Hr3(-1,0
     &    ,1)*dx - 72.D0*Hr3(0,-1,0) - 280.D0/3.D0*Hr3(0,-1,0)*x + 32.D0
     &    /3.D0*Hr3(0,-1,0)*x**2 - 32.D0/3.D0*Hr3(0,-1,0)*dm - 8.D0/3.D0
     &    *Hr3(0,0,0) + 32.D0*Hr3(0,0,0)*x - 32.D0/3.D0*Hr3(0,0,0)*x**2
     &     + 4.D0/3.D0*Hr3(0,0,1) - 212.D0/3.D0*Hr3(0,0,1)*x - 32.D0/3.D
     &    0*Hr3(0,0,1)*x**2 - 16.D0/3.D0*Hr3(0,0,1)*dm - 32.D0/3.D0*
     &    Hr3(0,1,0) - 32.D0/3.D0*Hr3(0,1,0)*x + 116.D0/3.D0*Hr3(1,0,0)
     &     - 100.D0/3.D0*Hr3(1,0,0)*x - 16.D0/3.D0*Hr3(1,0,0)*dm + 96.D0
     &    *Hr4(0,-1,-1,0) - 96.D0*Hr4(0,-1,-1,0)*x - 48.D0*Hr4(0,-1,0,0
     &    ) + 48.D0*Hr4(0,-1,0,0)*x - 16.D0*Hr4(0,0,-1,0) + 80.D0*Hr4(0
     &    ,0,-1,0)*x + 16.D0*Hr4(0,0,0,0) + 24.D0*Hr4(0,0,0,1) + 24.D0*
     &    Hr4(0,0,0,1)*x + 24.D0*Hr4(0,1,0,0) + 24.D0*Hr4(0,1,0,0)*x )
      dPgg2M = dPgg2M + nf*cf*ca * (  - 61046.D0/27.D0 + 62036.D0/27.D
     &    0*x - 172.D0*z4 - 332.D0*z4*x - 110.D0/3.D0*dm + 1096.D0/3.D0
     &    *z3 + 112.D0/3.D0*z3*x - 160.D0/3.D0*z3*x**2 + 32.D0*z3*dm + 
     &    3392.D0/9.D0*z2 + 1364.D0/9.D0*z2*x - 48.D0*z2*x**2 - 176.D0*
     &    Hr1(-1)*z2 - 176.D0*Hr1(-1)*z2*x + 32.D0*Hr1(-1)*z2*x**2 + 32.
     &    D0*Hr1(-1)*z2*dx - 33452.D0/27.D0*Hr1(0) - 13448.D0/27.D0*
     &    Hr1(0)*x + 144.D0*Hr1(0)*z3 - 16.D0*Hr1(0)*z3*x + 448.D0/3.D0
     &    *Hr1(0)*z2 - 344.D0/3.D0*Hr1(0)*z2*x - 128.D0/3.D0*Hr1(0)*z2*
     &    x**2 - 4508.D0/9.D0*Hr1(1) + 4508.D0/9.D0*Hr1(1)*x + 240.D0*
     &    Hr1(1)*z2 - 240.D0*Hr1(1)*z2*x - 32.D0/3.D0*Hr1(1)*z2*x**2 + 
     &    32.D0/3.D0*Hr1(1)*z2*dx + 488.D0/3.D0*Hr2(-1,0) + 488.D0/3.D0
     &    *Hr2(-1,0)*x - 48.D0*Hr2(-1,0)*x**2 - 48.D0*Hr2(-1,0)*dx - 
     &    128.D0*Hr2(0,-1)*z2 + 128.D0*Hr2(0,-1)*z2*x - 6892.D0/9.D0*
     &    Hr2(0,0) + 2768.D0/9.D0*Hr2(0,0)*x + 48.D0*Hr2(0,0)*x**2 + 96.
     &    D0*Hr2(0,0)*z2 + 96.D0*Hr2(0,0)*z2*x - 1928.D0/9.D0*Hr2(0,1)
     &     - 1364.D0/9.D0*Hr2(0,1)*x + 128.D0*Hr2(0,1)*z2 + 128.D0*Hr2(
     &    0,1)*z2*x - 208.D0/3.D0*Hr2(1,0) + 208.D0/3.D0*Hr2(1,0)*x - 
     &    260.D0/3.D0*Hr2(1,1) + 260.D0/3.D0*Hr2(1,1)*x - 480.D0*Hr3(-1
     &    ,-1,0) - 480.D0*Hr3(-1,-1,0)*x + 64.D0/3.D0*Hr3(-1,-1,0)*x**2
     &     + 64.D0/3.D0*Hr3(-1,-1,0)*dx + 48.D0*Hr3(-1,0,0) + 48.D0*
     &    Hr3(-1,0,0)*x - 64.D0/3.D0*Hr3(-1,0,0)*x**2 - 64.D0/3.D0*Hr3(
     &    -1,0,0)*dx - 64.D0*Hr3(-1,0,1) - 64.D0*Hr3(-1,0,1)*x - 64.D0/
     &    3.D0*Hr3(-1,0,1)*x**2 - 64.D0/3.D0*Hr3(-1,0,1)*dx + 64.D0*
     &    Hr3(0,-1,0) + 256.D0*Hr3(0,-1,0)*x - 64.D0/3.D0*Hr3(0,-1,0)*
     &    x**2 - 1448.D0/3.D0*Hr3(0,0,0) - 704.D0/3.D0*Hr3(0,0,0)*x + 
     &    64.D0/3.D0*Hr3(0,0,0)*x**2 - 256.D0/3.D0*Hr3(0,0,1) + 344.D0/
     &    3.D0*Hr3(0,0,1)*x + 64.D0/3.D0*Hr3(0,0,1)*x**2 - 304.D0/3.D0*
     &    Hr3(0,1,0) + 296.D0/3.D0*Hr3(0,1,0)*x - 152.D0/3.D0*Hr3(0,1,1
     &    ) + 88.D0/3.D0*Hr3(0,1,1)*x - 232.D0*Hr3(1,0,0) + 232.D0*Hr3(
     &    1,0,0)*x - 80.D0*Hr3(1,1,0) + 80.D0*Hr3(1,1,0)*x - 80.D0*Hr3(
     &    1,1,1) + 80.D0*Hr3(1,1,1)*x - 256.D0*Hr4(0,-1,-1,0) + 256.D0*
     &    Hr4(0,-1,-1,0)*x + 64.D0*Hr4(0,-1,0,0) - 64.D0*Hr4(0,-1,0,0)*
     &    x - 128.D0*Hr4(0,0,-1,0)*x - 192.D0*Hr4(0,0,0,0) + 128.D0*
     &    Hr4(0,0,0,0)*x - 96.D0*Hr4(0,0,0,1) - 96.D0*Hr4(0,0,0,1)*x - 
     &    64.D0*Hr4(0,0,1,0) - 64.D0*Hr4(0,0,1,0)*x - 16.D0*Hr4(0,0,1,1
     &    ) - 16.D0*Hr4(0,0,1,1)*x - 112.D0*Hr4(0,1,0,0) - 112.D0*Hr4(0
     &    ,1,0,0)*x - 32.D0*Hr4(0,1,1,0) - 32.D0*Hr4(0,1,1,0)*x - 32.D0
     &    *Hr4(0,1,1,1) - 32.D0*Hr4(0,1,1,1)*x )
      dPgg2M = dPgg2M + nf*cf**2 * ( 64.D0 - 64.D0*x + 180.D0*z4 + 
     &    260.D0*z4*x - 80.D0*z3 - 336.D0*z3*x + 64.D0*z2 - 76.D0*z2*x
     &     + 32.D0*z2*x**2 + 128.D0*Hr1(-1)*z2 + 128.D0*Hr1(-1)*z2*x - 
     &    100.D0*Hr1(0) + 40.D0*Hr1(0)*x - 16.D0*Hr1(0)*z3 + 112.D0*
     &    Hr1(0)*z3*x + 64.D0*Hr1(0)*z2 - 16.D0*Hr1(0)*z2*x + 164.D0*
     &    Hr1(1) - 164.D0*Hr1(1)*x - 128.D0*Hr1(1)*z2 + 128.D0*Hr1(1)*
     &    z2*x + 80.D0*Hr2(-1,0) + 80.D0*Hr2(-1,0)*x + 32.D0*Hr2(-1,0)*
     &    x**2 + 32.D0*Hr2(-1,0)*dx + 64.D0*Hr2(0,-1)*z2 - 64.D0*Hr2(0,
     &    -1)*z2*x - 184.D0*Hr2(0,0) - 12.D0*Hr2(0,0)*x - 32.D0*Hr2(0,0
     &    )*x**2 + 32.D0*Hr2(0,0)*z2 + 32.D0*Hr2(0,0)*z2*x + 16.D0*Hr2(
     &    0,1) + 76.D0*Hr2(0,1)*x - 64.D0*Hr2(0,1)*z2 - 64.D0*Hr2(0,1)*
     &    z2*x - 56.D0*Hr2(1,0) + 56.D0*Hr2(1,0)*x + 52.D0*Hr2(1,1) - 
     &    52.D0*Hr2(1,1)*x + 256.D0*Hr3(-1,-1,0) + 256.D0*Hr3(-1,-1,0)*
     &    x - 128.D0*Hr3(-1,0,0) - 128.D0*Hr3(-1,0,0)*x - 256.D0*Hr3(0,
     &    -1,0)*x - 112.D0*Hr3(0,0,0) + 48.D0*Hr3(0,0,0)*x - 64.D0*Hr3(
     &    0,0,1) + 16.D0*Hr3(0,0,1)*x - 48.D0*Hr3(0,1,0) - 64.D0*Hr3(0,
     &    1,0)*x + 32.D0*Hr3(0,1,1) - 48.D0*Hr3(0,1,1)*x + 88.D0*Hr3(1,
     &    0,0) - 88.D0*Hr3(1,0,0)*x + 80.D0*Hr3(1,1,0) - 80.D0*Hr3(1,1,
     &    0)*x + 80.D0*Hr3(1,1,1) - 80.D0*Hr3(1,1,1)*x + 128.D0*Hr4(0,
     &    -1,-1,0) - 128.D0*Hr4(0,-1,-1,0)*x - 64.D0*Hr4(0,-1,0,0) + 64.
     &    D0*Hr4(0,-1,0,0)*x + 128.D0*Hr4(0,0,-1,0)*x - 32.D0*Hr4(0,0,0
     &    ,0) - 32.D0*Hr4(0,0,0,0)*x - 32.D0*Hr4(0,0,0,1) - 32.D0*Hr4(0
     &    ,0,0,1)*x + 16.D0*Hr4(0,0,1,1) + 16.D0*Hr4(0,0,1,1)*x + 48.D0
     &    *Hr4(0,1,0,0) + 48.D0*Hr4(0,1,0,0)*x + 32.D0*Hr4(0,1,1,0) + 
     &    32.D0*Hr4(0,1,1,0)*x + 32.D0*Hr4(0,1,1,1) + 32.D0*Hr4(0,1,1,1
     &    )*x )
      dPgg2M = dPgg2M + nf2*ca * (  - 214.D0/27.D0 + 230.D0/27.D0*x
     &     - 16.D0/27.D0*dm + 32.D0/9.D0*z2 + 32.D0/9.D0*z2*x - 112.D0/
     &    27.D0*Hr1(0) + 134.D0/27.D0*Hr1(0)*x - 82.D0/9.D0*Hr1(1) + 82.
     &    D0/9.D0*Hr1(1)*x - 16.D0/9.D0*Hr2(0,0) - 16.D0/9.D0*Hr2(0,0)*
     &    x - 32.D0/9.D0*Hr2(0,1) - 32.D0/9.D0*Hr2(0,1)*x )
      dPgg2M = dPgg2M + nf2*cf * (  - 1168.D0/27.D0 + 1168.D0/27.D0
     &    *x - 16.D0/3.D0*z3 - 16.D0/3.D0*z3*x + 208.D0/9.D0*z2 - 128.D0
     &    /9.D0*z2*x - 520.D0/27.D0*Hr1(0) + 104.D0/27.D0*Hr1(0)*x + 32.
     &    D0/3.D0*Hr1(0)*z2 + 32.D0/3.D0*Hr1(0)*z2*x - 208.D0/9.D0*Hr1(
     &    1) + 208.D0/9.D0*Hr1(1)*x - 56.D0/9.D0*Hr2(0,0) + 184.D0/9.D0
     &    *Hr2(0,0)*x - 208.D0/9.D0*Hr2(0,1) + 128.D0/9.D0*Hr2(0,1)*x
     &     - 80.D0/3.D0*Hr2(1,0) + 80.D0/3.D0*Hr2(1,0)*x - 40.D0/3.D0*
     &    Hr2(1,1) + 40.D0/3.D0*Hr2(1,1)*x - 16.D0/3.D0*Hr3(0,0,0) - 16.
     &    D0/3.D0*Hr3(0,0,0)*x - 32.D0/3.D0*Hr3(0,0,1) - 32.D0/3.D0*
     &    Hr3(0,0,1)*x - 32.D0/3.D0*Hr3(0,1,0) - 32.D0/3.D0*Hr3(0,1,0)*
     &    x - 16.D0/3.D0*Hr3(0,1,1) - 16.D0/3.D0*Hr3(0,1,1)*x )
*
* ...The soft (`+'-distribution) part of the splitting function
*
       A3G =
     ,      ca**3    * ( 490.D0/3.D0 - 1072.D0/9.D0*z2 + 88.D0/3.D0*z3
     ,                 + 176.D0/5.D0*z2**2 )
     ,    + ca**2*nf * ( - 836./27.D0 + 160./9.D0*z2 - 112./3.D0*z3 )
     ,    + ca*cf*nf * ( - 110./3.D0 + 32.*z3 ) - ca*nf2 * 16./27.D0
*
       dPGG2S = DM * A3G
*
* ...The regular piece of the coefficient function
*
       DP2GGMA = dPgg2M - DPGG2S
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The singular (soft) piece of DeltaP_gg^(2)
*
       FUNCTION DP2GGB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P2GSOFT / A3G
*
       DP2GGB  = A3G/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
* ..The 'local' piece of DeltaP_gg^(2)
*
       FUNCTION DP2GGC (Y, NF)
*
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF, NF2
       COMMON / P2GSOFT / A3G
       PARAMETER ( Z2 = 1.6449 34066 84822 64365 D0,
     ,             Z3 = 1.2020 56903 15959 42854 D0,
     ,             Z4 = 1.0823 23233 71113 81915 D0,
     ,             Z5 = 1.0369 27755 14336 99263 D0 )
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
* ...The coefficient of delta(1-x)
*
       P2GDELT =
     ,    + 79.D0/2.D0*ca**3
     ,    + 8.D0/3.D0*z2*ca**3
     ,    + 22.D0/3.D0*z2**2*ca**3
     ,    + 536.D0/3.D0*z3*ca**3
     ,    - 16.D0*z2*z3*ca**3
     ,    - 80.D0*z5*ca**3
     ,    + cf**2*nf
     ,    - 233.D0/18*ca**2*nf
     ,    - 8.D0/3.D0*z2*ca**2*nf
     ,    - 80.D0/3.D0*z3*ca**2*nf
     ,    - 4.D0/3.D0*z2**2*ca**2*nf
     ,    - 241.D0/18.D0*ca*cf*nf
     ,    + 29.D0/18.D0*ca*nf2
     ,    + 11.D0/9.D0*cf*nf2
*
       DP2GGC = DLOG (1.D0-Y) * A3G + P2GDELT
*
       RETURN
       END
*
* =================================================================av==
