*
*--------------------------------------------------------------------------
*
      subroutine zztoregion(n,x,j,a,b)
      implicit real*8(a-h,o-z)
*
      dimension x(n)
*
      if(j.gt.1) then
         do i=1,j-1
             if(x(i).lt.0.d0.or.x(i).gt.1.d0) then
                stop
             endif
         enddo
      endif
*
      a= 0.d0
      b= 1.d0
*
      return
      end
*
*-------------------------------------------------------------------
*
      subroutine zztocorrqcd(scal,zm,als,bqm2,cqm2,rbm2,rcm2,rsm2)
      implicit real*8(a-h,o-z)
*
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      data rz3/1.20205690315959428540d0/,nf/5/
*
      scal2= scal*scal
      bqm= sqrt(bqm2)
      cqm= sqrt(cqm2)
*
*-----computes the running charm mass
*-----computes the running bottom mass
*
      scal1g= 1.d0
      als1g= zztoralphas(zm,scal1g,als,nf)/pi
      alsc4= zztoralphas(zm,cqm,als,nf)/pi
      alsc42= alsc4*alsc4
      als1g2= als1g*als1g
      alsb= zztoralphas(zm,bqm,als,nf)/pi
      alsb2= alsb*alsb
      alsz= zztoralphas(zm,scal,als,nf)/pi
      alsz2= alsz*alsz
*
*-----first the running of the s-quark mass
*     up to c/b-threshold
*
      fn= 3.d0
      b03= (11.d0-2.d0/3.d0*fn)/4.d0
      b13= (102.d0-38.d0/3.d0*fn)/16.d0
      b23= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0*
     #     fn*fn)/64.d0
      g03= 1.d0
      g13= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0
      g23= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn-
     #      140.d0/81.d0*fn*fn)/64.d0
      rex3= g03/b03
      b03s= b03*b03
      b03c= b03s*b03
      b13s= b13*b13
      sm1g= 0.189d0
      alsdf= alsc4-als1g
      cfm13= g13/b03-b13*g03/b03s
      cfm23= g23/b03-b13*g13/b03s-b23*g03/b03s+
     #          b13s*g03/b03c
      rsmc= sm1g*(alsc4/als1g)**rex3*(1.d0+cfm13*
     #      alsdf+0.5d0*cfm13*cfm13*alsdf*alsdf+0.5d0*cfm23*
     #      (alsc42-als1g2))
*
      fn= 4.d0
      b04= (11.d0-2.d0/3.d0*fn)/4.d0
      b14= (102.d0-38.d0/3.d0*fn)/16.d0
      b24= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0*
     #      fn*fn)/64.d0
      g04= 1.d0
      g14= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0
      g24= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn-
     #      140.d0/81.d0*fn*fn)/64.d0
      rex4= g04/b04
      b04s= b04*b04
      b04c= b04s*b04
      b14s= b14*b14
      alsdf= alsb-alsc4
      cfm14= g14/b04-b14*g04/b04s
      cfm24= g24/b04-b14*g14/b04s-b24*g04/b04s+
     #       b14s*g04/b04c
      rsmb= rsmc*(alsb/alsc4)**rex4*(1.d0+cfm14*
     #      alsdf+0.5d0*cfm14*cfm14*alsdf*alsdf+0.5d0*cfm24*
     #      (alsb2-alsc42))
*
*-----c quark mass at c-threshold
*
      zero= 0.d0
      x1= 0.5d0
      x2= 2.0d0
      xacc= 1.d-12
      fn= 4.d0
      cmm4= zztorrunm(x1,x2,xacc,cqm,alsc4,rsmc,zero,fn)
*
*-----c quark mass at b-threshold
*
      cmb= cmm4*(alsb/alsc4)**rex4*(1.d0+cfm14*
     #     (alsb-alsc4)+0.5d0*cfm14*cfm14*
     #     (alsb-alsc4)**2+0.5d0*cfm24*(alsb2-alsc42))
*
*-----running c mass
*
      fn= 5.d0
      b05= (11.d0-2.d0/3.d0*fn)/4.d0
      b15= (102.d0-38.d0/3.d0*fn)/16.d0
      b25= (2857.d0/2.d0-5033.d0/18.d0*fn+325.d0/54.d0*
     #      fn*fn)/64.d0
      g05= 1.d0
      g15= (202.d0/3.d0-20.d0/9.d0*fn)/16.d0
      g25= (1249.d0-(2216.d0/27.d0+160.d0/3.d0*rz3)*fn-
     #      140.d0/81.d0*fn*fn)/64.d0
      rex5= g05/b05
      b05s= b05*b05
      b05c= b05s*b05
      b15s= b15*b15
      cfm15= g15/b05-b15*g05/b05s
      cfm25= g25/b05-b15*g15/b05s-b25*g05/b05s+b15s*g05/b05c
      rcm= cmb*(alsz/alsb)**rex5*(1.d0+cfm15*
     #     (alsz-alsb)+0.5d0*cfm15*cfm15*
     #     (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2))
      rcm2= rcm*rcm
*
*-----running s mass
*
      rsm= rsmb*(alsz/alsb)**rex5*(1.d0+cfm15*
     #     (alsz-alsb)+0.5d0*cfm15*cfm15*
     #     (alsz-alsb)**2+0.5d0*cfm25*(alsz2-alsb2))
      rsm2= rsm*rsm
*
*-----b quark mass at b-threshold
*
      x1= 0.5d0
      x2= 6.0d0
      xacc= 1.d-12
      fn= 5.d0
      bmm5= zztorrunm(x1,x2,xacc,bqm,alsb,cmb,rsmb,fn)
*
*-----running b mass
*
      alsdf= alsz-alsb
      rbm= bmm5*(alsz/alsb)**rex5*(1.d0+cfm15*
     #     alsdf+0.5d0*cfm15*cfm15*alsdf*alsdf+0.5d0*cfm25*
     #     (alsz2-alsb2))
      rbm2= rbm*rbm
*
      return
      end
*
*-------------------------------------------------------------------
*
      real*8 function zztoralphas(rs0,rs,als,nf)
      implicit real*8(a-h,o-z)
*
      common/zztmass/em,rmm,tm,uqm,dqm,cqm,sqm,bqm
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension b0(5),b1(5),b2(5)
*
*-----limits for lambda-5 are 1 mev < lambda_5 < 10 gev
*
      if(als.eq.0.d0) then
         zztoralphas= 0.d0
      else
         x1= 0.001d0
         x2= 10.0d0
         xacc= 1.d-12
         qcdl= zztoqcdlam(nf,als,rs0,x1,x2,xacc)
         pqcdl5= qcdl
         do i=1,5
             b0(i)= (11.d0-2.d0/3.d0*i)/4.d0
             b1(i)= (102.d0-38.d0/3.d0*i)/16.d0
             b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0-
     #              325.d0/27.d0*i))/64.d0
         enddo
*
         if(rs.lt.bqm) then
            rat= bqm/qcdl
            rl= 2.d0*log(rat)
            rll= log(rl)
            rb= log(b0(5)/b0(4))
            fac= b1(5)/b0(5)-b1(4)/b0(4)
            facp= b2(5)/b0(5)-b2(4)/b0(4)
            fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4)
            rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+
     #           b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*(
     #           fac2-facp-7.d0/72.d0)
            rhs= rhs/b0(4)
            ratl2= exp(rhs)
            qcdl= qcdl/sqrt(ratl2)      
            pqcdl4= qcdl
            nfe= nf-1
            if(rs.lt.cqm) then
               rat= cqm/qcdl
               rl= 2.d0*log(rat)
               rll= log(rl)
               rb= log(b0(4)/b0(3))
               fac= b1(4)/b0(4)-b1(3)/b0(3)
               facp= b2(4)/b0(4)-b2(3)/b0(3)
               fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3)
               rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+
     #              b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*(
     #              fac2-facp-7.d0/72.d0)
               rhs= rhs/b0(3)
               ratl2= exp(rhs)
               qcdl= qcdl/sqrt(ratl2)      
               pqcdl3= qcdl
               nfe= nf-2
            endif
         else
            nfe= nf
         endif
*
         qcdb0= 11.d0-2.d0/3.d0*nfe
         qcdb1= 102.d0-38.d0/3.d0*nfe
         qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+
     #          325.d0/27.d0*nfe*nfe)
         qcda= 2.d0*log(rs/qcdl)
*
         zztoralphas= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
     #            log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
     #            0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))
*
      endif
*
      return
      end
*
*--------------------------------------------------------------------
*
      real*8 function zztorals(qcdl,rs,nf)
      implicit real*8(a-h,o-z)
*
      common/zztmass/em,rmm,tm,uqm,dqm,cqm,sqm,bqm
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension b0(5),b1(5),b2(5)
*
      pqcdl5= qcdl
      do i=1,5
          b0(i)= (11.d0-2.d0/3.d0*i)/4.d0
          b1(i)= (102.d0-38.d0/3.d0*i)/16.d0
          b2(i)= 0.5d0*(2857.d0-i*(5033.d0/9.d0-
     #           325.d0/27.d0*i))/64.d0
      enddo
*
      if(rs.lt.bqm) then
         rat= bqm/qcdl
         rl= 2.d0*log(rat)
         rll= log(rl)
         rb= log(b0(5)/b0(4))
         fac= b1(5)/b0(5)-b1(4)/b0(4)
         facp= b2(5)/b0(5)-b2(4)/b0(4)
         fac2= b1(5)*b1(5)/b0(5)/b0(5)-b1(4)*b1(4)/b0(4)/b0(4)
         rhs= (b0(5)-b0(4))*rl+fac*rll-b1(4)/b0(4)*rb+
     #        b1(5)/b0(5)/b0(5)*fac*rll/rl+1.d0/b0(5)/rl*(
     #        fac2-facp-7.d0/72.d0)
         rhs= rhs/b0(4)
         ratl2= exp(rhs)
         qcdl= qcdl/sqrt(ratl2)      
         pqcdl4= qcdl
         nfe= nf-1
         if(rs.lt.cqm) then
            rat= cqm/qcdl
            rl= 2.d0*log(rat)
            rll= log(rl)
            rb= log(b0(4)/b0(3))
            fac= b1(4)/b0(4)-b1(3)/b0(3)
            facp= b2(4)/b0(4)-b2(3)/b0(3)
            fac2= b1(4)*b1(4)/b0(4)/b0(4)-b1(3)*b1(3)/b0(3)/b0(3)
            rhs= (b0(4)-b0(3))*rl+fac*rll-b1(3)/b0(3)*rb+
     #           b1(4)/b0(4)/b0(4)*fac*rll/rl+1.d0/b0(4)/rl*(
     #           fac2-facp-7.d0/72.d0)
            rhs= rhs/b0(3)
            ratl2= exp(rhs)
            qcdl= qcdl/sqrt(ratl2)      
            pqcdl3= qcdl
            nfe= nf-2
         endif
      else
         nfe= nf
      endif
*
      qcdb0= 11.d0-2.d0/3.d0*nfe
      qcdb1= 102.d0-38.d0/3.d0*nfe
      qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nfe+
     #       325.d0/27.d0*nfe*nfe)
      qcda= 2.d0*log(rs/qcdl)
*
      zztorals= 4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
     #         log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
     #         0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0))
*
      return
      end
*
*----------------------------------------------------------
*
      real*8 function zztorrunm(x1,x2,xacc,qm,als,rm1,rm2,fn)
      implicit real*8(a-h,o-z)
*
      parameter (jmax=50)
*
      fmid= zztoqcdmass(qm,als,rm1,rm2,fn,x2)
      f= zztoqcdmass(qm,als,rm1,rm2,fn,x1)
      if (f*fmid.ge.0.d0) then
          print*,'root must be bracketed for bisection'
          print 1,qm
    1     format(/' error detected by rrunm ',/
     #            ' current value of quark mass  = ',e20.5)
          stop
      endif
      if (f.lt.0.d0) then
          zztorrunm= x1
          dx= x2-x1
      else
          zztorrunm= x2
          dx= x1-x2
      endif
      do j=1,jmax
          dx= dx*0.5d0
          xmid= zztorrunm+dx
          fmid= zztoqcdmass(qm,als,rm1,rm2,fn,xmid)
          if (fmid.le.0.d0) zztorrunm= xmid
          if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return
      enddo
      pause 'too many bisections'
      end       
*
*-----zztoqcdmass------------------------------------------------------
*
      real*8 function zztoqcdmass(qm,als,rm1,rm2,fn,x)
      implicit real*8(a-h,o-z)
*
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      rz2= pis/6.d0
      rz3= 1.20205690315959428540d0
      ra4= 0.5174790617d0
      rz5= 1.03692775514336992633d0
      rln= 2.d0*log(qm/x)
      rlns= rln*rln
      r1= rm1/x
      r2= rm2/x
      delta0= 3.d0/4.d0*rz2-3.d0/8.d0
      delta1= r1*(pis/8.d0+r1*(-0.597d0+0.230d0*r1))
      delta2= r2*(pis/8.d0+r2*(-0.597d0+0.230d0*r2))
*
      rhs= 1.d0+als*(4.d0/3.d0+rln+als*(3817.d0/288.d0-8.d0/3.d0+
     #     2.d0/3.d0*(2.d0+log(2.d0))*rz2-rz3/6.d0-fn/3.d0*(rz2+
     #     71.d0/48.d0)+4.d0/3.d0*(delta0+delta1+delta2)+(173.d0/
     #     24.d0-13.d0/36.d0*fn)*rln+(15.d0/8.d0-fn/12.d0)*rlns))
      zztoqcdmass= qm-x*rhs
*
      return
      end
*
*---------------------------------------------------------------------
*
      real*8 function zztoqcdlam(nf,als,rs,x1,x2,xacc)
      implicit real*8(a-h,o-z)
*
      parameter (jmax=50,nout=21)
*
      fmid= zztoqcdscale(nf,als,rs,x2)
      f= zztoqcdscale(nf,als,rs,x1)
      if (f*fmid.ge.0.d0) then
          print*,'root must be bracketed for bisection'
          print 1,als,x1,x2
    1     format(/' error detected by zztoqcdlam ',/
     #            ' current value of alpha_s  = ',e20.5,/
     #            ' interval   ',2e20.5)
          stop
      endif
      if (f.lt.0.d0) then
          zztoqcdlam= x1
          dx= x2-x1
      else
          zztoqcdlam= x2
          dx= x1-x2
      endif
      do j=1,jmax
          dx= dx*0.5d0
          xmid= zztoqcdlam+dx
          fmid= zztoqcdscale(nf,als,rs,xmid)
          if (fmid.le.0.d0) zztoqcdlam= xmid
          if(abs(dx).lt.xacc.or.fmid.eq.0.d0) return
      enddo
      pause 'too many bisections'
      end       
*
*------------------------------------------------------------------
*
      real*8 function zztoqcdscale(nf,als,rs,x)
      implicit real*8(a-h,o-z)
*
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      qcdb0= 11.d0-2.d0/3.d0*nf
      qcdb1= 102.d0-38.d0/3.d0*nf
      qcdb2= 0.5d0*(2857.d0-5033.d0/9.d0*nf+325.d0/27.d0*nf*nf)
      qcda= 2.d0*log(rs/x)
      zztoqcdscale= als-(4.d0*pi/qcdb0/qcda*(1.d0-qcdb1/qcdb0**2/qcda*
     #            log(qcda)+(qcdb1/qcdb0**2/qcda)**2*((log(qcda)-
     #            0.5d0)**2+qcdb2*qcdb0/qcdb1**2-5.d0/4.d0)))
*
      return
      end
*
*-------------------------------------------------------------
*
      subroutine zztofst(nt,st,fvect,iflag)
      implicit real*8(a-h,o-z)
*
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
      common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2
*
      dimension zs(2),zl(2)
      dimension st(nt),fvect(nt)
      dimension pggfz(2),pgglqz(2),p3qz(2),fzz(2),p3qw(2),fw(2),
     #          fzw(2),fwz(2)
*
      tqm2= st(1)
*
*-----alpha(zm)
*
      wm2= wm*wm
      zm2= zm*zm
      pzr= -zm2
      pzi= 0.d0
      call zztopoleg(pzr,pzi,pggfz,pgglqz,pggnpz)
      call zztopolez0(pggf0,fw0)
      x= alphai-0.25d0*(pggfz(1)-pggf0+pggnpz)/pi
      y= -0.25d0*(pggfz(2)+pgglqz(2))/pi
      alzr= x/(x*x+y*y)
      alzi= -y/(x*x+y*y)
*
*-----re(1/g^2(zm))
*
      apis= 16.d0*pis
      p2zr= -zm2
      p2zi= 0.d0
      call zztopole(p2zr,p2zi,p3qz,fzz,fwz)
      bx= 1.d0/8.d0/gf/zm2+(fw0-fzz(1))/zm2/apis
      xih= -p3qz(2)/apis
      ac= 4.d0*pi*alzr
      bc= -1.d0-8.d0*pi*xih*alzi
      cc= -4.d0*pi*xih*xih*alzr+bx
      ifail= 0
      call c02ajf(ac,bc,cc,zs,zl,ifail)
      g2z= zs(1)
*
      p2r= -wm2
      p2i= 0.d0
      call zztopole(p2r,p2i,p3qw,fzw,fw)
*
      fvect(1)= 1.d0+0.5d0*gf*wm2/pis*((fw0-fw(1))/wm2+p3qw(1))-
     #          8.d0*gf*wm2*(g2z+p3qz(1)/apis)
*
      return
      end
*
*--------------------------------------------------------------------
* 
      subroutine zztopoleg(p2r,p2i,pggf,pgglq,pggnp)  
      implicit real*8(a-h,o-z)  
*   
      common/zzted/eps,delta
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
      common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2
*
      dimension pggf(2),pgglq(2)
      dimension bfl(2),bfh(2),bflq(2)  
      dimension b0e(2),b1e(2),b21e(2)
      dimension b0m(2),b1m(2),b21m(2)
      dimension b0t(2),b1t(2),b21t(2)
      dimension b0h(2),b1h(2),b21h(2)
      dimension b0lq(2),b1lq(2),b21lq(2)
*
      call zztosbff(p2r,p2i,em2,em2,b0e,b1e,b21e)
      call zztosbff(p2r,p2i,rmm2,rmm2,b0m,b1m,b21m)
      call zztosbff(p2r,p2i,tm2,tm2,b0t,b1t,b21t)
      call zztosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h)
      call zztosbff(p2r,p2i,uqm2,uqm2,b0lq,b1lq,b21lq)
*
      do i=1,2
          bfl(i)= 2.d0*b21e(i)-b0e(i)+
     #            2.d0*b21m(i)-b0m(i)+
     #            2.d0*b21t(i)-b0t(i)   
          bfh(i)= 2.d0*b21h(i)-b0h(i)
          bflq(i)= 2.d0*b21lq(i)-b0lq(i)
      enddo
*
      do i=1,2
       pggf(i)= 4.d0*bfl(i)+16.d0/3.d0*bfh(i)
       pgglq(i)= 44.d0/3.d0*bflq(i)
      enddo  
*
      ap2x= abs(p2r)   
      px= sqrt(ap2x)
      if(px.lt.0.3d0) then
         a= 0.d0
         b= 0.00835d0
         c= 1.d0
      else if(px.lt.3.d0) then
         a= 0.d0
         b= 0.00238d0
         c= 3.927d0
      else if(px.lt.100.d0) then
         a= 0.00165d0
         b= 0.00299d0
         c= 1.d0
      else if(px.gt.100.d0) then
         a= 0.00221d0
         b= 0.00293d0
         c= 1.d0
      endif
*
      pggnp= 4.d0*pi*alphai*(a+b*log(1.d0+c*ap2x))
*
      return
      end
*
*-------------------------------------------------------------------------
*
      subroutine zztopolez0(pggf0,fw0)
      implicit real*8(a-h,o-z)
*
      common/zzted/eps,delta
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
      common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2
*
      call zztorbff0(em2,em2,b0ee0,b1ee0,b21ee0)
      call zztorbff0(rmm2,rmm2,b0mm0,b1mm0,b21mm0)
      call zztorbff0(tm2,tm2,b0tau0,b1tau0,b21tau0)
      call zztorbff0(tqm2,tqm2,b0tt0,b1tt0,b21tt0)
*
      bfe= 2.d0*b21ee0-b0ee0
      bfm= 2.d0*b21mm0-b0mm0
      bftau= 2.d0*b21tau0-b0tau0
      bft= 2.d0*b21tt0-b0tt0
*
      pggf0= 4.d0*(bfe+bfm+bftau)+16.d0/3.d0*bft
*
      call zztorbff0(tqm2,dmy2,b0tb0,b1tb0,b21tb0)
*
      fw0= -3.d0*tqm2*(b0tb0+b1tb0)
*
      return
      end
*
*-------------------------------------------------------------------
*
      subroutine zztopole(p2r,p2i,p3q,fz,fw)
      implicit real*8(a-h,o-z)
*
      common/zzted/eps,delta
      common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2
*
      dimension bfl(2),bfh(2)
      dimension p3q(2),fz(2),fw(2)
      dimension b0l(2),b1l(2),b21l(2)
      dimension b0h(2),b1h(2),b21h(2)
      dimension b0c(2),b1c(2),b21c(2)
*
      call zztosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l)
      call zztosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h)
      call zztosbff(p2r,p2i,tqm2,dmy2,b0c,b1c,b21c)
*
      do i=1,2
          bfl(i)= 2.d0*b21l(i)-b0l(i)
          bfh(i)= 2.d0*b21h(i)-b0h(i)
      enddo
      bfsr= bfh(1)-bfl(1)
      bfsi= bfh(2)-bfl(2)
*
      do i=1,2
       p3q(i)= 10.d0*bfl(i)+2.d0*bfh(i)
      enddo
*
      fz(1)= -0.5d0*(p2r*bfsr-p2i*bfsi)-3.d0/2.d0*tqm2*b0h(1)
      fz(2)= -0.5d0*(p2r*bfsi+p2i*bfsr)-3.d0/2.d0*tqm2*b0h(2)
*      
      auxr= 3.d0*(b21c(1)+b1c(1))-2.d0*b21h(1)+b0h(1)-b21l(1)-
     #      3.d0*b1l(1)-b0l(1)
      auxi= 3.d0*(b21c(2)+b1c(2))-2.d0*b21h(2)+b0h(2)-b21l(2)-
     #      3.d0*b1l(2)-b0l(2)
*
      fw(1)= 2.d0*(p2r*auxr-p2i*auxi)-3.d0*tqm2*(b1c(1)+b0c(1))
      fw(2)= 2.d0*(p2r*auxi+p2i*auxr)-3.d0*tqm2*(b1c(2)+b0c(2))
*
      return
      end
*
*------------------------------------------------------------------
*
      subroutine zztofsw(n,sw,fvecw,iflag)
      implicit real*8(a-h,o-z)
*
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension sw(n),fvecw(n)
      dimension p3qw(2),fzw(2),fw(2),p3q(2),fz(2),fcw(2)
*
      wm2= wm*wm
      swr= sw(1)
      swi= sw(2)
      wmu= sqrt(swr)
      wlg= -swi/wmu
*
      call zztopolez0(pggf0,fw0)
      p2r= -wm2
      p2i= 0.d0
      call zztopole(p2r,p2i,p3qw,fzw,fw)
      p2r= -swr
      p2i= -swi
      call zztopole(p2r,p2i,p3q,fz,fcw)
*
      x0= 1.d0/8.d0/gf/wm2
      xr= 1.d0+0.5d0*gf*wm2/pis*((fw0-fw(1))/wm2+p3qw(1)-p3q(1))
      xi= -0.5d0*gf*wm2/pis*p3q(2)
      xr= x0*xr
      xi= x0*xi
*
      fvecw(1)= gf*(swr*xr-swi*xi+(fcw(1)-fw0)/16.d0/pis)
     #          -1.d0/8.d0
      fvecw(2)= swr*xi+swi*xr+fcw(2)/16.d0/pis
*
      return
      end
*
*----------------------------------------------------------
*
      subroutine zztofsz(n,sz,fvecz,iflag)
      implicit real*8(a-h,o-z)
*
      common/zztafsz/g2z
      common/zztiaz/alzr,alzi,alzir,alzii
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension sz(n),fvecz(n)
      dimension p3qz(2),fzz(2),p3q(2),fz(2),fwz(2),fw(2)
      dimension pggfz(2),pgglqz(2),pggf(2),pgglq(2),evl(2)
*
      szr= sz(1)
      szi= sz(2)
*
      zm2= zm*zm
      apis= 16.d0*pis
      call zztopolez0(pggf0,fw0)
      pzr= -zm2
      pzi= 0.d0
      call zztopoleg(pzr,pzi,pggfz,pgglqz,pggnpz)
      call zztopole(pzr,pzi,p3qz,fzz,fwz)
      p2r= -szr
      p2i= -szi
      call zztopoleg(p2r,p2i,pggf,pgglq,pggnp)
      call zztopole(p2r,p2i,p3q,fz,fw)
      do i=1,2
          evl(i)= pggfz(i)-pggf(i)+pgglqz(i)-pgglq(i)
      enddo
      axr= alzir+0.25d0*evl(1)/pi
      axi= alzii+0.25d0*evl(2)/pi
      alr= axr/(axr*axr+axi*axi)
      ali= -axi/(axr*axr+axi*axi)
*
      xr= g2z+(p3qz(1)-p3q(1))/apis
      xi= -p3q(2)/apis
*
      yr= xr-4.d0*pi*(alr*(xr*xr-xi*xi)-2.d0*ali*xr*xi)
      yi= xi-4.d0*pi*(2.d0*alr*xr*xi+ali*(xr*xr-xi*xi))
*
      fvecz(1)= 8.d0*gf*(szr*yr-szi*yi+(fz(1)-fw0)/apis)-1.d0
      fvecz(2)= szr*yi+szi*yr+fz(2)/apis
*
      return
      end
*
*-------------------------------------------------------------------------
*
      subroutine zztosbff(p2r,p2i,rm12,rm22,b0,b1,b21)
      implicit real*8(a-h,o-z)
*
      common/zzted/eps,delta
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension cmp2(2)
      dimension clnmp2(2)
      dimension b0(2),b1(2),b21(2)
*
      if(rm12.ne.0.d0.and.rm12.eq.rm22.and.
     #   (p2r.gt.0.d0.and.sqrt(p2r).lt.(2.d-5*sqrt(rm12)))) then 
       call zztorbff0(rm12,rm22,ab0,ab1,ab21)
       b0(1)= ab0-p2r/6.d0/rm12
       b1(1)= ab1+p2r/12.d0/rm12
       b21(1)= ab21-p2r/20.d0/rm12
       b0(2)= 0.d0
       b1(2)= 0.d0
       b21(2)= 0.d0
       return
      endif
*
      p2m= p2r*p2r+p2i*p2i
*
      if(rm12.eq.0.d0.and.rm22.eq.0.d0) then
         if(p2i.eq.0) then
            cmp2(1)= p2r
            cmp2(2)= -eps
            call zztoclnx(cmp2,clnmp2)
            b0(1)= delta+2.d0-clnmp2(1)
            b0(2)= -clnmp2(2)
         else
            cmp2(1)= p2r
            cmp2(2)= p2i
            call zztoclnx(cmp2,clnmp2)
            b0(1)= delta+2.d0-clnmp2(1)
            b0(2)= -clnmp2(2)+2.d0*pi
         endif   
         b1(1)= -0.5d0*b0(1)
         b1(2)= -0.5d0*b0(2)
         b21(1)= 1.d0/18.d0+1.d0/3.d0*b0(1)
         b21(2)= 1.d0/3.d0*b0(2)
      else if(rm22.eq.0.d0.and.rm12.ne.0.d0) then
         a0= rm12*(-delta-1.d0+log(rm12))
         if(p2i.eq.0) then
            cmp2(1)= 1.d0+p2r/rm12  
            cmp2(2)= -eps/rm12 
            call zztoclnx(cmp2,clnmp2) 
            b0(1)= delta+2.d0-log(rm12)-(1.d0+rm12/p2r)*clnmp2(1) 
            b0(2)= -(1.d0+rm12/p2r)*clnmp2(2)   
            b1(1)= 0.5d0*(a0+(rm12-p2r)*b0(1))/p2r
            b1(2)= 0.5d0*(rm12-p2r)*b0(2)/p2r
            b21(1)= (3.d0*rm12+p2r)/18.d0/p2r+(rm12-p2r)*a0/3.d0/
     #              (p2r*p2r)+(1.d0-rm12/p2r*(1.d0-rm12/p2r))*
     #              b0(1)/3.d0
            b21(2)= (1.d0-rm12/p2r*(1.d0-rm12/p2r))*b0(2)/3.d0
         else 
            cmp2(1)= 1.d0+p2r/rm12  
            cmp2(2)= p2i/rm12
            call zztoclnx(cmp2,clnmp2) 
            xr= p2r/p2m
            xi= -p2i/p2m  
            fctr= 1.d0+rm12*xr  
            fcti= rm12*xi
            b0(1)= delta+2.d0-log(rm12)-fctr*clnmp2(1)+fcti*clnmp2(2)
            b0(2)= -fctr*clnmp2(2)-fcti*clnmp2(1)
            if(p2r.lt.(-rm12)) then
               b0(1)= b0(1)+2.d0*pi*rm12*p2i/p2m
               b0(2)= b0(2)+2.d0*pi*(1.d0+rm12*p2r/p2m)   
            endif
            x2r= xr*xr-xi*xi
            x2i= 2.d0*xr*xi
            b1(1)= 0.5d0*a0*xr-0.5d0*b0(1)+0.5d0*rm12*(xr*b0(1)-
     #             xi*b0(2))
            b1(2)= 0.5d0*a0*xi-0.5d0*b0(2)+0.5d0*rm12*(xr*b0(2)+
     #             xi*b0(1))
            b21(1)= 3.d0/18.d0*rm12*xr+1.d0/18.d0+a0/3.d0*
     #              (rm12*x2r-xr)+1.d0/3.d0*((1.d0-rm12*xr+
     #              rm12*rm12*x2r)*b0(1)-rm12*(-xi+rm12*x2i)*b0(2))
            b21(2)= 3.d0/18.d0*rm12*xi+a0/3.d0*
     #              (rm12*x2i-xi)+1.d0/3.d0*((1.d0-rm12*xr+
     #              rm12*rm12*x2r)*b0(2)+rm12*(-xi+rm12*x2i)*b0(1))
         endif
      else if(rm12.ne.0.d0.and.rm12.eq.rm22) then
         a0= rm12*(-delta-1.d0+log(rm12))
         if(p2i.eq.0) then  
            p2d= p2r*p2r+eps*eps
            b2r= 1.d0+4.d0*rm12*p2r/p2d 
            b2i= 4.d0*rm12*eps/p2d
            call a02aaf(b2r,b2i,br,bi)
            cmp2(1)= (br*br+bi*bi-1.d0)/((br-1.d0)*(br-1.d0)+bi*bi)
            cmp2(2)= -2.d0*bi/((br-1.d0)*(br-1.d0)+bi*bi) 
            call zztoclnx(cmp2,clnmp2) 
            b0(1)= delta+2.d0-log(rm12)-(br*clnmp2(1)-   
     #             bi*clnmp2(2))
            b0(2)= -(br*clnmp2(2)+bi*clnmp2(1))          
            b1(1)= -0.5d0*b0(1)
            b1(2)= -0.5d0*b0(2)
            b21(1)= (6.d0*rm12+p2r)/18.d0/p2r+a0/3.d0/p2r+
     #              1.d0/3.d0*(1.d0+rm12/p2r)*b0(1)
            b21(2)= 1.d0/3.d0*(1.d0+rm12/p2r)*b0(2)
         else   
            xr= p2r/p2m
            xi= -p2i/p2m
            b2r= 1.d0+4.d0*rm12*xr 
            b2i= 4.d0*rm12*xi 
            call a02aaf(b2r,b2i,br,bi)
            cmp2(1)= (br*br+bi*bi-1.d0)/((br-1.d0)*(br-1.d0)+bi*bi) 
            cmp2(2)= -2.d0*bi/((br-1.d0)*(br-1.d0)+bi*bi) 
            call zztoclnx(cmp2,clnmp2)
            b0(1)= delta+2.d0-log(rm12)-(br*clnmp2(1)-
     #             bi*clnmp2(2))
            b0(2)= -(br*clnmp2(2)+bi*clnmp2(1)) 
            if(p2r.lt.(-4.d0*rm12)) then
               b0(1)= b0(1)-2.d0*pi*bi
               b0(2)= b0(2)+2.d0*pi*br  
            endif 
            x2r= xr*xr-xi*xi
            x2i= 2.d0*xr*xi
            b1(1)= -0.5d0*b0(1)
            b1(2)= -0.5d0*b0(2)
            b21(1)= rm12*xr/3.d0+1.d0/18.d0+a0/3.d0*xr+
     #              1.d0/3.d0*((1.d0+rm12*xr)*b0(1)-rm12*xi*
     #              b0(2))
            b21(2)= rm12*xi/3.d0+a0/3.d0*xi+1.d0/3.d0*((1.d0+
     #              rm12*xr)*b0(2)+rm12*xi*b0(1))
         endif
      else
         a01= rm12*(-delta-1.d0+log(rm12))
         a02= rm22*(-delta-1.d0+log(rm22)) 
         if(p2i.eq.0) then
            b2r= (p2r-rm22+rm12)*(p2r-rm22+rm12)- 
     #           4.d0*rm12*rm22-eps*eps 
            b2i= -2.d0*(p2r+rm22+rm12)*eps
            call a02aaf(b2r,b2i,br,bi)
            cmp2(1)= 0.5d0*(-p2r-rm22-rm12+br)/sqrt(rm22*rm12)
            cmp2(2)= 0.5d0*(eps+bi)/sqrt(rm22*rm12) 
            call zztoclnx(cmp2,clnmp2)
            crr= (br*clnmp2(1)-bi*clnmp2(2))/p2r
            cri= (br*clnmp2(2)+bi*clnmp2(1))/p2r
            b0(1)= delta-crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/
     #             p2r*log(rm22/rm12)+2.d0  
            b0(2)= -cri 
            b1(1)= 0.5d0/p2r*(a01-a02+(rm12-rm22-p2r)*b0(1))
            b1(2)= 0.5d0/p2r*(rm12-rm22-p2r)*b0(2)
            b21(1)= (3.d0*(rm12+rm22)+p2r)/18.d0/p2r+
     #              (rm12-rm22-p2r)*a01/3.d0/(p2r*p2r)+
     #              (rm22-rm12+2.d0*p2r)*a02/3.d0/(p2r*p2r)+
     #              (1.d0+1.d0/p2r*(2.d0*rm22-rm12+(rm22-rm12)**2/
     #              p2r))*b0(1)
            b21(2)= (1.d0+1.d0/p2r*(2.d0*rm22-rm12+(rm22-rm12)**2/
     #              p2r))*b0(2)
          else
            xr= p2r/p2m
            xi= -p2i/p2m
            x2r= xr*xr-xi*xi
            x2i= 2.d0*xr*xi
            b2r= (p2r-rm22+rm12)*(p2r-rm22+rm12)-
     #           4.d0*rm12*rm22-p2i*p2i 
            b2i= 2.d0*(p2r+rm22+rm12)*p2i
            call a02aaf(b2r,b2i,br,bi)
            cmp2(1)= 0.5d0*(-p2r-rm22-rm12+br)/sqrt(rm22*rm12)
            cmp2(2)= 0.5d0*(-p2i+bi)/sqrt(rm22*rm12)
            call zztoclnx(cmp2,clnmp2)
            crr= (br*clnmp2(1)-bi*clnmp2(2))*xr-
     #           (br*clnmp2(2)+bi*clnmp2(1))*xi
            cri= (br*clnmp2(1)-bi*clnmp2(2))*xi+
     #           (br*clnmp2(2)+bi*clnmp2(1))*xr
            b0(1)= delta+crr-0.5d0*log(rm22*rm12)+0.5d0*(rm22-rm12)/
     #             p2r*log(rm22/rm12)+2.d0
            b0(2)= cri
            tsti= -(sqrt(rm12)+sqrt(rm22))*(sqrt(rm12)+sqrt(rm22))  
            if(p2r.lt.tsti) then
               b0(1)= b0(1)-2.d0*pi*(br*xi+bi*xr) 
               b0(2)= b0(2)+2.d0*pi*(br*xr-bi*xi)
            endif
            b1(1)= 0.5d0*xr*(a01-a02+(rm12-rm22)*b0(1))-
     #             0.5d0*xi*(rm12-rm22)*b0(2)-0.5d0*b0(1)
            b1(2)= 0.5d0*xi*(a01-a02+(rm12-rm22)*b0(1))+
     #             0.5d0*xr*(rm12-rm22)*b0(2)-0.5d0*b0(2)
            b21(1)= (rm12+rm22)/6.d0*xr+1.d0/18.d0+
     #              (rm12-rm22)/3.d0*x2r*a01-a01/3.d0*xr+
     #              (rm22-rm12)/3.d0*a02*x2r+2.d0/3.d0*a02*xr+
     #              b0(1)/3.d0+(2.d0*rm22-rm12)/3.d0*(xr*b0(1)-
     #              xi*b0(2))+(rm22-rm12)**2/3.d0*(x2r*b0(1)-x2i*b0(2))
            b21(2)= (rm12+rm22)/6.d0*xi+
     #              (rm12-rm22)/3.d0*x2i*a01-a01/3.d0*xi+
     #              (rm22-rm12)/3.d0*a02*x2i+2.d0/3.d0*a02*xi+
     #              b0(2)/3.d0+(2.d0*rm22-rm12)/3.d0*(xr*b0(2)+
     #              xi*b0(1))+(rm22-rm12)**2/3.d0*(x2r*b0(2)+x2i*b0(1))
         endif
      endif
*
      return
      end
*
*-------------------------------------------------------------------
*
      subroutine zztorbff0(rm12,rm22,b0,b1,b21)
      implicit real*8(a-h,o-z)
*
      common/zzted/eps,delta
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension arg(2),cln(2),fr(3)
*
      if(rm12.eq.0.d0) then
         arm12= 1.d-20
      else 
         arm12= rm12
      endif
      if(rm22.eq.0.d0) then
         arm22= 1.d-20
      else 
         arm22= rm22
      endif
      if (arm12.eq.arm22) then
          fact= delta-log(arm12)
          b0= fact
          b1= -0.5d0*fact
          b21= fact/3.d0
          return
      else
          n= 3
          yr= arm12/(arm12-arm22)
          omyr= -arm22/(arm12-arm22)
          yi= -eps/(arm12-arm22)
          call zztorcg(n,yr,yi,omyr,fr)
          arg(1)= omyr
          arg(2)= -yi
          call zztoclnx(arg,cln)
          f1r= fr(1)-cln(1)
          f2r= 2.d0*fr(2)-cln(1)
          f3r= 3.d0*fr(3)-cln(1)
          fact= delta-log(arm22)
          b0= fact-f1r
          b1= 0.5d0*(-fact+f2r)
          b21= 1.d0/3.d0*(fact-f3r)
          return
      endif
      end 
*
*----------------------------------------------------------
*
      subroutine zztoclnx(arg,res)
      implicit real*8(a-h,o-z)
*
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension arg(2),aarg(2),res(2)
*
      do i= 1,2
          aarg(i)= abs(arg(i))
      enddo
      azm2= (arg(1))**2+(arg(2))**2
      azm= sqrt(azm2)
      res(1)= log(azm)
      if (arg(1).eq.0.d0) then
          if (arg(2).gt.0.d0) then
              teta= pi/2.d0
          else
              teta= -pi/2.d0
          endif
          res(2)= teta
          return
      else if (arg(2).eq.0.d0) then 
               if (arg(1).gt.0.d0) then
                   teta= 0.d0
               else
                   teta= pi
               endif
          res(2)= teta
          return
      else
          tnteta= aarg(2)/aarg(1)
          teta= atan(tnteta)
          sr= arg(1)/aarg(1)
          si= arg(2)/aarg(2)
          if (sr.gt.0.d0) then
              res(2)= si*teta
          else
              res(2)= si*(pi-teta)
          endif
          return
      endif
      end
*
*---------------------------------------------------------
*
      subroutine zztoclnomx(arg,omarg,res)
      implicit real*8(a-h,o-z)
*
      dimension arg(2),omarg(2),res(2),ares(2),ct(10),sn(10)
*
      zr= arg(1)
      zi= arg(2)
      azm2= zr*zr+zi*zi
      azm= sqrt(azm2)
      if (azm.lt.1.d-7) then
          ct(1)= zr/azm
          sn(1)= zi/azm
          do n=2,10
              ct(n)= ct(1)*ct(n-1)-sn(1)*sn(n-1)
              sn(n)= sn(1)*ct(n-1)+ct(1)*sn(n-1)
          enddo
          ares(1)= ct(10)/10.d0
          ares(2)= sn(10)/10.d0
          do k=9,1,-1
              ares(1)= ares(1)*azm+ct(k)/k
              ares(2)= ares(2)*azm+sn(k)/k
          enddo
          ares(1)= -ares(1)*azm
          ares(2)= -ares(2)*azm
      else
          call zztoclnx(omarg,ares)
      endif
      do i= 1,2
          res(i)= ares(i)
      enddo
*
      return
      end
*
*-----------------------------------------------------
*
      subroutine zztorcg(n,zr,zi,omzr,gfr)
      implicit real*8(a-h,o-z)
*
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      dimension gfr(n)
      dimension a(4),ra(16),z(2),omz(2),oz(2),clnomz(2),clnoz(2)
      dimension ca(2,10),aux(2),zp(2),ct(16),sn(16)
      data ra/-0.2d0,0.666666666666667d-1,
     #            -0.952380952380952d-2,-0.396825396825397d-3,
     #             0.317460317460317d-3,-0.132275132275132d-4,
     #            -0.962000962000962d-5,0.105218855218855d-5,
     #             0.266488361726450d-6,-0.488745528428000d-7,
     #            -0.675397500794000d-8,0.190720263471000d-8,
     #             0.153663007690000d-9,-0.679697905790000d-10,
     #            -0.293683556000000d-11,0.228836696000000d-11/
*
      z(1)= zr
      z(2)= zi
      omz(1)= omzr
      omz(2)= -zi
      if (zr.eq.0.d0.and.zi.eq.0.d0) then
          do k=1,n
             gfr(k)= -1.d0/k/k
          enddo
          return
      else if (zr.eq.1.d0.and.zi.eq.0.d0) then
               a(1)= -1.d0
               do j=2,4
                   a(j)= ((j-1.d0)*a(j-1)-1.d0/j)/j
              enddo
              do k=1,n
                 gfr(k)= a(k)
              enddo
              return
      else       
          zmod2= zr*zr+zi*zi
          zmod= sqrt(zmod2)
          call zztoclnx(omz,clnomz)
* 
*        |z| < 4        
*      
         if (zmod.lt.4.d0) then
             oz(1)= -z(1)
             oz(2)= -z(2)      
             call zztoclnx(oz,clnoz)
             ca(1,1)= omz(1)*clnomz(1)-omz(2)*clnomz(2)+z(1)*clnoz(1)-
     #                z(2)*clnoz(2)-1.d0
             ca(2,1)= omz(1)*clnomz(2)+omz(2)*clnomz(1)+z(1)*clnoz(2)+
     #                z(2)*clnoz(1)
            if (n.eq.1) then
                gfr(1)= ca(1,1)
                return
            else
                do j= 2,n
                  jm1= j-1
                  ca(1,j)= ((j-1.d0)*(z(1)*ca(1,jm1)-z(2)*ca(2,jm1))+
     #                     omz(1)*clnomz(1)-omz(2)*clnomz(2)-1.d0/j)/j
                  ca(2,j)= ((j-1.d0)*(z(1)*ca(2,jm1)+z(2)*ca(1,jm1))+
     #                     omz(1)*clnomz(2)+omz(2)*clnomz(1))/j
                enddo
                do k=1,n
                    gfr(k)= ca(1,k)
                enddo
                return
            endif
*
*        |z| > 4
*
      else
         aux(1)= (-zr*omzr+zi**2)/zmod2
         aux(2)= zi/zmod2
         call zztoclnx(aux,zp)
         zpm2= zp(1)*zp(1)+zp(2)*zp(2)
         zpm= sqrt(zpm2)
         ct(1)= zp(1)/zpm
         sn(1)= zp(2)/zpm
         do k=2,16
             ct(k)= ct(1)*ct(k-1)-sn(1)*sn(k-1)
             sn(k)= sn(1)*ct(k-1)+ct(1)*sn(k-1)
         enddo
         ca(1,4)= ra(16)*ct(16)*zpm+ra(15)*ct(15)
         ca(2,4)= ra(16)*sn(16)*zpm+ra(15)*sn(15)
         do j=14,1,-1
             ca(1,4)= ca(1,4)*zpm+ra(j)*ct(j)
             ca(2,4)= ca(2,4)*zpm+ra(j)*sn(j)
         enddo
         ca(1,4)= ca(1,4)*zpm
         ca(2,4)= ca(2,4)*zpm
         do j= 3,1,-1
             jp1= j+1
             ca(1,j)= ((ca(1,jp1)+1.d0/jp1)*z(1)+ca(2,jp1)*z(2))/zmod2
             ca(2,j)= (ca(2,jp1)*z(1)-(ca(1,jp1)+1.d0/jp1)*z(2))/zmod2
         enddo
         do k=1,n
             gfr(k)= (ca(1,k)+clnomz(1))/k
         enddo
         return
      endif
      endif
      end
*
*----------------------------------------------------------
*
      subroutine zztocfl(p2r,p2i,p3qf,s3qb,pgg,s33,sww,r1,r2,r3)
      implicit real*8(a-h,o-z)
*
      common/zzted/eps,delta
      common/zttvbm/zm,wm,hm,zm2,wm2,hm2,zg,rs,alsz
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
      common/zztfmass2/em2,rmm2,tm2,uqm2,dqm2,cqm2,sqm2,bqm2,tqm2,dmy2
*
      dimension bfl(2),bfh(2),co(2)
      dimension b0l(2),b1l(2),b21l(2)
      dimension b0h(2),b1h(2),b21h(2)
      dimension b0c(2),b1c(2),b21c(2)
      dimension b0w(2),b1w(2),b21w(2)
      dimension b0z(2),b1z(2),b21z(2)
      dimension b0zh(2),b1zh(2),b21zh(2)
      dimension b0zw(2),b1zw(2),b21zw(2)
      dimension b0pw(2),b1pw(2),b21pw(2)
      dimension b0wh(2),b1wh(2),b21wh(2)
      dimension pggf(2),pggb(2),p3qf(2),s3qb(2),s3q(2),pgg(2),r1(2),
     #          r2(2),bcl(2),bch(2),s33f(2),s33(2),sr33(2),swwf(2),
     #          p33(2),pwwb(2),srwwb(2),sww(2),r3(2)
*
      co(1)= 1.d0
      co(2)= 0.d0
*
      call zztosbff(p2r,p2i,dmy2,dmy2,b0l,b1l,b21l)
      call zztosbff(p2r,p2i,tqm2,tqm2,b0h,b1h,b21h)
      call zztosbff(p2r,p2i,tqm2,dmy2,b0c,b1c,b21c)
*
      do i=1,2
          bfl(i)= 2.d0*b21l(i)-b0l(i)
          bfh(i)= 2.d0*b21h(i)-b0h(i)
      enddo
      do i=1,2
        pggf(i)= 80.d0/3.d0*bfl(i)+16.d0/3.d0*bfh(i)
       enddo  
*
      do i=1,2
       p3qf(i)= 10.d0*bfl(i)+2.d0*bfh(i)
      enddo
*
      call zztosbff(p2r,p2i,wm2,wm2,b0w,b1w,b21w)
      call zztosbff(p2r,p2i,zm2,zm2,b0z,b1z,b21z)
      call zztosbff(p2r,p2i,zm2,hm2,b0zh,b1zh,b21zh)
      call zztosbff(p2r,p2i,zm2,wm2,b0zw,b1zw,b21zw)
      call zztosbff(p2r,p2i,wm2,hm2,b0wh,b1wh,b21wh)
      call zztosbff(p2r,p2i,dmy2,wm2,b0pw,b1pw,b21pw)
*
      s3qb(1)= p2r*(2.d0/3.d0-10.d0*b21w(1)+13.d0/2.d0*b0w(1))-
     #         p2i*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))-
     #         2.d0*wm2*b0w(1)+2.d0*wm2*(delta-log(wm2))
      s3qb(2)= p2i*(2.d0/3.d0-10.d0*b21w(1)+13.d0/2.d0*b0w(1))+
     #         p2r*(-10.d0*b21w(2)+13.d0/2.d0*b0w(2))-
     #         2.d0*wm2*b0w(2)
*
      pggb(1)= 2.d0/3.d0-12.d0*b21w(1)+7.d0*b0w(1)
      pggb(2)= -12.d0*b21w(2)+7.d0*b0w(2)
*
      s3q(1)= p2r*p3qf(1)-p2i*p3qf(2)+s3qb(1)
      s3q(2)= p2r*p3qf(2)+p2i*p3qf(1)+s3qb(2)
*
      do i=1,2
       pgg(i)= pggf(i)+pggb(i)
      enddo
*
      r1(1)= s3q(1)-5.d0/46.d0*(p2r*pgg(1)-p2i*pgg(2))
      r1(2)= s3q(2)-5.d0/46.d0*(p2r*pgg(2)+p2i*pgg(1))
*
      do i=1,2
       p33(i)= 2.d0/3.d0*co(i)-9.d0*b21w(i)+25.d0/4.d0*b0w(i)-
     #         b21z(i)-b1z(i)-1.d0/4.d0*b0z(i)
       sr33(i)= -2.d0*wm2*b0w(i)+1.d0/2.d0*(zm2-hm2)*b1zh(i)+
     #          1.d0/4.d0*(5.d0*zm2-hm2)*b0zh(i)
      enddo
*
      s33f(1)= p2r*p3qf(1)-p2i*p3qf(2)-1/2*(p2r*(bfh(1)-bfl(1))-
     #         p2i*(bfh(2)-bfl(2)))-3.d0/2.d0*tqm2*b0h(1)
      s33f(2)= p2r*p3qf(2)+p2i*p3qf(1)-1/2*(p2r*(bfh(2)-bfl(2))+
     #         p2i*(bfh(1)-bfl(1)))-3.d0/2.d0*tqm2*b0h(2)
      s33(1)= p2r*p33(1)-p2i*p33(2)+sr33(1)+s33f(1)+4.d0*wm2*
     #        (delta-log(wm2))
      s33(2)= p2r*p33(2)+p2i*p33(1)+sr33(2)+s33f(2)
*
      r2(1)= s33(1)-5.d0/46.d0*(p2r*pgg(1)-p2i*pgg(2))
      r2(2)= s33(2)-5.d0/46.d0*(p2r*pgg(2)+p2i*pgg(1))
*
      do i=1,2
       bcl(i)= b21l(i)+b0l(i)
       bch(i)= b21c(i)+b0c(i)
      enddo
*
      swwf(1)= p2r*(12.d0*bfl(1)+6.d0*(bch(1)-bcl(1)))-
     #         p2i*(12.d0*bfl(2)+6.d0*(bch(2)-bcl(2)))-
     #         3.d0*tqm2*(b1c(1)+b0c(1))
      swwf(2)= p2r*(12.d0*bfl(2)+6.d0*(bch(2)-bcl(2)))+
     #         p2i*(12.d0*bfl(1)+6.d0*(bch(1)-bcl(1)))-
     #         3.d0*tqm2*(b1c(2)+b0c(2))
*
      do i=1,2
       pwwb(i)= 2.d0/3.d0*co(i)-9.d0*(b21zw(i)+b1zw(i))+7.d0/4.d0*
     #          b0zw(i)-b21wh(i)-b1wh(i)-1.d0/4.d0*b0wh(i)+sth2*(
     #          8.d0*b21zw(i)-2.d0*b0zw(i)+8.d0*b1zw(i)-
     #          8.d0*b21pw(i)-8.d0*b1pw(i)+2.d0*b0pw(i))
       srwwb(i)= 9.d0/2.d0*(zm2-wm2)*b1zw(i)+1.d0/4.d0*(13.d0*zm2-
     #           21.d0*wm2)*b0zw(i)+1.d0/2.d0*(wm2-hm2)*b1wh(i)+
     #           1.d0/4.d0*(5.d0*wm2-hm2)*b0wh(i)+sth2*(
     #           2.d0*(wm2-zm2)*(2.d0*b1zw(i)+b0zw(i))-2.d0*wm2*(
     #           2.d0*b1pw(i)+b0pw(i)))
      enddo
      sww(1)= p2r*pwwb(1)-p2i*pwwb(2)+srwwb(1)+swwf(1)+4.d0*wm2*
     #        (delta-log(wm2))
      sww(2)= p2r*pwwb(2)+p2i*pwwb(1)+srwwb(2)+swwf(2)
*
      r3(1)= sww(1)-5.d0/46.d0*(p2r*pgg(1)-p2i*pgg(2))
      r3(2)= sww(2)-5.d0/46.d0*(p2r*pgg(2)+p2i*pgg(1))
*
      return
      end
*
*----------------------------------------------------------
*
      subroutine zztohadr5(e,der,eder)
******************************************************************
*                                                                *
*      subroutine for the evaluation of the light hadron         *
*           contributions to delta_r  and  delta_g               *
*                    using fits to the                           *
*          qed vacuum polarization from e^+ e^- data             *
*                                                                *
*    f. jegerlehner, paul scherrer institute, ch-5232 villigen   *
*                                                                *
*    reference: f. jegerlehner, z. phys. c32 (1986) 195          *
*               h. burkhardt et al., z. phys. c42 (1989) 497     *
*               s. eidelman, f. jegerlehner, z. phys. c (1995)   *
*                                                                *
******************************************************************
*       version: 24/02/1995
*
*  notation: e energy ( momentum transfer ): e>0 timelike , e<0 spacelike
*            st2 is sin^2(theta); st2= 0.2322 is the reference value
*  the routine returns the hadronic contribution of 5 flavors (u,d,s,c,b)
*                 to   der= delta_r with hadronic error errder
*                and   deg= delta_g with hadronic error errdeg
*  the effective value of the fine structure constant alphaqed at energy
*  e is alphaqed(e)= alphaqed(0)/(1-delta_r) ,similarly for the su(2)
*  coupling alphasu2(e)= alphasu2(0)/(1-delta_g), where delta_r(g) is the
*  sum of leptonic, hadronic contributions (top to be added).
*
*  this program does not yet know how to compute delta r and delta g for
*  energies in the ranges  |e|>1tev and 2m_pi < e < 40(13) gev !!!!!!!!!
*
      implicit real*8(a-h,o-z)
      parameter(nf=9,ns=4)
*
      dimension res(ns)
      dimension ae(nf,ns),eu(nf),eo(nf),rm1(nf)
      dimension c1(nf,ns),c2(nf,ns),c3(nf,ns),c4(nf,ns),rl1(nf,ns)
*
       do i=1,nf
           eu(i)= 0.d0
           eo(i)= 0.d0
           rm1(i)= 0.d0
           do j=1,ns
               ae(i,j)= 0.d0
               c1(i,j)= 0.d0
               c2(i,j)= 0.d0
               c3(i,j)= 0.d0
               c4(i,j)= 0.d0
           enddo
       enddo
*
* #1# delta_r
* fit parameters spacelike  -1000 to  -200 gev
      eu(1)  = -1.d3
      eo(1)  = -2.d2
      rm1(1) = -1.d3
      c1(1,1)=  4.2069394d-02
      c2(1,1)=  2.9253566d-03
      c3(1,1)= -6.7782454d-04
      c4(1,1)=  9.3214130d-06
* fit parameters spacelike  -200 to  -20 gev
      eu(2)  = -2.d2
      eo(2)  =  -2.d1
      rm1(2) =  -1.d2
      c1(2,1)=  2.8526291d-02
      c2(2,1)=  2.9520725d-03
      c3(2,1)= -2.7906310d-03
      c4(2,1)=  6.4174528d-05
* fit parameters spacelike   -20 to   -2 gev
      eu(3)  = -2.d1
      eo(3)  = -2.d0
      rm1(3) = -2.d1
      rl1(3,1)= 9.3055d-03
      c1(3,1)=  2.8668314d-03
      c2(3,1)=  0.3514608d0
      c3(3,1)=  0.5496359d0
      c4(3,1)=  1.9892334d-04
      ae(3,1)=  3.d0
* fit parameters spacelike    -2 to    0.25 gev
      eu(4)  = -2.d0
      eo(4)  =  0.25d0
      rm1(4) = -2.d0
      rl1(4,1)= 9.3055d-03
      c1(4,1)=  2.2694240d-03
      c2(4,1)=  8.073429d0
      c3(4,1)=  0.1636393d0
      c4(4,1)= -3.3545541d-05
      ae(4,1)=  2.d0
* fit parameters timelike   0.25 to    2 gev
      eu(5)  =  0.25d0
      eo(5)  =  2.d0
* fit parameters timelike   2    to   40 gev
      eu(6)  =  2.d0
      eo(6)  =  4.d1
* fit parameters timelike     40 to   80 gev
      eu(7)  =  4.d1
      eo(7)  =  8.d1
      rm1(7) =  8.d1
      c1(7,1)=  2.7266588d-02
      c2(7,1)=  2.9285045d-03
      c3(7,1)= -4.7720564d-03
      c4(7,1)=  7.7295507d-04
* fit parameters timelike     80 to  250 gev
      eu(8)  =  8.d1
      eo(8)  =  2.5d2
      rm1(8) =  91.18880d0
      c1(8,1)=  2.8039809d-02
      c2(8,1)=  2.9373798d-03
      c3(8,1)= -2.8432352d-03
      c4(8,1)= -5.2537734d-04
* fit parameters timelike    250 to 1000 gev
      eu(9)  =  2.5d2
      eo(9)  =  1.d3
      rm1(9) =  1.d3
      c1(9,1)=  4.2092260d-02
      c2(9,1)=  2.9233438d-03
      c3(9,1)= -3.2966913d-04
      c4(9,1)=  3.4324117d-07
* #2# delta_g
* fit parameters spacelike  -1000 to  -200 gev
      c1(1,2)=  8.6415343d-02
      c2(1,2)=  6.0127582d-03
      c3(1,2)= -6.7379221d-04
      c4(1,2)=  9.0877611d-06
* fit parameters spacelike  -200 to  -20 gev
      c1(2,2)=  5.8580618d-02
      c2(2,2)=  6.0678599d-03
      c3(2,2)= -2.4153464d-03
      c4(2,2)=  6.1934326d-05
* fit parameters spacelike   -20 to   -2 gev
      rl1(3,2)= 1.9954d-02
      c1(3,2)=  5.7231588d-03
      c2(3,2)=  0.3588257d0
      c3(3,2)=  0.5532265d0
      c4(3,2)=  6.0730567d-04
      ae(3,2)=  3.d0
* fit parameters spacelike    -2 to    0.25 gev
      rl1(4,2)= 1.9954d-02
      c1(4,2)=  4.8065037d-03
      c2(4,2)=  8.255167d0
      c3(4,2)=  0.1599882d0
      c4(4,2)= -1.8624817d-04
      ae(3,2)=  2.d0
* fit parameters timelike     40 to   80 gev
      c1(7,2)=  5.5985276d-02
      c2(7,2)=  6.0203830d-03
      c3(7,2)= -5.0066952d-03
      c4(7,2)=  7.1363564d-04
* fit parameter timelike     80 to  250 gev
      c1(8,2)=  5.7575710d-02
      c2(8,2)=  6.0372148d-03
      c3(8,2)= -3.4556778d-03
      c4(8,2)= -4.9574347d-04
* fit parameters timelike    250 to 1000 gev
      c1(9,2)=  8.6462371d-02
      c2(9,2)=  6.0088057d-03
      c3(9,2)= -3.3235471d-04
      c4(9,2)=  5.9021050d-07
* #3# error delta_r
* fit parameters spacelike  -1000 to  -200 gev
      c1(1,3)=  6.3289929d-04
      c2(1,3)=  3.3592437d-06
      c3(1,3)=  0.d0
      c4(1,3)=  0.d0
* fit parameters spacelike  -200 to  -20 gev
      c1(2,3)=  6.2759849d-04
      c2(2,3)= -1.0816625d-06
      c3(2,3)=  5.050189d0
      c4(2,3)= -9.6505374d-02
      ae(2,3)=  1.d0
* fit parameters spacelike   -20 to   -2 gev
      rl1(3,3)= 2.0243d-04
      c1(3,3)=  1.0147886d-04
      c2(3,3)=  1.819327d0
      c3(3,3)= -0.1174904d0
      c4(3,3)= -1.2404939d-04
      ae(3,3)=  3.d0
* fit parameters spacelike    -2 to    0.25 gev
      rl1(4,3)= 2.0243d-04
      c1(4,3)= -7.1368617d-05
      c2(4,3)=  9.980347d-04
      c3(4,3)=  1.669151d0
      c4(4,3)=  3.5645600d-05
      ae(4,3)=  2.d0
* fit parameters timelike     40 to   80 gev
      c1(7,3)=  6.4947648d-04
      c2(7,3)=  4.9386853d-07
      c3(7,3)= -55.22332d0
      c4(7,3)=  26.13011d0
* fit parameters timelike     80 to  250 gev
      c1(8,3)=  6.4265809d-04
      c2(8,3)= -2.8453374d-07
      c3(8,3)= -23.38172d0
      c4(8,3)= -6.251794d0
* fit parameter timelike    250 to 1000 gev
      c1(9,3)=  6.3369947d-04
      c2(9,3)= -2.0898329d-07
      c3(9,3)=  0.d0
      c4(9,3)=  0.d0
* #4# error delta_g
* fit parameters spacelike  -1000 to  -200 gev
      c1(1,4)=  1.2999176d-03
      c2(1,4)=  7.4505529d-06
      c3(1,4)=  0.d0
      c4(1,4)=  0.d0
* fit parameters spacelike  -200 to  -20 gev
      c1(2,4)=  1.2883141d-03
      c2(2,4)= -1.3790827d-06
      c3(2,4)=  8.056159d0
      c4(2,4)= -0.1536313d0
      ae(2,4)=  1.d0
* fit parameters spacelike   -20 to   -2 gev
      rl1(3,4)= 4.3408d-04
      c1(3,4)=  2.0489733d-04
      c2(3,4)=  2.065011d0
      c3(3,4)= -0.6172962d0
      c4(3,4)= -2.5603661d-04
      ae(3,4)=  3.d0
* fit parameters spacelike    -2 to    0.25 gev
      rl1(4,4)= 4.3408d-04
      c1(4,4)= -1.5095409d-04
      c2(4,4)=  9.9847501d-04
      c3(4,4)=  1.636659d0
      c4(4,4)=  7.5892596d-05
      ae(4,4)=  2.d0
* fit parameters timelike     40 to   80 gev
      c1(7,4)=  1.3335156d-03
      c2(7,4)=  2.2939612d-07
      c3(7,4)= -246.4966d0
      c4(7,4)=  114.9956d0
* fit parameters timelike     80 to  250 gev
      c1(8,4)=  1.3196438d-03
      c2(8,4)=  2.8937683d-09
      c3(8,4)=  5449.778d0
      c4(8,4)=  930.3875d0
* fit parameters timelike    250 to 1000 gev
      c1(9,4)=  1.3016918d-03
      c2(9,4)= -3.6027674d-07
      c3(9,4)=  0.d0
      c4(9,4)=  0.d0
*
       se=  654.d0/643.d0      ! rescaling error to published version 1995
       s= e*e
       der= 0.d0
       deg= 0.d0
       eder= 0.d0
       edeg= 0.d0
       if((e.gt.1.e3).or.(e.lt.-1.e3)) goto 100
       if((e.lt.4.e1).and.(e.gt.0.25d0)) goto 100
       i= 1
       do while (e.ge.eo(i))
         i= i+1
       enddo
       if(e.eq.1.e3) i= 9
       if(e.eq.0.d0 ) goto 100
       s0= sign(1.d0,rm1(i))*rm1(i)**2
       s = sign(1.d0,e)*e*e
       x1= s0/s
       xi= 1.d0/x1
       x2= x1*x1
       if(ae(i,1).le.0.d0) then
          do j= 1,4
              xlar= xi+ae(i,j)*exp(-xi)
              xlog= log(xlar)
              res(j)= c1(i,j)+c2(i,j)*(xlog+c3(i,j)*(x1-1.d0)+
     #                c4(i,j)*(x2-1.0))
          enddo
       else if (ae(i,1).eq.2.d0) then
          hx= xi*xi
          do j= 1,2
              fx= 1.d0-c2(i,j)*s
              gx= c3(i,j)*s/(c3(i,j)-s)
              xx= log(abs(fx))+c2(i,j)*gx
              res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx
          enddo
          do j= 3,4
              u= abs(s)
              gx= -c3(i,j)*u/(c3(i,j)+u)
              xx= xi**3/(sqrt(abs(xi))**5+c2(i,j))
              res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx
          enddo
       else if (ae(i,1).eq.3.0) then
          hx= xi
          do j= 1,4
              fx= 1.d0-c2(i,j)*s
              gx= c3(i,j)*s/(c3(i,j)-s)
              xx= log(abs(fx))+c2(i,j)*gx
              res(j)= c1(i,j)*xx-rl1(i,j)*gx+c4(i,j)*hx
          enddo
       endif
       der= res(1)
       eder= res(3)*se
100    return
       end
*
*-----------------------------------------------------
*
      real*8 function zztoli2(xr,xi)
      implicit real*8(a,b,d-h,o-z)
      implicit complex*16 (c)
*
      common/zzted/eps,delta
      common/zztparam/pi,pis,pih,gf,g8,sth2,cth8,alphai
*
      data  b0/1.d0/,b1/-0.25d0/,b2/0.277777777777778d-1/,
     #      b4/-0.277777777777778d-3/,b6/0.472411186696901d-5/
*
      c0= cmplx(0.d0)
      pis= pi*pi
      pis6= pis/6.d0
      cx= cmplx(xr,xi)
      comx= 1.d0-cx
      if (xr.lt.0.d0) then
          cy= 1.d0-cx
          sign1= -1.d0
          clnx= log(cx)
          clnomx= log(comx)
          cadd1= pis6-clnx*clnomx
      else
          cy= cx
          sign1= 1.d0
          cadd1= c0
      endif
      comy= 1.d0-cy
      ym2= conjg(cy)*cy
      ym= sqrt(ym2)
      if (ym.gt.1.d0) then
          cz= conjg(cy)/ym2
          sign2= -1.d0
          coy= -cy
          clnoy= log(coy)
          cadd2= -pis6-0.5d0*clnoy*clnoy
      else
          cz= cy
          sign2= 1.d0
          cadd2= c0
      endif
      comz= 1.d0-cz
      zr= real(cz)
      if (zr.gt.0.5d0) then
          ct= 1.d0-cz
          comt= 1.d0-ct
          sign3= -1.d0
          clnz= log(cz)
          clnomz= log(comz)
          cadd3= pis6-clnz*clnomz
      else
          ct= cz
          comt= 1.d0-ct
          sign3= 1.d0
          cadd3= c0
      endif
      cpar= log(comt)
      cpar2= cpar*cpar
*      
      cres= -cpar*(b0-cpar*(b1-cpar*(b2+cpar2*(b4+b6*cpar2))))
*
      cli2= sign1*(sign2*(sign3*cres+cadd3)+cadd2)+cadd1
*
      zztoli2= real(cli2)
*
      return
      end
*
      subroutine zztocflag(oflag,oval)
      implicit real*8(a-h,o-z)
*
      character*(*) oflag,oval
      character*4 ochan
      character*1 out,ofl,otop,oqed,onu,oqedf,oqcdf,oscheme
*
      common/zztn/onu
      common/zztout/out
      common/zztch/ochan
      common/zztoqed/oqed
      common/zztqedfsr/oqedf
      common/zztqedfsr/oqcdf
      common/zztofl/ofl,otop
      common/zztosch/oscheme
*
      if(oflag.eq.'onu') then
         onu= oval
      else if(oflag.eq.'out') then
         out= oval
      else if(oflag.eq.'ochan') then
         ochan= oval
      else if(oflag.eq.'oqed') then
         oqed= oval
      else if(oflag.eq.'oqedf') then
         oqedf= oval
      else if(oflag.eq.'oqcdf') then
         oqcdf= oval
      else if(oflag.eq.'ofl') then
         ofl= oval
      else if(oflag.eq.'otop') then
         otop= oval
      else if(oflag.eq.'oscheme') then
         oscheme= oval
      endif
*
      return
      end