double complex function cdig(alpha,x)
c --- Written By Eric Kostlan & Dmitry Gokhman
c --- March  1986
      double complex alpha,x,cdh
      double complex re,one,p,q
      double precision xlim,zero,dnrm
      data re,one/(0.36787944117144232d0,0.),1./
      data xlim,zero/1.,0./ibuf/34/
c --- If x is near the negative real axis, then shift to x=1.
      if(dnrm(x).lt.xlim.or.dreal(x).lt.zero.and.
     *   dabs(dimag(x)).lt.xlim)then
         cdig=re/cdh(alpha,one)
         ilim=dreal(x/re)
         do 1 i=0,ibuf-ilim
            call term(alpha,x,i,p,q)
            cdig=cdig+p*q
    1    continue
      else
         cdig=cdexp(-x+alpha*cdlog(x))/cdh(alpha,x)
      endif
      return
      end
c
      double complex function cdh(alpha,x)
c --- Written By Eric Kostlan & Dmitry Gokhman
c --- March  1986
      double complex alpha,x,cdhs
      double complex one,term,sum,cn,alpha1
      data one/1./
c --- If Re(alpha-x) is too big, shift alpha.
      n=dreal(alpha-x)
      if(n.gt.0)then
         cn=n
         alpha1=alpha-cn
         term=one/x
         sum=term
         do 1 i=1,n-1
            cn=n-i
            term=term*(alpha1+cn)/x
            sum=term+sum
    1    continue
         sum=sum+term*alpha1/cdhs(alpha1,x)
         cdh=one/sum
      else
         cdh=cdhs(alpha,x)
      endif
      return
      end
c
      double complex function cdhs(alpha,x)
c --- Written By Eric Kostlan & Dmitry Gokhman
c --- March  1986
      double complex zero,half,one,alpha,x
      double complex p0,q0,p1,q1,r0,r1,ci,factor
      double precision tol1,tol2,error,dnrm
      data zero,half,one/0.,0.5,1./
      data tol1,tol2,error/1.d10,1.d-10,5.d-18/ilim/100000/
      q0=one
      q1=one
      p0=x
      p1=x+one-alpha
      do 1 i=1,ilim
         ci=i
         if(p0.ne.zero.and.q0.ne.zero.and.q1.ne.zero)then
            r0=p0/q0
            r1=p1/q1
            if(dnrm(r0-r1).le.dnrm(r1)*error)then
               cdhs=r1
               return
            endif
c --------- Occasionally renormalize the sequences to avoid over(under)flow.
            if(dnrm(p0).gt.tol1.or.dnrm(p0).lt.tol2.or.
     *         dnrm(q0).gt.tol1.or.dnrm(q0).lt.tol2)then
               factor=p0*q0
               p0=p0/factor
               q0=q0/factor
               p1=p1/factor
               q1=q1/factor
            endif
         endif
         p0=x*p1+ci*p0
         q0=x*q1+ci*q0
         p1=p0+(ci+one-alpha)*p1
         q1=q0+(ci+one-alpha)*q1
    1 continue
c --- If the peripheral routines are written correctly,
c --- the following four statements should never be executed.
      print*,'cdhs:  *** Warning: i >',ilim
      print*,'cdhs:  *** r0,r1= ',r0,r1
      cdhs=half*(r0+r1)
      return
      end
c
      subroutine term(alpha,x,i,p,q)
c --- Calculate p*q = -1**i(1-x**(alpha+i))/(alpha+i)i! carefully.
      double complex alpha,x,p,q,ci,alphai
      double complex zero,one,two,cdlx
      double precision tol,xlim,dnrm
      data zero,one,two/0.,1.,2./tol/3.d-7/xlim/39./
      if(i.eq.0)q=one
      ci=i
      alphai=alpha+ci
      if(x.eq.zero)then
         p=one/alphai
         if(i.ne.0)q=-q/ci
         return
      endif
      cdlx=cdlog(x)
c --- If (1-x**alphai)=-x**alphai on the computer,
c --- then change the inductive scheme to avoid overflow.
      if(dreal(alphai*cdlx).gt.xlim.and.i.ne.0)then
         p=p*(alphai-one)/alphai
         q=-q*x/ci
         return
      endif
      if(dnrm(alphai).gt.tol)then
         p=(one-x**alphai)/alphai
      else
         p=-cdlx*(one+cdlx*alphai/two)
      endif
      if(i.ne.0)q=-q/ci
      return
      end
c
      double precision function dnrm(z)
      double complex z
      dnrm=dabs(dreal(z))+dabs(dimag(z))
      return
      end