double precision function dbetai (x, pin, qin)
c august 1980 edition.  w. fullerton, c3, los alamos scientific lab.
c based on bosten and battiste, remark on algorithm 179, comm. acm,
c v 17, p 153, (1974).
c
c             input arguments --
c x      upper limit of integration.  x must be in (0,1) inclusive.
c p      first beta distribution parameter.  p must be gt 0.0.
c q      second beta distribution parameter.  q must be gt 0.0.
c betai  the incomplete beta function ratio is the probability that a
c        random variable from a beta distribution having parameters
c        p and q will be less than or equal to x.
c
      double precision x, pin, qin, alneps, alnsml, c, eps, finsum, p,
     1  ps, p1, q, sml, term, xb, xi, y, dint, d1mach, dlbeta,
     2  dexp, dlog
      external d1mach, dexp, dint, dlbeta, dlog
      data eps, alneps, sml, alnsml / 4*0.0d0 /
c
      if (eps.ne.0.0d0) go to 10
      eps = d1mach(3)
      alneps = dlog (eps)
      sml = d1mach(1)
      alnsml = dlog (sml)
c
 10   if (x.lt.0.d0 .or. x.gt.1.d0) call seteru (
     1  35hdbetai  x is not in the range (0,1), 35, 1, 2)
      if (pin.le.0.d0 .or. qin.le.0.d0) call seteru (
     1  29hdbetai  a and/or b is le zero, 29, 2, 2)
c
      y = x
      p = pin
      q = qin
      if (q.le.p .and. x.lt.0.8d0) go to 20
      if (x.lt.0.2d0) go to 20
      y = 1.0d0 - y
      p = qin
      q = pin
c
 20   if ((p+q)*y/(p+1.d0).lt.eps) go to 80
c
c evaluate the infinite sum first.  term will equal
c y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) .
c
      ps = q - dint(q)
      if (ps.eq.0.d0) ps = 1.0d0
      xb = p*dlog(y) - dlbeta(ps,p) - dlog(p)
      dbetai = 0.0d0
      if (xb.lt.alnsml) go to 40
c
      dbetai = dexp (xb)
      term = dbetai*p
      if (ps.eq.1.0d0) go to 40
      n = max1 (sngl(alneps/dlog(y)), 4.0)
      do 30 i=1,n
        xi = i
        term = term * (xi-ps)*y/xi
        dbetai = dbetai + term/(p+xi)
 30   continue
c
c now evaluate the finite sum, maybe.
c
 40   if (q.le.1.0d0) go to 70
c
      xb = p*dlog(y) + q*dlog(1.0d0-y) - dlbeta(p,q) - dlog(q)
      ib = max1 (sngl(xb/alnsml), 0.0)
      term = dexp (xb - dble(float(ib))*alnsml )
      c = 1.0d0/(1.d0-y)
      p1 = q*c/(p+q-1.d0)
c
      finsum = 0.0d0
      n = q
      if (q.eq.dble(float(n))) n = n - 1
      do 50 i=1,n
        if (p1.le.1.0d0 .and. term/eps.le.finsum) go to 60
        xi = i
        term = (q-xi+1.0d0)*c*term/(p+q-xi)
c
        if (term.gt.1.0d0) ib = ib - 1
        if (term.gt.1.0d0) term = term*sml
c
        if (ib.eq.0) finsum = finsum + term
 50   continue
c
 60   dbetai = dbetai + finsum
 70   if (y.ne.x .or. p.ne.pin) dbetai = 1.0d0 - dbetai
      dbetai = dmax1 (dmin1 (dbetai, 1.0d0), 0.0d0)
      return
c
 80   dbetai = 0.0d0
      xb = p*dlog(dmax1(y,sml)) - dlog(p) - dlbeta(p,q)
      if (xb.gt.alnsml .and. y.ne.0.0d0) dbetai = dexp(xb)
      if (y.ne.x .or. p.ne.pin) dbetai = 1.0d0 - dbetai
c
      return
      end