*###[ ffzli2:
	subroutine ffzli2(zdilog,zlog,cx,ier)
***#[*comment:***********************************************************
*									*
*	Computes the dilogarithm (Li2, Sp) for any (complex) cx 	*
*	to a precision precc.  It assumes that cx is already in the	*
*	area |cx|<=1, Re(cx)<=1/2.  As it is available it also returns	*
*	log(1-cx) = zlog.						*
*									*
*	Input:	cx	(complex)					*
*									*
*	Output:	zdilog	(complex) Li2(cx)				*
*		zlog	(complex) log(1-cx) = -Li1(cx)			*
*		ier	(integer) 0=OK,1=num,2=err			*
*									*
*	Calls:	log,zfflo1,(d/a)imag,real/dble				*
*									*
***#]*comment:*********************************************************** 
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	DOUBLE COMPLEX cx,zlog,zdilog
*
*	local variables
*
	DOUBLE PRECISION xprec,bdn02,bdn05,bdn10,bdn15,
     +		xi,xr,xdilog,xlog,absc,xa,a,ffbnd
	DOUBLE COMPLEX cc,cz,cz2,zfflo1
	external ffbnd,zfflo1
	save xprec,bdn02,bdn05,bdn10,bdn15
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*  #] declarations: 
*  #[ initialisations:
	data xprec /-1D0/
	if ( xprec .ne. precc ) then
	    xprec = precc
	    bdn02 = ffbnd(1,2,bf)
	    bdn05 = ffbnd(1,5,bf)
	    bdn10 = ffbnd(1,10,bf)
	    bdn15 = ffbnd(1,15,bf)
*		we don't have bf(21) ...
	endif
*  #] initialisations: 
*  #[ exceptional cases:
	xi  = DIMAG(cx)
	xr  = DBLE(cx)
	if ( xi .eq. 0) then
	    call ffxli2(xdilog,xlog,xr,ier)
	    zdilog = xdilog
	    zlog = xlog
	    return
	endif
	xa = abs(xi) + abs(xr)
	if ( xa .lt. precc ) then
	    zdilog = cx
	    zlog = -cx
	    return
	endif
*  #] exceptional cases: 
*  #[ get log,dilog:
	if ( xa .lt. xloss**2 ) then
	    zlog = zfflo1(cx,ier)
	else
	    zlog = log(1-cx)
	endif
	cz = -zlog
	if ( absc(cz) .lt. xclog2 ) then
	    zdilog = cz
	else
	cz2 = cz*cz
	a = xa**2
	if ( a .gt. bdn15 ) then
	    zdilog = cz2*(DBLE(bf(16)) + cz2*(DBLE(bf(17))
     +		   + cz2*(DBLE(bf(18)) + cz2*(DBLE(bf(19))
     +		   + cz2*(DBLE(bf(20)))))))
	else
	    zdilog = 0
	endif
	if ( a .gt. bdn10 ) then
	    zdilog = cz2*(DBLE(bf(11)) + cz2*(DBLE(bf(12))
     +		   + cz2*(DBLE(bf(13)) + cz2*(DBLE(bf(14))
     +		   + cz2*(DBLE(bf(15)) + zdilog)))))
	endif
	if ( a .gt. bdn05 ) then
	    zdilog = cz2*(DBLE(bf(6)) + cz2*(DBLE(bf(7))
     +		   + cz2*(DBLE(bf(8)) + cz2*(DBLE(bf(9))
     +		   + cz2*(DBLE(bf(10)) + zdilog)))))
	endif
	if ( a .gt. bdn02 ) then
	    zdilog = cz2*(DBLE(bf(3)) + cz2*(DBLE(bf(4))
     +		   + cz2*(DBLE(bf(5)) + zdilog)))
	endif
*	watch the powers of z.
	zdilog = cz + cz2*(DBLE(bf(1)) + cz*(DBLE(bf(2)) + zdilog))
	endif
*  #] get log,dilog: 
*###] ffzli2:
	end


*###[ ffzzdl:
	subroutine ffzzdl(zdilog,ipi12,zlog,cx,ier)
***#[*comment:***************************************************
*								*
*	Computes the dilogarithm (Li2, Sp) for any (complex) cx *
*	to about 15 significant figures. This can be improved	*
*	by adding more of the bf's. For real cx > 1 an error is	*
*	generated as the imaginary part is undefined then.	*
*	For use in ffcdbd zlog = log(1-cx) is also calculated	*
*								*
*	Input:	cx	(complex)				*
*								*
*	Output:	zdilog	(complex) Li2(cx) mod factors pi^2/12	*
*		ipi12	(integer) these factors			*
*		zlog	(complex) log(1-cx)			*
*								*
*	Calls:	log,zfflo1,(d/a)imag,real/dble			*
*								*
***#]*comment:*************************************************** 
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ipi12,ier
	DOUBLE COMPLEX zdilog,zlog,cx
*
*	local variables
*
	integer jsgn
	DOUBLE PRECISION xprec,bdn02,bdn05,bdn10,bdn15,
     +		xi,xr,s1,s2,xa,a,absc,ffbnd
	DOUBLE COMPLEX cfact,cx1,cy,cz,cz2,zfflo1,c
	external ffbnd,zfflo1
	save xprec,bdn02,bdn05,bdn10,bdn15
*
*	common blocks
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations: 
*  #[ initialisations:
	data xprec /-1D0/
	if ( xprec .ne. precc ) then
	    xprec = precc
	    bdn02 = ffbnd(1,2,bf)
	    bdn05 = ffbnd(1,5,bf)
	    bdn10 = ffbnd(1,10,bf)
	    bdn15 = ffbnd(1,15,bf)
	endif
*  #] initialisations: 
*  #[ exceptional cases:
	xi  = DIMAG(cx)
	xr  = DBLE(cx)
	if ( xi .eq. 0 ) then
	    if ( xr .gt. 1 ) call fferr(31,ier)
	    call ffzxdl(zdilog,ipi12,zlog,xr,1,ier)
	    return
	endif
	if ( abs(xi) .lt. xalog2 ) then
	    s1 = 0
	else
	    s1 = xi**2
	endif
	if ( abs(xr) .lt. xalog2 ) then
	    s2 = 0
	else
	    s2 = xr**2
	endif
	xa = sqrt(s1 + s2)
	if ( xa .lt. precc ) then
	    zdilog = cx
	    zlog = -cx
	    ipi12 = 0
	    return
	endif
*  #] exceptional cases: 
*  #[ transform to |x|<1, Re(x) < 0.5:
	if ( xr .le. .5D0) then
	    if (xa .gt. 1) then
		if ( 1/xa .lt. xalogm ) then
		    cfact = 0
		elseif ( 1/xa .lt. xclogm ) then
		    cx1 = cx*DBLE(1/xa)
		    cfact = log(-cx1) + log(DBLE(xa))
		else
		    cfact = log(-cx)
		endif
		cy = - cfact**2/2
		ipi12 = -2
		if ( xa*xloss**2 .gt. 1) then
		    if ( 1/xa .lt. xclogm ) then
			cx1 = cx*DBLE(1/xa)
			cx1 = 1/cx1
			cx1 = cx1*DBLE(1/xa)
		    else
			cx1 = 1/cx
		    endif
		    cz = -zfflo1(cx1,ier)
		else
		    cz = -log(1-1/cx)
		endif
		zlog = log(1-cx)
		jsgn = -1
	    else
		cy = 0
		ipi12 = 0
		if ( xa .lt. xloss**2 ) then
		    zlog = zfflo1(cx,ier)
		else
		    zlog = log(1-cx)
		endif
		cz = -zlog
		jsgn = 1
	    endif
	else
	    if (xa .le. sqrt(2*xr)) then
		cz = -log(cx)
		if ( abs(xr-1) + abs(xi) .lt. xclogm ) then
		    cy = 0
		else
		    zlog = log(1-cx)
		    cy = cz*zlog
		endif
		ipi12 = 2
		jsgn = -1
	    else
		if ( 1/xa .lt. xalogm ) then
		    cfact = 0
		elseif ( 1/xa .lt. xclogm ) then
		    cx1 = cx*DBLE(1/xa)
		    cfact = log(-cx1) + log(DBLE(xa))
		else
		    cfact = log(-cx)
		endif
		cy = - cfact**2/2
		ipi12 = -2
		if ( xa*xloss .gt. 1) then
		    if ( 1/xa .lt. xclogm ) then
			cx1 = cx*DBLE(1/xa)
			cx1 = 1/cx1
			cx1 = cx1*DBLE(1/xa)
		    else
			cx1 = 1/cx
		    endif
		    cz = -zfflo1(cx1,ier)
		else
		    cz = -log(1-1/cx)
		endif
		zlog = log(1-cx)
		jsgn = -1
	    endif
	endif
*  #] transform to |x|<1, Re(x) < 0.5: 
*  #[ get dilog:
	if ( absc(cz) .lt. xclogm ) then
	    zdilog = cz
	else
	cz2 = cz*cz
	a = DBLE(cz)**2 + DIMAG(cz)**2
	if ( a .gt. bdn15 ) then
	    zdilog = cz2*(DBLE(bf(16)) + cz2*(DBLE(bf(17))
     +		   + cz2*(DBLE(bf(18)) + cz2*(DBLE(bf(19))
     +		   + cz2*(DBLE(bf(20)))))))
	else
	    zdilog = 0
	endif
	if ( a .gt. bdn10 ) then
	    zdilog = cz2*(DBLE(bf(11)) + cz2*(DBLE(bf(12))
     +		   + cz2*(DBLE(bf(13)) + cz2*(DBLE(bf(14))
     +		   + cz2*(DBLE(bf(15)) + zdilog)))))
	endif
	if ( a .gt. bdn05 ) then
	    zdilog = cz2*(DBLE(bf(6)) + cz2*(DBLE(bf(7))
     +		   + cz2*(DBLE(bf(8)) + cz2*(DBLE(bf(9))
     +		   + cz2*(DBLE(bf(10)) + zdilog)))))
	endif
	if ( a .gt. bdn02 ) then
	    zdilog = cz2*(DBLE(bf(3)) + cz2*(DBLE(bf(4))
     +		   + cz2*(DBLE(bf(5)) + zdilog)))
	endif
*	watch the powers of z.
	zdilog = cz + cz2*(DBLE(bf(1)) + cz*(DBLE(bf(2)) + zdilog))
	endif
	if(jsgn.eq.1)then
	    zdilog =  zdilog + cy
	else
	    zdilog = -zdilog + cy
	endif
*  #] get dilog: 
*###] ffzzdl: 
	end


*###[ zfflog:
	DOUBLE COMPLEX function zfflog(cx,ieps,cy,ier)
***#[*comment:***********************************************************
*									*
*	Calculate the complex logarithm of cx.  The following cases	*
*	are treted separately:						*
*		|cx| too small:		give warning and return 0	*
*					(for Absoft, Apollo DN300)	*
*		Im(cx) = 0, Re(cx) < 0:	take sign according to ieps	*
*									*
***#]*comment:*********************************************************** 
*  #[ declarations:
*
*	arguments
*
	implicit none
	integer ieps,ier
	DOUBLE COMPLEX cx,cy
*
*	local variables
*
	DOUBLE COMPLEX c,ctroep
	DOUBLE PRECISION absc,xa,xlog1p
*
*	common blocks, statement function
*
	include 'ff.h'
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations: 
*  #[ calculations:
	xa = absc(cx)
	if ( xa .lt. xalogm ) then
	    if ( cx .ne. 0 ) call fferr(23,ier)
	    zfflog = 0
	elseif ( DBLE(cx) .lt. 0 .and. DIMAG(cx) .eq. 0 ) then
*     +		 abs(DIMAG(cx)) .lt. precc*abs(DBLE(cx)) ) then
	    xlog1p = log(-DBLE(cx))
*	    checked imaginary parts 19-May-1988
	    if ( abs(ieps) .eq. 1 ) then
		if ( ieps*DBLE(cy) .lt. 0 ) then
		    zfflog = DCMPLX(xlog1p,-pi)
		elseif ( ieps*DBLE(cy) .gt. 0 ) then
		    zfflog = DCMPLX(xlog1p,pi)
		else
		    call fferr(51,ier)
		    zfflog = DCMPLX(xlog1p,pi)
		endif
	    elseif ( ieps .ge. 2 .and. ieps .le. 3 ) then
		zfflog = DCMPLX(xlog1p,-pi)
	    elseif ( ieps .le. -2 .and. ieps .ge. -3 ) then
		zfflog = DCMPLX(xlog1p,pi)
	    else
		call fferr(51,ier)
		zfflog = DCMPLX(xlog1p,pi)
	    endif
	elseif ( xa .lt. xclogm .or. 1/xa .lt. xclogm ) then
	    ctroep = cx*DBLE(1/xa)
	    zfflog = log(ctroep) + DBLE(log(xa))
	else
*	    print *,'zfflog: neem log van ',cx
	    zfflog = log(cx)
	endif
*  #] calculations: 
*###] zfflog: 
	end


*###[ zfflo1:
	DOUBLE COMPLEX function zfflo1(cx,ier)
***#[*comment:***************************************************
*	calculates log(1-x) for |x|<.14 in a faster way to ~15	*
*	significant figures.					*
***#]*comment:*************************************************** 
*  #[ declarations:
	implicit none
	integer ier
	DOUBLE COMPLEX cx,c,zfflog
	DOUBLE PRECISION xprec,bdn01,bdn05,bdn10,bdn15,bdn19,
     +		absc,xa,ffbnd
	external zfflog,ffbnd
	save xprec,bdn01,bdn05,bdn10,bdn15,bdn19
	include 'ff.h'
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*  #] declarations: 
*  #[ initialisations:
	data xprec /-1D0/
	if ( precc .ne. xprec ) then
	    xprec = precc
*	    determine the boundaries for 1,5,10,15 terms
	    bdn01 = ffbnd(1,1,xninv)
	    bdn05 = ffbnd(1,5,xninv)
	    bdn10 = ffbnd(1,10,xninv)
	    bdn15 = ffbnd(1,15,xninv)
	    bdn19 = ffbnd(1,19,xninv)
	endif
*  #] initialisations: 
*  #[ calculations:
	xa = absc(cx)
	if ( xa .gt. bdn19 ) then
	    c = cx-1
	    xa = absc(c)
	    zfflo1 = zfflog(1-cx,0,czero,ier)
	    return
	endif
	if ( xa .gt. bdn15 ) then
	    zfflo1 = cx*( DBLE(xninv(16)) + cx*( DBLE(xninv(17))
     +		   + cx*( DBLE(xninv(18)) + cx*( DBLE(xninv(19))
     +		   + cx*( DBLE(xninv(20)) )))))
	else
	    zfflo1 = 0
	endif
	if ( xa .gt. bdn10 ) then
	    zfflo1 = cx*( DBLE(xninv(11)) + cx*( DBLE(xninv(12))
     +		   + cx*( DBLE(xninv(13)) + cx*( DBLE(xninv(14))
     +		   + cx*( DBLE(xninv(15)) + zfflo1 )))))
	endif
	if ( xa .gt. bdn05 ) then
	    zfflo1 = cx*( DBLE(xninv(6)) + cx*( DBLE(xninv(7))
     +		   + cx*( DBLE(xninv(8)) + cx*( DBLE(xninv(9))
     +		   + cx*( DBLE(xninv(10)) + zfflo1 )))))
	endif
	if ( xa .gt. bdn01 ) then
	    zfflo1 = cx*( DBLE(xninv(2)) + cx*( DBLE(xninv(3))
     +		   + cx*( DBLE(xninv(4)) + cx*( DBLE(xninv(5))
     +		   + zfflo1 ))))
	endif
	zfflo1 = - cx*( DBLE(xninv(1)) + zfflo1 )
*  #] calculations: 
*###] zfflo1: 
	end


*###[ zfflo2:
	DOUBLE COMPLEX function zfflo2(x,ier)
***#[*comment:***************************************************
*	calculates log(1-x)+x for |x|<.14 in a faster way to	*
*	~15 significant figures.				*
***#]*comment:*************************************************** 
*  #[ declarations:
	implicit none
	integer ier
	DOUBLE COMPLEX x,zfflo1,cc
	DOUBLE PRECISION bdn01,bdn05,bdn10,bdn15,bdn18,xprec,xa,
     +		ffbnd,absc
	external ffbnd,zfflo1
	save xprec,bdn01,bdn05,bdn10,bdn15,bdn18
	include 'ff.h'
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*  #] declarations: 
*  #[ initialisation:
	data xprec /-1D0/
	if ( xprec .ne. precc ) then
	    xprec = precx
	    precx = precc
*	    determine the boundaries for 1,5,10,15 terms
	    bdn01 = ffbnd(1,1,xninv(2))
	    bdn05 = ffbnd(1,5,xninv(2))
	    bdn10 = ffbnd(1,10,xninv(2))
	    bdn15 = ffbnd(1,15,xninv(2))
	    bdn18 = ffbnd(1,18,xninv(2))
	    precx = xprec
	    xprec = precc
	endif
*  #] initialisation: 
*  #[ calculations:
	xa = absc(x)
	if ( xa .gt. bdn18 ) then
	    zfflo2 = zfflo1(x,ier) + x
	    return
	endif
	if ( xa .gt. bdn15 ) then
	    zfflo2 = x*( DBLE(xninv(17)) + x*( DBLE(xninv(18)) +
     +		x*( DBLE(xninv(19)) + x*( DBLE(xninv(20)) ))))
	else
	    zfflo2 = 0
	endif
	if ( xa .gt. bdn10 ) then
	    zfflo2 = x*( DBLE(xninv(12)) + x*( DBLE(xninv(13)) +
     +		x*( DBLE(xninv(14)) + x*( DBLE(xninv(15)) +
     +		x*( DBLE(xninv(16)) + zfflo2 )))))
	endif
	if ( xa .gt. bdn05 ) then
	    zfflo2 = x*( DBLE(xninv(7)) + x*( DBLE(xninv(8)) +
     +		x*( DBLE(xninv(9)) +x*( DBLE(xninv(10)) +
     +		x*( DBLE(xninv(11)) + zfflo2 )))))
	endif
	if ( xa .gt. bdn01 ) then
	    zfflo2 = x*( DBLE(xninv(3)) + x*( DBLE(xninv(4)) +
     +		x*( DBLE(xninv(5)) + x*( DBLE(xninv(6)) + zfflo2 ))))
	endif
	zfflo2 = - x**2*( DBLE(xninv(2)) + zfflo2 )
*  #] calculations: 
*###] zfflo2: 
	end


*###[ zfflo3:
	DOUBLE COMPLEX function zfflo3(x,ier)
***#[*comment:***************************************************
*	calculates log(1-x)+x+x^2/2 for |x|<.14 in a faster 	*
*	way to ~15 significant figures.				*
***#]*comment:*************************************************** 
*  #[ declarations:
	implicit none
	integer ier
	DOUBLE COMPLEX x,zfflo2,cc
	DOUBLE PRECISION bdn01,bdn05,bdn10,bdn15,xprec,xa,ffbnd,
     +		absc
	external zfflo2,ffbnd
	save xprec,bdn01,bdn05,bdn10,bdn15
	include 'ff.h'
	absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
*  #] declarations: 
*  #[ initialisation:
	data xprec /-1D0/
	if ( xprec .ne. precx ) then
	    xprec = precx
	    precx = precc
*	    determine the boundaries for 1,5,10,15 terms
	    bdn01 = ffbnd(1,1,xninv(3))
	    bdn05 = ffbnd(1,5,xninv(3))
	    bdn10 = ffbnd(1,10,xninv(3))
	    bdn15 = ffbnd(1,15,xninv(3))
	    precx = xprec
	    xprec = precc
	endif
*  #] initialisation: 
*  #[ calculations:
	xa = absc(x)
	if ( xa .gt. bdn15 ) then
	    zfflo3 = zfflo2(x,ier) + x**2/2
	    return
	endif
	if ( xa .gt. bdn10 ) then
	    zfflo3 = x*( DBLE(xninv(13)) + x*( DBLE(xninv(14)) +
     +		x*( DBLE(xninv(15)) + x*( DBLE(xninv(16)) +
     +		x*( DBLE(xninv(17)) )))))
	else
	    zfflo3 = 0
	endif
	if ( xa .gt. bdn05 ) then
	    zfflo3 = x*( DBLE(xninv(8)) + x*( DBLE(xninv(9)) +
     +		x*( DBLE(xninv(10)) + x*( DBLE(xninv(11)) +
     +		x*( DBLE(xninv(12)) + zfflo3 )))))
	endif
	if ( xa .gt. bdn01 ) then
	    zfflo3 = x*( DBLE(xninv(4)) + x*( DBLE(xninv(5)) +
     +		x*( DBLE(xninv(6)) + x*( DBLE(xninv(7)) + zfflo3 ))))
	endif
	zfflo3 = - x**3*( DBLE(xninv(3)) + zfflo3 )
*  #] calculations: 
*###] zfflo3: 
	end