Calc beta incomleate function by Lentz's algorithm

betainc.ts

type i32 = number;
type f64 = number;

function isNeg(x: f64): boolean {
  return x < 0.0 && x == Math.floor(x);
}

// Returns the value ln(gamma(x)).
function logGamma(x: f64): f64 {
  let t: f64, r: f64;
  t  = x + 5.5;
  t -= (x + 0.5) * Math.log(t);

  r  = 1.000000000190015;
  r += 76.180091729471460e+0 / (x + 1.0);
  r -= 86.505320329416770e+0 / (x + 2.0);
  r += 24.014098240830910e+0 / (x + 3.0);
  r -= 1.2317395724501550e+0 / (x + 4.0);
  r += 0.1208650973866179e-2 / (x + 5.0);
  r -= 0.5395239384953000e-5 / (x + 6.0);

  return Math.log(2.5066282746310005 * r / x) - t;

  // let d = [
  //    2.48574089138753565546e-5,
  //    1.05142378581721974210e+0,
  //   -3.45687097222016235469e+0,
  //    4.51227709466894823700e+0,
  //   -2.98285225323576655721e+0,
  //    1.05639711577126713077e+0,
  //   -1.95428773191645869583e-1,
  //    1.70970543404441224307e-2,
  //   -5.71926117404305781283e-4,
  //    4.63399473359905636708e-6,
  //   -2.71994908488607703910e-9
  // ];

  // let s = d[0], z = x - 1.0;
  // for (let k = 1; k <= 10; k++) {
  //   s += d[k] / (z + k);
  // }
  // z += 0.5;
  // return Math.log(1.860382734205265717 * s) - z + z * Math.log1p(x + 9.400511);
}

function betainc(a: f64, b: f64, x: f64): f64 {
  if (x  < 0.0 || x > 1.0) return NaN;
  if (isNeg(a) || isNeg(b) || isNeg(a + b)) return NaN;

  const EPS = 1e-9;
  if (x < EPS) return 0;
  if (x > 1 - EPS) return 1;
  if (a < EPS && b < EPS) return 0.5;
  if (a < EPS) return 1;
  if (b < EPS) return 0;

  // According to https://dlmf.nist.gov/8.17#v, the continued fraction converges rapidly for x < (a + 1) / (a + b + 2)
  // If x>=(a + 1) / (a + b + 2), we use Ix(a, b) = I(1 - x, b, a)
  if (x > (a + 1.0) / (a + b + 2.0)) {
    return 1.0 - betainc(b, a, 1.0 - x);
  }

  // We use the Lentz's algorithm to calculate the continued fraction
  //   f = 1 + d[1] / (1 + d[2] / (1 + d[3] / (...)))
  // The implementation is based on the following formulas:
  //   u[0] = 1, v[0] = 0, f[0] = 1
  //   u[i] = 1 + d[i] / u[i - 1]
  //   v[i] = 1 / (1 + d[i] * v[i - 1])
  //   f[i] = f[i - 1] * u[i] * v[i]

  const MAX_ITERS = 400;
  const norm = (z: f64) => (
    Math.abs(z) < 1e-30 ? 1e-30 : z
  );

  function logBeta(a: f64, b: f64): f64 {
    return logGamma(a) + logGamma(b) - logGamma(a + b);
  }

  let u = 2.0, v = 1.0, f = 2.0, d = 1.0; // d[0]
  for (let i = 1; i <= MAX_ITERS; i++) {
    let m   = i >>> 1;
    let am  = a + m;
    let am2 = am + m;
    if ((i & 1) === 0) {
      d = m * (b - m) * x / ((am2 - 1.0) * am2); // d[2m]
    } else {
      d = -(am * (am + b) * x) / (am2 * (am2 + 1.0)); // d[2m + 1]
    }
    u = norm(1.0 + d / u);
    v = 1.0 / norm(1.0 + d * v);
    let uv = u * v;
    f *= uv;

    if (Math.abs(uv - 1.0) <= Number.EPSILON) break;
  }
  // Ix(a, b) = x^a * (1-x)^b / (a*B(a, b)) * 1 / (1 + d[1] / (1 + d[2] / (1 + d[3] / (...))))
  return Math.exp(Math.log(x) * a + Math.log1p(-x) * b - logBeta(a, b)) / a * (f - 1.0);
}

function test(a: f64, b: f64, x: f64, expected: f64) {
  let res = betainc(a, b, x);
  console.log(`\nactual   betainc(${a}, ${b}, ${x}) = ${res}`);
  console.log(`expected betainc(${a}, ${b}, ${x}) = ${expected}`);
  console.log(`|error| = ${ Math.abs((expected || 0) - (res || 0)).toExponential(1) }`);
}

test(0.5, 1.0, 0.0, 0.0);
test(0.5, 0.5, 1.0, 1.0);
test(1.0, 3.0, 0.5, 0.875);
test(1.4, 3.1, 0.5, 0.8148904036225296);
test(0.5, 5.0, 0.2, 0.8550723945959197);
test(0.5, 5.0, 0.5, 0.9898804402645662);
test(0.001, 0.3333, 0.1111, 0.9953402935460681);
test(0.999, 1.567, 1.46e-8, 2.3287434894343686e-08);
test(1.234e-8, 1.0, 0.1, 0.9999999715861003);
test(128.5, 3.0, 0.5, 4.4580244557709776e-36);
test(0.5, 4.0, 0.9999999, 0.9999999999999999);
test(20, 30, 0.001, 2.750687665855991e-47);
test(20, 30, 0.0001, 2.819953178893307e-67);
test(0.5, 33.0, 1.0000001, NaN);

/* Results:

actual   betainc(0.5, 1, 0) = 0
expected betainc(0.5, 1, 0) = 0
|error| = 0.0e+0

actual   betainc(0.5, 0.5, 1) = 1
expected betainc(0.5, 0.5, 1) = 1
|error| = 0.0e+0

actual   betainc(1, 3, 0.5) = 0.8749999999999999
expected betainc(1, 3, 0.5) = 0.875
|error| = 1.1e-16

actual   betainc(1.4, 3.1, 0.5) = 0.8148904036225282
expected betainc(1.4, 3.1, 0.5) = 0.8148904036225296
|error| = 1.4e-15

actual   betainc(0.5, 5, 0.2) = 0.8550723945958834
expected betainc(0.5, 5, 0.2) = 0.8550723945959197
|error| = 3.6e-14

actual   betainc(0.5, 5, 0.5) = 0.9898804402645667
expected betainc(0.5, 5, 0.5) = 0.9898804402645662
|error| = 4.4e-16

actual   betainc(0.001, 0.3333, 0.1111) = 0.9953402935460702
expected betainc(0.001, 0.3333, 0.1111) = 0.9953402935460681
|error| = 2.1e-15

actual   betainc(0.999, 1.567, 1.46e-8) = 2.3287434894343845e-8
expected betainc(0.999, 1.567, 1.46e-8) = 2.3287434894343686e-8
|error| = 1.6e-22

actual   betainc(1.234e-8, 1, 0.1) = 0.9999999715861008
expected betainc(1.234e-8, 1, 0.1) = 0.9999999715861003
|error| = 4.4e-16

actual   betainc(128.5, 3, 0.5) = 4.458024455776908e-36
expected betainc(128.5, 3, 0.5) = 4.4580244557709776e-36
|error| = 5.9e-48

actual   betainc(0.5, 4, 0.9999999) = 1
expected betainc(0.5, 4, 0.9999999) = 0.9999999999999999
|error| = 1.1e-16

actual   betainc(20, 30, 0.001) = 2.7506876659140094e-47
expected betainc(20, 30, 0.001) = 2.750687665855991e-47
|error| = 5.8e-58

actual   betainc(20, 30, 0.0001) = 2.8199531789528524e-67
expected betainc(20, 30, 0.0001) = 2.819953178893307e-67
|error| = 6.0e-78

actual   betainc(0.5, 33, 1.0000001) = NaN
expected betainc(0.5, 33, 1.0000001) = NaN
|error| = 0.0e+0

*/