*
* ..File: xpij2te.f 
*
*
* ..The exact 3-loop MS(bar) singlet splitting functions  P_ij^(2)T
*    for the evolution of unpolarized fragmentation densities at
*    mu_r = mu_f.  The expansion parameter is alpha_s/(4 pi).
*
* ..The function  X2QGTA  includes a parameter IMOD for estimating, 
*    via  IMOD = 1 and IMOD = -1, a residual uncertainty discussed in 
*    the AMV paper.  IMOD = 0  provides the main result.
* 
* ..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.
*
* ..References: S. Moch and A. Vogt, 
*               Phys. Lett. B659 (2008) 290, arXiv:0708.3899  (ps, gg)
*               A. Almasy, S. Moch and A. Vogt,
*               Nucl. Phys. B854 (2012) 133, arXiv:1107.2263  (qg, gq)
*
* =====================================================================
*
*
* ..The (regular) pure-singlet splitting functions P_ps^(2)T.
*    P_qq^(2) is obtained by adding the non-singlet quantity P_NS^(2)T+
*
       FUNCTION X2PSTA (X, NF)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N, IMOD
       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 )
*
* ...Colour factors and an abbreviation
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
       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
*
      PTqqps2 =
     &  + nf*cf*ca * ( 820.D0/3.D0 - 1172.D0/3.D0*x + 7496.D0/81.D0*
     &    x**2 + 2008.D0/81.D0*dx + 296.D0/3.D0*z3 - 1312.D0/3.D0*z3*x
     &     + 256.D0/3.D0*z3*x**2 + 2212.D0/9.D0*z2 + 4912.D0/9.D0*z2*x
     &     + 640.D0/9.D0*z2*x**2 + 512.D0/9.D0*z2*dx - 852.D0/5.D0*
     &    z2**2 - 812.D0/5.D0*z2**2*x + 40.D0*Hr1(-1)*z2 + 40.D0*Hr1(-1
     &    )*z2*x + 32.D0/3.D0*Hr1(-1)*z2*x**2 + 32.D0/3.D0*Hr1(-1)*z2*
     &    dx - 2380.D0/27.D0*Hr1(0) + 23540.D0/27.D0*Hr1(0)*x + 9172.D0/
     &    27.D0*Hr1(0)*x**2 + 1240.D0/27.D0*Hr1(0)*dx - 160.D0*Hr1(0)*
     &    z3 - 32.D0*Hr1(0)*z3*x + 196.D0/3.D0*Hr1(0)*z2 - 404.D0/3.D0*
     &    Hr1(0)*z2*x + 32.D0/3.D0*Hr1(0)*z2*x**2 + 64.D0/3.D0*Hr1(0)*
     &    z2*dx - 2596.D0/9.D0*Hr1(1) + 592.D0/9.D0*Hr1(1)*x + 7868.D0/
     &    27.D0*Hr1(1)*x**2 - 1856.D0/27.D0*Hr1(1)*dx - 72.D0*Hr1(1)*z2
     &     + 72.D0*Hr1(1)*z2*x + 96.D0*Hr1(1)*z2*x**2 - 96.D0*Hr1(1)*z2
     &    *dx + 688.D0/3.D0*Hr2(-1,0) + 400.D0/3.D0*Hr2(-1,0)*x - 992.D0
     &    /9.D0*Hr2(-1,0)*x**2 - 128.D0/9.D0*Hr2(-1,0)*dx + 16.D0*Hr2(0
     &    ,-1)*z2 )
      PTqqps2 = PTqqps2 + nf*cf*ca * (  - 16.D0*Hr2(0,-1)*z2*x - 2440.D0
     &    /9.D0*Hr2(0,0) - 10468.D0/9.D0*Hr2(0,0)*x - 2312.D0/9.D0*Hr2(
     &    0,0)*x**2 - 64.D0/9.D0*Hr2(0,0)*dx - 8.D0*Hr2(0,0)*z2 + 104.D0
     &    *Hr2(0,0)*z2*x - 1540.D0/9.D0*Hr2(0,1) - 3520.D0/9.D0*Hr2(0,1
     &    )*x - 2032.D0/9.D0*Hr2(0,1)*x**2 - 112.D0/9.D0*Hr2(0,1)*dx - 
     &    144.D0*Hr2(0,1)*z2 - 144.D0*Hr2(0,1)*z2*x + 4.D0/3.D0*Hr2(1,0
     &    ) + 140.D0/3.D0*Hr2(1,0)*x - 72.D0*Hr2(1,0)*x**2 + 24.D0*Hr2(
     &    1,0)*dx - 12.D0*Hr2(1,1) + 12.D0*Hr2(1,1)*x + 16.D0/9.D0*Hr2(
     &    1,1)*x**2 - 16.D0/9.D0*Hr2(1,1)*dx + 16.D0*Hr3(-1,-1,0) + 16.D
     &    0*Hr3(-1,-1,0)*x - 64.D0/3.D0*Hr3(-1,-1,0)*x**2 - 64.D0/3.D0*
     &    Hr3(-1,-1,0)*dx - 56.D0*Hr3(-1,0,0) - 56.D0*Hr3(-1,0,0)*x + 
     &    128.D0/3.D0*Hr3(-1,0,0)*x**2 + 128.D0/3.D0*Hr3(-1,0,0)*dx - 
     &    32.D0*Hr3(-1,0,1) - 32.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 + 120.D0*Hr3(0,-1,0) - 
     &    232.D0*Hr3(0,-1,0)*x + 256.D0/3.D0*Hr3(0,-1,0)*x**2 + 64.D0/3.
     &    D0*Hr3(0,-1,0)*dx )
      PTqqps2 = PTqqps2 + nf*cf*ca * (  - 112.D0/3.D0*Hr3(0,0,0) + 2288.
     &    D0/3.D0*Hr3(0,0,0)*x - 128.D0/3.D0*Hr3(0,0,0)*dx + 236.D0/3.D0
     &    *Hr3(0,0,1) + 380.D0/3.D0*Hr3(0,0,1)*x + 224.D0/3.D0*Hr3(0,0,
     &    1)*x**2 + 64.D0/3.D0*Hr3(0,0,1)*dx - 64.D0*Hr3(0,1,0) + 24.D0
     &    *Hr3(0,1,0)*x + 224.D0/3.D0*Hr3(0,1,0)*x**2 - 64.D0/3.D0*Hr3(
     &    0,1,0)*dx - 136.D0/3.D0*Hr3(0,1,1) - 88.D0/3.D0*Hr3(0,1,1)*x
     &     + 32.D0/3.D0*Hr3(0,1,1)*x**2 - 32.D0/3.D0*Hr3(0,1,1)*dx + 60.
     &    D0*Hr3(1,0,0) - 60.D0*Hr3(1,0,0)*x - 96.D0*Hr3(1,0,0)*x**2 + 
     &    96.D0*Hr3(1,0,0)*dx + 48.D0*Hr3(1,0,1) - 48.D0*Hr3(1,0,1)*x
     &     - 64.D0*Hr3(1,0,1)*x**2 + 64.D0*Hr3(1,0,1)*dx + 32.D0*Hr3(1,
     &    1,0) - 32.D0*Hr3(1,1,0)*x - 128.D0/3.D0*Hr3(1,1,0)*x**2 + 128.
     &    D0/3.D0*Hr3(1,1,0)*dx + 16.D0*Hr3(1,1,1) - 16.D0*Hr3(1,1,1)*x
     &     - 64.D0/3.D0*Hr3(1,1,1)*x**2 + 64.D0/3.D0*Hr3(1,1,1)*dx + 32.
     &    D0*Hr4(0,-1,-1,0) - 32.D0*Hr4(0,-1,-1,0)*x - 80.D0*Hr4(0,-1,0
     &    ,0) + 80.D0*Hr4(0,-1,0,0)*x - 48.D0*Hr4(0,0,-1,0) + 112.D0*
     &    Hr4(0,0,-1,0)*x )
      PTqqps2 = PTqqps2 + nf*cf*ca * ( 64.D0*Hr4(0,0,0,0) - 304.D0*Hr4(
     &    0,0,0,0)*x + 8.D0*Hr4(0,0,0,1) + 8.D0*Hr4(0,0,0,1)*x - 144.D0
     &    *Hr4(0,0,1,0) - 144.D0*Hr4(0,0,1,0)*x - 48.D0*Hr4(0,0,1,1) - 
     &    48.D0*Hr4(0,0,1,1)*x + 136.D0*Hr4(0,1,0,0) + 136.D0*Hr4(0,1,0
     &    ,0)*x + 96.D0*Hr4(0,1,0,1) + 96.D0*Hr4(0,1,0,1)*x + 64.D0*
     &    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 )
      PTqqps2 = PTqqps2 + nf*cf**2 * (  - 74.D0/27.D0 + 6266.D0/27.D0*x
     &     - 156.D0*x**2 - 220.D0/3.D0*dx + 288.D0*z3 + 352.D0*z3*x + 
     &    256.D0/3.D0*z3*x**2 + 64.D0*z3*dx - 200.D0/3.D0*z2 - 1444.D0/
     &    3.D0*z2*x - 2624.D0/9.D0*z2*x**2 + 984.D0/5.D0*z2**2 + 984.D0/
     &    5.D0*z2**2*x - 598.D0/9.D0*Hr1(0) - 4406.D0/9.D0*Hr1(0)*x - 
     &    9836.D0/27.D0*Hr1(0)*x**2 + 160.D0*Hr1(0)*z3 + 160.D0*Hr1(0)*
     &    z3*x + 80.D0*Hr1(0)*z2*x + 128.D0/3.D0*Hr1(0)*z2*x**2 + 856.D0
     &    /3.D0*Hr1(1) - 140.D0/3.D0*Hr1(1)*x - 8972.D0/27.D0*Hr1(1)*
     &    x**2 + 2528.D0/27.D0*Hr1(1)*dx + 80.D0*Hr1(1)*z2 - 80.D0*Hr1(
     &    1)*z2*x - 320.D0/3.D0*Hr1(1)*z2*x**2 + 320.D0/3.D0*Hr1(1)*z2*
     &    dx + 976.D0/3.D0*Hr2(0,0) + 1184.D0/3.D0*Hr2(0,0)*x + 608.D0/
     &    9.D0*Hr2(0,0)*x**2 - 16.D0*Hr2(0,0)*z2 - 16.D0*Hr2(0,0)*z2*x
     &     + 296.D0*Hr2(0,1) + 348.D0*Hr2(0,1)*x + 512.D0/3.D0*Hr2(0,1)
     &    *x**2 + 224.D0/9.D0*Hr2(0,1)*dx + 160.D0*Hr2(0,1)*z2 + 160.D0
     &    *Hr2(0,1)*z2*x + 428.D0/3.D0*Hr2(1,0) - 284.D0/3.D0*Hr2(1,0)*
     &    x )
      PTqqps2 = PTqqps2 + nf*cf**2 * (  - 1120.D0/9.D0*Hr2(1,0)*x**2 + 
     &    688.D0/9.D0*Hr2(1,0)*dx + 28.D0/3.D0*Hr2(1,1) - 28.D0/3.D0*
     &    Hr2(1,1)*x - 64.D0/3.D0*Hr2(1,1)*x**2 + 64.D0/3.D0*Hr2(1,1)*
     &    dx + 72.D0*Hr3(0,0,0) + 184.D0*Hr3(0,0,0)*x + 160.D0/3.D0*
     &    Hr3(0,0,0)*x**2 + 32.D0*Hr3(0,0,1) + 80.D0*Hr3(0,0,1)*x - 128.
     &    D0/3.D0*Hr3(0,0,1)*dx + 208.D0*Hr3(0,1,0) + 224.D0*Hr3(0,1,0)
     &    *x + 64.D0/3.D0*Hr3(0,1,0)*x**2 + 128.D0/3.D0*Hr3(0,1,0)*dx
     &     + 64.D0*Hr3(0,1,1) + 48.D0*Hr3(0,1,1)*x + 64.D0/3.D0*Hr3(0,1
     &    ,1)*dx - 88.D0*Hr3(1,0,0) + 88.D0*Hr3(1,0,0)*x + 352.D0/3.D0*
     &    Hr3(1,0,0)*x**2 - 352.D0/3.D0*Hr3(1,0,0)*dx - 48.D0*Hr3(1,0,1
     &    ) + 48.D0*Hr3(1,0,1)*x + 64.D0*Hr3(1,0,1)*x**2 - 64.D0*Hr3(1,
     &    0,1)*dx - 32.D0*Hr3(1,1,0) + 32.D0*Hr3(1,1,0)*x + 128.D0/3.D0
     &    *Hr3(1,1,0)*x**2 - 128.D0/3.D0*Hr3(1,1,0)*dx - 16.D0*Hr3(1,1,
     &    1) + 16.D0*Hr3(1,1,1)*x + 64.D0/3.D0*Hr3(1,1,1)*x**2 - 64.D0/
     &    3.D0*Hr3(1,1,1)*dx - 64.D0*Hr4(0,0,0,0) - 64.D0*Hr4(0,0,0,0)*
     &    x )
      PTqqps2 = PTqqps2 + nf*cf**2 * (  - 112.D0*Hr4(0,0,0,1) - 112.D0*
     &    Hr4(0,0,0,1)*x + 80.D0*Hr4(0,0,1,0) + 80.D0*Hr4(0,0,1,0)*x + 
     &    48.D0*Hr4(0,0,1,1) + 48.D0*Hr4(0,0,1,1)*x - 176.D0*Hr4(0,1,0,
     &    0) - 176.D0*Hr4(0,1,0,0)*x - 96.D0*Hr4(0,1,0,1) - 96.D0*Hr4(0
     &    ,1,0,1)*x - 64.D0*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 )
      PTqqps2 = PTqqps2 + nf2*cf * (  - 344.D0/27.D0 - 376.D0/27.D0*x
     &     + 752.D0/27.D0*x**2 - 32.D0/27.D0*dx + 64.D0/3.D0*z3 + 64.D0/
     &    3.D0*z3*x + 296.D0/9.D0*z2 + 440.D0/9.D0*z2*x + 32.D0/3.D0*z2
     &    *x**2 + 736.D0/27.D0*Hr1(0) + 304.D0/27.D0*Hr1(0)*x - 16.D0*
     &    Hr1(0)*x**2 - 16.D0/3.D0*Hr1(0)*z2 - 16.D0/3.D0*Hr1(0)*z2*x
     &     - 128.D0/3.D0*Hr1(1) + 80.D0/3.D0*Hr1(1)*x + 752.D0/27.D0*
     &    Hr1(1)*x**2 - 320.D0/27.D0*Hr1(1)*dx + 64.D0/9.D0*Hr2(0,0) + 
     &    208.D0/9.D0*Hr2(0,0)*x + 32.D0/3.D0*Hr2(0,0)*x**2 - 296.D0/9.D
     &    0*Hr2(0,1) - 440.D0/9.D0*Hr2(0,1)*x - 32.D0/3.D0*Hr2(0,1)*
     &    x**2 + 8.D0*Hr2(1,0) - 8.D0*Hr2(1,0)*x - 32.D0/3.D0*Hr2(1,0)*
     &    x**2 + 32.D0/3.D0*Hr2(1,0)*dx + 8.D0/3.D0*Hr2(1,1) - 8.D0/3.D0
     &    *Hr2(1,1)*x - 32.D0/9.D0*Hr2(1,1)*x**2 + 32.D0/9.D0*Hr2(1,1)*
     &    dx - 32.D0/3.D0*Hr3(0,0,0) - 32.D0/3.D0*Hr3(0,0,0)*x + 16.D0/
     &    3.D0*Hr3(0,0,1) + 16.D0/3.D0*Hr3(0,0,1)*x + 16.D0*Hr3(0,1,0)
     &     + 16.D0*Hr3(0,1,0)*x + 16.D0/3.D0*Hr3(0,1,1) + 16.D0/3.D0*
     &    Hr3(0,1,1)*x )

*
       X2PSTA = PTqqps2
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The quark-gluon splitting functions P_qg^(2)T
*
       FUNCTION X2QGTA (X, NF, IMOD)
*
       IMPLICIT REAL*8 (A - Z)
       COMPLEX*16 HC1, HC2, HC3, HC4 
       INTEGER NF, NF2, N1, N2, NW, I1, I2, I3, N, CORR, IMOD
       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 )
*
* ...Colour factors and an abbreviation
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
       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
*
      PTqg2 =
     &  + nf*ca**2 * ( 2057.D0/9.D0 - 672.D0*x + 53753.D0/81.D0*x**2 + 
     &    2092.D0/81.D0*dx - 196.D0*z3 - 2416.D0/3.D0*z3*x + 328.D0/3.D0
     &    *z3*x**2 - 310.D0/3.D0*z2 + 308.D0*z2*x - 536.D0/3.D0*z2*x**2
     &     + 544.D0/9.D0*z2*dx + 34.D0*z2**2 + 252.D0/5.D0*z2**2*x + 
     &    704.D0/5.D0*z2**2*x**2 + 136.D0*Hr1(-1)*z3 + 272.D0*Hr1(-1)*
     &    z3*x + 272.D0*Hr1(-1)*z3*x**2 + 512.D0/3.D0*Hr1(-1)*z2 + 760.D
     &    0/3.D0*Hr1(-1)*z2*x - 24.D0*Hr1(-1)*z2*x**2 + 32.D0/3.D0*Hr1(
     &    -1)*z2*dx - 2996.D0/27.D0*Hr1(0) - 4832.D0/27.D0*Hr1(0)*x + 
     &    5044.D0/9.D0*Hr1(0)*x**2 + 40.D0*Hr1(0)*dx + 16.D0*Hr1(0)*z3
     &     + 960.D0*Hr1(0)*z3*x + 96.D0*Hr1(0)*z2 - 484.D0/3.D0*Hr1(0)*
     &    z2*x + 824.D0/3.D0*Hr1(0)*z2*x**2 + 64.D0/3.D0*Hr1(0)*z2*dx
     &     - 3940.D0/27.D0*Hr1(1) - 8416.D0/27.D0*Hr1(1)*x + 7448.D0/27.
     &    D0*Hr1(1)*x**2 + 752.D0/27.D0*Hr1(1)*dx + 248.D0*Hr1(1)*z3 - 
     &    496.D0*Hr1(1)*z3*x + 496.D0*Hr1(1)*z3*x**2 + 440.D0/3.D0*Hr1(
     &    1)*z2 - 664.D0/3.D0*Hr1(1)*z2*x + 464.D0/3.D0*Hr1(1)*z2*x**2
     &     - 64.D0/3.D0*Hr1(1)*z2*dx )
      PTqg2 = PTqg2 + nf*ca**2 * (  - 208.D0*Hr2(-1,-1)*z2 - 416.D0*
     &    Hr2(-1,-1)*z2*x - 416.D0*Hr2(-1,-1)*z2*x**2 - 152.D0/3.D0*
     &    Hr2(-1,0) - 256.D0/3.D0*Hr2(-1,0)*x - 248.D0*Hr2(-1,0)*x**2
     &     - 32.D0/3.D0*Hr2(-1,0)*dx + 112.D0*Hr2(-1,0)*z2 + 224.D0*
     &    Hr2(-1,0)*z2*x + 224.D0*Hr2(-1,0)*z2*x**2 + 72.D0*Hr2(0,-1)*
     &    z2 - 176.D0*Hr2(0,-1)*z2*x + 96.D0*Hr2(0,-1)*z2*x**2 - 598.D0/
     &    9.D0*Hr2(0,0) - 3508.D0/9.D0*Hr2(0,0)*x - 7208.D0/9.D0*Hr2(0,
     &    0)*x**2 - 128.D0/9.D0*Hr2(0,0)*dx - 28.D0*Hr2(0,0)*z2 + 632.D0
     &    *Hr2(0,0)*z2*x - 96.D0*Hr2(0,0)*z2*x**2 - 1262.D0/9.D0*Hr2(0,
     &    1) + 4828.D0/9.D0*Hr2(0,1)*x - 2900.D0/3.D0*Hr2(0,1)*x**2 + 
     &    112.D0/9.D0*Hr2(0,1)*dx + 168.D0*Hr2(0,1)*z2 - 432.D0*Hr2(0,1
     &    )*z2*x + 384.D0*Hr2(0,1)*z2*x**2 - 562.D0/9.D0*Hr2(1,0) + 92.D
     &    0/9.D0*Hr2(1,0)*x + 704.D0/3.D0*Hr2(1,0)*x**2 - 8.D0/9.D0*
     &    Hr2(1,0)*dx + 48.D0*Hr2(1,0)*z2 - 96.D0*Hr2(1,0)*z2*x + 96.D0
     &    *Hr2(1,0)*z2*x**2 + 134.D0/9.D0*Hr2(1,1) - 700.D0/9.D0*Hr2(1,
     &    1)*x )
      PTqg2 = PTqg2 + nf*ca**2 * ( 652.D0/3.D0*Hr2(1,1)*x**2 - 128.D0/9.
     &    D0*Hr2(1,1)*dx - 16.D0*Hr2(1,1)*z2 + 32.D0*Hr2(1,1)*z2*x - 32.
     &    D0*Hr2(1,1)*z2*x**2 + 320.D0/3.D0*Hr3(-1,-1,0) + 112.D0/3.D0*
     &    Hr3(-1,-1,0)*x - 208.D0*Hr3(-1,-1,0)*x**2 - 64.D0/3.D0*Hr3(-1
     &    ,-1,0)*dx - 872.D0/3.D0*Hr3(-1,0,0) - 520.D0/3.D0*Hr3(-1,0,0)
     &    *x + 304.D0/3.D0*Hr3(-1,0,0)*x**2 + 128.D0/3.D0*Hr3(-1,0,0)*
     &    dx - 352.D0/3.D0*Hr3(-1,0,1) - 704.D0/3.D0*Hr3(-1,0,1)*x - 80.
     &    D0*Hr3(-1,0,1)*x**2 - 64.D0/3.D0*Hr3(-1,0,1)*dx + 184.D0/3.D0
     &    *Hr3(0,-1,0) - 344.D0/3.D0*Hr3(0,-1,0)*x + 832.D0/3.D0*Hr3(0,
     &    -1,0)*x**2 + 64.D0/3.D0*Hr3(0,-1,0)*dx + 64.D0/3.D0*Hr3(0,0,0
     &    ) + 2680.D0/3.D0*Hr3(0,0,0)*x - 192.D0*Hr3(0,0,0)*x**2 - 128.D
     &    0/3.D0*Hr3(0,0,0)*dx + 64.D0/3.D0*Hr3(0,0,1) + 292.D0*Hr3(0,0
     &    ,1)*x + 120.D0*Hr3(0,0,1)*x**2 - 64.D0/3.D0*Hr3(0,0,1)*dx - 
     &    100.D0*Hr3(0,1,0) - 312.D0*Hr3(0,1,0)*x + 72.D0*Hr3(0,1,0)*
     &    x**2 + 64.D0/3.D0*Hr3(0,1,0)*dx + 128.D0*Hr3(0,1,1) - 400.D0*
     &    Hr3(0,1,1)*x )
      PTqg2 = PTqg2 + nf*ca**2 * ( 640.D0/3.D0*Hr3(0,1,1)*x**2 + 32.D0/
     &    3.D0*Hr3(0,1,1)*dx + 728.D0/3.D0*Hr3(1,0,0) + 188.D0/3.D0*
     &    Hr3(1,0,0)*x - 304.D0/3.D0*Hr3(1,0,0)*x**2 + 224.D0/3.D0*Hr3(
     &    1,0,0)*dx + 160.D0/3.D0*Hr3(1,0,1) + 304.D0/3.D0*Hr3(1,0,1)*x
     &     - 472.D0/3.D0*Hr3(1,0,1)*x**2 + 32.D0*Hr3(1,0,1)*dx + 160.D0/
     &    3.D0*Hr3(1,1,0) + 304.D0/3.D0*Hr3(1,1,0)*x - 472.D0/3.D0*Hr3(
     &    1,1,0)*x**2 + 32.D0*Hr3(1,1,0)*dx - 196.D0/3.D0*Hr3(1,1,1) + 
     &    344.D0/3.D0*Hr3(1,1,1)*x - 400.D0/3.D0*Hr3(1,1,1)*x**2 + 32.D0
     &    /3.D0*Hr3(1,1,1)*dx - 96.D0*Hr4(-1,-1,-1,0) - 192.D0*Hr4(-1,
     &    -1,-1,0)*x - 192.D0*Hr4(-1,-1,-1,0)*x**2 + 64.D0*Hr4(-1,-1,0,
     &    0) + 128.D0*Hr4(-1,-1,0,0)*x + 128.D0*Hr4(-1,-1,0,0)*x**2 + 
     &    160.D0*Hr4(-1,-1,0,1) + 320.D0*Hr4(-1,-1,0,1)*x + 320.D0*Hr4(
     &    -1,-1,0,1)*x**2 + 112.D0*Hr4(-1,0,-1,0) + 224.D0*Hr4(-1,0,-1,
     &    0)*x + 224.D0*Hr4(-1,0,-1,0)*x**2 + 8.D0*Hr4(-1,0,0,0) + 16.D0
     &    *Hr4(-1,0,0,0)*x + 16.D0*Hr4(-1,0,0,0)*x**2 - 64.D0*Hr4(-1,0,
     &    1,0) )
      PTqg2 = PTqg2 + nf*ca**2 * (  - 128.D0*Hr4(-1,0,1,0)*x - 128.D0*
     &    Hr4(-1,0,1,0)*x**2 - 32.D0*Hr4(-1,0,1,1) - 64.D0*Hr4(-1,0,1,1
     &    )*x - 64.D0*Hr4(-1,0,1,1)*x**2 + 16.D0*Hr4(0,-1,-1,0) - 352.D0
     &    *Hr4(0,-1,-1,0)*x - 64.D0*Hr4(0,-1,-1,0)*x**2 - 104.D0*Hr4(0,
     &    -1,0,0) + 432.D0*Hr4(0,-1,0,0)*x - 32.D0*Hr4(0,-1,0,0)*x**2
     &     - 64.D0*Hr4(0,-1,0,1) - 128.D0*Hr4(0,-1,0,1)*x**2 - 56.D0*
     &    Hr4(0,0,-1,0) + 496.D0*Hr4(0,0,-1,0)*x + 64.D0*Hr4(0,0,0,0)
     &     - 1264.D0*Hr4(0,0,0,0)*x - 100.D0*Hr4(0,0,0,1) - 776.D0*Hr4(
     &    0,0,0,1)*x + 96.D0*Hr4(0,0,0,1)*x**2 + 64.D0*Hr4(0,0,1,0) - 
     &    96.D0*Hr4(0,0,1,0)*x + 160.D0*Hr4(0,0,1,0)*x**2 + 176.D0*Hr4(
     &    0,0,1,1) - 192.D0*Hr4(0,0,1,1)*x + 320.D0*Hr4(0,0,1,1)*x**2
     &     + 92.D0*Hr4(0,1,0,0) + 408.D0*Hr4(0,1,0,0)*x - 32.D0*Hr4(0,1
     &    ,0,0)*x**2 + 288.D0*Hr4(0,1,0,1)*x - 96.D0*Hr4(0,1,0,1)*x**2
     &     + 288.D0*Hr4(0,1,1,0)*x - 96.D0*Hr4(0,1,1,0)*x**2 - 80.D0*
     &    Hr4(0,1,1,1) + 256.D0*Hr4(0,1,1,1)*x - 192.D0*Hr4(0,1,1,1)*
     &    x**2 )
      PTqg2 = PTqg2 + nf*ca**2 * (  - 16.D0*Hr4(1,0,-1,0) + 32.D0*Hr4(1
     &    ,0,-1,0)*x - 32.D0*Hr4(1,0,-1,0)*x**2 - 120.D0*Hr4(1,0,0,0)
     &     + 240.D0*Hr4(1,0,0,0)*x - 240.D0*Hr4(1,0,0,0)*x**2 - 128.D0*
     &    Hr4(1,0,0,1) + 256.D0*Hr4(1,0,0,1)*x - 256.D0*Hr4(1,0,0,1)*
     &    x**2 - 32.D0*Hr4(1,0,1,1) + 64.D0*Hr4(1,0,1,1)*x - 64.D0*Hr4(
     &    1,0,1,1)*x**2 - 160.D0*Hr4(1,1,0,0) + 320.D0*Hr4(1,1,0,0)*x
     &     - 320.D0*Hr4(1,1,0,0)*x**2 - 32.D0*Hr4(1,1,0,1) + 64.D0*Hr4(
     &    1,1,0,1)*x - 64.D0*Hr4(1,1,0,1)*x**2 - 32.D0*Hr4(1,1,1,0) + 
     &    64.D0*Hr4(1,1,1,0)*x - 64.D0*Hr4(1,1,1,0)*x**2 + 32.D0*Hr4(1,
     &    1,1,1) - 64.D0*Hr4(1,1,1,1)*x + 64.D0*Hr4(1,1,1,1)*x**2 )
      PTqg2 = PTqg2 + nf*cf*ca * ( 17597.D0/36.D0 + 2659.D0/6.D0*x - 
     &    3092.D0/3.D0*x**2 - 220.D0/3.D0*dx + 1340.D0/3.D0*z3 + 4904.D0
     &    /3.D0*z3*x + 488.D0*z3*x**2 + 64.D0*z3*dx + 28.D0/9.D0*z2 + 
     &    4540.D0/9.D0*z2*x - 2620.D0/3.D0*z2*x**2 - 44.D0/5.D0*z2**2
     &     + 2632.D0/5.D0*z2**2*x - 128.D0*z2**2*x**2 - 40.D0*Hr1(-1)*
     &    z3 - 80.D0*Hr1(-1)*z3*x - 80.D0*Hr1(-1)*z3*x**2 - 96.D0*Hr1(
     &    -1)*z2 - 160.D0*Hr1(-1)*z2*x + 32.D0*Hr1(-1)*z2*x**2 - 2422.D0
     &    /27.D0*Hr1(0) + 30203.D0/27.D0*Hr1(0)*x - 17918.D0/27.D0*Hr1(
     &    0)*x**2 - 56.D0*Hr1(0)*z3 + 16.D0*Hr1(0)*z3*x - 160.D0*Hr1(0)
     &    *z3*x**2 - 260.D0/3.D0*Hr1(0)*z2 - 176.D0/3.D0*Hr1(0)*z2*x - 
     &    8.D0*Hr1(0)*z2*x**2 + 2788.D0/27.D0*Hr1(1) + 19660.D0/27.D0*
     &    Hr1(1)*x - 19970.D0/27.D0*Hr1(1)*x**2 - 80.D0/27.D0*Hr1(1)*dx
     &     - 792.D0*Hr1(1)*z3 + 1584.D0*Hr1(1)*z3*x - 1584.D0*Hr1(1)*z3
     &    *x**2 - 188.D0*Hr1(1)*z2 + 272.D0*Hr1(1)*z2*x - 120.D0*Hr1(1)
     &    *z2*x**2 + 32.D0*Hr1(1)*z2*dx + 80.D0*Hr2(-1,-1)*z2 + 160.D0*
     &    Hr2(-1,-1)*z2*x )
      PTqg2 = PTqg2 + nf*cf*ca * ( 160.D0*Hr2(-1,-1)*z2*x**2 + 672.D0*
     &    Hr2(-1,0) + 544.D0*Hr2(-1,0)*x - 16.D0*Hr2(-1,0)*x**2 - 64.D0
     &    *Hr2(0,-1)*z2 - 32.D0*Hr2(0,-1)*z2*x - 128.D0*Hr2(0,-1)*z2*
     &    x**2 - 1007.D0/9.D0*Hr2(0,0) - 6254.D0/9.D0*Hr2(0,0)*x + 2444.
     &    D0/9.D0*Hr2(0,0)*x**2 + 32.D0*Hr2(0,0)*z2 - 256.D0*Hr2(0,0)*
     &    z2*x + 256.D0*Hr2(0,0)*z2*x**2 + 2204.D0/9.D0*Hr2(0,1) - 3772.
     &    D0/9.D0*Hr2(0,1)*x + 3176.D0/3.D0*Hr2(0,1)*x**2 - 200.D0*Hr2(
     &    0,1)*z2 + 592.D0*Hr2(0,1)*z2*x - 512.D0*Hr2(0,1)*z2*x**2 - 
     &    212.D0*Hr2(1,0) + 2036.D0/3.D0*Hr2(1,0)*x - 3284.D0/3.D0*Hr2(
     &    1,0)*x**2 + 304.D0/3.D0*Hr2(1,0)*dx + 80.D0*Hr2(1,0)*z2 - 160.
     &    D0*Hr2(1,0)*z2*x + 160.D0*Hr2(1,0)*z2*x**2 - 488.D0/9.D0*Hr2(
     &    1,1) + 868.D0/9.D0*Hr2(1,1)*x - 1904.D0/9.D0*Hr2(1,1)*x**2 + 
     &    304.D0/9.D0*Hr2(1,1)*dx + 32.D0*Hr2(1,1)*z2 - 64.D0*Hr2(1,1)*
     &    z2*x + 64.D0*Hr2(1,1)*z2*x**2 - 288.D0*Hr3(-1,-1,0) - 448.D0*
     &    Hr3(-1,-1,0)*x - 64.D0*Hr3(-1,-1,0)*x**2 + 16.D0*Hr3(-1,0,0)
     &     - 80.D0*Hr3(-1,0,0)*x )
      PTqg2 = PTqg2 + nf*cf*ca * (  - 144.D0*Hr3(-1,0,0)*x**2 - 48.D0*
     &    Hr3(-1,0,1) - 64.D0*Hr3(-1,0,1)*x - 64.D0*Hr3(-1,0,1)*x**2 + 
     &    336.D0*Hr3(0,-1,0) + 16.D0*Hr3(0,-1,0)*x + 160.D0*Hr3(0,-1,0)
     &    *x**2 - 148.D0/3.D0*Hr3(0,0,0) + 1832.D0/3.D0*Hr3(0,0,0)*x + 
     &    640.D0/3.D0*Hr3(0,0,0)*x**2 + 388.D0/3.D0*Hr3(0,0,1) + 256.D0/
     &    3.D0*Hr3(0,0,1)*x - 40.D0/3.D0*Hr3(0,0,1)*x**2 + 28.D0/3.D0*
     &    Hr3(0,1,0) + 2032.D0/3.D0*Hr3(0,1,0)*x + 280.D0/3.D0*Hr3(0,1,
     &    0)*x**2 - 592.D0/3.D0*Hr3(0,1,1) + 1352.D0/3.D0*Hr3(0,1,1)*x
     &     - 192.D0*Hr3(0,1,1)*x**2 - 16.D0*Hr3(1,0,0) - 504.D0*Hr3(1,0
     &    ,0)*x + 368.D0*Hr3(1,0,0)*x**2 - 96.D0*Hr3(1,0,0)*dx + 268.D0/
     &    3.D0*Hr3(1,0,1) - 608.D0/3.D0*Hr3(1,0,1)*x + 728.D0/3.D0*Hr3(
     &    1,0,1)*x**2 - 32.D0*Hr3(1,0,1)*dx + 452.D0/3.D0*Hr3(1,1,0) - 
     &    1024.D0/3.D0*Hr3(1,1,0)*x + 1048.D0/3.D0*Hr3(1,1,0)*x**2 - 32.
     &    D0*Hr3(1,1,0)*dx + 424.D0/3.D0*Hr3(1,1,1) - 536.D0/3.D0*Hr3(1
     &    ,1,1)*x + 544.D0/3.D0*Hr3(1,1,1)*x**2 - 32.D0/3.D0*Hr3(1,1,1)
     &    *dx )
      PTqg2 = PTqg2 + nf*cf*ca * ( 96.D0*Hr4(-1,-1,-1,0) + 192.D0*Hr4(
     &    -1,-1,-1,0)*x + 192.D0*Hr4(-1,-1,-1,0)*x**2 + 64.D0*Hr4(-1,-1
     &    ,0,0) + 128.D0*Hr4(-1,-1,0,0)*x + 128.D0*Hr4(-1,-1,0,0)*x**2
     &     - 32.D0*Hr4(-1,-1,0,1) - 64.D0*Hr4(-1,-1,0,1)*x - 64.D0*Hr4(
     &    -1,-1,0,1)*x**2 - 80.D0*Hr4(-1,0,-1,0) - 160.D0*Hr4(-1,0,-1,0
     &    )*x - 160.D0*Hr4(-1,0,-1,0)*x**2 - 24.D0*Hr4(-1,0,0,0) - 48.D0
     &    *Hr4(-1,0,0,0)*x - 48.D0*Hr4(-1,0,0,0)*x**2 + 96.D0*Hr4(-1,0,
     &    1,0) + 192.D0*Hr4(-1,0,1,0)*x + 192.D0*Hr4(-1,0,1,0)*x**2 + 
     &    32.D0*Hr4(-1,0,1,1) + 64.D0*Hr4(-1,0,1,1)*x + 64.D0*Hr4(-1,0,
     &    1,1)*x**2 - 128.D0*Hr4(0,-1,-1,0) - 64.D0*Hr4(0,-1,-1,0)*x - 
     &    256.D0*Hr4(0,-1,-1,0)*x**2 - 96.D0*Hr4(0,-1,0,0)*x - 128.D0*
     &    Hr4(0,-1,0,0)*x**2 + 144.D0*Hr4(0,0,-1,0) + 32.D0*Hr4(0,0,-1,
     &    0)*x + 384.D0*Hr4(0,0,-1,0)*x**2 - 32.D0*Hr4(0,0,0,0) - 64.D0
     &    *Hr4(0,0,0,0)*x + 64.D0*Hr4(0,0,0,1) + 96.D0*Hr4(0,0,0,1)*x
     &     + 128.D0*Hr4(0,0,0,1)*x**2 - 160.D0*Hr4(0,0,1,0) - 160.D0*
     &    Hr4(0,0,1,0)*x )
      PTqg2 = PTqg2 + nf*cf*ca * (  - 320.D0*Hr4(0,0,1,0)*x**2 - 176.D0
     &    *Hr4(0,0,1,1) + 192.D0*Hr4(0,0,1,1)*x - 384.D0*Hr4(0,0,1,1)*
     &    x**2 + 136.D0*Hr4(0,1,0,0) - 1040.D0*Hr4(0,1,0,0)*x + 672.D0*
     &    Hr4(0,1,0,0)*x**2 + 168.D0*Hr4(0,1,0,1) - 624.D0*Hr4(0,1,0,1)
     &    *x + 448.D0*Hr4(0,1,0,1)*x**2 + 184.D0*Hr4(0,1,1,0) - 656.D0*
     &    Hr4(0,1,1,0)*x + 512.D0*Hr4(0,1,1,0)*x**2 + 136.D0*Hr4(0,1,1,
     &    1) - 368.D0*Hr4(0,1,1,1)*x + 320.D0*Hr4(0,1,1,1)*x**2 + 80.D0
     &    *Hr4(1,0,-1,0) - 160.D0*Hr4(1,0,-1,0)*x + 160.D0*Hr4(1,0,-1,0
     &    )*x**2 - 56.D0*Hr4(1,0,0,0) + 112.D0*Hr4(1,0,0,0)*x - 112.D0*
     &    Hr4(1,0,0,0)*x**2 + 112.D0*Hr4(1,0,0,1) - 224.D0*Hr4(1,0,0,1)
     &    *x + 224.D0*Hr4(1,0,0,1)*x**2 - 240.D0*Hr4(1,0,1,0) + 480.D0*
     &    Hr4(1,0,1,0)*x - 480.D0*Hr4(1,0,1,0)*x**2 - 32.D0*Hr4(1,0,1,1
     &    ) + 64.D0*Hr4(1,0,1,1)*x - 64.D0*Hr4(1,0,1,1)*x**2 + 208.D0*
     &    Hr4(1,1,0,0) - 416.D0*Hr4(1,1,0,0)*x + 416.D0*Hr4(1,1,0,0)*
     &    x**2 - 48.D0*Hr4(1,1,0,1) + 96.D0*Hr4(1,1,0,1)*x - 96.D0*Hr4(
     &    1,1,0,1)*x**2 )
      PTqg2 = PTqg2 + nf*cf*ca * (  - 16.D0*Hr4(1,1,1,0) + 32.D0*Hr4(1,
     &    1,1,0)*x - 32.D0*Hr4(1,1,1,0)*x**2 - 64.D0*Hr4(1,1,1,1) + 128.
     &    D0*Hr4(1,1,1,1)*x - 128.D0*Hr4(1,1,1,1)*x**2 )
      PTqg2 = PTqg2 + nf*cf**2 * (  - 227.D0/4.D0 - 1967.D0/2.D0*x + 
     &    1069.D0*x**2 + 156.D0*z3 - 816.D0*z3*x - 640.D0*z3*x**2 + 106.
     &    D0*z2 - 564.D0*z2*x + 564.D0*z2*x**2 + 56.D0*z2**2 - 1904.D0/
     &    5.D0*z2**2*x + 352.D0/5.D0*z2**2*x**2 + 128.D0*Hr1(-1)*z2 + 
     &    160.D0*Hr1(-1)*z2*x + 32.D0*Hr1(-1)*z2*x**2 + 92.D0*Hr1(0) - 
     &    1143.D0*Hr1(0)*x + 110.D0*Hr1(0)*x**2 + 40.D0*Hr1(0)*z3 - 304.
     &    D0*Hr1(0)*z3*x + 160.D0*Hr1(0)*z3*x**2 + 12.D0*Hr1(0)*z2 + 
     &    272.D0*Hr1(0)*z2*x - 176.D0*Hr1(0)*z2*x**2 + 180.D0*Hr1(1) - 
     &    480.D0*Hr1(1)*x + 318.D0*Hr1(1)*x**2 + 496.D0*Hr1(1)*z3 - 992.
     &    D0*Hr1(1)*z3*x + 992.D0*Hr1(1)*z3*x**2 + 108.D0*Hr1(1)*z2 - 
     &    64.D0*Hr1(1)*z2*x - 80.D0*Hr1(1)*z2*x**2 - 496.D0*Hr2(-1,0)
     &     - 240.D0*Hr2(-1,0)*x + 224.D0*Hr2(-1,0)*x**2 - 96.D0*Hr2(-1,
     &    0)*z2 - 192.D0*Hr2(-1,0)*z2*x - 192.D0*Hr2(-1,0)*z2*x**2 + 64.
     &    D0*Hr2(0,-1)*z2 + 64.D0*Hr2(0,-1)*z2*x + 128.D0*Hr2(0,-1)*z2*
     &    x**2 + 107.D0*Hr2(0,0) - 42.D0*Hr2(0,0)*x - 76.D0*Hr2(0,0)*
     &    x**2 )
      PTqg2 = PTqg2 + nf*cf**2 * (  - 40.D0*Hr2(0,0)*z2 + 144.D0*Hr2(0,
     &    0)*z2*x - 160.D0*Hr2(0,0)*z2*x**2 - 2.D0*Hr2(0,1) + 116.D0*
     &    Hr2(0,1)*x - 324.D0*Hr2(0,1)*x**2 + 104.D0*Hr2(0,1)*z2 - 144.D
     &    0*Hr2(0,1)*z2*x + 224.D0*Hr2(0,1)*z2*x**2 + 294.D0*Hr2(1,0)
     &     - 560.D0*Hr2(1,0)*x + 468.D0*Hr2(1,0)*x**2 - 112.D0*Hr2(1,0)
     &    *z2 + 224.D0*Hr2(1,0)*z2*x - 224.D0*Hr2(1,0)*z2*x**2 + 34.D0*
     &    Hr2(1,1) - 32.D0*Hr2(1,1)*x - 4.D0*Hr2(1,1)*x**2 - 16.D0*Hr2(
     &    1,1)*z2 + 32.D0*Hr2(1,1)*z2*x - 32.D0*Hr2(1,1)*z2*x**2 + 256.D
     &    0*Hr3(-1,-1,0) + 320.D0*Hr3(-1,-1,0)*x + 64.D0*Hr3(-1,-1,0)*
     &    x**2 + 64.D0*Hr3(-1,0,0) + 320.D0*Hr3(-1,0,0)*x + 256.D0*Hr3(
     &    -1,0,0)*x**2 - 240.D0*Hr3(0,-1,0) + 64.D0*Hr3(0,-1,0)*x - 256.
     &    D0*Hr3(0,-1,0)*x**2 + 36.D0*Hr3(0,0,0) - 584.D0*Hr3(0,0,0)*x
     &     - 192.D0*Hr3(0,0,0)*x**2 - 28.D0*Hr3(0,0,1) - 144.D0*Hr3(0,0
     &    ,1)*x + 16.D0*Hr3(0,0,1)*x**2 + 184.D0*Hr3(0,1,0) - 128.D0*
     &    Hr3(0,1,0)*x - 112.D0*Hr3(0,1,0)*x**2 + 104.D0*Hr3(0,1,1) - 
     &    80.D0*Hr3(0,1,1)*x )
      PTqg2 = PTqg2 + nf*cf**2 * ( 16.D0*Hr3(0,1,1)*x**2 - 288.D0*Hr3(1
     &    ,0,0) + 352.D0*Hr3(1,0,0)*x - 160.D0*Hr3(1,0,0)*x**2 - 124.D0
     &    *Hr3(1,0,1) + 64.D0*Hr3(1,0,1)*x - 48.D0*Hr3(1,0,1)*x**2 - 
     &    180.D0*Hr3(1,1,0) + 192.D0*Hr3(1,1,0)*x - 144.D0*Hr3(1,1,0)*
     &    x**2 - 76.D0*Hr3(1,1,1) + 64.D0*Hr3(1,1,1)*x - 48.D0*Hr3(1,1,
     &    1)*x**2 - 192.D0*Hr4(-1,-1,0,0) - 384.D0*Hr4(-1,-1,0,0)*x - 
     &    384.D0*Hr4(-1,-1,0,0)*x**2 - 96.D0*Hr4(-1,0,-1,0) - 192.D0*
     &    Hr4(-1,0,-1,0)*x - 192.D0*Hr4(-1,0,-1,0)*x**2 + 176.D0*Hr4(-1
     &    ,0,0,0) + 352.D0*Hr4(-1,0,0,0)*x + 352.D0*Hr4(-1,0,0,0)*x**2
     &     + 128.D0*Hr4(0,-1,-1,0) + 128.D0*Hr4(0,-1,-1,0)*x + 256.D0*
     &    Hr4(0,-1,-1,0)*x**2 + 32.D0*Hr4(0,-1,0,0) + 128.D0*Hr4(0,-1,0
     &    ,0)*x + 256.D0*Hr4(0,-1,0,0)*x**2 - 96.D0*Hr4(0,0,-1,0) + 64.D
     &    0*Hr4(0,0,-1,0)*x - 256.D0*Hr4(0,0,-1,0)*x**2 + 16.D0*Hr4(0,0
     &    ,0,0) - 256.D0*Hr4(0,0,0,0)*x + 128.D0*Hr4(0,0,0,0)*x**2 + 8.D
     &    0*Hr4(0,0,0,1) - 16.D0*Hr4(0,0,0,1)*x + 32.D0*Hr4(0,0,0,1)*
     &    x**2 )
      PTqg2 = PTqg2 + nf*cf**2 * ( 96.D0*Hr4(0,0,1,0) - 192.D0*Hr4(0,0,
     &    1,0)*x + 352.D0*Hr4(0,0,1,0)*x**2 + 64.D0*Hr4(0,0,1,1) - 128.D
     &    0*Hr4(0,0,1,1)*x + 192.D0*Hr4(0,0,1,1)*x**2 - 176.D0*Hr4(0,1,
     &    0,0) + 288.D0*Hr4(0,1,0,0)*x - 416.D0*Hr4(0,1,0,0)*x**2 - 104.
     &    D0*Hr4(0,1,0,1) + 208.D0*Hr4(0,1,0,1)*x - 224.D0*Hr4(0,1,0,1)
     &    *x**2 - 120.D0*Hr4(0,1,1,0) + 240.D0*Hr4(0,1,1,0)*x - 288.D0*
     &    Hr4(0,1,1,0)*x**2 - 56.D0*Hr4(0,1,1,1) + 112.D0*Hr4(0,1,1,1)*
     &    x - 128.D0*Hr4(0,1,1,1)*x**2 - 96.D0*Hr4(1,0,-1,0) + 192.D0*
     &    Hr4(1,0,-1,0)*x - 192.D0*Hr4(1,0,-1,0)*x**2 + 336.D0*Hr4(1,0,
     &    0,0) - 672.D0*Hr4(1,0,0,0)*x + 672.D0*Hr4(1,0,0,0)*x**2 + 80.D
     &    0*Hr4(1,0,0,1) - 160.D0*Hr4(1,0,0,1)*x + 160.D0*Hr4(1,0,0,1)*
     &    x**2 + 304.D0*Hr4(1,0,1,0) - 608.D0*Hr4(1,0,1,0)*x + 608.D0*
     &    Hr4(1,0,1,0)*x**2 + 64.D0*Hr4(1,0,1,1) - 128.D0*Hr4(1,0,1,1)*
     &    x + 128.D0*Hr4(1,0,1,1)*x**2 + 16.D0*Hr4(1,1,0,0) - 32.D0*
     &    Hr4(1,1,0,0)*x + 32.D0*Hr4(1,1,0,0)*x**2 + 80.D0*Hr4(1,1,0,1)
     &     - 160.D0*Hr4(1,1,0,1)*x )
      PTqg2 = PTqg2 + nf*cf**2 * ( 160.D0*Hr4(1,1,0,1)*x**2 + 48.D0*
     &    Hr4(1,1,1,0) - 96.D0*Hr4(1,1,1,0)*x + 96.D0*Hr4(1,1,1,0)*x**2
     &     + 32.D0*Hr4(1,1,1,1) - 64.D0*Hr4(1,1,1,1)*x + 64.D0*Hr4(1,1,
     &    1,1)*x**2 )
      PTqg2 = PTqg2 + nf2*ca * (  - 14.D0/9.D0 - 1216.D0/9.D0*x + 916.
     &    D0/9.D0*x**2 - 44.D0/9.D0*dx - 40.D0*z3 + 352.D0/3.D0*z3*x - 
     &    272.D0/3.D0*z3*x**2 + 100.D0/3.D0*z2 + 208.D0/3.D0*z2*x - 16.D
     &    0/9.D0*z2*x**2 + 32.D0/9.D0*z2*dx + 64.D0/3.D0*Hr1(-1)*z2 + 
     &    128.D0/3.D0*Hr1(-1)*z2*x + 128.D0/3.D0*Hr1(-1)*z2*x**2 + 44.D0
     &    /27.D0*Hr1(0) + 752.D0/27.D0*Hr1(0)*x - 2212.D0/27.D0*Hr1(0)*
     &    x**2 - 256.D0/27.D0*Hr1(0)*dx - 24.D0*Hr1(0)*z2 + 112.D0/3.D0
     &    *Hr1(0)*z2*x - 112.D0/3.D0*Hr1(0)*z2*x**2 + 592.D0/27.D0*Hr1(
     &    1) - 1232.D0/27.D0*Hr1(1)*x + 2632.D0/27.D0*Hr1(1)*x**2 - 320.
     &    D0/27.D0*Hr1(1)*dx - 32.D0/3.D0*Hr1(1)*z2 + 64.D0/3.D0*Hr1(1)
     &    *z2*x - 64.D0/3.D0*Hr1(1)*z2*x**2 - 16.D0*Hr2(-1,0) + 32.D0/3.
     &    D0*Hr2(-1,0)*x + 512.D0/9.D0*Hr2(-1,0)*x**2 + 32.D0/9.D0*Hr2(
     &    -1,0)*dx - 296.D0/9.D0*Hr2(0,0) - 416.D0/9.D0*Hr2(0,0)*x - 
     &    536.D0/9.D0*Hr2(0,0)*x**2 - 64.D0/9.D0*Hr2(0,0)*dx + 356.D0/9.
     &    D0*Hr2(0,1) - 1360.D0/9.D0*Hr2(0,1)*x + 632.D0/9.D0*Hr2(0,1)*
     &    x**2 )
      PTqg2 = PTqg2 + nf2*ca * ( 16.D0/9.D0*Hr2(1,0) + 496.D0/9.D0*
     &    Hr2(1,0)*x - 664.D0/9.D0*Hr2(1,0)*x**2 + 32.D0/3.D0*Hr2(1,0)*
     &    dx - 320.D0/9.D0*Hr2(1,1) + 544.D0/9.D0*Hr2(1,1)*x - 200.D0/3.
     &    D0*Hr2(1,1)*x**2 + 32.D0/9.D0*Hr2(1,1)*dx + 64.D0/3.D0*Hr3(-1
     &    ,-1,0) + 128.D0/3.D0*Hr3(-1,-1,0)*x + 128.D0/3.D0*Hr3(-1,-1,0
     &    )*x**2 + 32.D0/3.D0*Hr3(-1,0,0) + 64.D0/3.D0*Hr3(-1,0,0)*x + 
     &    64.D0/3.D0*Hr3(-1,0,0)*x**2 - 32.D0/3.D0*Hr3(-1,0,1) - 64.D0/
     &    3.D0*Hr3(-1,0,1)*x - 64.D0/3.D0*Hr3(-1,0,1)*x**2 - 64.D0/3.D0
     &    *Hr3(0,-1,0) + 128.D0/3.D0*Hr3(0,-1,0)*x - 64.D0/3.D0*Hr3(0,
     &    -1,0)*x**2 - 64.D0/3.D0*Hr3(0,0,0) - 640.D0/3.D0*Hr3(0,0,0)*x
     &     + 104.D0/3.D0*Hr3(0,0,1) - 16.D0*Hr3(0,0,1)*x + 176.D0/3.D0*
     &    Hr3(0,0,1)*x**2 + 8.D0*Hr3(0,1,0) + 80.D0*Hr3(0,1,0)*x - 16.D0
     &    *Hr3(0,1,0)*x**2 - 32.D0*Hr3(0,1,1) + 96.D0*Hr3(0,1,1)*x - 
     &    224.D0/3.D0*Hr3(0,1,1)*x**2 - 152.D0/3.D0*Hr3(1,0,0) + 304.D0/
     &    3.D0*Hr3(1,0,0)*x - 304.D0/3.D0*Hr3(1,0,0)*x**2 - 16.D0/3.D0*
     &    Hr3(1,0,1) )
      PTqg2 = PTqg2 + nf2*ca * ( 32.D0/3.D0*Hr3(1,0,1)*x - 32.D0/3.D0
     &    *Hr3(1,0,1)*x**2 - 16.D0/3.D0*Hr3(1,1,0) + 32.D0/3.D0*Hr3(1,1
     &    ,0)*x - 32.D0/3.D0*Hr3(1,1,0)*x**2 + 40.D0/3.D0*Hr3(1,1,1) - 
     &    80.D0/3.D0*Hr3(1,1,1)*x + 80.D0/3.D0*Hr3(1,1,1)*x**2 )
      PTqg2 = PTqg2 + nf2*cf * (  - 4847.D0/54.D0 - 2375.D0/27.D0*x
     &     + 4066.D0/27.D0*x**2 + 200.D0/27.D0*dx + 304.D0/3.D0*z3 - 
     &    608.D0/3.D0*z3*x + 160.D0*z3*x**2 + 272.D0/9.D0*z2 + 8.D0/9.D0
     &    *z2*x + 32.D0*z2*x**2 - 64.D0/9.D0*z2*dx + 1684.D0/27.D0*Hr1(
     &    0) - 4022.D0/27.D0*Hr1(0)*x - 4480.D0/27.D0*Hr1(0)*x**2 + 512.
     &    D0/27.D0*Hr1(0)*dx + 80.D0/3.D0*Hr1(0)*z2 - 352.D0/3.D0*Hr1(0
     &    )*z2*x + 48.D0*Hr1(0)*z2*x**2 - 808.D0/27.D0*Hr1(1) + 560.D0/
     &    27.D0*Hr1(1)*x - 248.D0/27.D0*Hr1(1)*x**2 - 8.D0*Hr1(1)*z2 + 
     &    16.D0*Hr1(1)*z2*x - 16.D0*Hr1(1)*z2*x**2 + 48.D0*Hr2(-1,0) - 
     &    496.D0/9.D0*Hr2(-1,0)*x**2 - 64.D0/9.D0*Hr2(-1,0)*dx + 902.D0/
     &    9.D0*Hr2(0,0) + 956.D0/9.D0*Hr2(0,0)*x + 1400.D0/9.D0*Hr2(0,0
     &    )*x**2 + 128.D0/9.D0*Hr2(0,0)*dx - 416.D0/9.D0*Hr2(0,1) + 376.
     &    D0/9.D0*Hr2(0,1)*x - 88.D0/3.D0*Hr2(0,1)*x**2 + 64.D0*Hr2(1,0
     &    ) - 248.D0/3.D0*Hr2(1,0)*x + 232.D0/3.D0*Hr2(1,0)*x**2 + 344.D
     &    0/9.D0*Hr2(1,1) - 376.D0/9.D0*Hr2(1,1)*x + 328.D0/9.D0*Hr2(1,
     &    1)*x**2 )
      PTqg2 = PTqg2 + nf2*cf * ( 32.D0*Hr3(0,-1,0) - 64.D0*Hr3(0,-1,0
     &    )*x + 64.D0/3.D0*Hr3(0,-1,0)*x**2 - 32.D0/3.D0*Hr3(0,0,0) + 
     &    112.D0/3.D0*Hr3(0,0,0)*x + 32.D0*Hr3(0,0,0)*x**2 - 64.D0/3.D0
     &    *Hr3(0,0,1) + 128.D0/3.D0*Hr3(0,0,1)*x - 176.D0/3.D0*Hr3(0,0,
     &    1)*x**2 + 80.D0/3.D0*Hr3(0,1,0) - 160.D0/3.D0*Hr3(0,1,0)*x + 
     &    176.D0/3.D0*Hr3(0,1,0)*x**2 + 64.D0/3.D0*Hr3(0,1,1) - 128.D0/
     &    3.D0*Hr3(0,1,1)*x + 48.D0*Hr3(0,1,1)*x**2 + 32.D0*Hr3(1,0,0)
     &     - 64.D0*Hr3(1,0,0)*x + 64.D0*Hr3(1,0,0)*x**2 - 40.D0/3.D0*
     &    Hr3(1,0,1) + 80.D0/3.D0*Hr3(1,0,1)*x - 80.D0/3.D0*Hr3(1,0,1)*
     &    x**2 - 56.D0/3.D0*Hr3(1,1,0) + 112.D0/3.D0*Hr3(1,1,0)*x - 112.
     &    D0/3.D0*Hr3(1,1,0)*x**2 - 40.D0/3.D0*Hr3(1,1,1) + 80.D0/3.D0*
     &    Hr3(1,1,1)*x - 80.D0/3.D0*Hr3(1,1,1)*x**2 - 32.D0*Hr4(0,0,0,0
     &    ) + 64.D0*Hr4(0,0,0,0)*x )
      PTqg2 = PTqg2 + nf**3 * (  - 8.D0/9.D0 + 16.D0/3.D0*x - 16.D0/3.D0
     &    *x**2 + 8.D0/3.D0*z2 - 16.D0/3.D0*z2*x + 16.D0/3.D0*z2*x**2
     &     + 40.D0/9.D0*Hr1(0) - 32.D0/9.D0*Hr1(0)*x + 32.D0/9.D0*Hr1(0
     &    )*x**2 - 40.D0/9.D0*Hr1(1) + 32.D0/9.D0*Hr1(1)*x - 32.D0/9.D0
     &    *Hr1(1)*x**2 + 8.D0/3.D0*Hr2(0,0) - 16.D0/3.D0*Hr2(0,0)*x + 
     &    16.D0/3.D0*Hr2(0,0)*x**2 - 8.D0/3.D0*Hr2(0,1) + 16.D0/3.D0*
     &    Hr2(0,1)*x - 16.D0/3.D0*Hr2(0,1)*x**2 - 8.D0/3.D0*Hr2(1,0) + 
     &    16.D0/3.D0*Hr2(1,0)*x - 16.D0/3.D0*Hr2(1,0)*x**2 + 8.D0/3.D0*
     &    Hr2(1,1) - 16.D0/3.D0*Hr2(1,1)*x + 16.D0/3.D0*Hr2(1,1)*x**2 )
      PTqg2 = PTqg2
     , + 2.*z2*(ca-cf)*b0* 2.*nf* (1.- 2.*x + 2.*x*x )
     ,     * ( 11.+24.*Hr1(0) ) * IMOD
*
       X2QGTA = PTqg2
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The gluon-quark splitting functions P_gq^(2)T
*
       FUNCTION X2GQTA (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 )
*
* ...Colour factors and an abbreviation
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
       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
*
      PTgq2 =
     &  + cf*ca**2 * (  - 17798.D0/27.D0 - 3436.D0/27.D0*x + 4136.D0/9.D
     &    0*x**2 + 3514.D0/9.D0*dx - 400.D0*z3 + 148.D0*z3*x - 640.D0/3.
     &    D0*z3*x**2 + 160.D0*z3*dx - 848.D0*z2 - 2800.D0/9.D0*z2*x - 
     &    352.D0/9.D0*z2*x**2 - 5176.D0/9.D0*z2*dx + 572.D0/5.D0*z2**2
     &     + 314.D0*z2**2*x + 1984.D0/5.D0*z2**2*dx - 16.D0*Hr1(-1)*z3
     &     - 8.D0*Hr1(-1)*z3*x - 16.D0*Hr1(-1)*z3*dx - 168.D0*Hr1(-1)*
     &    z2 - 160.D0*Hr1(-1)*z2*x - 32.D0/3.D0*Hr1(-1)*z2*x**2 - 176.D0
     &    /3.D0*Hr1(-1)*z2*dx - 32612.D0/27.D0*Hr1(0) - 23630.D0/27.D0*
     &    Hr1(0)*x - 12964.D0/27.D0*Hr1(0)*x**2 + 13384.D0/27.D0*Hr1(0)
     &    *dx + 352.D0*Hr1(0)*z3 + 240.D0*Hr1(0)*z3*x + 416.D0*Hr1(0)*
     &    z3*dx - 12.D0*Hr1(0)*z2 + 532.D0/3.D0*Hr1(0)*z2*x - 128.D0/3.D
     &    0*Hr1(0)*z2*x**2 - 752.D0/3.D0*Hr1(0)*z2*dx - 22270.D0/27.D0*
     &    Hr1(1) - 688.D0/27.D0*Hr1(1)*x - 1924.D0/9.D0*Hr1(1)*x**2 + 
     &    31504.D0/27.D0*Hr1(1)*dx - 656.D0*Hr1(1)*z3 + 328.D0*Hr1(1)*
     &    z3*x + 656.D0*Hr1(1)*z3*dx - 280.D0*Hr1(1)*z2 - 32.D0*Hr1(1)*
     &    z2*x )
      PTgq2 = PTgq2 + cf*ca**2 * (  - 64.D0*Hr1(1)*z2*x**2 + 344.D0*
     &    Hr1(1)*z2*dx - 32.D0*Hr2(-1,-1)*z2 - 16.D0*Hr2(-1,-1)*z2*x - 
     &    32.D0*Hr2(-1,-1)*z2*dx - 1280.D0/9.D0*Hr2(-1,0) + 1052.D0/9.D0
     &    *Hr2(-1,0)*x + 704.D0/9.D0*Hr2(-1,0)*x**2 - 248.D0/3.D0*Hr2(
     &    -1,0)*dx + 160.D0*Hr2(-1,0)*z2 + 80.D0*Hr2(-1,0)*z2*x + 160.D0
     &    *Hr2(-1,0)*z2*dx - 144.D0*Hr2(0,-1)*z2 - 8.D0*Hr2(0,-1)*z2*x
     &     - 64.D0*Hr2(0,-1)*z2*dx + 11848.D0/9.D0*Hr2(0,0) + 7280.D0/9.
     &    D0*Hr2(0,0)*x + 3688.D0/9.D0*Hr2(0,0)*x**2 + 3920.D0/9.D0*
     &    Hr2(0,0)*dx + 152.D0*Hr2(0,0)*z2 - 156.D0*Hr2(0,0)*z2*x - 128.
     &    D0*Hr2(0,0)*z2*dx + 6896.D0/9.D0*Hr2(0,1) + 4220.D0/9.D0*Hr2(
     &    0,1)*x + 1552.D0/9.D0*Hr2(0,1)*x**2 + 4664.D0/9.D0*Hr2(0,1)*
     &    dx + 48.D0*Hr2(0,1)*z2 + 248.D0*Hr2(0,1)*z2*x + 320.D0*Hr2(0,
     &    1)*z2*dx + 1160.D0/9.D0*Hr2(1,0) - 604.D0/9.D0*Hr2(1,0)*x + 
     &    760.D0/3.D0*Hr2(1,0)*x**2 - 4136.D0/9.D0*Hr2(1,0)*dx + 128.D0
     &    *Hr2(1,0)*z2 - 64.D0*Hr2(1,0)*z2*x - 128.D0*Hr2(1,0)*z2*dx + 
     &    808.D0/9.D0*Hr2(1,1) )
      PTgq2 = PTgq2 + cf*ca**2 * (  - 440.D0/9.D0*Hr2(1,1)*x + 560.D0/9.
     &    D0*Hr2(1,1)*x**2 - 532.D0/3.D0*Hr2(1,1)*dx - 96.D0*Hr2(1,1)*
     &    z2 + 48.D0*Hr2(1,1)*z2*x + 96.D0*Hr2(1,1)*z2*dx - 16.D0*Hr3(
     &    -1,-1,0) - 160.D0*Hr3(-1,-1,0)*x + 64.D0/3.D0*Hr3(-1,-1,0)*
     &    x**2 + 352.D0/3.D0*Hr3(-1,-1,0)*dx - 328.D0/3.D0*Hr3(-1,0,0)
     &     + 544.D0/3.D0*Hr3(-1,0,0)*x - 128.D0/3.D0*Hr3(-1,0,0)*x**2
     &     - 272.D0*Hr3(-1,0,0)*dx + 160.D0*Hr3(-1,0,1) + 80.D0*Hr3(-1,
     &    0,1)*x + 64.D0/3.D0*Hr3(-1,0,1)*x**2 + 352.D0/3.D0*Hr3(-1,0,1
     &    )*dx - 280.D0/3.D0*Hr3(0,-1,0) + 928.D0/3.D0*Hr3(0,-1,0)*x - 
     &    256.D0/3.D0*Hr3(0,-1,0)*x**2 - 752.D0/3.D0*Hr3(0,-1,0)*dx - 
     &    848.D0/3.D0*Hr3(0,0,0) - 3344.D0/3.D0*Hr3(0,0,0)*x - 224.D0/3.
     &    D0*Hr3(0,0,0)*x**2 + 1856.D0/3.D0*Hr3(0,0,0)*dx - 1292.D0/3.D0
     &    *Hr3(0,0,1) - 284.D0/3.D0*Hr3(0,0,1)*x - 128.D0/3.D0*Hr3(0,0,
     &    1)*x**2 - 32.D0/3.D0*Hr3(0,0,1)*dx - 576.D0*Hr3(0,1,0) - 424.D
     &    0*Hr3(0,1,0)*x - 512.D0/3.D0*Hr3(0,1,0)*x**2 + 1504.D0/3.D0*
     &    Hr3(0,1,0)*dx )
      PTgq2 = PTgq2 + cf*ca**2 * (  - 296.D0/3.D0*Hr3(0,1,1) - 260.D0/3.
     &    D0*Hr3(0,1,1)*x - 128.D0/3.D0*Hr3(0,1,1)*x**2 + 448.D0/3.D0*
     &    Hr3(0,1,1)*dx + 1300.D0/3.D0*Hr3(1,0,0) - 1004.D0/3.D0*Hr3(1,
     &    0,0)*x - 1024.D0/3.D0*Hr3(1,0,0)*dx + 192.D0*Hr3(1,0,1) - 72.D
     &    0*Hr3(1,0,1)*x + 32.D0*Hr3(1,0,1)*x**2 - 248.D0*Hr3(1,0,1)*dx
     &     + 368.D0/3.D0*Hr3(1,1,0) - 352.D0/3.D0*Hr3(1,1,0)*x + 32.D0/
     &    3.D0*Hr3(1,1,0)*x**2 - 424.D0/3.D0*Hr3(1,1,0)*dx + 104.D0/3.D0
     &    *Hr3(1,1,1) - 76.D0/3.D0*Hr3(1,1,1)*x + 32.D0/3.D0*Hr3(1,1,1)
     &    *x**2 - 160.D0/3.D0*Hr3(1,1,1)*dx - 192.D0*Hr4(-1,-1,-1,0) - 
     &    96.D0*Hr4(-1,-1,-1,0)*x - 192.D0*Hr4(-1,-1,-1,0)*dx + 192.D0*
     &    Hr4(-1,-1,0,0) + 96.D0*Hr4(-1,-1,0,0)*x + 192.D0*Hr4(-1,-1,0,
     &    0)*dx - 64.D0*Hr4(-1,-1,0,1) - 32.D0*Hr4(-1,-1,0,1)*x - 64.D0
     &    *Hr4(-1,-1,0,1)*dx + 352.D0*Hr4(-1,0,-1,0) + 176.D0*Hr4(-1,0,
     &    -1,0)*x + 352.D0*Hr4(-1,0,-1,0)*dx - 336.D0*Hr4(-1,0,0,0) - 
     &    168.D0*Hr4(-1,0,0,0)*x - 336.D0*Hr4(-1,0,0,0)*dx + 64.D0*Hr4(
     &    -1,0,0,1) )
      PTgq2 = PTgq2 + cf*ca**2 * ( 32.D0*Hr4(-1,0,0,1)*x + 64.D0*Hr4(-1
     &    ,0,0,1)*dx - 192.D0*Hr4(-1,0,1,0) - 96.D0*Hr4(-1,0,1,0)*x - 
     &    192.D0*Hr4(-1,0,1,0)*dx - 64.D0*Hr4(-1,0,1,1) - 32.D0*Hr4(-1,
     &    0,1,1)*x - 64.D0*Hr4(-1,0,1,1)*dx - 160.D0*Hr4(0,-1,-1,0) + 
     &    48.D0*Hr4(0,-1,-1,0)*x + 240.D0*Hr4(0,-1,0,0) - 136.D0*Hr4(0,
     &    -1,0,0)*x - 96.D0*Hr4(0,-1,0,0)*dx + 64.D0*Hr4(0,-1,0,1) + 32.
     &    D0*Hr4(0,-1,0,1)*x + 64.D0*Hr4(0,-1,0,1)*dx + 112.D0*Hr4(0,0,
     &    -1,0) - 248.D0*Hr4(0,0,-1,0)*x - 256.D0*Hr4(0,0,-1,0)*dx - 
     &    128.D0*Hr4(0,0,0,0) + 632.D0*Hr4(0,0,0,0)*x + 512.D0*Hr4(0,0,
     &    0,0)*dx - 152.D0*Hr4(0,0,0,1) - 92.D0*Hr4(0,0,0,1)*x - 128.D0
     &    *Hr4(0,0,0,1)*dx + 384.D0*Hr4(0,0,1,0) + 480.D0*Hr4(0,0,1,0)*
     &    x + 576.D0*Hr4(0,0,1,0)*dx + 128.D0*Hr4(0,0,1,1) + 160.D0*
     &    Hr4(0,0,1,1)*x + 192.D0*Hr4(0,0,1,1)*dx + 232.D0*Hr4(0,1,0,0)
     &     - 76.D0*Hr4(0,1,0,0)*x - 160.D0*Hr4(0,1,0,0)*dx + 64.D0*Hr4(
     &    0,1,0,1) - 176.D0*Hr4(0,1,0,1)*x - 256.D0*Hr4(0,1,0,1)*dx + 
     &    128.D0*Hr4(0,1,1,0) )
      PTgq2 = PTgq2 + cf*ca**2 * (  - 112.D0*Hr4(0,1,1,0)*x - 192.D0*
     &    Hr4(0,1,1,0)*dx + 64.D0*Hr4(0,1,1,1) - 80.D0*Hr4(0,1,1,1)*x
     &     - 128.D0*Hr4(0,1,1,1)*dx + 96.D0*Hr4(1,0,-1,0) - 48.D0*Hr4(1
     &    ,0,-1,0)*x - 96.D0*Hr4(1,0,-1,0)*dx - 560.D0*Hr4(1,0,0,0) + 
     &    280.D0*Hr4(1,0,0,0)*x + 560.D0*Hr4(1,0,0,0)*dx - 160.D0*Hr4(1
     &    ,0,0,1) + 80.D0*Hr4(1,0,0,1)*x + 160.D0*Hr4(1,0,0,1)*dx - 640.
     &    D0*Hr4(1,0,1,0) + 320.D0*Hr4(1,0,1,0)*x + 640.D0*Hr4(1,0,1,0)
     &    *dx - 64.D0*Hr4(1,0,1,1) + 32.D0*Hr4(1,0,1,1)*x + 64.D0*Hr4(1
     &    ,0,1,1)*dx - 160.D0*Hr4(1,1,0,0) + 80.D0*Hr4(1,1,0,0)*x + 160.
     &    D0*Hr4(1,1,0,0)*dx - 128.D0*Hr4(1,1,0,1) + 64.D0*Hr4(1,1,0,1)
     &    *x + 128.D0*Hr4(1,1,0,1)*dx - 192.D0*Hr4(1,1,1,0) + 96.D0*
     &    Hr4(1,1,1,0)*x + 192.D0*Hr4(1,1,1,0)*dx - 64.D0*Hr4(1,1,1,1)
     &     + 32.D0*Hr4(1,1,1,1)*x + 64.D0*Hr4(1,1,1,1)*dx )
      PTgq2 = PTgq2 + cf**2*ca * ( 1735.D0/6.D0 - 4811.D0/12.D0*x - 140.
     &    D0/3.D0*x**2 + 200.D0*dx + 552.D0*z3 + 424.D0*z3*x + 128.D0/3.
     &    D0*z3*x**2 + 264.D0*z3*dx + 3452.D0/3.D0*z2 + 436.D0/3.D0*z2*
     &    x + 1760.D0/9.D0*z2*x**2 + 528.D0*z2*dx - 608.D0/5.D0*z2**2
     &     - 2144.D0/5.D0*z2**2*x - 3664.D0/5.D0*z2**2*dx - 176.D0*Hr1(
     &    -1)*z3 - 88.D0*Hr1(-1)*z3*x - 176.D0*Hr1(-1)*z3*dx - 64.D0*
     &    Hr1(-1)*z2 - 32.D0*Hr1(-1)*z2*x - 24.D0*Hr1(-1)*z2*dx + 40301.
     &    D0/27.D0*Hr1(0) + 18911.D0/27.D0*Hr1(0)*x + 8932.D0/27.D0*
     &    Hr1(0)*x**2 + 244.D0*Hr1(0)*dx - 272.D0*Hr1(0)*z3 - 376.D0*
     &    Hr1(0)*z3*x - 832.D0*Hr1(0)*z3*dx + 72.D0*Hr1(0)*z2 + 128.D0*
     &    Hr1(0)*z2*x - 32.D0/3.D0*Hr1(0)*z2*x**2 + 192.D0*Hr1(0)*z2*dx
     &     + 26008.D0/27.D0*Hr1(1) + 3988.D0/27.D0*Hr1(1)*x + 5908.D0/
     &    27.D0*Hr1(1)*x**2 - 39092.D0/27.D0*Hr1(1)*dx + 816.D0*Hr1(1)*
     &    z3 - 408.D0*Hr1(1)*z3*x - 816.D0*Hr1(1)*z3*dx + 672.D0*Hr1(1)
     &    *z2 - 180.D0*Hr1(1)*z2*x + 224.D0/3.D0*Hr1(1)*z2*x**2 - 1760.D
     &    0/3.D0*Hr1(1)*z2*dx )
      PTgq2 = PTgq2 + cf**2*ca * ( 288.D0*Hr2(-1,-1)*z2 + 144.D0*Hr2(-1
     &    ,-1)*z2*x + 288.D0*Hr2(-1,-1)*z2*dx - 88.D0*Hr2(-1,0) - 268.D0
     &    *Hr2(-1,0)*x + 192.D0*Hr2(-1,0)*dx - 192.D0*Hr2(-1,0)*z2 - 96.
     &    D0*Hr2(-1,0)*z2*x - 192.D0*Hr2(-1,0)*z2*dx - 128.D0*Hr2(0,-1)
     &    *z2 - 16.D0*Hr2(0,-1)*z2*x + 32.D0*Hr2(0,-1)*z2*dx - 2338.D0/
     &    3.D0*Hr2(0,0) - 441.D0*Hr2(0,0)*x - 1408.D0/9.D0*Hr2(0,0)*
     &    x**2 + 80.D0*Hr2(0,0)*dx - 16.D0*Hr2(0,0)*z2 - 40.D0*Hr2(0,0)
     &    *z2*x - 128.D0*Hr2(0,0)*z2*dx - 8980.D0/9.D0*Hr2(0,1) - 4576.D
     &    0/9.D0*Hr2(0,1)*x - 352.D0/3.D0*Hr2(0,1)*x**2 - 5800.D0/9.D0*
     &    Hr2(0,1)*dx + 272.D0*Hr2(0,1)*z2 - 456.D0*Hr2(0,1)*z2*x - 608.
     &    D0*Hr2(0,1)*z2*dx + 488.D0/3.D0*Hr2(1,0) - 200.D0*Hr2(1,0)*x
     &     - 704.D0/9.D0*Hr2(1,0)*x**2 + 620.D0/9.D0*Hr2(1,0)*dx - 32.D0
     &    *Hr2(1,0)*z2 + 16.D0*Hr2(1,0)*z2*x + 32.D0*Hr2(1,0)*z2*dx - 
     &    160.D0/9.D0*Hr2(1,1) - 280.D0/9.D0*Hr2(1,1)*x - 352.D0/9.D0*
     &    Hr2(1,1)*x**2 + 848.D0/9.D0*Hr2(1,1)*dx + 192.D0*Hr2(1,1)*z2
     &     - 96.D0*Hr2(1,1)*z2*x )
      PTgq2 = PTgq2 + cf**2*ca * (  - 192.D0*Hr2(1,1)*z2*dx + 192.D0*
     &    Hr3(-1,-1,0) + 128.D0*Hr3(-1,-1,0)*x + 48.D0*Hr3(-1,-1,0)*dx
     &     - 192.D0*Hr3(-1,0,0) - 72.D0*Hr3(-1,0,0)*x - 48.D0*Hr3(-1,0,
     &    0)*dx + 160.D0*Hr3(-1,0,1) + 96.D0*Hr3(-1,0,1)*x + 48.D0*Hr3(
     &    -1,0,1)*dx - 144.D0*Hr3(0,-1,0) + 128.D0*Hr3(0,-1,0)*x + 96.D0
     &    *Hr3(0,-1,0)*dx + 368.D0/3.D0*Hr3(0,0,0) - 268.D0/3.D0*Hr3(0,
     &    0,0)*x + 64.D0/3.D0*Hr3(0,0,0)*x**2 + 136.D0/3.D0*Hr3(0,0,1)
     &     - 680.D0/3.D0*Hr3(0,0,1)*x - 32.D0*Hr3(0,0,1)*x**2 + 640.D0/
     &    3.D0*Hr3(0,0,1)*dx + 600.D0*Hr3(0,1,0) + 96.D0*Hr3(0,1,0)*x
     &     + 224.D0/3.D0*Hr3(0,1,0)*x**2 - 1280.D0/3.D0*Hr3(0,1,0)*dx
     &     + 424.D0/3.D0*Hr3(0,1,1) + 28.D0/3.D0*Hr3(0,1,1)*x + 32.D0*
     &    Hr3(0,1,1)*x**2 - 464.D0/3.D0*Hr3(0,1,1)*dx - 1760.D0/3.D0*
     &    Hr3(1,0,0) + 832.D0/3.D0*Hr3(1,0,0)*x - 64.D0/3.D0*Hr3(1,0,0)
     &    *x**2 + 336.D0*Hr3(1,0,0)*dx - 976.D0/3.D0*Hr3(1,0,1) + 404.D0
     &    /3.D0*Hr3(1,0,1)*x - 32.D0*Hr3(1,0,1)*x**2 + 904.D0/3.D0*Hr3(
     &    1,0,1)*dx )
      PTgq2 = PTgq2 + cf**2*ca * (  - 256.D0/3.D0*Hr3(1,1,0) + 332.D0/3.
     &    D0*Hr3(1,1,0)*x - 32.D0/3.D0*Hr3(1,1,0)*x**2 + 88.D0*Hr3(1,1,
     &    0)*dx - 152.D0/3.D0*Hr3(1,1,1) + 232.D0/3.D0*Hr3(1,1,1)*x - 
     &    32.D0/3.D0*Hr3(1,1,1)*x**2 + 160.D0/3.D0*Hr3(1,1,1)*dx + 192.D
     &    0*Hr4(-1,-1,-1,0) + 96.D0*Hr4(-1,-1,-1,0)*x + 192.D0*Hr4(-1,
     &    -1,-1,0)*dx + 64.D0*Hr4(-1,-1,0,0) + 32.D0*Hr4(-1,-1,0,0)*x
     &     + 64.D0*Hr4(-1,-1,0,0)*dx - 192.D0*Hr4(-1,-1,0,1) - 96.D0*
     &    Hr4(-1,-1,0,1)*x - 192.D0*Hr4(-1,-1,0,1)*dx - 160.D0*Hr4(-1,0
     &    ,-1,0) - 80.D0*Hr4(-1,0,-1,0)*x - 160.D0*Hr4(-1,0,-1,0)*dx - 
     &    16.D0*Hr4(-1,0,0,0) - 8.D0*Hr4(-1,0,0,0)*x - 16.D0*Hr4(-1,0,0
     &    ,0)*dx - 64.D0*Hr4(-1,0,0,1) - 32.D0*Hr4(-1,0,0,1)*x - 64.D0*
     &    Hr4(-1,0,0,1)*dx + 128.D0*Hr4(-1,0,1,0) + 64.D0*Hr4(-1,0,1,0)
     &    *x + 128.D0*Hr4(-1,0,1,0)*dx + 64.D0*Hr4(-1,0,1,1) + 32.D0*
     &    Hr4(-1,0,1,1)*x + 64.D0*Hr4(-1,0,1,1)*dx - 128.D0*Hr4(0,-1,-1
     &    ,0) - 96.D0*Hr4(0,-1,-1,0)*x - 64.D0*Hr4(0,-1,-1,0)*dx - 32.D0
     &    *Hr4(0,-1,0,0) )
      PTgq2 = PTgq2 + cf**2*ca * (  - 32.D0*Hr4(0,-1,0,0)*x + 64.D0*
     &    Hr4(0,-1,0,1) - 32.D0*Hr4(0,-1,0,1)*x - 64.D0*Hr4(0,-1,0,1)*
     &    dx - 32.D0*Hr4(0,0,-1,0) + 48.D0*Hr4(0,0,-1,0)*x + 64.D0*Hr4(
     &    0,0,0,0) + 64.D0*Hr4(0,0,0,0)*x + 336.D0*Hr4(0,0,0,1) + 376.D0
     &    *Hr4(0,0,0,1)*x + 512.D0*Hr4(0,0,0,1)*dx - 208.D0*Hr4(0,0,1,0
     &    ) - 344.D0*Hr4(0,0,1,0)*x - 448.D0*Hr4(0,0,1,0)*dx - 16.D0*
     &    Hr4(0,0,1,1) - 216.D0*Hr4(0,0,1,1)*x - 320.D0*Hr4(0,0,1,1)*dx
     &     - 320.D0*Hr4(0,1,0,0) + 208.D0*Hr4(0,1,0,0)*x + 256.D0*Hr4(0
     &    ,1,0,0)*dx - 336.D0*Hr4(0,1,0,1) + 312.D0*Hr4(0,1,0,1)*x + 
     &    448.D0*Hr4(0,1,0,1)*dx - 336.D0*Hr4(0,1,1,0) + 216.D0*Hr4(0,1
     &    ,1,0)*x + 384.D0*Hr4(0,1,1,0)*dx - 272.D0*Hr4(0,1,1,1) + 184.D
     &    0*Hr4(0,1,1,1)*x + 320.D0*Hr4(0,1,1,1)*dx - 96.D0*Hr4(1,0,-1,
     &    0) + 48.D0*Hr4(1,0,-1,0)*x + 96.D0*Hr4(1,0,-1,0)*dx + 48.D0*
     &    Hr4(1,0,0,0) - 24.D0*Hr4(1,0,0,0)*x - 48.D0*Hr4(1,0,0,0)*dx
     &     - 224.D0*Hr4(1,0,0,1) + 112.D0*Hr4(1,0,0,1)*x + 224.D0*Hr4(1
     &    ,0,0,1)*dx )
      PTgq2 = PTgq2 + cf**2*ca * ( 544.D0*Hr4(1,0,1,0) - 272.D0*Hr4(1,0
     &    ,1,0)*x - 544.D0*Hr4(1,0,1,0)*dx - 64.D0*Hr4(1,0,1,1) + 32.D0
     &    *Hr4(1,0,1,1)*x + 64.D0*Hr4(1,0,1,1)*dx - 160.D0*Hr4(1,1,0,0)
     &     + 80.D0*Hr4(1,1,0,0)*x + 160.D0*Hr4(1,1,0,0)*dx + 32.D0*Hr4(
     &    1,1,0,1) - 16.D0*Hr4(1,1,0,1)*x - 32.D0*Hr4(1,1,0,1)*dx + 224.
     &    D0*Hr4(1,1,1,0) - 112.D0*Hr4(1,1,1,0)*x - 224.D0*Hr4(1,1,1,0)
     &    *dx + 128.D0*Hr4(1,1,1,1) - 64.D0*Hr4(1,1,1,1)*x - 128.D0*
     &    Hr4(1,1,1,1)*dx )
      PTgq2 = PTgq2 + cf**3 * (  - 1915.D0/2.D0 + 731.D0/4.D0*x + 794.D0
     &    *dx - 96.D0*z3 - 196.D0*z3*x - 720.D0*z3*dx - 224.D0*z2 + 148.
     &    D0*z2*x - 256.D0*z2*dx - 776.D0/5.D0*z2**2 + 84.D0/5.D0*z2**2
     &    *x + 1216.D0/5.D0*z2**2*dx - 64.D0*Hr1(-1)*z2 - 32.D0*Hr1(-1)
     &    *z2*x + 197.D0*Hr1(0) + 29.D0*Hr1(0)*x + 208.D0*Hr1(0)*dx - 
     &    80.D0*Hr1(0)*z3 - 200.D0*Hr1(0)*z3*x + 128.D0*Hr1(0)*z3*dx - 
     &    24.D0*Hr1(0)*z2 - 20.D0*Hr1(0)*z2*x - 192.D0*Hr1(0)*z2*dx - 
     &    230.D0*Hr1(1) + 16.D0*Hr1(1)*x + 228.D0*Hr1(1)*dx - 64.D0*
     &    Hr1(1)*z3 + 32.D0*Hr1(1)*z3*x + 64.D0*Hr1(1)*z3*dx - 352.D0*
     &    Hr1(1)*z2 + 252.D0*Hr1(1)*z2*x + 144.D0*Hr1(1)*z2*dx - 80.D0*
     &    Hr2(-1,0) + 208.D0*Hr2(-1,0)*x - 256.D0*Hr2(-1,0)*dx + 128.D0
     &    *Hr2(0,-1)*z2 + 96.D0*Hr2(0,-1)*z2*x + 128.D0*Hr2(0,-1)*z2*dx
     &     - 290.D0*Hr2(0,0) - 251.D0*Hr2(0,0)*x - 64.D0*Hr2(0,0)*z2 - 
     &    64.D0*Hr2(0,0)*z2*x + 208.D0*Hr2(0,1) - 140.D0*Hr2(0,1)*x - 
     &    464.D0*Hr2(0,1)*z2 + 200.D0*Hr2(0,1)*z2*x + 320.D0*Hr2(0,1)*
     &    z2*dx )
      PTgq2 = PTgq2 + cf**3 * (  - 264.D0*Hr2(1,0) + 140.D0*Hr2(1,0)*x
     &     + 180.D0*Hr2(1,0)*dx - 128.D0*Hr2(1,0)*z2 + 64.D0*Hr2(1,0)*
     &    z2*x + 128.D0*Hr2(1,0)*z2*dx - 72.D0*Hr2(1,1) + 80.D0*Hr2(1,1
     &    )*x + 60.D0*Hr2(1,1)*dx - 96.D0*Hr2(1,1)*z2 + 48.D0*Hr2(1,1)*
     &    z2*x + 96.D0*Hr2(1,1)*z2*dx - 128.D0*Hr3(-1,-1,0) - 64.D0*
     &    Hr3(-1,-1,0)*x + 224.D0*Hr3(-1,0,0) + 96.D0*Hr3(-1,0,0)*x + 
     &    96.D0*Hr3(-1,0,0)*dx + 128.D0*Hr3(0,-1,0) - 112.D0*Hr3(0,-1,0
     &    )*x - 192.D0*Hr3(0,-1,0)*dx + 160.D0*Hr3(0,0,0) + 40.D0*Hr3(0
     &    ,0,0)*x + 248.D0*Hr3(0,0,1) - 28.D0*Hr3(0,0,1)*x - 104.D0*
     &    Hr3(0,1,0) + 40.D0*Hr3(0,1,0)*x - 72.D0*Hr3(0,1,1) + 56.D0*
     &    Hr3(0,1,1)*x + 128.D0*Hr3(1,0,0) - 36.D0*Hr3(1,0,0)*x + 48.D0
     &    *Hr3(1,0,0)*dx + 128.D0*Hr3(1,0,1) - 60.D0*Hr3(1,0,1)*x - 48.D
     &    0*Hr3(1,0,1)*dx - 32.D0*Hr3(1,1,0) + 4.D0*Hr3(1,1,0)*x + 48.D0
     &    *Hr3(1,1,0)*dx + 16.D0*Hr3(1,1,1) - 52.D0*Hr3(1,1,1)*x - 128.D
     &    0*Hr4(-1,-1,0,0) - 64.D0*Hr4(-1,-1,0,0)*x - 128.D0*Hr4(-1,-1,
     &    0,0)*dx )
      PTgq2 = PTgq2 + cf**3 * (  - 64.D0*Hr4(-1,0,-1,0) - 32.D0*Hr4(-1,
     &    0,-1,0)*x - 64.D0*Hr4(-1,0,-1,0)*dx + 32.D0*Hr4(-1,0,0,0) + 
     &    16.D0*Hr4(-1,0,0,0)*x + 32.D0*Hr4(-1,0,0,0)*dx + 256.D0*Hr4(0
     &    ,-1,-1,0) + 192.D0*Hr4(0,-1,-1,0)*x + 256.D0*Hr4(0,-1,-1,0)*
     &    dx - 64.D0*Hr4(0,-1,0,0) - 64.D0*Hr4(0,-1,0,0)*x - 128.D0*
     &    Hr4(0,-1,0,0)*dx - 64.D0*Hr4(0,0,-1,0) - 96.D0*Hr4(0,0,-1,0)*
     &    x - 32.D0*Hr4(0,0,0,0) + 96.D0*Hr4(0,0,0,0)*x - 128.D0*Hr4(0,
     &    0,0,1) + 64.D0*Hr4(0,0,0,1)*x - 176.D0*Hr4(0,0,1,0) + 88.D0*
     &    Hr4(0,0,1,0)*x + 128.D0*Hr4(0,0,1,0)*dx - 240.D0*Hr4(0,0,1,1)
     &     + 120.D0*Hr4(0,0,1,1)*x + 256.D0*Hr4(0,0,1,1)*dx - 16.D0*
     &    Hr4(0,1,0,0) + 40.D0*Hr4(0,1,0,0)*x + 128.D0*Hr4(0,1,0,0)*dx
     &     + 144.D0*Hr4(0,1,0,1) - 72.D0*Hr4(0,1,0,1)*x - 64.D0*Hr4(0,1
     &    ,0,1)*dx + 80.D0*Hr4(0,1,1,0) - 40.D0*Hr4(0,1,1,0)*x - 64.D0*
     &    Hr4(0,1,1,0)*dx + 208.D0*Hr4(0,1,1,1) - 104.D0*Hr4(0,1,1,1)*x
     &     - 192.D0*Hr4(0,1,1,1)*dx + 64.D0*Hr4(1,0,-1,0) - 32.D0*Hr4(1
     &    ,0,-1,0)*x )
      PTgq2 = PTgq2 + cf**3 * (  - 64.D0*Hr4(1,0,-1,0)*dx + 192.D0*Hr4(
     &    1,0,0,0) - 96.D0*Hr4(1,0,0,0)*x - 192.D0*Hr4(1,0,0,0)*dx + 
     &    256.D0*Hr4(1,0,0,1) - 128.D0*Hr4(1,0,0,1)*x - 256.D0*Hr4(1,0,
     &    0,1)*dx - 32.D0*Hr4(1,0,1,0) + 16.D0*Hr4(1,0,1,0)*x + 32.D0*
     &    Hr4(1,0,1,0)*dx + 128.D0*Hr4(1,0,1,1) - 64.D0*Hr4(1,0,1,1)*x
     &     - 128.D0*Hr4(1,0,1,1)*dx + 192.D0*Hr4(1,1,0,0) - 96.D0*Hr4(1
     &    ,1,0,0)*x - 192.D0*Hr4(1,1,0,0)*dx + 96.D0*Hr4(1,1,0,1) - 48.D
     &    0*Hr4(1,1,0,1)*x - 96.D0*Hr4(1,1,0,1)*dx - 32.D0*Hr4(1,1,1,0)
     &     + 16.D0*Hr4(1,1,1,0)*x + 32.D0*Hr4(1,1,1,0)*dx - 64.D0*Hr4(1
     &    ,1,1,1) + 32.D0*Hr4(1,1,1,1)*x + 64.D0*Hr4(1,1,1,1)*dx )
      PTgq2 = PTgq2 + nf*cf*ca * ( 332.D0/27.D0 - 2624.D0/27.D0*x + 
     &    1892.D0/81.D0*x**2 + 3778.D0/81.D0*dx + 64.D0*z3 - 40.D0*z3*x
     &     - 80.D0*z3*dx + 8.D0/3.D0*z2 - 176.D0/9.D0*z2*x - 160.D0/9.D0
     &    *z2*dx + 1280.D0/27.D0*Hr1(0) + 308.D0/27.D0*Hr1(0)*x + 64.D0/
     &    9.D0*Hr1(0)*x**2 + 1156.D0/27.D0*Hr1(0)*dx - 16.D0/3.D0*Hr1(0
     &    )*z2*x - 32.D0/3.D0*Hr1(0)*z2*dx + 448.D0/27.D0*Hr1(1) - 392.D
     &    0/27.D0*Hr1(1)*x - 448.D0/27.D0*Hr1(1)*dx - 32.D0*Hr1(1)*z2
     &     + 16.D0*Hr1(1)*z2*x + 32.D0*Hr1(1)*z2*dx - 160.D0/9.D0*Hr2(
     &    -1,0) - 104.D0/9.D0*Hr2(-1,0)*x - 160.D0/9.D0*Hr2(-1,0)*dx - 
     &    256.D0/9.D0*Hr2(0,0) - 464.D0/9.D0*Hr2(0,0)*x - 128.D0/9.D0*
     &    Hr2(0,0)*x**2 + 656.D0/9.D0*Hr2(0,0)*dx + 136.D0/9.D0*Hr2(0,1
     &    ) - 104.D0/9.D0*Hr2(0,1)*x - 160.D0/9.D0*Hr2(0,1)*dx - 320.D0/
     &    9.D0*Hr2(1,0) + 136.D0/9.D0*Hr2(1,0)*x + 320.D0/9.D0*Hr2(1,0)
     &    *dx - 160.D0/9.D0*Hr2(1,1) + 104.D0/9.D0*Hr2(1,1)*x + 160.D0/
     &    9.D0*Hr2(1,1)*dx - 32.D0/3.D0*Hr3(-1,0,0) - 16.D0/3.D0*Hr3(-1
     &    ,0,0)*x )
      PTgq2 = PTgq2 + nf*cf*ca * (  - 32.D0/3.D0*Hr3(-1,0,0)*dx - 32.D0/
     &    3.D0*Hr3(0,-1,0) - 16.D0/3.D0*Hr3(0,-1,0)*x - 32.D0/3.D0*Hr3(
     &    0,-1,0)*dx + 32.D0/3.D0*Hr3(0,0,0) + 128.D0/3.D0*Hr3(0,0,0)*x
     &     + 128.D0/3.D0*Hr3(0,0,0)*dx + 32.D0/3.D0*Hr3(0,0,1) - 16.D0/
     &    3.D0*Hr3(0,0,1)*x - 32.D0/3.D0*Hr3(0,0,1)*dx - 16.D0/3.D0*
     &    Hr3(0,1,1) + 8.D0/3.D0*Hr3(0,1,1)*x + 16.D0/3.D0*Hr3(0,1,1)*
     &    dx - 160.D0/3.D0*Hr3(1,0,0) + 80.D0/3.D0*Hr3(1,0,0)*x + 160.D0
     &    /3.D0*Hr3(1,0,0)*dx - 32.D0/3.D0*Hr3(1,1,0) + 16.D0/3.D0*Hr3(
     &    1,1,0)*x + 32.D0/3.D0*Hr3(1,1,0)*dx + 16.D0/3.D0*Hr3(1,1,1)
     &     - 8.D0/3.D0*Hr3(1,1,1)*x - 16.D0/3.D0*Hr3(1,1,1)*dx )
      PTgq2 = PTgq2 + nf*cf**2 * (  - 12803.D0/27.D0 + 30475.D0/54.D0*x
     &     - 21784.D0/81.D0*x**2 + 14008.D0/81.D0*dx - 2330.D0/27.D0*
     &    Hr1(0) + 2494.D0/27.D0*Hr1(0)*x + 3040.D0/27.D0*Hr1(0)*x**2
     &     + 200.D0/9.D0*Hr1(0)*dx - 484.D0/27.D0*Hr1(1) + 464.D0/27.D0
     &    *Hr1(1)*x + 448.D0/27.D0*Hr1(1)*dx + 16.D0*Hr1(1)*z2 - 8.D0*
     &    Hr1(1)*z2*x - 16.D0*Hr1(1)*z2*dx - 204.D0*Hr2(0,0) + 142.D0/3.
     &    D0*Hr2(0,0)*x - 256.D0/9.D0*Hr2(0,0)*x**2 - 128.D0/3.D0*Hr2(0
     &    ,0)*dx - 320.D0/9.D0*Hr2(0,1) + 208.D0/9.D0*Hr2(0,1)*x + 320.D
     &    0/9.D0*Hr2(0,1)*dx + 80.D0/3.D0*Hr2(1,0) - 40.D0/3.D0*Hr2(1,0
     &    )*x - 80.D0/3.D0*Hr2(1,0)*dx + 160.D0/9.D0*Hr2(1,1) - 104.D0/
     &    9.D0*Hr2(1,1)*x - 160.D0/9.D0*Hr2(1,1)*dx - 176.D0/3.D0*Hr3(0
     &    ,0,0) + 160.D0/3.D0*Hr3(0,0,0)*x - 128.D0/3.D0*Hr3(0,0,0)*dx
     &     - 64.D0/3.D0*Hr3(0,0,1) + 32.D0/3.D0*Hr3(0,0,1)*x + 64.D0/3.D
     &    0*Hr3(0,0,1)*dx + 32.D0/3.D0*Hr3(0,1,1) - 16.D0/3.D0*Hr3(0,1,
     &    1)*x - 32.D0/3.D0*Hr3(0,1,1)*dx + 80.D0/3.D0*Hr3(1,0,0) - 40.D
     &    0/3.D0*Hr3(1,0,0)*x )
      PTgq2 = PTgq2 + nf*cf**2 * (  - 80.D0/3.D0*Hr3(1,0,0)*dx + 16.D0/
     &    3.D0*Hr3(1,0,1) - 8.D0/3.D0*Hr3(1,0,1)*x - 16.D0/3.D0*Hr3(1,0
     &    ,1)*dx + 16.D0/3.D0*Hr3(1,1,0) - 8.D0/3.D0*Hr3(1,1,0)*x - 16.D
     &    0/3.D0*Hr3(1,1,0)*dx - 16.D0/3.D0*Hr3(1,1,1) + 8.D0/3.D0*Hr3(
     &    1,1,1)*x + 16.D0/3.D0*Hr3(1,1,1)*dx + 64.D0*Hr4(0,0,0,0) - 32.
     &    D0*Hr4(0,0,0,0)*x )
*
       X2GQTA = PTgq2
*
       RETURN
       END
* ---------------------------------------------------------------------
*
*
* ..The regular piece of the gluon-gluon splitting functions P_gg^(2)T
*
       FUNCTION X2GGTA (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,
     ,             Z5 = 1.0369 27755 14336 99263 D0 )
*
* ..The soft coefficient for use in X2GGTB and X2GGTC
*
       COMMON / P2GSOFT / A3G
*
* ...Colour factors
*
       CF  = 4./3.D0
       CA  = 3.D0
       NF2 = NF*NF
*
* ...Some abbreviations
*
       DX = 1.D0/X
       DM = 1.D0/(1.D0-X)
       DP = 1.D0/(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
*    (without the delta(1-x) part, but with the soft contribution)
*
      PTgg2 =
     &  + ca**3 * (  - 54088.D0/27.D0 + 49678.D0/27.D0*x - 146182.D0/81.
     &    D0*x**2 + 146182.D0/81.D0*dx + 490.D0/3.D0*dm + 1064.D0/3.D0*
     &    z3 + 704.D0/3.D0*z3*x + 1408.D0/3.D0*z3*x**2 - 1144.D0/3.D0*
     &    z3*dx + 176.D0/3.D0*z3*dm + 976.D0/3.D0*z2 - 2632.D0/9.D0*z2*
     &    x + 1072.D0/9.D0*z2*x**2 - 3112.D0/9.D0*z2*dx - 1072.D0/9.D0*
     &    z2*dp - 1072.D0/9.D0*z2*dm - 272.D0*z2**2 - 1152.D0/5.D0*
     &    z2**2*x - 416.D0/5.D0*z2**2*x**2 - 464.D0/5.D0*z2**2*dx + 88.D
     &    0*z2**2*dp - 24.D0/5.D0*z2**2*dm - 384.D0*Hr1(-1)*z3 - 192.D0
     &    *Hr1(-1)*z3*x - 192.D0*Hr1(-1)*z3*x**2 - 192.D0*Hr1(-1)*z3*dx
     &     + 192.D0*Hr1(-1)*z3*dp - 424.D0*Hr1(-1)*z2 - 424.D0*Hr1(-1)*
     &    z2*x - 88.D0*Hr1(-1)*z2*x**2 - 88.D0*Hr1(-1)*z2*dx + 1274.D0/
     &    9.D0*Hr1(0) + 9382.D0/9.D0*Hr1(0)*x + 5524.D0/27.D0*Hr1(0)*
     &    x**2 + 32608.D0/27.D0*Hr1(0)*dx + 8272.D0/27.D0*Hr1(0)*dm - 
     &    800.D0*Hr1(0)*z3*x - 288.D0*Hr1(0)*z3*dx + 144.D0*Hr1(0)*z3*
     &    dp - 144.D0*Hr1(0)*z3*dm + 808.D0/3.D0*Hr1(0)*z2 + 512.D0/3.D0
     &    *Hr1(0)*z2*x )
      PTgg2 = PTgg2 + ca**3 * ( 176.D0/3.D0*Hr1(0)*z2*x**2 - 880.D0/3.D0
     &    *Hr1(0)*z2*dx - 88.D0*Hr1(0)*z2*dp - 176.D0/3.D0*Hr1(0)*z2*dm
     &     + 124.D0/9.D0*Hr1(1) - 124.D0/9.D0*Hr1(1)*x + 1652.D0/27.D0*
     &    Hr1(1)*x**2 - 1652.D0/27.D0*Hr1(1)*dx + 192.D0*Hr1(1)*z3 - 96.
     &    D0*Hr1(1)*z3*x + 96.D0*Hr1(1)*z3*x**2 - 96.D0*Hr1(1)*z3*dx - 
     &    96.D0*Hr1(1)*z3*dm - 24.D0*Hr1(1)*z2 + 24.D0*Hr1(1)*z2*x + 88.
     &    D0*Hr1(1)*z2*x**2 - 88.D0*Hr1(1)*z2*dx + 512.D0*Hr2(-1,-1)*z2
     &     + 256.D0*Hr2(-1,-1)*z2*x + 256.D0*Hr2(-1,-1)*z2*x**2 + 256.D0
     &    *Hr2(-1,-1)*z2*dx - 256.D0*Hr2(-1,-1)*z2*dp - 872.D0/9.D0*
     &    Hr2(-1,0) - 3016.D0/9.D0*Hr2(-1,0)*x - 680.D0/3.D0*Hr2(-1,0)*
     &    x**2 - 680.D0/3.D0*Hr2(-1,0)*dx - 2144.D0/9.D0*Hr2(-1,0)*dp
     &     - 64.D0*Hr2(-1,0)*z2 - 32.D0*Hr2(-1,0)*z2*x - 32.D0*Hr2(-1,0
     &    )*z2*x**2 - 32.D0*Hr2(-1,0)*z2*dx + 32.D0*Hr2(-1,0)*z2*dp - 
     &    256.D0*Hr2(0,-1)*z2 + 160.D0*Hr2(0,-1)*z2*x - 96.D0*Hr2(0,-1)
     &    *z2*x**2 + 96.D0*Hr2(0,-1)*z2*dx + 96.D0*Hr2(0,-1)*z2*dm - 
     &    116.D0/9.D0*Hr2(0,0) )
      PTgg2 = PTgg2 + ca**3 * ( 9044.D0/9.D0*Hr2(0,0)*x + 2248.D0/9.D0*
     &    Hr2(0,0)*x**2 + 6224.D0/9.D0*Hr2(0,0)*dx + 1072.D0/9.D0*Hr2(0
     &    ,0)*dp - 1280.D0/9.D0*Hr2(0,0)*dm + 96.D0*Hr2(0,0)*z2 - 608.D0
     &    *Hr2(0,0)*z2*x - 256.D0*Hr2(0,0)*z2*dx + 128.D0*Hr2(0,0)*z2*
     &    dp - 128.D0*Hr2(0,0)*z2*dm + 1760.D0/9.D0*Hr2(0,1) - 784.D0/9.
     &    D0*Hr2(0,1)*x + 1072.D0/9.D0*Hr2(0,1)*x**2 - 1072.D0/9.D0*
     &    Hr2(0,1)*dx - 2144.D0/9.D0*Hr2(0,1)*dm - 256.D0*Hr2(0,1)*z2
     &     - 32.D0*Hr2(0,1)*z2*x - 96.D0*Hr2(0,1)*z2*x**2 + 32.D0*Hr2(0
     &    ,1)*z2*dx + 32.D0*Hr2(0,1)*z2*dp + 64.D0*Hr2(0,1)*z2*dm + 
     &    2344.D0/9.D0*Hr2(1,0) - 200.D0/9.D0*Hr2(1,0)*x + 1072.D0/9.D0
     &    *Hr2(1,0)*x**2 - 1072.D0/9.D0*Hr2(1,0)*dx - 2144.D0/9.D0*Hr2(
     &    1,0)*dm - 64.D0*Hr2(1,0)*z2 + 32.D0*Hr2(1,0)*z2*x - 32.D0*
     &    Hr2(1,0)*z2*x**2 + 32.D0*Hr2(1,0)*z2*dx + 32.D0*Hr2(1,0)*z2*
     &    dm - 48.D0*Hr3(-1,-1,0) - 48.D0*Hr3(-1,-1,0)*x + 176.D0*Hr3(
     &    -1,-1,0)*x**2 + 176.D0*Hr3(-1,-1,0)*dx + 928.D0/3.D0*Hr3(-1,0
     &    ,0) )
      PTgg2 = PTgg2 + ca**3 * ( 224.D0/3.D0*Hr3(-1,0,0)*x - 704.D0/3.D0
     &    *Hr3(-1,0,0)*x**2 - 704.D0/3.D0*Hr3(-1,0,0)*dx - 704.D0/3.D0*
     &    Hr3(-1,0,0)*dp + 400.D0*Hr3(-1,0,1) + 400.D0*Hr3(-1,0,1)*x + 
     &    176.D0*Hr3(-1,0,1)*x**2 + 176.D0*Hr3(-1,0,1)*dx + 64.D0*Hr3(0
     &    ,-1,0) - 320.D0/3.D0*Hr3(0,-1,0)*x - 1232.D0/3.D0*Hr3(0,-1,0)
     &    *dx - 352.D0/3.D0*Hr3(0,-1,0)*dp + 176.D0/3.D0*Hr3(0,-1,0)*dm
     &     + 176.D0*Hr3(0,0,0) - 1264.D0/3.D0*Hr3(0,0,0)*x + 352.D0*
     &    Hr3(0,0,0)*x**2 + 704.D0*Hr3(0,0,0)*dx + 176.D0*Hr3(0,0,0)*dp
     &     - 528.D0*Hr3(0,0,0)*dm + 376.D0/3.D0*Hr3(0,0,1) - 512.D0/3.D0
     &    *Hr3(0,0,1)*x + 176.D0/3.D0*Hr3(0,0,1)*x**2 + 704.D0/3.D0*
     &    Hr3(0,0,1)*dx - 616.D0/3.D0*Hr3(0,0,1)*dm + 944.D0/3.D0*Hr3(0
     &    ,1,0) - 16.D0/3.D0*Hr3(0,1,0)*x + 352.D0/3.D0*Hr3(0,1,0)*x**2
     &     + 352.D0/3.D0*Hr3(0,1,0)*dx - 704.D0/3.D0*Hr3(0,1,0)*dm + 
     &    632.D0/3.D0*Hr3(1,0,0) - 16.D0/3.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 - 616.D0/3.D0*
     &    Hr3(1,0,0)*dm )
      PTgg2 = PTgg2 + ca**3 * ( 256.D0*Hr4(-1,-1,0,0) + 128.D0*Hr4(-1,
     &    -1,0,0)*x + 128.D0*Hr4(-1,-1,0,0)*x**2 + 128.D0*Hr4(-1,-1,0,0
     &    )*dx - 128.D0*Hr4(-1,-1,0,0)*dp - 512.D0*Hr4(-1,-1,0,1) - 256.
     &    D0*Hr4(-1,-1,0,1)*x - 256.D0*Hr4(-1,-1,0,1)*x**2 - 256.D0*
     &    Hr4(-1,-1,0,1)*dx + 256.D0*Hr4(-1,-1,0,1)*dp + 256.D0*Hr4(-1,
     &    0,-1,0) + 128.D0*Hr4(-1,0,-1,0)*x + 128.D0*Hr4(-1,0,-1,0)*
     &    x**2 + 128.D0*Hr4(-1,0,-1,0)*dx - 128.D0*Hr4(-1,0,-1,0)*dp - 
     &    640.D0*Hr4(-1,0,0,0) - 320.D0*Hr4(-1,0,0,0)*x - 320.D0*Hr4(-1
     &    ,0,0,0)*x**2 - 320.D0*Hr4(-1,0,0,0)*dx + 320.D0*Hr4(-1,0,0,0)
     &    *dp - 128.D0*Hr4(-1,0,1,0) - 64.D0*Hr4(-1,0,1,0)*x - 64.D0*
     &    Hr4(-1,0,1,0)*x**2 - 64.D0*Hr4(-1,0,1,0)*dx + 64.D0*Hr4(-1,0,
     &    1,0)*dp + 320.D0*Hr4(0,-1,-1,0)*x + 64.D0*Hr4(0,-1,-1,0)*x**2
     &     + 192.D0*Hr4(0,-1,-1,0)*dx - 128.D0*Hr4(0,-1,-1,0)*dp + 64.D0
     &    *Hr4(0,-1,-1,0)*dm + 128.D0*Hr4(0,-1,0,0) - 544.D0*Hr4(0,-1,0
     &    ,0)*x - 96.D0*Hr4(0,-1,0,0)*x**2 - 224.D0*Hr4(0,-1,0,0)*dx + 
     &    160.D0*Hr4(0,-1,0,0)*dp )
      PTgg2 = PTgg2 + ca**3 * (  - 64.D0*Hr4(0,-1,0,0)*dm + 256.D0*Hr4(
     &    0,-1,0,1) + 128.D0*Hr4(0,-1,0,1)*x**2 - 64.D0*Hr4(0,-1,0,1)*
     &    dp - 64.D0*Hr4(0,-1,0,1)*dm - 64.D0*Hr4(0,0,-1,0) - 704.D0*
     &    Hr4(0,0,-1,0)*x - 128.D0*Hr4(0,0,-1,0)*x**2 - 256.D0*Hr4(0,0,
     &    -1,0)*dx + 192.D0*Hr4(0,0,-1,0)*dp - 64.D0*Hr4(0,0,-1,0)*dm
     &     - 256.D0*Hr4(0,0,0,0) + 1984.D0*Hr4(0,0,0,0)*x - 128.D0*Hr4(
     &    0,0,0,0)*x**2 + 512.D0*Hr4(0,0,0,0)*dx - 192.D0*Hr4(0,0,0,0)*
     &    dp + 320.D0*Hr4(0,0,0,0)*dm - 96.D0*Hr4(0,0,0,1) + 800.D0*
     &    Hr4(0,0,0,1)*x - 256.D0*Hr4(0,0,0,1)*x**2 + 384.D0*Hr4(0,0,0,
     &    1)*dx - 64.D0*Hr4(0,0,0,1)*dp + 320.D0*Hr4(0,0,0,1)*dm - 128.D
     &    0*Hr4(0,0,1,0) + 512.D0*Hr4(0,0,1,0)*x - 192.D0*Hr4(0,0,1,0)*
     &    x**2 + 256.D0*Hr4(0,0,1,0)*dx - 32.D0*Hr4(0,0,1,0)*dp + 224.D0
     &    *Hr4(0,0,1,0)*dm - 256.D0*Hr4(0,0,1,1) + 128.D0*Hr4(0,0,1,1)*
     &    x - 128.D0*Hr4(0,0,1,1)*x**2 + 128.D0*Hr4(0,0,1,1)*dx + 128.D0
     &    *Hr4(0,0,1,1)*dm - 288.D0*Hr4(0,1,0,0) + 384.D0*Hr4(0,1,0,0)*
     &    x )
      PTgg2 = PTgg2 + ca**3 * (  - 224.D0*Hr4(0,1,0,0)*x**2 + 224.D0*
     &    Hr4(0,1,0,0)*dx + 224.D0*Hr4(0,1,0,0)*dm - 256.D0*Hr4(0,1,0,1
     &    ) + 128.D0*Hr4(0,1,0,1)*x - 128.D0*Hr4(0,1,0,1)*x**2 + 128.D0
     &    *Hr4(0,1,0,1)*dx + 128.D0*Hr4(0,1,0,1)*dm - 256.D0*Hr4(0,1,1,
     &    0) + 128.D0*Hr4(0,1,1,0)*x - 128.D0*Hr4(0,1,1,0)*x**2 + 128.D0
     &    *Hr4(0,1,1,0)*dx + 128.D0*Hr4(0,1,1,0)*dm + 128.D0*Hr4(1,0,-1
     &    ,0) - 64.D0*Hr4(1,0,-1,0)*x + 64.D0*Hr4(1,0,-1,0)*x**2 - 64.D0
     &    *Hr4(1,0,-1,0)*dx - 64.D0*Hr4(1,0,-1,0)*dm - 640.D0*Hr4(1,0,0
     &    ,0) + 320.D0*Hr4(1,0,0,0)*x - 320.D0*Hr4(1,0,0,0)*x**2 + 320.D
     &    0*Hr4(1,0,0,0)*dx + 320.D0*Hr4(1,0,0,0)*dm - 256.D0*Hr4(1,0,0
     &    ,1) + 128.D0*Hr4(1,0,0,1)*x - 128.D0*Hr4(1,0,0,1)*x**2 + 128.D
     &    0*Hr4(1,0,0,1)*dx + 128.D0*Hr4(1,0,0,1)*dm - 256.D0*Hr4(1,0,1
     &    ,0) + 128.D0*Hr4(1,0,1,0)*x - 128.D0*Hr4(1,0,1,0)*x**2 + 128.D
     &    0*Hr4(1,0,1,0)*dx + 128.D0*Hr4(1,0,1,0)*dm - 256.D0*Hr4(1,1,0
     &    ,0) + 128.D0*Hr4(1,1,0,0)*x - 128.D0*Hr4(1,1,0,0)*x**2 + 128.D
     &    0*Hr4(1,1,0,0)*dx )
      PTgg2 = PTgg2 + ca**3 * ( 128.D0*Hr4(1,1,0,0)*dm )
      PTgg2 = PTgg2 + nf*ca**2 * (  - 2174.D0/9.D0 + 7358.D0/27.D0*x - 
     &    19264.D0/81.D0*x**2 + 19264.D0/81.D0*dx - 836.D0/27.D0*dm + 
     &    496.D0/3.D0*z3 + 112.D0/3.D0*z3*x + 128.D0/3.D0*z3*x**2 - 176.
     &    D0/3.D0*z3*dx - 128.D0/3.D0*z3*dm - 1700.D0/9.D0*z2 + 1076.D0/
     &    9.D0*z2*x - 176.D0/3.D0*z2*x**2 - 128.D0/3.D0*z2*dx + 160.D0/
     &    9.D0*z2*dp + 160.D0/9.D0*z2*dm + 276.D0/5.D0*z2**2 + 156.D0/5.
     &    D0*z2**2*x - 88.D0*Hr1(-1)*z2 - 88.D0*Hr1(-1)*z2*x + 16.D0/3.D
     &    0*Hr1(-1)*z2*x**2 + 16.D0/3.D0*Hr1(-1)*z2*dx + 1408.D0/9.D0*
     &    Hr1(0) + 202.D0/3.D0*Hr1(0)*x + 340.D0/3.D0*Hr1(0)*x**2 + 
     &    5912.D0/27.D0*Hr1(0)*dx - 2600.D0/27.D0*Hr1(0)*dm - 96.D0*
     &    Hr1(0)*z3*x - 196.D0/3.D0*Hr1(0)*z2 + 76.D0/3.D0*Hr1(0)*z2*x
     &     - 32.D0/3.D0*Hr1(0)*z2*x**2 - 32.D0/3.D0*Hr1(0)*z2*dx + 16.D0
     &    *Hr1(0)*z2*dp + 32.D0/3.D0*Hr1(0)*z2*dm + 1354.D0/9.D0*Hr1(1)
     &     - 1354.D0/9.D0*Hr1(1)*x - 820.D0/27.D0*Hr1(1)*x**2 + 820.D0/
     &    27.D0*Hr1(1)*dx - 56.D0*Hr1(1)*z2 + 56.D0*Hr1(1)*z2*x + 16.D0
     &    *Hr1(1)*z2*x**2 )
      PTgg2 = PTgg2 + nf*ca**2 * (  - 16.D0*Hr1(1)*z2*dx + 1208.D0/9.D0
     &    *Hr2(-1,0) + 1528.D0/9.D0*Hr2(-1,0)*x - 304.D0/3.D0*Hr2(-1,0)
     &    *x**2 - 304.D0/3.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 + 1264.D0/9.D0*
     &    Hr2(0,0) + 196.D0/3.D0*Hr2(0,0)*x + 560.D0/9.D0*Hr2(0,0)*x**2
     &     + 2176.D0/9.D0*Hr2(0,0)*dx - 160.D0/9.D0*Hr2(0,0)*dp - 224.D0
     &    /9.D0*Hr2(0,0)*dm - 24.D0*Hr2(0,0)*z2 - 40.D0*Hr2(0,0)*z2*x
     &     + 44.D0/3.D0*Hr2(0,1) + 460.D0/3.D0*Hr2(0,1)*x - 176.D0/3.D0
     &    *Hr2(0,1)*x**2 + 176.D0/3.D0*Hr2(0,1)*dx + 320.D0/9.D0*Hr2(0,
     &    1)*dm - 48.D0*Hr2(0,1)*z2 - 48.D0*Hr2(0,1)*z2*x - 784.D0/9.D0
     &    *Hr2(1,0) + 464.D0/9.D0*Hr2(1,0)*x - 176.D0/3.D0*Hr2(1,0)*
     &    x**2 + 176.D0/3.D0*Hr2(1,0)*dx + 320.D0/9.D0*Hr2(1,0)*dm - 
     &    112.D0*Hr3(-1,-1,0) - 112.D0*Hr3(-1,-1,0)*x + 32.D0*Hr3(-1,-1
     &    ,0)*x**2 + 32.D0*Hr3(-1,-1,0)*dx - 40.D0/3.D0*Hr3(-1,0,0) + 
     &    88.D0/3.D0*Hr3(-1,0,0)*x - 160.D0/3.D0*Hr3(-1,0,0)*x**2 - 160.
     &    D0/3.D0*Hr3(-1,0,0)*dx )
      PTgg2 = PTgg2 + nf*ca**2 * ( 128.D0/3.D0*Hr3(-1,0,0)*dp + 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 + 104.D0*Hr3(0,-1,0) - 136.D0
     &    /3.D0*Hr3(0,-1,0)*x - 128.D0/3.D0*Hr3(0,-1,0)*x**2 - 32.D0*
     &    Hr3(0,-1,0)*dx + 64.D0/3.D0*Hr3(0,-1,0)*dp - 32.D0/3.D0*Hr3(0
     &    ,-1,0)*dm - 32.D0*Hr3(0,0,0) + 1000.D0/3.D0*Hr3(0,0,0)*x - 64.
     &    D0*Hr3(0,0,0)*x**2 + 128.D0*Hr3(0,0,0)*dx - 32.D0*Hr3(0,0,0)*
     &    dp + 96.D0*Hr3(0,0,0)*dm - 124.D0/3.D0*Hr3(0,0,1) + 44.D0/3.D0
     &    *Hr3(0,0,1)*x - 160.D0/3.D0*Hr3(0,0,1)*x**2 + 128.D0/3.D0*
     &    Hr3(0,0,1)*dx + 112.D0/3.D0*Hr3(0,0,1)*dm - 224.D0/3.D0*Hr3(0
     &    ,1,0) + 160.D0/3.D0*Hr3(0,1,0)*x - 128.D0/3.D0*Hr3(0,1,0)*
     &    x**2 + 128.D0/3.D0*Hr3(0,1,0)*dx + 128.D0/3.D0*Hr3(0,1,0)*dm
     &     - 140.D0/3.D0*Hr3(1,0,0) + 28.D0/3.D0*Hr3(1,0,0)*x - 128.D0/
     &    3.D0*Hr3(1,0,0)*x**2 + 128.D0/3.D0*Hr3(1,0,0)*dx + 112.D0/3.D0
     &    *Hr3(1,0,0)*dm - 96.D0*Hr4(0,-1,-1,0) + 96.D0*Hr4(0,-1,-1,0)*
     &    x )
      PTgg2 = PTgg2 + nf*ca**2 * ( 48.D0*Hr4(0,-1,0,0) - 48.D0*Hr4(0,-1
     &    ,0,0)*x + 80.D0*Hr4(0,0,-1,0) - 16.D0*Hr4(0,0,-1,0)*x + 16.D0
     &    *Hr4(0,0,0,0)*x + 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 )
      PTgg2 = PTgg2 + nf*cf*ca * ( 16180.D0/27.D0 - 12022.D0/27.D0*x + 
     &    50210.D0/81.D0*x**2 - 59714.D0/81.D0*dx - 110.D0/3.D0*dm - 
     &    776.D0/3.D0*z3 + 352.D0/3.D0*z3*x - 640.D0/3.D0*z3*x**2 + 512.
     &    D0/3.D0*z3*dx + 32.D0*z3*dm - 1432.D0/9.D0*z2 - 7972.D0/9.D0*
     &    z2*x - 352.D0/9.D0*z2*x**2 + 608.D0/9.D0*z2*dx + 616.D0/5.D0*
     &    z2**2 + 936.D0/5.D0*z2**2*x + 208.D0*Hr1(-1)*z2 + 208.D0*Hr1(
     &    -1)*z2*x - 64.D0/3.D0*Hr1(-1)*z2*x**2 - 64.D0/3.D0*Hr1(-1)*z2
     &    *dx - 5636.D0/27.D0*Hr1(0) - 25544.D0/27.D0*Hr1(0)*x - 500.D0
     &    *Hr1(0)*x**2 - 14552.D0/27.D0*Hr1(0)*dx - 16.D0*Hr1(0)*dm + 
     &    160.D0*Hr1(0)*z3 + 512.D0*Hr1(0)*z3*x - 304.D0/3.D0*Hr1(0)*z2
     &     + 8.D0/3.D0*Hr1(0)*z2*x - 32.D0/3.D0*Hr1(0)*z2*x**2 + 128.D0/
     &    3.D0*Hr1(0)*z2*dx + 2080.D0/9.D0*Hr1(1) - 76.D0/9.D0*Hr1(1)*x
     &     - 8068.D0/27.D0*Hr1(1)*x**2 + 2056.D0/27.D0*Hr1(1)*dx + 224.D
     &    0*Hr1(1)*z2 - 224.D0*Hr1(1)*z2*x - 448.D0/3.D0*Hr1(1)*z2*x**2
     &     + 448.D0/3.D0*Hr1(1)*z2*dx - 1664.D0/3.D0*Hr2(-1,0) - 1376.D0
     &    /3.D0*Hr2(-1,0)*x )
      PTgg2 = PTgg2 + nf*cf*ca * ( 2528.D0/9.D0*Hr2(-1,0)*x**2 + 1664.D0
     &    /9.D0*Hr2(-1,0)*dx + 128.D0*Hr2(0,-1)*z2 - 128.D0*Hr2(0,-1)*
     &    z2*x + 904.D0/3.D0*Hr2(0,0) + 2000.D0/3.D0*Hr2(0,0)*x + 1544.D
     &    0/9.D0*Hr2(0,0)*x**2 - 3712.D0/9.D0*Hr2(0,0)*dx - 16.D0*Hr2(0
     &    ,0)*z2 + 112.D0*Hr2(0,0)*z2*x + 1672.D0/9.D0*Hr2(0,1) + 2740.D
     &    0/9.D0*Hr2(0,1)*x + 272.D0*Hr2(0,1)*x**2 - 176.D0/9.D0*Hr2(0,
     &    1)*dx + 288.D0*Hr2(0,1)*z2 + 288.D0*Hr2(0,1)*z2*x + 148.D0/3.D
     &    0*Hr2(1,0) - 292.D0/3.D0*Hr2(1,0)*x + 1000.D0/9.D0*Hr2(1,0)*
     &    x**2 - 568.D0/9.D0*Hr2(1,0)*dx + 12.D0*Hr2(1,1) - 12.D0*Hr2(1
     &    ,1)*x - 16.D0/9.D0*Hr2(1,1)*x**2 + 16.D0/9.D0*Hr2(1,1)*dx + 
     &    288.D0*Hr3(-1,-1,0) + 288.D0*Hr3(-1,-1,0)*x - 256.D0/3.D0*
     &    Hr3(-1,-1,0)*x**2 - 256.D0/3.D0*Hr3(-1,-1,0)*dx - 208.D0*Hr3(
     &    -1,0,0) - 208.D0*Hr3(-1,0,0)*x + 224.D0/3.D0*Hr3(-1,0,0)*x**2
     &     + 224.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 )
      PTgg2 = PTgg2 + nf*cf*ca * (  - 416.D0*Hr3(0,-1,0) + 160.D0*Hr3(0
     &    ,-1,0)*x + 256.D0/3.D0*Hr3(0,-1,0)*dx + 640.D0/3.D0*Hr3(0,0,0
     &    ) - 920.D0/3.D0*Hr3(0,0,0)*x - 640.D0/3.D0*Hr3(0,0,0)*dx - 
     &    128.D0*Hr3(0,0,1) - 88.D0*Hr3(0,0,1)*x - 32.D0*Hr3(0,0,1)*
     &    x**2 - 64.D0/3.D0*Hr3(0,0,1)*dx + 64.D0/3.D0*Hr3(0,1,0) - 104.
     &    D0/3.D0*Hr3(0,1,0)*x - 160.D0/3.D0*Hr3(0,1,0)*x**2 + 136.D0/3.
     &    D0*Hr3(0,1,1) + 88.D0/3.D0*Hr3(0,1,1)*x - 32.D0/3.D0*Hr3(0,1,
     &    1)*x**2 + 32.D0/3.D0*Hr3(0,1,1)*dx - 200.D0*Hr3(1,0,0) + 200.D
     &    0*Hr3(1,0,0)*x + 512.D0/3.D0*Hr3(1,0,0)*x**2 - 512.D0/3.D0*
     &    Hr3(1,0,0)*dx - 48.D0*Hr3(1,0,1) + 48.D0*Hr3(1,0,1)*x + 64.D0
     &    *Hr3(1,0,1)*x**2 - 64.D0*Hr3(1,0,1)*dx - 32.D0*Hr3(1,1,0) + 
     &    32.D0*Hr3(1,1,0)*x + 128.D0/3.D0*Hr3(1,1,0)*x**2 - 128.D0/3.D0
     &    *Hr3(1,1,0)*dx - 16.D0*Hr3(1,1,1) + 16.D0*Hr3(1,1,1)*x + 64.D0
     &    /3.D0*Hr3(1,1,1)*x**2 - 64.D0/3.D0*Hr3(1,1,1)*dx + 256.D0*
     &    Hr4(0,-1,-1,0) - 256.D0*Hr4(0,-1,-1,0)*x - 192.D0*Hr4(0,-1,0,
     &    0) )
      PTgg2 = PTgg2 + nf*cf*ca * ( 192.D0*Hr4(0,-1,0,0)*x - 256.D0*Hr4(
     &    0,0,-1,0) + 128.D0*Hr4(0,0,-1,0)*x + 128.D0*Hr4(0,0,0,0) - 
     &    448.D0*Hr4(0,0,0,0)*x - 112.D0*Hr4(0,0,0,1) - 112.D0*Hr4(0,0,
     &    0,1)*x + 80.D0*Hr4(0,0,1,0) + 80.D0*Hr4(0,0,1,0)*x + 48.D0*
     &    Hr4(0,0,1,1) + 48.D0*Hr4(0,0,1,1)*x - 304.D0*Hr4(0,1,0,0) - 
     &    304.D0*Hr4(0,1,0,0)*x - 96.D0*Hr4(0,1,0,1) - 96.D0*Hr4(0,1,0,
     &    1)*x - 64.D0*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 )
      PTgg2 = PTgg2 + nf*cf**2 * (  - 1682.D0/9.D0 - 382.D0/9.D0*x + 
     &    292.D0/3.D0*x**2 + 132.D0*dx - 288.D0*z3 - 224.D0*z3*x - 256.D
     &    0/3.D0*z3*x**2 - 320.D0/3.D0*z3*dx + 824.D0/3.D0*z2 + 2516.D0/
     &    3.D0*z2*x + 320.D0*z2*x**2 - 152.D0*z2**2 - 184.D0*z2**2*x - 
     &    64.D0*Hr1(-1)*z2 - 64.D0*Hr1(-1)*z2*x + 64.D0/3.D0*Hr1(-1)*z2
     &    *x**2 + 64.D0/3.D0*Hr1(-1)*z2*dx + 370.D0/3.D0*Hr1(0) + 1786.D
     &    0/3.D0*Hr1(0)*x + 3604.D0/9.D0*Hr1(0)*x**2 + 32.D0*Hr1(0)*dx
     &     - 160.D0*Hr1(0)*z3 - 288.D0*Hr1(0)*z3*x - 80.D0*Hr1(0)*z2*x
     &     - 128.D0/3.D0*Hr1(0)*z2*x**2 - 1432.D0/3.D0*Hr1(1) + 716.D0/
     &    3.D0*Hr1(1)*x + 3316.D0/9.D0*Hr1(1)*x**2 - 1168.D0/9.D0*Hr1(1
     &    )*dx - 144.D0*Hr1(1)*z2 + 144.D0*Hr1(1)*z2*x + 128.D0*Hr1(1)*
     &    z2*x**2 - 128.D0*Hr1(1)*z2*dx + 1072.D0/3.D0*Hr2(-1,0) + 1072.
     &    D0/3.D0*Hr2(-1,0)*x - 64.D0*Hr2(0,-1)*z2 + 64.D0*Hr2(0,-1)*z2
     &    *x - 760.D0/3.D0*Hr2(0,0) - 784.D0*Hr2(0,0)*x - 160.D0/3.D0*
     &    Hr2(0,0)*x**2 + 16.D0*Hr2(0,0)*z2 + 16.D0*Hr2(0,0)*z2*x - 296.
     &    D0*Hr2(0,1) )
      PTgg2 = PTgg2 + nf*cf**2 * (  - 556.D0*Hr2(0,1)*x - 512.D0/3.D0*
     &    Hr2(0,1)*x**2 - 160.D0/3.D0*Hr2(0,1)*dx - 224.D0*Hr2(0,1)*z2
     &     - 224.D0*Hr2(0,1)*z2*x - 116.D0/3.D0*Hr2(1,0) - 28.D0/3.D0*
     &    Hr2(1,0)*x + 416.D0/3.D0*Hr2(1,0)*x**2 - 272.D0/3.D0*Hr2(1,0)
     &    *dx - 28.D0/3.D0*Hr2(1,1) + 28.D0/3.D0*Hr2(1,1)*x + 64.D0/3.D0
     &    *Hr2(1,1)*x**2 - 64.D0/3.D0*Hr2(1,1)*dx - 128.D0*Hr3(-1,-1,0)
     &     - 128.D0*Hr3(-1,-1,0)*x + 128.D0/3.D0*Hr3(-1,-1,0)*x**2 + 
     &    128.D0/3.D0*Hr3(-1,-1,0)*dx + 64.D0*Hr3(-1,0,0) + 64.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 + 128.D0*Hr3(0,-1,0) - 128.D0/3.D0*Hr3(0,-1,0)*x**2 - 
     &    72.D0*Hr3(0,0,0) - 40.D0*Hr3(0,0,0)*x - 160.D0/3.D0*Hr3(0,0,0
     &    )*x**2 + 16.D0*Hr3(0,0,1)*x - 128.D0/3.D0*Hr3(0,0,1)*dx - 192.
     &    D0*Hr3(0,1,0) - 176.D0*Hr3(0,1,0)*x - 64.D0/3.D0*Hr3(0,1,0)*
     &    x**2 - 256.D0/3.D0*Hr3(0,1,0)*dx - 64.D0*Hr3(0,1,1) - 48.D0*
     &    Hr3(0,1,1)*x - 64.D0/3.D0*Hr3(0,1,1)*dx + 152.D0*Hr3(1,0,0)
     &     - 152.D0*Hr3(1,0,0)*x )
      PTgg2 = PTgg2 + nf*cf**2 * (  - 416.D0/3.D0*Hr3(1,0,0)*x**2 + 416.
     &    D0/3.D0*Hr3(1,0,0)*dx + 48.D0*Hr3(1,0,1) - 48.D0*Hr3(1,0,1)*x
     &     - 64.D0*Hr3(1,0,1)*x**2 + 64.D0*Hr3(1,0,1)*dx + 32.D0*Hr3(1,
     &    1,0) - 32.D0*Hr3(1,1,0)*x - 128.D0/3.D0*Hr3(1,1,0)*x**2 + 128.
     &    D0/3.D0*Hr3(1,1,0)*dx + 16.D0*Hr3(1,1,1) - 16.D0*Hr3(1,1,1)*x
     &     - 64.D0/3.D0*Hr3(1,1,1)*x**2 + 64.D0/3.D0*Hr3(1,1,1)*dx - 
     &    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) - 32.D
     &    0*Hr4(0,0,0,0) - 32.D0*Hr4(0,0,0,0)*x - 16.D0*Hr4(0,0,0,1) - 
     &    16.D0*Hr4(0,0,0,1)*x - 144.D0*Hr4(0,0,1,0) - 144.D0*Hr4(0,0,1
     &    ,0)*x - 48.D0*Hr4(0,0,1,1) - 48.D0*Hr4(0,0,1,1)*x + 240.D0*
     &    Hr4(0,1,0,0) + 240.D0*Hr4(0,1,0,0)*x + 96.D0*Hr4(0,1,0,1) + 
     &    96.D0*Hr4(0,1,0,1)*x + 64.D0*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 )
      PTgg2 = PTgg2 + nf2*ca * (  - 110.D0/27.D0 + 14.D0/3.D0*x - 472.
     &    D0/81.D0*x**2 + 472.D0/81.D0*dx - 16.D0/27.D0*dm + 32.D0/9.D0
     &    *z2 + 32.D0/9.D0*z2*x - 104.D0/9.D0*Hr1(0) + 130.D0/9.D0*Hr1(
     &    0)*x - 88.D0/9.D0*Hr1(0)*x**2 + 368.D0/27.D0*Hr1(0)*dx + 160.D
     &    0/27.D0*Hr1(0)*dm + 22.D0/9.D0*Hr1(1) - 22.D0/9.D0*Hr1(1)*x
     &     + 104.D0/27.D0*Hr1(1)*x**2 - 104.D0/27.D0*Hr1(1)*dx - 80.D0/
     &    9.D0*Hr2(0,0) + 112.D0/9.D0*Hr2(0,0)*x - 64.D0/9.D0*Hr2(0,0)*
     &    x**2 + 64.D0/9.D0*Hr2(0,0)*dx + 64.D0/9.D0*Hr2(0,0)*dm - 32.D0
     &    /9.D0*Hr2(0,1) - 32.D0/9.D0*Hr2(0,1)*x )
      PTgg2 = PTgg2 + nf2*cf * ( 544.D0/9.D0 - 304.D0/9.D0*x - 1216.D0
     &    /81.D0*x**2 - 944.D0/81.D0*dx - 64.D0/3.D0*z3 - 64.D0/3.D0*z3
     &    *x - 200.D0/9.D0*z2 - 440.D0/9.D0*z2*x - 224.D0/9.D0*z2*x**2
     &     - 208.D0/27.D0*Hr1(0) - 64.D0/27.D0*Hr1(0)*x + 16.D0*Hr1(0)*
     &    x**2 - 736.D0/27.D0*Hr1(0)*dx + 16.D0/3.D0*Hr1(0)*z2 + 16.D0/
     &    3.D0*Hr1(0)*z2*x - 32.D0/3.D0*Hr1(1) + 80.D0/3.D0*Hr1(1)*x - 
     &    1168.D0/27.D0*Hr1(1)*x**2 + 736.D0/27.D0*Hr1(1)*dx - 160.D0/3.
     &    D0*Hr2(0,0) - 64.D0*Hr2(0,0)*x - 160.D0/9.D0*Hr2(0,0)*x**2 - 
     &    128.D0/9.D0*Hr2(0,0)*dx + 296.D0/9.D0*Hr2(0,1) + 344.D0/9.D0*
     &    Hr2(0,1)*x + 32.D0/3.D0*Hr2(0,1)*x**2 + 128.D0/9.D0*Hr2(0,1)*
     &    dx - 8.D0/3.D0*Hr2(1,0) + 8.D0/3.D0*Hr2(1,0)*x + 32.D0/9.D0*
     &    Hr2(1,0)*x**2 - 32.D0/9.D0*Hr2(1,0)*dx - 8.D0/3.D0*Hr2(1,1)
     &     + 8.D0/3.D0*Hr2(1,1)*x + 32.D0/9.D0*Hr2(1,1)*x**2 - 32.D0/9.D
     &    0*Hr2(1,1)*dx - 64.D0/3.D0*Hr3(0,0,0) - 64.D0/3.D0*Hr3(0,0,0)
     &    *x + 16.D0*Hr3(0,0,1) + 16.D0*Hr3(0,0,1)*x - 16.D0/3.D0*Hr3(0
     &    ,1,0) )
      PTgg2 = PTgg2 + nf2*cf * (  - 16.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 + 88.D0/3.D0*z3 - 1072.D0/9.D0*z2
     ,                   + 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
*
       GGG2L = DM * A3G
*
* ...The regular piece of the splitting function
*
       X2GGTA = PTgg2 - GGG2L
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The singular (soft) piece of P_gg^(2)T  (as in the spacelike case)
*
       FUNCTION X2GGTB (Y, NF)
       IMPLICIT REAL*8 (A - Z)
       INTEGER NF
*
       COMMON / P2GSOFT / A3G
*
       X2GGTB  = A3G/(1.D0-Y)
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..The 'local' piece of P_gg^(2)T  (as in the spacelike case)
*
       FUNCTION X2GGTC (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,
     ,             Z5 = 1.0369 27755 14336 99263 D0 )
*
       COMMON / P2GSOFT / A3G
*
* ...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
*
       X2GGTC = LOG (1.D0-Y) * A3G + P2GDELT
*
       RETURN
       END
*
* =================================================================av==
