[RFC]: implementation modification to improve accuracy of `@stdlib/math/base/special/binomcoef` about stdlib HOT 3 CLOSED

@Andrew-Cottrell I did a bit more investigating. I was able to update the implementation to only compute the GCD when overflow is a risk. For larger binomial coefficients, this drops performance by 50%. Not great, but also not the end of the world.

For reference:

Thanks for opening this issue and helping point in the right direction. This resolves the failing tests in fibpoly, which have subsequently been updated.

from stdlib.

An alternative version of the modification avoids calling gcd when unnecessary, at the cost of increased code complexity.

	var res;
	var j;
	var s; // 's' for safe
	var d; // 'd' for divisor

        ...

	s = Math.floor(Number.MAX_SAFE_INTEGER / n);
	res = n;
	for ( j = 2; j <= k && res < s; j++ ) {
		res = ( res * ( n - j + 1 ) ) / j;
	}
	for ( ; j <= k; j++ ) {
		d = gcd( res, j );
		res = ( res / d ) * ( d * (n - j + 1) / j );
	}
	return res;

The next version swaps the order of the numerator values to be increasing rather than decreasing, which helps performance by sometimes evaluating more terms using the non-gcd version of the calculation.

	s = Math.floor(Number.MAX_SAFE_INTEGER / n);
	res = n - k + 1;
	n = res + 1;
	for ( j = 2; j <= k && res < s; j++, n++ ) {
		res = (res * n) / j;
	}
	for ( ; j <= k; j++, n++ ) {
		d = gcd( res, j );
		res = ( res / d ) * ( d * n / j );
	}
	return res;

In some non-scientific testing, this version appears to be about three times faster, than the simple modification, when computing the four example binomial coefficients mentioned above. And about seven times faster when computing smaller binomial coefficients, such as binomcoef( 54, 20 ), that don't need to call the gcd function.

from stdlib.

@Andrew-Cottrell Thanks for opening this issue and your proposed solutions. Based on my testing, performance drops by over 10x when moving to GCD, even for the improved version.

The accuracy does improve, which is a plus, but I'd prefer, if possible, to find a different algorithm which doesn't force an order of magnitude drop in performance, especially as binomcoef is used as a dependency by other functions within stdlib. Meaning, a significant drop in binomcoef perf will likely necessarily mean a drop in dependants perf.

I've done a bit of research to try and find other algorithms which improve accuracy without sacrificing perf; however, I haven't had much luck finding one which is universally good.

from stdlib.