*
* ..File: c2lmom.f
*
*
* ..The subroutines c21p, c22p, c23p, cl1p, cl2p and cl3p
*    for the (complex) moments of the parametrized non-singlet 
*    and singlet MS(bar) coefficient functions for the unpolarized 
*    deep-inelastic structure functions F_2 and F_L, mu_r = mu_f.
*
* ..The QCD colour factors have been hard-wired in the parametrizations.
*    The coupling constant is normalized as  a_s = alpha_s/(4*pi). 
*    
* ..The required routines for simple harmonic sums (expressed through 
*    Euler's psi function and its derivatives) are appended at the end 
*    of the file.
* 
* ..References: hep-ph/0411112 = Phys. Lett. B606 (2005) 123 
*               hep-ph/0504242 = Nucl. Phys. B724 (2005) 3 
*
*
* =====================================================================
c--------------------
      subroutine c21p(cn,nf,c2pls1n,c2g1n)
      IMPLICIT complex*16 (A-H,O-Z)

      integer nf
      real*8 dl,pi,ge,z2,z3,z4,z5,z6

      pi = 3.1415926535897932385d0
      dl = dlog(2.0d0)
      ge = 5.7721566490153286061d-1
      z2 = 1.6449340668482264365d0
      z3 = 1.2020569031595942854d0
      z4 = 1.0823232337111381916d0
      z5 = 1.0369277551433699263d0
      z6 = 1.0173430619844491398d0

      c2pls1n=+4.d0/3.d0
     &  *(-9.d0-2.d0*(cn+1)**(-2)+4.d0*(cn+1)**(-1)+2.d0*(cn)**(-2)+3.
     &  d0*(cn)**(-1)+2.d0*s1(1,1.d0+cn)*(cn+1)**(-1)+3.d0*s1(1,cn)-2.d0
     &    *s1(1,cn)*(cn)**(-1)+2.d0*s1(1,cn)**2-2.d0*s1(2,cn))

      c2g1n=
     & +nf*(4.d0*(cn+2)**(-2)-16.d0*(cn+2)**(-1)-4.d0*(cn+1)**(-2)+16.d0
     &    *(cn+1)**(-1)+2.d0*(cn)**(-2)-2.d0*(cn)**(-1)+4.d0*s1(1,1.d0+
     &    cn)*(cn+1)**(-1)-4.d0*s1(1,2.d0+cn)*(cn+2)**(-1)-2.d0*s1(1,cn
     &    )*(cn)**(-1))

      return
      end
* =====================================================================
c--------------------
      subroutine c22p(cn,nf,c2pls2n,c2ps2n,c2g2n)
      IMPLICIT complex*16 (A-H,O-Z)

      integer nf
      real*8 dl,pi,ge,z2,z3,z4,z5,z6

      pi = 3.1415926535897932385d0
      dl = dlog(2.0d0)
      ge = 5.7721566490153286061d-1
      z2 = 1.6449340668482264365d0
      z3 = 1.2020569031595942854d0
      z4 = 1.0823232337111381916d0
      z5 = 1.0369277551433699263d0
      z6 = 1.0173430619844491398d0

      c2pls2n = 
     &     + nf * ( + 216883.d0/4629.d0 - 8113.d0/1000.d0*z2*cn**(-1)
     &     + 40.d0/9.d0*cn**(-3) - 16.d0/3.d0*cn**(-2)
     &     - 73101.d0/50000.d0*cn**(-1) + 111.d0/100.d0*(cn+1.d0)**(-4)
     &     - 891.d0/50.d0*(cn+1.d0)**(-1)
     &     - 1297.d0/100.d0*(cn+2.d0)**(-1) - 79361.d0/12500.d0*s1(1,cn)
     &     + 8113.d0/1000.d0*s1(1,cn)*cn**(-2)
     &     - 43949.d0/2700.d0*s1(1,cn)*cn**(-1)
     &     - 116.d0/27.d0*s1(1,cn)**2 + 5.d0/18.d0*s1(1,cn)**2*cn**(-1)
     &     - 16.d0/27.d0*s1(1,cn)**3 - 16.d0/9.d0*s1(1,cn)*s1(2,cn)
     &     - 116.d0/27.d0*s1(2,cn) + 75517.d0/9000.d0*s1(2,cn)*cn**(-1)
     &     - 32.d0/27.d0*s1(3,cn) )
      c2pls2n = c2pls2n
     &     - 338513.d0/1000.d0 + 151.d0/2.d0*z2*cn**(-2)
     &     - 2187.d0/10.d0*z3*cn**(-1) + 160.d0/9.d0*cn**(-4)
     &     - 207.d0/5.d0*cn**(-3) + 3548.d0/125.d0*cn**(-2)
     &     + 7641.d0/1000.d0*cn**(-1) + 2157.d0/125.d0*(cn+1.d0)**(-5)
     &     - 8067.d0/10.d0*(cn+1.d0)**(-1) - 188641.d0/1000.d0*s1(1,cn)
     &     - 1471.d0/5.d0*s1(1,cn)*z2*cn**(-1)
     &     - 151.d0/2.d0*s1(1,cn)*cn**(-3)
     &     + 1649763.d0/2500.d0*s1(1,cn)*cn**(-1)
     &     - 77763.d0/5000.d0*s1(1,cn)**2
     &     + 1471.d0/10.d0*s1(1,cn)**2*cn**(-2)
      c2pls2n = c2pls2n
     &     + 818.d0/75.d0*s1(1,cn)**2*cn**(-1) + 184.d0/9.d0*s1(1,cn)**3
     &     + 1583.d0/450.d0*s1(1,cn)**3*cn**(-1)
     &     + 32.d0/9.d0*s1(1,cn)**4 + 64.d0/3.d0*s1(1,cn)**2*s1(2,cn)
     &     + 184.d0/3.d0*s1(1,cn)*s1(2,cn)
     &     + 45713.d0/150.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 256.d0/9.d0*s1(1,cn)*s1(3,cn) - 77763.d0/5000.d0*s1(2,cn)
     &     + 358.d0/5.d0*s1(2,cn)*cn**(-2)
     &     + 818.d0/75.d0*s1(2,cn)*cn**(-1) + 32.d0/3.d0*s1(2,cn)**2
     &     + 368.d0/9.d0*s1(3,cn) + 101581.d0/450.d0*s1(3,cn)*cn**(-1)
     &     + 64.d0/3.d0*s1(4,cn)

      c2ps2n =
     &     + nf * ( - 1013.d0/25.d0*z2*cn**(-2)
     &     + 3393.d0/50.d0*z2*(cn+1.d0)**(-2)
     &     - 1013.d0/25.d0*z3*cn**(-1)
     &     + 3393.d0/50.d0*z3*(cn+1.d0)**(-1) - 80.d0/3.d0*cn**(-4)
     &     - 1053.d0/500.d0*cn**(-3) - 1084.d0/25.d0*cn**(-2)
     &     + 88831.d0/10000.d0*cn**(-1)
     &     + 52903.d0/10000.d0*(cn-1.d0)**(-1)
     &     - 17001.d0/500.d0*(cn+1.d0)**(-4)
     &     + 1084.d0/25.d0*(cn+1.d0)**(-2)
     &     - 185767.d0/5000.d0*(cn+1.d0)**(-1)
     &     + 3639.d0/100.d0*(cn+2.d0)**(-1)
     &     - 1341.d0/100.d0*(cn+3.d0)**(-1) )
      c2ps2n = c2ps2n
     &     + nf * ( - 3393.d0/50.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-3)
     &     - 32.d0/3.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 8.d0/3.d0*s1(1,cn+1.d0)**2*(cn+1.d0)**(-1)
     &     + 4.d0/25.d0*s1(1,cn+1.d0)**3*(cn+1.d0)**(-1)
     &     + 12.d0/25.d0*s1(1,cn+1.d0)*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 2.d0/25.d0*s1(1,cn+2.d0)**3*(cn+2.d0)**(-1)
     &     - 6.d0/25.d0*s1(1,cn+2.d0)*s1(2,cn+2.d0)*(cn+2.d0)**(-1)
     &     + 1013.d0/25.d0*s1(1,cn)*cn**(-3)
     &     + 32.d0/3.d0*s1(1,cn)*cn**(-1)
     &     + 8.d0/3.d0*s1(1,cn)**2*cn**(-1)
     &     - 2.d0/25.d0*s1(1,cn)**3*cn**(-1)
     &     - 6.d0/25.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     - 3393.d0/50.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-2)
     &     - 8.d0/3.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-1) )
      c2ps2n = c2ps2n
     &     + nf * ( + 1013.d0/25.d0*s1(2,cn)*cn**(-2)
     &     + 8.d0/3.d0*s1(2,cn)*cn**(-1)
     &     - 3377.d0/50.d0*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 4.d0/25.d0*s1(3,cn+2.d0)*(cn+2.d0)**(-1)
     &     + 1009.d0/25.d0*s1(3,cn)*cn**(-1) )

      c2g2n = 
     &     + nf * ( + 2086.d0*z2*cn**(-2) - 1469.d0/5.d0*z2*cn**(-1)
     &     - 3556.d0/5.d0*z2*(cn+1.d0)**(-1) + 2086.d0*z3*cn**(-1)
     &     - 140.d0/3.d0*cn**(-4) - 7109.d0/500.d0*cn**(-3)
     &     - 578.d0/5.d0*cn**(-2) - 63973.d0/10000.d0*cn**(-1)
     &     + 119033.d0/10000.d0*(cn-1.d0)**(-1)
     &     + 2408.d0*(cn+1.d0)**(-3) - 36761.d0/50.d0*(cn+1.d0)**(-1)
     &     + 7603.d0/10.d0*(cn+2.d0)**(-1)
     &     + 3556.d0/5.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-2)
     &     - 593.d0/20.d0*s1(1,cn+1.d0)**3*(cn+1.d0)**(-1) )
      c2g2n = c2g2n
     &     + nf * ( - 1779.d0/20.d0*s1(1,cn+1.d0)*s1(2,cn+1.d0)*
     &    (cn+1.d0)**(-1) - 2086.d0*s1(1,cn)*cn**(-3)
     &     + 1469.d0/5.d0*s1(1,cn)*cn**(-2)
     &     + 872.d0/25.d0*s1(1,cn)*cn**(-1) - 24.d0*s1(1,cn)**2*cn**(-1)
     &     + 4177.d0/180.d0*s1(1,cn)**3*cn**(-1)
     &     + 4177.d0/60.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 3556.d0/5.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 2086.d0*s1(2,cn)*cn**(-2) + 1349.d0/5.d0*s1(2,cn)*cn**(-1)
     &     - 593.d0/10.d0*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 183563.d0/90.d0*s1(3,cn)*cn**(-1) )

      return
      end
* =====================================================================
c--------------------
      subroutine c23p(cn,nf,c2pls3n,c2ps3n,c2g3n)
      IMPLICIT complex*16 (A-H,O-Z)

      dimension fl11ns(7),fl11ps(7),fl11g(7)
      data fl11ns 
     &     /0.d0, 2.d0, 0.5d0, 0.d0, 0.5d0, 0.2d0, 0.5d0 /
      data fl11g  
     &     /0.d0, 1.d0, 0.1d0, 0.d0, 0.1d0, 0.018181818d0, 0.1d0 /
      data fl11ps 
     &     /0.d0,-1.d0,-0.4d0, 0.d0,-0.4d0,-0.181818182d0,-0.4d0 /

      integer nf
      real*8 dl,pi,ge,z2,z3,z4,z5,z6

      pi = 3.1415926535897932385d0
      dl = dlog(2.0d0)
      ge = 5.7721566490153286061d-1
      z2 = 1.6449340668482264365d0
      z3 = 1.2020569031595942854d0
      z4 = 1.0823232337111381916d0
      z5 = 1.0369277551433699263d0
      z6 = 1.0173430619844491398d0

      c2pls3n=
     &     + fl11ns(nf)*nf * ( - 1486.d0/125.d0
     &     + 640.d0/27.d0*(cn+1.d0)**(-4)
     &     - 3200.d0/81.d0*(cn+1.d0)**(-3)
     &     - 5959.d0/100.d0*(cn+1.d0)**(-2)
     &     + 99703.d0/1000.d0*(cn+1.d0)**(-1)
     &     - 1842.d0/25.d0*(cn+2.d0)**(-5)
     &     + 10437.d0/50.d0*(cn+2.d0)**(-4)
     &     - 1018.d0/5.d0*(cn+2.d0)**(-3)
     &     - 17671.d0/100.d0*(cn+2.d0)**(-1)
     &     + 7.d0/50.d0*(cn+3.d0)**(-1) + 1003.d0/20.d0*(cn+4.d0)**(-1)
     &     + 363.d0/40.d0*s1(1,cn+2.d0)*(cn+2.d0)**(-1)
     &     - 363.d0/40.d0*s1(1,cn+3.d0)*(cn+3.d0)**(-1) )
      c2pls3n = c2pls3n
     &     + nf * ( + 2574687.d0/1000.d0
     &     + 7162.d0/5.d0*z2*cn**(-2) - 4365.d0*z2*cn**(-1)
     &     - 52668.d0/5.d0*z3*cn**(-1) + 5824.d0/81.d0*cn**(-5)
     &     - 24448.d0/81.d0*cn**(-4) + 1873.d0/5.d0*cn**(-3)
     &     - 1378.d0/5.d0*cn**(-2) + 102241.d0/1000.d0*cn**(-1)
     &     + 12306.d0/125.d0*(cn+1.d0)**(-5) - 6752.d0*(cn+1.d0)**(-1)
     &     - 2778.d0*(cn+2.d0)**(-1)
     &     - 171.d0*s1(1,cn+1.d0)**4*(cn+1.d0)**(-1)
     &     - 1026.d0*s1(1,cn+1.d0)**2*s1(2,cn+1.d0)*(cn+1.d0)**(-1) )
      c2pls3n = c2pls3n
     &     + nf * (- 1368.d0*s1(1,cn+1.d0)*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 729359.d0/1000.d0*s1(1,cn) - 11966.d0*s1(1,cn)*z2*cn**(-1)
     &     - 7162.d0/5.d0*s1(1,cn)*cn**(-3) + 4365.d0*s1(1,cn)*cn**(-2)
     &     + 103097.d0/500.d0*s1(1,cn)*cn**(-1)
     &     + 147453.d0/1000.d0*s1(1,cn)**2
     &     + 5983.d0*s1(1,cn)**2*cn**(-2)
     &     - 608127.d0/1000.d0*s1(1,cn)**2*cn**(-1)
     &     - 220573.d0/3000.d0*s1(1,cn)**3
     &     - 11683.d0/162.d0*s1(1,cn)**3*cn**(-1)
     &     - 1648.d0/81.d0*s1(1,cn)**4 + 171.d0*s1(1,cn)**4*cn**(-1)
     &     - 128.d0/81.d0*s1(1,cn)**5 )
      c2pls3n = c2pls3n
     &     + nf * ( - 1280.d0/81.d0*s1(1,cn)**3*s1(2,cn)
     &     - 3296.d0/27.d0*s1(1,cn)**2*s1(2,cn)
     &     + 1026.d0*s1(1,cn)**2*s1(2,cn)*cn**(-1)
     &     - 2560.d0/81.d0*s1(1,cn)**2*s1(3,cn)
     &     - 220573.d0/1000.d0*s1(1,cn)*s1(2,cn)
     &     + 634481.d0/54.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     - 640.d0/27.d0*s1(1,cn)*s1(2,cn)**2
     &     - 13184.d0/81.d0*s1(1,cn)*s1(3,cn)
     &     + 1368.d0*s1(1,cn)*s1(3,cn)*cn**(-1)
     &     - 1280.d0/27.d0*s1(1,cn)*s1(4,cn)
     &     - 513.d0*s1(2,cn+1.d0)**2*(cn+1.d0)**(-1)
     &     + 147453.d0/1000.d0*s1(2,cn)
     &     + 22753.d0/5.d0*s1(2,cn)*cn**(-2)
     &     + 3756873.d0/1000.d0*s1(2,cn)*cn**(-1) )
      c2pls3n = c2pls3n
     &     + nf * ( - 1648.d0/27.d0*s1(2,cn)**2
     &     + 513.d0*s1(2,cn)**2*cn**(-1)
     &     - 2560.d0/81.d0*s1(2,cn)*s1(3,cn)
     &     - 220573.d0/1500.d0*s1(3,cn)
     &     + 4207693.d0/405.d0*s1(3,cn)*cn**(-1)
     &     - 1026.d0*s1(4,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 3296.d0/27.d0*s1(4,cn)
     &     + 1026.d0*s1(4,cn)*cn**(-1)
     &     - 1024.d0/27.d0*s1(5,cn) )
      c2pls3n = c2pls3n
     &     + nf**2 * ( - 516183.d0/5000.d0 - 623.d0/25.d0*z2*cn**(-2)
     &     + 9607.d0/100.d0*z2*cn**(-1) + 3671.d0/25.d0*z3*cn**(-1)
     &     + 736.d0/81.d0*cn**(-4) - 3968.d0/243.d0*cn**(-3)
     &     + 4021.d0/500.d0*cn**(-2) + 10083.d0/10000.d0*cn**(-1)
     &     + 646.d0/5.d0*(cn+1.d0)**(-1) + 205.d0/2.d0*(cn+2.d0)**(-1)
     &     - 10083.d0/10000.d0*s1(1,cn)
     &     + 4294.d0/25.d0*s1(1,cn)*z2*cn**(-1)
     &     + 623.d0/25.d0*s1(1,cn)*cn**(-3)
     &     - 9607.d0/100.d0*s1(1,cn)*cn**(-2) )
      c2pls3n = c2pls3n
     &     + nf**2 * ( + 228299.d0/20000.d0*s1(1,cn)*cn**(-1)
     &     + 153501.d0/40000.d0*s1(1,cn)**2
     &     - 2147.d0/25.d0*s1(1,cn)**2*cn**(-2)
     &     + 101101.d0/8100.d0*s1(1,cn)**2*cn**(-1)
     &     + 464.d0/243.d0*s1(1,cn)**3 + 16.d0/81.d0*s1(1,cn)**4
     &     + 32.d0/27.d0*s1(1,cn)**2*s1(2,cn)
     &     + 464.d0/81.d0*s1(1,cn)*s1(2,cn)
     &     - 4294.d0/25.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 128.d0/81.d0*s1(1,cn)*s1(3,cn)
     &     + 153501.d0/40000.d0*s1(2,cn)
     &     - 1524.d0/25.d0*s1(2,cn)*cn**(-2)
     &     - 338533.d0/4050.d0*s1(2,cn)*cn**(-1)
     &     + 16.d0/27.d0*s1(2,cn)**2 + 928.d0/243.d0*s1(3,cn)
     &     - 3671.d0/25.d0*s1(3,cn)*cn**(-1) + 32.d0/27.d0*s1(4,cn) )
      c2pls3n = c2pls3n
     &     - 458519.d0/50.d0 - 3672.d0*z2*cn**(-2) + 6158.d0*z2*cn**(-1)
     &     - 3672.d0*z3*cn**(-1) + 1280.d0/9.d0*cn**(-6)
     &     - 23456.d0/27.d0*cn**(-5) + 9273.d0/5.d0*cn**(-4)
     &     - 8996.d0/5.d0*cn**(-3) + 3879.d0/5.d0*cn**(-2)
     &     - 12081.d0/25.d0*cn**(-1) - 14157.d0/25.d0*(cn+1.d0)**(-6)
     &     + 7725.d0*(cn+1.d0)**(-1) + 57256.d0*(cn+2.d0)**(-1)
     &     + 12898.d0*(cn+3.d0)**(-1)
     &     + 56000.d0*s1(1,cn+1.d0)**2*(cn+1.d0)**(-1)
      c2pls3n = c2pls3n
     &     - 111069.d0/25.d0*s1(1,cn) + 3672.d0*s1(1,cn)*cn**(-3)
     &     - 6158.d0*s1(1,cn)*cn**(-2)
     &     - 230071.d0/20.d0*s1(1,cn)*cn**(-1)
     &     - 146569.d0/40.d0*s1(1,cn)**2
     &     - 2890323.d0/50.d0*s1(1,cn)**2*cn**(-1)
     &     - 58577.d0/150.d0*s1(1,cn)**3
     &     + 2866901.d0/1000.d0*s1(1,cn)**3*cn**(-1)
     &     + 501099.d0/4000.d0*s1(1,cn)**4
     &     + 896.d0/27.d0*s1(1,cn)**4*cn**(-1)
     &     + 1088.d0/27.d0*s1(1,cn)**5
     &     + 256.d0/81.d0*s1(1,cn)**6
     &     + 1280.d0/27.d0*s1(1,cn)**4*s1(2,cn)
     &     + 10880.d0/27.d0*s1(1,cn)**3*s1(2,cn)
     &     + 10240.d0/81.d0*s1(1,cn)**3*s1(3,cn)
      c2pls3n = c2pls3n
     &	   + 1503297.d0/2000.d0*s1(1,cn)**2*s1(2,cn)
     &     + 1792.d0/9.d0*s1(1,cn)**2*s1(2,cn)*cn**(-1)
     &     + 1280.d0/9.d0*s1(1,cn)**2*s1(2,cn)**2
     &     + 21760.d0/27.d0*s1(1,cn)**2*s1(3,cn)
     &     + 2560.d0/9.d0*s1(1,cn)**2*s1(4,cn)
     &     - 58577.d0/50.d0*s1(1,cn)*s1(2,cn)
     &     + 8600703.d0/1000.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 5440.d0/9.d0*s1(1,cn)*s1(2,cn)**2
     &     + 10240.d0/27.d0*s1(1,cn)*s1(2,cn)*s1(3,cn)
     &     + 501099.d0/500.d0*s1(1,cn)*s1(3,cn)
     &     + 7168.d0/27.d0*s1(1,cn)*s1(3,cn)*cn**(-1)
     &     + 10880.d0/9.d0*s1(1,cn)*s1(4,cn)
     &     + 4096.d0/9.d0*s1(1,cn)*s1(5,cn)
     &     + 56000.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 146569.d0/40.d0*s1(2,cn) + 3672.d0*s1(2,cn)*cn**(-2)
      c2pls3n = c2pls3n
     &     - 3198223.d0/50.d0*s1(2,cn)*cn**(-1)
     &     + 1503297.d0/4000.d0*s1(2,cn)**2
     &     + 896.d0/9.d0*s1(2,cn)**2*cn**(-1)
     &     + 1280.d0/27.d0*s1(2,cn)**3
     &     + 21760.d0/27.d0*s1(2,cn)*s1(3,cn)
     &     + 2560.d0/9.d0*s1(2,cn)*s1(4,cn) - 58577.d0/75.d0*s1(3,cn)
     &     + 4702901.d0/500.d0*s1(3,cn)*cn**(-1)
     &     + 10240.d0/81.d0*s1(3,cn)**2 + 1503297.d0/2000.d0*s1(4,cn)
     &     + 1792.d0/9.d0*s1(4,cn)*cn**(-1)
     &     + 8704.d0/9.d0*s1(5,cn) + 10240.d0/27.d0*s1(6,cn)

      c2ps3n=
     &     + fl11ps(nf)*nf * ( - 1486.d0/125.d0
     &     + 640.d0/27.d0*(cn+1.d0)**(-4)
     &     - 3200.d0/81.d0*(cn+1.d0)**(-3)
     &     - 5959.d0/100.d0*(cn+1.d0)**(-2)
     &     + 99703.d0/1000.d0*(cn+1.d0)**(-1)
     &     - 1842.d0/25.d0*(cn+2.d0)**(-5)
     &     + 10437.d0/50.d0*(cn+2.d0)**(-4)
     &     - 1018.d0/5.d0*(cn+2.d0)**(-3)
     &     - 17671.d0/100.d0*(cn+2.d0)**(-1)
     &     + 7.d0/50.d0*(cn+3.d0)**(-1)
     &     + 1003.d0/20.d0*(cn+4.d0)**(-1)
     &     + 363.d0/40.d0*s1(1,cn+2.d0)*(cn+2.d0)**(-1)
     &     - 363.d0/40.d0*s1(1,cn+3.d0)*(cn+3.d0)**(-1) )
      c2ps3n = c2ps3n
     &     + nf * ( + 17662.d0*z2*cn**(-2) - 4162.d0*z2*cn**(-1)
     &     + 4162.d0*z2*(cn+1.d0)**(-1) + 17662.d0*z3*cn**(-1)
     &     + 1600.d0/3.d0*cn**(-6) + 33952.d0/27.d0*cn**(-5)
     &     + 100119.d0/50.d0*cn**(-4) + 2416.d0*cn**(-3)
     &     - 1615.d0*cn**(-2) + 561641.d0/50.d0*cn**(-1)
     &     + 207131.d0/500.d0*(cn-1.d0)**(-2)
     &     - 136591.d0/50.d0*(cn-1.d0)**(-1)
     &     + 9264.d0/5.d0*(cn+1.d0)**(-6) + 6666.d0*(cn+1.d0)**(-3) )
      c2ps3n = c2ps3n
     &     + nf * ( - 13215.d0*(cn+1.d0)**(-1)
     &     + 9551.d0/2.d0*(cn+2.d0)**(-1) - 123.d0/2.d0*(cn+3.d0)**(-1)
     &     - 4162.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-2)
     &     - 542.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 13057.d0/100.d0*s1(1,cn+1.d0)**2*(cn+1.d0)**(-1)
     &     - 6032.d0/81.d0*s1(1,cn+1.d0)**3*(cn+1.d0)**(-1)
     &     - 856.d0/81.d0*s1(1,cn+1.d0)**4*(cn+1.d0)**(-1)
     &     - 1712.d0/27.d0*s1(1,cn+1.d0)**2*s1(2,cn+1.d0)*
     &    (cn+1.d0)**(-1)
     &     - 6032.d0/27.d0*s1(1,cn+1.d0)*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 6848.d0/81.d0*s1(1,cn+1.d0)*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 17662.d0*s1(1,cn)*cn**(-3) + 4162.d0*s1(1,cn)*cn**(-2) )
      c2ps3n = c2ps3n
     &     + nf * ( + 542.d0*s1(1,cn)*cn**(-1)
     &     + 13057.d0/100.d0*s1(1,cn)**2*cn**(-1)
     &     + 6032.d0/81.d0*s1(1,cn)**3*cn**(-1)
     &     + 856.d0/81.d0*s1(1,cn)**4*cn**(-1)
     &     + 1712.d0/27.d0*s1(1,cn)**2*s1(2,cn)*cn**(-1)
     &     + 6032.d0/27.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 6848.d0/81.d0*s1(1,cn)*s1(3,cn)*cn**(-1)
     &     - 429257.d0/100.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 856.d0/27.d0*s1(2,cn+1.d0)**2*(cn+1.d0)**(-1)
     &     - 17662.d0*s1(2,cn)*cn**(-2)
     &     + 429257.d0/100.d0*s1(2,cn)*cn**(-1)
     &     + 856.d0/27.d0*s1(2,cn)**2*cn**(-1)
     &     - 12064.d0/81.d0*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 1418558.d0/81.d0*s1(3,cn)*cn**(-1)
     &     - 1712.d0/27.d0*s1(4,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 1712.d0/27.d0*s1(4,cn)*cn**(-1) )
      c2ps3n = c2ps3n
     &     + nf**2 * ( - 1603.d0/5.d0*z2*cn**(-2)
     &     + 677.d0/5.d0*z2*cn**(-1) - 677.d0/5.d0*z2*(cn+1.d0)**(-1)
     &     - 1603.d0/5.d0*z3*cn**(-1) - 1472.d0/27.d0*cn**(-5)
     &     + 7136.d0/81.d0*cn**(-4) - 4519.d0/25.d0*cn**(-3)
     &     + 1049.d0/5.d0*cn**(-2) - 1302563.d0/5000.d0*cn**(-1)
     &     + 201213.d0/5000.d0*(cn-1.d0)**(-1)
     &     + 3621.d0/25.d0*(cn+1.d0)**(-4)
     &     - 2154.d0/5.d0*(cn+1.d0)**(-3)
     &     + 28007.d0/100.d0*(cn+1.d0)**(-1)
     &     - 1187.d0/5.d0*(cn+2.d0)**(-1)
     &     + 888.d0/5.d0*(cn+3.d0)**(-1)
     &     + 677.d0/5.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-2)
     &     + 2309.d0/100.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-1) )
      c2ps3n = c2ps3n
     &     + nf**2 * ( - 208.d0/81.d0*s1(1,cn+1.d0)**2*(cn+1.d0)**(-1)
     &     - 64.d0/81.d0*s1(1,cn+1.d0)**3*(cn+1.d0)**(-1)
     &     - 64.d0/27.d0*s1(1,cn+1.d0)*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 1603.d0/5.d0*s1(1,cn)*cn**(-3)
     &     - 677.d0/5.d0*s1(1,cn)*cn**(-2)
     &     - 2309.d0/100.d0*s1(1,cn)*cn**(-1)
     &     + 208.d0/81.d0*s1(1,cn)**2*cn**(-1)
     &     + 64.d0/81.d0*s1(1,cn)**3*cn**(-1)
     &     + 64.d0/27.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 53797.d0/405.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 1603.d0/5.d0*s1(2,cn)*cn**(-2) 
     &     - 53797.d0/405.d0*s1(2,cn)*cn**(-1)
     &     - 128.d0/81.d0*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 130483.d0/405.d0*s1(3,cn)*cn**(-1) )

      c2g3n = 
     &     + fl11g(nf)*nf * ( - 1862.d0/25.d0*z2*(cn+1.d0)**(-2)
     &     + 792.d0/25.d0*z2*(cn+1.d0)**(-1)
     &     - 1862.d0/25.d0*z3*(cn+1.d0)**(-1)
     &     + 160.d0/9.d0*(cn+1.d0)**(-5) - 1040.d0/27.d0*(cn+1.d0)**(-4)
     &     - 804.d0/25.d0*(cn+1.d0)**(-3) + 412.d0/5.d0*(cn+1.d0)**(-2)
     &     - 6447.d0/100.d0*(cn+1.d0)**(-1)
     &     - 3381.d0/50.d0*(cn+2.d0)**(-4) + 608.d0/5.d0*(cn+2.d0)**(-1)
     &     - 2291.d0/50.d0*(cn+3.d0)**(-1)
     &     + 1862.d0/25.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-3)
     &     - 792.d0/25.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-2)
     &     - 476.d0/25.d0*s1(1,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 623.d0/1000.d0*s1(1,cn+1.d0)**3*(cn+1.d0)**(-1)
     &     + 1869.d0/1000.d0*s1(1,cn+1.d0)*s1(2,cn+1.d0)*
     &    (cn+1.d0)**(-1) + 3211.d0/1000.d0*s1(1,cn)**2*cn**(-1) )
      c2g3n = c2g3n
     &     + fl11g(nf)*nf * ( - 623.d0/1000.d0*s1(1,cn)**3*cn**(-1)
     &     - 1869.d0/1000.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 1862.d0/25.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-2)
     &     - 792.d0/25.d0*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 3211.d0/1000.d0*s1(2,cn)*cn**(-1)
     &     + 37863.d0/500.d0*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 623.d0/500.d0*s1(3,cn)*cn**(-1) )
      c2g3n = c2g3n
     &     + nf * ( + 5.d0/8.d0 - 474020.d0*z2*cn**(-2)
     &     - 138230.d0*z2*cn**(-1) - 474020.d0*z3*cn**(-1)
     &     + 5360.d0/3.d0*cn**(-6) + 19844.d0/27.d0*cn**(-5)
     &     + 8640.d0*cn**(-4) - 7008.d0/5.d0*cn**(-3)
     &     + 11860.d0*cn**(-2) - 15729.d0*cn**(-1)
     &     + 932089.d0/1000.d0*(cn-1.d0)**(-2)
     &     - 318127.d0/50.d0*(cn-1.d0)**(-1) + 30396.d0*(cn+1.d0)**(-6)
     &     - 488300.d0*(cn+1.d0)**(-3) - 310510.d0*(cn+1.d0)**(-1)
     &     + 331570.d0*(cn+2.d0)**(-1) )
      c2g3n = c2g3n
     &     + nf * ( - 1372.d0*s1(1,cn+1.d0)**4*(cn+1.d0)**(-1)
     &     - 8232.d0*s1(1,cn+1.d0)**2*s1(2,cn+1.d0)*(cn+1.d0)**(-1)
     &     - 10976.d0*s1(1,cn+1.d0)*s1(3,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 474020.d0*s1(1,cn)*cn**(-3) + 138230.d0*s1(1,cn)*cn**(-2)
     &     + 2405.d0*s1(1,cn)*cn**(-1)
     &     + 4896.d0/5.d0*s1(1,cn)**2*cn**(-1)
     &     - 8931.d0/100.d0*s1(1,cn)**3*cn**(-1)
     &     + 22825.d0/18.d0*s1(1,cn)**4*cn**(-1)
     &     - 322.d0/27.d0*s1(1,cn)**5*cn**(-1)
     &     - 3220.d0/27.d0*s1(1,cn)**3*s1(2,cn)*cn**(-1)
     &     + 22825.d0/3.d0*s1(1,cn)**2*s1(2,cn)*cn**(-1) )
      c2g3n = c2g3n
     &     + nf * ( - 6440.d0/27.d0*s1(1,cn)**2*s1(3,cn)*cn**(-1)
     &     - 26793.d0/100.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     - 1610.d0/9.d0*s1(1,cn)*s1(2,cn)**2*cn**(-1)
     &     + 91300.d0/9.d0*s1(1,cn)*s1(3,cn)*cn**(-1)
     &     - 3220.d0/9.d0*s1(1,cn)*s1(4,cn)*cn**(-1)
     &     - 4116.d0*s1(2,cn+1.d0)**2*(cn+1.d0)**(-1)
     &     + 474020.d0*s1(2,cn)*cn**(-2)
     &     + 696046.d0/5.d0*s1(2,cn)*cn**(-1)
     &     + 22825.d0/6.d0*s1(2,cn)**2*cn**(-1)
     &     - 6440.d0/27.d0*s1(2,cn)*s1(3,cn)*cn**(-1)
     &     + 23692069.d0/50.d0*s1(3,cn)*cn**(-1)
     &     - 8232.d0*s1(4,cn+1.d0)*(cn+1.d0)**(-1)
     &     + 22825.d0/3.d0*s1(4,cn)*cn**(-1)
     &     - 2576.d0/9.d0*s1(5,cn)*cn**(-1) )
      c2g3n = c2g3n
     &     + nf**2 * ( - 9292.d0*z2*cn**(-2) + 1583.d0/2.d0*z2*cn**(-1)
     &     - 9292.d0*z3*cn**(-1) - 160.d0*cn**(-6)
     &     + 6472.d0/27.d0*cn**(-5) - 3129.d0/5.d0*cn**(-4)
     &     + 836.d0*cn**(-3) - 739.d0*cn**(-2)
     &     + 2381.d0/500.d0*cn**(-1) + 84423.d0/1000.d0*(cn-1.d0)**(-1)
     &     - 8038.d0*(cn+1.d0)**(-3) - 190.d0*(cn+1.d0)**(-1)
     &     - 4092.d0/5.d0*(cn+2.d0)**(-1) + 9292.d0*s1(1,cn)*cn**(-3) )
      c2g3n = c2g3n
     &     + nf**2 * ( - 1583.d0/2.d0*s1(1,cn)*cn**(-2)
     &     + 2261.d0/10.d0*s1(1,cn)*cn**(-1)
     &     + 3607.d0/1000.d0*s1(1,cn)**2*cn**(-1)
     &     + 368.d0/25.d0*s1(1,cn)**3*cn**(-1)
     &     + 131.d0/81.d0*s1(1,cn)**4*cn**(-1)
     &     + 262.d0/27.d0*s1(1,cn)**2*s1(2,cn)*cn**(-1)
     &     + 1104.d0/25.d0*s1(1,cn)*s1(2,cn)*cn**(-1)
     &     + 1048.d0/81.d0*s1(1,cn)*s1(3,cn)*cn**(-1)
     &     + 9292.d0*s1(2,cn)*cn**(-2)
     &     - 787893.d0/1000.d0*s1(2,cn)*cn**(-1)
     &     + 131.d0/27.d0*s1(2,cn)**2*cn**(-1)
     &     + 233036.d0/25.d0*s1(3,cn)*cn**(-1)
     &     + 262.d0/27.d0*s1(4,cn)*cn**(-1) )

      return
      end
* =====================================================================
c--------------------
      subroutine cl1p(cn,nf,clpls1n,clg1n)
      IMPLICIT complex*16 (A-H,O-Z)
      integer nf

      clpls1n=
     &  + 16.d0/3.d0*(cn+1)**(-1)

      clg1n=
     &  + nf * ( 8.d0*(cn+1)**(-1) - 8.d0*(cn+2)**(-1) )

      return
      end
* =====================================================================
c--------------------
      subroutine cl2p(cn,nf,clpls2n,clps2n,clg2n)
      IMPLICIT complex*16 (A-H,O-Z)

      integer nf
      real*8 dl,pi,ge,z2,z3,z4,z5,z6

      pi = 3.1415926535897932385d0
      dl = dlog(2.0d0)
      ge = 5.7721566490153286061d-1
      z2 = 1.6449340668482264365d0
      z3 = 1.2020569031595942854d0
      z4 = 1.0823232337111381916d0
      z5 = 1.0369277551433699263d0
      z6 = 1.0173430619844491398d0

      clpls2n=
     &  + nf * ( 32.d0/9.d0*cn**(-1) + 64.d0/9.d0*(cn+1)**(-2) - 400.d0/
     &    27.d0*(cn+1)**(-1) - 32.d0/9.d0*s1(1,1.d0 + cn)*(cn+1)**(-1)
     &     )
      clpls2n=clpls2n
     &     - 3.d0/250.d0 + 42047.d0/500.d0*z2*cn**(-1) + 128.d0
     &    /9.d0*cn**(-2) - 18669.d0/500.d0*cn**(-1) + 1841.d0/50.d0*
     &    (cn+1)**(-3) - 329.d0/10.d0*(cn+1)**(-2) + 8953.d0/100.d0*
     &    (cn+1)**(-1) + 1691.d0/50.d0*(cn+2)**(-1) + 93.d0/2.d0*s1(1,1.
     &    d0 + cn)*(cn+1)**(-1) + 128.d0/9.d0*s1(1,1.d0 + cn)**2*
     &    (cn+1)**(-1) - 42047.d0/500.d0*s1(1,cn)*cn**(-2) + 128.d0/9.d0
     &    *s1(2,1.d0 + cn)*(cn+1)**(-1) - 42047.d0/500.d0*s1(2,cn)*
     &    cn**(-1)

      clps2n=
     &  + nf * ( 421.d0/500.d0*cn**(-3) - 2809.d0/100.d0*cn**(-2) + 
     &    1319.d0/50.d0*cn**(-1) - 237.d0/100.d0*(cn-1)**(-1) + 76.d0/
     &    25.d0*(cn+1)**(-3) + 2809.d0/100.d0*(cn+1)**(-2) - 1623.d0/25.
     &    d0*(cn+1)**(-1) + 3009.d0/50.d0*(cn+2)**(-1) - 1927.d0/100.d0
     &    *(cn+3)**(-1) + 9273.d0/250.d0*s1(1,1.d0 + cn)*(cn+1)**(-1)
     &     - 6591.d0/250.d0*s1(1,2.d0 + cn)*(cn+2)**(-1) + 1303.d0/250.d
     &    0*s1(1,3.d0 + cn)*(cn+3)**(-1) - 797.d0/50.d0*s1(1,cn)*
     &    cn**(-1) )

      clg2n=
     &  + nf * (  - 1161.d0*z2*(cn+1)**(-1) - 1983.d0/50.d0*cn**(-2) + 
     &    5333.d0/1000.d0*cn**(-1) - 5333.d0/1000.d0*(cn-1)**(-1) + 
     &    3003.d0/25.d0*(cn+1)**(-3) + 1983.d0/50.d0*(cn+1)**(-2) + 
     &    1161.d0*s1(1,1.d0 + cn)*(cn+1)**(-2) + 4324.d0/5.d0*s1(1,1.d0
     &     + cn)*(cn+1)**(-1) - 7197.d0/50.d0*s1(1,1.d0 + cn)**2*
     &    (cn+1)**(-1) + 246.d0/5.d0*s1(1,2.d0 + cn)**2*(cn+2)**(-1) - 
     &    4324.d0/5.d0*s1(1,cn)*cn**(-1) + 4737.d0/50.d0*s1(1,cn)**2*
     &    cn**(-1) + 50853.d0/50.d0*s1(2,1.d0 + cn)*(cn+1)**(-1) + 246.d
     &    0/5.d0*s1(2,2.d0 + cn)*(cn+2)**(-1) + 4737.d0/50.d0*s1(2,cn)*
     &    cn**(-1) )

      return
      end
* =====================================================================
c--------------------
      subroutine cl3p(cn,nf,clpls3n,clps3n,clg3n)
      IMPLICIT complex*16 (A-H,O-Z)

      dimension fl11ns(7),fl11ps(7),fl11g(7)
      data fl11ns 
     &     /0.d0, 2.d0, 0.5d0, 0.d0, 0.5d0, 0.2d0, 0.5d0 /
      data fl11g  
     &     /0.d0, 1.d0, 0.1d0, 0.d0, 0.1d0, 0.018181818d0, 0.1d0 /
      data fl11ps 
     &     /0.d0,-1.d0,-0.4d0, 0.d0,-0.4d0,-0.181818182d0,-0.4d0 /

      integer nf
      real*8 dl,pi,ge,z2,z3,z4,z5,z6

      pi = 3.1415926535897932385d0
      dl = dlog(2.0d0)
      ge = 5.7721566490153286061d-1
      z2 = 1.6449340668482264365d0
      z3 = 1.2020569031595942854d0
      z4 = 1.0823232337111381916d0
      z5 = 1.0369277551433699263d0
      z6 = 1.0173430619844491398d0

      clpls3n=
     &  + fl11ns(nf)*nf * ( 640.d0/27.d0*(cn+1)**(-4) - 1280.d0/81.d0*
     &    (cn+1)**(-3) - 9773.d0/1000.d0*(cn+1)**(-2) + 80283.d0/1000.d0
     &    *(cn+1)**(-1) + 3416.d0/25.d0*(cn+2)**(-3) - 1819.d0/5.d0*
     &    (cn+2)**(-2) + 4281.d0/20.d0*(cn+2)**(-1) - 37567.d0/100.d0*
     &    (cn+3)**(-1) + 2731.d0/50.d0*(cn+4)**(-1) )
      clpls3n=clpls3n
     &     + nf * ( 3.d0/500.d0 + 589.d0*z2*cn**(-2) - 8873.d0/
     &    10.d0*z2*cn**(-1) + 11768.d0/25.d0*z3*cn**(-1) + 128.d0/3.d0*
     &    cn**(-3) - 20073.d0/100.d0*cn**(-2) + 2042.d0/5.d0*cn**(-1)
     &     + 3584.d0/27.d0*(cn+1)**(-4) - 1869.d0/200.d0*(cn+1)**(-1)
     &     - 9193.d0/10.d0*(cn+2)**(-1) - 2397.d0/10.d0*s1(1,1.d0 + cn)
     &    **2*(cn+1)**(-1) + 2063.d0/100.d0*s1(1,1.d0 + cn)**3*
     &    (cn+1)**(-1) + 6189.d0/100.d0*s1(1,1.d0 + cn)*s1(2,1.d0 + cn)
     &    *(cn+1)**(-1) - 2957.d0/25.d0*s1(1,cn)*z2*cn**(-1) - 589.d0*
     &    s1(1,cn)*cn**(-3) + 8873.d0/10.d0*s1(1,cn)*cn**(-2) - 3441.d0/
     &    10.d0*s1(1,cn)*cn**(-1) + 2957.d0/50.d0*s1(1,cn)**2*cn**(-2)
     &     + 2547.d0/20.d0*s1(1,cn)**2*cn**(-1) - 269503.d0/8100.d0*s1(
     &    1,cn)**3*cn**(-1) + 49853.d0/2700.d0*s1(1,cn)*s1(2,cn)*
     &    cn**(-1) - 2397.d0/10.d0*s1(2,1.d0 + cn)*(cn+1)**(-1) - 26493.
     &    d0/50.d0*s1(2,cn)*cn**(-2) + 20293.d0/20.d0*s1(2,cn)*cn**(-1)
     &     + 2063.d0/50.d0*s1(3,1.d0 + cn)*(cn+1)**(-1) - 2175919.d0/
     &    4050.d0*s1(3,cn)*cn**(-1) )
      clpls3n=clpls3n
     &     + nf**2 * ( 128.d0/27.d0*cn**(-2) - 1216.d0/81.d0*
     &    cn**(-1) + 128.d0/9.d0*(cn+1)**(-3) - 3200.d0/81.d0*
     &    (cn+1)**(-2) + 10144.d0/243.d0*(cn+1)**(-1) - 256.d0/27.d0*
     &    s1(1,1.d0 + cn)*(cn+1)**(-2) + 1600.d0/81.d0*s1(1,1.d0 + cn)*
     &    (cn+1)**(-1) + 64.d0/27.d0*s1(1,1.d0 + cn)**2*(cn+1)**(-1) - 
     &    128.d0/27.d0*s1(1,cn)*cn**(-1) - 64.d0/27.d0*s1(2,1.d0 + cn)*
     &    (cn+1)**(-1) )
      clpls3n=clpls3n
     &     + 113.d0/1000.d0 - 8447.d0/5.d0*z2*cn**(-2) - 2724.d
     &    0*z3*cn**(-1) + 2560.d0/27.d0*cn**(-4) - 14912.d0/27.d0*
     &    cn**(-3) + 13557.d0/10.d0*cn**(-2) - 4441.d0/2.d0*cn**(-1) - 
     &    6249.d0/5.d0*(cn+1)**(-4) - 7884.d0*(cn+1)**(-1) + 4168.d0*
     &    (cn+2)**(-1) - 1253.d0/10.d0*s1(1,1.d0 + cn)**3*(cn+1)**(-1)
     &     - 978.d0/5.d0*s1(1,1.d0 + cn)**4*(cn+1)**(-1) - 5868.d0/5.d0
     &    *s1(1,1.d0 + cn)**2*s1(2,1.d0 + cn)*(cn+1)**(-1) - 3759.d0/10.
     &    d0*s1(1,1.d0 + cn)*s1(2,1.d0 + cn)*(cn+1)**(-1) - 7824.d0/5.d0
     &    *s1(1,1.d0 + cn)*s1(3,1.d0 + cn)*(cn+1)**(-1) - 5173.d0/5.d0*
     &    s1(1,cn)*z2*cn**(-1) + 8447.d0/5.d0*s1(1,cn)*cn**(-3) + 2729.d
     &    0*s1(1,cn)*cn**(-1) + 5173.d0/10.d0*s1(1,cn)**2*cn**(-2) + 
     &    3253.d0/5.d0*s1(1,cn)**2*cn**(-1) + 3027.d0/10.d0*s1(1,cn)**3
     &    *cn**(-1) + 28966.d0/135.d0*s1(1,cn)**4*cn**(-1) + 57932.d0/
     &    45.d0*s1(1,cn)**2*s1(2,cn)*cn**(-1) + 19427.d0/10.d0*s1(1,cn)
     &    *s1(2,cn)*cn**(-1) + 231728.d0/135.d0*s1(1,cn)*s1(3,cn)*
     &    cn**(-1)
      clpls3n=clpls3n
     &     - 2934.d0/5.d0*s1(2,1.d0 + cn)**2*(cn+1)**(-1) + 
     &    22067.d0/10.d0*s1(2,cn)*cn**(-2) + 3253.d0/5.d0*s1(2,cn)*
     &    cn**(-1) + 28966.d0/45.d0*s1(2,cn)**2*cn**(-1) - 1253.d0/5.d0
     &    *s1(3,1.d0 + cn)*(cn+1)**(-1) + 16647.d0/5.d0*s1(3,cn)*
     &    cn**(-1) - 5868.d0/5.d0*s1(4,1.d0 + cn)*(cn+1)**(-1) + 57932.d
     &    0/45.d0*s1(4,cn)*cn**(-1)

      clps3n=
     &  + fl11ps(nf)*nf * ( 640.d0/27.d0*(cn+1)**(-4) - 1280.d0/81.d0*
     &    (cn+1)**(-3) - 9773.d0/1000.d0*(cn+1)**(-2) + 80283.d0/1000.d0
     &    *(cn+1)**(-1) + 3416.d0/25.d0*(cn+2)**(-3) - 1819.d0/5.d0*
     &    (cn+2)**(-2) + 4281.d0/20.d0*(cn+2)**(-1) - 37567.d0/100.d0*
     &    (cn+3)**(-1) + 2731.d0/50.d0*(cn+4)**(-1) )
      clps3n=clps3n
     &     + nf * ( 4368.d0*z2*cn**(-2) - 6059.d0*z2*cn**(-1)
     &     + 6059.d0*z2*(cn+1)**(-1) + 4368.d0*z3*cn**(-1) + 3200.d0/9.d
     &    0*cn**(-4) + 5696.d0/9.d0*cn**(-3) - 192.d0/5.d0*cn**(-2) + 
     &    48773.d0/50.d0*cn**(-1) + 182.d0*(cn-1)**(-2) - 88553.d0/100.d
     &    0*(cn-1)**(-1) - 718.d0/5.d0*(cn+1)**(-2) - 33033.d0/100.d0*
     &    (cn+1)**(-1) + 6382.d0/5.d0*(cn+2)**(-1) - 1036.d0*
     &    (cn+3)**(-1) - 6059.d0*s1(1,1.d0 + cn)*(cn+1)**(-2) + 10248.d0
     &    *s1(1,1.d0 + cn)*(cn+1)**(-1) + 7936.d0/9.d0*s1(1,1.d0 + cn)
     &    **2*(cn+1)**(-1) + 3136.d0/27.d0*s1(1,1.d0 + cn)**3*
     &    (cn+1)**(-1) + 3136.d0/9.d0*s1(1,1.d0 + cn)*s1(2,1.d0 + cn)*
     &    (cn+1)**(-1) - 5124.d0*s1(1,2.d0 + cn)*(cn+2)**(-1) - 3968.d0/
     &    9.d0*s1(1,2.d0 + cn)**2*(cn+2)**(-1) - 1568.d0/27.d0*s1(1,2.d0
     &     + cn)**3*(cn+2)**(-1) - 1568.d0/9.d0*s1(1,2.d0 + cn)*s1(2,2.d
     &    0 + cn)*(cn+2)**(-1) - 4368.d0*s1(1,cn)*cn**(-3) + 6059.d0*
     &    s1(1,cn)*cn**(-2) - 5124.d0*s1(1,cn)*cn**(-1) - 3968.d0/9.d0*
     &    s1(1,cn)**2*cn**(-1) )
      clps3n=clps3n
     &     + nf * (  - 1568.d0/27.d0*s1(1,cn)**3*cn**(-1) - 
     &    1568.d0/9.d0*s1(1,cn)*s1(2,cn)*cn**(-1) - 46595.d0/9.d0*s1(2,
     &    1.d0 + cn)*(cn+1)**(-1) - 3968.d0/9.d0*s1(2,2.d0 + cn)*
     &    (cn+2)**(-1) - 4368.d0*s1(2,cn)*cn**(-2) + 50563.d0/9.d0*s1(2
     &    ,cn)*cn**(-1) + 6272.d0/27.d0*s1(3,1.d0 + cn)*(cn+1)**(-1) - 
     &    3136.d0/27.d0*s1(3,2.d0 + cn)*(cn+2)**(-1) - 121072.d0/27.d0*
     &    s1(3,cn)*cn**(-1) )
      clps3n=clps3n
     &     + nf**2 * ( 1759.d0/25.d0*z2*cn**(-2) - 3653.d0/50.d
     &    0*z2*cn**(-1) + 3653.d0/50.d0*z2*(cn+1)**(-1) + 1759.d0/25.d0
     &    *z3*cn**(-1) - 256.d0/9.d0*cn**(-3) + 6941.d0/100.d0*cn**(-2)
     &     - 47319.d0/500.d0*cn**(-1) + 40239.d0/1000.d0*(cn-1)**(-1)
     &     - 1762.d0/25.d0*(cn+1)**(-3) - 6941.d0/100.d0*(cn+1)**(-2)
     &     + 138399.d0/1000.d0*(cn+1)**(-1) - 3846.d0/25.d0*
     &    (cn+2)**(-1) + 1746.d0/25.d0*(cn+3)**(-1) - 3653.d0/50.d0*s1(
     &    1,1.d0 + cn)*(cn+1)**(-2) + 1476.d0/25.d0*s1(1,1.d0 + cn)*
     &    (cn+1)**(-1) + 64.d0/9.d0*s1(1,1.d0 + cn)**2*(cn+1)**(-1) - 
     &    738.d0/25.d0*s1(1,2.d0 + cn)*(cn+2)**(-1) - 32.d0/9.d0*s1(1,2.
     &    d0 + cn)**2*(cn+2)**(-1) - 1759.d0/25.d0*s1(1,cn)*cn**(-3) + 
     &    3653.d0/50.d0*s1(1,cn)*cn**(-2) - 738.d0/25.d0*s1(1,cn)*
     &    cn**(-1) - 32.d0/9.d0*s1(1,cn)**2*cn**(-1) - 29677.d0/450.d0*
     &    s1(2,1.d0 + cn)*(cn+1)**(-1) - 32.d0/9.d0*s1(2,2.d0 + cn)*
     &    (cn+2)**(-1) - 1759.d0/25.d0*s1(2,cn)*cn**(-2) + 31277.d0/450.
     &    d0*s1(2,cn)*cn**(-1) )
      clps3n=clps3n
     &     + nf**2 * (  - 1759.d0/25.d0*s1(3,cn)*cn**(-1) )

      clg3n=
     &  + fl11g(nf)*nf**2 * ( 160.d0/9.d0*(cn+1)**(-5) - 560.d0/27.d0*
     &    (cn+1)**(-4) - 154.d0/5.d0*(cn+1)**(-3) + 3583.d0/50.d0*
     &    (cn+1)**(-2) - 13349.d0/200.d0*(cn+1)**(-1) - 339.d0/25.d0*
     &    (cn+2)**(-4) + 2201.d0/500.d0*(cn+2)**(-3) - 121.d0/1000.d0*
     &    (cn+2)**(-2) + 1684.d0/25.d0*(cn+2)**(-1) - 123.d0/200.d0*
     &    (cn+3)**(-1) - 493.d0/25.d0*s1(1,1.d0 + cn)*(cn+1)**(-1) - 31.
     &    d0/20.d0*s1(1,1.d0 + cn)**2*(cn+1)**(-1) - 21.d0/2000.d0*s1(1
     &    ,1.d0 + cn)**3*(cn+1)**(-1) - 63.d0/2000.d0*s1(1,1.d0 + cn)*
     &    s1(2,1.d0 + cn)*(cn+1)**(-1) + 493.d0/25.d0*s1(1,2.d0 + cn)*
     &    (cn+2)**(-1) + 31.d0/20.d0*s1(1,cn)**2*cn**(-1) + 21.d0/2000.d
     &    0*s1(1,cn)**3*cn**(-1) + 63.d0/2000.d0*s1(1,cn)*s1(2,cn)*
     &    cn**(-1) - 31.d0/20.d0*s1(2,1.d0 + cn)*(cn+1)**(-1) + 31.d0/
     &    20.d0*s1(2,cn)*cn**(-1) - 21.d0/1000.d0*s1(3,1.d0 + cn)*
     &    (cn+1)**(-1) + 21.d0/1000.d0*s1(3,cn)*cn**(-1) )
      clg3n=clg3n
     &     + nf * ( 176476.d0*z2*cn**(-2) - 72549.d0*z2*cn**(-1)
     &     + 176476.d0*z3*cn**(-1) + 10720.d0/9.d0*cn**(-4) + 4537.d0/
     &    50.d0*cn**(-3) + 1117.d0*cn**(-2) + 57537.d0/10.d0*cn**(-1)
     &     + 204753.d0/500.d0*(cn-1)**(-2) - 20447.d0/10.d0*
     &    (cn-1)**(-1) + 133546.d0*(cn+1)**(-3) - 37223.d0*(cn+1)**(-1)
     &     + 23981.d0*(cn+2)**(-1) + 9533.d0*(cn+3)**(-1) + 53160.d0*
     &    s1(1,1.d0 + cn)*(cn+1)**(-1) - 6319.d0*s1(1,1.d0 + cn)**2*
     &    (cn+1)**(-1) - 47024.d0/27.d0*s1(1,1.d0 + cn)**3*(cn+1)**(-1)
     &     - 144.d0*s1(1,1.d0 + cn)**4*(cn+1)**(-1) - 864.d0*s1(1,1.d0
     &     + cn)**2*s1(2,1.d0 + cn)*(cn+1)**(-1) - 47024.d0/9.d0*s1(1,1.
     &    d0 + cn)*s1(2,1.d0 + cn)*(cn+1)**(-1) - 1152.d0*s1(1,1.d0 + 
     &    cn)*s1(3,1.d0 + cn)*(cn+1)**(-1) - 176476.d0*s1(1,cn)*
     &    cn**(-3) + 72549.d0*s1(1,cn)*cn**(-2) - 53160.d0*s1(1,cn)*
     &    cn**(-1) + 6319.d0*s1(1,cn)**2*cn**(-1) + 47024.d0/27.d0*s1(1
     &    ,cn)**3*cn**(-1) + 144.d0*s1(1,cn)**4*cn**(-1) + 864.d0*s1(1,
     &    cn)**2*s1(2,cn)*cn**(-1) )
      clg3n=clg3n
     &     + nf * ( 47024.d0/9.d0*s1(1,cn)*s1(2,cn)*cn**(-1) + 
     &    1152.d0*s1(1,cn)*s1(3,cn)*cn**(-1) - 6319.d0*s1(2,1.d0 + cn)*
     &    (cn+1)**(-1) - 432.d0*s1(2,1.d0 + cn)**2*(cn+1)**(-1) - 
     &    176476.d0*s1(2,cn)*cn**(-2) + 78868.d0*s1(2,cn)*cn**(-1) + 
     &    432.d0*s1(2,cn)**2*cn**(-1) - 94048.d0/27.d0*s1(3,1.d0 + cn)*
     &    (cn+1)**(-1) - 4670804.d0/27.d0*s1(3,cn)*cn**(-1) - 864.d0*
     &    s1(4,1.d0 + cn)*(cn+1)**(-1) + 864.d0*s1(4,cn)*cn**(-1) )
      clg3n=clg3n
     &     + nf**2 * ( 712.d0/25.d0*z2*cn**(-2) - 3113.d0/10.d0*
     &    z2*cn**(-1) + 712.d0/25.d0*z3*cn**(-1) - 320.d0/3.d0*cn**(-4)
     &     + 2000.d0/9.d0*cn**(-3) - 366.d0*cn**(-2) + 9907901.d0/20270.
     &    d0*cn**(-1) + 179397.d0/2027.d0*(cn-1)**(-1) - 4617.d0/25.d0*
     &    (cn+1)**(-4) - 13063.d0/10.d0*(cn+1)**(-1) + 729.d0*
     &    (cn+2)**(-1) - 21033.d0/10.d0*s1(1,1.d0 + cn)*(cn+1)**(-1) + 
     &    1216.d0/9.d0*s1(1,1.d0 + cn)**2*(cn+1)**(-1) + 32.d0/3.d0*s1(
     &    1,1.d0 + cn)**3*(cn+1)**(-1) + 32.d0*s1(1,1.d0 + cn)*s1(2,1.d0
     &     + cn)*(cn+1)**(-1) + 1511.d0*s1(1,2.d0 + cn)*(cn+2)**(-1) - 
     &    712.d0/25.d0*s1(1,cn)*cn**(-3) + 3113.d0/10.d0*s1(1,cn)*
     &    cn**(-2) + 5923.d0/10.d0*s1(1,cn)*cn**(-1) - 1216.d0/9.d0*s1(
     &    1,cn)**2*cn**(-1) - 32.d0/3.d0*s1(1,cn)**3*cn**(-1) - 32.d0*
     &    s1(1,cn)*s1(2,cn)*cn**(-1) + 1216.d0/9.d0*s1(2,1.d0 + cn)*
     &    (cn+1)**(-1) - 712.d0/25.d0*s1(2,cn)*cn**(-2) + 15857.d0/90.d0
     &    *s1(2,cn)*cn**(-1) + 64.d0/3.d0*s1(3,1.d0 + cn)*(cn+1)**(-1)
     &     - 3736.d0/75.d0*s1(3,cn)*cn**(-1) )

      return
      end
* =====================================================================
*
* ..harmonic sums and the like
*
c--------------
      complex*16 function s1(k,cn)
      IMPLICIT complex*16 (A-H,O-Z)
      integer k

      integer nf
      real*8 dl,pi,ge,z2,z3,z4,z5,z6

      pi = 3.1415926535897932385d0
      dl = dlog(2.0d0)
      ge = 5.7721566490153286061d-1
      z2 = 1.6449340668482264365d0
      z3 = 1.2020569031595942854d0
      z4 = 1.0823232337111381916d0
      z5 = 1.0369277551433699263d0
      z6 = 1.0173430619844491398d0

      if (k.eq.1) then
        cn1=cn+1
        call psi0(cn1,res)
        s1=ge + res
      end if
      if (k.eq.2) then
        cn1=cn+1
        call psi1(cn1,res)
        s1=z2 - res
      end if
      if (k.eq.3) then
        cn1=cn+1
        call psi2(cn1,res)
        s1=z3 + res/2.d0
      end if
      if (k.eq.4) then
        cn1=cn+1
        call psi3(cn1,res)
        s1=z4 - res/6.d0
      end if
      if (k.eq.5) then
        cn1=cn+1
        call psi4(cn1,res)
        s1=z5 + res/24.d0
      end if
      if (k.eq.6) then
        cn1=cn+1
        call psi5(cn1,res)
        s1=z6 - res/120.d0
      end if
      return
      end

*
* ..psi(cn) for complex argument
*      
c--------------
      subroutine psi0(cn,res)
      real*8 r
      complex*16 cn,t,t0,y,y2,res

      t=dcmplx(0.d0,0.d0)
2     r=dsqrt(dreal(cn)**2+dimag(cn)**2)
      if(r.gt.16.d0) goto 1
      t=t-1.d0/cn
      cn=cn+1.d0
      goto 2

1     y=1.d0/cn
      y2=y*y

      t0 = 
     & (-1.d0/2.d0+(-1.d0/12.d0+(1.d0/120.d0+(-1.d0/252.d0+(1.d0/240.d0
     & +(-1.d0/132.d0+(691.d0/32760.d0+(-1.d0/12.d0+(3617.d0/8160.d0
     & +(-43867.d0/14364.d0+(174611.d0/6600.d0+(-77683.d0/276.d0
     & +(236364091.d0/65520.d0+(-657931.d0/12.d0+3392780147.d0/3480.d0
     & *y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y)*y-log(y)

      res=t+t0

      return
      end
*
* ..psi'(cn) (1st derivative) for complex argument
*      
c--------------
      subroutine psi1(cn,res)
      real*8 r
      complex*16 cn,t,t0,y,y2,res

      t=dcmplx(0.d0,0.d0)
2     r=dsqrt(dreal(cn)**2+dimag(cn)**2)
      if(r.gt.16.d0) goto 1
      t=t+1.d0/cn**2
      cn=cn+1.d0
      goto 2
1     y=1.d0/cn
      y2=y*y

      t0 = 
     & (1.d0+(1.d0/2.d0+(1.d0/6.d0+(-1.d0/30.d0+(1.d0/42.d0
     & +(-1.d0/30.d0+(5.d0/66.d0+(-691.d0/2730.d0+(7.d0/6.d0
     & +(-3617.d0/510.d0+(43867.d0/798.d0+(-174611.d0/330.d0
     & +(854513.d0/138.d0+(-236364091.d0/2730.d0+(8553103.d0/6
     & -23749461029.d0/870*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)
     & *y2)*y2)*y)*y)*y

      res=t+t0

      return
      end
*
* ..psi''(cn) (2nd derivative) for complex argument
*      
c--------------
      subroutine psi2(cn,res)
      real*8 r
      complex*16 cn,y,y2,t,t0,res

      t=dcmplx(0.d0,0.d0)
2     r=dsqrt(dreal(cn)**2+dimag(cn)**2)
      if(r.gt.16.d0) goto 1
      t=t-2.d0/cn**3
      cn=cn+1.d0
      goto 2
1     y=1.d0/cn
      y2=y*y

      t0 = 
     & (-1.d0+(-1.d0+(-1.d0/2.d0+(1.d0/6.d0+(-1.d0/6.d0+(3.d0/10.d0
     & +(-5.d0/6.d0+(691.d0/210.d0+(-35.d0/2.d0+(3617.d0/30.d0
     & +(-43867.d0/42.d0+(1222277.d0/110.d0+(-854513.d0/6.d0
     & +(1181820455.d0/546-76977927.d0/2*y2)*y2)*y2)*y2)*y2)*y2)*y2)
     & *y2)*y2)*y2)*y2)*y2)*y)*y)*y2;

      res=t+t0

      return
      end
*
* ..psi'''(cn) (3rd derivative) for complex argument
*      
c--------------
      subroutine psi3(cn,res)
      real*8 r
      complex*16 cn,y,y2,t,t0,res

      t=dcmplx(0.d0,0.d0)
2     r=dsqrt(dreal(cn)**2+dimag(cn)**2)
      if(r.gt.16.d0) goto 1
      t=t+6.d0/cn**4
      cn=cn+1.d0
      goto 2
1     y=1.d0/cn
      y2=y*y

      t0 = 
     & (2.d0+(3.d0+(2.d0+(-1.d0+(4.d0/3.d0+(-3.d0+(10.d0+(-691.d0/15.d0
     & +(280.d0+(-10851.d0/5.d0+(438670.d0/21.d0+(-1222277.d0/5.d0
     & +(3418052.d0+(-1181820455.d0/21.d0+1077690978*y2)*y2)*y2)*y2)
     & *y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y)*y)*y2*y;

      res=t+t0

      return
      end

*
* ..psi''''(cn) (4th derivative) for complex argument
*      
c--------------
      subroutine psi4(cn,res)
      real*8 r
      complex*16 cn,y,y2,t,t0,res

      t=dcmplx(0.d0,0.d0)
2     r=dsqrt(dreal(cn)**2+dimag(cn)**2)
      if(r.gt.16.d0) goto 1
      t=t-24.d0/cn**5
      cn=cn+1.d0
      goto 2
1     y=1.d0/cn
      y2=y*y

      t0 = 
     & (-6.d0+(-12.d0+(-10.d0+(7.d0+(-12.d0+(33.d0+(-130.d0+(691.d0
     & +(-4760.d0+(206169.d0/5.d0+(-438670.d0+(28112371.d0/5.d0
     & +(-85451300.d0+10636384095.d0/7.d0*y2)*y2)*y2)*y2)*y2)*y2)
     & *y2)*y2)*y2)*y2)*y2)*y)*y)*y2*y2;

      res=t+t0

      return
      end

*
* ..psi'''''(cn) (5th derivative) for complex argument
*      
c--------------
      subroutine psi5(cn,res)
      real*8 r
      complex*16 cn,y,y2,t,t0,res

      t=dcmplx(0.d0,0.d0)
2     r=dsqrt(dreal(cn)**2+dimag(cn)**2)
      if(r.gt.16.d0) goto 1
      t=t+120.d0/cn**6
      cn=cn+1.d0
      goto 2
1     y=1.d0/cn
      y2=y*y

      t0 = 
     & (24.d0+(60.d0+(60.d0+(-56.d0+(120.d0+(-396.d0+(1820.d0
     & +(-11056.d0+(85680.d0+(-824676.d0+(9650740.d0
     & +(-674696904.d0/5.d0+(2221733800.d0-42545536380.d0*y2)*y2)*y2)
     & *y2)*y2)*y2)*y2)*y2)*y2)*y2)*y2)*y)*y)*y2*y2*y;

      res=t+t0

      return
      end

*
* ================================================================csm==












