code: def betacf(a, b, x):
      
      ITMAX = 200
      EPS = 3.0e-7
  
      bm = az = am = 1.0
      qab = a + b
      qap = a + 1.0
      qam = a - 1.0
      bz = 1.0 - qab * x / qap
      for i in range(ITMAX + 1):
          em = float(i + 1)
          tem = em + em
          d = em * (b - em) * x / ((qam + tem) * (a + tem))
          ap = az + d * am
          bp = bz + d * bm
          d = -(a + em) * (qab + em) * x / ((qap + tem) * (a + tem))
          app = ap + d * az
          bpp = bp + d * bz
          aold = az
          am = ap / bpp
          bm = bp / bpp
          az = app / bpp
          bz = 1.0
          if (abs(az - aold) < (EPS * abs(az))):
              return az
      print('a or b too big, or ITMAX too small in Betacf.')
docs: This function evaluates the continued fraction form of the incomplete
      Beta function, betai.  (Adapted from: Numerical Recipies in C.)
  
      Usage:   lbetacf(a,b,x)