docs: def _gcf(a,x):
      "Continued fraction representation of Gamma. NumRec sect 6.1"
      ITMAX=100
      EPS=3.e-7
      FPMIN=1.e-30
  
      gln=lgamma(a)
      b=x+1.-a
      c=1./FPMIN
      d=1./b
      h=d
      for i in range(1,ITMAX+1):
          an=-i*(i-a)
          b=b+2.
          d=an*d+b
          if abs(d) < FPMIN: d=FPMIN
          c=b+an/c
          if abs(c) < FPMIN: c=FPMIN
          d=1./d
          delt=d*c
          h=h*delt
          if abs(delt-1.) < EPS: break
      else:
          print('a too large, ITMAX too small in gcf')
      gammcf=np.exp(-x+a*np.log(x)-gln)*h
      return gammcf,gln