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