      function beta_i(x,a,b) ! bernd berg, apr 4 2000.
c incomplete beta function. 
      include '../../libs/fortran/implicit.sta'
      include '../../libs/fortran/constants.par'
      parameter(iter_max=200,eps=1.d-10)
 
      if(x.lt.zero .or. x.gt.one) stop 'bad argument x in beta_i'
      if(x.eq.zero .or. x.eq.one) then
        bt=zero
      else
        bt=exp(gamma_ln(a+b)-gamma_ln(a)-gamma_ln(b)
     &    +a*log(x)+b*log(one-x))
      endif
 
      if(x.lt.(a+one)/(a+b+two)) then
        xx=x
        aa=a
        bb=b
      else
        xx=one-x
        aa=b
        bb=a
      end if

      apb=aa+bb
      ap1=aa+one
      am1=aa-one
 
      bcfm=one
      bm=one
      bcf=one
      bz=one-apb*xx/ap1
      do iter=1,iter_max
        xiter=iter*one
        two_iter=xiter+xiter
        c1=xiter*(bb-xiter)*xx/((am1+two_iter)*(aa+two_iter))
        c2=-(aa+xiter)*(apb+xiter)*xx/((aa+two_iter)*(ap1+two_iter))
        bcfp=bcf+c1*bcfm
        bp=bz+c1*bm
        bpp=bp+c2*bz
        bcfold=bcf
        bcfm=bcfp/bpp
        bm=bp/bpp
        bcf=(bcfp+c2*bcf)/bpp
        bz=one
        if(abs(bcf-bcfold).lt.eps*abs(bcf)) go to 1
      end do
      stop 'beta_i: a or b too big, or iter_max too small'
c
1     if(x.lt.(a+one)/(a+b+two)) then
        beta_i=bt*bcf/a
      else
        beta_i=one-bt*bcf/b
      endif
      return
      end
