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,  dlbeta
        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