update efficiency of calculating binomial coefficients

[blazeshomida]

I can make a formal PR/Issue, but just sitting on my phone talking with chatGPT trying to figure out how to efficiently calculate binomial coefficients,I was looking at how you solve binomial coefficients in this package; @stdlib/math-base-special-binomcoef, and was wondering if a solution such as this would be better?

function binomialCoefficientLog(n, k) {
  if (k > n) return 0;
  if (k === 0 || k === n) return 1;

  let logSum = 0;
  for (let i = 0; i < k; i++) {
    logSum += Math.log(n - i) - Math.log(i + 1);
  }

  return Math.exp(logSum);
}

This basically breaks it down to it's most simplest form (with log for maintaining precision) by stopping the loop once we reach k, so only calculating necessary values. Basically O(k) time complexity and O(1) space complexity.

Obviously error handling would need to be added but in a nutshell this is what I was thinking.

[stdlib-bot]

:wave: Hi there! :wave:

And thank you for opening your first issue! We will get back to you shortly. :runner: :dash:

[kgryte]

@blazeshomida Thanks for opening this issue. The main concern I have for your proposed implementation is precision loss. Especially for large n and k, you will likely have significant error accumulation. The current implementation attempts to minimize error.

[blazeshomida]

@kgryte So I've just done some simple testing in codesandbox, and I see your point with precision loss because of the use of logarithms. However, with a slight revision, you can achieve a result closer to what you're currently using, although I have seen slight differences near the end with very large numbers.

This approach sacrifices a bit of precision when reaching large numbers, but any result smaller than that is accurate. It's not so much an update to make it better, just a slightly different approach, trading off precision for efficiency.

function binomialCoefficient(n: number, k: number) {
  if (k > n) return 0;
  if (k === 0 || k === n) return 1;
  if (n - k < k) k = n - k;

  let result = 1;
  for (let i = 0; i < k; i++) {
    result *= (n - i) / (i + 1);
  }

  // Precision check
  if (Math.abs(result - Math.round(result)) > Number.EPSILON) {
    console.warn("Result may suffer from precision loss.");
  }

  return result;
}

This function leverages symmetry and calculates only necessary values up to k, achieving O( k ) time complexity and O( 1 ) space complexity. While it may introduce slight precision loss for very large numbers, it offers an efficient alternative for most practical purposes.

If you want to close this issue, feel free. After digging into it a bit more I understand the reasonings for the route you chose.

[kgryte]

Some related discussion: https://github.com/stdlib-js/stdlib/issues/1155.

There may be an opportunity for adding a "fast" and less accurate version to @stdlib/math/base/special/fast, which is where we house other similar functions. In that case, I'd be slightly more inlined to simply use our older implementation (see https://github.com/stdlib-js/stdlib/blob/d54143e5d29e8c320914df3e989dab7f923a5b31/lib/node_modules/%40stdlib/math/base/special/binomcoef/lib/main.js).