*
* ..File: xpij2n.f   UPDATE 7/2000
*
*
* ..Parametrisation of the 3-loop MS(bar) splitting functions P_ij^(2) 
*    for the evolution of unpolarized singlet partons at mu_r = mu_f.
*    The expansion parameter is alpha_s/(4 pi).
*
* ..The two sets providing the error estimate are called via  IMOD = 1
*    and  IMOD = 2.  Any other value of IMOD invokes their average.
*
* ..The distributions (in the mathematical sense) are given as in eq.
*   (B.27) 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 results are based on the lowest six even-integer moments 
*    calculated by Larin et al. and  Retey and Vermaseren, supplemented
*    by the leading small-x terms obtained by Catani and Hautmann, 
*    Fadin and Lipatov, and Camici and Ciafaloni. 
*
* ..Reference: W.L. van Neerven and A. Vogt, 
*              hep-ph/9907472  (non-singlet, needed for P_qq),
*              hep-ph/0006154, and  hep-ph/0007362
*    It is appropriate to cite also the sources of information, i.e, 
*    refs. [2-9] of hep-ph/0007362.
*
* 
* =====================================================================
*
*
* ..This is the (regular) pure-singlet splitting functions P_ps^(2).
*    P_qq^(2) is obtained by adding the non-singlet quantity P_NS^(2)+.
*
       FUNCTION P2PSA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       DIMENSION CL3(2)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       G2PSA1 = (229.497 * DL1 + 722.99 * Y**2 - 2678.77 + 560.20/Y)
     1          * (1.-Y) - 2008.61 * DL - 998.15 * DL**2
       G2PSA2 = (-73.845 * DL1**2 - 305.988 * DL1 - 2063.19 * Y
     1          + 387.95/Y ) * (1.-Y) - 1999.35 * Y*DL + 732.68 * DL  
*
       IF (IMOD .EQ. 1) THEN
         G2PSA = G2PSA1
       ELSE IF (IMOD .EQ. 2) THEN
         G2PSA = G2PSA2
       ELSE
         G2PSA = 0.5 * (G2PSA1 + G2PSA2)
       END IF
*
       G2PSA = G2PSA + 3584./(27.* Y) * DL 
     1         + NF * ( (7.282 * DL1 + 38.779 * Y**2 - 32.022 * Y 
     2                  + 6.252 - 1.767 / Y ) * (1.-Y) -7.453 * DL**2 ) 
       P2PSA = - NF * G2PSA
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the gluon->quark splitting functions P_qg^(2).
*
       FUNCTION P2QGA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       DIMENSION CL3(2)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       G2QGA1 =   31.830 * DL1**3 - 1252.267 * DL1 - 1722.47
     1          - 1999.89 * Y - 1223.43 * DL**2 + 1334.61 / Y
       G2QGA2 = - 19.428 * DL1**4 - 159.833 * DL1**3 - 309.384 * DL1**2
     1          - 2631.00 * (1.-Y) + 67.25 * DL**2 + 776.793 / Y
*
       IF (IMOD .EQ. 1) THEN
         G2QGA = G2QGA1
       ELSE IF (IMOD .EQ. 2) THEN
         G2QGA = G2QGA2
       ELSE
         G2QGA = 0.5 * (G2QGA1 + G2QGA2)
       END IF
*
       G2QGA = G2QGA + 896./ (3.*Y) * DL 
     1         + NF * ( 0.9085 * DL1**2 + 35.803 * DL1 + 128.023
     2                  - 200.929 * (1.-Y) - 40.542 * DL - 3.284 / Y)
       P2QGA = - NF * G2QGA
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the quark->gluon splitting functions P_gq^(2).
*
       FUNCTION P2GQA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       G2GQA1 = -13.1212 * DL1**4 - 126.665 * DL1**3 - 308.536 * DL1**2
     1          - 361.21 + 2113.45 * DL + 17.965 * DL/Y
     2   + NF* (- 2.4427 * DL1**4 - 27.763 * DL1**3 - 80.548 * DL1**2 
     3          +227.135 + 151.04 * DL**2 - 65.91  * DL/Y )
       G2GQA2 = 4.5107 * DL1**4 + 66.618 * DL1**3 + 231.535 * DL1**2
     1          + 1224.22* (1.-Y) - 240.08 * DL**2 - 379.60 * (DL+4.)/Y
     2   + NF* (1.4028 * DL1**4 + 11.638 * DL1**3 - 164.963 * DL1
     3          + 1066.78* (1.-Y) + 182.08 * DL**2 - 138.54 * (DL+2.)/Y)
*
       IF (IMOD .EQ. 1) THEN
         G2GQA = G2GQA1
       ELSE IF (IMOD .EQ. 2) THEN
         G2GQA = G2GQA2
       ELSE
         G2GQA = 0.5 * (G2GQA1 + G2GQA2)
       END IF
*
       G2GQA = G2GQA 
     1         + NF**2* (- 1.9361 * DL1**2 - 11.178 * DL1 - 11.632 
     2                   + 15.145 * (1.-Y) - 3.354 * DL + 2.133 / Y )
       P2GQA = - G2GQA
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the regular piece of the gg splitting function P_gg^(2).
*
       FUNCTION P2GGA (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       INTEGER IMOD, NF
       DL  = LOG (Y)
       DL1 = LOG (1.-Y)
*
       G2GGA1 = 732.715 * DL1**2 + 20640.069 * Y + 15428.58 * (1.-Y**2) 
     1          + 15213.6 * DL**2 - 16700.88 * 1./Y
     2   + NF* (425.708 * DL1 - 914.548  * Y**2 + 1122.86
     3          + 444.21 * DL**2 - 376.98 * 1./Y )
       G2GGA2 = - 3748.934 * DL1 + 35974.45 * (1.+Y**2) - 60879.62 * Y
     2          - 2002.96 * DL**2 - 9762.09 * 1./Y
     3   + NF* (- 62.630 * DL1**2 - 801.90 - 1891.40 * DL
     4          - 813.78 * DL**2 - 1.360 * 1./Y ) 
*
       IF (IMOD .EQ. 1) THEN
         G2GGA = G2GGA1
       ELSE IF (IMOD .EQ. 2) THEN
         G2GGA = G2GGA2
       ELSE
         G2GGA = 0.5 * (G2GGA1 + G2GGA2)
       END IF
*
       G2GGA = G2GGA - (2675.85 + NF * 157.18) * DL / Y 
     1         + NF**2 * (- 37.6417 * Y**2 + 72.926 * Y - 32.349
     2                    + 0.991 * DL**2 - 2.818 / Y )
       P2GGA = - G2GGA
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the singular piece of the gg splitting function P_gg^(2).
*
       FUNCTION P2GGB (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       INTEGER IMOD, NF
       D1 = 1./(1.-Y)
*
       G2GGB1 = (-2626.38 + NF * 415.71) * D1
       G2GGB2 = (-2678.22 + NF * 412.00) * D1
*
       IF (IMOD .EQ. 1) THEN
         G2GGB = G2GGB1
       ELSE IF (IMOD .EQ. 2) THEN
         G2GGB = G2GGB2
       ELSE
         G2GGB = 0.5 * (G2GGB1 + G2GGB2)
       END IF
*
       G2GGB =  G2GGB + NF**2 * 16./9.D0 * D1 
       P2GGB = -G2GGB
*
       RETURN
       END
*
* ---------------------------------------------------------------------
*
*
* ..This is the 'local' piece of the gg splitting function P_gg^(2).
*
       FUNCTION P2GGC (Y, NF, IMOD)
*
       IMPLICIT REAL*8 (A-Z)
       INTEGER IMOD, NF
       DL1 = LOG (1.-Y)
*
       G2GGC1 = (-2626.38 + NF * 415.71)* DL1 - 4424.168 + NF * 548.569
       G2GGC2 = (-2678.22 + NF * 412.00)* DL1 - 4590.570 + NF * 534.951
*
       IF (IMOD .EQ. 1) THEN
         G2GGC = G2GGC1
       ELSE IF (IMOD .EQ. 2) THEN
         G2GGC = G2GGC2
       ELSE
         G2GGC = 0.5 * (G2GGC1 + G2GGC2)
       END IF
*
       G2GGC =  G2GGC + NF**2 * (16./9. * DL1 - 6.4882)
       P2GGC = -G2GGC
*
       RETURN
       END
*
* =================================================================av==
