export default function logbeta(a: number, b: number): number {
  // Log Beta Function using the relationship with Gamma function
  // ln(Beta(a,b)) = ln(Gamma(a)) + ln(Gamma(b)) - ln(Gamma(a+b))
  return lnGamma(a) + lnGamma(b) - lnGamma(a + b);
}

// Efficient implementation of ln(Gamma(x)) using Lanczos approximation
function lnGamma(x: number): number {
  const p: number[] = [
    676.5203681218851, -1259.1392167224028, 771.32342877765313,
    -176.61502916214059, 12.507343278686905, -0.13857109526572012,
    9.9843695780195716e-6, 1.5056327351493116e-7,
  ];

  if (x < 0.5) {
    return Math.log(Math.PI) - Math.log(Math.sin(Math.PI * x)) - lnGamma(1 - x);
  }

  // eslint-disable-next-line no-param-reassign
  x -= 1;
  const y: number = x + 7.5;
  let sum = 0.99999999999980993;
  // eslint-disable-next-line no-plusplus
  for (let i = 0; i < p.length; i++) {
    sum += p[i] / (x + i + 1);
  }

  return (
    Math.log(Math.sqrt(2 * Math.PI)) +
    (x + 0.5) * Math.log(y) -
    y +
    Math.log(sum)
  );
}