function betai (x, pin, qin)
c august 1980 version.  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
      external aint, albeta, alog, exp, r1mach
      data eps, alneps, sml, alnsml / 4*0.0 /
c
      if (eps.ne.0.) go to 10
      eps = r1mach(3)
      alneps = alog(eps)
      sml = r1mach(1)
      alnsml = alog(sml)
c
 10   if (x.lt.0. .or. x.gt.1.0) call seteru (
     1  35hbetai   x is not in the range (0,1), 35, 1, 2)
      if (pin.le.0. .or. qin.le.0.) call seteru (
     1  29hbetai   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.8) go to 20
      if (x.lt.0.2) go to 20
      y = 1.0 - y
      p = qin
      q = pin
c
 20   if ((p+q)*y/(p+1.).lt.eps) go to 80
c
c evaluate the infinite sum first.
c term will equal y**p/beta(ps,p) * (1.-ps)i * y**i / fac(i)
c
      ps = q - aint(q)
      if (ps.eq.0.) ps = 1.0
      xb = p*alog(y) -  albeta(ps, p) - alog(p)
      betai = 0.0
      if (xb.lt.alnsml) go to 40
c
      betai = exp (xb)
      term = betai*p
      if (ps.eq.1.0) go to 40
c
      n = amax1 (alneps/alog(y), 4.0)
      do 30 i=1,n
        term = term*(float(i)-ps)*y/float(i)
        betai = betai + term/(p+float(i))
 30   continue
c
c now evaluate the finite sum, maybe.
c
 40   if (q.le.1.0) go to 70
c
      xb = p*alog(y) + q*alog(1.0-y) - albeta(p,q) - alog(q)
      ib = amax1 (xb/alnsml, 0.0)
      term = exp (xb - float(ib)*alnsml)
      c = 1.0/(1.0-y)
      p1 = q*c/(p+q-1.)
c
      finsum = 0.0
      n = q
      if (q.eq.float(n)) n = n - 1
      do 50 i=1,n
        if (p1.le.1.0 .and. term/eps.le.finsum) go to 60
        term = (q-float(i-1))*c*term/(p+q-float(i))
c
        if (term.gt.1.0) ib = ib - 1
        if (term.gt.1.0) term = term*sml
c
        if (ib.eq.0) finsum = finsum + term
 50   continue
c
 60   betai = betai + finsum
 70   if (y.ne.x .or. p.ne.pin) betai = 1.0 - betai
      betai = amax1 (amin1 (betai, 1.0), 0.0)
      return
c
 80   betai = 0.0
      xb = p*alog(amax1(y,sml)) - alog(p) - albeta(p,q)
      if (xb.gt.alnsml .and. y.ne.0.) betai = exp (xb)
      if (y.ne.x .or. p.ne.pin) betai = 1.0 - betai
      return
c
      end