code: func BetaInc(x, a, b float64) float64 {
  	// Based on Numerical Recipes in C, section 6.4. This uses the
  	// continued fraction definition of I:
  	//
  	//  (xᵃ*(1-x)ᵇ)/(a*B(a,b)) * (1/(1+(d₁/(1+(d₂/(1+...))))))
  	//
  	// where B(a,b) is the beta function and
  	//
  	//  d_{2m+1} = -(a+m)(a+b+m)x/((a+2m)(a+2m+1))
  	//  d_{2m}   = m(b-m)x/((a+2m-1)(a+2m))
  	if x < 0 || x > 1 {
  		return math.NaN()
  	}
  	bt := 0.0
  	if 0 < x && x < 1 {
  		// Compute the coefficient before the continued
  		// fraction.
  		bt = math.Exp(lgamma(a+b) - lgamma(a) - lgamma(b) +
  			a*math.Log(x) + b*math.Log(1-x))
  	}
  	if x < (a+1)/(a+b+2) {
  		// Compute continued fraction directly.
  		return bt * betacf(x, a, b) / a
  	} else {
  		// Compute continued fraction after symmetry transform.
  		return 1 - bt*betacf(1-x, b, a)/b
  	}
  }
docs: // BetaInc returns the value of the regularized incomplete beta
  // function Iₓ(a, b).
  //
  // This is not to be confused with the "incomplete beta function",
  // which can be computed as BetaInc(x, a, b)*Beta(a, b).
  //
  // If x < 0 or x > 1, returns NaN.