c
      real function riembi(p,q,x)
c implementation of Lentz's continued fraction 
c based on ideas from Numerical Recipes (Press et. al., 1986)
      implicit none
      real p,q,x,gamlog,rbfrac,zero,one,two,sav
      parameter (zero=0.,one=1.,two=2.)
c Abramovitz & Stegun 1971 (26.5.8)
      if (x.lt.zero.or.x.gt.one) then
        pause 'riembi: argument X out of range'
      elseif (x.eq.zero.or.x.eq.one) then
        sav=zero
      else
c Abramovitz & Stegun 1971 (6.2.2)
        sav=exp(gamlog(p+q)
     $-gamlog(p)-gamlog(q)+p*log(x)+q*log(one-x))
      endif
c Abramovitz & Stegun 1971 (26.5.8)
      if ((p+one)/(p+q+two).gt.x) then
        riembi=sav*rbfrac(p,q,x)/p
      else
c Abramovitz & Stegun 1971 (26.5.2)
        riembi=one-sav*rbfrac(q,p,one-x)/q
      endif
      end
c
c                                   
c